1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 #ifndef LIBMESH_INF_FE_H
21 #define LIBMESH_INF_FE_H
22
23 #include "libmesh/libmesh_config.h"
24
25 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
26
27 // Local includes
28 #include "libmesh/fe_base.h"
29 #include "libmesh/inf_fe_map.h"
30
31 // C++ includes
32 #include <cstddef>
33
34 namespace libMesh
35 {
36
37
38 // forward declarations
39 class Elem;
40 class FEComputeData;
41
42
43 /**
44 * Infinite elements are in some sense directional, compared
45 * to conventional finite elements. All methods related
46 * to the radial part, which extends perpendicular from the base,
47 * are collected in this class. This class offers static methods for
48 * use in \p InfFE and \p InfFEMap
49 *
50 * \author Daniel Dreyer
51 * \date 2003
52 */
53 class InfFERadial
54 {
55 private:
56
57 /**
58 * Never use an object of this type.
59 */
InfFERadial()60 InfFERadial() {}
61
62 public:
63
64 /**
65 * \returns The decay in the radial direction of
66 * the \p Dim dimensional infinite element.
67 */
68 static Real decay (const unsigned int dim, const Real v);
69
70 /**
71 * \returns The first (local) derivative of the
72 * decay in radial direction of the infinite element.
73 *
74 * This is only valid for 3D?? - RHS
75 */
decay_deriv(const Real)76 static Real decay_deriv (const Real) { return -.5; }
77
78 /**
79 * \returns The radial weight D, used as an additional weight
80 * for the test function, evaluated at local radial coordinate \p v.
81 */
D(const Real v)82 static Real D (const Real v) { return (1.-v)*(1.-v)/4.; }
83
84 /**
85 * \returns The first (local) radial derivative of the radial weight D.
86 */
D_deriv(const Real v)87 static Real D_deriv (const Real v) { return (v-1.)/2.; }
88
89 /**
90 * \returns The Order of the mapping functions
91 * in the radial direction. Currently, this is always \p FIRST.
92 */
mapping_order()93 static Order mapping_order() { return FIRST; }
94
95 /**
96 * \returns The number of shape functions in radial direction
97 * associated with this infinite element.
98 * Either way, if the modes are stored as nodal dofs (\p n_dofs_at_node)
99 * or as element dofs (\p n_dofs_per_elem), in each case we have the
100 * same number of modes in radial direction.
101 *
102 * \note For the case of 1D infinite elements, in the base the
103 * dof-per-node scheme is used.
104 *
105 * From the formulation of the infinite elements, we have
106 * 1 mode, when \p o_radial=CONST.
107 * Therefore, we have a total of \p o_radial+1 modes in radial direction.
108 */
n_dofs(const Order o_radial)109 static unsigned int n_dofs (const Order o_radial)
110 { return static_cast<unsigned int>(o_radial)+1; }
111
112 /**
113 * \returns The number of dofs in radial direction on "onion slice"
114 * \p n (either 0 or 1) for an infinite element of type \p inf_elem_type and
115 * radial order \p o_radial.
116 *
117 * Currently, the first radial mode is associated with the nodes in
118 * the base. All higher radial modes are associated with
119 * the physically existing nodes further out.
120 */
121 static unsigned int n_dofs_at_node (const Order o_radial,
122 const unsigned int n_onion);
123
124 /**
125 * \returns The number of modes in radial direction interior to the element,
126 * not associated with any interior nodes.
127 *
128 * \note These modes are a discontinuous approximation, therefore
129 * we have no special formulation for coupling in the base, like in the
130 * case of associating (possibly) multiple dofs per (outer) node.
131 */
n_dofs_per_elem(const Order o_radial)132 static unsigned int n_dofs_per_elem (const Order o_radial)
133 { return static_cast<unsigned int>(o_radial)+1; }
134
135 };
136
137
138
139 /**
140 * This nested class contains most of the static methods related
141 * to the base part of an infinite element. Only static members
142 * are provided, for use within \p InfFE and \p InfFEMap.
143 *
144 * \author Daniel Dreyer
145 * \date 2003
146 */
147 class InfFEBase
148 {
149 private:
150
151 /**
152 * Never use an object of this type.
153 */
InfFEBase()154 InfFEBase() {}
155
156 public:
157
158 /**
159 * Build the base element of an infinite element. Be careful,
160 * this method allocates memory! So be sure to delete the
161 * new element afterward.
162 */
163 static Elem * build_elem (const Elem * inf_elem);
164
165 /**
166 * \returns The base element associated to
167 * \p type. This is, for example, \p TRI3 for
168 * \p INFPRISM6.
169 */
170 static ElemType get_elem_type (const ElemType type);
171
172 /**
173 * \returns The number of shape functions used in the
174 * mapping in the base element of type \p base_elem_type
175 * mapped with order \p base_mapping_order
176 */
177 static unsigned int n_base_mapping_sf (const Elem & base_elem,
178 const Order base_mapping_order);
179
180 };
181
182
183
184 /**
185 * A specific instantiation of the \p FEBase class. This
186 * class is templated, and specific template instantiations
187 * will result in different Infinite Element families, similar
188 * to the \p FE class. \p InfFE builds a \p FE<Dim-1,T_base>,
189 * and most of the requests related to the base are handed over
190 * to this object. All methods related to the radial part
191 * are collected in the class \p InfFERadial. Similarly,
192 * most of the static methods concerning base approximation
193 * are contained in \p InfFEBase.
194 *
195 * Having different shape approximation families in radial direction
196 * introduces the requirement for an additional \p Order in this
197 * class. Therefore, the \p FEType internals change when infinite
198 * elements are enabled.
199 * When the specific infinite element type is not known at compile
200 * time, use the \p FEBase::build() member to create abstract
201 * (but still optimized) infinite elements at run time.
202 *
203 * The node numbering scheme is the one from the current
204 * infinite element. Each node in the base holds exactly
205 * the same number of dofs as an adjacent conventional \p FE
206 * would contain. The nodes further out hold the additional
207 * dof necessary for radial approximation. The order of the outer nodes'
208 * components is such that the radial shapes have highest
209 * priority, followed by the base shapes.
210 *
211 * \author Daniel Dreyer
212 * \date 2003
213 * \brief Base class for all the infinite geometric element types.
214 */
215 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
216 class InfFE : public FEBase
217 {
218
219 /*
220 * Protect the nested class
221 */
222 protected:
223
224
225 public:
226
227 // InfFE continued
228
229 /**
230 * Constructor and empty destructor.
231 * Initializes some data structures. Builds a \p FE<Dim-1,T_base>
232 * object to handle approximation in the base, so that
233 * there is no need to template \p InfFE<Dim,T_radial,T_map> also with
234 * respect to the base approximation \p T_base.
235 *
236 * The same remarks concerning compile-time optimization for
237 * \p FE also hold for \p InfFE. Use the
238 * \p FEBase::build_InfFE(const unsigned int, const FEType &)
239 * method to build specific instantiations of \p InfFE at
240 * run time.
241 */
242 explicit
243 InfFE(const FEType & fet);
~InfFE()244 ~InfFE() {}
245
246 // The static public members for access from FEInterface etc
247
248 /**
249 * \returns The value of the \f$ i^{th} \f$ shape function at
250 * point \p p. This method lets you specify the relevant
251 * data directly, and is therefore allowed to be static.
252 *
253 * \note This class member is not as efficient as its counterpart in
254 * \p FE<Dim,T>, and is not employed in the \p reinit() cycle.
255 *
256 * \note This method does not return physically correct shapes,
257 * instead use \p compute_data(). The \p shape() methods should
258 * only be used for mapping.
259 */
260 static Real shape(const FEType & fet,
261 const ElemType t,
262 const unsigned int i,
263 const Point & p);
264
265 /**
266 * \returns The value of the \f$ i^{th} \f$ shape function at
267 * point \p p. This method lets you specify the relevant
268 * data directly, and is therefore allowed to be static.
269 *
270 * \note This class member is not as efficient as its counterpart in
271 * \p FE<Dim,T>, and is not employed in the \p reinit() cycle.
272 *
273 * \note This method does not return physically correct shapes,
274 * instead use \p compute_data(). The \p shape() methods should
275 * only be used for mapping.
276 */
277 static Real shape(const FEType & fet,
278 const Elem * elem,
279 const unsigned int i,
280 const Point & p);
281
282 /**
283 * \returns The value of the \f$ i^{th} \f$ shape function at
284 * point \p p. This method lets you specify the relevant
285 * data directly, and is therefore allowed to be static.
286 *
287 * \note This class member is not as efficient as its counterpart in
288 * \p FE<Dim,T>, and is not employed in the \p reinit() cycle.
289 *
290 * \note This method does not return physically correct shapes,
291 * instead use \p compute_data(). The \p shape() methods should
292 * only be used for mapping.
293 */
294 static Real shape(const FEType fet,
295 const Elem * elem,
296 const unsigned int i,
297 const Point & p,
298 const bool add_p_level);
299
300 /**
301 * \returns The \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
302 * shape function at point \p p. This method lets you specify the relevant
303 * data directly, and is therefore allowed to be static.
304 *
305 * \note This class member is not as efficient as its counterpart in
306 * \p FE<Dim,T>, and is not employed in the \p reinit() cycle.
307 *
308 * \note This method does not return physically correct shape gradients,
309 * instead use \p compute_data(). The \p shape_deriv() methods should
310 * only be used for mapping.
311 */
312 static Real shape_deriv (const FEType & fet,
313 const Elem * inf_elem,
314 const unsigned int i,
315 const unsigned int j,
316 const Point & p);
317
318 /**
319 * \returns The \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
320 * shape function at point \p p. This method lets you specify the relevant
321 * data directly, and is therefore allowed to be static.
322 *
323 * \note This class member is not as efficient as its counterpart in
324 * \p FE<Dim,T>, and is not employed in the \p reinit() cycle.
325 *
326 * \note This method does not return physically correct shape gradients,
327 * instead use \p compute_data(). The \p shape_deriv() methods should
328 * only be used for mapping.
329 */
330 static Real shape_deriv (const FEType fet,
331 const Elem * inf_elem,
332 const unsigned int i,
333 const unsigned int j,
334 const Point & p,
335 const bool add_p_level);
336
337 /**
338 * \returns The \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
339 * shape function at point \p p. This method lets you specify the relevant
340 * data directly, and is therefore allowed to be static.
341 *
342 * \note This class member is not as efficient as its counterpart in
343 * \p FE<Dim,T>, and is not employed in the \p reinit() cycle.
344 *
345 * \note This method does not return physically correct shape gradients,
346 * instead use \p compute_data(). The \p shape_deriv() methods should
347 * only be used for mapping.
348 */
349 static Real shape_deriv (const FEType & fet,
350 const ElemType inf_elem_type,
351 const unsigned int i,
352 const unsigned int j,
353 const Point & p);
354
355 /**
356 * Generalized version of \p shape(), takes an \p Elem *. The \p data
357 * contains both input and output parameters. For frequency domain
358 * simulations, the complex-valued shape is returned. In time domain
359 * both the computed shape, and the phase is returned.
360 *
361 * \note The phase (proportional to the distance of the \p Point \p
362 * data.p from the envelope) is actually a measure how far into the
363 * future the results are.
364 */
365 static void compute_data(const FEType & fe_t,
366 const Elem * inf_elem,
367 FEComputeData & data);
368
369 /**
370 * \returns The number of shape functions associated with
371 * a finite element of type \p t and approximation order \p o.
372 *
373 * \deprecated Call the version of this function that takes a
374 * pointer-to-Elem instead.
375 */
n_shape_functions(const FEType & fet,const ElemType t)376 static unsigned int n_shape_functions (const FEType & fet,
377 const ElemType t)
378 { return n_dofs(fet, t); }
379
n_shape_functions(const FEType & fet,const Elem * inf_elem)380 static unsigned int n_shape_functions (const FEType & fet,
381 const Elem * inf_elem)
382 { return n_dofs(fet, inf_elem); }
383
384 /**
385 * \returns The number of shape functions associated with this
386 * infinite element. Currently, we have \p o_radial+1 modes in
387 * radial direction, and \code FE<Dim-1,T>::n_dofs(...) \endcode
388 * in the base.
389 *
390 * \deprecated Call the version of this function that takes an Elem*
391 * instead for consistency with other FEInterface::n_dofs() methods.
392 */
393 static unsigned int n_dofs(const FEType & fet,
394 const ElemType inf_elem_type);
395
396 static unsigned int n_dofs(const FEType & fet,
397 const Elem * inf_elem);
398
399 /**
400 * \returns The number of dofs at infinite element node \p n
401 * (not dof!) for an element of type \p t and order \p o.
402 *
403 * \deprecated Call the version of this function that takes an Elem*
404 * instead for consistency with other FEInterface::n_dofs() methods.
405 */
406 static unsigned int n_dofs_at_node(const FEType & fet,
407 const ElemType inf_elem_type,
408 const unsigned int n);
409
410 static unsigned int n_dofs_at_node(const FEType & fet,
411 const Elem * inf_elem,
412 const unsigned int n);
413
414 /**
415 * \returns The number of dofs interior to the element,
416 * not associated with any interior nodes.
417 *
418 * \deprecated Call the version of this function that takes an Elem*
419 * instead for consistency with other FEInterface::n_dofs() methods.
420 */
421 static unsigned int n_dofs_per_elem(const FEType & fet,
422 const ElemType inf_elem_type);
423
424 static unsigned int n_dofs_per_elem(const FEType & fet,
425 const Elem * inf_elem);
426
427 /**
428 * \returns The continuity of the element.
429 */
get_continuity()430 virtual FEContinuity get_continuity() const override
431 { return C_ZERO; } // FIXME - is this true??
432
433 /**
434 * \returns \p true if the element's higher order shape functions are
435 * hierarchic
436 */
is_hierarchic()437 virtual bool is_hierarchic() const override
438 { return false; } // FIXME - Inf FEs don't handle p elevation yet
439
440 /**
441 * Usually, this method would build the nodal soln from the
442 * element soln. But infinite elements require additional
443 * simulation-specific data to compute physically correct
444 * results. Use \p compute_data() to compute results. For
445 * compatibility an empty vector is returned.
446 */
447 static void nodal_soln(const FEType & fet,
448 const Elem * elem,
449 const std::vector<Number> & elem_soln,
450 std::vector<Number> & nodal_soln);
451
452
map(const Elem * inf_elem,const Point & reference_point)453 static Point map (const Elem * inf_elem,
454 const Point & reference_point) {
455 // libmesh_deprecated(); // soon
456 return InfFEMap::map(Dim, inf_elem, reference_point);
457 }
458
459
460 static Point inverse_map (const Elem * elem,
461 const Point & p,
462 const Real tolerance = TOLERANCE,
463 const bool secure = true) {
464 // libmesh_deprecated(); // soon
465 return InfFEMap::inverse_map(Dim, elem, p, tolerance, secure);
466 }
467
468
469 static void inverse_map (const Elem * elem,
470 const std::vector<Point> & physical_points,
471 std::vector<Point> & reference_points,
472 const Real tolerance = TOLERANCE,
473 const bool secure = true) {
474 // libmesh_deprecated(); // soon
475 return InfFEMap::inverse_map(Dim, elem, physical_points,
476 reference_points, tolerance, secure);
477 }
478
479
480 // The workhorses of InfFE. These are often used during matrix assembly.
481
482 /**
483 * This is at the core of this class. Use this for each
484 * new element in the mesh. Reinitializes all the physical
485 * element-dependent data based on the current element
486 * \p elem.
487 * \note: pts need to be in reference space coordinates, not physical ones.
488 */
489 virtual void reinit (const Elem * elem,
490 const std::vector<Point> * const pts = nullptr,
491 const std::vector<Real> * const weights = nullptr) override;
492
493 /**
494 * Reinitializes all the physical
495 * element-dependent data based on the \p side of an infinite
496 * element.
497 */
498 virtual void reinit (const Elem * inf_elem,
499 const unsigned int s,
500 const Real tolerance = TOLERANCE,
501 const std::vector<Point> * const pts = nullptr,
502 const std::vector<Real> * const weights = nullptr) override;
503
504 /**
505 * Not implemented yet. Reinitializes all the physical
506 * element-dependent data based on the \p edge of an infinite
507 * element.
508 */
509 virtual void edge_reinit (const Elem * elem,
510 const unsigned int edge,
511 const Real tolerance = TOLERANCE,
512 const std::vector<Point> * const pts = nullptr,
513 const std::vector<Real> * const weights = nullptr) override;
514
515 /**
516 * Computes the reference space quadrature points on the side of
517 * an element based on the side quadrature points.
518 */
side_map(const Elem *,const Elem *,const unsigned int,const std::vector<Point> &,std::vector<Point> &)519 virtual void side_map (const Elem * /* elem */,
520 const Elem * /* side */,
521 const unsigned int /* s */,
522 const std::vector<Point> & /* reference_side_points */,
523 std::vector<Point> & /* reference_points */) override
524 {
525 libmesh_not_implemented();
526 }
527
528 /**
529 * The use of quadrature rules with the \p InfFE class is somewhat
530 * different from the approach of the \p FE class. While the
531 * \p FE class requires an appropriately initialized quadrature
532 * rule object, and simply uses it, the \p InfFE class requires only
533 * the quadrature rule object of the current \p FE class.
534 * From this \p QBase *, it determines the necessary data,
535 * and builds two appropriate quadrature classes, one for radial,
536 * and another for base integration, using the convenient
537 * \p QBase::build() method.
538 */
539 virtual void attach_quadrature_rule (QBase * q) override;
540
541 /**
542 * \returns The number of shape functions associated with
543 * this infinite element.
544 */
n_shape_functions()545 virtual unsigned int n_shape_functions () const override
546 { return _n_total_approx_sf; }
547
548 /**
549 * \returns The total number of quadrature points. Call this
550 * to get an upper bound for the \p for loop in your simulation
551 * for matrix assembly of the current element.
552 */
n_quadrature_points()553 virtual unsigned int n_quadrature_points () const override
554 { libmesh_assert(radial_qrule); return _n_total_qp; }
555
556
557 protected:
558
559 // static members used by the workhorses
560
561 /**
562 * \returns The value of the \f$ i^{th} \f$ polynomial evaluated
563 * at \p v. This method provides the approximation
564 * in radial direction for the overall shape functions,
565 * which is defined in \p InfFE::shape().
566 * This method is allowed to be static, since it is independent
567 * of dimension and base_family. It is templated, though,
568 * w.r.t. to radial \p FEFamily.
569 *
570 * \returns The value of the \f$ i^{th} \f$ mapping shape function
571 * in radial direction evaluated at \p v when T_radial ==
572 * INFINITE_MAP. Currently, only one specific mapping shape is
573 * used. Namely the one by Marques JMMC, Owen DRJ: Infinite
574 * elements in quasi-static materially nonlinear problems, Computers
575 * and Structures, 1984.
576 */
577 static Real eval(Real v,
578 Order o_radial,
579 unsigned int i);
580
581 /**
582 * \returns The value of the first derivative of the
583 * \f$ i^{th} \f$ polynomial at coordinate \p v.
584 * See \p eval for details.
585 */
586 static Real eval_deriv(Real v,
587 Order o_radial,
588 unsigned int i);
589
590
591
592 // Non-static members used by the workhorses
593
594 /**
595 * Updates the protected member \p base_elem to the appropriate base element
596 * for the given \p inf_elem.
597 */
598 void update_base_elem (const Elem * inf_elem);
599
600 /**
601 * Do not use this derived member in \p InfFE<Dim,T_radial,T_map>.
602 */
init_base_shape_functions(const std::vector<Point> &,const Elem *)603 virtual void init_base_shape_functions(const std::vector<Point> &,
604 const Elem *) override
605 { libmesh_not_implemented(); }
606
607 /**
608 * Some of the member data only depend on the radial part of the
609 * infinite element. The parts that only change when the radial
610 * order changes, are initialized here.
611 */
612 void init_radial_shape_functions(const Elem * inf_elem,
613 const std::vector<Point> * radial_pts = nullptr);
614
615 /**
616 * Initialize all the data fields like \p weight, \p mode,
617 * \p phi, \p dphidxi, \p dphideta, \p dphidzeta, etc.
618 * for the current element. This method prepares the data
619 * related to the base part, and some of the combined fields.
620 */
621 void init_shape_functions(const std::vector<Point> & radial_qp,
622 const std::vector<Point> & base_qp,
623 const Elem * inf_elem);
624
625 /**
626 * Initialize all the data fields like \p weight,
627 * \p phi, etc for the side \p s.
628 */
629 void init_face_shape_functions (const std::vector<Point> &,
630 const Elem * inf_side);
631
632 /**
633 * Combines the shape functions, which were formed in
634 * \p init_shape_functions(Elem *), with geometric data.
635 * Has to be called every time the geometric configuration
636 * changes. Afterward, the fields are ready to be used
637 * to compute global derivatives, the jacobian etc, see
638 * \p FEAbstract::compute_map().
639 */
640 void combine_base_radial(const Elem * inf_elem);
641
642 /**
643 * After having updated the jacobian and the transformation
644 * from local to global coordinates in FEAbstract::compute_map(),
645 * the first derivatives of the shape functions are
646 * transformed to global coordinates, giving \p dphi,
647 * \p dphidx/y/z, \p dphasedx/y/z, \p dweight. This method
648 * should barely be re-defined in derived classes, but
649 * still should be usable for children. Therefore, keep
650 * it protected.
651 */
652 virtual void compute_shape_functions(const Elem *, const std::vector<Point> &) override;
653
654
655
656 // Miscellaneous static members
657
658 /**
659 * Computes the indices in the base \p base_node and in radial
660 * direction \p radial_node (either 0 or 1) associated to the
661 * node \p outer_node_index of an infinite element of type
662 * \p inf_elem_type.
663 */
664 static void compute_node_indices (const ElemType inf_elem_type,
665 const unsigned int outer_node_index,
666 unsigned int & base_node,
667 unsigned int & radial_node);
668
669 /**
670 * Does the same as \p compute_node_indices(), but stores
671 * the maps for the current element type. Provided the
672 * infinite element type changes seldom, this is probably
673 * faster than using \p compute_node_indices () alone.
674 * This is possible since the number of nodes is not likely
675 * to change.
676 */
677 static void compute_node_indices_fast (const ElemType inf_elem_type,
678 const unsigned int outer_node_index,
679 unsigned int & base_node,
680 unsigned int & radial_node);
681
682 /**
683 * Computes the indices of shape functions in the base \p base_shape and
684 * in radial direction \p radial_shape (0 in the base, \f$ \ge 1 \f$ further
685 * out) associated to the shape with global index \p i of an infinite element
686 * of type \p inf_elem_type.
687 *
688 * \deprecated Call the version of this function that takes an Elem * instead.
689 */
690 static void compute_shape_indices (const FEType & fet,
691 const ElemType inf_elem_type,
692 const unsigned int i,
693 unsigned int & base_shape,
694 unsigned int & radial_shape);
695
696 static void compute_shape_indices (const FEType & fet,
697 const Elem * inf_elem,
698 const unsigned int i,
699 unsigned int & base_shape,
700 unsigned int & radial_shape);
701
702 /**
703 * the radial distance of the base nodes from the origin
704 */
705 std::vector<Real> dist;
706
707 /**
708 * the additional radial weight \f$ 1/{r^2} \f$ in local coordinates,
709 * over all quadrature points. The weight does not vary in base
710 * direction. However, for uniform access to the data fields from the
711 * outside, this data field is expanded to all quadrature points.
712 */
713 std::vector<Real> dweightdv;
714
715 /**
716 * the radial decay \f$ 1/r \f$ in local coordinates. Needed when
717 * setting up the overall shape functions.
718 *
719 * \note It is this decay which ensures that the Sommerfeld
720 * radiation condition is satisfied in advance.
721 */
722 std::vector<Real> som;
723 /**
724 * the first local derivative of the radial decay \f$ 1/r \f$ in local
725 * coordinates. Needed when setting up the overall shape functions.
726 */
727 std::vector<Real> dsomdv;
728
729 /**
730 * the radial approximation shapes in local coordinates
731 * Needed when setting up the overall shape functions.
732 */
733 std::vector<std::vector<Real>> mode;
734
735 /**
736 * the first local derivative of the radial approximation shapes.
737 * Needed when setting up the overall shape functions.
738 */
739 std::vector<std::vector<Real>> dmodedv;
740
741 /**
742 * the radial mapping shapes in local coordinates
743 */
744 std::vector<std::vector<Real>> radial_map;
745
746
747 /**
748 * the first local derivative of the radial mapping shapes
749 */
750 std::vector<std::vector<Real>> dradialdv_map;
751
752 /**
753 * the first local derivative (for 3D, the first in the base)
754 * of the phase term in local coordinates.
755 * Needed in the overall weak form of infinite element formulations.
756 */
757 std::vector<Real> dphasedxi;
758
759 /**
760 * the second local derivative (for 3D, the second in the base)
761 * of the phase term in local coordinates.
762 * Needed in the overall weak form of infinite element formulations.
763 */
764 std::vector<Real> dphasedeta;
765
766 /**
767 * the third local derivative (for 3D, the derivative in radial
768 * direction) of the phase term in local coordinates.
769 * Needed in the overall weak form of infinite element formulations.
770 */
771 std::vector<Real> dphasedzeta;
772
773
774
775
776 // numbering scheme maps
777
778 /**
779 * The internal structure of the \p InfFE
780 * -- tensor product of base element times radial
781 * nodes -- has to be determined from the node numbering
782 * of the current infinite element. This vector
783 * maps the infinite \p Elem node number to the
784 * radial node (either 0 or 1).
785 */
786 std::vector<unsigned int> _radial_node_index;
787
788 /**
789 * The internal structure of the \p InfFE
790 * -- tensor product of base element times radial
791 * nodes -- has to be determined from the node numbering
792 * of the current element. This vector
793 * maps the infinite \p Elem node number to the
794 * associated node in the base element.
795 */
796 std::vector<unsigned int> _base_node_index;
797
798 /**
799 * The internal structure of the \p InfFE
800 * -- tensor product of base element shapes times radial
801 * shapes -- has to be determined from the dof numbering
802 * scheme of the current infinite element. This vector
803 * maps the infinite \p Elem dof index to the radial
804 * \p InfFE shape index (\p 0..radial_order+1 ).
805 */
806 std::vector<unsigned int> _radial_shape_index;
807
808 /**
809 * The internal structure of the \p InfFE
810 * -- tensor product of base element shapes times radial
811 * shapes -- has to be determined from the dof numbering
812 * scheme of the current infinite element. This vector
813 * maps the infinite \p Elem dof index to the associated
814 * dof in the base \p FE.
815 */
816 std::vector<unsigned int> _base_shape_index;
817
818
819
820
821 // some more protected members
822
823 /**
824 * The number of total approximation shape functions for
825 * the current configuration
826 */
827 unsigned int _n_total_approx_sf;
828
829 /**
830 * The total number of quadrature points
831 * for the current configuration
832 */
833 unsigned int _n_total_qp;
834
835 /**
836 * this vector contains the combined integration weights, so
837 * that \p FEAbstract::compute_map() can still be used
838 */
839 std::vector<Real> _total_qrule_weights;
840
841 /**
842 * The quadrature rule for the base element associated
843 * with the current infinite element
844 */
845 std::unique_ptr<QBase> base_qrule;
846
847 /**
848 * The quadrature rule for the base element associated
849 * with the current infinite element
850 */
851 std::unique_ptr<QBase> radial_qrule;
852
853 /**
854 * The base element associated with the
855 * current infinite element
856 */
857 std::unique_ptr<Elem> base_elem;
858
859 /**
860 * Have a \p FE<Dim-1,T_base> handy for base approximation.
861 * Since this one is created using the \p FEBase::build() method,
862 * the \p InfFE class is not required to be templated w.r.t.
863 * to the base approximation shape.
864 */
865 std::unique_ptr<FEBase> base_fe;
866
867 /**
868 * This \p FEType stores the characteristics for which the data
869 * structures \p phi, \p phi_map etc are currently initialized.
870 * This avoids re-initializing the radial part.
871 *
872 * \note Currently only \p order may change, both the FE families
873 * and \p base_order must remain constant.
874 */
875 FEType current_fe_type;
876
877
878 private:
879
880 /**
881 * \returns \p false, currently not required.
882 */
883 virtual bool shapes_need_reinit() const override;
884
885 /**
886 * When \p compute_node_indices_fast() is used, this static
887 * variable remembers the element type for which the
888 * static variables in \p compute_node_indices_fast()
889 * are currently set. Using a class member for the
890 * element type helps initializing it to a default value.
891 */
892 static ElemType _compute_node_indices_fast_current_elem_type;
893
894
895 #ifdef DEBUG
896
897 /**
898 * static members that are used to issue warning messages only once.
899 */
900 static bool _warned_for_nodal_soln;
901 static bool _warned_for_shape;
902 static bool _warned_for_dshape;
903
904 #endif
905
906 /**
907 * Make all \p InfFE<Dim,T_radial,T_map> classes
908 * friends of each other, so that the protected
909 * \p eval() may be accessed.
910 */
911 template <unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
912 friend class InfFE;
913
914 friend class InfFEMap;
915 };
916
917
918
919 // InfFERadial class inline members
920 inline
decay(const unsigned int dim,const Real v)921 Real InfFERadial::decay(const unsigned int dim, const Real v)
922 {
923 switch (dim)
924 //TODO:[DD] What decay do i have in 2D and 1D?
925 {
926 case 3:
927 return (1.-v)/2.;
928
929 case 2:
930 return 0.;
931
932 case 1:
933 return 0.;
934
935 default:
936 libmesh_error_msg("Invalid dim = " << dim);
937 }
938 }
939
940 } // namespace libMesh
941
942
943 #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
944
945
946 #endif // LIBMESH_INF_FE_H
947