libMesh
fe_interface.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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/fe_interface.h"
22 
23 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
24 #include "libmesh/fe_interface_macros.h"
25 #include "libmesh/inf_fe.h"
26 #endif
27 
28 #include "libmesh/elem.h"
29 #include "libmesh/fe.h"
30 #include "libmesh/fe_compute_data.h"
31 #include "libmesh/dof_map.h"
32 #include "libmesh/enum_fe_family.h"
33 #include "libmesh/enum_order.h"
34 #include "libmesh/enum_elem_type.h"
35 #include "libmesh/enum_to_string.h"
36 
37 namespace libMesh
38 {
39 
40 //------------------------------------------------------------
41 //FEInterface class members
43 {
44  libmesh_error_msg("ERROR: Do not define an object of this type.");
45 }
46 
47 
48 
49 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
50 
51 bool
53 {
54  return false;
55 }
56 
57 #else
58 
59 bool
61 {
62 
63  switch (et)
64  {
65  case INFEDGE2:
66  case INFQUAD4:
67  case INFQUAD6:
68  case INFHEX8:
69  case INFHEX16:
70  case INFHEX18:
71  case INFPRISM6:
72  case INFPRISM12:
73  return true;
74 
75  default:
76  return false;
77  }
78 }
79 
80 #endif //ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
81 
82 #define fe_family_case_func(family, dim, func_and_args, prefix, suffix) \
83  case family: \
84  prefix FE<dim,family>::func_and_args suffix \
85 
86 #define fe_family_case(family) \
87  case family:
88 
89 #define fe_family_scalar_case_func(dim, func_and_args, prefix, suffix) \
90  fe_family_case_func(CLOUGH, dim, func_and_args, prefix, suffix) \
91  fe_family_case_func(HERMITE, dim, func_and_args, prefix, suffix) \
92  fe_family_case_func(HIERARCHIC, dim, func_and_args, prefix, suffix) \
93  fe_family_case_func(L2_HIERARCHIC, dim, func_and_args, prefix, suffix) \
94  fe_family_case_func(SIDE_HIERARCHIC, dim, func_and_args, prefix, suffix) \
95  fe_family_case_func(LAGRANGE, dim, func_and_args, prefix, suffix) \
96  fe_family_case_func(L2_LAGRANGE, dim, func_and_args, prefix, suffix) \
97  fe_family_case_func(MONOMIAL, dim, func_and_args, prefix, suffix) \
98  fe_family_case_func(SCALAR, dim, func_and_args, prefix, suffix) \
99  fe_family_case_func(XYZ, dim, func_and_args, prefix, suffix) \
100  fe_family_case_func(SUBDIVISION, 2, func_and_args, \
101  libmesh_assert_equal_to (dim, 2); prefix, suffix)
102 
103 #define fe_family_scalar_case() \
104  fe_family_case(CLOUGH) \
105  fe_family_case(HERMITE) \
106  fe_family_case(HIERARCHIC) \
107  fe_family_case(L2_HIERARCHIC) \
108  fe_family_case(SIDE_HIERARCHIC) \
109  fe_family_case(LAGRANGE) \
110  fe_family_case(L2_LAGRANGE) \
111  fe_family_case(MONOMIAL) \
112  fe_family_case(SCALAR) \
113  fe_family_case(XYZ) \
114  fe_family_case(SUBDIVISION) \
115 
116 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
117 
118 #define fe_family_horder_case_func(dim, func_and_args, prefix, suffix) \
119  fe_family_case_func(BERNSTEIN, dim, func_and_args, prefix, suffix) \
120  fe_family_case_func(SZABAB, dim, func_and_args, prefix, suffix) \
121  fe_family_case_func(RATIONAL_BERNSTEIN, dim, func_and_args, prefix, suffix) \
122 
123 #define fe_family_horder_case() \
124  fe_family_case(BERNSTEIN) \
125  fe_family_case(SZABAB) \
126  fe_family_case(RATIONAL_BERNSTEIN) \
127 
128 #else
129 
130 #define fe_family_horder_case_func(dim, func_and_args, prefix, suffix)
131 #define fe_family_horder_case()
132 
133 #endif
134 
135 #define fe_family_vector_case_func(dim, func_and_args, prefix, suffix) \
136  fe_family_case_func(HIERARCHIC_VEC, dim, func_and_args, prefix, suffix) \
137  fe_family_case_func(L2_HIERARCHIC_VEC, dim, func_and_args, prefix, suffix) \
138  fe_family_case_func(LAGRANGE_VEC, dim, func_and_args, prefix, suffix) \
139  fe_family_case_func(L2_LAGRANGE_VEC, dim, func_and_args, prefix, suffix) \
140  fe_family_case_func(MONOMIAL_VEC, dim, func_and_args, prefix, suffix) \
141  fe_family_case_func(NEDELEC_ONE, dim, func_and_args, prefix, suffix) \
142  fe_family_case_func(RAVIART_THOMAS, dim, func_and_args, prefix, suffix) \
143  fe_family_case_func(L2_RAVIART_THOMAS, dim, func_and_args, prefix, suffix) \
144 
145 #define fe_family_vector_case() \
146  fe_family_case(HIERARCHIC_VEC) \
147  fe_family_case(L2_HIERARCHIC_VEC) \
148  fe_family_case(LAGRANGE_VEC) \
149  fe_family_case(L2_LAGRANGE_VEC) \
150  fe_family_case(MONOMIAL_VEC) \
151  fe_family_case(NEDELEC_ONE) \
152  fe_family_case(RAVIART_THOMAS) \
153  fe_family_case(L2_RAVIART_THOMAS) \
154 
155 #define fe_family_switch(dim, func_and_args, prefix, suffix) \
156  do { \
157  switch (fe_t.family) \
158  { \
159  fe_family_scalar_case_func(dim, func_and_args, prefix, suffix) \
160  fe_family_horder_case_func(dim, func_and_args, prefix, suffix) \
161  default: \
162  libmesh_error_msg("Unsupported family = " << Utility::enum_to_string(fe_t.family)); \
163  } \
164  } while (0)
165 
166 #define fe_family_with_vec_switch(dim, func_and_args, prefix, suffix) \
167  do { \
168  switch (fe_t.family) \
169  { \
170  fe_family_scalar_case_func(dim, func_and_args, prefix, suffix) \
171  fe_family_horder_case_func(dim, func_and_args, prefix, suffix) \
172  fe_family_vector_case_func(dim, func_and_args, prefix, suffix) \
173  default: \
174  libmesh_error_msg("Unsupported family = " << Utility::enum_to_string(fe_t.family)); \
175  } \
176  } while (0)
177 
178 #define fe_scalar_vec_error_switch(dim, func_and_args, prefix, suffix) \
179  do { \
180  switch (fe_t.family) \
181  { \
182  fe_family_scalar_case_func(dim, func_and_args, prefix, suffix) \
183  fe_family_horder_case_func(dim, func_and_args, prefix, suffix) \
184  fe_family_vector_case() \
185  libmesh_error_msg("Error: Can only request scalar valued elements for Real FEInterface::func_and_args"); \
186  default: \
187  libmesh_error_msg("Unsupported family = " << Utility::enum_to_string(fe_t.family)); \
188  } \
189  } while (0)
190 
191 
192 #define fe_vector_scalar_error_switch(dim, func_and_args, prefix, suffix) \
193  do { \
194  switch (fe_t.family) \
195  { \
196  fe_family_vector_case_func(dim, func_and_args, prefix, suffix) \
197  fe_family_scalar_case() \
198  fe_family_horder_case() \
199  libmesh_error_msg("Error: Can only request vector valued elements for RealGradient FEInterface::func_and_args"); \
200  default: \
201  libmesh_error_msg("Unsupported family = " << Utility::enum_to_string(fe_t.family)); \
202  } \
203  } while (0)
204 
205 
206 #define fe_switch(func_and_args) \
207  do { \
208  switch (dim) \
209  { \
210  /* 0D */ \
211  case 0: \
212  fe_family_switch (0, func_and_args, return, ;); \
213  /* 1D */ \
214  case 1: \
215  fe_family_switch (1, func_and_args, return, ;); \
216  /* 2D */ \
217  case 2: \
218  fe_family_switch (2, func_and_args, return, ;); \
219  /* 3D */ \
220  case 3: \
221  fe_family_switch (3, func_and_args, return, ;); \
222  default: \
223  libmesh_error_msg("Invalid dim = " << dim); \
224  } \
225  } while (0)
226 
227 #define fe_with_vec_switch(func_and_args) \
228  do { \
229  switch (dim) \
230  { \
231  /* 0D */ \
232  case 0: \
233  fe_family_with_vec_switch (0, func_and_args, return, ;); \
234  /* 1D */ \
235  case 1: \
236  fe_family_with_vec_switch (1, func_and_args, return, ;); \
237  /* 2D */ \
238  case 2: \
239  fe_family_with_vec_switch (2, func_and_args, return, ;); \
240  /* 3D */ \
241  case 3: \
242  fe_family_with_vec_switch (3, func_and_args, return, ;); \
243  default: \
244  libmesh_error_msg("Invalid dim = " << dim); \
245  } \
246  } while (0)
247 
248 #define void_fe_with_vec_switch(func_and_args) \
249  do { \
250  switch (dim) \
251  { \
252  /* 0D */ \
253  case 0: \
254  fe_family_with_vec_switch (0, func_and_args, , ; return;); \
255  /* 1D */ \
256  case 1: \
257  fe_family_with_vec_switch (1, func_and_args, , ; return;); \
258  /* 2D */ \
259  case 2: \
260  fe_family_with_vec_switch (2, func_and_args, , ; return;); \
261  /* 3D */ \
262  case 3: \
263  fe_family_with_vec_switch (3, func_and_args, , ; return;); \
264  default: \
265  libmesh_error_msg("Invalid dim = " << dim); \
266  } \
267  } while (0)
268 
269 
270 
271 #ifdef LIBMESH_ENABLE_DEPRECATED
272 unsigned int FEInterface::n_shape_functions(const unsigned int dim,
273  const FEType & fe_t,
274  const ElemType t)
275 {
276  libmesh_deprecated();
277 
278 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
279  /*
280  * Since the FEType, stored in DofMap/(some System child), has to
281  * be the _same_ for InfFE and FE, we have to catch calls
282  * to infinite elements through the element type.
283  */
284 
285  if (is_InfFE_elem(t))
286  return ifem_n_shape_functions(dim, fe_t, t);
287 
288 #endif
289 
290  const Order o = fe_t.order;
291 
292  fe_with_vec_switch(n_shape_functions(t, o));
293 }
294 #endif // LIBMESH_ENABLE_DEPRECATED
295 
296 
297 
298 unsigned int
300  const Elem * elem,
301  const bool add_p_level)
302 {
303  // dim is required by the fe_with_vec_switch macro
304  auto dim = elem->dim();
305 
306 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
307  /*
308  * Since the FEType, stored in DofMap/(some System child), has to
309  * be the _same_ for InfFE and FE, we have to catch calls
310  * to infinite elements through the element type.
311  */
312 
313  if (is_InfFE_elem(elem->type()))
314  return ifem_n_shape_functions(fe_t, elem);
315 
316 #endif
317 
318  // Account for Elem::p_level() when computing total_order
319  auto total_order = fe_t.order + add_p_level*elem->p_level();
320 
321  fe_with_vec_switch(n_dofs(elem, total_order));
322 }
323 
324 
325 
326 unsigned int
328  const int extra_order,
329  const Elem * elem)
330 {
331  // dim is required by the fe_with_vec_switch macro
332  auto dim = elem->dim();
333 
334 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
335  /*
336  * Since the FEType, stored in DofMap/(some System child), has to
337  * be the _same_ for InfFE and FE, we have to catch calls
338  * to infinite elements through the element type.
339  */
340 
341  if (is_InfFE_elem(elem->type()))
342  return ifem_n_shape_functions(fe_t, elem);
343 
344 #endif
345 
346  // Ignore Elem::p_level() and instead use extra_order to compute total_order.
347  auto total_order = fe_t.order + extra_order;
348 
349  fe_with_vec_switch(n_dofs(elem, total_order));
350 }
351 
352 
353 
354 #ifdef LIBMESH_ENABLE_DEPRECATED
355 unsigned int FEInterface::n_dofs(const unsigned int dim,
356  const FEType & fe_t,
357  const ElemType t)
358 {
359  libmesh_deprecated();
360 
361 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
362 
363  if (is_InfFE_elem(t))
364  return ifem_n_dofs(dim, fe_t, t);
365 
366 #endif
367 
368  const Order o = fe_t.order;
369 
370  fe_with_vec_switch(n_dofs(t, o));
371 }
372 
373 
374 
375 
376 unsigned int
377 FEInterface::n_dofs (const unsigned int dim,
378  const FEType & fe_t,
379  const Elem * elem)
380 {
381  libmesh_deprecated();
382 
383  fe_with_vec_switch(n_dofs(elem, fe_t.order + elem->p_level()));
384 }
385 #endif // LIBMESH_ENABLE_DEPRECATED
386 
387 
388 
389 unsigned int
391  const Elem * elem,
392  const bool add_p_level)
393 {
394  // dim is required by the fe_with_vec_switch macro
395  auto dim = elem->dim();
396 
397 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
398 
399  // InfElems currently don't support p_level()
400  if (is_InfFE_elem(elem->type()))
401  return ifem_n_dofs(fe_t, elem);
402 
403 #endif
404 
405  // Account for Elem::p_level() when computing total_order
406  fe_with_vec_switch(n_dofs(elem, fe_t.order + add_p_level*elem->p_level()));
407 }
408 
409 
410 
411 unsigned int
413  int extra_order,
414  const Elem * elem)
415 {
416  // dim is required by the fe_with_vec_switch macro
417  auto dim = elem->dim();
418 
419 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
420 
421  // InfElems currently don't support p_level()
422  if (is_InfFE_elem(elem->type()))
423  return ifem_n_dofs(fe_t, elem);
424 
425 #endif
426 
427  // Elem::p_level() is ignored, extra_order is used instead.
428  auto total_order = fe_t.order + extra_order;
429 
430  fe_with_vec_switch(n_dofs(elem, total_order));
431 }
432 
433 
434 
435 #ifdef LIBMESH_ENABLE_DEPRECATED
436 unsigned int FEInterface::n_dofs_at_node(const unsigned int dim,
437  const FEType & fe_t,
438  const ElemType t,
439  const unsigned int n)
440 {
441  libmesh_deprecated();
442 
443 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
444 
445  if (is_InfFE_elem(t))
446  return ifem_n_dofs_at_node(dim, fe_t, t, n);
447 
448 #endif
449 
450  const Order o = fe_t.order;
451 
452  fe_with_vec_switch(n_dofs_at_node(t, o, n));
453 }
454 
455 
456 
459  const FEType & fe_t)
460 {
461  libmesh_deprecated();
462 
463  fe_with_vec_switch(n_dofs_at_node);
464 }
465 #endif // LIBMESH_ENABLE_DEPRECATED
466 
467 
468 
471  const Elem * elem)
472 {
473  // dim is required by the fe_with_vec_switch macro
474  auto dim = elem->dim();
475 
476  fe_with_vec_switch(n_dofs_at_node);
477 }
478 
479 
480 
481 unsigned int
483  const Elem * elem,
484  const unsigned int n,
485  const bool add_p_level)
486 {
487  // dim is required by the fe_with_vec_switch macro
488  auto dim = elem->dim();
489 
490 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
491 
492  if (is_InfFE_elem(elem->type()))
493  return ifem_n_dofs_at_node(fe_t, elem, n);
494 
495 #endif
496 
497  // Account for Elem::p_level() when computing total_order
498  auto total_order = fe_t.order + add_p_level*elem->p_level();
499 
500  fe_with_vec_switch(n_dofs_at_node(*elem, total_order, n));
501 }
502 
503 
504 
505 unsigned int
507  const int extra_order,
508  const Elem * elem,
509  const unsigned int n)
510 {
511  // dim is required by the fe_with_vec_switch macro
512  auto dim = elem->dim();
513 
514 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
515 
516  if (is_InfFE_elem(elem->type()))
517  return ifem_n_dofs_at_node(fe_t, elem, n);
518 
519 #endif
520 
521  // Ignore Elem::p_level() and instead use extra_order to compute total_order.
522  auto total_order = fe_t.order + extra_order;
523 
524  fe_with_vec_switch(n_dofs_at_node(*elem, total_order, n));
525 }
526 
527 
528 
529 #ifdef LIBMESH_ENABLE_DEPRECATED
530 unsigned int FEInterface::n_dofs_per_elem(const unsigned int dim,
531  const FEType & fe_t,
532  const ElemType t)
533 {
534  libmesh_deprecated();
535 
536 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
537 
538  if (is_InfFE_elem(t))
539  return ifem_n_dofs_per_elem(dim, fe_t, t);
540 
541 #endif
542 
543  const Order o = fe_t.order;
544 
545  fe_with_vec_switch(n_dofs_per_elem(t, o));
546 }
547 #endif // LIBMESH_ENABLE_DEPRECATED
548 
549 
550 
551 unsigned int
553  const Elem * elem,
554  const bool add_p_level)
555 {
556  // dim is required by the fe_with_vec_switch macro
557  auto dim = elem->dim();
558 
559 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
560 
561  if (is_InfFE_elem(elem->type()))
562  return ifem_n_dofs_per_elem(fe_t, elem);
563 
564 #endif
565 
566  // Account for Elem::p_level() when computing total_order
567  auto total_order = fe_t.order + add_p_level*elem->p_level();
568 
569  fe_with_vec_switch(n_dofs_per_elem(*elem, total_order));
570 }
571 
572 
573 
574 unsigned int
576  const int extra_order,
577  const Elem * elem)
578 {
579  // dim is required by the fe_with_vec_switch macro
580  auto dim = elem->dim();
581 
582 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
583 
584  if (is_InfFE_elem(elem->type()))
585  return ifem_n_dofs_per_elem(fe_t, elem);
586 
587 #endif
588 
589  // Ignore Elem::p_level() and instead use extra_order to compute total_order.
590  auto total_order = fe_t.order + extra_order;
591 
592  fe_with_vec_switch(n_dofs_per_elem(*elem, total_order));
593 }
594 
595 
596 
597 void FEInterface::dofs_on_side(const Elem * const elem,
598  const unsigned int dim,
599  const FEType & fe_t,
600  unsigned int s,
601  std::vector<unsigned int> & di,
602  const bool add_p_level)
603 {
604  const Order o = fe_t.order;
605 
606  void_fe_with_vec_switch(dofs_on_side(elem, o, s, di, add_p_level));
607 }
608 
609 
610 
611 void FEInterface::dofs_on_edge(const Elem * const elem,
612  const unsigned int dim,
613  const FEType & fe_t,
614  unsigned int e,
615  std::vector<unsigned int> & di,
616  const bool add_p_level)
617 {
618  const Order o = fe_t.order;
619 
620  void_fe_with_vec_switch(dofs_on_edge(elem, o, e, di, add_p_level));
621 }
622 
623 
624 
625 void FEInterface::nodal_soln(const unsigned int dim,
626  const FEType & fe_t,
627  const Elem * elem,
628  const std::vector<Number> & elem_soln,
629  std::vector<Number> & nodal_soln,
630  const bool add_p_level,
631  const unsigned int vdim)
632 {
633 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
634 
635  if (is_InfFE_elem(elem->type()))
636  {
637  ifem_nodal_soln(dim, fe_t, elem, elem_soln, nodal_soln);
638  return;
639  }
640 
641 #endif
642 
643  const Order order = fe_t.order;
644 
645  void_fe_with_vec_switch(nodal_soln(elem, order, elem_soln, nodal_soln, add_p_level, vdim));
646 }
647 
648 
649 
651  const Elem * elem,
652  const unsigned int side,
653  const std::vector<Number> & elem_soln,
654  std::vector<Number> & nodal_soln,
655  const bool add_p_level,
656  const unsigned int vdim)
657 {
658 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
659 
660  if (is_InfFE_elem(elem->type()))
661  {
662  libmesh_not_implemented();
663  return;
664  }
665 
666 #endif
667 
668  const Order order = fe_t.order;
669  const unsigned int dim = elem->dim();
670 
671  void_fe_with_vec_switch(side_nodal_soln(elem, order, side, elem_soln, nodal_soln, add_p_level, vdim));
672 }
673 
674 
675 
676 #ifdef LIBMESH_ENABLE_DEPRECATED
678  const FEType & fe_t,
679  const Elem * elem,
680  const Point & p)
681 {
682  libmesh_deprecated();
683 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
684  if (is_InfFE_elem(elem->type()))
685  return ifem_map(dim, fe_t, elem, p);
686 #endif
687  fe_with_vec_switch(map(elem, p));
688 }
689 
690 
691 
692 Point FEInterface::inverse_map (const unsigned int dim,
693  const FEType & fe_t,
694  const Elem * elem,
695  const Point & p,
696  const Real tolerance,
697  const bool secure)
698 {
699  libmesh_deprecated();
700 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
701 
702  if (is_InfFE_elem(elem->type()))
703  return ifem_inverse_map(dim, fe_t, elem, p,tolerance, secure);
704 
705 #endif
706 
707  fe_with_vec_switch(inverse_map(elem, p, tolerance, secure));
708 }
709 
710 
711 
712 void FEInterface::inverse_map (const unsigned int dim,
713  const FEType & fe_t,
714  const Elem * elem,
715  const std::vector<Point> & physical_points,
716  std::vector<Point> & reference_points,
717  const Real tolerance,
718  const bool secure)
719 {
720  libmesh_deprecated();
721 
722  const std::size_t n_pts = physical_points.size();
723 
724  // Resize the vector
725  reference_points.resize(n_pts);
726 
727  if (n_pts == 0)
728  {
729  libMesh::err << "WARNING: empty vector physical_points!"
730  << std::endl;
731  libmesh_here();
732  return;
733  }
734 
735 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
736 
737  if (is_InfFE_elem(elem->type()))
738  {
739  ifem_inverse_map(dim, fe_t, elem, physical_points, reference_points, tolerance, secure);
740  return;
741  // libmesh_not_implemented();
742  }
743 
744 #endif
745 
746  void_fe_with_vec_switch(inverse_map(elem, physical_points, reference_points, tolerance, secure));
747 }
748 
749 
750 
752  const ElemType t,
753  const Real eps)
754 {
755  return FEBase::on_reference_element(p,t,eps);
756 }
757 
758 
759 
760 Real FEInterface::shape(const unsigned int dim,
761  const FEType & fe_t,
762  const ElemType t,
763  const unsigned int i,
764  const Point & p)
765 {
766  libmesh_deprecated();
767 
768 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
769 
770  if (is_InfFE_elem(t))
771  return ifem_shape(dim, fe_t, t, i, p);
772 
773 #endif
774 
775  const Order o = fe_t.order;
776 
777  fe_switch(shape(t,o,i,p));
778 }
779 
780 Real FEInterface::shape(const unsigned int dim,
781  const FEType & fe_t,
782  const Elem * elem,
783  const unsigned int i,
784  const Point & p)
785 {
786  libmesh_deprecated();
787 
788 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
789 
790  if (elem && is_InfFE_elem(elem->type()))
791  return ifem_shape(fe_t, elem, i, p);
792 
793 #endif
794 
795  const Order o = fe_t.order;
796 
797  fe_switch(shape(elem,o,i,p));
798 }
799 #endif // LIBMESH_ENABLE_DEPRECATED
800 
801 
802 
803 Real
805  const Elem * elem,
806  const unsigned int i,
807  const Point & p,
808  const bool add_p_level)
809 {
810  // dim is required by the fe_switch macro
811  auto dim = elem->dim();
812 
813 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
814 
815  if (elem && is_InfFE_elem(elem->type()))
816  return ifem_shape(fe_t, elem, i, p);
817 
818 #endif
819 
820  // We are calling
821  //
822  // FE<X,Y>::shape(Elem *, Order, unsigned, Point, true)
823  //
824  // See fe.h for more details.
825  fe_switch(shape(elem, fe_t.order, i, p, add_p_level));
826 }
827 
828 
829 
830 Real
832  int extra_order,
833  const Elem * elem,
834  const unsigned int i,
835  const Point & p)
836 {
837  // dim is required by the fe_switch macro
838  auto dim = elem->dim();
839 
840 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
841 
842  if (elem && is_InfFE_elem(elem->type()))
843  return ifem_shape(fe_t, elem, i, p);
844 
845 #endif
846 
847  // We are calling
848  //
849  // FE<X,Y>::shape(Elem *, Order, unsigned, Point, false)
850  //
851  // with the last parameter set to "false" so that the
852  // Elem::p_level() is not used internally and the "total_order" that
853  // we compute is used instead. See fe.h for more details.
854  auto total_order = fe_t.order + extra_order;
855 
856  fe_switch(shape(elem, total_order, i, p, false));
857 }
858 
859 
860 
861 #ifdef LIBMESH_ENABLE_DEPRECATED
862 template<>
863 void FEInterface::shape<Real>(const unsigned int dim,
864  const FEType & fe_t,
865  const ElemType t,
866  const unsigned int i,
867  const Point & p,
868  Real & phi)
869 {
870  libmesh_deprecated();
871 
872 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
873 
874  if (is_InfFE_elem(t))
875  {
876  phi = ifem_shape(dim, fe_t, t, i, p);
877  return;
878  }
879 
880 #endif
881 
882  const Order o = fe_t.order;
883 
884  switch(dim)
885  {
886  case 0:
887  fe_scalar_vec_error_switch(0, shape(t,o,i,p), phi = , ; break;);
888  break;
889  case 1:
890  fe_scalar_vec_error_switch(1, shape(t,o,i,p), phi = , ; break;);
891  break;
892  case 2:
893  fe_scalar_vec_error_switch(2, shape(t,o,i,p), phi = , ; break;);
894  break;
895  case 3:
896  fe_scalar_vec_error_switch(3, shape(t,o,i,p), phi = , ; break;);
897  break;
898  default:
899  libmesh_error_msg("Invalid dimension = " << dim);
900  }
901 
902  return;
903 }
904 
905 
906 
907 template<>
908 void FEInterface::shape<Real>(const unsigned int dim,
909  const FEType & fe_t,
910  const Elem * elem,
911  const unsigned int i,
912  const Point & p,
913  Real & phi)
914 {
915  libmesh_deprecated();
916 
917 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
918 
919  if (elem && is_InfFE_elem(elem->type()))
920  {
921  phi = ifem_shape(fe_t, elem, i, p);
922  return;
923  }
924 #endif
925 
926  const Order o = fe_t.order;
927 
928  switch(dim)
929  {
930  case 0:
931  fe_scalar_vec_error_switch(0, shape(elem,o,i,p), phi = , ; break;);
932  break;
933  case 1:
934  fe_scalar_vec_error_switch(1, shape(elem,o,i,p), phi = , ; break;);
935  break;
936  case 2:
937  fe_scalar_vec_error_switch(2, shape(elem,o,i,p), phi = , ; break;);
938  break;
939  case 3:
940  fe_scalar_vec_error_switch(3, shape(elem,o,i,p), phi = , ; break;);
941  break;
942  default:
943  libmesh_error_msg("Invalid dimension = " << dim);
944  }
945 
946  return;
947 }
948 #endif // LIBMESH_ENABLE_DEPRECATED
949 
950 
951 
952 template<>
953 void FEInterface::shape<Real>(const FEType & fe_t,
954  const Elem * elem,
955  const unsigned int i,
956  const Point & p,
957  Real & phi)
958 {
959  // dim is required by the fe_switch macro
960  auto dim = elem->dim();
961 
962 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
963 
964  if (is_InfFE_elem(elem->type()))
965  {
966  phi = ifem_shape(fe_t, elem, i, p);
967  return;
968  }
969 
970 #endif
971 
972  // Below we call FE<X,Y>::shape(Elem *, Order, unsigned, Point, true)
973  // so that the Elem::p_level() is accounted for internally.
974  switch(dim)
975  {
976  case 0:
977  fe_scalar_vec_error_switch(0, shape(elem, fe_t.order, i, p, true), phi = , ; break;);
978  break;
979  case 1:
980  fe_scalar_vec_error_switch(1, shape(elem, fe_t.order, i, p, true), phi = , ; break;);
981  break;
982  case 2:
983  fe_scalar_vec_error_switch(2, shape(elem, fe_t.order, i, p, true), phi = , ; break;);
984  break;
985  case 3:
986  fe_scalar_vec_error_switch(3, shape(elem, fe_t.order, i, p, true), phi = , ; break;);
987  break;
988  default:
989  libmesh_error_msg("Invalid dimension = " << dim);
990  }
991 }
992 
993 
994 
995 template<>
996 void FEInterface::shape<Real>(const FEType & fe_t,
997  int extra_order,
998  const Elem * elem,
999  const unsigned int i,
1000  const Point & p,
1001  Real & phi)
1002 {
1003  // dim is required by the fe_switch macro
1004  auto dim = elem->dim();
1005 
1006 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1007 
1008  if (is_InfFE_elem(elem->type()))
1009  {
1010  phi = ifem_shape(fe_t, elem, i, p);
1011  return;
1012  }
1013 
1014 #endif
1015 
1016  // Ignore Elem::p_level() and instead use extra_order to compute total_order
1017  auto total_order = fe_t.order + extra_order;
1018 
1019  // Below we call
1020  //
1021  // FE<X,Y>::shape(Elem *, Order, unsigned, Point, false)
1022  //
1023  // so that the Elem::p_level() is ignored and the total_order that
1024  // we compute is used instead.
1025  switch(dim)
1026  {
1027  case 0:
1028  fe_scalar_vec_error_switch(0, shape(elem, total_order, i, p, false), phi = , ; break;);
1029  break;
1030  case 1:
1031  fe_scalar_vec_error_switch(1, shape(elem, total_order, i, p, false), phi = , ; break;);
1032  break;
1033  case 2:
1034  fe_scalar_vec_error_switch(2, shape(elem, total_order, i, p, false), phi = , ; break;);
1035  break;
1036  case 3:
1037  fe_scalar_vec_error_switch(3, shape(elem, total_order, i, p, false), phi = , ; break;);
1038  break;
1039  default:
1040  libmesh_error_msg("Invalid dimension = " << dim);
1041  }
1042 }
1043 
1044 
1045 
1046 template<>
1047 void FEInterface::shapes<Real>(const unsigned int dim,
1048  const FEType & fe_t,
1049  const Elem * elem,
1050  const unsigned int i,
1051  const std::vector<Point> & p,
1052  std::vector<Real> & phi,
1053  const bool add_p_level)
1054 {
1055 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1056 
1057  if (elem && is_InfFE_elem(elem->type()))
1058  {
1059  FEType elevated = fe_t;
1060  elevated.order = fe_t.order + add_p_level * elem->p_level();
1061  for (auto qpi : index_range(p))
1062  phi[qpi] = ifem_shape(elevated, elem, i, p[qpi]);
1063  return;
1064  }
1065 #endif
1066 
1067  const Order o = fe_t.order;
1068 
1069  switch(dim)
1070  {
1071  case 0:
1072  fe_scalar_vec_error_switch(0, shapes(elem,o,i,p,phi,add_p_level), , ; return;);
1073  break;
1074  case 1:
1075  fe_scalar_vec_error_switch(1, shapes(elem,o,i,p,phi,add_p_level), , ; return;);
1076  break;
1077  case 2:
1078  fe_scalar_vec_error_switch(2, shapes(elem,o,i,p,phi,add_p_level), , ; return;);
1079  break;
1080  case 3:
1081  fe_scalar_vec_error_switch(3, shapes(elem,o,i,p,phi,add_p_level), , ; return;);
1082  break;
1083  default:
1084  libmesh_error_msg("Invalid dimension = " << dim);
1085  }
1086 
1087  return;
1088 }
1089 
1090 
1091 template<>
1092 void FEInterface::all_shapes<Real>(const unsigned int dim,
1093  const FEType & fe_t,
1094  const Elem * elem,
1095  const std::vector<Point> & p,
1096  std::vector<std::vector<Real>> & phi,
1097  const bool add_p_level)
1098 {
1099 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1100 
1101  if (elem && is_InfFE_elem(elem->type()))
1102  {
1103  for (auto i : index_range(phi))
1104  FEInterface::shapes<Real>(dim, fe_t, elem, i, p, phi[i], add_p_level);
1105  return;
1106  }
1107 #endif
1108 
1109  const Order o = fe_t.order;
1110 
1111  switch(dim)
1112  {
1113  case 0:
1114  fe_scalar_vec_error_switch(0, all_shapes(elem,o,p,phi,add_p_level), , ; return;);
1115  break;
1116  case 1:
1117  fe_scalar_vec_error_switch(1, all_shapes(elem,o,p,phi,add_p_level), , ; return;);
1118  break;
1119  case 2:
1120  fe_scalar_vec_error_switch(2, all_shapes(elem,o,p,phi,add_p_level), , ; return;);
1121  break;
1122  case 3:
1123  fe_scalar_vec_error_switch(3, all_shapes(elem,o,p,phi,add_p_level), , ; return;);
1124  break;
1125  default:
1126  libmesh_error_msg("Invalid dimension = " << dim);
1127  }
1128 
1129  return;
1130 }
1131 
1132 
1133 
1134 #ifdef LIBMESH_ENABLE_DEPRECATED
1135 template<>
1136 void FEInterface::shape<RealGradient>(const unsigned int dim,
1137  const FEType & fe_t,
1138  const ElemType t,
1139  const unsigned int i,
1140  const Point & p,
1141  RealGradient & phi)
1142 {
1143  libmesh_deprecated();
1144 
1145  // This API does not currently support infinite elements.
1146 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1147  if (is_InfFE_elem(t))
1148  {
1149  libmesh_not_implemented();
1150  }
1151 #endif
1152  const Order o = fe_t.order;
1153 
1154  switch(dim)
1155  {
1156  case 0:
1157  fe_vector_scalar_error_switch(0, shape(t,o,i,p), phi = , ; break;);
1158  break;
1159  case 1:
1160  fe_vector_scalar_error_switch(1, shape(t,o,i,p), phi = , ; break;);
1161  break;
1162  case 2:
1163  fe_vector_scalar_error_switch(2, shape(t,o,i,p), phi = , ; break;);
1164  break;
1165  case 3:
1166  fe_vector_scalar_error_switch(3, shape(t,o,i,p), phi = , ; break;);
1167  break;
1168  default:
1169  libmesh_error_msg("Invalid dimension = " << dim);
1170  }
1171 }
1172 #endif // LIBMESH_ENABLE_DEPRECATED
1173 
1174 
1175 
1176 template<>
1177 void FEInterface::shape<RealGradient>(const FEType & fe_t,
1178  const Elem * elem,
1179  const unsigned int i,
1180  const Point & p,
1181  RealGradient & phi)
1182 {
1183  // This API does not currently support infinite elements.
1184 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1185  if (is_InfFE_elem(elem->type()))
1186  {
1187  libmesh_not_implemented();
1188  }
1189 #endif
1190 
1191  auto dim = elem->dim();
1192 
1193  // We are calling
1194  //
1195  // FE<X,Y>::shape(Elem *, Order, unsigned, Point, true)
1196  //
1197  // with the last parameter set to "true" so that the Elem::p_level()
1198  // is accounted for internally. See fe.h for more details.
1199 
1200  switch(dim)
1201  {
1202  case 0:
1203  fe_vector_scalar_error_switch(0, shape(elem, fe_t.order, i, p, true), phi = , ; break;);
1204  break;
1205  case 1:
1206  fe_vector_scalar_error_switch(1, shape(elem, fe_t.order, i, p, true), phi = , ; break;);
1207  break;
1208  case 2:
1209  fe_vector_scalar_error_switch(2, shape(elem, fe_t.order, i, p, true), phi = , ; break;);
1210  break;
1211  case 3:
1212  fe_vector_scalar_error_switch(3, shape(elem, fe_t.order, i, p, true), phi = , ; break;);
1213  break;
1214  default:
1215  libmesh_error_msg("Invalid dimension = " << dim);
1216  }
1217 }
1218 
1219 
1220 
1221 template<>
1222 void FEInterface::shape<RealGradient>(const FEType & fe_t,
1223  int extra_order,
1224  const Elem * elem,
1225  const unsigned int i,
1226  const Point & p,
1227  RealGradient & phi)
1228 {
1229  // This API does not currently support infinite elements.
1230 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1231  if (is_InfFE_elem(elem->type()))
1232  {
1233  libmesh_not_implemented();
1234  }
1235 #endif
1236 
1237  auto dim = elem->dim();
1238 
1239  // We are calling
1240  //
1241  // FE<X,Y>::shape(Elem *, Order, unsigned, Point, false)
1242  //
1243  // with the last parameter set to "false" so that the
1244  // Elem::p_level() is not used internally and the "total_order" that
1245  // we compute is used instead. See fe.h for more details.
1246  auto total_order = fe_t.order + extra_order;
1247 
1248  switch(dim)
1249  {
1250  case 0:
1251  fe_vector_scalar_error_switch(0, shape(elem, total_order, i, p, false), phi = , ; break;);
1252  break;
1253  case 1:
1254  fe_vector_scalar_error_switch(1, shape(elem, total_order, i, p, false), phi = , ; break;);
1255  break;
1256  case 2:
1257  fe_vector_scalar_error_switch(2, shape(elem, total_order, i, p, false), phi = , ; break;);
1258  break;
1259  case 3:
1260  fe_vector_scalar_error_switch(3, shape(elem, total_order, i, p, false), phi = , ; break;);
1261  break;
1262  default:
1263  libmesh_error_msg("Invalid dimension = " << dim);
1264  }
1265 }
1266 
1267 
1268 
1269 template<>
1270 void FEInterface::shapes<RealGradient>(const unsigned int dim,
1271  const FEType & fe_t,
1272  const Elem * elem,
1273  const unsigned int i,
1274  const std::vector<Point> & p,
1275  std::vector<RealGradient> & phi,
1276  const bool add_p_level)
1277 {
1278  // This is actually an issue for infinite elements: They require type 'Gradient'!
1279  if (elem->infinite())
1280  libmesh_not_implemented();
1281 
1282  const Order o = fe_t.order;
1283 
1284  switch(dim)
1285  {
1286  case 0:
1287  fe_vector_scalar_error_switch(0, shapes(elem,o,i,p,phi,add_p_level), , ; return;);
1288  break;
1289  case 1:
1290  fe_vector_scalar_error_switch(1, shapes(elem,o,i,p,phi,add_p_level), , ; return;);
1291  break;
1292  case 2:
1293  fe_vector_scalar_error_switch(2, shapes(elem,o,i,p,phi,add_p_level), , ; return;);
1294  break;
1295  case 3:
1296  fe_vector_scalar_error_switch(3, shapes(elem,o,i,p,phi,add_p_level), , ; return;);
1297  break;
1298  default:
1299  libmesh_error_msg("Invalid dimension = " << dim);
1300  }
1301 
1302  return;
1303 }
1304 
1305 
1306 
1307 template<>
1308 void FEInterface::all_shapes<RealGradient>(const unsigned int dim,
1309  const FEType & fe_t,
1310  const Elem * elem,
1311  const std::vector<Point> & p,
1312  std::vector<std::vector<RealGradient>> & phi,
1313  const bool add_p_level)
1314 {
1315  // This is actually an issue for infinite elements: They require type 'Gradient'!
1316  if (elem->infinite())
1317  libmesh_not_implemented();
1318 
1319  const Order o = fe_t.order;
1320 
1321  switch(dim)
1322  {
1323  case 0:
1324  fe_vector_scalar_error_switch(0, all_shapes(elem,o,p,phi,add_p_level), , ; return;);
1325  break;
1326  case 1:
1327  fe_vector_scalar_error_switch(1, all_shapes(elem,o,p,phi,add_p_level), , ; return;);
1328  break;
1329  case 2:
1330  fe_vector_scalar_error_switch(2, all_shapes(elem,o,p,phi,add_p_level), , ; return;);
1331  break;
1332  case 3:
1333  fe_vector_scalar_error_switch(3, all_shapes(elem,o,p,phi,add_p_level), , ; return;);
1334  break;
1335  default:
1336  libmesh_error_msg("Invalid dimension = " << dim);
1337  }
1338 
1339  return;
1340 }
1341 
1342 
1343 #ifdef LIBMESH_ENABLE_DEPRECATED
1345 FEInterface::shape_function(const unsigned int dim,
1346  const FEType & fe_t,
1347  const ElemType libmesh_inf_var(t))
1348 {
1349  libmesh_deprecated();
1350 
1351 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1352  if (is_InfFE_elem(t))
1353  {
1354  inf_fe_switch(shape);
1355  }
1356 #endif
1357  fe_switch(shape);
1358 }
1359 #endif // LIBMESH_ENABLE_DEPRECATED
1360 
1361 
1362 
1365  const Elem * elem)
1366 {
1367  // dim is needed by the fe_switch macros below
1368  auto dim = elem->dim();
1369 
1370 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1371  if (is_InfFE_elem(elem->type()))
1372  inf_fe_switch(shape);
1373 #endif
1374  fe_switch(shape);
1375 }
1376 
1377 
1378 
1379 #ifdef LIBMESH_ENABLE_DEPRECATED
1380 Real FEInterface::shape_deriv(const unsigned int dim,
1381  const FEType & fe_t,
1382  const ElemType t,
1383  const unsigned int i,
1384  const unsigned int j,
1385  const Point & p)
1386 {
1387  libmesh_deprecated();
1388 
1389  libmesh_assert_greater (dim,j);
1390 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1391 
1392  if (is_InfFE_elem(t)){
1393  return ifem_shape_deriv(dim, fe_t, t, i, j, p);
1394  }
1395 
1396 #endif
1397 
1398  const Order o = fe_t.order;
1399  fe_switch(shape_deriv(t,o,i,j,p));
1400  return 0;
1401 }
1402 
1403 
1404 Real FEInterface::shape_deriv(const unsigned int dim,
1405  const FEType & fe_t,
1406  const Elem * elem,
1407  const unsigned int i,
1408  const unsigned int j,
1409  const Point & p)
1410 {
1411  libmesh_deprecated();
1412 
1413  libmesh_assert_greater (dim,j);
1414 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1415 
1416  if (elem->infinite()){
1417  return ifem_shape_deriv(fe_t, elem, i, j, p);
1418  }
1419 
1420 #endif
1421 
1422  const Order o = fe_t.order;
1423 
1424  switch(dim)
1425  {
1426  case 0:
1427  fe_family_switch (0, shape_deriv(elem, o, i, j, p), return , ;);
1428  case 1:
1429  fe_family_switch (1, shape_deriv(elem, o, i, j, p), return , ;);
1430  case 2:
1431  fe_family_switch (2, shape_deriv(elem, o, i, j, p), return , ;);
1432  case 3:
1433  fe_family_switch (3, shape_deriv(elem, o, i, j, p), return , ;);
1434  default:
1435  libmesh_error_msg("Invalid dimension = " << dim);
1436  }
1437  return 0;
1438 }
1439 #endif // LIBMESH_ENABLE_DEPRECATED
1440 
1441 
1442 
1444  const Elem * elem,
1445  const unsigned int i,
1446  const unsigned int j,
1447  const Point & p)
1448 {
1449  auto dim = elem->dim();
1450 
1451  libmesh_assert_greater (dim, j);
1452 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1453 
1454  if (is_InfFE_elem(elem->type())){
1455  return ifem_shape_deriv(fe_t, elem, i, j, p);
1456  }
1457 
1458 #endif
1459 
1460  // Account for Elem::p_level() when computing total order. Note: we are calling
1461  // FE::shape_deriv() with the final argument == true so that the Elem::p_level()
1462  // is accounted for automatically internally.
1463  fe_switch(shape_deriv(elem, fe_t.order, i, j, p, true));
1464 
1465  // We'll never get here
1466  return 0;
1467 }
1468 
1469 
1470 
1472  int extra_order,
1473  const Elem * elem,
1474  const unsigned int i,
1475  const unsigned int j,
1476  const Point & p)
1477 {
1478  auto dim = elem->dim();
1479 
1480  libmesh_assert_greater (dim, j);
1481 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1482 
1483  if (elem->infinite()){
1484  return ifem_shape_deriv(fe_t, elem, i, j, p);
1485  }
1486 
1487 #endif
1488 
1489  // Ignore Elem::p_level() when computing total order, use
1490  // extra_order instead.
1491  auto total_order = fe_t.order + extra_order;
1492 
1493  // We call shape_deriv() with the final argument == false so that
1494  // the Elem::p_level() is ignored internally.
1495  switch(dim)
1496  {
1497  case 0:
1498  fe_family_switch (0, shape_deriv(elem, total_order, i, j, p, false), return , ;);
1499  case 1:
1500  fe_family_switch (1, shape_deriv(elem, total_order, i, j, p, false), return , ;);
1501  case 2:
1502  fe_family_switch (2, shape_deriv(elem, total_order, i, j, p, false), return , ;);
1503  case 3:
1504  fe_family_switch (3, shape_deriv(elem, total_order, i, j, p, false), return , ;);
1505  default:
1506  libmesh_error_msg("Invalid dimension = " << dim);
1507  }
1508 
1509  // We'll never get here
1510  return 0;
1511 }
1512 
1513 
1514 
1515 template<>
1516 void FEInterface::shape_derivs<Real>(const FEType & fe_t,
1517  const Elem * elem,
1518  const unsigned int i,
1519  const unsigned int j,
1520  const std::vector<Point> & p,
1521  std::vector<Real> & dphi,
1522  const bool add_p_level)
1523 {
1524  const Order o = fe_t.order;
1525 
1526  switch(elem->dim())
1527  {
1528  case 0:
1529  fe_scalar_vec_error_switch(0, shape_derivs(elem,o,i,j,p,dphi,add_p_level), , ; break;);
1530  break;
1531  case 1:
1532  fe_scalar_vec_error_switch(1, shape_derivs(elem,o,i,j,p,dphi,add_p_level), , ; break;);
1533  break;
1534  case 2:
1535  fe_scalar_vec_error_switch(2, shape_derivs(elem,o,i,j,p,dphi,add_p_level), , ; break;);
1536  break;
1537  case 3:
1538  fe_scalar_vec_error_switch(3, shape_derivs(elem,o,i,j,p,dphi,add_p_level), , ; break;);
1539  break;
1540  default:
1541  libmesh_error_msg("Invalid dimension = " << elem->dim());
1542  }
1543 
1544  return;
1545 }
1546 
1547 
1548 
1549 // Only instantiating Real because we don't expect vector-valued
1550 // mapping functions, and that's what's using this interface.
1551 template<>
1552 void FEInterface::all_shape_derivs<Real>(const unsigned int dim,
1553  const FEType & fe_t,
1554  const Elem * elem,
1555  const std::vector<Point> & p,
1556  std::vector<std::vector<Real>> * comps[3],
1557  const bool add_p_level)
1558 {
1559 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1560 
1561  if (elem && is_InfFE_elem(elem->type()))
1562  {
1563  for (auto j : make_range(dim))
1564  for (auto i : index_range(*comps[j]))
1565  FEInterface::shape_derivs<Real>(fe_t, elem, i, j, p, (*comps[j])[i], add_p_level);
1566  return;
1567  }
1568 #endif
1569 
1570  const Order o = fe_t.order;
1571 
1572  switch(dim)
1573  {
1574  case 0:
1575  fe_scalar_vec_error_switch(0, all_shape_derivs(elem,o,p,comps,add_p_level), , ; return;);
1576  break;
1577  case 1:
1578  fe_scalar_vec_error_switch(1, all_shape_derivs(elem,o,p,comps,add_p_level), , ; return;);
1579  break;
1580  case 2:
1581  fe_scalar_vec_error_switch(2, all_shape_derivs(elem,o,p,comps,add_p_level), , ; return;);
1582  break;
1583  case 3:
1584  fe_scalar_vec_error_switch(3, all_shape_derivs(elem,o,p,comps,add_p_level), , ; return;);
1585  break;
1586  default:
1587  libmesh_error_msg("Invalid dimension = " << dim);
1588  }
1589 
1590  return;
1591 }
1592 
1593 
1594 
1595 template<>
1596 void FEInterface::shape_derivs<RealGradient>(const FEType & fe_t,
1597  const Elem * elem,
1598  const unsigned int i,
1599  const unsigned int j,
1600  const std::vector<Point> & p,
1601  std::vector<RealGradient> & dphi,
1602  const bool add_p_level)
1603 {
1604  const Order o = fe_t.order;
1605 
1606  switch(elem->dim())
1607  {
1608  case 0:
1609  fe_vector_scalar_error_switch(0, shape_derivs(elem,o,i,j,p,dphi,add_p_level), , ; break;);
1610  break;
1611  case 1:
1612  fe_vector_scalar_error_switch(1, shape_derivs(elem,o,i,j,p,dphi,add_p_level), , ; break;);
1613  break;
1614  case 2:
1615  fe_vector_scalar_error_switch(2, shape_derivs(elem,o,i,j,p,dphi,add_p_level), , ; break;);
1616  break;
1617  case 3:
1618  fe_vector_scalar_error_switch(3, shape_derivs(elem,o,i,j,p,dphi,add_p_level), , ; break;);
1619  break;
1620  default:
1621  libmesh_error_msg("Invalid dimension = " << elem->dim());
1622  }
1623 
1624  return;
1625 }
1626 
1627 
1628 
1629 #ifdef LIBMESH_ENABLE_DEPRECATED
1631 FEInterface::shape_deriv_function(const unsigned int dim,
1632  const FEType & fe_t,
1633  const ElemType libmesh_inf_var(t))
1634 {
1635  libmesh_deprecated();
1636 
1637 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1638  if (is_InfFE_elem(t))
1639  {
1640  inf_fe_switch(shape_deriv);
1641  }
1642 #endif
1643  fe_switch(shape_deriv);
1644 }
1645 #endif // LIBMESH_ENABLE_DEPRECATED
1646 
1647 
1648 
1651  const Elem * elem)
1652 {
1653  // dim is needed by the fe_switch macros below
1654  auto dim = elem->dim();
1655 
1656 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1657  if (is_InfFE_elem(elem->type()))
1658  {
1659  inf_fe_switch(shape_deriv);
1660  }
1661 #endif
1662  fe_switch(shape_deriv);
1663 }
1664 
1665 
1666 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1667 #ifdef LIBMESH_ENABLE_DEPRECATED
1669  const FEType & fe_t,
1670  const ElemType t,
1671  const unsigned int i,
1672  const unsigned int j,
1673  const Point & p)
1674 {
1675  libmesh_deprecated();
1676 
1677  libmesh_assert_greater_equal (dim*(dim-1),j);
1678 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1679  if (is_InfFE_elem(t))
1680  libmesh_not_implemented();
1681 #endif
1682 
1683  const Order o = fe_t.order;
1684 
1685  switch(dim)
1686  {
1687  case 0:
1688  fe_family_switch (0, shape_second_deriv(t, o, i, j, p), return , ;);
1689  case 1:
1690  fe_family_switch (1, shape_second_deriv(t, o, i, j, p), return , ;);
1691  case 2:
1692  fe_family_switch (2, shape_second_deriv(t, o, i, j, p), return , ;);
1693  case 3:
1694  fe_family_switch (3, shape_second_deriv(t, o, i, j, p), return , ;);
1695  default:
1696  libmesh_error_msg("Invalid dimension = " << dim);
1697  }
1698  return 0;
1699 }
1700 
1701 
1703  const FEType & fe_t,
1704  const Elem * elem,
1705  const unsigned int i,
1706  const unsigned int j,
1707  const Point & p)
1708 {
1709  libmesh_deprecated();
1710 
1711  libmesh_assert_greater_equal (dim*(dim-1),j);
1712 
1713  if (elem->infinite())
1714  libmesh_not_implemented();
1715 
1716  const Order o = fe_t.order;
1717 
1718  switch(dim)
1719  {
1720  case 0:
1721  fe_family_switch (0, shape_second_deriv(elem, o, i, j, p), return , ;);
1722  case 1:
1723  fe_family_switch (1, shape_second_deriv(elem, o, i, j, p), return , ;);
1724  case 2:
1725  fe_family_switch (2, shape_second_deriv(elem, o, i, j, p), return , ;);
1726  case 3:
1727  fe_family_switch (3, shape_second_deriv(elem, o, i, j, p), return , ;);
1728  default:
1729  libmesh_error_msg("Invalid dimension = " << dim);
1730  }
1731  return 0;
1732 }
1733 #endif // LIBMESH_ENABLE_DEPRECATED
1734 
1735 
1736 
1738  const Elem * elem,
1739  const unsigned int i,
1740  const unsigned int j,
1741  const Point & p)
1742 {
1743  auto dim = elem->dim();
1744 
1745  libmesh_assert_greater_equal (dim*(dim-1),j);
1746 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1747  if (is_InfFE_elem(elem->type()))
1748  libmesh_not_implemented();
1749 #endif
1750 
1751  // We are calling FE::shape_second_deriv() with the final argument
1752  // == true so that the Elem::p_level() is accounted for
1753  // automatically internally.
1754  switch(dim)
1755  {
1756  case 0:
1757  fe_family_switch (0, shape_second_deriv(elem, fe_t.order, i, j, p, true), return , ;);
1758  case 1:
1759  fe_family_switch (1, shape_second_deriv(elem, fe_t.order, i, j, p, true), return , ;);
1760  case 2:
1761  fe_family_switch (2, shape_second_deriv(elem, fe_t.order, i, j, p, true), return , ;);
1762  case 3:
1763  fe_family_switch (3, shape_second_deriv(elem, fe_t.order, i, j, p, true), return , ;);
1764  default:
1765  libmesh_error_msg("Invalid dimension = " << dim);
1766  }
1767 
1768  // We'll never get here
1769  return 0;
1770 }
1771 
1772 
1773 
1775  int extra_order,
1776  const Elem * elem,
1777  const unsigned int i,
1778  const unsigned int j,
1779  const Point & p)
1780 {
1781  auto dim = elem->dim();
1782 
1783  libmesh_assert_greater_equal (dim*(dim-1),j);
1784 
1785  if (elem->infinite())
1786  libmesh_not_implemented();
1787 
1788  // Ignore Elem::p_level() when computing total order, use
1789  // extra_order instead.
1790  auto total_order = fe_t.order + extra_order;
1791 
1792  // We are calling FE::shape_second_deriv() with the final argument
1793  // == false so that the Elem::p_level() is ignored and the
1794  // total_order we compute is used instead.
1795  switch(dim)
1796  {
1797  case 0:
1798  fe_family_switch (0, shape_second_deriv(elem, total_order, i, j, p, false), return , ;);
1799  case 1:
1800  fe_family_switch (1, shape_second_deriv(elem, total_order, i, j, p, false), return , ;);
1801  case 2:
1802  fe_family_switch (2, shape_second_deriv(elem, total_order, i, j, p, false), return , ;);
1803  case 3:
1804  fe_family_switch (3, shape_second_deriv(elem, total_order, i, j, p, false), return , ;);
1805  default:
1806  libmesh_error_msg("Invalid dimension = " << dim);
1807  }
1808 
1809  // We'll never get here
1810  return 0;
1811 }
1812 
1813 
1814 
1815 #ifdef LIBMESH_ENABLE_DEPRECATED
1818  const FEType & fe_t,
1819  const ElemType libmesh_inf_var(t))
1820 {
1821  libmesh_deprecated();
1822 
1823 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1824  if (is_InfFE_elem(t))
1825  libmesh_not_implemented();
1826 #endif
1827  fe_switch(shape_second_deriv);
1828 }
1829 #endif // LIBMESH_ENABLE_DEPRECATED
1830 
1831 
1832 
1835  const Elem * elem)
1836 {
1837  auto dim = elem->dim();
1838 
1839 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1840  if (is_InfFE_elem(elem->type()))
1841  libmesh_not_implemented();
1842 #endif
1843  fe_switch(shape_second_deriv);
1844 }
1845 
1846 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
1847 
1848 
1849 #ifdef LIBMESH_ENABLE_DEPRECATED
1850 template<>
1851 void FEInterface::shape<RealGradient>(const unsigned int dim,
1852  const FEType & fe_t,
1853  const Elem * elem,
1854  const unsigned int i,
1855  const Point & p,
1856  RealGradient & phi)
1857 {
1858  libmesh_deprecated();
1859 
1860  // This is actually an issue for infinite elements: They require type 'Gradient'!
1861  if (elem->infinite())
1862  libmesh_not_implemented();
1863 
1864  const Order o = fe_t.order;
1865 
1866  switch(dim)
1867  {
1868  case 0:
1869  fe_vector_scalar_error_switch(0, shape(elem,o,i,p), phi = , ; break;);
1870  break;
1871  case 1:
1872  fe_vector_scalar_error_switch(1, shape(elem,o,i,p), phi = , ; break;);
1873  break;
1874  case 2:
1875  fe_vector_scalar_error_switch(2, shape(elem,o,i,p), phi = , ; break;);
1876  break;
1877  case 3:
1878  fe_vector_scalar_error_switch(3, shape(elem,o,i,p), phi = , ; break;);
1879  break;
1880  default:
1881  libmesh_error_msg("Invalid dimension = " << dim);
1882  }
1883 
1884  return;
1885 }
1886 #endif // LIBMESH_ENABLE_DEPRECATED
1887 
1888 
1889 void FEInterface::compute_data(const unsigned int dim,
1890  const FEType & fe_t,
1891  const Elem * elem,
1892  FEComputeData & data)
1893 {
1894 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1895 
1896  if (elem && is_InfFE_elem(elem->type()))
1897  {
1898  data.init();
1899  ifem_compute_data(dim, fe_t, elem, data);
1900  return;
1901  }
1902 
1903 #endif
1904 
1905  const unsigned int n_dof = n_dofs (fe_t, elem);
1906  const Point & p = data.p;
1907  data.shape.resize(n_dof);
1908 
1909  if (data.need_derivative())
1910  {
1911  data.dshape.resize(n_dof);
1912  data.local_transform.resize(dim);
1913 
1914  for (unsigned int d=0; d<dim; d++)
1915  data.local_transform[d].resize(dim);
1916 
1917  auto fe = FEBase::build(dim, fe_t);
1918  std::vector<Point> pt = {p};
1919  fe->get_dphideta(); // to compute the map
1920  fe->reinit(elem, &pt);
1921 
1922  // compute the reference->physical map.
1923  data.local_transform[0][0] = fe->get_dxidx()[0];
1924  if (dim > 1)
1925  {
1926  data.local_transform[1][0] = fe->get_detadx()[0];
1927  data.local_transform[1][1] = fe->get_detady()[0];
1928  data.local_transform[0][1] = fe->get_dxidy()[0];
1929  if (dim > 2)
1930  {
1931  data.local_transform[2][0] = fe->get_dzetadx()[0];
1932  data.local_transform[2][1] = fe->get_dzetady()[0];
1933  data.local_transform[2][2] = fe->get_dzetadz()[0];
1934  data.local_transform[1][2] = fe->get_detadz()[0];
1935  data.local_transform[0][2] = fe->get_dxidz()[0];
1936  }
1937  }
1938  }
1939 
1940  // set default values for all the output fields
1941  data.init();
1942 
1943  for (unsigned int n=0; n<n_dof; n++)
1944  {
1945  // Here we pass the original fe_t object. Additional p-levels
1946  // (if any) are handled internally by the shape() and
1947  // shape_deriv() functions since they have access to the elem
1948  // pointer. Note that we are already using the n_dof value
1949  // appropriate to the elevated p-level.
1950  data.shape[n] = shape(fe_t, elem, n, p);
1951  if (data.need_derivative())
1952  {
1953  for (unsigned int j=0; j<dim; j++)
1954  data.dshape[n](j) = shape_deriv(fe_t, elem, n, j, p);
1955  }
1956  }
1957 }
1958 
1959 
1960 
1961 #ifdef LIBMESH_ENABLE_AMR
1962 
1964  DofMap & dof_map,
1965  const unsigned int variable_number,
1966  const Elem * elem)
1967 {
1968  libmesh_assert(elem);
1969 
1970  const FEType & fe_t = dof_map.variable_type(variable_number);
1971 
1972  switch (elem->dim())
1973  {
1974  case 0:
1975  case 1:
1976  {
1977  // No constraints in 0D/1D.
1978  return;
1979  }
1980 
1981 
1982  case 2:
1983  {
1984  switch (fe_t.family)
1985  {
1986  case CLOUGH:
1987  FE<2,CLOUGH>::compute_constraints (constraints,
1988  dof_map,
1989  variable_number,
1990  elem); return;
1991 
1992  case HERMITE:
1994  dof_map,
1995  variable_number,
1996  elem); return;
1997 
1998  case LAGRANGE:
2000  dof_map,
2001  variable_number,
2002  elem); return;
2003 
2004  case HIERARCHIC:
2006  dof_map,
2007  variable_number,
2008  elem); return;
2009 
2010  case HIERARCHIC_VEC:
2012  dof_map,
2013  variable_number,
2014  elem); return;
2015 
2016  case SIDE_HIERARCHIC:
2018  dof_map,
2019  variable_number,
2020  elem); return;
2021 
2022  case LAGRANGE_VEC:
2024  dof_map,
2025  variable_number,
2026  elem); return;
2027 
2028  fe_family_horder_case_func(2, compute_constraints (constraints,
2029  dof_map,
2030  variable_number,
2031  elem), , ; return;)
2032  default:
2033  return;
2034  }
2035  }
2036 
2037 
2038  case 3:
2039  {
2040  switch (fe_t.family)
2041  {
2042  case HERMITE:
2044  dof_map,
2045  variable_number,
2046  elem); return;
2047 
2048  case LAGRANGE:
2050  dof_map,
2051  variable_number,
2052  elem); return;
2053 
2054  case HIERARCHIC:
2056  dof_map,
2057  variable_number,
2058  elem); return;
2059 
2060  case SIDE_HIERARCHIC:
2062  dof_map,
2063  variable_number,
2064  elem); return;
2065 
2066  case LAGRANGE_VEC:
2068  dof_map,
2069  variable_number,
2070  elem); return;
2071 
2072  case HIERARCHIC_VEC:
2074  dof_map,
2075  variable_number,
2076  elem); return;
2077 
2078  fe_family_horder_case_func(3, compute_constraints (constraints,
2079  dof_map,
2080  variable_number,
2081  elem), , ; return;)
2082  default:
2083  return;
2084  }
2085  }
2086 
2087 
2088  default:
2089  libmesh_error_msg("Invalid dimension = " << elem->dim());
2090  }
2091 }
2092 
2093 #endif // #ifdef LIBMESH_ENABLE_AMR
2094 
2095 
2096 
2097 #ifdef LIBMESH_ENABLE_PERIODIC
2098 
2100  DofMap & dof_map,
2101  const PeriodicBoundaries & boundaries,
2102  const MeshBase & mesh,
2103  const PointLocatorBase * point_locator,
2104  const unsigned int variable_number,
2105  const Elem * elem)
2106 {
2107  const FEType & fe_t = dof_map.variable_type(variable_number);
2108  // No element-specific optimizations currently exist, although
2109  // we do have to select the right compute_periodic_constraints
2110  // for the OutputType of this FEType.
2111  switch (field_type(fe_t))
2112  {
2113  case TYPE_SCALAR:
2115  constraints, dof_map, boundaries, mesh, point_locator, variable_number, elem);
2116  break;
2117  case TYPE_VECTOR:
2119  constraints, dof_map, boundaries, mesh, point_locator, variable_number, elem);
2120  break;
2121  default:
2122  libmesh_error_msg(
2123  "compute_periodic_constraints only set up for vector or scalar FEFieldTypes");
2124  }
2125 }
2126 
2127 #endif // #ifdef LIBMESH_ENABLE_PERIODIC
2128 
2129 
2130 
2131 unsigned int FEInterface::max_order(const FEType & fe_t,
2132  const ElemType & el_t)
2133 {
2134  // Yeah, I know, infinity is much larger than 11, but our
2135  // solvers don't seem to like high degree polynomials, and our
2136  // quadrature rules and number_lookups tables
2137  // need to go up higher.
2138  const unsigned int unlimited = 11;
2139 
2140  // If we used 0 as a default, then elements missing from this
2141  // table (e.g. infinite elements) would be considered broken.
2142  const unsigned int unknown = unlimited;
2143 
2144  switch (fe_t.family)
2145  {
2146  case LAGRANGE:
2147  case LAGRANGE_VEC:
2148  switch (el_t)
2149  {
2150  case EDGE2:
2151  case EDGE3:
2152  case EDGE4:
2153  return 3;
2154  case TRI3:
2155  case TRISHELL3:
2156  case C0POLYGON:
2157  return 1;
2158  case TRI6:
2159  return 2;
2160  case TRI7:
2161  return 3;
2162  case QUAD4:
2163  case QUADSHELL4:
2164  return 1;
2165  case QUAD8:
2166  case QUADSHELL8:
2167  case QUAD9:
2168  case QUADSHELL9:
2169  return 2;
2170  case TET4:
2171  return 1;
2172  case TET10:
2173  return 2;
2174  case TET14:
2175  return 3;
2176  case HEX8:
2177  return 1;
2178  case HEX20:
2179  case HEX27:
2180  return 2;
2181  case PRISM6:
2182  return 1;
2183  case PRISM15:
2184  case PRISM18:
2185  return 2;
2186  case PRISM20:
2187  case PRISM21:
2188  return 3;
2189  case PYRAMID5:
2190  return 1;
2191  case PYRAMID13:
2192  case PYRAMID14:
2193  return 2;
2194  case PYRAMID18:
2195  return 3;
2196  case C0POLYHEDRON:
2197  return 1;
2198  default:
2199  return unknown;
2200  }
2201  break;
2202  case MONOMIAL:
2203  case L2_LAGRANGE:
2204  case L2_LAGRANGE_VEC:
2205  case L2_HIERARCHIC:
2206  case L2_HIERARCHIC_VEC:
2207  case MONOMIAL_VEC:
2208  switch (el_t)
2209  {
2210  case EDGE2:
2211  case EDGE3:
2212  case EDGE4:
2213  case C0POLYGON:
2214  case TRI3:
2215  case TRISHELL3:
2216  case TRI6:
2217  case TRI7:
2218  case QUAD4:
2219  case QUADSHELL4:
2220  case QUAD8:
2221  case QUADSHELL8:
2222  case QUAD9:
2223  case QUADSHELL9:
2224  case TET4:
2225  case TET10:
2226  case TET14:
2227  case HEX8:
2228  case HEX20:
2229  case HEX27:
2230  case PRISM6:
2231  case PRISM15:
2232  case PRISM18:
2233  case PRISM20:
2234  case PRISM21:
2235  case PYRAMID5:
2236  case PYRAMID13:
2237  case PYRAMID14:
2238  case PYRAMID18:
2239  case C0POLYHEDRON:
2240  return unlimited;
2241  default:
2242  return unknown;
2243  }
2244  break;
2245 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
2246  case BERNSTEIN:
2247  case RATIONAL_BERNSTEIN:
2248  switch (el_t)
2249  {
2250  case EDGE2:
2251  case EDGE3:
2252  case EDGE4:
2253  return unlimited;
2254  case TRI3:
2255  case TRISHELL3:
2256  return 1;
2257  case TRI6:
2258  case TRI7:
2259  return 6;
2260  case QUAD4:
2261  case QUADSHELL4:
2262  return 1;
2263  case QUAD8:
2264  case QUADSHELL8:
2265  case QUAD9:
2266  case QUADSHELL9:
2267  return unlimited;
2268  case TET4:
2269  return 1;
2270  case TET10:
2271  case TET14:
2272  return 2;
2273  case HEX8:
2274  return 1;
2275  case HEX20:
2276  return 2;
2277  case HEX27:
2278  return 4;
2279  case C0POLYGON:
2280  case PRISM6:
2281  case PRISM15:
2282  case PRISM18:
2283  case PRISM20:
2284  case PRISM21:
2285  case PYRAMID5:
2286  case PYRAMID13:
2287  case PYRAMID14:
2288  case PYRAMID18:
2289  case C0POLYHEDRON:
2290  return 0;
2291  default:
2292  return unknown;
2293  }
2294  break;
2295  case SZABAB:
2296  switch (el_t)
2297  {
2298  case EDGE2:
2299  return 1;
2300  case EDGE3:
2301  case EDGE4:
2302  return 7;
2303  case TRI3:
2304  case TRISHELL3:
2305  return 1;
2306  case TRI6:
2307  case TRI7:
2308  return 7;
2309  case QUAD4:
2310  case QUADSHELL4:
2311  return 1;
2312  case QUAD8:
2313  case QUADSHELL8:
2314  case QUAD9:
2315  case QUADSHELL9:
2316  return 7;
2317  case C0POLYGON:
2318  case TET4:
2319  case TET10:
2320  case TET14:
2321  case HEX8:
2322  case HEX20:
2323  case HEX27:
2324  case PRISM6:
2325  case PRISM15:
2326  case PRISM18:
2327  case PRISM20:
2328  case PRISM21:
2329  case PYRAMID5:
2330  case PYRAMID13:
2331  case PYRAMID14:
2332  case PYRAMID18:
2333  case C0POLYHEDRON:
2334  return 0;
2335  default:
2336  return unknown;
2337  }
2338  break;
2339 #endif
2340  case XYZ:
2341  switch (el_t)
2342  {
2343  case EDGE2:
2344  case EDGE3:
2345  case EDGE4:
2346  case C0POLYGON:
2347  case TRI3:
2348  case TRISHELL3:
2349  case TRI6:
2350  case TRI7:
2351  case QUAD4:
2352  case QUADSHELL4:
2353  case QUAD8:
2354  case QUADSHELL8:
2355  case QUAD9:
2356  case QUADSHELL9:
2357  case TET4:
2358  case TET10:
2359  case TET14:
2360  case HEX8:
2361  case HEX20:
2362  case HEX27:
2363  case PRISM6:
2364  case PRISM15:
2365  case PRISM18:
2366  case PRISM20:
2367  case PRISM21:
2368  case PYRAMID5:
2369  case PYRAMID13:
2370  case PYRAMID14:
2371  case PYRAMID18:
2372  case C0POLYHEDRON:
2373  return unlimited;
2374  default:
2375  return unknown;
2376  }
2377  break;
2378  case CLOUGH:
2379  switch (el_t)
2380  {
2381  case EDGE2:
2382  case EDGE3:
2383  return 3;
2384  case EDGE4:
2385  case TRI3:
2386  case TRISHELL3:
2387  return 0;
2388  case TRI6:
2389  case TRI7:
2390  return 3;
2391  case C0POLYGON:
2392  case QUAD4:
2393  case QUADSHELL4:
2394  case QUAD8:
2395  case QUADSHELL8:
2396  case QUAD9:
2397  case QUADSHELL9:
2398  case TET4:
2399  case TET10:
2400  case TET14:
2401  case HEX8:
2402  case HEX20:
2403  case HEX27:
2404  case PRISM6:
2405  case PRISM15:
2406  case PRISM18:
2407  case PRISM20:
2408  case PRISM21:
2409  case PYRAMID5:
2410  case PYRAMID13:
2411  case PYRAMID14:
2412  case PYRAMID18:
2413  case C0POLYHEDRON:
2414  return 0;
2415  default:
2416  return unknown;
2417  }
2418  break;
2419  case HERMITE:
2420  switch (el_t)
2421  {
2422  case EDGE2:
2423  case EDGE3:
2424  return unlimited;
2425  case EDGE4:
2426  case TRI3:
2427  case TRISHELL3:
2428  case TRI6:
2429  case TRI7:
2430  return 0;
2431  case QUAD4:
2432  case QUADSHELL4:
2433  return 3;
2434  case QUAD8:
2435  case QUADSHELL8:
2436  case QUAD9:
2437  case QUADSHELL9:
2438  return unlimited;
2439  case TET4:
2440  case TET10:
2441  case TET14:
2442  return 0;
2443  case HEX8:
2444  return 3;
2445  case HEX20:
2446  case HEX27:
2447  return unlimited;
2448  case C0POLYGON:
2449  case PRISM6:
2450  case PRISM15:
2451  case PRISM18:
2452  case PRISM20:
2453  case PRISM21:
2454  case PYRAMID5:
2455  case PYRAMID13:
2456  case PYRAMID14:
2457  case PYRAMID18:
2458  case C0POLYHEDRON:
2459  return 0;
2460  default:
2461  return unknown;
2462  }
2463  break;
2464  case HIERARCHIC:
2465  case HIERARCHIC_VEC:
2466  switch (el_t)
2467  {
2468  case EDGE2:
2469  case EDGE3:
2470  case EDGE4:
2471  return unlimited;
2472  case TRI3:
2473  case TRISHELL3:
2474  return 1;
2475  case TRI6:
2476  case TRI7:
2477  return unlimited;
2478  case QUAD4:
2479  case QUADSHELL4:
2480  return 1;
2481  case QUAD8:
2482  case QUADSHELL8:
2483  case QUAD9:
2484  case QUADSHELL9:
2485  return unlimited;
2486  case TET4:
2487  return 1;
2488  case TET10:
2489  return 2;
2490  case TET14:
2491  return unlimited;
2492  case HEX8:
2493  case HEX20:
2494  return 1;
2495  case HEX27:
2496  return unlimited;
2497  case PRISM6:
2498  case PRISM15:
2499  case PRISM18:
2500  case PRISM20:
2501  case PRISM21:
2502  return unlimited;
2503  case C0POLYGON:
2504  case PYRAMID5:
2505  case PYRAMID13:
2506  case PYRAMID14:
2507  case PYRAMID18:
2508  case C0POLYHEDRON:
2509  return 0;
2510  default:
2511  return unknown;
2512  }
2513  break;
2514  case SIDE_HIERARCHIC:
2515  switch (el_t)
2516  {
2517  case EDGE2:
2518  case EDGE3:
2519  case EDGE4:
2520  return unlimited; // although it's all the same as 0...
2521  case TRI3:
2522  case TRISHELL3:
2523  return 0;
2524  case TRI6:
2525  case TRI7:
2526  return unlimited;
2527  case C0POLYGON:
2528  case QUAD4:
2529  case QUADSHELL4:
2530  return 0;
2531  case QUAD8:
2532  case QUADSHELL8:
2533  case QUAD9:
2534  case QUADSHELL9:
2535  return unlimited;
2536  case TET4:
2537  case TET10:
2538  return 0;
2539  case TET14:
2540  return unlimited;
2541  case HEX8:
2542  case HEX20:
2543  return 0;
2544  case HEX27:
2545  return unlimited;
2546  case PRISM6:
2547  case PRISM15:
2548  case PRISM18:
2549  return 0;
2550  case PRISM20:
2551  case PRISM21:
2552  return unlimited;
2553  case PYRAMID5:
2554  case PYRAMID13:
2555  case PYRAMID14:
2556  case PYRAMID18:
2557  case C0POLYHEDRON:
2558  return 0;
2559  default:
2560  return unknown;
2561  }
2562  break;
2563  case SUBDIVISION:
2564  switch (el_t)
2565  {
2566  case TRI3SUBDIVISION:
2567  return unlimited;
2568  default:
2569  return unknown;
2570  }
2571  break;
2572  case NEDELEC_ONE:
2573  switch (el_t)
2574  {
2575  case TRI6:
2576  case TRI7:
2577  case QUAD8:
2578  case QUAD9:
2579  return 5;
2580  case TET10:
2581  case TET14:
2582  case HEX20:
2583  case HEX27:
2584  return 1;
2585  default:
2586  return 0;
2587  }
2588  break;
2589  case RAVIART_THOMAS:
2590  case L2_RAVIART_THOMAS:
2591  switch (el_t)
2592  {
2593  case TRI6:
2594  case TRI7:
2595  case QUAD8:
2596  case QUAD9:
2597  return 5;
2598  case TET14:
2599  case HEX27:
2600  return 1;
2601  default:
2602  return 0;
2603  }
2604  break;
2605  default:
2606  return 0;
2607  break;
2608  }
2609 }
2610 
2611 
2612 
2614 {
2615  switch (fe_t.family)
2616  {
2617  case LAGRANGE:
2618  case L2_LAGRANGE:
2619  case MONOMIAL:
2620  case MONOMIAL_VEC:
2621  case L2_HIERARCHIC:
2622  case SIDE_HIERARCHIC:
2623  case XYZ:
2624  case SUBDIVISION:
2625  case LAGRANGE_VEC:
2626  case L2_LAGRANGE_VEC:
2627  case NEDELEC_ONE:
2628  case RAVIART_THOMAS:
2629  case L2_RAVIART_THOMAS:
2630  return false;
2631  case CLOUGH:
2632  case HERMITE:
2633  case HIERARCHIC:
2634  case HIERARCHIC_VEC:
2635  case L2_HIERARCHIC_VEC:
2636  fe_family_horder_case()
2637  default:
2638  return true;
2639  }
2640 }
2641 
2643 {
2644  return FEInterface::field_type(fe_type.family);
2645 }
2646 
2648 {
2649  switch (fe_family)
2650  {
2651  fe_family_vector_case()
2652  return TYPE_VECTOR;
2653  default:
2654  return TYPE_SCALAR;
2655  }
2656 }
2657 
2659 {
2660  switch (fe_family)
2661  {
2662  case HIERARCHIC:
2663  case L2_HIERARCHIC:
2664  case HIERARCHIC_VEC:
2665  case L2_HIERARCHIC_VEC:
2666  case BERNSTEIN:
2667  case RATIONAL_BERNSTEIN:
2668  case SZABAB:
2669  case NEDELEC_ONE:
2670  case RAVIART_THOMAS:
2671  return true;
2672  default:
2673  return false;
2674  }
2675 }
2676 
2677 unsigned int FEInterface::n_vec_dim (const MeshBase & mesh,
2678  const FEType & fe_type)
2679 {
2680  // We assume the number of vector components is the mesh spatial dimension.
2681  return field_type(fe_type.family) == TYPE_VECTOR ? mesh.spatial_dimension() : 1;
2682 }
2683 
2684 
2685 
2686 
2687 bool
2689 {
2690  switch (fe_type.family)
2691  {
2692  case HERMITE:
2693  case HIERARCHIC:
2694  case HIERARCHIC_VEC:
2695  case L2_HIERARCHIC:
2696  case L2_HIERARCHIC_VEC:
2697  case MONOMIAL:
2698  case MONOMIAL_VEC:
2699  case SIDE_HIERARCHIC:
2700  case SZABAB:
2701  case XYZ:
2702  return true;
2703 
2704  case BERNSTEIN:
2705  case CLOUGH: // maybe some day?
2706  case LAGRANGE:
2707  case LAGRANGE_VEC:
2708  case L2_LAGRANGE:
2709  case L2_LAGRANGE_VEC:
2710  case NEDELEC_ONE:
2711  case RATIONAL_BERNSTEIN:
2712  case RAVIART_THOMAS:
2713  case L2_RAVIART_THOMAS:
2714  case SCALAR:
2715  case SUBDIVISION:
2716  return false;
2717 
2718  default:
2719  libmesh_error_msg("Unknown FE Family " << Utility::enum_to_string(fe_type.family));
2720  }
2721 }
2722 
2723 
2724 
2726 {
2727  switch (fe_type.family)
2728  {
2729  // Discontinuous elements
2730  case MONOMIAL:
2731  case MONOMIAL_VEC:
2732  case L2_HIERARCHIC:
2733  case L2_LAGRANGE:
2734  case XYZ:
2735  case SCALAR:
2736  case L2_RAVIART_THOMAS:
2737  case L2_HIERARCHIC_VEC:
2738  case L2_LAGRANGE_VEC:
2739  return DISCONTINUOUS;
2740 
2741  // C0 elements
2742  case LAGRANGE:
2743  case HIERARCHIC:
2744  case BERNSTEIN:
2745  case SZABAB:
2746  case RATIONAL_BERNSTEIN:
2747  case INFINITE_MAP:
2748  case JACOBI_20_00:
2749  case JACOBI_30_00:
2750  case LEGENDRE:
2751  case LAGRANGE_VEC:
2752  case HIERARCHIC_VEC:
2753  return C_ZERO;
2754 
2755  // C1 elements
2756  case CLOUGH:
2757  case HERMITE:
2758  case SUBDIVISION:
2759  return C_ONE;
2760 
2761  case NEDELEC_ONE:
2762  return H_CURL;
2763 
2764  case RAVIART_THOMAS:
2765  return H_DIV;
2766 
2767  // Side elements
2768  case SIDE_HIERARCHIC:
2769  return SIDE_DISCONTINUOUS;
2770 
2771  default:
2772  libmesh_error_msg("Unknown FE Family " << Utility::enum_to_string(fe_type.family));
2773  }
2774 }
2775 
2776 
2777 } // namespace libMesh
Real(* shape_deriv_ptr)(const FEType fet, const Elem *elem, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level)
Typedef for pointer to a function that returns FE shape function derivative values.
Definition: fe_interface.h:648
static unsigned int n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:530
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
OStreamProxy err
FEFamily family
The type of finite element.
Definition: fe_type.h:221
static Point map(unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p)
This is now deprecated; use FEMap::map instead.
Definition: fe_interface.C:677
bool need_derivative()
Check whether derivatives should be computed or not.
ElemType
Defines an enum for geometric element types.
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:355
static shape_ptr shape_function(const unsigned int dim, const FEType &fe_t, const ElemType t)
static unsigned int ifem_n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
static Real shape_deriv(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const unsigned int j, const Point &p)
static Real ifem_shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
unsigned int dim
class FEComputeData hides arbitrary data to be passed to and from children of FEBase through the FEIn...
We&#39;re using a class instead of a typedef to allow forward declarations and future flexibility...
const Point & p
Holds the point where the data are to be computed.
static void dofs_on_edge(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int e, std::vector< unsigned int > &di, const bool add_p_level=true)
Fills the vector di with the local degree of freedom indices associated with edge e of element elem A...
Definition: fe_interface.C:611
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2220
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
static FEFieldType field_type(const FEType &fe_type)
std::vector< Gradient > dshape
Storage for the computed shape derivative values.
unsigned int p_level() const
Definition: elem.h:3108
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:215
std::vector< std::vector< Real > > local_transform
Storage for local to global mapping at p.
The libMesh namespace provides an interface to certain functionality in the library.
static Point ifem_inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
void init()
Inits the output data to default values, provided the fields are correctly resized.
static unsigned int max_order(const FEType &fe_t, const ElemType &el_t)
static bool on_reference_element(const Point &p, const ElemType t, const Real eps=TOLERANCE)
Definition: fe_abstract.C:637
static void side_nodal_soln(const FEType &fe_t, const Elem *elem, const unsigned int side, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln, const bool add_p_level=true, const unsigned int vdim=1)
Build the nodal soln on one side from the (full) element soln.
Definition: fe_interface.C:650
This is the MeshBase class.
Definition: mesh_base.h:75
static bool orientation_dependent(const FEFamily &fe_family)
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to var...
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
FEInterface()
Empty constructor.
Definition: fe_interface.C:42
static void ifem_compute_data(const unsigned int dim, const FEType &fe_t, const Elem *elem, FEComputeData &data)
static bool extra_hanging_dofs(const FEType &fe_t)
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:272
static Real shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
Definition: fe_interface.C:760
static void compute_data(const unsigned int dim, const FEType &fe_t, const Elem *elem, FEComputeData &data)
Lets the appropriate child of FEBase compute the requested data for the input specified in data...
static bool is_hierarchic(const FEType &fe_type)
Returns whether or not the input FEType&#39;s higher-order shape functions are always hierarchic...
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libmesh_assert(ctx)
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
This is now deprecated; use FEMap::inverse_map instead.
Definition: fe_interface.C:692
static Point ifem_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p)
This is the base class for point locators.
static void ifem_nodal_soln(const unsigned int dim, const FEType &fe_t, const Elem *elem, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
static bool on_reference_element(const Point &p, const ElemType t, const Real eps=TOLERANCE)
Definition: fe_interface.C:751
static shape_deriv_ptr shape_deriv_function(const unsigned int dim, const FEType &fe_t, const ElemType t)
static unsigned int ifem_n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
static unsigned int n_vec_dim(const MeshBase &mesh, const FEType &fe_type)
std::vector< Number > shape
Storage for the computed shape function values.
static n_dofs_at_node_ptr n_dofs_at_node_function(const unsigned int dim, const FEType &fe_t)
Definition: fe_interface.C:458
static void compute_periodic_constraints(DofConstraints &constraints, DofMap &dof_map, const PeriodicBoundaries &boundaries, const MeshBase &mesh, const PointLocatorBase *point_locator, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for meshes with periodic boundary conditions) correspon...
Definition: fe_base.C:1878
std::string enum_to_string(const T e)
static FEContinuity get_continuity(const FEType &fe_type)
Returns the input FEType&#39;s FEContinuity based on the underlying FEFamily and potentially the Order...
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for non-conforming adapted meshes) corresponding to var...
static unsigned int n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
Definition: fe_interface.C:436
static Real shape_second_deriv(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const unsigned int j, const Point &p)
static unsigned int ifem_n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Real(* shape_ptr)(const FEType fe_t, const Elem *elem, const unsigned int i, const Point &p, const bool add_p_level)
Typedef for pointer to a function that returns FE shape function values.
Definition: fe_interface.h:541
unsigned int(* n_dofs_at_node_ptr)(const ElemType, const Order, const unsigned int)
Definition: fe_interface.h:151
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static shape_second_deriv_ptr shape_second_deriv_function(const unsigned int dim, const FEType &fe_t, const ElemType t)
virtual unsigned short dim() const =0
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
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:140
static bool is_InfFE_elem(const ElemType et)
Definition: fe_interface.C:52
Real(* shape_second_deriv_ptr)(const FEType fet, const Elem *elem, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level)
Typedef for pointer to a function that returns FE shape function second derivative values...
Definition: fe_interface.h:753
virtual bool infinite() const =0
FEFamily
defines an enum for finite element families.
static void nodal_soln(const unsigned int dim, const FEType &fe_t, const Elem *elem, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln, const bool add_p_level=true, const unsigned int vdim=1)
Build the nodal soln from the element soln.
Definition: fe_interface.C:625
static void compute_periodic_constraints(DofConstraints &constraints, DofMap &dof_map, const PeriodicBoundaries &boundaries, const MeshBase &mesh, const PointLocatorBase *point_locator, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for periodic boundary conditions) corresponding to vari...
static Real ifem_shape_deriv(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const unsigned int j, const Point &p)
virtual ElemType type() const =0
static void dofs_on_side(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int s, std::vector< unsigned int > &di, const bool add_p_level=true)
Fills the vector di with the local degree of freedom indices associated with side s of element elem A...
Definition: fe_interface.C:597
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
The constraint matrix storage format.
Definition: dof_map.h:108
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:117
static unsigned int ifem_n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
FEFieldType
defines an enum for finite element field types - i.e.