21 #include "libmesh/libmesh_config.h" 23 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 25 #include "libmesh/inf_fe.h" 26 #include "libmesh/fe.h" 27 #include "libmesh/quadrature_gauss.h" 28 #include "libmesh/elem.h" 29 #include "libmesh/libmesh_logging.h" 30 #include "libmesh/int_range.h" 31 #include "libmesh/type_tensor.h" 32 #include "libmesh/fe_interface.h" 42 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
46 calculate_map_scaled(false),
47 calculate_phi_scaled(false),
48 calculate_dphi_scaled(false),
51 _n_total_approx_sf (0),
61 current_fe_type (
FEType(fet.order,
79 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
86 const Order radial_int_order =
static_cast<Order>(2 * (
static_cast<unsigned int>(fe_type.radial_order.get_order()) + 1) +2);
87 const unsigned int qrule_dim = q->
get_dim();
93 base_fe->attach_quadrature_rule(base_qrule.get());
97 radial_qrule = std::make_unique<QGauss>(1, radial_int_order);
109 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_base>
120 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
122 const std::vector<Point> *
const pts,
123 const std::vector<Real> *
const weights)
130 this->determine_calculations();
135 libmesh_assert_equal_to (base_fe->qrule, base_qrule.get());
138 bool init_shape_functions_required =
false;
141 if (current_fe_type.radial_order != fe_type.radial_order)
143 current_fe_type.radial_order = fe_type.radial_order;
148 radial_qrule->init(
EDGE2);
151 this->init_radial_shape_functions(inf_elem);
153 init_shape_functions_required=
true;
157 bool update_base_elem_required=
true;
164 ((this->get_type() != inf_elem->
type()) ||
165 (base_fe->shapes_need_reinit())))
170 elem_type = inf_elem->
type();
171 this->update_base_elem(inf_elem);
172 update_base_elem_required=
false;
175 base_qrule->init(base_elem->type());
176 init_shape_functions_required=
true;
185 if (update_base_elem_required)
186 this->update_base_elem(inf_elem);
188 if (calculate_phi_scaled || calculate_dphi_scaled || calculate_phi || calculate_dphi)
190 base_fe->init_base_shape_functions(base_fe->qrule->get_points(),
195 base_fe->_fe_map->compute_map (base_fe->dim, base_fe->qrule->get_weights(),
196 base_elem.get(), base_fe->calculate_d2phi);
197 base_fe->compute_shape_functions(base_elem.get(), base_fe->qrule->get_points());
201 if (init_shape_functions_required)
202 this->init_shape_functions (radial_qrule->get_points(),
203 base_fe->qrule->get_points(),
208 this->compute_shape_functions (inf_elem,
209 base_fe->qrule->get_points(),
210 radial_qrule->get_points()
218 elem_type = inf_elem->
type();
225 std::vector<Point> radial_pts;
229 radial_pts.push_back(
radius);
230 unsigned int n_radial_pts=1;
231 unsigned int n_angular_pts=1;
234 radius = (*pts)[p](Dim-1);
240 if (p == (n_radial_pts)*n_angular_pts)
242 radial_pts.push_back(
radius);
247 libmesh_error_msg(
"We assumed that the "<<pts->size()
248 <<
" points are of tensor-product type with " 249 <<n_radial_pts<<
" radial points and " 250 <<n_angular_pts<<
" angular points."<<std::endl
251 <<
"But apparently point "<<p+1
252 <<
" does not fit that scheme: Its radius is " 253 <<
radius <<
"but should have " 254 <<radial_pts[n_radial_pts*n_angular_pts-p]<<
".");
260 else if (n_radial_pts == 1)
271 libmesh_error_msg(
"Calling reinit() with an empty point list is prohibited.\n");
274 const std::size_t radial_pts_size = radial_pts.size();
275 const std::size_t base_pts_size = pts->size() / radial_pts_size;
277 libmesh_assert_equal_to
278 (base_pts_size * radial_pts_size, pts->size());
281 std::vector<Point> base_pts;
282 base_pts.reserve(base_pts_size);
283 for (std::size_t p=0; p != base_pts_size; ++p)
285 Point pt = (*pts)[p];
287 base_pts.push_back(pt);
291 this->init_radial_shape_functions(inf_elem, &radial_pts);
294 this->update_base_elem(inf_elem);
300 this->determine_calculations();
302 base_fe->reinit( base_elem.get(), &base_pts);
304 this->init_shape_functions (radial_pts, base_pts, inf_elem);
307 if (weights !=
nullptr)
309 this->compute_shape_functions (inf_elem,base_pts,radial_pts);
313 this->compute_shape_functions (inf_elem, base_pts, radial_pts);
318 if (this->calculate_dual)
319 libmesh_not_implemented_msg(
"Dual shape support for infinite elements is " 320 "not currently implemented");
325 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
328 this->calculations_started =
true;
334 #ifdef LIBMESH_ENABLE_DEPRECATED 335 if (!this->calculate_nothing &&
336 !this->calculate_phi && !this->calculate_dphi &&
337 !this->calculate_dphiref &&
338 !this->calculate_phi_scaled && !this->calculate_dphi_scaled &&
339 !this->calculate_xyz && !this->calculate_jxw &&
340 !this->calculate_map_scaled && !this->calculate_map &&
341 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
342 !this->calculate_d2phi &&
344 !this->calculate_curl_phi && !this->calculate_div_phi)
346 libmesh_deprecated();
347 this->calculate_phi = this->calculate_dphi = this->calculate_jxw =
true;
348 this->calculate_dphiref =
true;
349 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 350 this->calculate_d2phi =
true;
352 this->calculate_phi_scaled = this->calculate_dphi_scaled = this->calculate_xyz =
true;
355 this->calculate_curl_phi =
true;
356 this->calculate_div_phi =
true;
359 #else //LIBMESH_ENABLE_DEPRECATED 362 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 363 libmesh_assert (this->calculate_nothing || this->calculate_d2phi ||
364 this->calculate_phi || this->calculate_dphi ||
365 this->calculate_dphiref ||
366 this->calculate_phi_scaled || this->calculate_dphi_scaled ||
367 this->calculate_xyz || this->calculate_jxw ||
368 this->calculate_map_scaled || this->calculate_map ||
369 this->calculate_curl_phi || this->calculate_div_phi);
372 this->calculate_phi || this->calculate_dphi ||
373 this->calculate_dphiref ||
374 this->calculate_phi_scaled || this->calculate_dphi_scaled ||
375 this->calculate_xyz || this->calculate_jxw ||
376 this->calculate_map_scaled || this->calculate_map ||
377 this->calculate_curl_phi || this->calculate_div_phi);
379 #endif // LIBMESH_ENABLE_DEPRECATED 383 this->calculate_map =
true;
384 if (this->calculate_dphi)
385 this->calculate_map =
true;
386 if (this->calculate_dphi_scaled)
387 this->calculate_map_scaled =
true;
388 if (calculate_xyz && !calculate_map)
391 this->calculate_map_scaled =
true;
392 base_fe->calculate_phi = this->calculate_phi || this->calculate_phi_scaled
393 || this->calculate_dphi || this->calculate_dphi_scaled;
394 base_fe->calculate_dphi = this->calculate_dphi || this->calculate_dphi_scaled;
395 if (this->calculate_map || this->calculate_map_scaled
396 || this->calculate_dphiref)
398 base_fe->calculate_dphiref =
true;
402 base_fe->determine_calculations();
404 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 405 if (this->calculate_d2phi)
406 libmesh_not_implemented();
407 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES 412 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
416 const std::vector<Point> * radial_pts)
422 LOG_SCOPE(
"init_radial_shape_functions()",
"InfFE");
425 const Order radial_approx_order = fe_type.radial_order;
426 const unsigned int n_radial_approx_shape_functions =
429 const std::size_t n_radial_qp =
430 radial_pts ? radial_pts->size() : radial_qrule->n_points();
431 const std::vector<Point> & radial_qp =
432 radial_pts ? *radial_pts : radial_qrule->get_points();
435 if (calculate_phi || calculate_dphi || calculate_phi_scaled || calculate_dphi_scaled)
437 mode.resize (n_radial_approx_shape_functions);
438 for (
unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
439 mode[i].resize (n_radial_qp);
442 for (
unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
443 for (std::size_t p=0; p<n_radial_qp; ++p)
447 if (calculate_dphi || calculate_dphi_scaled)
449 dmodedv.resize (n_radial_approx_shape_functions);
450 for (
unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
451 dmodedv[i].resize (n_radial_qp);
454 for (
unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
455 for (std::size_t p=0; p<n_radial_qp; ++p)
460 if (calculate_phi || calculate_phi_scaled || calculate_dphi || calculate_dphi_scaled)
462 som.resize (n_radial_qp);
464 for (std::size_t p=0; p<n_radial_qp; ++p)
467 if (calculate_dphi || calculate_dphi_scaled)
469 dsomdv.resize (n_radial_qp);
471 for (std::size_t p=0; p<n_radial_qp; ++p)
478 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
480 const std::vector<Point> & base_qp,
481 const Elem * inf_elem)
486 LOG_SCOPE(
"init_shape_functions()",
"InfFE");
490 const std::size_t n_radial_qp = radial_qp.size();
492 if (calculate_phi || calculate_dphi || calculate_phi_scaled || calculate_dphi_scaled)
493 libmesh_assert_equal_to(n_radial_approx_sf, mode.size());
494 if (calculate_phi || calculate_dphi || calculate_dphi_scaled)
495 libmesh_assert_equal_to(som.size(), n_radial_qp);
511 unsigned int n_base_approx_shape_functions;
513 n_base_approx_shape_functions = base_fe->n_shape_functions();
515 n_base_approx_shape_functions = 1;
520 n_radial_approx_sf * n_base_approx_shape_functions;
524 const unsigned int n_base_qp = cast_int<unsigned int>(base_qp.size());
527 _n_total_qp = n_radial_qp * n_base_qp;
534 _radial_shape_index.resize(_n_total_approx_sf);
535 _base_shape_index.resize(_n_total_approx_sf);
538 for (
unsigned int n=0; n<_n_total_approx_sf; ++n)
540 compute_shape_indices (this->fe_type,
543 _base_shape_index[n],
544 _radial_shape_index[n]);
545 libmesh_assert_less (_base_shape_index[n], n_base_approx_shape_functions);
546 libmesh_assert_less (_radial_shape_index[n], n_radial_approx_sf);
564 if (calculate_phi || calculate_dphi)
565 weight.resize (_n_total_qp);
566 if (calculate_phi_scaled || calculate_dphi_scaled)
567 weightxr_sq.resize (_n_total_qp);
568 if (calculate_dphi || calculate_dphi_scaled)
569 dweightdv.resize (n_radial_qp);
571 dweight.resize (_n_total_qp);
572 if (calculate_dphi_scaled)
573 dweightxr_sq.resize(_n_total_qp);
575 if (calculate_dphi || calculate_dphi_scaled)
576 dphase.resize (_n_total_qp);
580 _total_qrule_weights.resize(_n_total_qp,1);
585 if (calculate_map_scaled)
586 JxWxdecay.resize(_n_total_qp);
588 JxW.resize(_n_total_qp);
589 if (calculate_map_scaled || calculate_map)
591 xyz.resize(_n_total_qp);
592 dxidx_map_scaled.resize(_n_total_qp);
593 dxidy_map_scaled.resize(_n_total_qp);
594 dxidz_map_scaled.resize(_n_total_qp);
595 detadx_map_scaled.resize(_n_total_qp);
596 detady_map_scaled.resize(_n_total_qp);
597 detadz_map_scaled.resize(_n_total_qp);
598 dzetadx_map_scaled.resize(_n_total_qp);
599 dzetady_map_scaled.resize(_n_total_qp);
600 dzetadz_map_scaled.resize(_n_total_qp);
604 dxidx_map.resize(_n_total_qp);
605 dxidy_map.resize(_n_total_qp);
606 dxidz_map.resize(_n_total_qp);
607 detadx_map.resize(_n_total_qp);
608 detady_map.resize(_n_total_qp);
609 detadz_map.resize(_n_total_qp);
610 dzetadx_map.resize(_n_total_qp);
611 dzetady_map.resize(_n_total_qp);
612 dzetadz_map.resize(_n_total_qp);
615 phi.resize (_n_total_approx_sf);
616 if (calculate_phi_scaled)
617 phixr.resize (_n_total_approx_sf);
620 dphi.resize (_n_total_approx_sf);
621 dphidx.resize (_n_total_approx_sf);
622 dphidy.resize (_n_total_approx_sf);
623 dphidz.resize (_n_total_approx_sf);
626 if (calculate_dphi_scaled)
628 dphixr.resize (_n_total_approx_sf);
629 dphixr_sq.resize(_n_total_approx_sf);
631 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 635 libmesh_not_implemented();
636 d2phi.resize (_n_total_approx_sf);
637 d2phidx2.resize (_n_total_approx_sf);
638 d2phidxdy.resize (_n_total_approx_sf);
639 d2phidxdz.resize (_n_total_approx_sf);
640 d2phidy2.resize (_n_total_approx_sf);
641 d2phidydz.resize (_n_total_approx_sf);
642 d2phidz2.resize (_n_total_approx_sf);
643 d2phidxi2.resize (_n_total_approx_sf);
647 d2phidxideta.resize(_n_total_approx_sf);
648 d2phideta2.resize(_n_total_approx_sf);
653 d2phidetadzeta.resize(_n_total_approx_sf);
654 d2phidxidzeta.resize(_n_total_approx_sf);
655 d2phidzeta2.resize(_n_total_approx_sf);
658 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 660 if (calculate_dphi || calculate_dphi_scaled)
662 dphidxi.resize (_n_total_approx_sf);
665 dphideta.resize(_n_total_approx_sf);
668 dphidzeta.resize(_n_total_approx_sf);
677 for (
unsigned int i=0; i<_n_total_approx_sf; ++i)
678 phi[i].resize (_n_total_qp);
681 for (
unsigned int i=0; i<_n_total_approx_sf; ++i)
683 dphi[i].resize (_n_total_qp);
684 dphidx[i].resize (_n_total_qp);
685 dphidy[i].resize (_n_total_qp);
686 dphidz[i].resize (_n_total_qp);
689 if (calculate_phi_scaled)
690 for (
unsigned int i=0; i<_n_total_approx_sf; ++i)
692 phixr[i].resize (_n_total_qp);
694 if (calculate_dphi_scaled)
695 for (
unsigned int i=0; i<_n_total_approx_sf; ++i)
697 dphixr[i].resize(_n_total_qp);
698 dphixr_sq[i].resize(_n_total_qp);
700 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 702 for (
unsigned int i=0; i<_n_total_approx_sf; ++i)
704 d2phi[i].resize (_n_total_qp);
705 d2phidx2[i].resize (_n_total_qp);
706 d2phidxdy[i].resize (_n_total_qp);
707 d2phidxdz[i].resize (_n_total_qp);
708 d2phidy2[i].resize (_n_total_qp);
709 d2phidydz[i].resize (_n_total_qp);
710 d2phidy2[i].resize (_n_total_qp);
711 d2phidxi2[i].resize (_n_total_qp);
715 d2phidxideta[i].resize (_n_total_qp);
716 d2phideta2[i].resize (_n_total_qp);
720 d2phidxidzeta[i].resize (_n_total_qp);
721 d2phidetadzeta[i].resize (_n_total_qp);
722 d2phidzeta2[i].resize (_n_total_qp);
725 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 727 if (calculate_dphi || calculate_dphi_scaled)
728 for (
unsigned int i=0; i<_n_total_approx_sf; ++i)
730 dphidxi[i].resize (_n_total_qp);
733 dphideta[i].resize (_n_total_qp);
736 dphidzeta[i].resize (_n_total_qp);
746 libmesh_assert_equal_to (radial_qp.size(), n_radial_qp);
748 if (radial_qrule && base_qrule)
750 const std::vector<Real> & radial_qw = radial_qrule->get_weights();
751 const std::vector<Real> & base_qw = base_qrule->get_weights();
752 libmesh_assert_equal_to (radial_qw.size(), n_radial_qp);
753 libmesh_assert_equal_to (base_qw.size(), n_base_qp);
755 for (
unsigned int rp=0; rp<n_radial_qp; ++rp)
756 for (
unsigned int bp=0; bp<n_base_qp; ++bp)
757 _total_qrule_weights[bp + rp*n_base_qp] = radial_qw[rp] * base_qw[bp];
761 for (
unsigned int rp=0; rp<n_radial_qp; ++rp)
763 if (calculate_phi || calculate_dphi)
764 for (
unsigned int bp=0; bp<n_base_qp; ++bp)
767 if ( calculate_phi_scaled)
768 for (
unsigned int bp=0; bp<n_base_qp; ++bp)
771 if ( calculate_dphi || calculate_dphi_scaled)
779 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
781 const std::vector<Point> & base_qp,
782 const std::vector<Point> & radial_qp
788 libmesh_assert_equal_to (base_elem->type(),
792 LOG_SCOPE(
"compute_shape_functions()",
"InfFE");
796 const std::size_t n_radial_qp = radial_qp.size();
797 const unsigned int n_base_qp = base_qp.size();
799 libmesh_assert_equal_to (_n_total_qp, n_radial_qp*n_base_qp);
800 libmesh_assert_equal_to (_n_total_qp, _total_qrule_weights.size());
803 libmesh_assert_equal_to(n_radial_qp, som.size());
805 if (this->calculate_map || this->calculate_map_scaled)
809 const std::vector<std::vector<Real>> & S_map = (base_fe->get_fe_map()).get_phi_map();
810 if (S_map[0].size() > 0)
811 libmesh_assert_equal_to(n_base_qp, S_map[0].size());
814 libmesh_assert_equal_to(n_radial_qp, radial_qrule->n_points());
816 libmesh_assert_equal_to(n_base_qp, base_qrule->n_points());
817 libmesh_assert_equal_to(_n_total_qp % n_radial_qp, 0);
822 base_fe->n_shape_functions();
832 unsigned int elem_dim = inf_elem->
dim();
839 libmesh_not_implemented();
844 std::vector<std::vector<Real>> S (0);
845 std::vector<std::vector<Real>> Ss(0);
846 std::vector<std::vector<Real>> St(0);
848 std::vector<Real> base_dxidx (0);
849 std::vector<Real> base_dxidy (0);
850 std::vector<Real> base_dxidz (0);
851 std::vector<Real> base_detadx(0);
852 std::vector<Real> base_detady(0);
853 std::vector<Real> base_detadz(0);
855 std::vector<Point> base_xyz (0);
857 if (calculate_phi || calculate_phi_scaled ||
858 calculate_dphi || calculate_dphi_scaled)
862 if (calculate_map || calculate_map_scaled)
864 Ss = base_fe->dphidxi;
865 St = base_fe->dphideta;
867 base_dxidx = base_fe->get_dxidx();
868 base_dxidy = base_fe->get_dxidy();
869 base_dxidz = base_fe->get_dxidz();
870 base_detadx = base_fe->get_detadx();
871 base_detady = base_fe->get_detady();
872 base_detadz = base_fe->get_detadz();
874 base_xyz = base_fe->get_xyz();
877 ElemType base_type= base_elem->type();
881 libmesh_assert_equal_to (phi.size(), _n_total_approx_sf);
884 libmesh_assert_equal_to (dphidxi.size(), _n_total_approx_sf);
885 libmesh_assert_equal_to (dphideta.size(), _n_total_approx_sf);
886 libmesh_assert_equal_to (dphidzeta.size(), _n_total_approx_sf);
891 for (
unsigned int rp=0; rp<n_radial_qp; ++rp)
892 for (
unsigned int bp=0; bp<n_base_qp; ++bp)
902 if (calculate_map || calculate_map_scaled)
904 xyz[tp] =
InfFEMap::map(elem_dim, inf_elem,
Point(base_qp[bp](0),base_qp[bp](1),radial_qp[rp](0)));
906 const Point r(xyz[tp]-origin);
907 a=(base_xyz[bp]-origin).
norm();
913 libmesh_assert_less(
std::abs(som[rp] -a/r_norm) , 1e-7);
919 Point e_xi(base_dxidx[bp],
922 Point e_eta(base_detadx[bp],
929 grad_a_scaled=unit_r - normal/(normal*unit_r);
931 const Real dxi_er=base_dxidx[bp]* unit_r(0) + base_dxidy[bp] *unit_r(1) + base_dxidz[bp] *unit_r(2);
932 const Real deta_er=base_detadx[bp]*unit_r(0) + base_detady[bp]*unit_r(1) + base_detadz[bp]*unit_r(2);
936 if (!base_elem->has_affine_map())
946 const unsigned int n_sf = base_elem->n_nodes();
948 for (
unsigned int i=0; i< n_sf; ++i)
951 base_elem->default_order(),
952 i, 0, base_qp[bp]) * e_xi
954 base_elem->default_order(),
955 i, 1, base_qp[bp]) * e_eta);
957 tmp += (base_elem->node_ref(i) -origin).
norm()* dL_da_i;
961 grad_a_scaled = ( tmp - (tmp*unit_r)*unit_r ) / ( 1. - tmp*unit_r);
966 dxidx_map_scaled[tp] = (grad_a_scaled(0) - unit_r(0))*dxi_er +base_dxidx[bp];
967 dxidy_map_scaled[tp] = (grad_a_scaled(1) - unit_r(1))*dxi_er +base_dxidy[bp];
968 dxidz_map_scaled[tp] = (grad_a_scaled(2) - unit_r(2))*dxi_er +base_dxidz[bp];
971 detadx_map_scaled[tp] = (grad_a_scaled(0) - unit_r(0))*deta_er + base_detadx[bp];
972 detady_map_scaled[tp] = (grad_a_scaled(1) - unit_r(1))*deta_er + base_detady[bp];
973 detadz_map_scaled[tp] = (grad_a_scaled(2) - unit_r(2))*deta_er + base_detadz[bp];
976 dzetadx_map_scaled[tp] =-2./a*(grad_a_scaled(0) - unit_r(0));
977 dzetady_map_scaled[tp] =-2./a*(grad_a_scaled(1) - unit_r(1));
978 dzetadz_map_scaled[tp] =-2./a*(grad_a_scaled(2) - unit_r(2));
984 dxidx_map[tp] = a/r_norm * dxidx_map_scaled[tp];
985 dxidy_map[tp] = a/r_norm * dxidy_map_scaled[tp];
986 dxidz_map[tp] = a/r_norm * dxidz_map_scaled[tp];
988 detadx_map[tp] = a/r_norm * detadx_map_scaled[tp];
989 detady_map[tp] = a/r_norm * detady_map_scaled[tp];
990 detadz_map[tp] = a/r_norm * detadz_map_scaled[tp];
994 dzetadx_map[tp] =-2.*a/(r_norm*r_norm)*(grad_a_scaled(0) - unit_r(0));
995 dzetady_map[tp] =-2.*a/(r_norm*r_norm)*(grad_a_scaled(1) - unit_r(1));
996 dzetadz_map[tp] =-2.*a/(r_norm*r_norm)*(grad_a_scaled(2) - unit_r(2));
1000 Real inv_jac = (dxidx_map[tp]*( detady_map[tp]*dzetadz_map[tp]- dzetady_map[tp]*detadz_map[tp]) +
1001 detadx_map[tp]*(dzetady_map[tp]* dxidz_map[tp]- dxidy_map[tp]*dzetadz_map[tp]) +
1002 dzetadx_map[tp]*( dxidy_map[tp]*detadz_map[tp]- detady_map[tp]* dxidz_map[tp]));
1004 if (inv_jac <= 1e-10)
1006 libmesh_error_msg(
"ERROR: negative inverse Jacobian " \
1015 JxW[tp] = _total_qrule_weights[tp]/inv_jac;
1019 if (calculate_map_scaled)
1021 Real inv_jacxR_pow4 = (dxidx_map_scaled[tp] *( detady_map_scaled[tp]*dzetadz_map_scaled[tp]
1022 - dzetady_map_scaled[tp]* detadz_map_scaled[tp]) +
1023 detadx_map_scaled[tp] *( dzetady_map_scaled[tp]* dxidz_map_scaled[tp]
1024 - dxidy_map_scaled[tp]*dzetadz_map_scaled[tp]) +
1025 dzetadx_map_scaled[tp]*( dxidy_map_scaled[tp]* detadz_map_scaled[tp]
1026 -detady_map_scaled[tp]* dxidz_map_scaled[tp]));
1027 if (inv_jacxR_pow4 <= 1e-7)
1029 libmesh_error_msg(
"ERROR: negative weighted inverse Jacobian " \
1037 JxWxdecay[tp] = _total_qrule_weights[tp]/inv_jacxR_pow4;
1043 if (calculate_dphi || calculate_dphi_scaled)
1044 dphase[tp] = unit_r - grad_a_scaled*a/r_norm;
1048 dweight[tp](0) = dweightdv[rp] * dzetadx_map[tp];
1049 dweight[tp](1) = dweightdv[rp] * dzetady_map[tp];
1050 dweight[tp](2) = dweightdv[rp] * dzetadz_map[tp];
1052 if (calculate_dphi_scaled)
1054 dweightxr_sq[tp](0) = dweightdv[rp] * dzetadx_map_scaled[tp];
1055 dweightxr_sq[tp](1) = dweightdv[rp] * dzetady_map_scaled[tp];
1056 dweightxr_sq[tp](2) = dweightdv[rp] * dzetadz_map_scaled[tp];
1059 if (calculate_phi || calculate_phi_scaled || calculate_dphi || calculate_dphi_scaled)
1061 for (
unsigned int i=0; i <_n_total_approx_sf ; ++i)
1064 unsigned int bi = _base_shape_index [i];
1065 unsigned int ri = _radial_shape_index[i];
1067 phi [i][tp] = S [bi][bp] * mode[ri][rp] * som[rp];
1069 if (calculate_phi_scaled)
1070 phixr [i][tp] = S [bi][bp] * mode[ri][rp];
1072 if (calculate_dphi || calculate_dphi_scaled)
1074 dphidxi [i][tp] = Ss[bi][bp] * mode[ri][rp] * som[rp];
1075 dphideta [i][tp] = St[bi][bp] * mode[ri][rp] * som[rp];
1076 dphidzeta[i][tp] = S [bi][bp]
1077 * (dmodedv[ri][rp] * som[rp] + mode[ri][rp] * dsomdv[rp]);
1085 dphidx[i][tp] = (dphidxi[i][tp]*dxidx_map[tp] +
1086 dphideta[i][tp]*detadx_map[tp] +
1087 dphidzeta[i][tp]*dzetadx_map[tp]);
1091 dphidy[i][tp] = (dphidxi[i][tp]*dxidy_map[tp] +
1092 dphideta[i][tp]*detady_map[tp] +
1093 dphidzeta[i][tp]*dzetady_map[tp]);
1097 dphidz[i][tp] = (dphidxi[i][tp]*dxidz_map[tp] +
1098 dphideta[i][tp]*detadz_map[tp] +
1099 dphidzeta[i][tp]*dzetadz_map[tp]);
1102 if (calculate_dphi_scaled)
1105 dphixr[i][tp](0)= (dphidxi[i][tp]*dxidx_map_scaled[tp] +
1106 dphideta[i][tp]*detadx_map_scaled[tp] +
1107 dphidzeta[i][tp]*dzetadx_map_scaled[tp]*som[rp]);
1109 dphixr[i][tp](1) = (dphidxi[i][tp]*dxidy_map_scaled[tp] +
1110 dphideta[i][tp]*detady_map_scaled[tp] +
1111 dphidzeta[i][tp]*dzetady_map_scaled[tp]*som[rp]);
1113 dphixr[i][tp](2) = (dphidxi[i][tp]*dxidz_map_scaled[tp] +
1114 dphideta[i][tp]*detadz_map_scaled[tp] +
1115 dphidzeta[i][tp]*dzetadz_map_scaled[tp]*som[rp]);
1117 const Real dphidxixr = Ss[bi][bp] * mode[ri][rp];
1118 const Real dphidetaxr= St[bi][bp] * mode[ri][rp];
1120 dphixr_sq[i][tp](0)= (dphidxixr*dxidx_map_scaled[tp] +
1121 dphidetaxr*detadx_map_scaled[tp] +
1122 dphidzeta[i][tp]*dzetadx_map_scaled[tp]);
1124 dphixr_sq[i][tp](1) = (dphidxixr*dxidy_map_scaled[tp] +
1125 dphidetaxr*detady_map_scaled[tp] +
1126 dphidzeta[i][tp]*dzetady_map_scaled[tp]);
1128 dphixr_sq[i][tp](2) = (dphidxixr*dxidz_map_scaled[tp] +
1129 dphidetaxr*detadz_map_scaled[tp] +
1130 dphidzeta[i][tp]*dzetadz_map_scaled[tp]);
1141 libmesh_error_msg(
"Unsupported dim = " <<
dim);
1147 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
1151 libmesh_not_implemented();
1159 #include "libmesh/inf_fe_instantiate_1D.h" 1160 #include "libmesh/inf_fe_instantiate_2D.h" 1161 #include "libmesh/inf_fe_instantiate_3D.h" 1165 #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
virtual void determine_calculations() override
Determine which values are to be calculated, for both the FE itself and for the FEMap.
ElemType
Defines an enum for geometric element types.
static ElemType get_elem_type(const ElemType type)
Order
defines an enum for polynomial orders.
A specific instantiation of the FEBase class.
static unsigned int n_dofs(const Order o_radial)
auto norm() const -> decltype(std::norm(T()))
virtual Point origin() const
static Point map(const unsigned int dim, const Elem *inf_elem, const Point &reference_point)
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
static Real decay_deriv(const unsigned int dim, const Real)
static Real Dxr_sq(const Real)
void init_shape_functions(const std::vector< Point > &radial_qp, const std::vector< Point > &base_qp, const Elem *inf_elem)
Initialize all the data fields like weight, mode, phi, dphidxi, dphideta, dphidzeta, etc.
This is the base class from which all geometric element types are derived.
static FEFieldType field_type(const FEType &fe_type)
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The libMesh namespace provides an interface to certain functionality in the library.
static std::unique_ptr< const Elem > build_elem(const Elem *inf_elem)
Build the base element of an infinite element.
virtual QuadratureType type() const =0
static Real D(const Real v)
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
A specific instantiation of the FEBase class.
friend class InfFE
Make all InfFE<Dim,T_radial,T_map> classes friends of each other, so that the protected eval() may be...
std::unique_ptr< FEBase > base_fe
Have a FE<Dim-1,T_base> handy for base approximation.
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
static Real decay(const unsigned int dim, const Real v)
unsigned int get_dim() const
InfMapType inf_map
The coordinate mapping type of the infinite element.
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &v) const
virtual bool shapes_need_reinit() const override
void init_radial_shape_functions(const Elem *inf_elem, const std::vector< Point > *radial_pts=nullptr)
Some of the member data only depend on the radial part of the infinite element.
FEFamily radial_family
The type of approximation in radial direction.
void compute_shape_functions(const Elem *inf_elem, const std::vector< Point > &base_qp, const std::vector< Point > &radial_qp)
After having updated the jacobian and the transformation from local to global coordinates in FEAbstra...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned short dim() const =0
static Real D_deriv(const Real v)
static std::unique_ptr< QBase > build(std::string_view name, const unsigned int dim, const Order order=INVALID_ORDER)
Builds a specific quadrature rule based on the name string.
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr) override
This is at the core of this class.
virtual void attach_quadrature_rule(QBase *q) override
The use of quadrature rules with the InfFE class is somewhat different from the approach of the FE cl...
FEType fe_type
The finite element type for this object.
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
The QBase class provides the basic functionality from which various quadrature rules can be derived...
This class forms the foundation from which generic finite elements may be derived.
void update_base_elem(const Elem *inf_elem)
Updates the protected member base_elem to the appropriate base element for the given inf_elem...