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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo