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