21 #include "libmesh/elem.h" 22 #include "libmesh/fe.h" 23 #include "libmesh/fe_interface.h" 24 #include "libmesh/fe_macro.h" 25 #include "libmesh/libmesh_logging.h" 26 #include "libmesh/quadrature.h" 27 #include "libmesh/tensor_value.h" 28 #include "libmesh/enum_elem_type.h" 29 #include "libmesh/quadrature_gauss.h" 35 void nonlagrange_dual_warning () {
36 libmesh_warning(
"dual calculations have only been verified for the LAGRANGE family");
47 template <
unsigned int Dim, FEFamily T>
56 libmesh_assert_equal_to (T, this->get_family());
60 template <
unsigned int Dim, FEFamily T>
64 static_cast<Order>(this->fe_type.order + this->_p_level));
68 template <
unsigned int Dim, FEFamily T>
79 template <
unsigned int Dim, FEFamily T>
83 return this->qrule->n_points();
87 template <
unsigned int Dim, FEFamily T>
91 std::vector<unsigned int> & di,
92 const bool add_p_level)
95 libmesh_assert_less (s, elem->
n_sides());
98 unsigned int nodenum = 0;
100 for (
unsigned int n = 0; n !=
n_nodes; ++n)
102 const unsigned int n_dofs =
103 n_dofs_at_node(elem->
type(),
static_cast<Order>(o + add_p_level*elem->
p_level()), n);
105 for (
unsigned int i = 0; i != n_dofs; ++i)
106 di.push_back(nodenum++);
114 template <
unsigned int Dim, FEFamily T>
118 std::vector<unsigned int> & di,
119 const bool add_p_level)
122 libmesh_assert_less (e, elem->
n_edges());
125 unsigned int nodenum = 0;
127 for (
unsigned int n = 0; n !=
n_nodes; ++n)
129 const unsigned int n_dofs =
130 n_dofs_at_node(elem->
type(),
static_cast<Order>(o + add_p_level*elem->
p_level()), n);
132 for (
unsigned int i = 0; i != n_dofs; ++i)
133 di.push_back(nodenum++);
141 template <
unsigned int Dim, FEFamily T>
143 const std::vector<Point> *
const pts,
144 const std::vector<Real> *
const weights)
151 this->determine_calculations();
155 bool cached_nodes_still_fit =
false;
165 this->elem_type = elem->
type();
166 this->_elem_p_level = elem->
p_level();
167 this->_p_level = this->_add_p_level_in_reinit * elem->
p_level();
170 this->_fe_map->template init_reference_to_physical_map<Dim>
172 this->init_shape_functions (*pts, elem);
175 this->shapes_on_quadrature =
false;
190 if (this->qrule->shapes_need_reinit())
191 this->shapes_on_quadrature =
false;
195 if (this->elem_type != elem->
type() ||
196 this->_elem_p_level != elem->
p_level() ||
197 !this->shapes_on_quadrature ||
201 this->elem_type = elem->
type();
202 this->_elem_p_level = elem->
p_level();
203 this->_p_level = this->_add_p_level_in_reinit * elem->
p_level();
205 this->_fe_map->template init_reference_to_physical_map<Dim>
206 (this->qrule->get_points(), elem);
207 this->init_shape_functions (this->qrule->get_points(), elem);
209 if (this->shapes_need_reinit())
211 cached_nodes.resize(elem->
n_nodes());
213 cached_nodes[n] = elem->
point(n);
220 cached_nodes_still_fit =
true;
221 if (cached_nodes.size() != elem->
n_nodes())
222 cached_nodes_still_fit =
false;
226 if (!(elem->
point(n) - elem->
point(0)).relative_fuzzy_equals
227 ((cached_nodes[n] - cached_nodes[0]), 1e-13))
229 cached_nodes_still_fit =
false;
234 if (this->shapes_need_reinit() && !cached_nodes_still_fit)
236 this->_fe_map->template init_reference_to_physical_map<Dim>
237 (this->qrule->get_points(), elem);
238 this->init_shape_functions (this->qrule->get_points(), elem);
239 cached_nodes.resize(elem->
n_nodes());
241 cached_nodes[n] = elem->
point(n);
246 this->shapes_on_quadrature =
true;
254 this->_elem_p_level = 0;
261 this->qrule->get_points() =
262 std::vector<Point>(1,
Point(0));
264 this->qrule->get_weights() =
265 std::vector<Real>(1,1);
269 this->qrule->get_points().clear();
270 this->qrule->get_weights().clear();
273 this->init_shape_functions (this->qrule->get_points(), elem);
276 this->init_shape_functions (*pts, elem);
282 if (weights !=
nullptr)
284 this->_fe_map->compute_map (this->
dim, *weights, elem, this->calculate_d2phi);
288 std::vector<Real> dummy_weights (pts->size(), 1.);
289 this->_fe_map->compute_map (this->
dim, dummy_weights, elem, this->calculate_d2phi);
294 this->_fe_map->compute_map (this->
dim, this->qrule->get_weights(), elem, this->calculate_d2phi);
299 if (!cached_nodes_still_fit)
302 this->compute_shape_functions (elem,*pts);
304 this->compute_shape_functions(elem,this->qrule->get_points());
305 if (this->calculate_dual)
308 nonlagrange_dual_warning();
316 if (elem && this->calculate_default_dual_coeff)
317 this->reinit_default_dual_shape_coeffs(elem);
320 this->compute_dual_shape_functions();
325 template <
unsigned int Dim, FEFamily T>
327 const std::vector<Point> & pts,
328 const std::vector<Real> & JxW)
331 this->elem_type = elem->
type();
332 this->_elem_p_level = elem->
p_level();
333 this->_p_level = this->_add_p_level_in_reinit * elem->
p_level();
335 const unsigned int n_shapes =
336 this->n_shape_functions(this->get_type(),
339 std::vector<std::vector<OutputShape>> phi_vals;
340 phi_vals.resize(n_shapes);
341 for (
const auto i :
make_range(phi_vals.size()))
342 phi_vals[i].resize(pts.size());
344 all_shapes(elem, this->get_order(), pts, phi_vals);
345 this->compute_dual_shape_coeffs(JxW, phi_vals);
348 template <
unsigned int Dim, FEFamily T>
353 FEType default_fe_type(this->get_order(), T);
359 this->reinit_dual_shape_coeffs(elem, default_qrule.get_points(), default_qrule.get_weights());
361 this->set_calculate_default_dual_coeff(
false);
365 template <
unsigned int Dim, FEFamily T>
368 if (!this->calculate_dual)
371 libmesh_assert_msg(this->calculate_phi,
372 "dual shape function calculation relies on " 373 "primal shape functions being calculated");
375 this->dual_phi.resize(n_shapes);
376 if (this->calculate_dphi)
377 this->dual_dphi.resize(n_shapes);
378 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 379 if (this->calculate_d2phi)
380 this->dual_d2phi.resize(n_shapes);
385 this->dual_phi[i].resize(n_qp);
386 if (this->calculate_dphi)
387 this->dual_dphi[i].resize(n_qp);
388 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 389 if (this->calculate_d2phi)
390 this->dual_d2phi[i].resize(n_qp);
395 template <
unsigned int Dim, FEFamily T>
400 LOG_SCOPE(
"init_shape_functions()",
"FE");
403 const unsigned int n_qp = cast_int<unsigned int>(qp.size());
407 const unsigned int n_approx_shape_functions =
408 this->n_shape_functions(this->get_type(),
414 unsigned int old_n_qp = 0;
420 if (this->calculate_phi)
422 if (this->phi.size() == n_approx_shape_functions)
424 old_n_qp = n_approx_shape_functions ? this->phi[0].size() : 0;
427 this->phi.resize (n_approx_shape_functions);
429 if (this->calculate_dphi)
431 if (this->dphi.size() == n_approx_shape_functions)
433 old_n_qp = n_approx_shape_functions ? this->dphi[0].size() : 0;
436 this->dphi.resize (n_approx_shape_functions);
437 this->dphidx.resize (n_approx_shape_functions);
438 this->dphidy.resize (n_approx_shape_functions);
439 this->dphidz.resize (n_approx_shape_functions);
442 if (this->calculate_dphiref)
446 if (this->dphidxi.size() == n_approx_shape_functions)
448 old_n_qp = n_approx_shape_functions ? this->dphidxi[0].size() : 0;
451 this->dphidxi.resize (n_approx_shape_functions);
455 this->dphideta.resize (n_approx_shape_functions);
458 this->dphidzeta.resize (n_approx_shape_functions);
461 if (this->calculate_curl_phi && (FEInterface::field_type(T) ==
TYPE_VECTOR))
462 this->curl_phi.resize(n_approx_shape_functions);
464 if (this->calculate_div_phi && (FEInterface::field_type(T) ==
TYPE_VECTOR))
465 this->div_phi.resize(n_approx_shape_functions);
467 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 468 if (this->calculate_d2phi)
470 if (this->d2phi.size() == n_approx_shape_functions)
472 old_n_qp = n_approx_shape_functions ? this->d2phi[0].size() : 0;
476 this->d2phi.resize (n_approx_shape_functions);
477 this->d2phidx2.resize (n_approx_shape_functions);
478 this->d2phidxdy.resize (n_approx_shape_functions);
479 this->d2phidxdz.resize (n_approx_shape_functions);
480 this->d2phidy2.resize (n_approx_shape_functions);
481 this->d2phidydz.resize (n_approx_shape_functions);
482 this->d2phidz2.resize (n_approx_shape_functions);
485 this->d2phidxi2.resize (n_approx_shape_functions);
489 this->d2phidxideta.resize (n_approx_shape_functions);
490 this->d2phideta2.resize (n_approx_shape_functions);
494 this->d2phidxidzeta.resize (n_approx_shape_functions);
495 this->d2phidetadzeta.resize (n_approx_shape_functions);
496 this->d2phidzeta2.resize (n_approx_shape_functions);
499 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 503 if (old_n_qp != n_qp)
504 for (
unsigned int i=0; i<n_approx_shape_functions; i++)
506 if (this->calculate_phi)
507 this->phi[i].resize (n_qp);
509 if (this->calculate_dphi)
511 this->dphi[i].resize (n_qp);
512 this->dphidx[i].resize (n_qp);
513 this->dphidy[i].resize (n_qp);
514 this->dphidz[i].resize (n_qp);
517 if (this->calculate_dphiref)
520 this->dphidxi[i].resize(n_qp);
523 this->dphideta[i].resize(n_qp);
526 this->dphidzeta[i].resize(n_qp);
529 if (this->calculate_curl_phi && (FEInterface::field_type(T) ==
TYPE_VECTOR))
530 this->curl_phi[i].resize(n_qp);
532 if (this->calculate_div_phi && (FEInterface::field_type(T) ==
TYPE_VECTOR))
533 this->div_phi[i].resize(n_qp);
535 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 536 if (this->calculate_d2phi)
538 this->d2phi[i].resize (n_qp);
539 this->d2phidx2[i].resize (n_qp);
540 this->d2phidxdy[i].resize (n_qp);
541 this->d2phidxdz[i].resize (n_qp);
542 this->d2phidy2[i].resize (n_qp);
543 this->d2phidydz[i].resize (n_qp);
544 this->d2phidz2[i].resize (n_qp);
546 this->d2phidxi2[i].resize (n_qp);
549 this->d2phidxideta[i].resize (n_qp);
550 this->d2phideta2[i].resize (n_qp);
554 this->d2phidxidzeta[i].resize (n_qp);
555 this->d2phidetadzeta[i].resize (n_qp);
556 this->d2phidzeta2[i].resize (n_qp);
559 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 563 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 571 if (this->calculate_phi || this->calculate_dphi)
573 this->
weight.resize (n_qp);
574 for (
unsigned int p=0; p<n_qp; p++)
578 if (this->calculate_dphi)
580 this->dweight.resize (n_qp);
581 this->dphase.resize (n_qp);
582 for (
unsigned int p=0; p<n_qp; p++)
584 this->dweight[p].zero();
585 this->dphase[p].zero();
589 #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 592 if (this->calculate_dphiref && Dim > 0)
594 std::vector<std::vector<OutputShape>> * comps[3]
595 { &this->dphidxi, &this->dphideta, &this->dphidzeta };
613 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 615 if (this->calculate_d2phi)
616 for (
unsigned int i=0; i<n_approx_shape_functions; i++)
617 for (
unsigned int p=0; p<n_qp; p++)
619 elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
620 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 631 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 633 if (this->calculate_d2phi)
634 for (
unsigned int i=0; i<n_approx_shape_functions; i++)
635 for (
unsigned int p=0; p<n_qp; p++)
638 elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
640 elem, this->fe_type.order, i, 1, qp[p], this->_add_p_level_in_reinit);
642 elem, this->fe_type.order, i, 2, qp[p], this->_add_p_level_in_reinit);
644 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 656 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 658 if (this->calculate_d2phi)
659 for (
unsigned int i=0; i<n_approx_shape_functions; i++)
660 for (
unsigned int p=0; p<n_qp; p++)
663 elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
665 elem, this->fe_type.order, i, 1, qp[p], this->_add_p_level_in_reinit);
667 elem, this->fe_type.order, i, 2, qp[p], this->_add_p_level_in_reinit);
669 elem, this->fe_type.order, i, 3, qp[p], this->_add_p_level_in_reinit);
671 elem, this->fe_type.order, i, 4, qp[p], this->_add_p_level_in_reinit);
673 elem, this->fe_type.order, i, 5, qp[p], this->_add_p_level_in_reinit);
675 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 682 libmesh_error_msg(
"Invalid dimension Dim = " << Dim);
685 if (this->calculate_dual)
686 this->init_dual_shape_functions(n_approx_shape_functions, n_qp);
689 template <
unsigned int Dim, FEFamily T>
693 const std::vector<Point> & p,
694 std::vector<std::vector<OutputShape>> * comps[3],
695 const bool add_p_level)
697 for (
unsigned int d=0; d != Dim; ++d)
699 auto & comps_d = *comps[d];
702 (elem,o,i,d,p,comps_d[i],add_p_level);
707 template <
unsigned int Dim, FEFamily T>
710 const unsigned int side,
711 const std::vector<Number> & elem_soln,
712 std::vector<Number> & nodal_soln_on_side,
713 const bool add_p_level)
715 std::vector<Number> full_nodal_soln;
716 nodal_soln(elem, o, elem_soln, full_nodal_soln, add_p_level);
717 const std::vector<unsigned int> side_nodes =
720 std::size_t n_side_nodes = side_nodes.size();
721 nodal_soln_on_side.resize(n_side_nodes);
723 nodal_soln_on_side[n] = full_nodal_soln[side_nodes[n]];
729 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS 731 template <
unsigned int Dim, FEFamily T>
735 this->elem_type = e->
type();
736 this->_fe_map->template init_reference_to_physical_map<Dim>(qp, e);
737 init_shape_functions(qp, e);
740 #endif // LIBMESH_ENABLE_INFINITE_ELEMENTS 743 template <
typename OutputShape>
746 const unsigned int i,
747 const unsigned int j,
749 const bool add_p_level,
750 OutputShape(*shape_func)
752 const unsigned int,
const Point &,
757 libmesh_assert_less (j, LIBMESH_DIM);
760 const Real eps = 1.e-6;
761 Point pp = p, pm = p;
790 libmesh_error_msg(
"Invalid derivative index j = " << j);
793 return (shape_func(elem, order, i, pp, add_p_level) -
794 shape_func(elem, order, i, pm, add_p_level))/2./eps;
798 template <
typename OutputShape>
801 const unsigned int i,
802 const unsigned int j,
804 OutputShape(*shape_func)
806 const unsigned int,
const Point &))
808 libmesh_assert_less (j, LIBMESH_DIM);
811 const Real eps = 1.e-6;
812 Point pp = p, pm = p;
841 libmesh_error_msg(
"Invalid derivative index j = " << j);
844 return (shape_func(type, order, i, pp) -
845 shape_func(type, order, i, pm))/2./eps;
849 template <
typename OutputShape>
853 const unsigned int i,
854 const unsigned int j,
856 const bool add_p_level,
857 OutputShape(*deriv_func)
859 const unsigned int,
const unsigned int,
860 const Point &,
const bool))
863 const Real eps = 1.e-5;
864 Point pp = p, pm = p;
865 unsigned int deriv_j = 0;
924 libmesh_error_msg(
"Invalid shape function derivative j = " << j);
927 return (deriv_func(elem, order, i, deriv_j, pp, add_p_level) -
928 deriv_func(elem, order, i, deriv_j, pm, add_p_level))/2./eps;
932 template <
typename OutputShape>
936 const unsigned int i,
937 const unsigned int j,
939 OutputShape(*deriv_func)
941 const unsigned int,
const unsigned int,
945 const Real eps = 1.e-5;
946 Point pp = p, pm = p;
947 unsigned int deriv_j = 0;
1006 libmesh_error_msg(
"Invalid shape function derivative j = " << j);
1009 return (deriv_func(type, order, i, deriv_j, pp) -
1010 deriv_func(type, order, i, deriv_j, pm))/2./eps;
1015 const FEType underlying_fe_type,
1016 std::vector<std::vector<Real>> & shapes,
1017 const std::vector<Point> & p,
1018 const bool add_p_level)
1020 const int extra_order = add_p_level * elem->
p_level();
1022 const int dim = elem->
dim();
1024 const unsigned int n_sf =
1025 FEInterface::n_shape_functions(underlying_fe_type, extra_order,
1028 libmesh_assert_equal_to (n_sf, elem->
n_nodes());
1030 std::vector<Real> node_weights(n_sf);
1032 const unsigned char datum_index = elem->
mapping_data();
1033 for (
unsigned int n=0; n<n_sf; n++)
1037 const std::size_t n_p = p.size();
1039 shapes.resize(n_sf);
1040 for (
unsigned int i=0; i != n_sf; ++i)
1042 auto & shapes_i = shapes[i];
1043 shapes_i.resize(n_p, 0);
1044 FEInterface::shapes(
dim, underlying_fe_type, elem, i, p,
1045 shapes_i, add_p_level);
1046 for (
auto & s : shapes_i)
1047 s *= node_weights[i];
1054 std::vector<std::vector<Real>> & shapes,
1055 std::vector<std::vector<std::vector<Real>>> & derivs,
1056 const std::vector<Point> & p,
1057 const bool add_p_level)
1059 const int extra_order = add_p_level * elem->
p_level();
1060 const unsigned int dim = elem->
dim();
1062 const unsigned int n_sf =
1063 FEInterface::n_shape_functions(fe_type, extra_order, elem);
1065 libmesh_assert_equal_to (n_sf, elem->
n_nodes());
1067 libmesh_assert_equal_to (
dim, derivs.size());
1068 for (
unsigned int d = 0; d !=
dim; ++d)
1069 derivs[d].resize(n_sf);
1071 std::vector<Real> node_weights(n_sf);
1073 const unsigned char datum_index = elem->
mapping_data();
1074 for (
unsigned int n=0; n<n_sf; n++)
1078 const std::size_t n_p = p.size();
1080 shapes.resize(n_sf);
1081 for (
unsigned int i=0; i != n_sf; ++i)
1082 shapes[i].resize(n_p, 0);
1084 FEInterface::all_shapes(
dim, fe_type, elem, p, shapes, add_p_level);
1086 for (
unsigned int i=0; i != n_sf; ++i)
1088 auto & shapes_i = shapes[i];
1090 for (
auto & s : shapes_i)
1091 s *= node_weights[i];
1093 for (
unsigned int d = 0; d !=
dim; ++d)
1095 auto & derivs_di = derivs[d][i];
1096 derivs_di.resize(n_p);
1097 FEInterface::shape_derivs(fe_type, elem, i, d, p,
1098 derivs_di, add_p_level);
1099 for (
auto & dip : derivs_di)
1100 dip *= node_weights[i];
1107 const FEType underlying_fe_type,
1108 const unsigned int i,
1110 const bool add_p_level)
1112 int extra_order = add_p_level * elem.
p_level();
1114 const unsigned int n_sf =
1115 FEInterface::n_shape_functions(underlying_fe_type, extra_order, &elem);
1117 libmesh_assert_equal_to (n_sf, elem.
n_nodes());
1119 std::vector<Real> node_weights(n_sf);
1123 Real weighted_shape_i = 0, weighted_sum = 0;
1125 for (
unsigned int sf=0; sf<n_sf; sf++)
1129 Real weighted_shape = node_weight *
1130 FEInterface::shape(underlying_fe_type, extra_order, &elem, sf, p);
1131 weighted_sum += weighted_shape;
1133 weighted_shape_i = weighted_shape;
1136 return weighted_shape_i / weighted_sum;
1141 const FEType underlying_fe_type,
1142 const unsigned int i,
1143 const unsigned int j,
1145 const bool add_p_level)
1147 libmesh_assert_less(j, elem.
dim());
1149 int extra_order = add_p_level * elem.
p_level();
1151 const unsigned int n_sf =
1152 FEInterface::n_shape_functions(underlying_fe_type, extra_order, &elem);
1155 libmesh_assert_equal_to (n_sf,
n_nodes);
1157 std::vector<Real> node_weights(
n_nodes);
1160 for (
unsigned int n=0; n<
n_nodes; n++)
1164 Real weighted_shape_i = 0, weighted_sum = 0,
1165 weighted_grad_i = 0, weighted_grad_sum = 0;
1167 for (
unsigned int sf=0; sf<n_sf; sf++)
1169 Real weighted_shape = node_weights[sf] *
1170 FEInterface::shape(underlying_fe_type, extra_order, &elem, sf, p);
1171 Real weighted_grad = node_weights[sf] *
1172 FEInterface::shape_deriv(underlying_fe_type, extra_order, &elem, sf, j, p);
1173 weighted_sum += weighted_shape;
1174 weighted_grad_sum += weighted_grad;
1177 weighted_shape_i = weighted_shape;
1178 weighted_grad_i = weighted_grad;
1182 return (weighted_sum * weighted_grad_i - weighted_shape_i * weighted_grad_sum) /
1183 weighted_sum / weighted_sum;
1187 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 1190 const FEType underlying_fe_type,
1191 const unsigned int i,
1192 const unsigned int j,
1194 const bool add_p_level)
1196 unsigned int j1, j2;
1230 int extra_order = add_p_level * elem.
p_level();
1232 const unsigned int n_sf =
1233 FEInterface::n_shape_functions(underlying_fe_type, extra_order,
1237 libmesh_assert_equal_to (n_sf,
n_nodes);
1239 std::vector<Real> node_weights(
n_nodes);
1242 for (
unsigned int n=0; n<
n_nodes; n++)
1246 Real weighted_shape_i = 0, weighted_sum = 0,
1247 weighted_grada_i = 0, weighted_grada_sum = 0,
1248 weighted_gradb_i = 0, weighted_gradb_sum = 0,
1249 weighted_hess_i = 0, weighted_hess_sum = 0;
1251 for (
unsigned int sf=0; sf<n_sf; sf++)
1253 Real weighted_shape = node_weights[sf] *
1254 FEInterface::shape(underlying_fe_type, extra_order, &elem, sf,
1256 Real weighted_grada = node_weights[sf] *
1257 FEInterface::shape_deriv(underlying_fe_type, extra_order,
1259 Real weighted_hess = node_weights[sf] *
1260 FEInterface::shape_second_deriv(underlying_fe_type,
1261 extra_order, &elem, sf, j, p);
1262 weighted_sum += weighted_shape;
1263 weighted_grada_sum += weighted_grada;
1264 Real weighted_gradb = weighted_grada;
1267 weighted_gradb = (j1 == j2) ? weighted_grada :
1269 FEInterface::shape_deriv(underlying_fe_type, extra_order,
1271 weighted_grada_sum += weighted_grada;
1273 weighted_hess_sum += weighted_hess;
1276 weighted_shape_i = weighted_shape;
1277 weighted_grada_i = weighted_grada;
1278 weighted_gradb_i = weighted_gradb;
1279 weighted_hess_i = weighted_hess;
1284 weighted_gradb_sum = weighted_grada_sum;
1286 return (weighted_sum * weighted_hess_i - weighted_grada_i * weighted_gradb_sum -
1287 weighted_shape_i * weighted_hess_sum - weighted_gradb_i * weighted_grada_sum +
1288 2 * weighted_grada_sum * weighted_shape_i * weighted_gradb_sum / weighted_sum) /
1289 weighted_sum / weighted_sum;
1292 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES 1296 const FEType underlying_fe_type,
1297 const std::vector<Point> & p,
1298 std::vector<std::vector<Real>> & v,
1299 const bool add_p_level)
1301 std::vector<std::vector<Real>> shapes;
1306 std::vector<Real> shape_sums(p.size(), 0);
1310 libmesh_assert_equal_to ( p.size(), shapes[i].size() );
1312 shape_sums[j] += shapes[i][j];
1317 libmesh_assert_equal_to ( p.size(), v[i].size() );
1319 v[i][j] = shapes[i][j] / shape_sums[j];
1324 template <
typename OutputShape>
1326 const FEType underlying_fe_type,
1327 const std::vector<Point> & p,
1328 std::vector<std::vector<OutputShape>> * comps[3],
1329 const bool add_p_level)
1331 const int my_dim = elem.
dim();
1333 std::vector<std::vector<Real>> shapes;
1334 std::vector<std::vector<std::vector<Real>>> derivs(my_dim);
1337 shapes, derivs, p, add_p_level);
1339 std::vector<Real> shape_sums(p.size(), 0);
1340 std::vector<std::vector<Real>> shape_deriv_sums(my_dim);
1341 for (
int d=0; d != my_dim; ++d)
1342 shape_deriv_sums[d].resize(p.size());
1346 libmesh_assert_equal_to ( p.size(), shapes[i].size() );
1348 shape_sums[j] += shapes[i][j];
1350 for (
int d=0; d != my_dim; ++d)
1352 shape_deriv_sums[d][j] += derivs[d][i][j];
1355 for (
int d=0; d != my_dim; ++d)
1357 auto & comps_d = *comps[d];
1358 libmesh_assert_equal_to(comps_d.size(), elem.
n_nodes());
1362 auto & comps_di = comps_d[i];
1363 auto & derivs_di = derivs[d][i];
1366 comps_di[j] = (shape_sums[j] * derivs_di[j] -
1367 shapes[i][j] * shape_deriv_sums[d][j]) /
1368 shape_sums[j] / shape_sums[j];
1377 const unsigned int,
const Point &,
const bool,
1379 (
const Elem *,
const Order,
const unsigned int,
1380 const Point &,
const bool));
1384 const unsigned int,
const Point &,
1392 const unsigned int,
const Point &,
const bool,
1394 (
const Elem *,
const Order,
const unsigned int,
1395 const Point &,
const bool));
1400 const unsigned int,
const Point &,
1403 const unsigned int,
const Point &));
1408 const unsigned int,
const Point &,
const bool,
1410 (
const Elem *,
const Order,
const unsigned int,
1411 const unsigned int,
const Point &,
const bool));
1416 const unsigned int,
const Point &,
const bool,
1418 (
const Elem *,
const Order,
const unsigned int,
1419 const unsigned int,
const Point &,
const bool));
1437 template LIBMESH_EXPORT
void rational_all_shape_derivs<double> (
const Elem & elem,
const FEType underlying_fe_type,
const std::vector<Point> & p, std::vector<std::vector<Real>> * comps[3],
const bool add_p_level);
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
virtual void init(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Initializes the data structures for a quadrature rule for an element of type type.
unsigned char mapping_data() const
ElemType
Defines an enum for geometric element types.
RealVectorValue RealGradient
Order
defines an enum for polynomial orders.
void rational_fe_weighted_shapes_derivs(const Elem *elem, const FEType fe_type, std::vector< std::vector< Real >> &shapes, std::vector< std::vector< std::vector< Real >>> &derivs, const std::vector< Point > &p, const bool add_p_level)
INSTANTIATE_SUBDIVISION_FE
This is the base class from which all geometric element types are derived.
template Real fe_fdm_second_deriv< Real >(const Elem *, const Order, const unsigned int, const unsigned int, const Point &, const bool, Real(*shape_func)(const Elem *, const Order, const unsigned int, const unsigned int, const Point &, const bool))
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
FE(const FEType &fet)
Constructor.
Order default_quadrature_order() const
unsigned int p_level() const
The libMesh namespace provides an interface to certain functionality in the library.
OutputShape fe_fdm_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*shape_func)(const Elem *, const Order, const unsigned int, const Point &, const bool))
Helper functions for finite differenced derivatives in cases where analytical calculations haven't be...
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const =0
A specific instantiation of the FEBase class.
template RealGradient fe_fdm_deriv< RealGradient >(const Elem *, const Order, const unsigned int, const unsigned int, const Point &, const bool, RealGradient(*shape_func)(const Elem *, const Order, const unsigned int, const Point &, const bool))
const dof_id_type n_nodes
ElemMappingType mapping_type() const
const Node & node_ref(const unsigned int i) const
virtual unsigned int n_nodes() const =0
void rational_all_shape_derivs(const Elem &elem, const FEType underlying_fe_type, const std::vector< Point > &p, std::vector< std::vector< OutputShape >> *comps[3], const bool add_p_level)
virtual unsigned int n_edges() const =0
OutputShape fe_fdm_second_deriv(const Elem *elem, const Order order, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level, OutputShape(*deriv_func)(const Elem *, const Order, const unsigned int, const unsigned int, const Point &, const bool))
virtual unsigned int n_sides() const =0
Real rational_fe_shape(const Elem &elem, const FEType underlying_fe_type, const unsigned int i, const Point &p, const bool add_p_level)
Real rational_fe_shape_second_deriv(const Elem &elem, const FEType underlying_fe_type, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level)
template RealGradient fe_fdm_second_deriv< RealGradient >(const Elem *, const Order, const unsigned int, const unsigned int, const Point &, const bool, RealGradient(*shape_func)(const Elem *, const Order, const unsigned int, const unsigned int, const Point &, const bool))
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned short dim() const =0
template LIBMESH_EXPORT void rational_all_shape_derivs< double >(const Elem &elem, const FEType underlying_fe_type, const std::vector< Point > &p, std::vector< std::vector< Real >> *comps[3], const bool add_p_level)
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
This class implements specific orders of Gauss quadrature.
IntRange< unsigned short > node_index_range() const
T get_extra_datum(const unsigned int index) const
Gets the value on this object of the extra datum associated with index, which should have been obtain...
void rational_all_shapes(const Elem &elem, const FEType underlying_fe_type, const std::vector< Point > &p, std::vector< std::vector< Real >> &v, const bool add_p_level)
template Real fe_fdm_deriv< Real >(const ElemType, const Order, const unsigned int, const unsigned int, const Point &, Real(*shape_func)(const ElemType, const Order, const unsigned int, const Point &))
void rational_fe_weighted_shapes(const Elem *elem, const FEType underlying_fe_type, std::vector< std::vector< Real >> &shapes, const std::vector< Point > &p, const bool add_p_level)
Helper functions for rational basis functions.
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
const Point & point(const unsigned int i) const
The QBase class provides the basic functionality from which various quadrature rules can be derived...
void ErrorVector unsigned int
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
This class forms the foundation from which generic finite elements may be derived.
virtual std::vector< unsigned int > nodes_on_side(const unsigned int) const =0
Most finite element types in libMesh are scalar-valued.
Real rational_fe_shape_deriv(const Elem &elem, const FEType underlying_fe_type, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level)