libMesh
fe_interface.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2019 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/elem.h"
23 #include "libmesh/fe.h"
24 #include "libmesh/fe_compute_data.h"
25 #include "libmesh/dof_map.h"
26 #include "libmesh/enum_fe_family.h"
27 #include "libmesh/enum_order.h"
28 #include "libmesh/enum_elem_type.h"
29 #include "libmesh/enum_to_string.h"
30 
31 namespace libMesh
32 {
33 
34 //------------------------------------------------------------
35 //FEInterface class members
36 FEInterface::FEInterface()
37 {
38  libmesh_error_msg("ERROR: Do not define an object of this type.");
39 }
40 
41 
42 
43 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
44 
45 bool
46 FEInterface::is_InfFE_elem(const ElemType)
47 {
48  return false;
49 }
50 
51 #else
52 
53 bool
54 FEInterface::is_InfFE_elem(const ElemType et)
55 {
56 
57  switch (et)
58  {
59  case INFEDGE2:
60  case INFQUAD4:
61  case INFQUAD6:
62  case INFHEX8:
63  case INFHEX16:
64  case INFHEX18:
65  case INFPRISM6:
66  case INFPRISM12:
67  {
68  return true;
69  }
70 
71  default:
72  {
73  return false;
74  }
75  }
76 }
77 
78 #endif //ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
79 
80 
81 
82 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
83 #define fe_family_switch(dim, func_and_args, prefix, suffix) \
84  do { \
85  switch (fe_t.family) \
86  { \
87  case CLOUGH: \
88  prefix FE<dim,CLOUGH>::func_and_args suffix \
89  case HERMITE: \
90  prefix FE<dim,HERMITE>::func_and_args suffix \
91  case HIERARCHIC: \
92  prefix FE<dim,HIERARCHIC>::func_and_args suffix \
93  case L2_HIERARCHIC: \
94  prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \
95  case LAGRANGE: \
96  prefix FE<dim,LAGRANGE>::func_and_args suffix \
97  case L2_LAGRANGE: \
98  prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \
99  case MONOMIAL: \
100  prefix FE<dim,MONOMIAL>::func_and_args suffix \
101  case SCALAR: \
102  prefix FE<dim,SCALAR>::func_and_args suffix \
103  case BERNSTEIN: \
104  prefix FE<dim,BERNSTEIN>::func_and_args suffix \
105  case SZABAB: \
106  prefix FE<dim,SZABAB>::func_and_args suffix \
107  case RATIONAL_BERNSTEIN: \
108  prefix FE<dim,RATIONAL_BERNSTEIN>::func_and_args suffix \
109  case XYZ: \
110  prefix FEXYZ<dim>::func_and_args suffix \
111  case SUBDIVISION: \
112  libmesh_assert_equal_to (dim, 2); \
113  prefix FE<2,SUBDIVISION>::func_and_args suffix \
114  default: \
115  libmesh_error_msg("Unsupported family = " << fe_t.family); \
116  } \
117  } while (0)
118 
119 #define fe_family_with_vec_switch(dim, func_and_args, prefix, suffix) \
120  do { \
121  switch (fe_t.family) \
122  { \
123  case CLOUGH: \
124  prefix FE<dim,CLOUGH>::func_and_args suffix \
125  case HERMITE: \
126  prefix FE<dim,HERMITE>::func_and_args suffix \
127  case HIERARCHIC: \
128  prefix FE<dim,HIERARCHIC>::func_and_args suffix \
129  case L2_HIERARCHIC: \
130  prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \
131  case LAGRANGE: \
132  prefix FE<dim,LAGRANGE>::func_and_args suffix \
133  case LAGRANGE_VEC: \
134  prefix FELagrangeVec<dim>::func_and_args suffix \
135  case L2_LAGRANGE: \
136  prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \
137  case MONOMIAL: \
138  prefix FE<dim,MONOMIAL>::func_and_args suffix \
139  case MONOMIAL_VEC: \
140  prefix FEMonomialVec<dim>::func_and_args suffix \
141  case SCALAR: \
142  prefix FE<dim,SCALAR>::func_and_args suffix \
143  case BERNSTEIN: \
144  prefix FE<dim,BERNSTEIN>::func_and_args suffix \
145  case SZABAB: \
146  prefix FE<dim,SZABAB>::func_and_args suffix \
147  case RATIONAL_BERNSTEIN: \
148  prefix FE<dim,RATIONAL_BERNSTEIN>::func_and_args suffix \
149  case XYZ: \
150  prefix FEXYZ<dim>::func_and_args suffix \
151  case SUBDIVISION: \
152  libmesh_assert_equal_to (dim, 2); \
153  prefix FE<2,SUBDIVISION>::func_and_args suffix \
154  case NEDELEC_ONE: \
155  prefix FENedelecOne<dim>::func_and_args suffix \
156  default: \
157  libmesh_error_msg("Unsupported family = " << fe_t.family); \
158  } \
159  } while (0)
160 
161 #define fe_scalar_vec_error_switch(dim, func_and_args, prefix, suffix) \
162  do { \
163  switch (fe_t.family) \
164  { \
165  case CLOUGH: \
166  prefix FE<dim,CLOUGH>::func_and_args suffix \
167  case HERMITE: \
168  prefix FE<dim,HERMITE>::func_and_args suffix \
169  case HIERARCHIC: \
170  prefix FE<dim,HIERARCHIC>::func_and_args suffix \
171  case L2_HIERARCHIC: \
172  prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \
173  case LAGRANGE: \
174  prefix FE<dim,LAGRANGE>::func_and_args suffix \
175  case L2_LAGRANGE: \
176  prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \
177  case MONOMIAL: \
178  prefix FE<dim,MONOMIAL>::func_and_args suffix \
179  case SCALAR: \
180  prefix FE<dim,SCALAR>::func_and_args suffix \
181  case BERNSTEIN: \
182  prefix FE<dim,BERNSTEIN>::func_and_args suffix \
183  case RATIONAL_BERNSTEIN: \
184  prefix FE<dim,RATIONAL_BERNSTEIN>::func_and_args suffix \
185  case SZABAB: \
186  prefix FE<dim,SZABAB>::func_and_args suffix \
187  case XYZ: \
188  prefix FEXYZ<dim>::func_and_args suffix \
189  case SUBDIVISION: \
190  libmesh_assert_equal_to (dim, 2); \
191  prefix FE<2,SUBDIVISION>::func_and_args suffix \
192  case LAGRANGE_VEC: \
193  case NEDELEC_ONE: \
194  case MONOMIAL_VEC: \
195  libmesh_error_msg("Error: Can only request scalar valued elements for Real FEInterface::func_and_args"); \
196  default: \
197  libmesh_error_msg("Unsupported family = " << fe_t.family); \
198  } \
199  } while (0)
200 
201 
202 #define fe_vector_scalar_error_switch(dim, func_and_args, prefix, suffix) \
203  do { \
204  switch (fe_t.family) \
205  { \
206  case LAGRANGE_VEC: \
207  prefix FELagrangeVec<dim>::func_and_args suffix \
208  case NEDELEC_ONE: \
209  prefix FENedelecOne<dim>::func_and_args suffix \
210  case MONOMIAL_VEC: \
211  prefix FEMonomialVec<dim>::func_and_args suffix \
212  case HERMITE: \
213  case HIERARCHIC: \
214  case L2_HIERARCHIC: \
215  case LAGRANGE: \
216  case L2_LAGRANGE: \
217  case MONOMIAL: \
218  case SCALAR: \
219  case BERNSTEIN: \
220  case SZABAB: \
221  case RATIONAL_BERNSTEIN: \
222  case XYZ: \
223  case SUBDIVISION: \
224  libmesh_error_msg("Error: Can only request vector valued elements for RealGradient FEInterface::shape"); \
225  default: \
226  libmesh_error_msg("Unsupported family = " << fe_t.family); \
227  } \
228  } while (0)
229 
230 #else
231 #define fe_family_switch(dim, func_and_args, prefix, suffix) \
232  do { \
233  switch (fe_t.family) \
234  { \
235  case CLOUGH: \
236  prefix FE<dim,CLOUGH>::func_and_args suffix \
237  case HERMITE: \
238  prefix FE<dim,HERMITE>::func_and_args suffix \
239  case HIERARCHIC: \
240  prefix FE<dim,HIERARCHIC>::func_and_args suffix \
241  case L2_HIERARCHIC: \
242  prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \
243  case LAGRANGE: \
244  prefix FE<dim,LAGRANGE>::func_and_args suffix \
245  case L2_LAGRANGE: \
246  prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \
247  case MONOMIAL: \
248  prefix FE<dim,MONOMIAL>::func_and_args suffix \
249  case SCALAR: \
250  prefix FE<dim,SCALAR>::func_and_args suffix \
251  case XYZ: \
252  prefix FEXYZ<dim>::func_and_args suffix \
253  case SUBDIVISION: \
254  libmesh_assert_equal_to (dim, 2); \
255  prefix FE<2,SUBDIVISION>::func_and_args suffix \
256  default: \
257  libmesh_error_msg("Unsupported family = " << fe_t.family); \
258  } \
259  } while (0)
260 
261 #define fe_family_with_vec_switch(dim, func_and_args, prefix, suffix) \
262  do { \
263  switch (fe_t.family) \
264  { \
265  case CLOUGH: \
266  prefix FE<dim,CLOUGH>::func_and_args suffix \
267  case HERMITE: \
268  prefix FE<dim,HERMITE>::func_and_args suffix \
269  case HIERARCHIC: \
270  prefix FE<dim,HIERARCHIC>::func_and_args suffix \
271  case L2_HIERARCHIC: \
272  prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \
273  case LAGRANGE: \
274  prefix FE<dim,LAGRANGE>::func_and_args suffix \
275  case LAGRANGE_VEC: \
276  prefix FELagrangeVec<dim>::func_and_args suffix \
277  case L2_LAGRANGE: \
278  prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \
279  case MONOMIAL: \
280  prefix FE<dim,MONOMIAL>::func_and_args suffix \
281  case MONOMIAL_VEC: \
282  prefix FEMonomialVec<dim>::func_and_args suffix \
283  case SCALAR: \
284  prefix FE<dim,SCALAR>::func_and_args suffix \
285  case XYZ: \
286  prefix FEXYZ<dim>::func_and_args suffix \
287  case SUBDIVISION: \
288  libmesh_assert_equal_to (dim, 2); \
289  prefix FE<2,SUBDIVISION>::func_and_args suffix \
290  case NEDELEC_ONE: \
291  prefix FENedelecOne<dim>::func_and_args suffix \
292  default: \
293  libmesh_error_msg("Unsupported family = " << fe_t.family); \
294  } \
295  } while (0)
296 
297 #define fe_scalar_vec_error_switch(dim, func_and_args, prefix, suffix) \
298  do { \
299  switch (fe_t.family) \
300  { \
301  case CLOUGH: \
302  prefix FE<dim,CLOUGH>::func_and_args suffix \
303  case HERMITE: \
304  prefix FE<dim,HERMITE>::func_and_args suffix \
305  case HIERARCHIC: \
306  prefix FE<dim,HIERARCHIC>::func_and_args suffix \
307  case L2_HIERARCHIC: \
308  prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \
309  case LAGRANGE: \
310  prefix FE<dim,LAGRANGE>::func_and_args suffix \
311  case L2_LAGRANGE: \
312  prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \
313  case MONOMIAL: \
314  prefix FE<dim,MONOMIAL>::func_and_args suffix \
315  case SCALAR: \
316  prefix FE<dim,SCALAR>::func_and_args suffix \
317  case XYZ: \
318  prefix FEXYZ<dim>::func_and_args suffix \
319  case SUBDIVISION: \
320  libmesh_assert_equal_to (dim, 2); \
321  prefix FE<2,SUBDIVISION>::func_and_args suffix \
322  case LAGRANGE_VEC: \
323  case NEDELEC_ONE: \
324  case MONOMIAL_VEC: \
325  libmesh_error_msg("Error: Can only request scalar valued elements for Real FEInterface::func_and_args"); \
326  default: \
327  libmesh_error_msg("Unsupported family = " << fe_t.family); \
328  } \
329  } while (0)
330 
331 
332 #define fe_vector_scalar_error_switch(dim, func_and_args, prefix, suffix) \
333  do { \
334  switch (fe_t.family) \
335  { \
336  case LAGRANGE_VEC: \
337  prefix FELagrangeVec<dim>::func_and_args suffix \
338  case NEDELEC_ONE: \
339  prefix FENedelecOne<dim>::func_and_args suffix \
340  case MONOMIAL_VEC: \
341  prefix FEMonomialVec<dim>::func_and_args suffix \
342  case HERMITE: \
343  case HIERARCHIC: \
344  case L2_HIERARCHIC: \
345  case LAGRANGE: \
346  case L2_LAGRANGE: \
347  case MONOMIAL: \
348  case SCALAR: \
349  case XYZ: \
350  case SUBDIVISION: \
351  libmesh_error_msg("Error: Can only request vector valued elements for RealGradient FEInterface::func_and_args"); \
352  default: \
353  libmesh_error_msg("Unsupported family = " << fe_t.family); \
354  } \
355  } while (0)
356 #endif
357 
358 
359 #define fe_switch(func_and_args) \
360  do { \
361  switch (dim) \
362  { \
363  /* 0D */ \
364  case 0: \
365  fe_family_switch (0, func_and_args, return, ;); \
366  /* 1D */ \
367  case 1: \
368  fe_family_switch (1, func_and_args, return, ;); \
369  /* 2D */ \
370  case 2: \
371  fe_family_switch (2, func_and_args, return, ;); \
372  /* 3D */ \
373  case 3: \
374  fe_family_switch (3, func_and_args, return, ;); \
375  default: \
376  libmesh_error_msg("Invalid dim = " << dim); \
377  } \
378  } while (0)
379 
380 #define fe_with_vec_switch(func_and_args) \
381  do { \
382  switch (dim) \
383  { \
384  /* 0D */ \
385  case 0: \
386  fe_family_with_vec_switch (0, func_and_args, return, ;); \
387  /* 1D */ \
388  case 1: \
389  fe_family_with_vec_switch (1, func_and_args, return, ;); \
390  /* 2D */ \
391  case 2: \
392  fe_family_with_vec_switch (2, func_and_args, return, ;); \
393  /* 3D */ \
394  case 3: \
395  fe_family_with_vec_switch (3, func_and_args, return, ;); \
396  default: \
397  libmesh_error_msg("Invalid dim = " << dim); \
398  } \
399  } while (0)
400 
401 
402 #define void_fe_switch(func_and_args) \
403  do { \
404  switch (dim) \
405  { \
406  /* 0D */ \
407  case 0: \
408  fe_family_switch (0, func_and_args, ;, ; return;); \
409  /* 1D */ \
410  case 1: \
411  fe_family_switch (1, func_and_args, ;, ; return;); \
412  /* 2D */ \
413  case 2: \
414  fe_family_switch (2, func_and_args, ;, ; return;); \
415  /* 3D */ \
416  case 3: \
417  fe_family_switch (3, func_and_args, ;, ; return;); \
418  default: \
419  libmesh_error_msg("Invalid dim = " << dim); \
420  } \
421  } while (0)
422 
423 #define void_fe_with_vec_switch(func_and_args) \
424  do { \
425  switch (dim) \
426  { \
427  /* 0D */ \
428  case 0: \
429  fe_family_with_vec_switch (0, func_and_args, ;, ; return;); \
430  /* 1D */ \
431  case 1: \
432  fe_family_with_vec_switch (1, func_and_args, ;, ; return;); \
433  /* 2D */ \
434  case 2: \
435  fe_family_with_vec_switch (2, func_and_args, ;, ; return;); \
436  /* 3D */ \
437  case 3: \
438  fe_family_with_vec_switch (3, func_and_args, ;, ; return;); \
439  default: \
440  libmesh_error_msg("Invalid dim = " << dim); \
441  } \
442  } while (0)
443 
444 
445 
446 unsigned int FEInterface::n_shape_functions(const unsigned int dim,
447  const FEType & fe_t,
448  const ElemType t)
449 {
450 
451 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
452  /*
453  * Since the FEType, stored in DofMap/(some System child), has to
454  * be the _same_ for InfFE and FE, we have to catch calls
455  * to infinite elements through the element type.
456  */
457 
458  if (is_InfFE_elem(t))
459  return ifem_n_shape_functions(dim, fe_t, t);
460 
461 #endif
462 
463  const Order o = fe_t.order;
464 
465  fe_with_vec_switch(n_shape_functions(t, o));
466 }
467 
468 
469 
470 
471 
472 unsigned int FEInterface::n_dofs(const unsigned int dim,
473  const FEType & fe_t,
474  const ElemType t)
475 {
476 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
477 
478  if (is_InfFE_elem(t))
479  return ifem_n_dofs(dim, fe_t, t);
480 
481 #endif
482 
483  const Order o = fe_t.order;
484 
485  fe_with_vec_switch(n_dofs(t, o));
486 }
487 
488 
489 
490 
491 unsigned int
492 FEInterface::n_dofs (const unsigned int dim,
493  const FEType & fe_t,
494  const Elem * elem)
495 {
496  FEType p_refined_fe_t = fe_t;
497  p_refined_fe_t.order = static_cast<Order>(p_refined_fe_t.order + elem->p_level());
498  return FEInterface::n_dofs(dim, p_refined_fe_t, elem->type());
499 }
500 
501 
502 
503 unsigned int FEInterface::n_dofs_at_node(const unsigned int dim,
504  const FEType & fe_t,
505  const ElemType t,
506  const unsigned int n)
507 {
508 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
509 
510  if (is_InfFE_elem(t))
511  return ifem_n_dofs_at_node(dim, fe_t, t, n);
512 
513 #endif
514 
515  const Order o = fe_t.order;
516 
517  fe_with_vec_switch(n_dofs_at_node(t, o, n));
518 }
519 
520 
521 
522 FEInterface::n_dofs_at_node_ptr
523 FEInterface::n_dofs_at_node_function(const unsigned int dim,
524  const FEType & fe_t)
525 {
526  fe_with_vec_switch(n_dofs_at_node);
527 }
528 
529 
530 
531 
532 
533 
534 unsigned int FEInterface::n_dofs_per_elem(const unsigned int dim,
535  const FEType & fe_t,
536  const ElemType t)
537 {
538 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
539 
540  if (is_InfFE_elem(t))
541  return ifem_n_dofs_per_elem(dim, fe_t, t);
542 
543 #endif
544 
545  const Order o = fe_t.order;
546 
547  fe_with_vec_switch(n_dofs_per_elem(t, o));
548 }
549 
550 
551 
552 
553 void FEInterface::dofs_on_side(const Elem * const elem,
554  const unsigned int dim,
555  const FEType & fe_t,
556  unsigned int s,
557  std::vector<unsigned int> & di)
558 {
559  const Order o = fe_t.order;
560 
561  void_fe_with_vec_switch(dofs_on_side(elem, o, s, di));
562 }
563 
564 
565 
566 void FEInterface::dofs_on_edge(const Elem * const elem,
567  const unsigned int dim,
568  const FEType & fe_t,
569  unsigned int e,
570  std::vector<unsigned int> & di)
571 {
572  const Order o = fe_t.order;
573 
574  void_fe_with_vec_switch(dofs_on_edge(elem, o, e, di));
575 }
576 
577 
578 
579 
580 void FEInterface::nodal_soln(const unsigned int dim,
581  const FEType & fe_t,
582  const Elem * elem,
583  const std::vector<Number> & elem_soln,
584  std::vector<Number> & nodal_soln)
585 {
586 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
587 
588  if (is_InfFE_elem(elem->type()))
589  {
590  ifem_nodal_soln(dim, fe_t, elem, elem_soln, nodal_soln);
591  return;
592  }
593 
594 #endif
595 
596  const Order order = fe_t.order;
597 
598  void_fe_with_vec_switch(nodal_soln(elem, order, elem_soln, nodal_soln));
599 }
600 
601 
602 
603 
604 Point FEInterface::map(unsigned int dim,
605  const FEType & fe_t,
606  const Elem * elem,
607  const Point & p)
608 {
609 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
610  if (is_InfFE_elem(elem->type()))
611  return ifem_map(dim, fe_t, elem, p);
612 #endif
613  fe_with_vec_switch(map(elem, p));
614 }
615 
616 
617 
618 
619 
620 Point FEInterface::inverse_map (const unsigned int dim,
621  const FEType & fe_t,
622  const Elem * elem,
623  const Point & p,
624  const Real tolerance,
625  const bool secure)
626 {
627 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
628 
629  if (is_InfFE_elem(elem->type()))
630  return ifem_inverse_map(dim, fe_t, elem, p,tolerance, secure);
631 
632 #endif
633 
634  fe_with_vec_switch(inverse_map(elem, p, tolerance, secure));
635 }
636 
637 
638 
639 
640 void FEInterface::inverse_map (const unsigned int dim,
641  const FEType & fe_t,
642  const Elem * elem,
643  const std::vector<Point> & physical_points,
644  std::vector<Point> & reference_points,
645  const Real tolerance,
646  const bool secure)
647 {
648  const std::size_t n_pts = physical_points.size();
649 
650  // Resize the vector
651  reference_points.resize(n_pts);
652 
653  if (n_pts == 0)
654  {
655  libMesh::err << "WARNING: empty vector physical_points!"
656  << std::endl;
657  libmesh_here();
658  return;
659  }
660 
661 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
662 
663  if (is_InfFE_elem(elem->type()))
664  {
665  ifem_inverse_map(dim, fe_t, elem, physical_points, reference_points, tolerance, secure);
666  return;
667  // libmesh_not_implemented();
668  }
669 
670 #endif
671 
672  void_fe_with_vec_switch(inverse_map(elem, physical_points, reference_points, tolerance, secure));
673 }
674 
675 
676 
677 bool FEInterface::on_reference_element(const Point & p,
678  const ElemType t,
679  const Real eps)
680 {
681  return FEBase::on_reference_element(p,t,eps);
682 }
683 
684 
685 
686 
687 Real FEInterface::shape(const unsigned int dim,
688  const FEType & fe_t,
689  const ElemType t,
690  const unsigned int i,
691  const Point & p)
692 {
693 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
694 
695  if (is_InfFE_elem(t))
696  return ifem_shape(dim, fe_t, t, i, p);
697 
698 #endif
699 
700  const Order o = fe_t.order;
701 
702  fe_switch(shape(t,o,i,p));
703 }
704 
705 Real FEInterface::shape(const unsigned int dim,
706  const FEType & fe_t,
707  const Elem * elem,
708  const unsigned int i,
709  const Point & p)
710 {
711 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
712 
713  if (elem && is_InfFE_elem(elem->type()))
714  return ifem_shape(dim, fe_t, elem, i, p);
715 
716 #endif
717 
718  const Order o = fe_t.order;
719 
720  fe_switch(shape(elem,o,i,p));
721 }
722 
723 template<>
724 void FEInterface::shape<Real>(const unsigned int dim,
725  const FEType & fe_t,
726  const ElemType t,
727  const unsigned int i,
728  const Point & p,
729  Real & phi)
730 {
731 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
732 
733  if (is_InfFE_elem(t))
734  {
735  phi = ifem_shape(dim, fe_t, t, i, p);
736  return;
737  }
738 
739 #endif
740 
741  const Order o = fe_t.order;
742 
743  switch(dim)
744  {
745  case 0:
746  fe_scalar_vec_error_switch(0, shape(t,o,i,p), phi = , ; break;);
747  break;
748  case 1:
749  fe_scalar_vec_error_switch(1, shape(t,o,i,p), phi = , ; break;);
750  break;
751  case 2:
752  fe_scalar_vec_error_switch(2, shape(t,o,i,p), phi = , ; break;);
753  break;
754  case 3:
755  fe_scalar_vec_error_switch(3, shape(t,o,i,p), phi = , ; break;);
756  break;
757  default:
758  libmesh_error_msg("Invalid dimension = " << dim);
759  }
760 
761  return;
762 }
763 
764 template<>
765 void FEInterface::shape<Real>(const unsigned int dim,
766  const FEType & fe_t,
767  const Elem * elem,
768  const unsigned int i,
769  const Point & p,
770  Real & phi)
771 {
772 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
773 
774  if (elem && is_InfFE_elem(elem->type()))
775  {
776  phi = ifem_shape(dim, fe_t, elem, i, p);
777  return;
778  }
779 #endif
780 
781  const Order o = fe_t.order;
782 
783  switch(dim)
784  {
785  case 0:
786  fe_scalar_vec_error_switch(0, shape(elem,o,i,p), phi = , ; break;);
787  break;
788  case 1:
789  fe_scalar_vec_error_switch(1, shape(elem,o,i,p), phi = , ; break;);
790  break;
791  case 2:
792  fe_scalar_vec_error_switch(2, shape(elem,o,i,p), phi = , ; break;);
793  break;
794  case 3:
795  fe_scalar_vec_error_switch(3, shape(elem,o,i,p), phi = , ; break;);
796  break;
797  default:
798  libmesh_error_msg("Invalid dimension = " << dim);
799  }
800 
801  return;
802 }
803 
804 template<>
805 void FEInterface::shape<RealGradient>(const unsigned int dim,
806  const FEType & fe_t,
807  const ElemType t,
808  const unsigned int i,
809  const Point & p,
810  RealGradient & phi)
811 {
812  // This even does not handle infinite elements at all!?
813  const Order o = fe_t.order;
814 
815  switch(dim)
816  {
817  case 0:
818  fe_vector_scalar_error_switch(0, shape(t,o,i,p), phi = , ; break;);
819  break;
820  case 1:
821  fe_vector_scalar_error_switch(1, shape(t,o,i,p), phi = , ; break;);
822  break;
823  case 2:
824  fe_vector_scalar_error_switch(2, shape(t,o,i,p), phi = , ; break;);
825  break;
826  case 3:
827  fe_vector_scalar_error_switch(3, shape(t,o,i,p), phi = , ; break;);
828  break;
829  default:
830  libmesh_error_msg("Invalid dimension = " << dim);
831  }
832 
833  return;
834 }
835 
836 
837 FEInterface::shape_ptr
838 FEInterface::shape_function(const unsigned int dim,
839  const FEType & fe_t)
840 {
841  fe_switch(shape);
842 }
843 
844 
845 Real FEInterface::shape_deriv(const unsigned int dim,
846  const FEType & fe_t,
847  const ElemType t,
848  const unsigned int i,
849  const unsigned int j,
850  const Point & p)
851 {
852  libmesh_assert_greater (dim,j);
853 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
854 
855  if (is_InfFE_elem(t)){
856  return ifem_shape_deriv(dim, fe_t, t, i, j, p);
857  }
858 
859 #endif
860 
861  const Order o = fe_t.order;
862 
863  switch(dim)
864  {
865  case 0:
866  fe_family_switch (0, shape_deriv(t, o, i, j, p), return , ;);
867  case 1:
868  fe_family_switch (1, shape_deriv(t, o, i, j, p), return , ;);
869  case 2:
870  fe_family_switch (2, shape_deriv(t, o, i, j, p), return , ;);
871  case 3:
872  fe_family_switch (3, shape_deriv(t, o, i, j, p), return , ;);
873  default:
874  libmesh_error_msg("Invalid dimension = " << dim);
875  }
876  return 0;
877 }
878 
879 
880 Real FEInterface::shape_deriv(const unsigned int dim,
881  const FEType & fe_t,
882  const Elem * elem,
883  const unsigned int i,
884  const unsigned int j,
885  const Point & p)
886 {
887  libmesh_assert_greater (dim,j);
888 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
889 
890  if (elem->infinite()){
891  return ifem_shape_deriv(dim, fe_t, elem, i, j, p);
892  }
893 
894 #endif
895 
896  const Order o = fe_t.order;
897 
898  switch(dim)
899  {
900  case 0:
901  fe_family_switch (0, shape_deriv(elem, o, i, j, p), return , ;);
902  case 1:
903  fe_family_switch (1, shape_deriv(elem, o, i, j, p), return , ;);
904  case 2:
905  fe_family_switch (2, shape_deriv(elem, o, i, j, p), return , ;);
906  case 3:
907  fe_family_switch (3, shape_deriv(elem, o, i, j, p), return , ;);
908  default:
909  libmesh_error_msg("Invalid dimension = " << dim);
910  }
911  return 0;
912 }
913 
914 
915 FEInterface::shape_deriv_ptr
916 FEInterface::shape_deriv_function(const unsigned int dim,
917  const FEType & fe_t)
918 {
919  fe_switch(shape_deriv);
920 }
921 
922 
923 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
924 Real FEInterface::shape_second_deriv(const unsigned int dim,
925  const FEType & fe_t,
926  const ElemType t,
927  const unsigned int i,
928  const unsigned int j,
929  const Point & p)
930 {
931  libmesh_assert_greater (dim*(dim-1),j);
932 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
933  if (is_InfFE_elem(t))
934  libmesh_not_implemented();
935 #endif
936 
937  const Order o = fe_t.order;
938 
939  switch(dim)
940  {
941  case 0:
942  fe_family_switch (0, shape_second_deriv(t, o, i, j, p), return , ;);
943  case 1:
944  fe_family_switch (1, shape_second_deriv(t, o, i, j, p), return , ;);
945  case 2:
946  fe_family_switch (2, shape_second_deriv(t, o, i, j, p), return , ;);
947  case 3:
948  fe_family_switch (3, shape_second_deriv(t, o, i, j, p), return , ;);
949  default:
950  libmesh_error_msg("Invalid dimension = " << dim);
951  }
952  return 0;
953 }
954 
955 
956 Real FEInterface::shape_second_deriv(const unsigned int dim,
957  const FEType & fe_t,
958  const Elem * elem,
959  const unsigned int i,
960  const unsigned int j,
961  const Point & p)
962 {
963  libmesh_assert_greater (dim,j);
964 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
965  if (elem->infinite())
966  libmesh_not_implemented();
967 #endif
968 
969  const Order o = fe_t.order;
970 
971  switch(dim)
972  {
973  case 0:
974  fe_family_switch (0, shape_second_deriv(elem, o, i, j, p), return , ;);
975  case 1:
976  fe_family_switch (1, shape_second_deriv(elem, o, i, j, p), return , ;);
977  case 2:
978  fe_family_switch (2, shape_second_deriv(elem, o, i, j, p), return , ;);
979  case 3:
980  fe_family_switch (3, shape_second_deriv(elem, o, i, j, p), return , ;);
981  default:
982  libmesh_error_msg("Invalid dimension = " << dim);
983  }
984  return 0;
985 }
986 
987 
988 FEInterface::shape_second_deriv_ptr
989 FEInterface::shape_second_deriv_function(const unsigned int dim,
990  const FEType & fe_t)
991 {
992  fe_switch(shape_second_deriv);
993 }
994 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
995 
996 
997 template<>
998 void FEInterface::shape<RealGradient>(const unsigned int dim,
999  const FEType & fe_t,
1000  const Elem * elem,
1001  const unsigned int i,
1002  const Point & p,
1003  RealGradient & phi)
1004 {
1005  const Order o = fe_t.order;
1006 
1007  switch(dim)
1008  {
1009  case 0:
1010  fe_vector_scalar_error_switch(0, shape(elem,o,i,p), phi = , ; break;);
1011  break;
1012  case 1:
1013  fe_vector_scalar_error_switch(1, shape(elem,o,i,p), phi = , ; break;);
1014  break;
1015  case 2:
1016  fe_vector_scalar_error_switch(2, shape(elem,o,i,p), phi = , ; break;);
1017  break;
1018  case 3:
1019  fe_vector_scalar_error_switch(3, shape(elem,o,i,p), phi = , ; break;);
1020  break;
1021  default:
1022  libmesh_error_msg("Invalid dimension = " << dim);
1023  }
1024 
1025  return;
1026 }
1027 
1028 void FEInterface::compute_data(const unsigned int dim,
1029  const FEType & fe_t,
1030  const Elem * elem,
1031  FEComputeData & data)
1032 {
1033 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1034 
1035  if (elem && is_InfFE_elem(elem->type()))
1036  {
1037  data.init();
1038  ifem_compute_data(dim, fe_t, elem, data);
1039  return;
1040  }
1041 
1042 #endif
1043 
1044  const unsigned int n_dof = n_dofs (dim, fe_t, elem);
1045  const Point & p = data.p;
1046  data.shape.resize(n_dof);
1047 
1048  if (data.need_derivative())
1049  {
1050  data.dshape.resize(n_dof);
1051  data.local_transform.resize(dim);
1052 
1053  for (unsigned int d=0; d<dim; d++)
1054  data.local_transform[d].resize(dim);
1055 
1056  UniquePtr<FEBase> fe (FEBase::build(dim, fe_t));
1057  std::vector<Point> pt(1);
1058  pt[0]=p;
1059  fe->get_dphideta(); // to compute the map
1060  fe->reinit(elem, &pt);
1061 
1062  // compute the reference->physical map.
1063  data.local_transform[0][0] = fe->get_dxidx()[0];
1064  if (dim > 1)
1065  {
1066  data.local_transform[1][0] = fe->get_detadx()[0];
1067  data.local_transform[1][1] = fe->get_detady()[0];
1068  data.local_transform[0][1] = fe->get_dxidy()[0];
1069  if (dim > 2)
1070  {
1071  data.local_transform[2][0] = fe->get_dzetadx()[0];
1072  data.local_transform[2][1] = fe->get_dzetady()[0];
1073  data.local_transform[2][2] = fe->get_dzetadz()[0];
1074  data.local_transform[1][2] = fe->get_detadz()[0];
1075  data.local_transform[0][2] = fe->get_dxidz()[0];
1076  }
1077  }
1078  }
1079 
1080  // set default values for all the output fields
1081  data.init();
1082 
1083  for (unsigned int n=0; n<n_dof; n++)
1084  {
1085  // Here we pass the original fe_t object. Additional p-levels
1086  // (if any) are handled internally by the shape() and
1087  // shape_deriv() functions since they have access to the elem
1088  // pointer. Note that we are already using the n_dof value
1089  // appropriate to the elevated p-level.
1090  data.shape[n] = shape(dim, fe_t, elem, n, p);
1091  if (data.need_derivative())
1092  {
1093  for (unsigned int j=0; j<dim; j++)
1094  data.dshape[n](j) = shape_deriv(dim, fe_t, elem, n, j, p);
1095  }
1096  }
1097 }
1098 
1099 
1100 
1101 #ifdef LIBMESH_ENABLE_AMR
1102 
1103 void FEInterface::compute_constraints (DofConstraints & constraints,
1104  DofMap & dof_map,
1105  const unsigned int variable_number,
1106  const Elem * elem)
1107 {
1108  libmesh_assert(elem);
1109 
1110  const FEType & fe_t = dof_map.variable_type(variable_number);
1111 
1112  switch (elem->dim())
1113  {
1114  case 0:
1115  case 1:
1116  {
1117  // No constraints in 0D/1D.
1118  return;
1119  }
1120 
1121 
1122  case 2:
1123  {
1124  switch (fe_t.family)
1125  {
1126  case CLOUGH:
1127  FE<2,CLOUGH>::compute_constraints (constraints,
1128  dof_map,
1129  variable_number,
1130  elem); return;
1131 
1132  case HERMITE:
1133  FE<2,HERMITE>::compute_constraints (constraints,
1134  dof_map,
1135  variable_number,
1136  elem); return;
1137 
1138  case LAGRANGE:
1139  FE<2,LAGRANGE>::compute_constraints (constraints,
1140  dof_map,
1141  variable_number,
1142  elem); return;
1143 
1144  case HIERARCHIC:
1145  FE<2,HIERARCHIC>::compute_constraints (constraints,
1146  dof_map,
1147  variable_number,
1148  elem); return;
1149 
1150  case L2_HIERARCHIC:
1151  FE<2,L2_HIERARCHIC>::compute_constraints (constraints,
1152  dof_map,
1153  variable_number,
1154  elem); return;
1155 
1156  case LAGRANGE_VEC:
1157  FE<2,LAGRANGE_VEC>::compute_constraints (constraints,
1158  dof_map,
1159  variable_number,
1160  elem); return;
1161 
1162 
1163 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1164  case SZABAB:
1165  FE<2,SZABAB>::compute_constraints (constraints,
1166  dof_map,
1167  variable_number,
1168  elem); return;
1169 
1170  case BERNSTEIN:
1171  FE<2,BERNSTEIN>::compute_constraints (constraints,
1172  dof_map,
1173  variable_number,
1174  elem); return;
1175 
1176  case RATIONAL_BERNSTEIN:
1177  FE<2,RATIONAL_BERNSTEIN>::compute_constraints (constraints,
1178  dof_map,
1179  variable_number,
1180  elem); return;
1181 
1182 #endif
1183  default:
1184  return;
1185  }
1186  }
1187 
1188 
1189  case 3:
1190  {
1191  switch (fe_t.family)
1192  {
1193  case HERMITE:
1194  FE<3,HERMITE>::compute_constraints (constraints,
1195  dof_map,
1196  variable_number,
1197  elem); return;
1198 
1199  case LAGRANGE:
1200  FE<3,LAGRANGE>::compute_constraints (constraints,
1201  dof_map,
1202  variable_number,
1203  elem); return;
1204 
1205  case HIERARCHIC:
1206  FE<3,HIERARCHIC>::compute_constraints (constraints,
1207  dof_map,
1208  variable_number,
1209  elem); return;
1210 
1211  case L2_HIERARCHIC:
1212  FE<3,L2_HIERARCHIC>::compute_constraints (constraints,
1213  dof_map,
1214  variable_number,
1215  elem); return;
1216 
1217  case LAGRANGE_VEC:
1218  FE<3,LAGRANGE_VEC>::compute_constraints (constraints,
1219  dof_map,
1220  variable_number,
1221  elem); return;
1222 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1223  case SZABAB:
1224  FE<3,SZABAB>::compute_constraints (constraints,
1225  dof_map,
1226  variable_number,
1227  elem); return;
1228 
1229  case BERNSTEIN:
1230  FE<3,BERNSTEIN>::compute_constraints (constraints,
1231  dof_map,
1232  variable_number,
1233  elem); return;
1234 
1235  case RATIONAL_BERNSTEIN:
1236  FE<3,RATIONAL_BERNSTEIN>::compute_constraints (constraints,
1237  dof_map,
1238  variable_number,
1239  elem); return;
1240 
1241 #endif
1242  default:
1243  return;
1244  }
1245  }
1246 
1247 
1248  default:
1249  libmesh_error_msg("Invalid dimension = " << elem->dim());
1250  }
1251 }
1252 
1253 #endif // #ifdef LIBMESH_ENABLE_AMR
1254 
1255 
1256 
1257 #ifdef LIBMESH_ENABLE_PERIODIC
1258 
1259 void FEInterface::compute_periodic_constraints (DofConstraints & constraints,
1260  DofMap & dof_map,
1261  const PeriodicBoundaries & boundaries,
1262  const MeshBase & mesh,
1263  const PointLocatorBase * point_locator,
1264  const unsigned int variable_number,
1265  const Elem * elem)
1266 {
1267  // No element-specific optimizations currently exist
1268  FEBase::compute_periodic_constraints (constraints,
1269  dof_map,
1270  boundaries,
1271  mesh,
1272  point_locator,
1273  variable_number,
1274  elem);
1275 }
1276 
1277 #endif // #ifdef LIBMESH_ENABLE_PERIODIC
1278 
1279 
1280 
1281 unsigned int FEInterface::max_order(const FEType & fe_t,
1282  const ElemType & el_t)
1283 {
1284  // Yeah, I know, infinity is much larger than 11, but our
1285  // solvers don't seem to like high degree polynomials, and our
1286  // quadrature rules and number_lookups tables
1287  // need to go up higher.
1288  const unsigned int unlimited = 11;
1289 
1290  // If we used 0 as a default, then elements missing from this
1291  // table (e.g. infinite elements) would be considered broken.
1292  const unsigned int unknown = unlimited;
1293 
1294  switch (fe_t.family)
1295  {
1296  case LAGRANGE:
1297  case L2_LAGRANGE: // TODO: L2_LAGRANGE can have higher "max_order" than LAGRANGE
1298  case LAGRANGE_VEC:
1299  switch (el_t)
1300  {
1301  case EDGE2:
1302  case EDGE3:
1303  case EDGE4:
1304  return 3;
1305  case TRI3:
1306  case TRISHELL3:
1307  return 1;
1308  case TRI6:
1309  return 2;
1310  case QUAD4:
1311  case QUADSHELL4:
1312  return 1;
1313  case QUAD8:
1314  case QUADSHELL8:
1315  case QUAD9:
1316  return 2;
1317  case TET4:
1318  return 1;
1319  case TET10:
1320  return 2;
1321  case HEX8:
1322  return 1;
1323  case HEX20:
1324  case HEX27:
1325  return 2;
1326  case PRISM6:
1327  return 1;
1328  case PRISM15:
1329  case PRISM18:
1330  return 2;
1331  case PYRAMID5:
1332  return 1;
1333  case PYRAMID13:
1334  case PYRAMID14:
1335  return 2;
1336  default:
1337  return unknown;
1338  }
1339  break;
1340  case MONOMIAL:
1341  case MONOMIAL_VEC:
1342  switch (el_t)
1343  {
1344  case EDGE2:
1345  case EDGE3:
1346  case EDGE4:
1347  case TRI3:
1348  case TRISHELL3:
1349  case TRI6:
1350  case QUAD4:
1351  case QUADSHELL4:
1352  case QUAD8:
1353  case QUADSHELL8:
1354  case QUAD9:
1355  case TET4:
1356  case TET10:
1357  case HEX8:
1358  case HEX20:
1359  case HEX27:
1360  case PRISM6:
1361  case PRISM15:
1362  case PRISM18:
1363  case PYRAMID5:
1364  case PYRAMID13:
1365  case PYRAMID14:
1366  return unlimited;
1367  default:
1368  return unknown;
1369  }
1370  break;
1371 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1372  case BERNSTEIN:
1373  case RATIONAL_BERNSTEIN:
1374  switch (el_t)
1375  {
1376  case EDGE2:
1377  case EDGE3:
1378  case EDGE4:
1379  return unlimited;
1380  case TRI3:
1381  case TRISHELL3:
1382  return 1;
1383  case TRI6:
1384  return 6;
1385  case QUAD4:
1386  case QUADSHELL4:
1387  return 1;
1388  case QUAD8:
1389  case QUADSHELL8:
1390  case QUAD9:
1391  return unlimited;
1392  case TET4:
1393  return 1;
1394  case TET10:
1395  return 2;
1396  case HEX8:
1397  return 1;
1398  case HEX20:
1399  return 2;
1400  case HEX27:
1401  return 4;
1402  case PRISM6:
1403  case PRISM15:
1404  case PRISM18:
1405  case PYRAMID5:
1406  case PYRAMID13:
1407  case PYRAMID14:
1408  return 0;
1409  default:
1410  return unknown;
1411  }
1412  break;
1413  case SZABAB:
1414  switch (el_t)
1415  {
1416  case EDGE2:
1417  return 1;
1418  case EDGE3:
1419  case EDGE4:
1420  return 7;
1421  case TRI3:
1422  case TRISHELL3:
1423  return 1;
1424  case TRI6:
1425  return 7;
1426  case QUAD4:
1427  case QUADSHELL4:
1428  return 1;
1429  case QUAD8:
1430  case QUADSHELL8:
1431  case QUAD9:
1432  return 7;
1433  case TET4:
1434  case TET10:
1435  case HEX8:
1436  case HEX20:
1437  case HEX27:
1438  case PRISM6:
1439  case PRISM15:
1440  case PRISM18:
1441  case PYRAMID5:
1442  case PYRAMID13:
1443  case PYRAMID14:
1444  return 0;
1445  default:
1446  return unknown;
1447  }
1448  break;
1449 #endif
1450  case XYZ:
1451  switch (el_t)
1452  {
1453  case EDGE2:
1454  case EDGE3:
1455  case EDGE4:
1456  case TRI3:
1457  case TRISHELL3:
1458  case TRI6:
1459  case QUAD4:
1460  case QUADSHELL4:
1461  case QUAD8:
1462  case QUADSHELL8:
1463  case QUAD9:
1464  case TET4:
1465  case TET10:
1466  case HEX8:
1467  case HEX20:
1468  case HEX27:
1469  case PRISM6:
1470  case PRISM15:
1471  case PRISM18:
1472  case PYRAMID5:
1473  case PYRAMID13:
1474  case PYRAMID14:
1475  return unlimited;
1476  default:
1477  return unknown;
1478  }
1479  break;
1480  case CLOUGH:
1481  switch (el_t)
1482  {
1483  case EDGE2:
1484  case EDGE3:
1485  return 3;
1486  case EDGE4:
1487  case TRI3:
1488  case TRISHELL3:
1489  return 0;
1490  case TRI6:
1491  return 3;
1492  case QUAD4:
1493  case QUADSHELL4:
1494  case QUAD8:
1495  case QUADSHELL8:
1496  case QUAD9:
1497  case TET4:
1498  case TET10:
1499  case HEX8:
1500  case HEX20:
1501  case HEX27:
1502  case PRISM6:
1503  case PRISM15:
1504  case PRISM18:
1505  case PYRAMID5:
1506  case PYRAMID13:
1507  case PYRAMID14:
1508  return 0;
1509  default:
1510  return unknown;
1511  }
1512  break;
1513  case HERMITE:
1514  switch (el_t)
1515  {
1516  case EDGE2:
1517  case EDGE3:
1518  return unlimited;
1519  case EDGE4:
1520  case TRI3:
1521  case TRISHELL3:
1522  case TRI6:
1523  return 0;
1524  case QUAD4:
1525  case QUADSHELL4:
1526  return 3;
1527  case QUAD8:
1528  case QUADSHELL8:
1529  case QUAD9:
1530  return unlimited;
1531  case TET4:
1532  case TET10:
1533  return 0;
1534  case HEX8:
1535  return 3;
1536  case HEX20:
1537  case HEX27:
1538  return unlimited;
1539  case PRISM6:
1540  case PRISM15:
1541  case PRISM18:
1542  case PYRAMID5:
1543  case PYRAMID13:
1544  case PYRAMID14:
1545  return 0;
1546  default:
1547  return unknown;
1548  }
1549  break;
1550  case HIERARCHIC:
1551  switch (el_t)
1552  {
1553  case EDGE2:
1554  case EDGE3:
1555  case EDGE4:
1556  return unlimited;
1557  case TRI3:
1558  case TRISHELL3:
1559  return 1;
1560  case TRI6:
1561  return unlimited;
1562  case QUAD4:
1563  case QUADSHELL4:
1564  return 1;
1565  case QUAD8:
1566  case QUADSHELL8:
1567  case QUAD9:
1568  return unlimited;
1569  case TET4:
1570  case TET10:
1571  return 0;
1572  case HEX8:
1573  case HEX20:
1574  return 1;
1575  case HEX27:
1576  return unlimited;
1577  case PRISM6:
1578  case PRISM15:
1579  case PRISM18:
1580  case PYRAMID5:
1581  case PYRAMID13:
1582  case PYRAMID14:
1583  return 0;
1584  default:
1585  return unknown;
1586  }
1587  break;
1588  case L2_HIERARCHIC:
1589  switch (el_t)
1590  {
1591  case EDGE2:
1592  case EDGE3:
1593  case EDGE4:
1594  return unlimited;
1595  case TRI3:
1596  case TRISHELL3:
1597  return 1;
1598  case TRI6:
1599  return unlimited;
1600  case QUAD4:
1601  case QUADSHELL4:
1602  return 1;
1603  case QUAD8:
1604  case QUADSHELL8:
1605  case QUAD9:
1606  return unlimited;
1607  case TET4:
1608  case TET10:
1609  return 0;
1610  case HEX8:
1611  case HEX20:
1612  return 1;
1613  case HEX27:
1614  return unlimited;
1615  case PRISM6:
1616  case PRISM15:
1617  case PRISM18:
1618  case PYRAMID5:
1619  case PYRAMID13:
1620  case PYRAMID14:
1621  return 0;
1622  default:
1623  return unknown;
1624  }
1625  break;
1626  case SUBDIVISION:
1627  switch (el_t)
1628  {
1629  case TRI3SUBDIVISION:
1630  return unlimited;
1631  default:
1632  return unknown;
1633  }
1634  break;
1635  case NEDELEC_ONE:
1636  switch (el_t)
1637  {
1638  case TRI6:
1639  case QUAD8:
1640  case QUAD9:
1641  case HEX20:
1642  case HEX27:
1643  return 1;
1644  default:
1645  return 0;
1646  }
1647  break;
1648  default:
1649  return 0;
1650  break;
1651  }
1652 }
1653 
1654 
1655 
1656 bool FEInterface::extra_hanging_dofs(const FEType & fe_t)
1657 {
1658  switch (fe_t.family)
1659  {
1660  case LAGRANGE:
1661  case L2_LAGRANGE:
1662  case MONOMIAL:
1663  case MONOMIAL_VEC:
1664  case L2_HIERARCHIC:
1665  case XYZ:
1666  case SUBDIVISION:
1667  case LAGRANGE_VEC:
1668  case NEDELEC_ONE:
1669  return false;
1670  case CLOUGH:
1671  case HERMITE:
1672  case HIERARCHIC:
1673 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1674  case BERNSTEIN:
1675  case SZABAB:
1676  case RATIONAL_BERNSTEIN:
1677 #endif
1678  default:
1679  return true;
1680  }
1681 }
1682 
1683 FEFieldType FEInterface::field_type(const FEType & fe_type)
1684 {
1685  return FEInterface::field_type(fe_type.family);
1686 }
1687 
1688 FEFieldType FEInterface::field_type (const FEFamily & fe_family)
1689 {
1690  switch (fe_family)
1691  {
1692  case LAGRANGE_VEC:
1693  case NEDELEC_ONE:
1694  case MONOMIAL_VEC:
1695  return TYPE_VECTOR;
1696  default:
1697  return TYPE_SCALAR;
1698  }
1699 }
1700 
1701 unsigned int FEInterface::n_vec_dim (const MeshBase & mesh,
1702  const FEType & fe_type)
1703 {
1704  switch (fe_type.family)
1705  {
1706  //FIXME: We currently assume that the number of vector components is tied
1707  // to the mesh dimension. This will break for mixed-dimension meshes.
1708  case LAGRANGE_VEC:
1709  case NEDELEC_ONE:
1710  case MONOMIAL_VEC:
1711  return mesh.mesh_dimension();
1712  default:
1713  return 1;
1714  }
1715 }
1716 
1717 
1718 
1719 FEContinuity FEInterface::get_continuity(const FEType & fe_type)
1720 {
1721  switch (fe_type.family)
1722  {
1723  // Discontinuous elements
1724  case MONOMIAL:
1725  case MONOMIAL_VEC:
1726  case L2_HIERARCHIC:
1727  case L2_LAGRANGE:
1728  case XYZ:
1729  case SCALAR:
1730  return DISCONTINUOUS;
1731 
1732  // C0 elements
1733  case LAGRANGE:
1734  case HIERARCHIC:
1735  case BERNSTEIN:
1736  case SZABAB:
1737  case RATIONAL_BERNSTEIN:
1738  case INFINITE_MAP:
1739  case JACOBI_20_00:
1740  case JACOBI_30_00:
1741  case LEGENDRE:
1742  case LAGRANGE_VEC:
1743  return C_ZERO;
1744 
1745  // C1 elements
1746  case CLOUGH:
1747  case HERMITE:
1748  case SUBDIVISION:
1749  return C_ONE;
1750 
1751  case NEDELEC_ONE:
1752  return H_CURL;
1753 
1754  default:
1755  libmesh_error_msg("Unknown FE Family " << Utility::enum_to_string(fe_type.family));
1756  }
1757 }
1758 
1759 
1760 } // namespace libMesh
libMesh::DofConstraints
The constraint matrix storage format.
Definition: dof_map.h:105
libMesh::C_ONE
Definition: enum_fe_family.h:80
libMesh::HEX20
Definition: enum_elem_type.h:48
libMesh::FEInterface::shape_function
static shape_ptr shape_function(const unsigned int dim, const FEType &fe_t)
Definition: fe_interface.C:838
libMesh::PRISM6
Definition: enum_elem_type.h:50
libMesh::PeriodicBoundaries
We're using a class instead of a typedef to allow forward declarations and future flexibility.
Definition: periodic_boundaries.h:51
libMesh::FEGenericBase::compute_periodic_constraints
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:1692
libMesh::FEInterface::on_reference_element
static bool on_reference_element(const Point &p, const ElemType t, const Real eps=TOLERANCE)
Definition: fe_interface.C:677
libMesh::FEInterface::extra_hanging_dofs
static bool extra_hanging_dofs(const FEType &fe_t)
Definition: fe_interface.C:1656
libMesh::FEComputeData
class FEComputeData hides arbitrary data to be passed to and from children of FEBase through the FEIn...
Definition: fe_compute_data.h:51
libMesh::FEInterface::n_shape_functions
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:446
libMesh::FEInterface::get_continuity
static FEContinuity get_continuity(const FEType &fe_type)
Returns the input FEType's FEContinuity based on the underlying FEFamily and potentially the Order,...
Definition: fe_interface.C:1719
libMesh::CLOUGH
Definition: enum_fe_family.h:53
libMesh::FEInterface::n_vec_dim
static unsigned int n_vec_dim(const MeshBase &mesh, const FEType &fe_type)
Definition: fe_interface.C:1701
libMesh::FEType::family
FEFamily family
The type of finite element.
Definition: fe_type.h:203
libMesh::FEInterface::inverse_map
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
This is now deprecated; use FEMap::inverse_map instead.
Definition: fe_interface.C:620
libMesh::FEInterface::ifem_nodal_soln
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)
Definition: fe_interface_inf_fe.C:229
libMesh::L2_HIERARCHIC
Definition: enum_fe_family.h:40
libMesh::HEX8
Definition: enum_elem_type.h:47
libMesh::INFINITE_MAP
Definition: enum_fe_family.h:48
libMesh::EDGE4
Definition: enum_elem_type.h:37
libMesh::INFHEX8
Definition: enum_elem_type.h:60
libMesh::FEInterface::n_dofs_at_node_function
static n_dofs_at_node_ptr n_dofs_at_node_function(const unsigned int dim, const FEType &fe_t)
Definition: fe_interface.C:523
libMesh::INFQUAD4
Definition: enum_elem_type.h:58
libMesh::FEInterface::ifem_n_dofs
static unsigned int ifem_n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface_inf_fe.C:134
libMesh::HERMITE
Definition: enum_fe_family.h:54
libMesh::TYPE_SCALAR
Definition: enum_fe_family.h:93
libMesh::FEInterface::n_dofs_at_node
static unsigned int n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
Definition: fe_interface.C:503
libMesh::FEInterface::shape_second_deriv
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)
Definition: fe_interface.C:924
libMesh::FEInterface::ifem_n_dofs_per_elem
static unsigned int ifem_n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface_inf_fe.C:198
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
libMesh::Elem::dim
virtual unsigned short dim() const =0
libMesh::JACOBI_20_00
Definition: enum_fe_family.h:49
libMesh::TET10
Definition: enum_elem_type.h:46
libMesh::XYZ
Definition: enum_fe_family.h:46
libMesh::Order
Order
Definition: enum_order.h:40
libMesh::Elem::infinite
virtual bool infinite() const =0
libMesh::FEInterface::shape_deriv_function
static shape_deriv_ptr shape_deriv_function(const unsigned int dim, const FEType &fe_t)
Definition: fe_interface.C:916
libMesh::MONOMIAL_VEC
Definition: enum_fe_family.h:62
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2512
libMesh::FEInterface::is_InfFE_elem
static bool is_InfFE_elem(const ElemType et)
Definition: fe_interface.C:46
mesh
MeshBase & mesh
Definition: mesh_communication.C:1257
libMesh::FEInterface::shape_ptr
Real(* shape_ptr)(const Elem *elem, const Order o, const unsigned int i, const Point &p, const bool add_p_level)
Definition: fe_interface.h:296
libMesh::FEInterface::n_dofs_at_node_ptr
unsigned int(* n_dofs_at_node_ptr)(const ElemType, const Order, const unsigned int)
Definition: fe_interface.h:126
libMesh::INFEDGE2
Definition: enum_elem_type.h:57
dim
unsigned int dim
Definition: adaptivity_ex3.C:113
libMesh::FEInterface::nodal_soln
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)
Build the nodal soln from the element soln.
Definition: fe_interface.C:580
libMesh::FEInterface::compute_periodic_constraints
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...
Definition: fe_interface.C:1259
libMesh::VectorValue< Real >
libMesh::TET4
Definition: enum_elem_type.h:45
libMesh::FEInterface::ifem_shape_deriv
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)
Definition: fe_interface_inf_fe.C:893
libMesh::TypeVector::size
auto size() const -> decltype(std::norm(T()))
Definition: type_vector.h:944
libMesh::PRISM15
Definition: enum_elem_type.h:51
libMesh::SZABAB
Definition: enum_fe_family.h:44
libMesh::INFPRISM6
Definition: enum_elem_type.h:63
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::MeshBase
This is the MeshBase class.
Definition: mesh_base.h:78
libMesh::JACOBI_30_00
Definition: enum_fe_family.h:50
libMesh::HEX27
Definition: enum_elem_type.h:49
libMesh::UniquePtr
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:33
libMesh::FEInterface::n_dofs_per_elem
static unsigned int n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:534
libMesh::FEGenericBase::build
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libMesh::INFHEX18
Definition: enum_elem_type.h:62
libMesh::QUAD4
Definition: enum_elem_type.h:41
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
libMesh::FEInterface::dofs_on_edge
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)
Fills the vector di with the local degree of freedom indices associated with edge e of element elem A...
Definition: fe_interface.C:566
libMesh::DISCONTINUOUS
Definition: enum_fe_family.h:78
libMesh::TRI3
Definition: enum_elem_type.h:39
libMesh::H_CURL
Definition: enum_fe_family.h:81
libMesh::C_ZERO
Definition: enum_fe_family.h:79
libMesh::LEGENDRE
Definition: enum_fe_family.h:51
libMesh::FE::compute_constraints
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...
libMesh::FEInterface::ifem_n_dofs_at_node
static unsigned int ifem_n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
Definition: fe_interface_inf_fe.C:165
libMesh::FEInterface::n_dofs
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:472
libMesh::DofMap::variable_type
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1924
libMesh::RATIONAL_BERNSTEIN
Definition: enum_fe_family.h:64
libMesh::BERNSTEIN
Definition: enum_fe_family.h:43
libMesh::INFHEX16
Definition: enum_elem_type.h:61
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::FEInterface::shape_deriv_ptr
Real(* shape_deriv_ptr)(const Elem *elem, const Order o, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level)
Definition: fe_interface.h:345
libMesh::FEInterface::ifem_n_shape_functions
static unsigned int ifem_n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface_inf_fe.C:102
libMesh::FEInterface::map
static Point map(unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p)
This is now deprecated; use FEMap::map instead.
Definition: fe_interface.C:604
libMesh::HIERARCHIC
Definition: enum_fe_family.h:37
libMesh::MONOMIAL
Definition: enum_fe_family.h:39
libMesh::INFPRISM12
Definition: enum_elem_type.h:64
libMesh::FEInterface::ifem_shape
static Real ifem_shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
Definition: fe_interface_inf_fe.C:675
libMesh::FEInterface::dofs_on_side
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)
Fills the vector di with the local degree of freedom indices associated with side s of element elem A...
Definition: fe_interface.C:553
libMesh::FEAbstract::on_reference_element
static bool on_reference_element(const Point &p, const ElemType t, const Real eps=TOLERANCE)
Definition: fe_abstract.C:606
libMesh::TRI6
Definition: enum_elem_type.h:40
libMesh::FEInterface::shape_deriv
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)
Definition: fe_interface.C:845
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
libMesh::FEType::order
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:197
libMesh::FEFamily
FEFamily
Definition: enum_fe_family.h:34
libMesh::INFQUAD6
Definition: enum_elem_type.h:59
libMesh::QUADSHELL8
Definition: enum_elem_type.h:73
libMesh::PYRAMID5
Definition: enum_elem_type.h:53
libMesh::DofMap
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:176
libMesh::FEInterface::field_type
static FEFieldType field_type(const FEType &fe_type)
Definition: fe_interface.C:1683
libMesh::NEDELEC_ONE
Definition: enum_fe_family.h:61
libMesh::TYPE_VECTOR
Definition: enum_fe_family.h:94
libMesh::EDGE3
Definition: enum_elem_type.h:36
libMesh::QUADSHELL4
Definition: enum_elem_type.h:72
libMesh::FEInterface::shape_second_deriv_ptr
Real(* shape_second_deriv_ptr)(const Elem *elem, const Order o, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level)
Definition: fe_interface.h:409
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:100
libMesh::FEInterface::compute_constraints
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...
Definition: fe_interface.C:1103
data
IterBase * data
Ideally this private member data should have protected access.
Definition: variant_filter_iterator.h:337
libMesh::L2_LAGRANGE
Definition: enum_fe_family.h:41
libMesh::QUAD9
Definition: enum_elem_type.h:43
libMesh::FEInterface::ifem_inverse_map
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)
Definition: fe_interface_inf_fe.C:510
libMesh::LAGRANGE_VEC
Definition: enum_fe_family.h:60
libMesh::err
OStreamProxy err
libMesh::LAGRANGE
Definition: enum_fe_family.h:36
libMesh::SCALAR
Definition: enum_fe_family.h:58
libMesh::SUBDIVISION
Definition: enum_fe_family.h:55
libMesh::TRI3SUBDIVISION
Definition: enum_elem_type.h:69
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::FEContinuity
FEContinuity
Definition: enum_fe_family.h:77
libMesh::PRISM18
Definition: enum_elem_type.h:52
libMesh::TRISHELL3
Definition: enum_elem_type.h:71
libMesh::FEInterface::shape_second_deriv_function
static shape_second_deriv_ptr shape_second_deriv_function(const unsigned int dim, const FEType &fe_t)
Definition: fe_interface.C:989
libMesh::PointLocatorBase
This is the base class for point locators.
Definition: point_locator_base.h:62
libMesh::FEInterface::ifem_compute_data
static void ifem_compute_data(const unsigned int dim, const FEType &fe_t, const Elem *elem, FEComputeData &data)
Definition: fe_interface_inf_fe.C:918
libMesh::FEInterface::max_order
static unsigned int max_order(const FEType &fe_t, const ElemType &el_t)
Definition: fe_interface.C:1281
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::FEFieldType
FEFieldType
Definition: enum_fe_family.h:92
libMesh::PYRAMID13
Definition: enum_elem_type.h:54
libMesh::PYRAMID14
Definition: enum_elem_type.h:55
libMesh::EDGE2
Definition: enum_elem_type.h:35
libMesh::QUAD8
Definition: enum_elem_type.h:42
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::FEInterface::FEInterface
FEInterface()
Empty constructor.
Definition: fe_interface.C:36
libMesh::FEInterface::compute_data
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,...
Definition: fe_interface.C:1028
libMesh::FEInterface::shape
static Real shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
Definition: fe_interface.C:687
libMesh::FEInterface::ifem_map
static Point ifem_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p)
Definition: fe_interface_inf_fe.C:479
Generated on Sat Jan 25 2020 12:06:57 for libMesh by   doxygen 1.8.16