25 #include "reference_element_aux.icc"
51 _bkm1_node_internal_d()
66 const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&
69 base::_initialize_data_guard (hat_K);
70 return _hat_node [hat_K.
variant()];
73 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
76 base::_initialize_data_guard (hat_K);
80 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
83 base::_initialize_data_guard (hat_K);
84 return _inv_vdm [hat_K.
variant()];
91 for (
size_type map_d = 0; map_d < 4; ++map_d) {
92 base::_ndof_on_subgeo_internal [map_d].fill (0);
109 base::_helper_make_discontinuous_ndof_on_subgeo (base::is_continuous(), base::_ndof_on_subgeo_internal, base::_ndof_on_subgeo);
110 base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_ndof_on_subgeo_internal, base::_first_idof_by_dimension_internal);
111 base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_ndof_on_subgeo, base::_first_idof_by_dimension);
121 qopt.set_order (2*(k+2));
123 for (
size_type map_d = 0; map_d < 4; ++map_d) {
124 base::_nnod_on_subgeo_internal [map_d].fill (0);
137 base::_helper_make_discontinuous_ndof_on_subgeo (base::is_continuous(), base::_nnod_on_subgeo_internal, base::_nnod_on_subgeo);
138 base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_nnod_on_subgeo_internal, base::_first_inod_by_dimension_internal);
139 base::_helper_initialize_first_ixxx_by_dimension_from_nxxx_on_subgeo (base::_nnod_on_subgeo, base::_first_inod_by_dimension);
156 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& hat_node = _hat_node [hat_K.
variant()];
157 hat_node.resize (base::_first_inod_by_dimension [
variant][
d+1]);
161 for (
size_type loc_isid = 0, loc_nsid = hat_K.
n_subgeo(
d-1); loc_isid < loc_nsid; ++loc_isid) {
165 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> hat_node_sid;
167 switch (base::_sopt.get_node()) {
176 error_macro (
"unsupported node set: " << base::_sopt.get_node_name());
178 point loc_vertex [4];
179 for (
size_type loc_jsidvert = 0; loc_jsidvert < loc_nsidvert; ++loc_jsidvert) {
181 loc_vertex[loc_jsidvert] = hat_K.
vertex (loc_jvertex);
183 for (
size_type loc_inod_sid = 0, loc_nnod_sid = hat_node_sid.size(); loc_inod_sid < loc_nnod_sid; ++loc_inod_sid) {
184 check_macro (loc_inod <
size_t(hat_node.size()),
"invalid loc_inod");
185 T xi0 = hat_node_sid [loc_inod_sid][0],
186 xi1 = hat_node_sid [loc_inod_sid][1];
187 if (loc_nsidvert == 4) {
192 hat_node[loc_inod][
alpha] = loc_vertex [0][
alpha];
194 hat_node[loc_inod][
alpha] += xi0*(loc_vertex[1][
alpha] - loc_vertex[0][
alpha]);
197 hat_node[loc_inod][
alpha] += xi1*(loc_vertex[loc_nsidvert-1][
alpha] - loc_vertex[0][
alpha]);
207 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> hat_node_internal;
209 switch (base::_sopt.get_node()) {
218 error_macro (
"unsupported node set: " << base::_sopt.get_node_name());
220 trace_macro (
"hat_node_internal.size="<<hat_node_internal.size());
221 for (
size_type loc_inod_int = 0; loc_inod_int <
size_type(hat_node_internal.size()); ++loc_inod_int) {
222 trace_macro (
"hat_node_internal["<<loc_inod_int<<
"]="<<hat_node_internal[loc_inod_int]);
223 hat_node[loc_inod] = hat_node_internal [loc_inod_int];
228 last_q = _quad.end(hat_K); iter_q != last_q; iter_q++) {
229 hat_node[loc_inod] = (*iter_q).x;
234 check_macro (loc_inod == loc_nnod,
"invalid node count: loc_inod="<<loc_inod<<
" and loc_nnod="<<loc_nnod);
244 size_type dim_Pkp1 = _b_pre_kp1.ndof (hat_K);
246 for (
size_type loc_isid = 0, loc_nsid = hat_K.
n_subgeo(
d-1); loc_isid < loc_nsid; ++loc_isid) {
249 loc_ndof_sid_tot += base::_nnod_on_subgeo_internal [
d][hat_S.
variant()];
252 size_type loc_nnod_sid_tot = loc_ndof_sid_tot;
256 warning_macro (
"first_inod(hat_K,d) ="<<base::first_inod(hat_K,
d));
257 warning_macro (
"first_inod(hat_K,d+1)="<<base::first_inod(hat_K,
d+1));
263 warning_macro (
"first_idof(hat_K,d-1)="<<base::first_idof(hat_K,
d-1));
264 warning_macro (
"first_idof(hat_K,d) ="<<base::first_idof(hat_K,
d));
265 warning_macro (
"first_idof(hat_K,d+1)="<<base::first_idof(hat_K,
d+1));
268 warning_macro (
"a_tilde(loc_ndof)="<<loc_ndof<<
",d*dim(Pkp1)="<<
d*dim_Pkp1<<
")");
275 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> tilde_a (loc_ndof,
d*dim_Pkp1);
281 trace_macro (
"basis(1): (q,0) and (0,r), q,r in basis(Pk)...");
282 for (
size_type loc_comp_idof = 0; loc_comp_idof < dim_Pk; ++loc_comp_idof) {
286 trace_macro (
"loc_idof="<<loc_idof<<
", loc_jpkp1="<<loc_jpkp1);
287 tilde_a (loc_idof, loc_jpkp1) = 1;
291 trace_macro (
"basis(2): (x0*p,x1*p), p in basis(P_k)\basis(P_{k-1}) [exactly k degree]");
292 std::vector<size_type> inod2ideg_kp1;
293 build_inod2ideg (hat_K, k+1, inod2ideg_kp1);
302 std::vector<point_basic<size_type> > power_index;
303 make_power_indexes_sorted_by_degrees (hat_K, k, power_index);
304 for (
size_type loc_ideg = dim_Pkm1, loc_ndeg = power_index.size(); loc_ideg < loc_ndeg; ++loc_ideg) {
307 size_type loc_idof =
d*dim_Pk + (loc_ideg - dim_Pkm1);
309 size_type loc_inod_kp1 = ilat2loc_inod (hat_K, k+1, ilat);
310 size_type loc_ideg_kp1 = inod2ideg_kp1 [loc_inod_kp1];
312 trace_macro (
"loc_idof="<<loc_idof<<
", loc_jpkp1d="<<loc_jpkp1d);
313 tilde_a (loc_idof, loc_jpkp1d) = 1;
328 std::vector<point_basic<size_type> > power_index_sub;
329 make_power_indexes_sorted_by_degrees (hat_sub, k, power_index_sub);
330 for (
size_type loc_ideg = 0, loc_ndeg = power_index_sub.size(); loc_ideg < loc_ndeg; ++loc_ideg) {
336 ilat [(
alpha+1)%
d] = ilat_sub[0];
337 if (
d == 3) ilat [(
alpha+2)%
d] = ilat_sub[1];
338 size_type loc_inod_kp1 = ilat2loc_inod (hat_K, k+1, ilat);
339 size_type loc_ideg_kp1 = inod2ideg_kp1 [loc_inod_kp1];
341 trace_macro (
"loc_idof="<<loc_idof<<
", loc_jpkp1d="<<loc_jpkp1d);
342 tilde_a (loc_idof, loc_jpkp1d) = 1;
348 default:
error_macro (
"unexpected element `"<<hat_K.
name()<<
"' (HINT: see BerDur-2013)");
352 cout << setprecision(std::numeric_limits<T>::digits10)
353 <<
"tilde_a=[" << endl << tilde_a <<
"]"<<endl
354 <<
"[u,s,vt]=svd(tilde_a)" << endl
355 <<
"id=eye(size(s))" << endl
356 <<
"a=u*id*vt" << endl
366 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
367 u (tilde_a.rows(), tilde_a.rows()),
368 vt (tilde_a.cols(), tilde_a.cols());
369 Eigen::Matrix<T,Eigen::Dynamic,1> s (tilde_a.rows());
370 Eigen::JacobiSVD<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >
371 svd (tilde_a, Eigen::ComputeFullU | Eigen::ComputeFullV);
373 s = svd.singularValues();
379 rank_s += (abs(s[loc_idof]) > eps) ? 1 : 0;
382 "invalid polynomial space dimension = " << rank_s <<
" < loc_ndof = " << loc_ndof);
383 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> id (tilde_a.rows(), tilde_a.cols());
385 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
a =
id*vt.transpose();
388 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> s1 (tilde_a.rows(), tilde_a.cols());
390 for (
size_type iloc = 0; iloc < size_t(tilde_a.rows()); iloc++) {
391 s1(iloc,iloc) = s(iloc);
393 cout <<
"s1 = ["<< endl << s1 <<
"];" << endl
394 <<
"u1 = ["<< endl <<
u <<
"];"<<endl
395 <<
"vt1 = ["<< endl << vt <<
"];"<<endl
396 <<
"id1 = ["<< endl <<
id <<
"];"<<endl
397 <<
"a1 = ["<< endl <<
a <<
"]"<<endl
398 <<
"err=norm(tilde_a-u*s*vt')"<<endl
399 <<
"err1=norm(tilde_a-u1*s1*vt1')"<<endl
400 <<
"err_a=abs(a-a1)"<<endl
401 <<
"err_svd=max(max(abs(a-a1)))"<<endl
402 <<
"err_u=max(max(abs(u-u1)))"<<endl
403 <<
"err_v=max(max(abs(vt-vt1)))"<<endl
404 <<
"err_s=max(diag(s1)-diag(s))"<<endl
406 T err_svd = (tilde_a -
u*s1*vt.transpose()).
norm();
407 cout <<
"err_svd = " << err_svd << endl;
418 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
420 bkm1_node_internal_d = _bkm1_node_internal_d [
variant];
423 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> hat_node_internal (loc_nnod_int);
424 for (
size_type loc_inod_int = 0; loc_inod_int < loc_nnod_int; ++loc_inod_int) {
425 size_type loc_inod = loc_nnod_sid_tot + loc_inod_int;
426 hat_node_internal [loc_inod_int] = hat_node [loc_inod];
428 string basis_dof_name;
429 switch (base::_sopt.get_raw_polynomial()) {
433 default:
error_macro (
"unsupported raw polynomial basis `"<<base::_sopt.get_raw_polynomial_name()<<
"'");
445 bkm1_node_internal_d [
alpha] = bkm1_node_internal_d [0];
453 b_dof_k (basis_dof_name+std::to_string(k));
454 std::array<Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>,3> hat_node_internal_comp;
456 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
458 bkm1_node_internal_comp,
459 bk_node_internal_comp;
462 hat_node_internal_comp [
alpha].resize (hat_node_internal.size());
463 for (
size_type loc_inod_int = 0, loc_nnod_int = hat_node_internal.size(); loc_inod_int < loc_nnod_int; ++loc_inod_int) {
465 hat_node_internal_comp [
alpha][loc_inod_int][0] = (1+hat_node_internal [loc_inod_int][
alpha])/2;
470 size_type loc_ndof_int = bkm1_node_internal_comp[0].cols()*
pow(bk_node_internal_comp[0].cols(),
d-1);
474 bkm1_node_internal_d [
alpha].resize (hat_node_internal.size(), loc_ndof_int);
475 for (
size_type loc_inod_int = 0, loc_nnod_int = hat_node_internal.size(); loc_inod_int < loc_nnod_int; ++loc_inod_int) {
480 bkm1_node_internal_d [
alpha ] (loc_inod_int,loc_idof_int)
481 = bkm1_node_internal_comp [
alpha ] (loc_inod_int,i)
482 * bk_node_internal_comp [alpha2] (loc_inod_int,j);
489 bkm1_node_internal_d [
alpha ] (loc_inod_int,loc_idof_int)
490 = bkm1_node_internal_comp [
alpha ] (loc_inod_int,i)
491 * bk_node_internal_comp [alpha2] (loc_inod_int,j)
492 * bk_node_internal_comp [alpha3] (loc_inod_int,k);
507 size_type loc_n_bkp1 = _b_pre_kp1.ndof (hat_K);
508 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> bkp1_node (loc_nnod, loc_n_bkp1);
512 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> bkp1d_dof (loc_ndof,
d*loc_n_bkp1);
513 bkp1d_dof.fill (std::numeric_limits<T>::max());
514 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> bkp1d_j_node (loc_nnod);
515 Eigen::Matrix<T,Eigen::Dynamic,1> bkp1d_j_dof (loc_ndof);
516 for (
size_type loc_j_bkp1 = 0; loc_j_bkp1 < loc_n_bkp1; ++loc_j_bkp1) {
518 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
520 bkp1d_j_node [loc_inod][
alpha] = bkp1_node(loc_inod,loc_j_bkp1);
522 bkp1d_j_dof.fill (std::numeric_limits<T>::max());
523 _compute_dofs (hat_K, bkp1d_j_node, bkp1d_j_dof);
525 check_macro (bkp1d_dof.rows() == bkp1d_j_dof.size(),
"invalid sizes");
526 bkp1d_dof.col (loc_j_bkp1d) = bkp1d_j_dof;
535 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& vdm = _vdm [hat_K.
variant()];
536 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& inv_vdm = _inv_vdm[hat_K.
variant()];
537 vdm = bkp1d_dof*
a.transpose();
538 bool invert_ok =
invert(vdm, inv_vdm);
540 "unisolvence failed for " <<
base::name() <<
"(" << hat_K.
name() <<
") basis");
544 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& bar_a = _bar_a [hat_K.
variant()];
545 bar_a = inv_vdm.transpose()*
a;
547 cout <<
"bkp1d_dof=[" << endl << bkp1d_dof <<
"]"<<endl
548 <<
"vdm=[" << endl << vdm <<
"]"<<endl
549 <<
"det_vdm=" << vdm.determinant() <<endl
550 <<
"cond_vdm=" <<
cond(vdm) <<endl
551 <<
"bar_a1=inv(vdm)'*a;"<<endl
552 <<
"bar_a=[" << endl << bar_a <<
"]"<<endl;
563 Eigen::Matrix<value_type,Eigen::Dynamic,1>&
value)
const
565 base::_initialize_data_guard (hat_K);
571 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& bar_a = _bar_a [hat_K.
variant()];
572 Eigen::Matrix<T,Eigen::Dynamic,1> bkp1;
573 _b_pre_kp1.evaluate (hat_K, hat_x, bkp1);
574 value.resize (loc_ndof);
585 for (
size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
587 for (
size_type loc_jdof_bkp1 = 0, loc_ndof_bkp1 = bkp1.size(); loc_jdof_bkp1 < loc_ndof_bkp1; ++loc_jdof_bkp1) {
590 value[loc_idof][
alpha] += bar_a(loc_idof,loc_jdof_bkp1d)*bkp1[loc_jdof_bkp1];
602 base::_initialize_data_guard (hat_K);
608 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& bar_a = _bar_a [hat_K.
variant()];
609 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> grad_bkp1;
610 _b_pre_kp1.grad_evaluate (hat_K, hat_x, grad_bkp1);
611 value.resize (loc_ndof);
623 for (
size_type loc_idof = 0; loc_idof < loc_ndof; ++loc_idof) {
625 for (
size_type loc_jdof_bkp1 = 0, loc_ndof_bkp1 = grad_bkp1.size(); loc_jdof_bkp1 < loc_ndof_bkp1; ++loc_jdof_bkp1) {
629 value[loc_idof](
alpha,
beta) += bar_a(loc_idof,loc_jdof_bkp1d)*grad_bkp1[loc_jdof_bkp1][
beta];
651 const Eigen::Matrix<value_type,Eigen::Dynamic,1>& f_xnod,
652 Eigen::Matrix<T,Eigen::Dynamic,1>& dof)
const
655 base::_initialize_data_guard (hat_K);
660 dof.resize (loc_ndof);
666 for (
size_type loc_isid = 0, loc_nsid = hat_K.
n_subgeo(
d-1); loc_isid < loc_nsid; ++loc_isid) {
672 for (
size_type loc_idof_sid = 0; loc_idof_sid < loc_ndof_sid; ++loc_idof_sid) {
673 dof[loc_idof] = hat_S_meas*
dot(f_xnod[loc_inod],hat_n);
681 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
683 bkm1_node_internal_d = _bkm1_node_internal_d [hat_K.
variant()];
686 size_type loc_ndof_int_d =
d*bkm1_node_internal_d[0].cols();
688 check_macro (loc_ndof == loc_ndof_boundary + loc_ndof_int_d,
689 "invalid internal dof count: loc_ndof="<<loc_ndof
690 <<
", loc_ndof_boundary="<<loc_ndof_boundary
691 <<
", loc_ndof_int_d="<<loc_ndof_int_d);
696 for (
size_type loc_idof_int = 0, loc_ndof_int = bkm1_node_internal_d[0].cols(); loc_idof_int < loc_ndof_int; ++loc_idof_int) {
697 size_type loc_inod = first_loc_inod_int + loc_idof_int;
699 dof [loc_idof] = f_xnod [loc_inod][
alpha];
704 for (
size_type loc_idof_int = 0, loc_ndof_int = bkm1_node_internal_d[0].cols(); loc_idof_int < loc_ndof_int; ++loc_idof_int) {
706 loc_inod = first_loc_inod_int;
710 last_q = _quad.end(hat_K); iter_q != last_q; iter_q++, inod_q++, ++loc_inod) {
711 sum += f_xnod [loc_inod][
alpha]
712 *bkm1_node_internal_d[
alpha] (inod_q, loc_idof_int)
716 dof [loc_idof] = sum;
721 check_macro (loc_idof == loc_ndof,
"invalid dof count");
727 #define _RHEOLEF_instanciation(T) \
728 template class basis_fem_RTk<T>;
see the Float page for the full documentation
see the point page for the full documentation
void grad_evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< tensor_basic< T >, Eigen::Dynamic, 1 > &value) const
void _compute_dofs(reference_element hat_K, const Eigen::Matrix< value_type, Eigen::Dynamic, 1 > &f_xnod, Eigen::Matrix< T, Eigen::Dynamic, 1 > &dof) const
void evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< value_type, Eigen::Dynamic, 1 > &value) const
const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > & vdm(reference_element hat_K) const
void _initialize_cstor_sizes() const
basis_raw_basic< T > _b_pre_kp1
basis_fem_RTk(size_type k, const basis_option &sopt)
const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > & hat_node(reference_element hat_K) const
std::string family_name() const
const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > & inv_vdm(reference_element hat_K) const
void _initialize_data(reference_element hat_K) const
see the basis_option page for the full documentation
reference_element::size_type size_type
piola_fem< T > _piola_fem
static std::string standard_naming(std::string family_name, size_t degree, const basis_option &sopt)
void set_family(family_type type)
rep::const_iterator const_iterator
see the reference_element page for the full documentation
void side_normal(size_type loc_isid, point_basic< Float > &hat_n) const
const point_basic< Float > & vertex(size_type iloc) const
static const variant_type H
size_type n_vertex() const
static const variant_type q
static const variant_type e
reference_element subgeo(size_type subgeo_dim, size_type loc_isid) const
size_type dimension() const
static const variant_type p
Float side_measure(size_type loc_isid) const
variant_type variant() const
size_type subgeo_local_vertex(size_type subgeo_dim, size_type loc_isid, size_type loc_jsidvert) const
static size_type n_node(variant_type variant, size_type order)
size_type n_subgeo(size_type subgeo_dim) const
static const variant_type T
static const variant_type P
static const variant_type t
integrate_option quadrature_option
#define trace_macro(message)
#define error_macro(message)
#define warning_macro(message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
Float alpha[pmax+1][pmax+1]
void basis_on_pointset_evaluate(const Basis &b, const reference_element &hat_K, const Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_x, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &vdm)
size_type ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
This file is part of Rheolef.
void pointset_lagrange_equispaced(reference_element hat_K, size_t order_in, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_xnod, size_t internal=0)
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
void pointset_lagrange_warburton(reference_element hat_K, size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_xnod, bool map_on_reference_element=true)
T cond(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &a)
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)