libMesh
fe.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // Local includes
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"
30 
31 namespace {
32  // Put this outside a templated class, so we only get 1 warning
33  // during our unit tests, not 1 warning for each of the zillion FE
34  // specializations we test.
35  void nonlagrange_dual_warning () {
36  libmesh_warning("dual calculations have only been verified for the LAGRANGE family");
37  }
38 }
39 
40 
41 namespace libMesh
42 {
43 
44 
45 // ------------------------------------------------------------
46 // FE class members
47 template <unsigned int Dim, FEFamily T>
48 FE<Dim,T>::FE (const FEType & fet) :
49  FEGenericBase<typename FEOutputType<T>::type> (Dim,fet),
50  last_side(INVALID_ELEM),
51  last_edge(INVALID_ELEM)
52 {
53  // Sanity check. Make sure the
54  // Family specified in the template instantiation
55  // matches the one in the FEType object
56  libmesh_assert_equal_to (T, this->get_family());
57 }
58 
59 
60 template <unsigned int Dim, FEFamily T>
61 unsigned int FE<Dim,T>::n_shape_functions () const
62 {
63  return FE<Dim,T>::n_dofs (this->elem_type,
64  static_cast<Order>(this->fe_type.order + this->_p_level));
65 }
66 
67 
68 template <unsigned int Dim, FEFamily T>
70 {
71  libmesh_assert(q);
72  this->qrule = q;
73  // make sure we don't cache results from a previous quadrature rule
74  this->elem_type = INVALID_ELEM;
75  return;
76 }
77 
78 
79 template <unsigned int Dim, FEFamily T>
80 unsigned int FE<Dim,T>::n_quadrature_points () const
81 {
82  libmesh_assert(this->qrule);
83  return this->qrule->n_points();
84 }
85 
86 
87 template <unsigned int Dim, FEFamily T>
88 void FE<Dim,T>::dofs_on_side(const Elem * const elem,
89  const Order o,
90  unsigned int s,
91  std::vector<unsigned int> & di,
92  const bool add_p_level)
93 {
94  libmesh_assert(elem);
95  libmesh_assert_less (s, elem->n_sides());
96 
97  di.clear();
98  unsigned int nodenum = 0;
99  const unsigned int n_nodes = elem->n_nodes();
100  for (unsigned int n = 0; n != n_nodes; ++n)
101  {
102  const unsigned int n_dofs =
103  n_dofs_at_node(elem->type(), static_cast<Order>(o + add_p_level*elem->p_level()), n);
104  if (elem->is_node_on_side(n, s))
105  for (unsigned int i = 0; i != n_dofs; ++i)
106  di.push_back(nodenum++);
107  else
108  nodenum += n_dofs;
109  }
110 }
111 
112 
113 
114 template <unsigned int Dim, FEFamily T>
115 void FE<Dim,T>::dofs_on_edge(const Elem * const elem,
116  const Order o,
117  unsigned int e,
118  std::vector<unsigned int> & di,
119  const bool add_p_level)
120 {
121  libmesh_assert(elem);
122  libmesh_assert_less (e, elem->n_edges());
123 
124  di.clear();
125  unsigned int nodenum = 0;
126  const unsigned int n_nodes = elem->n_nodes();
127  for (unsigned int n = 0; n != n_nodes; ++n)
128  {
129  const unsigned int n_dofs =
130  n_dofs_at_node(elem->type(), static_cast<Order>(o + add_p_level*elem->p_level()), n);
131  if (elem->is_node_on_edge(n, e))
132  for (unsigned int i = 0; i != n_dofs; ++i)
133  di.push_back(nodenum++);
134  else
135  nodenum += n_dofs;
136  }
137 }
138 
139 
140 
141 template <unsigned int Dim, FEFamily T>
142 void FE<Dim,T>::reinit(const Elem * elem,
143  const std::vector<Point> * const pts,
144  const std::vector<Real> * const weights)
145 {
146  // We can be called with no element. If we're evaluating SCALAR
147  // dofs we'll still have work to do.
148  // libmesh_assert(elem);
149 
150  // We're calculating now! Time to determine what.
151  this->determine_calculations();
152 
153  // Try to avoid calling init_shape_functions
154  // even when shapes_need_reinit
155  bool cached_nodes_still_fit = false;
156 
157  // Most of the hard work happens when we have an actual element
158  if (elem)
159  {
160  // Initialize the shape functions at the user-specified
161  // points
162  if (pts != nullptr)
163  {
164  // Set the type and p level for this element
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();
168 
169  // Initialize the shape functions
170  this->_fe_map->template init_reference_to_physical_map<Dim>
171  (*pts, elem);
172  this->init_shape_functions (*pts, elem);
173 
174  // The shape functions do not correspond to the qrule
175  this->shapes_on_quadrature = false;
176  }
177 
178  // If there are no user specified points, we use the
179  // quadrature rule
180 
181  // update the type in accordance to the current cell
182  // and reinit if the cell type has changed or (as in
183  // the case of the hierarchics) the shape functions need
184  // reinit, since they depend on the particular element shape
185  else
186  {
187  libmesh_assert(this->qrule);
188  this->qrule->init(elem->type(), elem->p_level());
189 
190  if (this->qrule->shapes_need_reinit())
191  this->shapes_on_quadrature = false;
192 
193  // We're not going to bother trying to cache nodal
194  // points *and* weights for fancier mapping types.
195  if (this->elem_type != elem->type() ||
196  this->_elem_p_level != elem->p_level() ||
197  !this->shapes_on_quadrature ||
198  elem->mapping_type() != LAGRANGE_MAP)
199  {
200  // Set the type and p level for this element
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();
204  // Initialize the shape functions
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);
208 
209  if (this->shapes_need_reinit())
210  {
211  cached_nodes.resize(elem->n_nodes());
212  for (auto n : elem->node_index_range())
213  cached_nodes[n] = elem->point(n);
214  }
215  }
216  else
217  {
218  // libmesh_assert_greater (elem->n_nodes(), 1);
219 
220  cached_nodes_still_fit = true;
221  if (cached_nodes.size() != elem->n_nodes())
222  cached_nodes_still_fit = false;
223  else
224  for (auto n : make_range(1u, elem->n_nodes()))
225  {
226  if (!(elem->point(n) - elem->point(0)).relative_fuzzy_equals
227  ((cached_nodes[n] - cached_nodes[0]), 1e-13))
228  {
229  cached_nodes_still_fit = false;
230  break;
231  }
232  }
233 
234  if (this->shapes_need_reinit() && !cached_nodes_still_fit)
235  {
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());
240  for (auto n : elem->node_index_range())
241  cached_nodes[n] = elem->point(n);
242  }
243  }
244 
245  // The shape functions correspond to the qrule
246  this->shapes_on_quadrature = true;
247  }
248  }
249  else // With no defined elem, so mapping or caching to
250  // be done, and our "quadrature rule" is one point for nonlocal
251  // (SCALAR) variables and zero points for local variables.
252  {
253  this->elem_type = INVALID_ELEM;
254  this->_elem_p_level = 0;
255  this->_p_level = 0;
256 
257  if (!pts)
258  {
259  if (T == SCALAR)
260  {
261  this->qrule->get_points() =
262  std::vector<Point>(1,Point(0));
263 
264  this->qrule->get_weights() =
265  std::vector<Real>(1,1);
266  }
267  else
268  {
269  this->qrule->get_points().clear();
270  this->qrule->get_weights().clear();
271  }
272 
273  this->init_shape_functions (this->qrule->get_points(), elem);
274  }
275  else
276  this->init_shape_functions (*pts, elem);
277  }
278 
279  // Compute the map for this element.
280  if (pts != nullptr)
281  {
282  if (weights != nullptr)
283  {
284  this->_fe_map->compute_map (this->dim, *weights, elem, this->calculate_d2phi);
285  }
286  else
287  {
288  std::vector<Real> dummy_weights (pts->size(), 1.);
289  this->_fe_map->compute_map (this->dim, dummy_weights, elem, this->calculate_d2phi);
290  }
291  }
292  else
293  {
294  this->_fe_map->compute_map (this->dim, this->qrule->get_weights(), elem, this->calculate_d2phi);
295  }
296 
297  // Compute the shape functions and the derivatives at all of the
298  // quadrature points.
299  if (!cached_nodes_still_fit)
300  {
301  if (pts != nullptr)
302  this->compute_shape_functions (elem,*pts);
303  else
304  this->compute_shape_functions(elem,this->qrule->get_points());
305  if (this->calculate_dual)
306  {
307  if (T != LAGRANGE)
308  nonlagrange_dual_warning();
309  // Check if we need to calculate the dual coefficients based on the default QRule
310  // We keep the default dual coeff calculation for the initial stage of the simulation
311  // and in the middel of the simulation when a customized QRule is not provided.
312  // This is used in MOOSE mortar-based contact. Currently, we re-compute dual_coeff
313  // for all the elements on the mortar segment mesh by setting `calculate_default_dual_coeff' = false
314  // in MOOSE (in `Assembly::reinitDual`) and use the customized QRule for calculating the dual shape coefficients
315  // This is to be improved in the future
316  if (elem && this->calculate_default_dual_coeff)
317  this->reinit_default_dual_shape_coeffs(elem);
318  // The dual shape functions relies on the customized shape functions
319  // and the coefficient matrix, \p dual_coeff
320  this->compute_dual_shape_functions();
321  }
322  }
323 }
324 
325 template <unsigned int Dim, FEFamily T>
327  const std::vector<Point> & pts,
328  const std::vector<Real> & JxW)
329 {
330  // Set the type and p level for this element
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();
334 
335  const unsigned int n_shapes =
336  this->n_shape_functions(this->get_type(),
337  this->get_order());
338 
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());
343 
344  all_shapes(elem, this->get_order(), pts, phi_vals);
345  this->compute_dual_shape_coeffs(JxW, phi_vals);
346 }
347 
348 template <unsigned int Dim, FEFamily T>
350 {
351  libmesh_assert(elem);
352 
353  FEType default_fe_type(this->get_order(), T);
354  QGauss default_qrule(elem->dim(), default_fe_type.default_quadrature_order());
355  default_qrule.init(elem->type(), elem->p_level());
356  // In preparation of computing dual_coeff, we compute the default shape
357  // function values and use these to compute the dual shape coefficients.
358  // The TRUE dual_phi values are computed in compute_dual_shape_functions()
359  this->reinit_dual_shape_coeffs(elem, default_qrule.get_points(), default_qrule.get_weights());
360  // we do not compute default dual coeff many times as this can be expensive
361  this->set_calculate_default_dual_coeff(false);
362 }
363 
364 
365 template <unsigned int Dim, FEFamily T>
366 void FE<Dim,T>::init_dual_shape_functions(const unsigned int n_shapes, const unsigned int n_qp)
367 {
368  if (!this->calculate_dual)
369  return;
370 
371  libmesh_assert_msg(this->calculate_phi,
372  "dual shape function calculation relies on "
373  "primal shape functions being calculated");
374 
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);
381 #endif
382 
383  for (auto i : index_range(this->dual_phi))
384  {
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);
391 #endif
392  }
393 }
394 
395 template <unsigned int Dim, FEFamily T>
396 void FE<Dim,T>::init_shape_functions(const std::vector<Point> & qp,
397  const Elem * elem)
398 {
399  // Start logging the shape function initialization
400  LOG_SCOPE("init_shape_functions()", "FE");
401 
402  // The number of quadrature points.
403  const unsigned int n_qp = cast_int<unsigned int>(qp.size());
404 
405  // Number of shape functions in the finite element approximation
406  // space.
407  const unsigned int n_approx_shape_functions =
408  this->n_shape_functions(this->get_type(),
409  this->get_order());
410 
411  // Maybe we already have correctly-sized data? Check data sizes,
412  // and get ready to break out of a "loop" if all these resize()
413  // calls are redundant.
414  unsigned int old_n_qp = 0;
415  do
416  {
417  // resize the vectors to hold current data
418  // Phi are the shape functions used for the FE approximation
419  // Phi_map are the shape functions used for the FE mapping
420  if (this->calculate_phi)
421  {
422  if (this->phi.size() == n_approx_shape_functions)
423  {
424  old_n_qp = n_approx_shape_functions ? this->phi[0].size() : 0;
425  break;
426  }
427  this->phi.resize (n_approx_shape_functions);
428  }
429  if (this->calculate_dphi)
430  {
431  if (this->dphi.size() == n_approx_shape_functions)
432  {
433  old_n_qp = n_approx_shape_functions ? this->dphi[0].size() : 0;
434  break;
435  }
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);
440  }
441 
442  if (this->calculate_dphiref)
443  {
444  if (Dim > 0)
445  {
446  if (this->dphidxi.size() == n_approx_shape_functions)
447  {
448  old_n_qp = n_approx_shape_functions ? this->dphidxi[0].size() : 0;
449  break;
450  }
451  this->dphidxi.resize (n_approx_shape_functions);
452  }
453 
454  if (Dim > 1)
455  this->dphideta.resize (n_approx_shape_functions);
456 
457  if (Dim > 2)
458  this->dphidzeta.resize (n_approx_shape_functions);
459  }
460 
461  if (this->calculate_curl_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
462  this->curl_phi.resize(n_approx_shape_functions);
463 
464  if (this->calculate_div_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
465  this->div_phi.resize(n_approx_shape_functions);
466 
467 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
468  if (this->calculate_d2phi)
469  {
470  if (this->d2phi.size() == n_approx_shape_functions)
471  {
472  old_n_qp = n_approx_shape_functions ? this->d2phi[0].size() : 0;
473  break;
474  }
475 
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);
483 
484  if (Dim > 0)
485  this->d2phidxi2.resize (n_approx_shape_functions);
486 
487  if (Dim > 1)
488  {
489  this->d2phidxideta.resize (n_approx_shape_functions);
490  this->d2phideta2.resize (n_approx_shape_functions);
491  }
492  if (Dim > 2)
493  {
494  this->d2phidxidzeta.resize (n_approx_shape_functions);
495  this->d2phidetadzeta.resize (n_approx_shape_functions);
496  this->d2phidzeta2.resize (n_approx_shape_functions);
497  }
498  }
499 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
500  }
501  while (false);
502 
503  if (old_n_qp != n_qp)
504  for (unsigned int i=0; i<n_approx_shape_functions; i++)
505  {
506  if (this->calculate_phi)
507  this->phi[i].resize (n_qp);
508 
509  if (this->calculate_dphi)
510  {
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);
515  }
516 
517  if (this->calculate_dphiref)
518  {
519  if (Dim > 0)
520  this->dphidxi[i].resize(n_qp);
521 
522  if (Dim > 1)
523  this->dphideta[i].resize(n_qp);
524 
525  if (Dim > 2)
526  this->dphidzeta[i].resize(n_qp);
527  }
528 
529  if (this->calculate_curl_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
530  this->curl_phi[i].resize(n_qp);
531 
532  if (this->calculate_div_phi && (FEInterface::field_type(T) == TYPE_VECTOR))
533  this->div_phi[i].resize(n_qp);
534 
535 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
536  if (this->calculate_d2phi)
537  {
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);
545  if (Dim > 0)
546  this->d2phidxi2[i].resize (n_qp);
547  if (Dim > 1)
548  {
549  this->d2phidxideta[i].resize (n_qp);
550  this->d2phideta2[i].resize (n_qp);
551  }
552  if (Dim > 2)
553  {
554  this->d2phidxidzeta[i].resize (n_qp);
555  this->d2phidetadzeta[i].resize (n_qp);
556  this->d2phidzeta2[i].resize (n_qp);
557  }
558  }
559 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
560  }
561 
562 
563 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
564  //------------------------------------------------------------
565  // Initialize the data fields, which should only be used for infinite
566  // elements, to some sensible values, so that using a FE with the
567  // variational formulation of an InfFE, correct element matrices are
568  // returned
569 
570  {
571  if (this->calculate_phi || this->calculate_dphi)
572  {
573  this->weight.resize (n_qp);
574  for (unsigned int p=0; p<n_qp; p++)
575  this->weight[p] = 1.;
576  }
577 
578  if (this->calculate_dphi)
579  {
580  this->dweight.resize (n_qp);
581  this->dphase.resize (n_qp);
582  for (unsigned int p=0; p<n_qp; p++)
583  {
584  this->dweight[p].zero();
585  this->dphase[p].zero();
586  }
587  }
588  }
589 #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
590 
591  // Compute the values of the shape function derivatives
592  if (this->calculate_dphiref && Dim > 0)
593  {
594  std::vector<std::vector<OutputShape>> * comps[3]
595  { &this->dphidxi, &this->dphideta, &this->dphidzeta };
596  FE<Dim,T>::all_shape_derivs(elem, this->fe_type.order, qp, comps, this->_add_p_level_in_reinit);
597  }
598 
599  switch (Dim)
600  {
601 
602  //------------------------------------------------------------
603  // 0D
604  case 0:
605  {
606  break;
607  }
608 
609  //------------------------------------------------------------
610  // 1D
611  case 1:
612  {
613 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
614  // Compute the value of shape function i Hessians at quadrature point p
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++)
618  this->d2phidxi2[i][p] = FE<Dim, T>::shape_second_deriv(
619  elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
620 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
621 
622  break;
623  }
624 
625 
626 
627  //------------------------------------------------------------
628  // 2D
629  case 2:
630  {
631 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
632  // Compute the value of shape function i Hessians at quadrature point p
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++)
636  {
637  this->d2phidxi2[i][p] = FE<Dim, T>::shape_second_deriv(
638  elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
639  this->d2phidxideta[i][p] = FE<Dim, T>::shape_second_deriv(
640  elem, this->fe_type.order, i, 1, qp[p], this->_add_p_level_in_reinit);
641  this->d2phideta2[i][p] = FE<Dim, T>::shape_second_deriv(
642  elem, this->fe_type.order, i, 2, qp[p], this->_add_p_level_in_reinit);
643  }
644 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
645 
646 
647  break;
648  }
649 
650 
651 
652  //------------------------------------------------------------
653  // 3D
654  case 3:
655  {
656 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
657  // Compute the value of shape function i Hessians at quadrature point p
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++)
661  {
662  this->d2phidxi2[i][p] = FE<Dim, T>::shape_second_deriv(
663  elem, this->fe_type.order, i, 0, qp[p], this->_add_p_level_in_reinit);
664  this->d2phidxideta[i][p] = FE<Dim, T>::shape_second_deriv(
665  elem, this->fe_type.order, i, 1, qp[p], this->_add_p_level_in_reinit);
666  this->d2phideta2[i][p] = FE<Dim, T>::shape_second_deriv(
667  elem, this->fe_type.order, i, 2, qp[p], this->_add_p_level_in_reinit);
668  this->d2phidxidzeta[i][p] = FE<Dim, T>::shape_second_deriv(
669  elem, this->fe_type.order, i, 3, qp[p], this->_add_p_level_in_reinit);
670  this->d2phidetadzeta[i][p] = FE<Dim, T>::shape_second_deriv(
671  elem, this->fe_type.order, i, 4, qp[p], this->_add_p_level_in_reinit);
672  this->d2phidzeta2[i][p] = FE<Dim, T>::shape_second_deriv(
673  elem, this->fe_type.order, i, 5, qp[p], this->_add_p_level_in_reinit);
674  }
675 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
676 
677  break;
678  }
679 
680 
681  default:
682  libmesh_error_msg("Invalid dimension Dim = " << Dim);
683  }
684 
685  if (this->calculate_dual)
686  this->init_dual_shape_functions(n_approx_shape_functions, n_qp);
687 }
688 
689 template <unsigned int Dim, FEFamily T>
690 void
692  const Order o,
693  const std::vector<Point> & p,
694  std::vector<std::vector<OutputShape>> * comps[3],
695  const bool add_p_level)
696 {
697  for (unsigned int d=0; d != Dim; ++d)
698  {
699  auto & comps_d = *comps[d];
700  for (auto i : index_range(comps_d))
702  (elem,o,i,d,p,comps_d[i],add_p_level);
703  }
704 }
705 
706 
707 template <unsigned int Dim, FEFamily T>
708 void
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)
714 {
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 =
718  elem->nodes_on_side(side);
719 
720  std::size_t n_side_nodes = side_nodes.size();
721  nodal_soln_on_side.resize(n_side_nodes);
722  for (auto n : make_range(n_side_nodes))
723  nodal_soln_on_side[n] = full_nodal_soln[side_nodes[n]];
724 }
725 
726 
727 
728 
729 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
730 
731 template <unsigned int Dim, FEFamily T>
732 void FE<Dim,T>::init_base_shape_functions(const std::vector<Point> & qp,
733  const Elem * e)
734 {
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);
738 }
739 
740 #endif // LIBMESH_ENABLE_INFINITE_ELEMENTS
741 
742 
743 template <typename OutputShape>
744 OutputShape fe_fdm_deriv(const Elem * elem,
745  const Order order,
746  const unsigned int i,
747  const unsigned int j,
748  const Point & p,
749  const bool add_p_level,
750  OutputShape(*shape_func)
751  (const Elem *, const Order,
752  const unsigned int, const Point &,
753  const bool))
754 {
755  libmesh_assert(elem);
756 
757  libmesh_assert_less (j, LIBMESH_DIM);
758 
759  // cheat by using finite difference approximations:
760  const Real eps = 1.e-6;
761  Point pp = p, pm = p;
762 
763  switch (j)
764  {
765  // d()/dxi
766  case 0:
767  {
768  pp(0) += eps;
769  pm(0) -= eps;
770  break;
771  }
772 
773  // d()/deta
774  case 1:
775  {
776  pp(1) += eps;
777  pm(1) -= eps;
778  break;
779  }
780 
781  // d()/dzeta
782  case 2:
783  {
784  pp(2) += eps;
785  pm(2) -= eps;
786  break;
787  }
788 
789  default:
790  libmesh_error_msg("Invalid derivative index j = " << j);
791  }
792 
793  return (shape_func(elem, order, i, pp, add_p_level) -
794  shape_func(elem, order, i, pm, add_p_level))/2./eps;
795 }
796 
797 
798 template <typename OutputShape>
799 OutputShape fe_fdm_deriv(const ElemType type,
800  const Order order,
801  const unsigned int i,
802  const unsigned int j,
803  const Point & p,
804  OutputShape(*shape_func)
805  (const ElemType, const Order,
806  const unsigned int, const Point &))
807 {
808  libmesh_assert_less (j, LIBMESH_DIM);
809 
810  // cheat by using finite difference approximations:
811  const Real eps = 1.e-6;
812  Point pp = p, pm = p;
813 
814  switch (j)
815  {
816  // d()/dxi
817  case 0:
818  {
819  pp(0) += eps;
820  pm(0) -= eps;
821  break;
822  }
823 
824  // d()/deta
825  case 1:
826  {
827  pp(1) += eps;
828  pm(1) -= eps;
829  break;
830  }
831 
832  // d()/dzeta
833  case 2:
834  {
835  pp(2) += eps;
836  pm(2) -= eps;
837  break;
838  }
839 
840  default:
841  libmesh_error_msg("Invalid derivative index j = " << j);
842  }
843 
844  return (shape_func(type, order, i, pp) -
845  shape_func(type, order, i, pm))/2./eps;
846 }
847 
848 
849 template <typename OutputShape>
850 OutputShape
852  const Order order,
853  const unsigned int i,
854  const unsigned int j,
855  const Point & p,
856  const bool add_p_level,
857  OutputShape(*deriv_func)
858  (const Elem *, const Order,
859  const unsigned int, const unsigned int,
860  const Point &, const bool))
861 {
862  // cheat by using finite difference approximations:
863  const Real eps = 1.e-5;
864  Point pp = p, pm = p;
865  unsigned int deriv_j = 0;
866 
867  switch (j)
868  {
869  // d^2() / dxi^2
870  case 0:
871  {
872  pp(0) += eps;
873  pm(0) -= eps;
874  deriv_j = 0;
875  break;
876  }
877 
878  // d^2() / dxi deta
879  case 1:
880  {
881  pp(1) += eps;
882  pm(1) -= eps;
883  deriv_j = 0;
884  break;
885  }
886 
887  // d^2() / deta^2
888  case 2:
889  {
890  pp(1) += eps;
891  pm(1) -= eps;
892  deriv_j = 1;
893  break;
894  }
895 
896  // d^2()/dxidzeta
897  case 3:
898  {
899  pp(2) += eps;
900  pm(2) -= eps;
901  deriv_j = 0;
902  break;
903  } // d^2()/deta^2
904 
905  // d^2()/detadzeta
906  case 4:
907  {
908  pp(2) += eps;
909  pm(2) -= eps;
910  deriv_j = 1;
911  break;
912  }
913 
914  // d^2()/dzeta^2
915  case 5:
916  {
917  pp(2) += eps;
918  pm(2) -= eps;
919  deriv_j = 2;
920  break;
921  }
922 
923  default:
924  libmesh_error_msg("Invalid shape function derivative j = " << j);
925  }
926 
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;
929 }
930 
931 
932 template <typename OutputShape>
933 OutputShape
935  const Order order,
936  const unsigned int i,
937  const unsigned int j,
938  const Point & p,
939  OutputShape(*deriv_func)
940  (const ElemType, const Order,
941  const unsigned int, const unsigned int,
942  const Point &))
943 {
944  // cheat by using finite difference approximations:
945  const Real eps = 1.e-5;
946  Point pp = p, pm = p;
947  unsigned int deriv_j = 0;
948 
949  switch (j)
950  {
951  // d^2() / dxi^2
952  case 0:
953  {
954  pp(0) += eps;
955  pm(0) -= eps;
956  deriv_j = 0;
957  break;
958  }
959 
960  // d^2() / dxi deta
961  case 1:
962  {
963  pp(1) += eps;
964  pm(1) -= eps;
965  deriv_j = 0;
966  break;
967  }
968 
969  // d^2() / deta^2
970  case 2:
971  {
972  pp(1) += eps;
973  pm(1) -= eps;
974  deriv_j = 1;
975  break;
976  }
977 
978  // d^2()/dxidzeta
979  case 3:
980  {
981  pp(2) += eps;
982  pm(2) -= eps;
983  deriv_j = 0;
984  break;
985  } // d^2()/deta^2
986 
987  // d^2()/detadzeta
988  case 4:
989  {
990  pp(2) += eps;
991  pm(2) -= eps;
992  deriv_j = 1;
993  break;
994  }
995 
996  // d^2()/dzeta^2
997  case 5:
998  {
999  pp(2) += eps;
1000  pm(2) -= eps;
1001  deriv_j = 2;
1002  break;
1003  }
1004 
1005  default:
1006  libmesh_error_msg("Invalid shape function derivative j = " << j);
1007  }
1008 
1009  return (deriv_func(type, order, i, deriv_j, pp) -
1010  deriv_func(type, order, i, deriv_j, pm))/2./eps;
1011 }
1012 
1013 
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)
1019 {
1020  const int extra_order = add_p_level * elem->p_level();
1021 
1022  const int dim = elem->dim();
1023 
1024  const unsigned int n_sf =
1025  FEInterface::n_shape_functions(underlying_fe_type, extra_order,
1026  elem);
1027 
1028  libmesh_assert_equal_to (n_sf, elem->n_nodes());
1029 
1030  std::vector<Real> node_weights(n_sf);
1031 
1032  const unsigned char datum_index = elem->mapping_data();
1033  for (unsigned int n=0; n<n_sf; n++)
1034  node_weights[n] =
1035  elem->node_ref(n).get_extra_datum<Real>(datum_index);
1036 
1037  const std::size_t n_p = p.size();
1038 
1039  shapes.resize(n_sf);
1040  for (unsigned int i=0; i != n_sf; ++i)
1041  {
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];
1048  }
1049 }
1050 
1051 
1053  const FEType fe_type,
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)
1058 {
1059  const int extra_order = add_p_level * elem->p_level();
1060  const unsigned int dim = elem->dim();
1061 
1062  const unsigned int n_sf =
1063  FEInterface::n_shape_functions(fe_type, extra_order, elem);
1064 
1065  libmesh_assert_equal_to (n_sf, elem->n_nodes());
1066 
1067  libmesh_assert_equal_to (dim, derivs.size());
1068  for (unsigned int d = 0; d != dim; ++d)
1069  derivs[d].resize(n_sf);
1070 
1071  std::vector<Real> node_weights(n_sf);
1072 
1073  const unsigned char datum_index = elem->mapping_data();
1074  for (unsigned int n=0; n<n_sf; n++)
1075  node_weights[n] =
1076  elem->node_ref(n).get_extra_datum<Real>(datum_index);
1077 
1078  const std::size_t n_p = p.size();
1079 
1080  shapes.resize(n_sf);
1081  for (unsigned int i=0; i != n_sf; ++i)
1082  shapes[i].resize(n_p, 0);
1083 
1084  FEInterface::all_shapes(dim, fe_type, elem, p, shapes, add_p_level);
1085 
1086  for (unsigned int i=0; i != n_sf; ++i)
1087  {
1088  auto & shapes_i = shapes[i];
1089 
1090  for (auto & s : shapes_i)
1091  s *= node_weights[i];
1092 
1093  for (unsigned int d = 0; d != dim; ++d)
1094  {
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];
1101  }
1102  }
1103 }
1104 
1105 
1107  const FEType underlying_fe_type,
1108  const unsigned int i,
1109  const Point & p,
1110  const bool add_p_level)
1111 {
1112  int extra_order = add_p_level * elem.p_level();
1113 
1114  const unsigned int n_sf =
1115  FEInterface::n_shape_functions(underlying_fe_type, extra_order, &elem);
1116 
1117  libmesh_assert_equal_to (n_sf, elem.n_nodes());
1118 
1119  std::vector<Real> node_weights(n_sf);
1120 
1121  const unsigned char datum_index = elem.mapping_data();
1122 
1123  Real weighted_shape_i = 0, weighted_sum = 0;
1124 
1125  for (unsigned int sf=0; sf<n_sf; sf++)
1126  {
1127  Real node_weight =
1128  elem.node_ref(sf).get_extra_datum<Real>(datum_index);
1129  Real weighted_shape = node_weight *
1130  FEInterface::shape(underlying_fe_type, extra_order, &elem, sf, p);
1131  weighted_sum += weighted_shape;
1132  if (sf == i)
1133  weighted_shape_i = weighted_shape;
1134  }
1135 
1136  return weighted_shape_i / weighted_sum;
1137 }
1138 
1139 
1141  const FEType underlying_fe_type,
1142  const unsigned int i,
1143  const unsigned int j,
1144  const Point & p,
1145  const bool add_p_level)
1146 {
1147  libmesh_assert_less(j, elem.dim());
1148 
1149  int extra_order = add_p_level * elem.p_level();
1150 
1151  const unsigned int n_sf =
1152  FEInterface::n_shape_functions(underlying_fe_type, extra_order, &elem);
1153 
1154  const unsigned int n_nodes = elem.n_nodes();
1155  libmesh_assert_equal_to (n_sf, n_nodes);
1156 
1157  std::vector<Real> node_weights(n_nodes);
1158 
1159  const unsigned char datum_index = elem.mapping_data();
1160  for (unsigned int n=0; n<n_nodes; n++)
1161  node_weights[n] =
1162  elem.node_ref(n).get_extra_datum<Real>(datum_index);
1163 
1164  Real weighted_shape_i = 0, weighted_sum = 0,
1165  weighted_grad_i = 0, weighted_grad_sum = 0;
1166 
1167  for (unsigned int sf=0; sf<n_sf; sf++)
1168  {
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;
1175  if (sf == i)
1176  {
1177  weighted_shape_i = weighted_shape;
1178  weighted_grad_i = weighted_grad;
1179  }
1180  }
1181 
1182  return (weighted_sum * weighted_grad_i - weighted_shape_i * weighted_grad_sum) /
1183  weighted_sum / weighted_sum;
1184 }
1185 
1186 
1187 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1188 
1190  const FEType underlying_fe_type,
1191  const unsigned int i,
1192  const unsigned int j,
1193  const Point & p,
1194  const bool add_p_level)
1195 {
1196  unsigned int j1, j2;
1197  switch (j)
1198  {
1199  case 0:
1200  // j = 0 ==> d^2 phi / dxi^2
1201  j1 = j2 = 0;
1202  break;
1203  case 1:
1204  // j = 1 ==> d^2 phi / dxi deta
1205  j1 = 0;
1206  j2 = 1;
1207  break;
1208  case 2:
1209  // j = 2 ==> d^2 phi / deta^2
1210  j1 = j2 = 1;
1211  break;
1212  case 3:
1213  // j = 3 ==> d^2 phi / dxi dzeta
1214  j1 = 0;
1215  j2 = 2;
1216  break;
1217  case 4:
1218  // j = 4 ==> d^2 phi / deta dzeta
1219  j1 = 1;
1220  j2 = 2;
1221  break;
1222  case 5:
1223  // j = 5 ==> d^2 phi / dzeta^2
1224  j1 = j2 = 2;
1225  break;
1226  default:
1227  libmesh_error();
1228  }
1229 
1230  int extra_order = add_p_level * elem.p_level();
1231 
1232  const unsigned int n_sf =
1233  FEInterface::n_shape_functions(underlying_fe_type, extra_order,
1234  &elem);
1235 
1236  const unsigned int n_nodes = elem.n_nodes();
1237  libmesh_assert_equal_to (n_sf, n_nodes);
1238 
1239  std::vector<Real> node_weights(n_nodes);
1240 
1241  const unsigned char datum_index = elem.mapping_data();
1242  for (unsigned int n=0; n<n_nodes; n++)
1243  node_weights[n] =
1244  elem.node_ref(n).get_extra_datum<Real>(datum_index);
1245 
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;
1250 
1251  for (unsigned int sf=0; sf<n_sf; sf++)
1252  {
1253  Real weighted_shape = node_weights[sf] *
1254  FEInterface::shape(underlying_fe_type, extra_order, &elem, sf,
1255  p);
1256  Real weighted_grada = node_weights[sf] *
1257  FEInterface::shape_deriv(underlying_fe_type, extra_order,
1258  &elem, sf, j1, p);
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;
1265  if (j1 != j2)
1266  {
1267  weighted_gradb = (j1 == j2) ? weighted_grada :
1268  node_weights[sf] *
1269  FEInterface::shape_deriv(underlying_fe_type, extra_order,
1270  &elem, sf, j2, p);
1271  weighted_grada_sum += weighted_grada;
1272  }
1273  weighted_hess_sum += weighted_hess;
1274  if (sf == i)
1275  {
1276  weighted_shape_i = weighted_shape;
1277  weighted_grada_i = weighted_grada;
1278  weighted_gradb_i = weighted_gradb;
1279  weighted_hess_i = weighted_hess;
1280  }
1281  }
1282 
1283  if (j1 == j2)
1284  weighted_gradb_sum = weighted_grada_sum;
1285 
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;
1290 }
1291 
1292 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1293 
1294 
1295 void rational_all_shapes (const Elem & elem,
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)
1300 {
1301  std::vector<std::vector<Real>> shapes;
1302 
1303  rational_fe_weighted_shapes(&elem, underlying_fe_type, shapes, p,
1304  add_p_level);
1305 
1306  std::vector<Real> shape_sums(p.size(), 0);
1307 
1308  for (auto i : index_range(v))
1309  {
1310  libmesh_assert_equal_to ( p.size(), shapes[i].size() );
1311  for (auto j : index_range(p))
1312  shape_sums[j] += shapes[i][j];
1313  }
1314 
1315  for (auto i : index_range(v))
1316  {
1317  libmesh_assert_equal_to ( p.size(), v[i].size() );
1318  for (auto j : index_range(v[i]))
1319  v[i][j] = shapes[i][j] / shape_sums[j];
1320  }
1321 }
1322 
1323 
1324 template <typename OutputShape>
1325 void rational_all_shape_derivs (const Elem & elem,
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)
1330 {
1331  const int my_dim = elem.dim();
1332 
1333  std::vector<std::vector<Real>> shapes;
1334  std::vector<std::vector<std::vector<Real>>> derivs(my_dim);
1335 
1336  rational_fe_weighted_shapes_derivs(&elem, underlying_fe_type,
1337  shapes, derivs, p, add_p_level);
1338 
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());
1343 
1344  for (auto i : index_range(shapes))
1345  {
1346  libmesh_assert_equal_to ( p.size(), shapes[i].size() );
1347  for (auto j : index_range(p))
1348  shape_sums[j] += shapes[i][j];
1349 
1350  for (int d=0; d != my_dim; ++d)
1351  for (auto j : index_range(p))
1352  shape_deriv_sums[d][j] += derivs[d][i][j];
1353  }
1354 
1355  for (int d=0; d != my_dim; ++d)
1356  {
1357  auto & comps_d = *comps[d];
1358  libmesh_assert_equal_to(comps_d.size(), elem.n_nodes());
1359 
1360  for (auto i : index_range(comps_d))
1361  {
1362  auto & comps_di = comps_d[i];
1363  auto & derivs_di = derivs[d][i];
1364 
1365  for (auto j : index_range(comps_di))
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];
1369  }
1370  }
1371 }
1372 
1373 
1374 
1375 template
1376 Real fe_fdm_deriv<Real>(const Elem *, const Order, const unsigned int,
1377  const unsigned int, const Point &, const bool,
1378  Real(*shape_func)
1379  (const Elem *, const Order, const unsigned int,
1380  const Point &, const bool));
1381 
1382 template
1383 Real fe_fdm_deriv<Real>(const ElemType, const Order, const unsigned int,
1384  const unsigned int, const Point &,
1385  Real(*shape_func)
1386  (const ElemType, const Order, const unsigned int,
1387  const Point &));
1388 
1389 template
1391 fe_fdm_deriv<RealGradient>(const Elem *, const Order, const unsigned int,
1392  const unsigned int, const Point &, const bool,
1393  RealGradient(*shape_func)
1394  (const Elem *, const Order, const unsigned int,
1395  const Point &, const bool));
1396 
1397 template
1398 Real
1399 fe_fdm_second_deriv<Real>(const ElemType, const Order, const unsigned int,
1400  const unsigned int, const Point &,
1401  Real(*shape_func)
1402  (const ElemType, const Order, const unsigned int,
1403  const unsigned int, const Point &));
1404 
1405 template
1406 Real
1407 fe_fdm_second_deriv<Real>(const Elem *, const Order, const unsigned int,
1408  const unsigned int, const Point &, const bool,
1409  Real(*shape_func)
1410  (const Elem *, const Order, const unsigned int,
1411  const unsigned int, const Point &, const bool));
1412 
1413 template
1415 fe_fdm_second_deriv<RealGradient>(const Elem *, const Order, const unsigned int,
1416  const unsigned int, const Point &, const bool,
1417  RealGradient(*shape_func)
1418  (const Elem *, const Order, const unsigned int,
1419  const unsigned int, const Point &, const bool));
1420 
1421 
1422 
1423 
1424 //--------------------------------------------------------------
1425 // Explicit instantiations using macro from fe_macro.h
1426 
1427 INSTANTIATE_FE(0);
1428 
1429 INSTANTIATE_FE(1);
1430 
1431 INSTANTIATE_FE(2);
1432 
1433 INSTANTIATE_FE(3);
1434 
1436 
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);
1438 } // namespace libMesh
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
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.
Definition: quadrature.C:64
unsigned char mapping_data() const
Definition: elem.h:2973
ElemType
Defines an enum for geometric element types.
RealVectorValue RealGradient
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
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)
Definition: fe.C:1052
INSTANTIATE_SUBDIVISION_FE
Definition: fe.C:1435
unsigned int dim
INSTANTIATE_FE(3)
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
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.
Definition: fe.C:48
Order default_quadrature_order() const
Definition: fe_type.h:357
unsigned int p_level() const
Definition: elem.h:2945
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&#39;t be...
Definition: fe.C:744
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const =0
A specific instantiation of the FEBase class.
Definition: fe.h:127
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
Definition: tecplot_io.C:67
ElemMappingType mapping_type() const
Definition: elem.h:2957
const Node & node_ref(const unsigned int i) const
Definition: elem.h:2353
virtual unsigned int n_nodes() const =0
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:436
libmesh_assert(ctx)
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)
Definition: fe.C:1325
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))
Definition: fe.C:851
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)
Definition: fe.C:1106
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)
Definition: fe.C:1189
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...
Definition: int_range.h:134
This class implements specific orders of Gauss quadrature.
IntRange< unsigned short > node_index_range() const
Definition: elem.h:2492
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...
Definition: dof_object.h:1139
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)
Definition: fe.C:1295
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.
Definition: fe.C:1014
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
const Point & point(const unsigned int i) const
Definition: elem.h:2277
The QBase class provides the basic functionality from which various quadrature rules can be derived...
Definition: quadrature.h:53
void ErrorVector unsigned int
Definition: adjoints_ex3.C:360
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:111
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.
Definition: fe.h:51
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)
Definition: fe.C:1140