fe_interface.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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"
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 
30 namespace libMesh
31 {
32 
33 //------------------------------------------------------------
34 //FEInterface class members
36 {
37  libmesh_error_msg("ERROR: Do not define an object of this type.");
38 }
39 
40 
41 
42 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
43 
44 bool
46 {
47  return false;
48 }
49 
50 #else
51 
52 bool
54 {
55 
56  switch (et)
57  {
58  case INFEDGE2:
59  case INFQUAD4:
60  case INFQUAD6:
61  case INFHEX8:
62  case INFHEX16:
63  case INFHEX18:
64  case INFPRISM6:
65  case INFPRISM12:
66  {
67  return true;
68  }
69 
70  default:
71  {
72  return false;
73  }
74  }
75 }
76 
77 #endif //ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
78 
79 
80 
81 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
82 #define fe_family_switch(dim, func_and_args, prefix, suffix) \
83  do { \
84  switch (fe_t.family) \
85  { \
86  case CLOUGH: \
87  prefix FE<dim,CLOUGH>::func_and_args suffix \
88  case HERMITE: \
89  prefix FE<dim,HERMITE>::func_and_args suffix \
90  case HIERARCHIC: \
91  prefix FE<dim,HIERARCHIC>::func_and_args suffix \
92  case L2_HIERARCHIC: \
93  prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \
94  case LAGRANGE: \
95  prefix FE<dim,LAGRANGE>::func_and_args suffix \
96  case L2_LAGRANGE: \
97  prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \
98  case MONOMIAL: \
99  prefix FE<dim,MONOMIAL>::func_and_args suffix \
100  case SCALAR: \
101  prefix FE<dim,SCALAR>::func_and_args suffix \
102  case BERNSTEIN: \
103  prefix FE<dim,BERNSTEIN>::func_and_args suffix \
104  case SZABAB: \
105  prefix FE<dim,SZABAB>::func_and_args suffix \
106  case XYZ: \
107  prefix FEXYZ<dim>::func_and_args suffix \
108  case SUBDIVISION: \
109  libmesh_assert_equal_to (dim, 2); \
110  prefix FE<2,SUBDIVISION>::func_and_args suffix \
111  default: \
112  libmesh_error_msg("Unsupported family = " << fe_t.family); \
113  } \
114  } while (0)
115 
116 #define fe_family_with_vec_switch(dim, func_and_args, prefix, suffix) \
117  do { \
118  switch (fe_t.family) \
119  { \
120  case CLOUGH: \
121  prefix FE<dim,CLOUGH>::func_and_args suffix \
122  case HERMITE: \
123  prefix FE<dim,HERMITE>::func_and_args suffix \
124  case HIERARCHIC: \
125  prefix FE<dim,HIERARCHIC>::func_and_args suffix \
126  case L2_HIERARCHIC: \
127  prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \
128  case LAGRANGE: \
129  prefix FE<dim,LAGRANGE>::func_and_args suffix \
130  case LAGRANGE_VEC: \
131  prefix FELagrangeVec<dim>::func_and_args suffix \
132  case L2_LAGRANGE: \
133  prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \
134  case MONOMIAL: \
135  prefix FE<dim,MONOMIAL>::func_and_args suffix \
136  case SCALAR: \
137  prefix FE<dim,SCALAR>::func_and_args suffix \
138  case BERNSTEIN: \
139  prefix FE<dim,BERNSTEIN>::func_and_args suffix \
140  case SZABAB: \
141  prefix FE<dim,SZABAB>::func_and_args suffix \
142  case XYZ: \
143  prefix FEXYZ<dim>::func_and_args suffix \
144  case SUBDIVISION: \
145  libmesh_assert_equal_to (dim, 2); \
146  prefix FE<2,SUBDIVISION>::func_and_args suffix \
147  case NEDELEC_ONE: \
148  prefix FENedelecOne<dim>::func_and_args suffix \
149  default: \
150  libmesh_error_msg("Unsupported family = " << fe_t.family); \
151  } \
152  } while (0)
153 
154 #define fe_scalar_vec_error_switch(dim, func_and_args, prefix, suffix) \
155  do { \
156  switch (fe_t.family) \
157  { \
158  case CLOUGH: \
159  prefix FE<dim,CLOUGH>::func_and_args suffix \
160  case HERMITE: \
161  prefix FE<dim,HERMITE>::func_and_args suffix \
162  case HIERARCHIC: \
163  prefix FE<dim,HIERARCHIC>::func_and_args suffix \
164  case L2_HIERARCHIC: \
165  prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \
166  case LAGRANGE: \
167  prefix FE<dim,LAGRANGE>::func_and_args suffix \
168  case L2_LAGRANGE: \
169  prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \
170  case MONOMIAL: \
171  prefix FE<dim,MONOMIAL>::func_and_args suffix \
172  case SCALAR: \
173  prefix FE<dim,SCALAR>::func_and_args suffix \
174  case BERNSTEIN: \
175  prefix FE<dim,BERNSTEIN>::func_and_args suffix \
176  case SZABAB: \
177  prefix FE<dim,SZABAB>::func_and_args suffix \
178  case XYZ: \
179  prefix FEXYZ<dim>::func_and_args suffix \
180  case SUBDIVISION: \
181  libmesh_assert_equal_to (dim, 2); \
182  prefix FE<2,SUBDIVISION>::func_and_args suffix \
183  case LAGRANGE_VEC: \
184  case NEDELEC_ONE: \
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 = " << 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  case LAGRANGE_VEC: \
197  prefix FELagrangeVec<dim>::func_and_args suffix \
198  case NEDELEC_ONE: \
199  prefix FENedelecOne<dim>::func_and_args suffix \
200  case HERMITE: \
201  case HIERARCHIC: \
202  case L2_HIERARCHIC: \
203  case LAGRANGE: \
204  case L2_LAGRANGE: \
205  case MONOMIAL: \
206  case SCALAR: \
207  case BERNSTEIN: \
208  case SZABAB: \
209  case XYZ: \
210  case SUBDIVISION: \
211  libmesh_error_msg("Error: Can only request vector valued elements for RealGradient FEInterface::shape"); \
212  default: \
213  libmesh_error_msg("Unsupported family = " << fe_t.family); \
214  } \
215  } while (0)
216 
217 #else
218 #define fe_family_switch(dim, func_and_args, prefix, suffix) \
219  do { \
220  switch (fe_t.family) \
221  { \
222  case CLOUGH: \
223  prefix FE<dim,CLOUGH>::func_and_args suffix \
224  case HERMITE: \
225  prefix FE<dim,HERMITE>::func_and_args suffix \
226  case HIERARCHIC: \
227  prefix FE<dim,HIERARCHIC>::func_and_args suffix \
228  case L2_HIERARCHIC: \
229  prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \
230  case LAGRANGE: \
231  prefix FE<dim,LAGRANGE>::func_and_args suffix \
232  case L2_LAGRANGE: \
233  prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \
234  case MONOMIAL: \
235  prefix FE<dim,MONOMIAL>::func_and_args suffix \
236  case SCALAR: \
237  prefix FE<dim,SCALAR>::func_and_args suffix \
238  case XYZ: \
239  prefix FEXYZ<dim>::func_and_args suffix \
240  case SUBDIVISION: \
241  libmesh_assert_equal_to (dim, 2); \
242  prefix FE<2,SUBDIVISION>::func_and_args suffix \
243  default: \
244  libmesh_error_msg("Unsupported family = " << fe_t.family); \
245  } \
246  } while (0)
247 
248 #define fe_family_with_vec_switch(dim, func_and_args, prefix, suffix) \
249  do { \
250  switch (fe_t.family) \
251  { \
252  case CLOUGH: \
253  prefix FE<dim,CLOUGH>::func_and_args suffix \
254  case HERMITE: \
255  prefix FE<dim,HERMITE>::func_and_args suffix \
256  case HIERARCHIC: \
257  prefix FE<dim,HIERARCHIC>::func_and_args suffix \
258  case L2_HIERARCHIC: \
259  prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \
260  case LAGRANGE: \
261  prefix FE<dim,LAGRANGE>::func_and_args suffix \
262  case LAGRANGE_VEC: \
263  prefix FELagrangeVec<dim>::func_and_args suffix \
264  case L2_LAGRANGE: \
265  prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \
266  case MONOMIAL: \
267  prefix FE<dim,MONOMIAL>::func_and_args suffix \
268  case SCALAR: \
269  prefix FE<dim,SCALAR>::func_and_args suffix \
270  case XYZ: \
271  prefix FEXYZ<dim>::func_and_args suffix \
272  case SUBDIVISION: \
273  libmesh_assert_equal_to (dim, 2); \
274  prefix FE<2,SUBDIVISION>::func_and_args suffix \
275  case NEDELEC_ONE: \
276  prefix FENedelecOne<dim>::func_and_args suffix \
277  default: \
278  libmesh_error_msg("Unsupported family = " << fe_t.family); \
279  } \
280  } while (0)
281 
282 #define fe_scalar_vec_error_switch(dim, func_and_args, prefix, suffix) \
283  do { \
284  switch (fe_t.family) \
285  { \
286  case CLOUGH: \
287  prefix FE<dim,CLOUGH>::func_and_args suffix \
288  case HERMITE: \
289  prefix FE<dim,HERMITE>::func_and_args suffix \
290  case HIERARCHIC: \
291  prefix FE<dim,HIERARCHIC>::func_and_args suffix \
292  case L2_HIERARCHIC: \
293  prefix FE<dim,L2_HIERARCHIC>::func_and_args suffix \
294  case LAGRANGE: \
295  prefix FE<dim,LAGRANGE>::func_and_args suffix \
296  case L2_LAGRANGE: \
297  prefix FE<dim,L2_LAGRANGE>::func_and_args suffix \
298  case MONOMIAL: \
299  prefix FE<dim,MONOMIAL>::func_and_args suffix \
300  case SCALAR: \
301  prefix FE<dim,SCALAR>::func_and_args suffix \
302  case XYZ: \
303  prefix FEXYZ<dim>::func_and_args suffix \
304  case SUBDIVISION: \
305  libmesh_assert_equal_to (dim, 2); \
306  prefix FE<2,SUBDIVISION>::func_and_args suffix \
307  case LAGRANGE_VEC: \
308  case NEDELEC_ONE: \
309  libmesh_error_msg("Error: Can only request scalar valued elements for Real FEInterface::func_and_args"); \
310  default: \
311  libmesh_error_msg("Unsupported family = " << fe_t.family); \
312  } \
313  } while (0)
314 
315 
316 #define fe_vector_scalar_error_switch(dim, func_and_args, prefix, suffix) \
317  do { \
318  switch (fe_t.family) \
319  { \
320  case LAGRANGE_VEC: \
321  prefix FELagrangeVec<dim>::func_and_args suffix \
322  case NEDELEC_ONE: \
323  prefix FENedelecOne<dim>::func_and_args suffix \
324  case HERMITE: \
325  case HIERARCHIC: \
326  case L2_HIERARCHIC: \
327  case LAGRANGE: \
328  case L2_LAGRANGE: \
329  case MONOMIAL: \
330  case SCALAR: \
331  case XYZ: \
332  case SUBDIVISION: \
333  libmesh_error_msg("Error: Can only request vector valued elements for RealGradient FEInterface::func_and_args"); \
334  default: \
335  libmesh_error_msg("Unsupported family = " << fe_t.family); \
336  } \
337  } while (0)
338 #endif
339 
340 
341 #define fe_switch(func_and_args) \
342  do { \
343  switch (dim) \
344  { \
345  /* 0D */ \
346  case 0: \
347  fe_family_switch (0, func_and_args, return, ;); \
348  /* 1D */ \
349  case 1: \
350  fe_family_switch (1, func_and_args, return, ;); \
351  /* 2D */ \
352  case 2: \
353  fe_family_switch (2, func_and_args, return, ;); \
354  /* 3D */ \
355  case 3: \
356  fe_family_switch (3, func_and_args, return, ;); \
357  default: \
358  libmesh_error_msg("Invalid dim = " << dim); \
359  } \
360  } while (0)
361 
362 #define fe_with_vec_switch(func_and_args) \
363  do { \
364  switch (dim) \
365  { \
366  /* 0D */ \
367  case 0: \
368  fe_family_with_vec_switch (0, func_and_args, return, ;); \
369  /* 1D */ \
370  case 1: \
371  fe_family_with_vec_switch (1, func_and_args, return, ;); \
372  /* 2D */ \
373  case 2: \
374  fe_family_with_vec_switch (2, func_and_args, return, ;); \
375  /* 3D */ \
376  case 3: \
377  fe_family_with_vec_switch (3, func_and_args, return, ;); \
378  default: \
379  libmesh_error_msg("Invalid dim = " << dim); \
380  } \
381  } while (0)
382 
383 
384 #define void_fe_switch(func_and_args) \
385  do { \
386  switch (dim) \
387  { \
388  /* 0D */ \
389  case 0: \
390  fe_family_switch (0, func_and_args, ;, ; return;); \
391  /* 1D */ \
392  case 1: \
393  fe_family_switch (1, func_and_args, ;, ; return;); \
394  /* 2D */ \
395  case 2: \
396  fe_family_switch (2, func_and_args, ;, ; return;); \
397  /* 3D */ \
398  case 3: \
399  fe_family_switch (3, func_and_args, ;, ; return;); \
400  default: \
401  libmesh_error_msg("Invalid dim = " << dim); \
402  } \
403  } while (0)
404 
405 #define void_fe_with_vec_switch(func_and_args) \
406  do { \
407  switch (dim) \
408  { \
409  /* 0D */ \
410  case 0: \
411  fe_family_with_vec_switch (0, func_and_args, ;, ; return;); \
412  /* 1D */ \
413  case 1: \
414  fe_family_with_vec_switch (1, func_and_args, ;, ; return;); \
415  /* 2D */ \
416  case 2: \
417  fe_family_with_vec_switch (2, func_and_args, ;, ; return;); \
418  /* 3D */ \
419  case 3: \
420  fe_family_with_vec_switch (3, func_and_args, ;, ; return;); \
421  default: \
422  libmesh_error_msg("Invalid dim = " << dim); \
423  } \
424  } while (0)
425 
426 
427 
428 unsigned int FEInterface::n_shape_functions(const unsigned int dim,
429  const FEType & fe_t,
430  const ElemType t)
431 {
432 
433 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
434  /*
435  * Since the FEType, stored in DofMap/(some System child), has to
436  * be the _same_ for InfFE and FE, we have to catch calls
437  * to infinite elements through the element type.
438  */
439 
440  if (is_InfFE_elem(t))
441  return ifem_n_shape_functions(dim, fe_t, t);
442 
443 #endif
444 
445  const Order o = fe_t.order;
446 
447  fe_with_vec_switch(n_shape_functions(t, o));
448 }
449 
450 
451 
452 
453 
454 unsigned int FEInterface::n_dofs(const unsigned int dim,
455  const FEType & fe_t,
456  const ElemType t)
457 {
458 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
459 
460  if (is_InfFE_elem(t))
461  return ifem_n_dofs(dim, fe_t, t);
462 
463 #endif
464 
465  const Order o = fe_t.order;
466 
467  fe_with_vec_switch(n_dofs(t, o));
468 }
469 
470 
471 
472 
473 unsigned int FEInterface::n_dofs_at_node(const unsigned int dim,
474  const FEType & fe_t,
475  const ElemType t,
476  const unsigned int n)
477 {
478 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
479 
480  if (is_InfFE_elem(t))
481  return ifem_n_dofs_at_node(dim, fe_t, t, n);
482 
483 #endif
484 
485  const Order o = fe_t.order;
486 
487  fe_with_vec_switch(n_dofs_at_node(t, o, n));
488 }
489 
490 
491 
494  const FEType & fe_t)
495 {
496  fe_with_vec_switch(n_dofs_at_node);
497 }
498 
499 
500 
501 
502 
503 
504 unsigned int FEInterface::n_dofs_per_elem(const unsigned int dim,
505  const FEType & fe_t,
506  const ElemType t)
507 {
508 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
509 
510  if (is_InfFE_elem(t))
511  return ifem_n_dofs_per_elem(dim, fe_t, t);
512 
513 #endif
514 
515  const Order o = fe_t.order;
516 
517  fe_with_vec_switch(n_dofs_per_elem(t, o));
518 }
519 
520 
521 
522 
523 void FEInterface::dofs_on_side(const Elem * const elem,
524  const unsigned int dim,
525  const FEType & fe_t,
526  unsigned int s,
527  std::vector<unsigned int> & di)
528 {
529  const Order o = fe_t.order;
530 
531  void_fe_with_vec_switch(dofs_on_side(elem, o, s, di));
532 }
533 
534 
535 
536 void FEInterface::dofs_on_edge(const Elem * const elem,
537  const unsigned int dim,
538  const FEType & fe_t,
539  unsigned int e,
540  std::vector<unsigned int> & di)
541 {
542  const Order o = fe_t.order;
543 
544  void_fe_with_vec_switch(dofs_on_edge(elem, o, e, di));
545 }
546 
547 
548 
549 
550 void FEInterface::nodal_soln(const unsigned int dim,
551  const FEType & fe_t,
552  const Elem * elem,
553  const std::vector<Number> & elem_soln,
554  std::vector<Number> & nodal_soln)
555 {
556 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
557 
558  if (is_InfFE_elem(elem->type()))
559  {
560  ifem_nodal_soln(dim, fe_t, elem, elem_soln, nodal_soln);
561  return;
562  }
563 
564 #endif
565 
566  const Order order = fe_t.order;
567 
568  void_fe_with_vec_switch(nodal_soln(elem, order, elem_soln, nodal_soln));
569 }
570 
571 
572 
573 
574 Point FEInterface::map(unsigned int dim,
575  const FEType & fe_t,
576  const Elem * elem,
577  const Point & p)
578 {
579 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
580  if (is_InfFE_elem(elem->type()))
581  return ifem_map(dim, fe_t, elem, p);
582 #endif
583  fe_with_vec_switch(map(elem, p));
584 }
585 
586 
587 
588 
589 
590 Point FEInterface::inverse_map (const unsigned int dim,
591  const FEType & fe_t,
592  const Elem * elem,
593  const Point & p,
594  const Real tolerance,
595  const bool secure)
596 {
597 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
598 
599  if (is_InfFE_elem(elem->type()))
600  return ifem_inverse_map(dim, fe_t, elem, p,tolerance, secure);
601 
602 #endif
603 
604  fe_with_vec_switch(inverse_map(elem, p, tolerance, secure));
605 }
606 
607 
608 
609 
610 void FEInterface::inverse_map (const unsigned int dim,
611  const FEType & fe_t,
612  const Elem * elem,
613  const std::vector<Point> & physical_points,
614  std::vector<Point> & reference_points,
615  const Real tolerance,
616  const bool secure)
617 {
618  const std::size_t n_pts = physical_points.size();
619 
620  // Resize the vector
621  reference_points.resize(n_pts);
622 
623  if (n_pts == 0)
624  {
625  libMesh::err << "WARNING: empty vector physical_points!"
626  << std::endl;
627  libmesh_here();
628  return;
629  }
630 
631 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
632 
633  if (is_InfFE_elem(elem->type()))
634  {
635  ifem_inverse_map(dim, fe_t, elem, physical_points, reference_points, tolerance, secure);
636  return;
637  // libmesh_not_implemented();
638  }
639 
640 #endif
641 
642  void_fe_with_vec_switch(inverse_map(elem, physical_points, reference_points, tolerance, secure));
643 }
644 
645 
646 
648  const ElemType t,
649  const Real eps)
650 {
651  return FEBase::on_reference_element(p,t,eps);
652 }
653 
654 
655 
656 
657 Real FEInterface::shape(const unsigned int dim,
658  const FEType & fe_t,
659  const ElemType t,
660  const unsigned int i,
661  const Point & p)
662 {
663 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
664 
665  if (is_InfFE_elem(t))
666  return ifem_shape(dim, fe_t, t, i, p);
667 
668 #endif
669 
670  const Order o = fe_t.order;
671 
672  fe_switch(shape(t,o,i,p));
673 }
674 
675 Real FEInterface::shape(const unsigned int dim,
676  const FEType & fe_t,
677  const Elem * elem,
678  const unsigned int i,
679  const Point & p)
680 {
681 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
682 
683  if (elem && is_InfFE_elem(elem->type()))
684  return ifem_shape(dim, fe_t, elem, i, p);
685 
686 #endif
687 
688  const Order o = fe_t.order;
689 
690  fe_switch(shape(elem,o,i,p));
691 }
692 
693 template<>
694 void FEInterface::shape<Real>(const unsigned int dim,
695  const FEType & fe_t,
696  const ElemType t,
697  const unsigned int i,
698  const Point & p,
699  Real & phi)
700 {
701 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
702 
703  if (is_InfFE_elem(t))
704  {
705  phi = ifem_shape(dim, fe_t, t, i, p);
706  return;
707  }
708 
709 #endif
710 
711  const Order o = fe_t.order;
712 
713  switch(dim)
714  {
715  case 0:
716  fe_scalar_vec_error_switch(0, shape(t,o,i,p), phi = , ; break;);
717  break;
718  case 1:
719  fe_scalar_vec_error_switch(1, shape(t,o,i,p), phi = , ; break;);
720  break;
721  case 2:
722  fe_scalar_vec_error_switch(2, shape(t,o,i,p), phi = , ; break;);
723  break;
724  case 3:
725  fe_scalar_vec_error_switch(3, shape(t,o,i,p), phi = , ; break;);
726  break;
727  default:
728  libmesh_error_msg("Invalid dimension = " << dim);
729  }
730 
731  return;
732 }
733 
734 template<>
735 void FEInterface::shape<Real>(const unsigned int dim,
736  const FEType & fe_t,
737  const Elem * elem,
738  const unsigned int i,
739  const Point & p,
740  Real & phi)
741 {
742 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
743 
744  if (elem && is_InfFE_elem(elem->type()))
745  {
746  phi = ifem_shape(dim, fe_t, elem, i, p);
747  return;
748  }
749 #endif
750 
751  const Order o = fe_t.order;
752 
753  switch(dim)
754  {
755  case 0:
756  fe_scalar_vec_error_switch(0, shape(elem,o,i,p), phi = , ; break;);
757  break;
758  case 1:
759  fe_scalar_vec_error_switch(1, shape(elem,o,i,p), phi = , ; break;);
760  break;
761  case 2:
762  fe_scalar_vec_error_switch(2, shape(elem,o,i,p), phi = , ; break;);
763  break;
764  case 3:
765  fe_scalar_vec_error_switch(3, shape(elem,o,i,p), phi = , ; break;);
766  break;
767  default:
768  libmesh_error_msg("Invalid dimension = " << dim);
769  }
770 
771  return;
772 }
773 
774 template<>
775 void FEInterface::shape<RealGradient>(const unsigned int dim,
776  const FEType & fe_t,
777  const ElemType t,
778  const unsigned int i,
779  const Point & p,
780  RealGradient & phi)
781 {
782  // This even does not handle infinite elements at all!?
783  const Order o = fe_t.order;
784 
785  switch(dim)
786  {
787  case 0:
788  fe_vector_scalar_error_switch(0, shape(t,o,i,p), phi = , ; break;);
789  break;
790  case 1:
791  fe_vector_scalar_error_switch(1, shape(t,o,i,p), phi = , ; break;);
792  break;
793  case 2:
794  fe_vector_scalar_error_switch(2, shape(t,o,i,p), phi = , ; break;);
795  break;
796  case 3:
797  fe_vector_scalar_error_switch(3, shape(t,o,i,p), phi = , ; break;);
798  break;
799  default:
800  libmesh_error_msg("Invalid dimension = " << dim);
801  }
802 
803  return;
804 }
805 
806 template<>
807 void FEInterface::shape<RealGradient>(const unsigned int dim,
808  const FEType & fe_t,
809  const Elem * elem,
810  const unsigned int i,
811  const Point & p,
812  RealGradient & phi)
813 {
814  const Order o = fe_t.order;
815 
816  switch(dim)
817  {
818  case 0:
819  fe_vector_scalar_error_switch(0, shape(elem,o,i,p), phi = , ; break;);
820  break;
821  case 1:
822  fe_vector_scalar_error_switch(1, shape(elem,o,i,p), phi = , ; break;);
823  break;
824  case 2:
825  fe_vector_scalar_error_switch(2, shape(elem,o,i,p), phi = , ; break;);
826  break;
827  case 3:
828  fe_vector_scalar_error_switch(3, shape(elem,o,i,p), phi = , ; break;);
829  break;
830  default:
831  libmesh_error_msg("Invalid dimension = " << dim);
832  }
833 
834  return;
835 }
836 
837 void FEInterface::compute_data(const unsigned int dim,
838  const FEType & fe_t,
839  const Elem * elem,
841 {
842 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
843 
844  if (elem && is_InfFE_elem(elem->type()))
845  {
846  data.init();
847  ifem_compute_data(dim, fe_t, elem, data);
848  return;
849  }
850 
851 #endif
852 
853  FEType p_refined = fe_t;
854  p_refined.order = static_cast<Order>(p_refined.order + elem->p_level());
855 
856  const unsigned int n_dof = n_dofs (dim, p_refined, elem->type());
857  const Point & p = data.p;
858  data.shape.resize(n_dof);
859 
860  // set default values for all the output fields
861  data.init();
862 
863  for (unsigned int n=0; n<n_dof; n++)
864  data.shape[n] = shape(dim, p_refined, elem, n, p);
865 
866  return;
867 }
868 
869 
870 
871 #ifdef LIBMESH_ENABLE_AMR
872 
874  DofMap & dof_map,
875  const unsigned int variable_number,
876  const Elem * elem)
877 {
878  libmesh_assert(elem);
879 
880  const FEType & fe_t = dof_map.variable_type(variable_number);
881 
882  switch (elem->dim())
883  {
884  case 0:
885  case 1:
886  {
887  // No constraints in 0D/1D.
888  return;
889  }
890 
891 
892  case 2:
893  {
894  switch (fe_t.family)
895  {
896  case CLOUGH:
898  dof_map,
899  variable_number,
900  elem); return;
901 
902  case HERMITE:
904  dof_map,
905  variable_number,
906  elem); return;
907 
908  case LAGRANGE:
910  dof_map,
911  variable_number,
912  elem); return;
913 
914  case HIERARCHIC:
916  dof_map,
917  variable_number,
918  elem); return;
919 
920  case L2_HIERARCHIC:
922  dof_map,
923  variable_number,
924  elem); return;
925 
926  case LAGRANGE_VEC:
928  dof_map,
929  variable_number,
930  elem); return;
931 
932 
933 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
934  case SZABAB:
936  dof_map,
937  variable_number,
938  elem); return;
939 
940  case BERNSTEIN:
942  dof_map,
943  variable_number,
944  elem); return;
945 
946 #endif
947  default:
948  return;
949  }
950  }
951 
952 
953  case 3:
954  {
955  switch (fe_t.family)
956  {
957  case HERMITE:
959  dof_map,
960  variable_number,
961  elem); return;
962 
963  case LAGRANGE:
965  dof_map,
966  variable_number,
967  elem); return;
968 
969  case HIERARCHIC:
971  dof_map,
972  variable_number,
973  elem); return;
974 
975  case L2_HIERARCHIC:
977  dof_map,
978  variable_number,
979  elem); return;
980 
981  case LAGRANGE_VEC:
983  dof_map,
984  variable_number,
985  elem); return;
986 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
987  case SZABAB:
989  dof_map,
990  variable_number,
991  elem); return;
992 
993  case BERNSTEIN:
995  dof_map,
996  variable_number,
997  elem); return;
998 
999 #endif
1000  default:
1001  return;
1002  }
1003  }
1004 
1005 
1006  default:
1007  libmesh_error_msg("Invalid dimension = " << elem->dim());
1008  }
1009 }
1010 
1011 #endif // #ifdef LIBMESH_ENABLE_AMR
1012 
1013 
1014 
1015 #ifdef LIBMESH_ENABLE_PERIODIC
1016 
1018  DofMap & dof_map,
1019  const PeriodicBoundaries & boundaries,
1020  const MeshBase & mesh,
1021  const PointLocatorBase * point_locator,
1022  const unsigned int variable_number,
1023  const Elem * elem)
1024 {
1025  // No element-specific optimizations currently exist
1027  dof_map,
1028  boundaries,
1029  mesh,
1030  point_locator,
1031  variable_number,
1032  elem);
1033 }
1034 
1035 #endif // #ifdef LIBMESH_ENABLE_PERIODIC
1036 
1037 
1038 
1039 unsigned int FEInterface::max_order(const FEType & fe_t,
1040  const ElemType & el_t)
1041 {
1042  // Yeah, I know, infinity is much larger than 11, but our
1043  // solvers don't seem to like high degree polynomials, and our
1044  // quadrature rules and number_lookups tables
1045  // need to go up higher.
1046  const unsigned int unlimited = 11;
1047 
1048  // If we used 0 as a default, then elements missing from this
1049  // table (e.g. infinite elements) would be considered broken.
1050  const unsigned int unknown = unlimited;
1051 
1052  switch (fe_t.family)
1053  {
1054  case LAGRANGE:
1055  case L2_LAGRANGE: // TODO: L2_LAGRANGE can have higher "max_order" than LAGRANGE
1056  case LAGRANGE_VEC:
1057  switch (el_t)
1058  {
1059  case EDGE2:
1060  case EDGE3:
1061  case EDGE4:
1062  return 3;
1063  case TRI3:
1064  case TRISHELL3:
1065  return 1;
1066  case TRI6:
1067  return 2;
1068  case QUAD4:
1069  case QUADSHELL4:
1070  return 1;
1071  case QUAD8:
1072  case QUADSHELL8:
1073  case QUAD9:
1074  return 2;
1075  case TET4:
1076  return 1;
1077  case TET10:
1078  return 2;
1079  case HEX8:
1080  return 1;
1081  case HEX20:
1082  case HEX27:
1083  return 2;
1084  case PRISM6:
1085  return 1;
1086  case PRISM15:
1087  case PRISM18:
1088  return 2;
1089  case PYRAMID5:
1090  return 1;
1091  case PYRAMID13:
1092  case PYRAMID14:
1093  return 2;
1094  default:
1095  return unknown;
1096  }
1097  break;
1098  case MONOMIAL:
1099  switch (el_t)
1100  {
1101  case EDGE2:
1102  case EDGE3:
1103  case EDGE4:
1104  case TRI3:
1105  case TRISHELL3:
1106  case TRI6:
1107  case QUAD4:
1108  case QUADSHELL4:
1109  case QUAD8:
1110  case QUADSHELL8:
1111  case QUAD9:
1112  case TET4:
1113  case TET10:
1114  case HEX8:
1115  case HEX20:
1116  case HEX27:
1117  case PRISM6:
1118  case PRISM15:
1119  case PRISM18:
1120  case PYRAMID5:
1121  case PYRAMID13:
1122  case PYRAMID14:
1123  return unlimited;
1124  default:
1125  return unknown;
1126  }
1127  break;
1128 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1129  case BERNSTEIN:
1130  switch (el_t)
1131  {
1132  case EDGE2:
1133  case EDGE3:
1134  case EDGE4:
1135  return unlimited;
1136  case TRI3:
1137  case TRISHELL3:
1138  return 1;
1139  case TRI6:
1140  return 6;
1141  case QUAD4:
1142  case QUADSHELL4:
1143  return 1;
1144  case QUAD8:
1145  case QUADSHELL8:
1146  case QUAD9:
1147  return unlimited;
1148  case TET4:
1149  return 1;
1150  case TET10:
1151  return 2;
1152  case HEX8:
1153  return 1;
1154  case HEX20:
1155  return 2;
1156  case HEX27:
1157  return 4;
1158  case PRISM6:
1159  case PRISM15:
1160  case PRISM18:
1161  case PYRAMID5:
1162  case PYRAMID13:
1163  case PYRAMID14:
1164  return 0;
1165  default:
1166  return unknown;
1167  }
1168  break;
1169  case SZABAB:
1170  switch (el_t)
1171  {
1172  case EDGE2:
1173  return 1;
1174  case EDGE3:
1175  case EDGE4:
1176  return 7;
1177  case TRI3:
1178  case TRISHELL3:
1179  return 1;
1180  case TRI6:
1181  return 7;
1182  case QUAD4:
1183  case QUADSHELL4:
1184  return 1;
1185  case QUAD8:
1186  case QUADSHELL8:
1187  case QUAD9:
1188  return 7;
1189  case TET4:
1190  case TET10:
1191  case HEX8:
1192  case HEX20:
1193  case HEX27:
1194  case PRISM6:
1195  case PRISM15:
1196  case PRISM18:
1197  case PYRAMID5:
1198  case PYRAMID13:
1199  case PYRAMID14:
1200  return 0;
1201  default:
1202  return unknown;
1203  }
1204  break;
1205 #endif
1206  case XYZ:
1207  switch (el_t)
1208  {
1209  case EDGE2:
1210  case EDGE3:
1211  case EDGE4:
1212  case TRI3:
1213  case TRISHELL3:
1214  case TRI6:
1215  case QUAD4:
1216  case QUADSHELL4:
1217  case QUAD8:
1218  case QUADSHELL8:
1219  case QUAD9:
1220  case TET4:
1221  case TET10:
1222  case HEX8:
1223  case HEX20:
1224  case HEX27:
1225  case PRISM6:
1226  case PRISM15:
1227  case PRISM18:
1228  case PYRAMID5:
1229  case PYRAMID13:
1230  case PYRAMID14:
1231  return unlimited;
1232  default:
1233  return unknown;
1234  }
1235  break;
1236  case CLOUGH:
1237  switch (el_t)
1238  {
1239  case EDGE2:
1240  case EDGE3:
1241  return 3;
1242  case EDGE4:
1243  case TRI3:
1244  case TRISHELL3:
1245  return 0;
1246  case TRI6:
1247  return 3;
1248  case QUAD4:
1249  case QUADSHELL4:
1250  case QUAD8:
1251  case QUADSHELL8:
1252  case QUAD9:
1253  case TET4:
1254  case TET10:
1255  case HEX8:
1256  case HEX20:
1257  case HEX27:
1258  case PRISM6:
1259  case PRISM15:
1260  case PRISM18:
1261  case PYRAMID5:
1262  case PYRAMID13:
1263  case PYRAMID14:
1264  return 0;
1265  default:
1266  return unknown;
1267  }
1268  break;
1269  case HERMITE:
1270  switch (el_t)
1271  {
1272  case EDGE2:
1273  case EDGE3:
1274  return unlimited;
1275  case EDGE4:
1276  case TRI3:
1277  case TRISHELL3:
1278  case TRI6:
1279  return 0;
1280  case QUAD4:
1281  case QUADSHELL4:
1282  return 3;
1283  case QUAD8:
1284  case QUADSHELL8:
1285  case QUAD9:
1286  return unlimited;
1287  case TET4:
1288  case TET10:
1289  return 0;
1290  case HEX8:
1291  return 3;
1292  case HEX20:
1293  case HEX27:
1294  return unlimited;
1295  case PRISM6:
1296  case PRISM15:
1297  case PRISM18:
1298  case PYRAMID5:
1299  case PYRAMID13:
1300  case PYRAMID14:
1301  return 0;
1302  default:
1303  return unknown;
1304  }
1305  break;
1306  case HIERARCHIC:
1307  switch (el_t)
1308  {
1309  case EDGE2:
1310  case EDGE3:
1311  case EDGE4:
1312  return unlimited;
1313  case TRI3:
1314  case TRISHELL3:
1315  return 1;
1316  case TRI6:
1317  return unlimited;
1318  case QUAD4:
1319  case QUADSHELL4:
1320  return 1;
1321  case QUAD8:
1322  case QUADSHELL8:
1323  case QUAD9:
1324  return unlimited;
1325  case TET4:
1326  case TET10:
1327  return 0;
1328  case HEX8:
1329  case HEX20:
1330  return 1;
1331  case HEX27:
1332  return unlimited;
1333  case PRISM6:
1334  case PRISM15:
1335  case PRISM18:
1336  case PYRAMID5:
1337  case PYRAMID13:
1338  case PYRAMID14:
1339  return 0;
1340  default:
1341  return unknown;
1342  }
1343  break;
1344  case L2_HIERARCHIC:
1345  switch (el_t)
1346  {
1347  case EDGE2:
1348  case EDGE3:
1349  case EDGE4:
1350  return unlimited;
1351  case TRI3:
1352  case TRISHELL3:
1353  return 1;
1354  case TRI6:
1355  return unlimited;
1356  case QUAD4:
1357  case QUADSHELL4:
1358  return 1;
1359  case QUAD8:
1360  case QUADSHELL8:
1361  case QUAD9:
1362  return unlimited;
1363  case TET4:
1364  case TET10:
1365  return 0;
1366  case HEX8:
1367  case HEX20:
1368  return 1;
1369  case HEX27:
1370  return unlimited;
1371  case PRISM6:
1372  case PRISM15:
1373  case PRISM18:
1374  case PYRAMID5:
1375  case PYRAMID13:
1376  case PYRAMID14:
1377  return 0;
1378  default:
1379  return unknown;
1380  }
1381  break;
1382  case SUBDIVISION:
1383  switch (el_t)
1384  {
1385  case TRI3SUBDIVISION:
1386  return unlimited;
1387  default:
1388  return unknown;
1389  }
1390  break;
1391  case NEDELEC_ONE:
1392  switch (el_t)
1393  {
1394  case TRI6:
1395  case QUAD8:
1396  case QUAD9:
1397  case HEX20:
1398  case HEX27:
1399  return 1;
1400  default:
1401  return 0;
1402  }
1403  break;
1404  default:
1405  return 0;
1406  break;
1407  }
1408 }
1409 
1410 
1411 
1413 {
1414  switch (fe_t.family)
1415  {
1416  case LAGRANGE:
1417  case L2_LAGRANGE:
1418  case MONOMIAL:
1419  case L2_HIERARCHIC:
1420  case XYZ:
1421  case SUBDIVISION:
1422  case LAGRANGE_VEC:
1423  case NEDELEC_ONE:
1424  return false;
1425  case CLOUGH:
1426  case HERMITE:
1427  case HIERARCHIC:
1428 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1429  case BERNSTEIN:
1430  case SZABAB:
1431 #endif
1432  default:
1433  return true;
1434  }
1435 }
1436 
1438 {
1439  return FEInterface::field_type(fe_type.family);
1440 }
1441 
1443 {
1444  switch (fe_family)
1445  {
1446  case LAGRANGE_VEC:
1447  case NEDELEC_ONE:
1448  return TYPE_VECTOR;
1449  default:
1450  return TYPE_SCALAR;
1451  }
1452 }
1453 
1454 unsigned int FEInterface::n_vec_dim (const MeshBase & mesh,
1455  const FEType & fe_type)
1456 {
1457  switch (fe_type.family)
1458  {
1459  //FIXME: We currently assume that the number of vector components is tied
1460  // to the mesh dimension. This will break for mixed-dimension meshes.
1461  case LAGRANGE_VEC:
1462  case NEDELEC_ONE:
1463  return mesh.mesh_dimension();
1464  default:
1465  return 1;
1466  }
1467 }
1468 
1469 } // namespace libMesh
static unsigned int n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:504
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)
Definition: fe_interface.C:536
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:179
FEFamily family
Definition: fe_type.h:204
static Point map(unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p)
Definition: fe_interface.C:574
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:454
static unsigned int ifem_n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
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)
Definition: fe_interface.C:523
Helper class used with FEInterface::compute_data().
Maps between boundary ids and PeriodicBoundaryBase objects.
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1792
The base class for all geometric element types.
Definition: elem.h:100
MeshBase & mesh
static FEFieldType field_type(const FEType &fe_type)
unsigned int p_level() const
Definition: elem.h:2555
OrderWrapper order
Definition: fe_type.h:198
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 unsigned int max_order(const FEType &fe_t, const ElemType &el_t)
static bool on_reference_element(const Point &p, const ElemType t, const Real eps=TOLERANCE)
Definition: fe_abstract.C:581
Base class for Mesh.
Definition: mesh_base.h:77
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
Definition: fe_interface.C:873
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:176
static void ifem_compute_data(const unsigned int dim, const FEType &fe_t, const Elem *elem, FEComputeData &data)
static bool extra_hanging_dofs(const FEType &fe_t)
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:428
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:657
static void compute_data(const unsigned int dim, const FEType &fe_t, const Elem *elem, FEComputeData &data)
Definition: fe_interface.C:837
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)
Definition: fe_interface.C:590
static Point ifem_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p)
OStreamProxy err(std::cerr)
static void ifem_nodal_soln(const unsigned int dim, const FEType &fe_t, const Elem *elem, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
static bool on_reference_element(const Point &p, const ElemType t, const Real eps=TOLERANCE)
Definition: fe_interface.C:647
static unsigned int ifem_n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
static unsigned int n_vec_dim(const MeshBase &mesh, const FEType &fe_type)
static n_dofs_at_node_ptr n_dofs_at_node_function(const unsigned int dim, const FEType &fe_t)
Definition: fe_interface.C:493
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)
Definition: fe_base.C:1650
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
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:473
static unsigned int ifem_n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
unsigned int(* n_dofs_at_node_ptr)(const ElemType, const Order, const unsigned int)
Definition: fe_interface.h:114
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned short dim() const =0
static bool is_InfFE_elem(const ElemType et)
Definition: fe_interface.C:45
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)
Definition: fe_interface.C:550
IterBase * data
Real size() const
Definition: type_vector.h:901
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)
virtual ElemType type() const =0
A geometric point in (x,y,z) space.
Definition: point.h:38
static unsigned int ifem_n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)