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