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),
60 current_fe_type (
FEType(fet.order,
78 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
85 const Order radial_int_order =
static_cast<Order>(2 * (
static_cast<unsigned int>(fe_type.radial_order.get_order()) + 1) +2);
86 const unsigned int qrule_dim = q->
get_dim();
92 base_fe->attach_quadrature_rule(base_qrule.get());
96 radial_qrule = std::make_unique<QGauss>(1, radial_int_order);
108 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_base>
119 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
121 const std::vector<Point> *
const pts,
122 const std::vector<Real> *
const weights)
129 this->determine_calculations();
134 libmesh_assert_equal_to (base_fe->qrule, base_qrule.get());
137 bool init_shape_functions_required =
false;
140 if (current_fe_type.radial_order != fe_type.radial_order)
142 current_fe_type.radial_order = fe_type.radial_order;
147 radial_qrule->init(
EDGE2);
150 this->init_radial_shape_functions(inf_elem);
152 init_shape_functions_required=
true;
156 bool update_base_elem_required=
true;
163 ((this->get_type() != inf_elem->
type()) ||
164 (base_fe->shapes_need_reinit())))
169 this->_elem = inf_elem;
170 this->_elem_type = inf_elem->
type();
171 this->update_base_elem(inf_elem);
172 update_base_elem_required=
false;
175 base_qrule->init(*base_elem);
176 init_shape_functions_required=
true;
180 this->_elem = inf_elem;
187 if (update_base_elem_required)
188 this->update_base_elem(inf_elem);
190 if (calculate_phi_scaled || calculate_dphi_scaled || calculate_phi || calculate_dphi)
192 base_fe->init_base_shape_functions(base_fe->qrule->get_points(),
197 base_fe->_fe_map->compute_map (base_fe->dim, base_fe->qrule->get_weights(),
198 base_elem.get(), base_fe->calculate_d2phi);
199 base_fe->compute_shape_functions(base_elem.get(), base_fe->qrule->get_points());
203 if (init_shape_functions_required)
204 this->init_shape_functions (radial_qrule->get_points(),
205 base_fe->qrule->get_points(),
210 this->compute_shape_functions (inf_elem,
211 base_fe->qrule->get_points(),
212 radial_qrule->get_points()
220 this->_elem = inf_elem;
221 this->_elem_type = inf_elem->
type();
228 std::vector<Point> radial_pts;
232 radial_pts.push_back(
radius);
233 unsigned int n_radial_pts=1;
234 unsigned int n_angular_pts=1;
237 radius = (*pts)[p](Dim-1);
240 if (std::abs(radial_pts[n_radial_pts-1](0) -
radius) > 1e-4)
243 if (p == (n_radial_pts)*n_angular_pts)
245 radial_pts.push_back(
radius);
250 libmesh_error_msg(
"We assumed that the "<<pts->size()
251 <<
" points are of tensor-product type with " 252 <<n_radial_pts<<
" radial points and " 253 <<n_angular_pts<<
" angular points."<<std::endl
254 <<
"But apparently point "<<p+1
255 <<
" does not fit that scheme: Its radius is " 256 <<
radius <<
"but should have " 257 <<radial_pts[n_radial_pts*n_angular_pts-p]<<
".");
263 else if (n_radial_pts == 1)
274 libmesh_error_msg(
"Calling reinit() with an empty point list is prohibited.\n");
277 const std::size_t radial_pts_size = radial_pts.size();
278 const std::size_t base_pts_size = pts->size() / radial_pts_size;
280 libmesh_assert_equal_to
281 (base_pts_size * radial_pts_size, pts->size());
284 std::vector<Point> base_pts;
285 base_pts.reserve(base_pts_size);
286 for (std::size_t p=0; p != base_pts_size; ++p)
288 Point pt = (*pts)[p];
290 base_pts.push_back(pt);
294 this->init_radial_shape_functions(inf_elem, &radial_pts);
297 this->update_base_elem(inf_elem);
303 this->determine_calculations();
305 base_fe->reinit( base_elem.get(), &base_pts);
307 this->init_shape_functions (radial_pts, base_pts, inf_elem);
310 if (weights !=
nullptr)
312 this->compute_shape_functions (inf_elem,base_pts,radial_pts);
316 this->compute_shape_functions (inf_elem, base_pts, radial_pts);
321 if (this->calculate_dual)
322 libmesh_not_implemented_msg(
"Dual shape support for infinite elements is " 323 "not currently implemented");
328 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
331 this->calculations_started =
true;
337 #ifdef LIBMESH_ENABLE_DEPRECATED 338 if (!this->calculate_nothing &&
339 !this->calculate_phi && !this->calculate_dphi &&
340 !this->calculate_dphiref &&
341 !this->calculate_phi_scaled && !this->calculate_dphi_scaled &&
342 !this->calculate_xyz && !this->calculate_jxw &&
343 !this->calculate_map_scaled && !this->calculate_map &&
344 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
345 !this->calculate_d2phi &&
347 !this->calculate_curl_phi && !this->calculate_div_phi)
349 libmesh_deprecated();
350 this->calculate_phi = this->calculate_dphi = this->calculate_jxw =
true;
351 this->calculate_dphiref =
true;
352 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 353 this->calculate_d2phi =
true;
355 this->calculate_phi_scaled = this->calculate_dphi_scaled = this->calculate_xyz =
true;
358 this->calculate_curl_phi =
true;
359 this->calculate_div_phi =
true;
362 #else //LIBMESH_ENABLE_DEPRECATED 365 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 366 libmesh_assert (this->calculate_nothing || this->calculate_d2phi ||
367 this->calculate_phi || this->calculate_dphi ||
368 this->calculate_dphiref ||
369 this->calculate_phi_scaled || this->calculate_dphi_scaled ||
370 this->calculate_xyz || this->calculate_jxw ||
371 this->calculate_map_scaled || this->calculate_map ||
372 this->calculate_curl_phi || this->calculate_div_phi);
375 this->calculate_phi || this->calculate_dphi ||
376 this->calculate_dphiref ||
377 this->calculate_phi_scaled || this->calculate_dphi_scaled ||
378 this->calculate_xyz || this->calculate_jxw ||
379 this->calculate_map_scaled || this->calculate_map ||
380 this->calculate_curl_phi || this->calculate_div_phi);
382 #endif // LIBMESH_ENABLE_DEPRECATED 386 this->calculate_map =
true;
387 if (this->calculate_dphi)
388 this->calculate_map =
true;
389 if (this->calculate_dphi_scaled)
390 this->calculate_map_scaled =
true;
391 if (calculate_xyz && !calculate_map)
394 this->calculate_map_scaled =
true;
395 base_fe->calculate_phi = this->calculate_phi || this->calculate_phi_scaled
396 || this->calculate_dphi || this->calculate_dphi_scaled;
397 base_fe->calculate_dphi = this->calculate_dphi || this->calculate_dphi_scaled;
398 if (this->calculate_map || this->calculate_map_scaled
399 || this->calculate_dphiref)
401 base_fe->calculate_dphiref =
true;
405 base_fe->determine_calculations();
407 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 408 if (this->calculate_d2phi)
409 libmesh_not_implemented();
410 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES 415 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
419 const std::vector<Point> * radial_pts)
425 LOG_SCOPE(
"init_radial_shape_functions()",
"InfFE");
428 const Order radial_approx_order = fe_type.radial_order;
429 const unsigned int n_radial_approx_shape_functions =
432 const std::size_t n_radial_qp =
433 radial_pts ? radial_pts->size() : radial_qrule->n_points();
434 const std::vector<Point> & radial_qp =
435 radial_pts ? *radial_pts : radial_qrule->get_points();
438 if (calculate_phi || calculate_dphi || calculate_phi_scaled || calculate_dphi_scaled)
440 mode.resize (n_radial_approx_shape_functions);
441 for (
unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
442 mode[i].resize (n_radial_qp);
445 for (
unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
446 for (std::size_t p=0; p<n_radial_qp; ++p)
450 if (calculate_dphi || calculate_dphi_scaled)
452 dmodedv.resize (n_radial_approx_shape_functions);
453 for (
unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
454 dmodedv[i].resize (n_radial_qp);
457 for (
unsigned int i=0; i<n_radial_approx_shape_functions; ++i)
458 for (std::size_t p=0; p<n_radial_qp; ++p)
463 if (calculate_phi || calculate_phi_scaled || calculate_dphi || calculate_dphi_scaled)
465 som.resize (n_radial_qp);
467 for (std::size_t p=0; p<n_radial_qp; ++p)
470 if (calculate_dphi || calculate_dphi_scaled)
472 dsomdv.resize (n_radial_qp);
474 for (std::size_t p=0; p<n_radial_qp; ++p)
481 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
483 const std::vector<Point> & base_qp,
484 const Elem * inf_elem)
489 LOG_SCOPE(
"init_shape_functions()",
"InfFE");
493 const std::size_t n_radial_qp = radial_qp.size();
495 if (calculate_phi || calculate_dphi || calculate_phi_scaled || calculate_dphi_scaled)
496 libmesh_assert_equal_to(n_radial_approx_sf, mode.size());
497 if (calculate_phi || calculate_dphi || calculate_dphi_scaled)
498 libmesh_assert_equal_to(som.size(), n_radial_qp);
514 unsigned int n_base_approx_shape_functions;
516 n_base_approx_shape_functions =
519 n_base_approx_shape_functions = 1;
524 n_radial_approx_sf * n_base_approx_shape_functions;
528 const unsigned int n_base_qp = cast_int<unsigned int>(base_qp.size());
531 _n_total_qp = n_radial_qp * n_base_qp;
538 _radial_shape_index.resize(_n_total_approx_sf);
539 _base_shape_index.resize(_n_total_approx_sf);
542 for (
unsigned int n=0; n<_n_total_approx_sf; ++n)
544 compute_shape_indices (this->fe_type,
547 _base_shape_index[n],
548 _radial_shape_index[n]);
549 libmesh_assert_less (_base_shape_index[n], n_base_approx_shape_functions);
550 libmesh_assert_less (_radial_shape_index[n], n_radial_approx_sf);
568 if (calculate_phi || calculate_dphi)
569 weight.resize (_n_total_qp);
570 if (calculate_phi_scaled || calculate_dphi_scaled)
571 weightxr_sq.resize (_n_total_qp);
572 if (calculate_dphi || calculate_dphi_scaled)
573 dweightdv.resize (n_radial_qp);
575 dweight.resize (_n_total_qp);
576 if (calculate_dphi_scaled)
577 dweightxr_sq.resize(_n_total_qp);
579 if (calculate_dphi || calculate_dphi_scaled)
580 dphase.resize (_n_total_qp);
584 _total_qrule_weights.resize(_n_total_qp,1);
589 if (calculate_map_scaled)
590 JxWxdecay.resize(_n_total_qp);
592 JxW.resize(_n_total_qp);
593 if (calculate_map_scaled || calculate_map)
595 xyz.resize(_n_total_qp);
596 dxidx_map_scaled.resize(_n_total_qp);
597 dxidy_map_scaled.resize(_n_total_qp);
598 dxidz_map_scaled.resize(_n_total_qp);
599 detadx_map_scaled.resize(_n_total_qp);
600 detady_map_scaled.resize(_n_total_qp);
601 detadz_map_scaled.resize(_n_total_qp);
602 dzetadx_map_scaled.resize(_n_total_qp);
603 dzetady_map_scaled.resize(_n_total_qp);
604 dzetadz_map_scaled.resize(_n_total_qp);
608 dxidx_map.resize(_n_total_qp);
609 dxidy_map.resize(_n_total_qp);
610 dxidz_map.resize(_n_total_qp);
611 detadx_map.resize(_n_total_qp);
612 detady_map.resize(_n_total_qp);
613 detadz_map.resize(_n_total_qp);
614 dzetadx_map.resize(_n_total_qp);
615 dzetady_map.resize(_n_total_qp);
616 dzetadz_map.resize(_n_total_qp);
619 phi.resize (_n_total_approx_sf);
620 if (calculate_phi_scaled)
621 phixr.resize (_n_total_approx_sf);
624 dphi.resize (_n_total_approx_sf);
625 dphidx.resize (_n_total_approx_sf);
626 dphidy.resize (_n_total_approx_sf);
627 dphidz.resize (_n_total_approx_sf);
630 if (calculate_dphi_scaled)
632 dphixr.resize (_n_total_approx_sf);
633 dphixr_sq.resize(_n_total_approx_sf);
635 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 639 libmesh_not_implemented();
640 d2phi.resize (_n_total_approx_sf);
641 d2phidx2.resize (_n_total_approx_sf);
642 d2phidxdy.resize (_n_total_approx_sf);
643 d2phidxdz.resize (_n_total_approx_sf);
644 d2phidy2.resize (_n_total_approx_sf);
645 d2phidydz.resize (_n_total_approx_sf);
646 d2phidz2.resize (_n_total_approx_sf);
647 d2phidxi2.resize (_n_total_approx_sf);
651 d2phidxideta.resize(_n_total_approx_sf);
652 d2phideta2.resize(_n_total_approx_sf);
657 d2phidetadzeta.resize(_n_total_approx_sf);
658 d2phidxidzeta.resize(_n_total_approx_sf);
659 d2phidzeta2.resize(_n_total_approx_sf);
662 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 664 if (calculate_dphi || calculate_dphi_scaled)
666 dphidxi.resize (_n_total_approx_sf);
669 dphideta.resize(_n_total_approx_sf);
672 dphidzeta.resize(_n_total_approx_sf);
681 for (
unsigned int i=0; i<_n_total_approx_sf; ++i)
682 phi[i].resize (_n_total_qp);
685 for (
unsigned int i=0; i<_n_total_approx_sf; ++i)
687 dphi[i].resize (_n_total_qp);
688 dphidx[i].resize (_n_total_qp);
689 dphidy[i].resize (_n_total_qp);
690 dphidz[i].resize (_n_total_qp);
693 if (calculate_phi_scaled)
694 for (
unsigned int i=0; i<_n_total_approx_sf; ++i)
696 phixr[i].resize (_n_total_qp);
698 if (calculate_dphi_scaled)
699 for (
unsigned int i=0; i<_n_total_approx_sf; ++i)
701 dphixr[i].resize(_n_total_qp);
702 dphixr_sq[i].resize(_n_total_qp);
704 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 706 for (
unsigned int i=0; i<_n_total_approx_sf; ++i)
708 d2phi[i].resize (_n_total_qp);
709 d2phidx2[i].resize (_n_total_qp);
710 d2phidxdy[i].resize (_n_total_qp);
711 d2phidxdz[i].resize (_n_total_qp);
712 d2phidy2[i].resize (_n_total_qp);
713 d2phidydz[i].resize (_n_total_qp);
714 d2phidy2[i].resize (_n_total_qp);
715 d2phidxi2[i].resize (_n_total_qp);
719 d2phidxideta[i].resize (_n_total_qp);
720 d2phideta2[i].resize (_n_total_qp);
724 d2phidxidzeta[i].resize (_n_total_qp);
725 d2phidetadzeta[i].resize (_n_total_qp);
726 d2phidzeta2[i].resize (_n_total_qp);
729 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 731 if (calculate_dphi || calculate_dphi_scaled)
732 for (
unsigned int i=0; i<_n_total_approx_sf; ++i)
734 dphidxi[i].resize (_n_total_qp);
737 dphideta[i].resize (_n_total_qp);
740 dphidzeta[i].resize (_n_total_qp);
750 libmesh_assert_equal_to (radial_qp.size(), n_radial_qp);
752 if (radial_qrule && base_qrule)
754 const std::vector<Real> & radial_qw = radial_qrule->get_weights();
755 const std::vector<Real> & base_qw = base_qrule->get_weights();
756 libmesh_assert_equal_to (radial_qw.size(), n_radial_qp);
757 libmesh_assert_equal_to (base_qw.size(), n_base_qp);
759 for (
unsigned int rp=0; rp<n_radial_qp; ++rp)
760 for (
unsigned int bp=0; bp<n_base_qp; ++bp)
761 _total_qrule_weights[bp + rp*n_base_qp] = radial_qw[rp] * base_qw[bp];
765 for (
unsigned int rp=0; rp<n_radial_qp; ++rp)
767 if (calculate_phi || calculate_dphi)
768 for (
unsigned int bp=0; bp<n_base_qp; ++bp)
771 if ( calculate_phi_scaled)
772 for (
unsigned int bp=0; bp<n_base_qp; ++bp)
775 if ( calculate_dphi || calculate_dphi_scaled)
783 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
785 const std::vector<Point> & base_qp,
786 const std::vector<Point> & radial_qp
792 libmesh_assert_equal_to (base_elem->type(),
796 LOG_SCOPE(
"compute_shape_functions()",
"InfFE");
800 const std::size_t n_radial_qp = radial_qp.size();
801 const unsigned int n_base_qp = base_qp.size();
803 libmesh_assert_equal_to (_n_total_qp, n_radial_qp*n_base_qp);
804 libmesh_assert_equal_to (_n_total_qp, _total_qrule_weights.size());
807 libmesh_assert_equal_to(n_radial_qp, som.size());
809 if (this->calculate_map || this->calculate_map_scaled)
813 const std::vector<std::vector<Real>> & S_map = (base_fe->get_fe_map()).get_phi_map();
814 if (S_map[0].size() > 0)
815 libmesh_assert_equal_to(n_base_qp, S_map[0].size());
818 libmesh_assert_equal_to(n_radial_qp, radial_qrule->n_points());
820 libmesh_assert_equal_to(n_base_qp, base_qrule->n_points());
821 libmesh_assert_equal_to(_n_total_qp % n_radial_qp, 0);
836 unsigned int elem_dim = inf_elem->
dim();
843 libmesh_not_implemented();
848 std::vector<std::vector<Real>> S (0);
849 std::vector<std::vector<Real>> Ss(0);
850 std::vector<std::vector<Real>> St(0);
852 std::vector<Real> base_dxidx (0);
853 std::vector<Real> base_dxidy (0);
854 std::vector<Real> base_dxidz (0);
855 std::vector<Real> base_detadx(0);
856 std::vector<Real> base_detady(0);
857 std::vector<Real> base_detadz(0);
859 std::vector<Point> base_xyz (0);
861 if (calculate_phi || calculate_phi_scaled ||
862 calculate_dphi || calculate_dphi_scaled)
866 if (calculate_map || calculate_map_scaled)
868 Ss = base_fe->dphidxi;
869 St = base_fe->dphideta;
871 base_dxidx = base_fe->get_dxidx();
872 base_dxidy = base_fe->get_dxidy();
873 base_dxidz = base_fe->get_dxidz();
874 base_detadx = base_fe->get_detadx();
875 base_detady = base_fe->get_detady();
876 base_detadz = base_fe->get_detadz();
878 base_xyz = base_fe->get_xyz();
881 ElemType base_type= base_elem->type();
885 libmesh_assert_equal_to (phi.size(), _n_total_approx_sf);
888 libmesh_assert_equal_to (dphidxi.size(), _n_total_approx_sf);
889 libmesh_assert_equal_to (dphideta.size(), _n_total_approx_sf);
890 libmesh_assert_equal_to (dphidzeta.size(), _n_total_approx_sf);
895 for (
unsigned int rp=0; rp<n_radial_qp; ++rp)
896 for (
unsigned int bp=0; bp<n_base_qp; ++bp)
906 if (calculate_map || calculate_map_scaled)
908 xyz[tp] =
InfFEMap::map(elem_dim, inf_elem,
Point(base_qp[bp](0),base_qp[bp](1),radial_qp[rp](0)));
910 const Point r(xyz[tp]-origin);
911 a=(base_xyz[bp]-origin).
norm();
917 libmesh_assert_less(std::abs(som[rp] -a/r_norm) , 1e-7);
923 Point e_xi(base_dxidx[bp],
926 Point e_eta(base_detadx[bp],
933 grad_a_scaled=unit_r - normal/(normal*unit_r);
935 const Real dxi_er=base_dxidx[bp]* unit_r(0) + base_dxidy[bp] *unit_r(1) + base_dxidz[bp] *unit_r(2);
936 const Real deta_er=base_detadx[bp]*unit_r(0) + base_detady[bp]*unit_r(1) + base_detadz[bp]*unit_r(2);
940 if (!base_elem->has_affine_map())
950 const unsigned int n_sf = base_elem->n_nodes();
952 for (
unsigned int i=0; i< n_sf; ++i)
955 base_elem->default_order(),
956 i, 0, base_qp[bp]) * e_xi
958 base_elem->default_order(),
959 i, 1, base_qp[bp]) * e_eta);
961 tmp += (base_elem->node_ref(i) -origin).
norm()* dL_da_i;
965 grad_a_scaled = ( tmp - (tmp*unit_r)*unit_r ) / ( 1. - tmp*unit_r);
970 dxidx_map_scaled[tp] = (grad_a_scaled(0) - unit_r(0))*dxi_er +base_dxidx[bp];
971 dxidy_map_scaled[tp] = (grad_a_scaled(1) - unit_r(1))*dxi_er +base_dxidy[bp];
972 dxidz_map_scaled[tp] = (grad_a_scaled(2) - unit_r(2))*dxi_er +base_dxidz[bp];
975 detadx_map_scaled[tp] = (grad_a_scaled(0) - unit_r(0))*deta_er + base_detadx[bp];
976 detady_map_scaled[tp] = (grad_a_scaled(1) - unit_r(1))*deta_er + base_detady[bp];
977 detadz_map_scaled[tp] = (grad_a_scaled(2) - unit_r(2))*deta_er + base_detadz[bp];
980 dzetadx_map_scaled[tp] =-2./a*(grad_a_scaled(0) - unit_r(0));
981 dzetady_map_scaled[tp] =-2./a*(grad_a_scaled(1) - unit_r(1));
982 dzetadz_map_scaled[tp] =-2./a*(grad_a_scaled(2) - unit_r(2));
988 dxidx_map[tp] = a/r_norm * dxidx_map_scaled[tp];
989 dxidy_map[tp] = a/r_norm * dxidy_map_scaled[tp];
990 dxidz_map[tp] = a/r_norm * dxidz_map_scaled[tp];
992 detadx_map[tp] = a/r_norm * detadx_map_scaled[tp];
993 detady_map[tp] = a/r_norm * detady_map_scaled[tp];
994 detadz_map[tp] = a/r_norm * detadz_map_scaled[tp];
998 dzetadx_map[tp] =-2.*a/(r_norm*r_norm)*(grad_a_scaled(0) - unit_r(0));
999 dzetady_map[tp] =-2.*a/(r_norm*r_norm)*(grad_a_scaled(1) - unit_r(1));
1000 dzetadz_map[tp] =-2.*a/(r_norm*r_norm)*(grad_a_scaled(2) - unit_r(2));
1004 Real inv_jac = (dxidx_map[tp]*( detady_map[tp]*dzetadz_map[tp]- dzetady_map[tp]*detadz_map[tp]) +
1005 detadx_map[tp]*(dzetady_map[tp]* dxidz_map[tp]- dxidy_map[tp]*dzetadz_map[tp]) +
1006 dzetadx_map[tp]*( dxidy_map[tp]*detadz_map[tp]- detady_map[tp]* dxidz_map[tp]));
1008 if (inv_jac <= 1e-10)
1010 libmesh_error_msg(
"ERROR: negative inverse Jacobian " \
1019 JxW[tp] = _total_qrule_weights[tp]/inv_jac;
1023 if (calculate_map_scaled)
1025 Real inv_jacxR_pow4 = (dxidx_map_scaled[tp] *( detady_map_scaled[tp]*dzetadz_map_scaled[tp]
1026 - dzetady_map_scaled[tp]* detadz_map_scaled[tp]) +
1027 detadx_map_scaled[tp] *( dzetady_map_scaled[tp]* dxidz_map_scaled[tp]
1028 - dxidy_map_scaled[tp]*dzetadz_map_scaled[tp]) +
1029 dzetadx_map_scaled[tp]*( dxidy_map_scaled[tp]* detadz_map_scaled[tp]
1030 -detady_map_scaled[tp]* dxidz_map_scaled[tp]));
1031 if (inv_jacxR_pow4 <= 1e-7)
1033 libmesh_error_msg(
"ERROR: negative weighted inverse Jacobian " \
1041 JxWxdecay[tp] = _total_qrule_weights[tp]/inv_jacxR_pow4;
1047 if (calculate_dphi || calculate_dphi_scaled)
1048 dphase[tp] = unit_r - grad_a_scaled*a/r_norm;
1052 dweight[tp](0) = dweightdv[rp] * dzetadx_map[tp];
1053 dweight[tp](1) = dweightdv[rp] * dzetady_map[tp];
1054 dweight[tp](2) = dweightdv[rp] * dzetadz_map[tp];
1056 if (calculate_dphi_scaled)
1058 dweightxr_sq[tp](0) = dweightdv[rp] * dzetadx_map_scaled[tp];
1059 dweightxr_sq[tp](1) = dweightdv[rp] * dzetady_map_scaled[tp];
1060 dweightxr_sq[tp](2) = dweightdv[rp] * dzetadz_map_scaled[tp];
1063 if (calculate_phi || calculate_phi_scaled || calculate_dphi || calculate_dphi_scaled)
1065 for (
unsigned int i=0; i <_n_total_approx_sf ; ++i)
1068 unsigned int bi = _base_shape_index [i];
1069 unsigned int ri = _radial_shape_index[i];
1071 phi [i][tp] = S [bi][bp] * mode[ri][rp] * som[rp];
1073 if (calculate_phi_scaled)
1074 phixr [i][tp] = S [bi][bp] * mode[ri][rp];
1076 if (calculate_dphi || calculate_dphi_scaled)
1078 dphidxi [i][tp] = Ss[bi][bp] * mode[ri][rp] * som[rp];
1079 dphideta [i][tp] = St[bi][bp] * mode[ri][rp] * som[rp];
1080 dphidzeta[i][tp] = S [bi][bp]
1081 * (dmodedv[ri][rp] * som[rp] + mode[ri][rp] * dsomdv[rp]);
1089 dphidx[i][tp] = (dphidxi[i][tp]*dxidx_map[tp] +
1090 dphideta[i][tp]*detadx_map[tp] +
1091 dphidzeta[i][tp]*dzetadx_map[tp]);
1095 dphidy[i][tp] = (dphidxi[i][tp]*dxidy_map[tp] +
1096 dphideta[i][tp]*detady_map[tp] +
1097 dphidzeta[i][tp]*dzetady_map[tp]);
1101 dphidz[i][tp] = (dphidxi[i][tp]*dxidz_map[tp] +
1102 dphideta[i][tp]*detadz_map[tp] +
1103 dphidzeta[i][tp]*dzetadz_map[tp]);
1106 if (calculate_dphi_scaled)
1109 dphixr[i][tp](0)= (dphidxi[i][tp]*dxidx_map_scaled[tp] +
1110 dphideta[i][tp]*detadx_map_scaled[tp] +
1111 dphidzeta[i][tp]*dzetadx_map_scaled[tp]*som[rp]);
1113 dphixr[i][tp](1) = (dphidxi[i][tp]*dxidy_map_scaled[tp] +
1114 dphideta[i][tp]*detady_map_scaled[tp] +
1115 dphidzeta[i][tp]*dzetady_map_scaled[tp]*som[rp]);
1117 dphixr[i][tp](2) = (dphidxi[i][tp]*dxidz_map_scaled[tp] +
1118 dphideta[i][tp]*detadz_map_scaled[tp] +
1119 dphidzeta[i][tp]*dzetadz_map_scaled[tp]*som[rp]);
1121 const Real dphidxixr = Ss[bi][bp] * mode[ri][rp];
1122 const Real dphidetaxr= St[bi][bp] * mode[ri][rp];
1124 dphixr_sq[i][tp](0)= (dphidxixr*dxidx_map_scaled[tp] +
1125 dphidetaxr*detadx_map_scaled[tp] +
1126 dphidzeta[i][tp]*dzetadx_map_scaled[tp]);
1128 dphixr_sq[i][tp](1) = (dphidxixr*dxidy_map_scaled[tp] +
1129 dphidetaxr*detady_map_scaled[tp] +
1130 dphidzeta[i][tp]*dzetady_map_scaled[tp]);
1132 dphixr_sq[i][tp](2) = (dphidxixr*dxidz_map_scaled[tp] +
1133 dphidetaxr*detadz_map_scaled[tp] +
1134 dphidzeta[i][tp]*dzetadz_map_scaled[tp]);
1145 libmesh_error_msg(
"Unsupported dim = " <<
dim);
1151 template <
unsigned int Dim, FEFamily T_radial, InfMapType T_map>
1155 libmesh_not_implemented();
1163 #include "libmesh/inf_fe_instantiate_1D.h" 1164 #include "libmesh/inf_fe_instantiate_2D.h" 1165 #include "libmesh/inf_fe_instantiate_3D.h" 1169 #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 unsigned int dim, const FEType &fe_t, const ElemType t)
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)
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...