libMesh
inf_fe_static.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 // Local includes
20 #include "libmesh/libmesh_config.h"
21 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
22 #include "libmesh/inf_fe.h"
23 #include "libmesh/inf_fe_macro.h"
24 #include "libmesh/fe.h"
25 #include "libmesh/fe_interface.h"
26 #include "libmesh/fe_compute_data.h"
27 #include "libmesh/elem.h"
28 
29 namespace libMesh
30 {
31 
32 
33 // ------------------------------------------------------------
34 // InfFE class static member initialization
35 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
37 
38 #ifdef DEBUG
39 
40 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
42 
43 
44 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
46 
47 #endif
48 
49 
50 
51 
52 // ------------------------------------------------------------
53 // InfFE static class members
54 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
55 unsigned int InfFE<Dim,T_radial,T_map>::n_dofs (const FEType & fet,
56  const ElemType inf_elem_type)
57 {
58  const ElemType base_et (Base::get_elem_type(inf_elem_type));
59 
60  if (Dim > 1)
61  return FEInterface::n_dofs(Dim-1, fet, base_et) * Radial::n_dofs(fet.radial_order);
62  else
63  return Radial::n_dofs(fet.radial_order);
64 }
65 
66 
67 
68 
69 
70 
71 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
73  const ElemType inf_elem_type,
74  const unsigned int n)
75 {
76  const ElemType base_et (Base::get_elem_type(inf_elem_type));
77 
78  unsigned int n_base, n_radial;
79  compute_node_indices(inf_elem_type, n, n_base, n_radial);
80 
81  // libMesh::out << "elem_type=" << inf_elem_type
82  // << ", fet.radial_order=" << fet.radial_order
83  // << ", n=" << n
84  // << ", n_radial=" << n_radial
85  // << ", n_base=" << n_base
86  // << std::endl;
87 
88  if (Dim > 1)
89  return FEInterface::n_dofs_at_node(Dim-1, fet, base_et, n_base)
90  * Radial::n_dofs_at_node(fet.radial_order, n_radial);
91  else
92  return Radial::n_dofs_at_node(fet.radial_order, n_radial);
93 }
94 
95 
96 
97 
98 
99 
100 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
102  const ElemType inf_elem_type)
103 {
104  const ElemType base_et (Base::get_elem_type(inf_elem_type));
105 
106  if (Dim > 1)
107  return FEInterface::n_dofs_per_elem(Dim-1, fet, base_et)
108  * Radial::n_dofs_per_elem(fet.radial_order);
109  else
110  return Radial::n_dofs_per_elem(fet.radial_order);
111 }
112 
113 
114 
115 
116 
117 
118 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
120  const Elem * /* elem */,
121  const std::vector<Number> & /* elem_soln */,
122  std::vector<Number> & nodal_soln)
123 {
124 #ifdef DEBUG
125  if (!_warned_for_nodal_soln)
126  {
127  libMesh::err << "WARNING: nodal_soln(...) does _not_ work for infinite elements." << std::endl
128  << " Will return an empty nodal solution. Use " << std::endl
129  << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!" << std::endl;
130  _warned_for_nodal_soln = true;
131  }
132 #endif
133 
134  /*
135  * In the base the infinite element couples to
136  * conventional finite elements. To not destroy
137  * the values there, clear \p nodal_soln. This
138  * indicates to the user of \p nodal_soln to
139  * not use this result.
140  */
141  nodal_soln.clear();
142  libmesh_assert (nodal_soln.empty());
143  return;
144 }
145 
146 
147 
148 
149 
150 
151 
152 
153 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
155  const ElemType inf_elem_type,
156  const unsigned int i,
157  const Point & p)
158 {
159  libmesh_assert_not_equal_to (Dim, 0);
160 
161 #ifdef DEBUG
162  // this makes only sense when used for mapping
163  if ((T_radial != INFINITE_MAP) && !_warned_for_shape)
164  {
165  libMesh::err << "WARNING: InfFE<Dim,T_radial,T_map>::shape(...) does _not_" << std::endl
166  << " return the correct trial function! Use " << std::endl
167  << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!"
168  << std::endl;
169  _warned_for_shape = true;
170  }
171 #endif
172 
173  const ElemType base_et (Base::get_elem_type(inf_elem_type));
174  const Order o_radial (fet.radial_order);
175  const Real v (p(Dim-1));
176 
177  unsigned int i_base, i_radial;
178  compute_shape_indices(fet, inf_elem_type, i, i_base, i_radial);
179 
180  //TODO:[SP/DD] exp(ikr) is still missing here!
181  // but is it intended? It would be probably somehow nice, but than it would be Number, not Real !
182  // --> thus it would destroy the interface...
183  if (Dim > 1)
184  return FEInterface::shape(Dim-1, fet, base_et, i_base, p)
185  * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
187  else
188  return InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
190 }
191 
192 
193 
194 
195 
196 
197 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
199  const Elem * inf_elem,
200  const unsigned int i,
201  const Point & p)
202 {
203  libmesh_assert(inf_elem);
204  libmesh_assert_not_equal_to (Dim, 0);
205 
206 #ifdef DEBUG
207  // this makes only sense when used for mapping
208  if ((T_radial != INFINITE_MAP) && !_warned_for_shape)
209  {
210  libMesh::err << "WARNING: InfFE<Dim,T_radial,T_map>::shape(...) does _not_" << std::endl
211  << " return the correct trial function! Use " << std::endl
212  << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!"
213  << std::endl;
214  _warned_for_shape = true;
215  }
216 #endif
217 
218  const Order o_radial (fet.radial_order);
219  const Real v (p(Dim-1));
220  UniquePtr<const Elem> base_el (inf_elem->build_side_ptr(0));
221 
222  unsigned int i_base, i_radial;
223  compute_shape_indices(fet, inf_elem->type(), i, i_base, i_radial);
224 
225  if (Dim > 1)
226  return FEInterface::shape(Dim-1, fet, base_el.get(), i_base, p)
227  * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
229  else
230  return InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
232 }
233 
234 
235 
236 
237 
238 
239 
240 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
242  const Elem * inf_elem,
244 {
245  libmesh_assert(inf_elem);
246  libmesh_assert_not_equal_to (Dim, 0);
247 
248  const Order o_radial (fet.radial_order);
249  const Order radial_mapping_order (Radial::mapping_order());
250  const Point & p (data.p);
251  const Real v (p(Dim-1));
252  UniquePtr<const Elem> base_el (inf_elem->build_side_ptr(0));
253 
254  /*
255  * compute \p interpolated_dist containing the mapping-interpolated
256  * distance of the base point to the origin. This is the same
257  * for all shape functions. Set \p interpolated_dist to 0, it
258  * is added to.
259  */
260  Real interpolated_dist = 0.;
261  switch (Dim)
262  {
263  case 1:
264  {
265  libmesh_assert_equal_to (inf_elem->type(), INFEDGE2);
266  interpolated_dist = Point(inf_elem->point(0) - inf_elem->point(1)).norm();
267  break;
268  }
269 
270  case 2:
271  {
272  const unsigned int n_base_nodes = base_el->n_nodes();
273 
274  const Point origin = inf_elem->origin();
275  const Order base_mapping_order (base_el->default_order());
276  const ElemType base_mapping_elem_type (base_el->type());
277 
278  // interpolate the base nodes' distances
279  for (unsigned int n=0; n<n_base_nodes; n++)
280  interpolated_dist += Point(base_el->point(n) - origin).norm()
281  * FE<1,LAGRANGE>::shape (base_mapping_elem_type, base_mapping_order, n, p);
282  break;
283  }
284 
285  case 3:
286  {
287  const unsigned int n_base_nodes = base_el->n_nodes();
288 
289  const Point origin = inf_elem->origin();
290  const Order base_mapping_order (base_el->default_order());
291  const ElemType base_mapping_elem_type (base_el->type());
292 
293  // interpolate the base nodes' distances
294  for (unsigned int n=0; n<n_base_nodes; n++)
295  interpolated_dist += Point(base_el->point(n) - origin).norm()
296  * FE<2,LAGRANGE>::shape (base_mapping_elem_type, base_mapping_order, n, p);
297  break;
298  }
299 
300  default:
301  libmesh_error_msg("Unknown Dim = " << Dim);
302  }
303 
304 
305 
306 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
307 
308  // assumption on time-harmonic behavior
309  Number sign_i (0,1.);
310 
311  // the wave number
312  const Number wavenumber = 2. * libMesh::pi * data.frequency / data.speed;
313 
314  // the exponent for time-harmonic behavior
315  const Number exponent = sign_i /* imaginary unit */
316  * wavenumber /* k (can be complex) */
317  * interpolated_dist /* together with next line: */
318  * InfFE<Dim,INFINITE_MAP,T_map>::eval(v, radial_mapping_order, 1); /* phase(s,t,v) */
319 
320  const Number time_harmonic = exp(exponent); /* e^(sign*i*k*phase(s,t,v)) */
321 
322  /*
323  * compute \p shape for all dof in the element
324  */
325  if (Dim > 1)
326  {
327  const unsigned int n_dof = n_dofs (fet, inf_elem->type());
328  data.shape.resize(n_dof);
329 
330  for (unsigned int i=0; i<n_dof; i++)
331  {
332  // compute base and radial shape indices
333  unsigned int i_base, i_radial;
334  compute_shape_indices(fet, inf_elem->type(), i, i_base, i_radial);
335 
336  data.shape[i] = (InfFE<Dim,T_radial,T_map>::Radial::decay(v) /* (1.-v)/2. in 3D */
337  * FEInterface::shape(Dim-1, fet, base_el.get(), i_base, p) /* S_n(s,t) */
338  * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)) /* L_n(v) */
339  * time_harmonic; /* e^(sign*i*k*phase(s,t,v) */
340  }
341  }
342  else
343  libmesh_error_msg("compute_data() for 1-dimensional InfFE not implemented.");
344 
345 #else
346 
347  const Real speed = data.speed;
348 
349  /*
350  * This is quite weird: the phase is actually
351  * a measure how advanced the pressure is that
352  * we compute. In other words: the further away
353  * the node \p data.p is, the further we look into
354  * the future...
355  */
356  data.phase = interpolated_dist /* phase(s,t,v)/c */
357  * InfFE<Dim,INFINITE_MAP,T_map>::eval(v, radial_mapping_order, 1) / speed;
358 
359  if (Dim > 1)
360  {
361  const unsigned int n_dof = n_dofs (fet, inf_elem->type());
362  data.shape.resize(n_dof);
363 
364  for (unsigned int i=0; i<n_dof; i++)
365  {
366  // compute base and radial shape indices
367  unsigned int i_base, i_radial;
368  compute_shape_indices(fet, inf_elem->type(), i, i_base, i_radial);
369 
370  data.shape[i] = InfFE<Dim,T_radial,T_map>::Radial::decay(v) /* (1.-v)/2. in 3D */
371  * FEInterface::shape(Dim-1, fet, base_el.get(), i_base, p) /* S_n(s,t) */
372  * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial); /* L_n(v) */
373  }
374  }
375  else
376  libmesh_error_msg("compute_data() for 1-dimensional InfFE not implemented.");
377 
378 #endif
379 }
380 
381 
382 
383 
384 
385 
386 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
388  const unsigned int outer_node_index,
389  unsigned int & base_node,
390  unsigned int & radial_node)
391 {
392  switch (inf_elem_type)
393  {
394  case INFEDGE2:
395  {
396  libmesh_assert_less (outer_node_index, 2);
397  base_node = 0;
398  radial_node = outer_node_index;
399  return;
400  }
401 
402 
403  // linear base approximation, easy to determine
404  case INFQUAD4:
405  {
406  libmesh_assert_less (outer_node_index, 4);
407  base_node = outer_node_index % 2;
408  radial_node = outer_node_index / 2;
409  return;
410  }
411 
412  case INFPRISM6:
413  {
414  libmesh_assert_less (outer_node_index, 6);
415  base_node = outer_node_index % 3;
416  radial_node = outer_node_index / 3;
417  return;
418  }
419 
420  case INFHEX8:
421  {
422  libmesh_assert_less (outer_node_index, 8);
423  base_node = outer_node_index % 4;
424  radial_node = outer_node_index / 4;
425  return;
426  }
427 
428 
429  // higher order base approximation, more work necessary
430  case INFQUAD6:
431  {
432  switch (outer_node_index)
433  {
434  case 0:
435  case 1:
436  {
437  radial_node = 0;
438  base_node = outer_node_index;
439  return;
440  }
441 
442  case 2:
443  case 3:
444  {
445  radial_node = 1;
446  base_node = outer_node_index-2;
447  return;
448  }
449 
450  case 4:
451  {
452  radial_node = 0;
453  base_node = 2;
454  return;
455  }
456 
457  case 5:
458  {
459  radial_node = 1;
460  base_node = 2;
461  return;
462  }
463 
464  default:
465  libmesh_error_msg("Unrecognized outer_node_index = " << outer_node_index);
466  }
467  }
468 
469 
470  case INFHEX16:
471  case INFHEX18:
472  {
473  switch (outer_node_index)
474  {
475  case 0:
476  case 1:
477  case 2:
478  case 3:
479  {
480  radial_node = 0;
481  base_node = outer_node_index;
482  return;
483  }
484 
485  case 4:
486  case 5:
487  case 6:
488  case 7:
489  {
490  radial_node = 1;
491  base_node = outer_node_index-4;
492  return;
493  }
494 
495  case 8:
496  case 9:
497  case 10:
498  case 11:
499  {
500  radial_node = 0;
501  base_node = outer_node_index-4;
502  return;
503  }
504 
505  case 12:
506  case 13:
507  case 14:
508  case 15:
509  {
510  radial_node = 1;
511  base_node = outer_node_index-8;
512  return;
513  }
514 
515  case 16:
516  {
517  libmesh_assert_equal_to (inf_elem_type, INFHEX18);
518  radial_node = 0;
519  base_node = 8;
520  return;
521  }
522 
523  case 17:
524  {
525  libmesh_assert_equal_to (inf_elem_type, INFHEX18);
526  radial_node = 1;
527  base_node = 8;
528  return;
529  }
530 
531  default:
532  libmesh_error_msg("Unrecognized outer_node_index = " << outer_node_index);
533  }
534  }
535 
536 
537  case INFPRISM12:
538  {
539  switch (outer_node_index)
540  {
541  case 0:
542  case 1:
543  case 2:
544  {
545  radial_node = 0;
546  base_node = outer_node_index;
547  return;
548  }
549 
550  case 3:
551  case 4:
552  case 5:
553  {
554  radial_node = 1;
555  base_node = outer_node_index-3;
556  return;
557  }
558 
559  case 6:
560  case 7:
561  case 8:
562  {
563  radial_node = 0;
564  base_node = outer_node_index-3;
565  return;
566  }
567 
568  case 9:
569  case 10:
570  case 11:
571  {
572  radial_node = 1;
573  base_node = outer_node_index-6;
574  return;
575  }
576 
577  default:
578  libmesh_error_msg("Unrecognized outer_node_index = " << outer_node_index);
579  }
580  }
581 
582 
583  default:
584  libmesh_error_msg("ERROR: Bad infinite element type=" << inf_elem_type << ", node=" << outer_node_index);
585  }
586 }
587 
588 
589 
590 
591 
592 
593 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
595  const unsigned int outer_node_index,
596  unsigned int & base_node,
597  unsigned int & radial_node)
598 {
599  libmesh_assert_not_equal_to (inf_elem_type, INVALID_ELEM);
600 
601  static std::vector<unsigned int> _static_base_node_index;
602  static std::vector<unsigned int> _static_radial_node_index;
603 
604  /*
605  * fast counterpart to compute_node_indices(), uses local static buffers
606  * to store the index maps. The class member
607  * \p _compute_node_indices_fast_current_elem_type remembers
608  * the current element type.
609  *
610  * Note that there exist non-static members storing the
611  * same data. However, you never know what element type
612  * is currently used by the \p InfFE object, and what
613  * request is currently directed to the static \p InfFE
614  * members (which use \p compute_node_indices_fast()).
615  * So separate these.
616  *
617  * check whether the work for this elemtype has already
618  * been done. If so, use this index. Otherwise, refresh
619  * the buffer to this element type.
620  */
621  if (inf_elem_type==_compute_node_indices_fast_current_elem_type)
622  {
623  base_node = _static_base_node_index [outer_node_index];
624  radial_node = _static_radial_node_index[outer_node_index];
625  return;
626  }
627  else
628  {
629  // store the map for _all_ nodes for this element type
630  _compute_node_indices_fast_current_elem_type = inf_elem_type;
631 
632  unsigned int n_nodes = libMesh::invalid_uint;
633 
634  switch (inf_elem_type)
635  {
636  case INFEDGE2:
637  {
638  n_nodes = 2;
639  break;
640  }
641  case INFQUAD4:
642  {
643  n_nodes = 4;
644  break;
645  }
646  case INFQUAD6:
647  {
648  n_nodes = 6;
649  break;
650  }
651  case INFHEX8:
652  {
653  n_nodes = 8;
654  break;
655  }
656  case INFHEX16:
657  {
658  n_nodes = 16;
659  break;
660  }
661  case INFHEX18:
662  {
663  n_nodes = 18;
664  break;
665  }
666  case INFPRISM6:
667  {
668  n_nodes = 6;
669  break;
670  }
671  case INFPRISM12:
672  {
673  n_nodes = 12;
674  break;
675  }
676  default:
677  libmesh_error_msg("ERROR: Bad infinite element type=" << inf_elem_type << ", node=" << outer_node_index);
678  }
679 
680 
681  _static_base_node_index.resize (n_nodes);
682  _static_radial_node_index.resize(n_nodes);
683 
684  for (unsigned int n=0; n<n_nodes; n++)
685  compute_node_indices (inf_elem_type,
686  n,
687  _static_base_node_index [outer_node_index],
688  _static_radial_node_index[outer_node_index]);
689 
690  // and return for the specified node
691  base_node = _static_base_node_index [outer_node_index];
692  radial_node = _static_radial_node_index[outer_node_index];
693  return;
694  }
695 }
696 
697 
698 
699 
700 
701 
702 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
704  const ElemType inf_elem_type,
705  const unsigned int i,
706  unsigned int & base_shape,
707  unsigned int & radial_shape)
708 {
709 
710  /*
711  * An example is provided: the numbers in comments refer to
712  * a fictitious InfHex18. The numbers are chosen as exemplary
713  * values. There is currently no base approximation that
714  * requires this many dof's at nodes, sides, faces and in the element.
715  *
716  * the order of the shape functions is heavily related with the
717  * order the dofs are assigned in \p DofMap::distributed_dofs().
718  * Due to the infinite elements with higher-order base approximation,
719  * some more effort is necessary.
720  *
721  * numbering scheme:
722  * 1. all vertices in the base, assign node->n_comp() dofs to each vertex
723  * 2. all vertices further out: innermost loop: radial shapes,
724  * then the base approximation shapes
725  * 3. all side nodes in the base, assign node->n_comp() dofs to each side node
726  * 4. all side nodes further out: innermost loop: radial shapes,
727  * then the base approximation shapes
728  * 5. (all) face nodes in the base, assign node->n_comp() dofs to each face node
729  * 6. (all) face nodes further out: innermost loop: radial shapes,
730  * then the base approximation shapes
731  * 7. element-associated dof in the base
732  * 8. element-associated dof further out
733  */
734 
735  const unsigned int radial_order = static_cast<unsigned int>(fet.radial_order.get_order()); // 4
736  const unsigned int radial_order_p_one = radial_order+1; // 5
737 
738  const ElemType base_elem_type (Base::get_elem_type(inf_elem_type)); // QUAD9
739 
740  // assume that the number of dof is the same for all vertices
741  unsigned int n_base_vertices = libMesh::invalid_uint; // 4
742  const unsigned int n_base_vertex_dof = FEInterface::n_dofs_at_node (Dim-1, fet, base_elem_type, 0);// 2
743 
744  unsigned int n_base_side_nodes = libMesh::invalid_uint; // 4
745  unsigned int n_base_side_dof = libMesh::invalid_uint; // 3
746 
747  unsigned int n_base_face_nodes = libMesh::invalid_uint; // 1
748  unsigned int n_base_face_dof = libMesh::invalid_uint; // 5
749 
750  const unsigned int n_base_elem_dof = FEInterface::n_dofs_per_elem (Dim-1, fet, base_elem_type);// 9
751 
752 
753  switch (inf_elem_type)
754  {
755  case INFEDGE2:
756  {
757  n_base_vertices = 1;
758  n_base_side_nodes = 0;
759  n_base_face_nodes = 0;
760  n_base_side_dof = 0;
761  n_base_face_dof = 0;
762  break;
763  }
764 
765  case INFQUAD4:
766  {
767  n_base_vertices = 2;
768  n_base_side_nodes = 0;
769  n_base_face_nodes = 0;
770  n_base_side_dof = 0;
771  n_base_face_dof = 0;
772  break;
773  }
774 
775  case INFQUAD6:
776  {
777  n_base_vertices = 2;
778  n_base_side_nodes = 1;
779  n_base_face_nodes = 0;
780  n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
781  n_base_face_dof = 0;
782  break;
783  }
784 
785  case INFHEX8:
786  {
787  n_base_vertices = 4;
788  n_base_side_nodes = 0;
789  n_base_face_nodes = 0;
790  n_base_side_dof = 0;
791  n_base_face_dof = 0;
792  break;
793  }
794 
795  case INFHEX16:
796  {
797  n_base_vertices = 4;
798  n_base_side_nodes = 4;
799  n_base_face_nodes = 0;
800  n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
801  n_base_face_dof = 0;
802  break;
803  }
804 
805  case INFHEX18:
806  {
807  n_base_vertices = 4;
808  n_base_side_nodes = 4;
809  n_base_face_nodes = 1;
810  n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
811  n_base_face_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, 8);
812  break;
813  }
814 
815 
816  case INFPRISM6:
817  {
818  n_base_vertices = 3;
819  n_base_side_nodes = 0;
820  n_base_face_nodes = 0;
821  n_base_side_dof = 0;
822  n_base_face_dof = 0;
823  break;
824  }
825 
826  case INFPRISM12:
827  {
828  n_base_vertices = 3;
829  n_base_side_nodes = 3;
830  n_base_face_nodes = 0;
831  n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
832  n_base_face_dof = 0;
833  break;
834  }
835 
836  default:
837  libmesh_error_msg("Unrecognized inf_elem_type = " << inf_elem_type);
838  }
839 
840 
841  {
842  // these are the limits describing the intervals where the shape function lies
843  const unsigned int n_dof_at_base_vertices = n_base_vertices*n_base_vertex_dof; // 8
844  const unsigned int n_dof_at_all_vertices = n_dof_at_base_vertices*radial_order_p_one; // 40
845 
846  const unsigned int n_dof_at_base_sides = n_base_side_nodes*n_base_side_dof; // 12
847  const unsigned int n_dof_at_all_sides = n_dof_at_base_sides*radial_order_p_one; // 60
848 
849  const unsigned int n_dof_at_base_face = n_base_face_nodes*n_base_face_dof; // 5
850  const unsigned int n_dof_at_all_faces = n_dof_at_base_face*radial_order_p_one; // 25
851 
852 
853  // start locating the shape function
854  if (i < n_dof_at_base_vertices) // range of i: 0..7
855  {
856  // belongs to vertex in the base
857  radial_shape = 0;
858  base_shape = i;
859  }
860 
861  else if (i < n_dof_at_all_vertices) // range of i: 8..39
862  {
863  /* belongs to vertex in the outer shell
864  *
865  * subtract the number of dof already counted,
866  * so that i_offset contains only the offset for the base
867  */
868  const unsigned int i_offset = i - n_dof_at_base_vertices; // 0..31
869 
870  // first the radial dof are counted, then the base dof
871  radial_shape = (i_offset % radial_order) + 1;
872  base_shape = i_offset / radial_order;
873  }
874 
875  else if (i < n_dof_at_all_vertices+n_dof_at_base_sides) // range of i: 40..51
876  {
877  // belongs to base, is a side node
878  radial_shape = 0;
879  base_shape = i - radial_order * n_dof_at_base_vertices; // 8..19
880  }
881 
882  else if (i < n_dof_at_all_vertices+n_dof_at_all_sides) // range of i: 52..99
883  {
884  // belongs to side node in the outer shell
885  const unsigned int i_offset = i - (n_dof_at_all_vertices
886  + n_dof_at_base_sides); // 0..47
887  radial_shape = (i_offset % radial_order) + 1;
888  base_shape = (i_offset / radial_order) + n_dof_at_base_vertices;
889  }
890 
891  else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_base_face) // range of i: 100..104
892  {
893  // belongs to the node in the base face
894  radial_shape = 0;
895  base_shape = i - radial_order*(n_dof_at_base_vertices
896  + n_dof_at_base_sides); // 20..24
897  }
898 
899  else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_all_faces) // range of i: 105..124
900  {
901  // belongs to the node in the outer face
902  const unsigned int i_offset = i - (n_dof_at_all_vertices
903  + n_dof_at_all_sides
904  + n_dof_at_base_face); // 0..19
905  radial_shape = (i_offset % radial_order) + 1;
906  base_shape = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides;
907  }
908 
909  else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_all_faces+n_base_elem_dof) // range of i: 125..133
910  {
911  // belongs to the base and is an element associated shape
912  radial_shape = 0;
913  base_shape = i - (n_dof_at_all_vertices
914  + n_dof_at_all_sides
915  + n_dof_at_all_faces); // 0..8
916  }
917 
918  else // range of i: 134..169
919  {
920  libmesh_assert_less (i, n_dofs(fet, inf_elem_type));
921  // belongs to the outer shell and is an element associated shape
922  const unsigned int i_offset = i - (n_dof_at_all_vertices
923  + n_dof_at_all_sides
924  + n_dof_at_all_faces
925  + n_base_elem_dof); // 0..19
926  radial_shape = (i_offset % radial_order) + 1;
927  base_shape = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides + n_dof_at_base_face;
928  }
929  }
930 
931  return;
932 }
933 
934 
935 
936 //--------------------------------------------------------------
937 // Explicit instantiations
938 //#include "libmesh/inf_fe_instantiate_1D.h"
939 //#include "libmesh/inf_fe_instantiate_2D.h"
940 //#include "libmesh/inf_fe_instantiate_3D.h"
941 
942 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs(const FEType &, const ElemType));
943 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs(const FEType &, const ElemType));
944 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs(const FEType &, const ElemType));
945 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const ElemType));
946 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const ElemType));
947 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const ElemType));
948 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const ElemType, const unsigned int));
949 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const ElemType, const unsigned int));
950 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const ElemType, const unsigned int));
951 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_shape_indices(const FEType &, const ElemType, const unsigned int, unsigned int &, unsigned int &));
952 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_shape_indices(const FEType &, const ElemType, const unsigned int, unsigned int &, unsigned int &));
953 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_shape_indices(const FEType &, const ElemType, const unsigned int, unsigned int &, unsigned int &));
954 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_node_indices(const ElemType, const unsigned int, unsigned int &, unsigned int &));
955 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_node_indices(const ElemType, const unsigned int, unsigned int &, unsigned int &));
956 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_node_indices(const ElemType, const unsigned int, unsigned int &, unsigned int &));
957 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape(const FEType &, const Elem *, const unsigned int, const Point & p));
958 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape(const FEType &, const Elem *, const unsigned int, const Point & p));
959 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape(const FEType &, const Elem *, const unsigned int, const Point & p));
960 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
961 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
962 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
963 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_data(const FEType &, const Elem *, FEComputeData &));
964 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_data(const FEType &, const Elem *, FEComputeData &));
965 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_data(const FEType &, const Elem *, FEComputeData &));
966 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, nodal_soln(const FEType &, const Elem *, const std::vector<Number> &, std::vector<Number> &));
967 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, nodal_soln(const FEType &, const Elem *, const std::vector<Number> &, std::vector<Number> &));
968 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, nodal_soln(const FEType &, const Elem *, const std::vector<Number> &, std::vector<Number> &));
969 
970 } // namespace libMesh
971 
972 #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
static void compute_node_indices_fast(const ElemType inf_elem_type, const unsigned int outer_node_index, unsigned int &base_node, unsigned int &radial_node)
Does the same as compute_node_indices(), but stores the maps for the current element type...
static unsigned int n_dofs_per_elem(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:473
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
OStreamProxy err
int get_order() const
Explicitly request the order as an int.
Definition: fe_type.h:77
static unsigned int n_dofs_per_elem(const FEType &fet, const ElemType inf_elem_type)
virtual UniquePtr< Elem > build_side_ptr(const unsigned int i, bool proxy=true)=0
static bool _warned_for_shape
Definition: inf_fe.h:807
A specific instantiation of the FEBase class.
Definition: fe.h:40
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:414
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value.
Definition: libmesh.h:184
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
virtual ElemType type() const =0
static unsigned int n_dofs_at_node(const FEType &fet, const ElemType inf_elem_type, const unsigned int n)
Definition: inf_fe_static.C:72
class FEComputeData hides arbitrary data to be passed to and from children of FEBase through the FEIn...
const Point & p
Holds the point where the data are to be computed.
ElemType
Defines an enum for geometric element types.
OrderWrapper radial_order
The approximation order in the base of the infinite element.
Definition: fe_type.h:236
Real phase
Storage for the computed phase lag.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Elem *, Base::build_elem(const Elem *))
static Real shape(const FEType &fet, const ElemType t, const unsigned int i, const Point &p)
The libMesh namespace provides an interface to certain functionality in the library.
static void compute_node_indices(const ElemType inf_elem_type, const unsigned int outer_node_index, unsigned int &base_node, unsigned int &radial_node)
Computes the indices in the base base_node and in radial direction radial_node (either 0 or 1) associ...
Real speed
The wave speed.
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
const dof_id_type n_nodes
Definition: tecplot_io.C:67
static unsigned int n_dofs(const FEType &fet, const ElemType inf_elem_type)
Definition: inf_fe_static.C:55
static Real decay(const Real v)
Definition: inf_fe.h:826
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:641
static Real eval(Real v, Order o_radial, unsigned int i)
Number frequency
The frequency to evaluate shape functions including the wave number depending terms.
std::vector< Number > shape
Storage for the computed shape function values.
virtual Point origin() const
Definition: elem.h:1480
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:436
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
static void compute_shape_indices(const FEType &fet, const ElemType inf_elem_type, const unsigned int i, unsigned int &base_shape, unsigned int &radial_shape)
Computes the indices of shape functions in the base base_shape and in radial direction radial_shape (...
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
static bool _warned_for_nodal_soln
static members that are used to issue warning messages only once.
Definition: inf_fe.h:806
static void compute_data(const FEType &fe_t, const Elem *inf_elem, FEComputeData &data)
Generalized version of shape(), takes an Elem *.
static void nodal_soln(const FEType &fet, const Elem *elem, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
Usually, this method would build the nodal soln from the element soln.
IterBase * data
Ideally this private member data should have protected access.
static ElemType _compute_node_indices_fast_current_elem_type
When compute_node_indices_fast() is used, this static variable remembers the element type for which t...
Definition: inf_fe.h:798
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
const Real pi
.
Definition: libmesh.h:172