inf_fe_static.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 // 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"
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>
36 ElemType InfFE<Dim,T_radial,T_map>::_compute_node_indices_fast_current_elem_type = INVALID_ELEM;
37 
38 #ifdef DEBUG
39 
40 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
41 bool InfFE<Dim,T_radial,T_map>::_warned_for_nodal_soln = false;
42 
43 
44 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
45 bool InfFE<Dim,T_radial,T_map>::_warned_for_shape = false;
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  if (Dim > 1)
182  return FEInterface::shape(Dim-1, fet, base_et, i_base, p)
183  * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
185  else
186  return InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
188 }
189 
190 
191 
192 
193 
194 
195 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
197  const Elem* inf_elem,
198  const unsigned int i,
199  const Point& p)
200 {
201  libmesh_assert(inf_elem);
202  libmesh_assert_not_equal_to (Dim, 0);
203 
204 #ifdef DEBUG
205  // this makes only sense when used for mapping
206  if ((T_radial != INFINITE_MAP) && !_warned_for_shape)
207  {
208  libMesh::err << "WARNING: InfFE<Dim,T_radial,T_map>::shape(...) does _not_" << std::endl
209  << " return the correct trial function! Use " << std::endl
210  << " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!"
211  << std::endl;
212  _warned_for_shape = true;
213  }
214 #endif
215 
216  const Order o_radial (fet.radial_order);
217  const Real v (p(Dim-1));
218  AutoPtr<Elem> base_el (inf_elem->build_side(0));
219 
220  unsigned int i_base, i_radial;
221  compute_shape_indices(fet, inf_elem->type(), i, i_base, i_radial);
222 
223  if (Dim > 1)
224  return FEInterface::shape(Dim-1, fet, base_el.get(), i_base, p)
225  * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
227  else
228  return InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
230 }
231 
232 
233 
234 
235 
236 
237 
238 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
240  const Elem* inf_elem,
242 {
243  libmesh_assert(inf_elem);
244  libmesh_assert_not_equal_to (Dim, 0);
245 
246  const Order o_radial (fet.radial_order);
247  const Order radial_mapping_order (Radial::mapping_order());
248  const Point& p (data.p);
249  const Real v (p(Dim-1));
250  AutoPtr<Elem> base_el (inf_elem->build_side(0));
251 
252  /*
253  * compute \p interpolated_dist containing the mapping-interpolated
254  * distance of the base point to the origin. This is the same
255  * for all shape functions. Set \p interpolated_dist to 0, it
256  * is added to.
257  */
258  Real interpolated_dist = 0.;
259  switch (Dim)
260  {
261  case 1:
262  {
263  libmesh_assert_equal_to (inf_elem->type(), INFEDGE2);
264  interpolated_dist = Point(inf_elem->point(0) - inf_elem->point(1)).size();
265  break;
266  }
267 
268  case 2:
269  {
270  const unsigned int n_base_nodes = base_el->n_nodes();
271 
272  const Point origin = inf_elem->origin();
273  const Order base_mapping_order (base_el->default_order());
274  const ElemType base_mapping_elem_type (base_el->type());
275 
276  // interpolate the base nodes' distances
277  for (unsigned int n=0; n<n_base_nodes; n++)
278  interpolated_dist += Point(base_el->point(n) - origin).size()
279  * FE<1,LAGRANGE>::shape (base_mapping_elem_type, base_mapping_order, n, p);
280  break;
281  }
282 
283  case 3:
284  {
285  const unsigned int n_base_nodes = base_el->n_nodes();
286 
287  const Point origin = inf_elem->origin();
288  const Order base_mapping_order (base_el->default_order());
289  const ElemType base_mapping_elem_type (base_el->type());
290 
291  // interpolate the base nodes' distances
292  for (unsigned int n=0; n<n_base_nodes; n++)
293  interpolated_dist += Point(base_el->point(n) - origin).size()
294  * FE<2,LAGRANGE>::shape (base_mapping_elem_type, base_mapping_order, n, p);
295  break;
296  }
297 #ifdef DEBUG
298  default:
299  libmesh_error();
300 #endif
301  }
302 
303 
304 
305 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
306 
307  // assumption on time-harmonic behavior
308  const short int sign (-1);
309 
310  // the wave number
311  const Real wavenumber = 2. * libMesh::pi * data.frequency / data.speed;
312 
313  // the exponent for time-harmonic behavior
314  const Real exponent = sign /* +1. or -1. */
315  * wavenumber /* k */
316  * interpolated_dist /* together with next line: */
317  * InfFE<Dim,INFINITE_MAP,T_map>::eval(v, radial_mapping_order, 1); /* phase(s,t,v) */
318 
319  const Number time_harmonic = Number(cos(exponent), sin(exponent)); /* e^(sign*i*k*phase(s,t,v)) */
320 
321  /*
322  * compute \p shape for all dof in the element
323  */
324  if (Dim > 1)
325  {
326  const unsigned int n_dof = n_dofs (fet, inf_elem->type());
327  data.shape.resize(n_dof);
328 
329  for (unsigned int i=0; i<n_dof; i++)
330  {
331  // compute base and radial shape indices
332  unsigned int i_base, i_radial;
333  compute_shape_indices(fet, inf_elem->type(), i, i_base, i_radial);
334 
335  data.shape[i] = (InfFE<Dim,T_radial,T_map>::Radial::decay(v) /* (1.-v)/2. in 3D */
336  * FEInterface::shape(Dim-1, fet, base_el.get(), i_base, p) /* S_n(s,t) */
337  * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)) /* L_n(v) */
338  * time_harmonic; /* e^(sign*i*k*phase(s,t,v) */
339  }
340  }
341  else
342  {
343  libMesh::err << "compute_data() for 1-dimensional InfFE not implemented." << std::endl;
344  libmesh_error();
345  }
346 
347 #else
348 
349  const Real speed = data.speed;
350 
351  /*
352  * This is quite weird: the phase is actually
353  * a measure how @e advanced the pressure is that
354  * we compute. In other words: the further away
355  * the node \p data.p is, the further we look into
356  * the future...
357  */
358  data.phase = interpolated_dist /* phase(s,t,v)/c */
359  * InfFE<Dim,INFINITE_MAP,T_map>::eval(v, radial_mapping_order, 1) / speed;
360 
361  if (Dim > 1)
362  {
363  const unsigned int n_dof = n_dofs (fet, inf_elem->type());
364  data.shape.resize(n_dof);
365 
366  for (unsigned int i=0; i<n_dof; i++)
367  {
368  // compute base and radial shape indices
369  unsigned int i_base, i_radial;
370  compute_shape_indices(fet, inf_elem->type(), i, i_base, i_radial);
371 
372  data.shape[i] = InfFE<Dim,T_radial,T_map>::Radial::decay(v) /* (1.-v)/2. in 3D */
373  * FEInterface::shape(Dim-1, fet, base_el.get(), i_base, p) /* S_n(s,t) */
374  * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial); /* L_n(v) */
375  }
376  }
377  else
378  {
379  libMesh::err << "compute_data() for 1-dimensional InfFE not implemented." << std::endl;
380  libmesh_error();
381  }
382 
383 #endif
384 }
385 
386 
387 
388 
389 
390 
391 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
393  const unsigned int outer_node_index,
394  unsigned int& base_node,
395  unsigned int& radial_node)
396 {
397  switch (inf_elem_type)
398  {
399  case INFEDGE2:
400  {
401  libmesh_assert_less (outer_node_index, 2);
402  base_node = 0;
403  radial_node = outer_node_index;
404  return;
405  }
406 
407 
408  // linear base approximation, easy to determine
409  case INFQUAD4:
410  {
411  libmesh_assert_less (outer_node_index, 4);
412  base_node = outer_node_index % 2;
413  radial_node = outer_node_index / 2;
414  return;
415  }
416 
417  case INFPRISM6:
418  {
419  libmesh_assert_less (outer_node_index, 6);
420  base_node = outer_node_index % 3;
421  radial_node = outer_node_index / 3;
422  return;
423  }
424 
425  case INFHEX8:
426  {
427  libmesh_assert_less (outer_node_index, 8);
428  base_node = outer_node_index % 4;
429  radial_node = outer_node_index / 4;
430  return;
431  }
432 
433 
434  // higher order base approximation, more work necessary
435  case INFQUAD6:
436  {
437  switch (outer_node_index)
438  {
439  case 0:
440  case 1:
441  {
442  radial_node = 0;
443  base_node = outer_node_index;
444  return;
445  }
446 
447  case 2:
448  case 3:
449  {
450  radial_node = 1;
451  base_node = outer_node_index-2;
452  return;
453  }
454 
455  case 4:
456  {
457  radial_node = 0;
458  base_node = 2;
459  return;
460  }
461 
462  case 5:
463  {
464  radial_node = 1;
465  base_node = 2;
466  return;
467  }
468 
469  default:
470  {
471  libmesh_error();
472  return;
473  }
474  }
475  }
476 
477 
478  case INFHEX16:
479  case INFHEX18:
480  {
481  switch (outer_node_index)
482  {
483  case 0:
484  case 1:
485  case 2:
486  case 3:
487  {
488  radial_node = 0;
489  base_node = outer_node_index;
490  return;
491  }
492 
493  case 4:
494  case 5:
495  case 6:
496  case 7:
497  {
498  radial_node = 1;
499  base_node = outer_node_index-4;
500  return;
501  }
502 
503  case 8:
504  case 9:
505  case 10:
506  case 11:
507  {
508  radial_node = 0;
509  base_node = outer_node_index-4;
510  return;
511  }
512 
513  case 12:
514  case 13:
515  case 14:
516  case 15:
517  {
518  radial_node = 1;
519  base_node = outer_node_index-8;
520  return;
521  }
522 
523  case 16:
524  {
525  libmesh_assert_equal_to (inf_elem_type, INFHEX18);
526  radial_node = 0;
527  base_node = 8;
528  return;
529  }
530 
531  case 17:
532  {
533  libmesh_assert_equal_to (inf_elem_type, INFHEX18);
534  radial_node = 1;
535  base_node = 8;
536  return;
537  }
538 
539  default:
540  {
541  libmesh_error();
542  return;
543  }
544  }
545  }
546 
547 
548  case INFPRISM12:
549  {
550  switch (outer_node_index)
551  {
552  case 0:
553  case 1:
554  case 2:
555  {
556  radial_node = 0;
557  base_node = outer_node_index;
558  return;
559  }
560 
561  case 3:
562  case 4:
563  case 5:
564  {
565  radial_node = 1;
566  base_node = outer_node_index-3;
567  return;
568  }
569 
570  case 6:
571  case 7:
572  case 8:
573  {
574  radial_node = 0;
575  base_node = outer_node_index-3;
576  return;
577  }
578 
579  case 9:
580  case 10:
581  case 11:
582  {
583  radial_node = 1;
584  base_node = outer_node_index-6;
585  return;
586  }
587 
588  default:
589  {
590  libmesh_error();
591  return;
592  }
593  }
594  }
595 
596 
597  default:
598  {
599  libMesh::err << "ERROR: Bad infinite element type=" << inf_elem_type
600  << ", node=" << outer_node_index << std::endl;
601  libmesh_error();
602  return;
603  }
604  }
605 }
606 
607 
608 
609 
610 
611 
612 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
614  const unsigned int outer_node_index,
615  unsigned int& base_node,
616  unsigned int& radial_node)
617 {
618  libmesh_assert_not_equal_to (inf_elem_type, INVALID_ELEM);
619 
620  static std::vector<unsigned int> _static_base_node_index;
621  static std::vector<unsigned int> _static_radial_node_index;
622 
623  /*
624  * fast counterpart to compute_node_indices(), uses local static buffers
625  * to store the index maps. The class member
626  * \p _compute_node_indices_fast_current_elem_type remembers
627  * the current element type.
628  *
629  * Note that there exist non-static members storing the
630  * same data. However, you never know what element type
631  * is currently used by the \p InfFE object, and what
632  * request is currently directed to the static \p InfFE
633  * members (which use \p compute_node_indices_fast()).
634  * So separate these.
635  *
636  * check whether the work for this elemtype has already
637  * been done. If so, use this index. Otherwise, refresh
638  * the buffer to this element type.
639  */
640  if (inf_elem_type==_compute_node_indices_fast_current_elem_type)
641  {
642  base_node = _static_base_node_index [outer_node_index];
643  radial_node = _static_radial_node_index[outer_node_index];
644  return;
645  }
646  else
647  {
648  // store the map for _all_ nodes for this element type
649  _compute_node_indices_fast_current_elem_type = inf_elem_type;
650 
651  unsigned int n_nodes = libMesh::invalid_uint;
652 
653  switch (inf_elem_type)
654  {
655  case INFEDGE2:
656  {
657  n_nodes = 2;
658  break;
659  }
660  case INFQUAD4:
661  {
662  n_nodes = 4;
663  break;
664  }
665  case INFQUAD6:
666  {
667  n_nodes = 6;
668  break;
669  }
670  case INFHEX8:
671  {
672  n_nodes = 8;
673  break;
674  }
675  case INFHEX16:
676  {
677  n_nodes = 16;
678  break;
679  }
680  case INFHEX18:
681  {
682  n_nodes = 18;
683  break;
684  }
685  case INFPRISM6:
686  {
687  n_nodes = 6;
688  break;
689  }
690  case INFPRISM12:
691  {
692  n_nodes = 12;
693  break;
694  }
695  default:
696  {
697  libMesh::err << "ERROR: Bad infinite element type=" << inf_elem_type
698  << ", node=" << outer_node_index << std::endl;
699  libmesh_error();
700  break;
701  }
702  }
703 
704 
705  _static_base_node_index.resize (n_nodes);
706  _static_radial_node_index.resize(n_nodes);
707 
708  for (unsigned int n=0; n<n_nodes; n++)
709  compute_node_indices (inf_elem_type,
710  n,
711  _static_base_node_index [outer_node_index],
712  _static_radial_node_index[outer_node_index]);
713 
714  // and return for the specified node
715  base_node = _static_base_node_index [outer_node_index];
716  radial_node = _static_radial_node_index[outer_node_index];
717  return;
718  }
719 }
720 
721 
722 
723 
724 
725 
726 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
728  const ElemType inf_elem_type,
729  const unsigned int i,
730  unsigned int& base_shape,
731  unsigned int& radial_shape)
732 {
733 
734  /*
735  * An example is provided: the numbers in comments refer to
736  * a fictitious InfHex18. The numbers are chosen as exemplary
737  * values. There is currently no base approximation that
738  * requires this many dof's at nodes, sides, faces and in the element.
739  *
740  * the order of the shape functions is heavily related with the
741  * order the dofs are assigned in \p DofMap::distributed_dofs().
742  * Due to the infinite elements with higher-order base approximation,
743  * some more effort is necessary.
744  *
745  * numbering scheme:
746  * 1. all vertices in the base, assign node->n_comp() dofs to each vertex
747  * 2. all vertices further out: innermost loop: radial shapes,
748  * then the base approximation shapes
749  * 3. all side nodes in the base, assign node->n_comp() dofs to each side node
750  * 4. all side nodes further out: innermost loop: radial shapes,
751  * then the base approximation shapes
752  * 5. (all) face nodes in the base, assign node->n_comp() dofs to each face node
753  * 6. (all) face nodes further out: innermost loop: radial shapes,
754  * then the base approximation shapes
755  * 7. element-associated dof in the base
756  * 8. element-associated dof further out
757  */
758 
759  const unsigned int radial_order = static_cast<unsigned int>(fet.radial_order); // 4
760  const unsigned int radial_order_p_one = radial_order+1; // 5
761 
762  const ElemType base_elem_type (Base::get_elem_type(inf_elem_type)); // QUAD9
763 
764  // assume that the number of dof is the same for all vertices
765  unsigned int n_base_vertices = libMesh::invalid_uint; // 4
766  const unsigned int n_base_vertex_dof = FEInterface::n_dofs_at_node (Dim-1, fet, base_elem_type, 0);// 2
767 
768  unsigned int n_base_side_nodes = libMesh::invalid_uint; // 4
769  unsigned int n_base_side_dof = libMesh::invalid_uint; // 3
770 
771  unsigned int n_base_face_nodes = libMesh::invalid_uint; // 1
772  unsigned int n_base_face_dof = libMesh::invalid_uint; // 5
773 
774  const unsigned int n_base_elem_dof = FEInterface::n_dofs_per_elem (Dim-1, fet, base_elem_type);// 9
775 
776 
777  switch (inf_elem_type)
778  {
779  case INFEDGE2:
780  {
781  n_base_vertices = 1;
782  n_base_side_nodes = 0;
783  n_base_face_nodes = 0;
784  n_base_side_dof = 0;
785  n_base_face_dof = 0;
786  break;
787  }
788 
789  case INFQUAD4:
790  {
791  n_base_vertices = 2;
792  n_base_side_nodes = 0;
793  n_base_face_nodes = 0;
794  n_base_side_dof = 0;
795  n_base_face_dof = 0;
796  break;
797  }
798 
799  case INFQUAD6:
800  {
801  n_base_vertices = 2;
802  n_base_side_nodes = 1;
803  n_base_face_nodes = 0;
804  n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
805  n_base_face_dof = 0;
806  break;
807  }
808 
809  case INFHEX8:
810  {
811  n_base_vertices = 4;
812  n_base_side_nodes = 0;
813  n_base_face_nodes = 0;
814  n_base_side_dof = 0;
815  n_base_face_dof = 0;
816  break;
817  }
818 
819  case INFHEX16:
820  {
821  n_base_vertices = 4;
822  n_base_side_nodes = 4;
823  n_base_face_nodes = 0;
824  n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
825  n_base_face_dof = 0;
826  break;
827  }
828 
829  case INFHEX18:
830  {
831  n_base_vertices = 4;
832  n_base_side_nodes = 4;
833  n_base_face_nodes = 1;
834  n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
835  n_base_face_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, 8);
836  break;
837  }
838 
839 
840  case INFPRISM6:
841  {
842  n_base_vertices = 3;
843  n_base_side_nodes = 0;
844  n_base_face_nodes = 0;
845  n_base_side_dof = 0;
846  n_base_face_dof = 0;
847  break;
848  }
849 
850  case INFPRISM12:
851  {
852  n_base_vertices = 3;
853  n_base_side_nodes = 3;
854  n_base_face_nodes = 0;
855  n_base_side_dof = FEInterface::n_dofs_at_node (Dim-1, fet,base_elem_type, n_base_vertices);
856  n_base_face_dof = 0;
857  break;
858  }
859 
860  default:
861  libmesh_error();
862  }
863 
864 
865  {
866  // these are the limits describing the intervals where the shape function lies
867  const unsigned int n_dof_at_base_vertices = n_base_vertices*n_base_vertex_dof; // 8
868  const unsigned int n_dof_at_all_vertices = n_dof_at_base_vertices*radial_order_p_one; // 40
869 
870  const unsigned int n_dof_at_base_sides = n_base_side_nodes*n_base_side_dof; // 12
871  const unsigned int n_dof_at_all_sides = n_dof_at_base_sides*radial_order_p_one; // 60
872 
873  const unsigned int n_dof_at_base_face = n_base_face_nodes*n_base_face_dof; // 5
874  const unsigned int n_dof_at_all_faces = n_dof_at_base_face*radial_order_p_one; // 25
875 
876 
877  // start locating the shape function
878  if (i < n_dof_at_base_vertices) // range of i: 0..7
879  {
880  // belongs to vertex in the base
881  radial_shape = 0;
882  base_shape = i;
883  }
884 
885  else if (i < n_dof_at_all_vertices) // range of i: 8..39
886  {
887  /* belongs to vertex in the outer shell
888  *
889  * subtract the number of dof already counted,
890  * so that i_offset contains only the offset for the base
891  */
892  const unsigned int i_offset = i - n_dof_at_base_vertices; // 0..31
893 
894  // first the radial dof are counted, then the base dof
895  radial_shape = (i_offset % radial_order) + 1;
896  base_shape = i_offset / radial_order;
897  }
898 
899  else if (i < n_dof_at_all_vertices+n_dof_at_base_sides) // range of i: 40..51
900  {
901  // belongs to base, is a side node
902  radial_shape = 0;
903  base_shape = i - radial_order * n_dof_at_base_vertices; // 8..19
904  }
905 
906  else if (i < n_dof_at_all_vertices+n_dof_at_all_sides) // range of i: 52..99
907  {
908  // belongs to side node in the outer shell
909  const unsigned int i_offset = i - (n_dof_at_all_vertices
910  + n_dof_at_base_sides); // 0..47
911  radial_shape = (i_offset % radial_order) + 1;
912  base_shape = (i_offset / radial_order) + n_dof_at_base_vertices;
913  }
914 
915  else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_base_face) // range of i: 100..104
916  {
917  // belongs to the node in the base face
918  radial_shape = 0;
919  base_shape = i - radial_order*(n_dof_at_base_vertices
920  + n_dof_at_base_sides); // 20..24
921  }
922 
923  else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_all_faces) // range of i: 105..124
924  {
925  // belongs to the node in the outer face
926  const unsigned int i_offset = i - (n_dof_at_all_vertices
927  + n_dof_at_all_sides
928  + n_dof_at_base_face); // 0..19
929  radial_shape = (i_offset % radial_order) + 1;
930  base_shape = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides;
931  }
932 
933  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
934  {
935  // belongs to the base and is an element associated shape
936  radial_shape = 0;
937  base_shape = i - (n_dof_at_all_vertices
938  + n_dof_at_all_sides
939  + n_dof_at_all_faces); // 0..8
940  }
941 
942  else // range of i: 134..169
943  {
944  libmesh_assert_less (i, n_dofs(fet, inf_elem_type));
945  // belongs to the outer shell and is an element associated shape
946  const unsigned int i_offset = i - (n_dof_at_all_vertices
947  + n_dof_at_all_sides
948  + n_dof_at_all_faces
949  + n_base_elem_dof); // 0..19
950  radial_shape = (i_offset % radial_order) + 1;
951  base_shape = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides + n_dof_at_base_face;
952  }
953  }
954 
955  return;
956 }
957 
958 
959 
960 //--------------------------------------------------------------
961 // Explicit instantiations
962 //#include "libmesh/inf_fe_instantiate_1D.h"
963 //#include "libmesh/inf_fe_instantiate_2D.h"
964 //#include "libmesh/inf_fe_instantiate_3D.h"
965 
966 INSTANTIATE_INF_FE_MBRF(1,CARTESIAN,unsigned int,n_dofs(const FEType&,const ElemType));
967 INSTANTIATE_INF_FE_MBRF(2,CARTESIAN,unsigned int,n_dofs(const FEType&,const ElemType));
968 INSTANTIATE_INF_FE_MBRF(3,CARTESIAN,unsigned int,n_dofs(const FEType&,const ElemType));
969 INSTANTIATE_INF_FE_MBRF(1,CARTESIAN,unsigned int,n_dofs_per_elem(const FEType&,const ElemType));
970 INSTANTIATE_INF_FE_MBRF(2,CARTESIAN,unsigned int,n_dofs_per_elem(const FEType&,const ElemType));
971 INSTANTIATE_INF_FE_MBRF(3,CARTESIAN,unsigned int,n_dofs_per_elem(const FEType&,const ElemType));
972 INSTANTIATE_INF_FE_MBRF(1,CARTESIAN,unsigned int,n_dofs_at_node(const FEType&,const ElemType,const unsigned int));
973 INSTANTIATE_INF_FE_MBRF(2,CARTESIAN,unsigned int,n_dofs_at_node(const FEType&,const ElemType,const unsigned int));
974 INSTANTIATE_INF_FE_MBRF(3,CARTESIAN,unsigned int,n_dofs_at_node(const FEType&,const ElemType,const unsigned int));
975 INSTANTIATE_INF_FE_MBRF(1,CARTESIAN,void,compute_shape_indices(const FEType&,const ElemType,const unsigned int,unsigned int&,unsigned int&));
976 INSTANTIATE_INF_FE_MBRF(2,CARTESIAN,void,compute_shape_indices(const FEType&,const ElemType,const unsigned int,unsigned int&,unsigned int&));
977 INSTANTIATE_INF_FE_MBRF(3,CARTESIAN,void,compute_shape_indices(const FEType&,const ElemType,const unsigned int,unsigned int&,unsigned int&));
978 INSTANTIATE_INF_FE_MBRF(1,CARTESIAN,void,compute_node_indices(const ElemType,const unsigned int,unsigned int&,unsigned int&));
979 INSTANTIATE_INF_FE_MBRF(2,CARTESIAN,void,compute_node_indices(const ElemType,const unsigned int,unsigned int&,unsigned int&));
980 INSTANTIATE_INF_FE_MBRF(3,CARTESIAN,void,compute_node_indices(const ElemType,const unsigned int,unsigned int&,unsigned int&));
981 INSTANTIATE_INF_FE_MBRF(1,CARTESIAN,Real,shape(const FEType&,const Elem*,const unsigned int,const Point& p));
982 INSTANTIATE_INF_FE_MBRF(2,CARTESIAN,Real,shape(const FEType&,const Elem*,const unsigned int,const Point& p));
983 INSTANTIATE_INF_FE_MBRF(3,CARTESIAN,Real,shape(const FEType&,const Elem*,const unsigned int,const Point& p));
984 INSTANTIATE_INF_FE_MBRF(1,CARTESIAN,Real,shape(const FEType&,const ElemType,const unsigned int,const Point&));
985 INSTANTIATE_INF_FE_MBRF(2,CARTESIAN,Real,shape(const FEType&,const ElemType,const unsigned int,const Point&));
986 INSTANTIATE_INF_FE_MBRF(3,CARTESIAN,Real,shape(const FEType&,const ElemType,const unsigned int,const Point&));
987 INSTANTIATE_INF_FE_MBRF(1,CARTESIAN,void,compute_data(const FEType&,const Elem*,FEComputeData&));
988 INSTANTIATE_INF_FE_MBRF(2,CARTESIAN,void,compute_data(const FEType&,const Elem*,FEComputeData&));
989 INSTANTIATE_INF_FE_MBRF(3,CARTESIAN,void,compute_data(const FEType&,const Elem*,FEComputeData&));
990 INSTANTIATE_INF_FE_MBRF(1,CARTESIAN,void,nodal_soln(const FEType&,const Elem*,const std::vector<Number>&,std::vector<Number>&));
991 INSTANTIATE_INF_FE_MBRF(2,CARTESIAN,void,nodal_soln(const FEType&,const Elem*,const std::vector<Number>&,std::vector<Number>&));
992 INSTANTIATE_INF_FE_MBRF(3,CARTESIAN,void,nodal_soln(const FEType&,const Elem*,const std::vector<Number>&,std::vector<Number>&));
993 
994 } // namespace libMesh
995 
996 #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS

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

Hosted By:
SourceForge.net Logo