Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 : */
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 : static Real decay_deriv (const unsigned int dim, const Real);
75 :
76 : /**
77 : * \returns The radial weight D, used as an additional weight
78 : * for the test function, evaluated at local radial coordinate \p v.
79 : * FIXME: This is only valid for 3D!!
80 : * And also makes assumptions on the mapping
81 : */
82 2961 : static Real D (const Real v) { return (1.-v)*(1.-v)/4.; }
83 :
84 30648 : static Real Dxr_sq (const Real /*v*/) { return 1.; }
85 :
86 : /**
87 : * \returns The first (local) radial derivative of the radial weight D.
88 : * FIXME: this is valid only in 3D.
89 : */
90 4942 : static Real D_deriv (const Real v) { return (v-1.)/2.; }
91 :
92 :
93 : /**
94 : * \returns The Order of the mapping functions
95 : * in the radial direction. Currently, this is always \p FIRST.
96 : */
97 0 : static Order mapping_order() { return FIRST; }
98 :
99 : /**
100 : * \returns The number of shape functions in radial direction
101 : * associated with this infinite element.
102 : * Either way, if the modes are stored as nodal dofs (\p n_dofs_at_node)
103 : * or as element dofs (\p n_dofs_per_elem), in each case we have the
104 : * same number of modes in radial direction.
105 : *
106 : * \note For the case of 1D infinite elements, in the base the
107 : * dof-per-node scheme is used.
108 : *
109 : * From the formulation of the infinite elements, we have
110 : * 1 mode, when \p o_radial=CONST.
111 : * Therefore, we have a total of \p o_radial+1 modes in radial direction.
112 : */
113 80716 : static unsigned int n_dofs (const Order o_radial)
114 234830 : { return static_cast<unsigned int>(o_radial)+1; }
115 :
116 : /**
117 : * \returns The number of dofs in radial direction on "onion slice"
118 : * \p n (either 0 or 1) for an infinite element of type \p inf_elem_type and
119 : * radial order \p o_radial.
120 : *
121 : * Currently, the first radial mode is associated with the nodes in
122 : * the base. All higher radial modes are associated with
123 : * the physically existing nodes further out.
124 : */
125 : static unsigned int n_dofs_at_node (const Order o_radial,
126 : const unsigned int n_onion);
127 :
128 : /**
129 : * \returns The number of modes in radial direction interior to the element,
130 : * not associated with any interior nodes.
131 : *
132 : * \note These modes are a discontinuous approximation, therefore
133 : * we have no special formulation for coupling in the base, like in the
134 : * case of associating (possibly) multiple dofs per (outer) node.
135 : */
136 8515 : static unsigned int n_dofs_per_elem (const Order o_radial)
137 22750 : { return static_cast<unsigned int>(o_radial)+1; }
138 :
139 : };
140 :
141 :
142 :
143 : /**
144 : * This nested class contains most of the static methods related
145 : * to the base part of an infinite element. Only static members
146 : * are provided, for use within \p InfFE and \p InfFEMap.
147 : *
148 : * \author Daniel Dreyer
149 : * \date 2003
150 : */
151 : class InfFEBase
152 : {
153 : private:
154 :
155 : /**
156 : * Never use an object of this type.
157 : */
158 : InfFEBase() {}
159 :
160 : public:
161 :
162 : /**
163 : * Build the base element of an infinite element. Here "base" refers
164 : * to the non-infinite/traditional geometric element face which is
165 : * associated with this infinite element. A const base Elem is built
166 : * from a const input InfElem.
167 : */
168 : static std::unique_ptr<const Elem> build_elem (const Elem * inf_elem);
169 :
170 : /**
171 : * \returns The base element associated to
172 : * \p type. This is, for example, \p TRI3 for
173 : * \p INFPRISM6.
174 : */
175 : static ElemType get_elem_type (const ElemType type);
176 :
177 : /**
178 : * \returns The number of shape functions used in the
179 : * mapping in the base element of type \p base_elem_type
180 : * mapped with order \p base_mapping_order
181 : */
182 : static unsigned int n_base_mapping_sf (const Elem & base_elem,
183 : const Order base_mapping_order);
184 :
185 : };
186 :
187 :
188 :
189 : /**
190 : * A specific instantiation of the \p FEBase class. This
191 : * class is templated, and specific template instantiations
192 : * will result in different Infinite Element families, similar
193 : * to the \p FE class. \p InfFE builds a \p FE<Dim-1,T_base>,
194 : * and most of the requests related to the base are handed over
195 : * to this object. All methods related to the radial part
196 : * are collected in the class \p InfFERadial. Similarly,
197 : * most of the static methods concerning base approximation
198 : * are contained in \p InfFEBase.
199 : *
200 : * Having different shape approximation families in radial direction
201 : * introduces the requirement for an additional \p Order in this
202 : * class. Therefore, the \p FEType internals change when infinite
203 : * elements are enabled.
204 : * When the specific infinite element type is not known at compile
205 : * time, use the \p FEBase::build() member to create abstract
206 : * (but still optimized) infinite elements at run time.
207 : *
208 : * The node numbering scheme is the one from the current
209 : * infinite element. Each node in the base holds exactly
210 : * the same number of dofs as an adjacent conventional \p FE
211 : * would contain. The nodes further out hold the additional
212 : * dof necessary for radial approximation. The order of the outer nodes'
213 : * components is such that the radial shapes have highest
214 : * priority, followed by the base shapes.
215 : *
216 : * For the derivation and some background of the implemented algorithm
217 : * see https://arxiv.org/abs/2501.05568.
218 : *
219 : * \author Daniel Dreyer
220 : * \date 2003
221 : * \brief Base class for all the infinite geometric element types.
222 : */
223 : template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
224 : class InfFE : public FEBase
225 : {
226 :
227 : /*
228 : * Protect the nested class
229 : */
230 : protected:
231 :
232 :
233 : public:
234 :
235 : // InfFE continued
236 :
237 : /**
238 : * Constructor and empty destructor.
239 : * Initializes some data structures. Builds a \p FE<Dim-1,T_base>
240 : * object to handle approximation in the base, so that
241 : * there is no need to template \p InfFE<Dim,T_radial,T_map> also with
242 : * respect to the base approximation \p T_base.
243 : *
244 : * The same remarks concerning compile-time optimization for
245 : * \p FE also hold for \p InfFE. Use the
246 : * \p FEBase::build_InfFE(const unsigned int, const FEType &)
247 : * method to build specific instantiations of \p InfFE at
248 : * run time.
249 : */
250 : explicit
251 : InfFE(const FEType & fet);
252 2140 : ~InfFE() = default;
253 :
254 : // The static public members for access from FEInterface etc
255 :
256 : /**
257 : * \returns The value of the \f$ i^{th} \f$ shape function at
258 : * point \p p. This method lets you specify the relevant
259 : * data directly, and is therefore allowed to be static.
260 : *
261 : * \note This class member is not as efficient as its counterpart in
262 : * \p FE<Dim,T>, and is not employed in the \p reinit() cycle.
263 : *
264 : * \note This method does not return physically correct shapes,
265 : * instead use \p compute_data(). The \p shape() methods should
266 : * only be used for mapping.
267 : */
268 : static Real shape(const FEType & fet,
269 : const ElemType t,
270 : const unsigned int i,
271 : const Point & p);
272 :
273 : /**
274 : * \returns The value of the \f$ i^{th} \f$ shape function at
275 : * point \p p. This method lets you specify the relevant
276 : * data directly, and is therefore allowed to be static.
277 : *
278 : * \note This class member is not as efficient as its counterpart in
279 : * \p FE<Dim,T>, and is not employed in the \p reinit() cycle.
280 : *
281 : * \note This method does not return physically correct shapes,
282 : * instead use \p compute_data(). The \p shape() methods should
283 : * only be used for mapping.
284 : */
285 : static Real shape(const FEType & fet,
286 : const Elem * elem,
287 : const unsigned int i,
288 : const Point & p);
289 :
290 : /**
291 : * \returns The value of the \f$ i^{th} \f$ shape function at
292 : * point \p p. This method lets you specify the relevant
293 : * data directly, and is therefore allowed to be static.
294 : *
295 : * \note This class member is not as efficient as its counterpart in
296 : * \p FE<Dim,T>, and is not employed in the \p reinit() cycle.
297 : *
298 : * \note This method does not return physically correct shapes,
299 : * instead use \p compute_data(). The \p shape() methods should
300 : * only be used for mapping.
301 : */
302 : static Real shape(const FEType fet,
303 : const Elem * elem,
304 : const unsigned int i,
305 : const Point & p,
306 : const bool add_p_level);
307 :
308 : /**
309 : * \returns The \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
310 : * shape function at point \p p. This method lets you specify the relevant
311 : * data directly, and is therefore allowed to be static.
312 : *
313 : * \note This class member is not as efficient as its counterpart in
314 : * \p FE<Dim,T>, and is not employed in the \p reinit() cycle.
315 : *
316 : * \note This method does not return physically correct shape gradients,
317 : * instead use \p compute_data(). The \p shape_deriv() methods should
318 : * only be used for mapping.
319 : */
320 : static Real shape_deriv (const FEType & fet,
321 : const Elem * inf_elem,
322 : const unsigned int i,
323 : const unsigned int j,
324 : const Point & p);
325 :
326 : /**
327 : * \returns The \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
328 : * shape function at point \p p. This method lets you specify the relevant
329 : * data directly, and is therefore allowed to be static.
330 : *
331 : * \note This class member is not as efficient as its counterpart in
332 : * \p FE<Dim,T>, and is not employed in the \p reinit() cycle.
333 : *
334 : * \note This method does not return physically correct shape gradients,
335 : * instead use \p compute_data(). The \p shape_deriv() methods should
336 : * only be used for mapping.
337 : */
338 : static Real shape_deriv (const FEType fet,
339 : const Elem * inf_elem,
340 : const unsigned int i,
341 : const unsigned int j,
342 : const Point & p,
343 : const bool add_p_level);
344 :
345 : /**
346 : * \returns The \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
347 : * shape function at point \p p. This method lets you specify the relevant
348 : * data directly, and is therefore allowed to be static.
349 : *
350 : * \note This class member is not as efficient as its counterpart in
351 : * \p FE<Dim,T>, and is not employed in the \p reinit() cycle.
352 : *
353 : * \note This method does not return physically correct shape gradients,
354 : * instead use \p compute_data(). The \p shape_deriv() methods should
355 : * only be used for mapping.
356 : */
357 : static Real shape_deriv (const FEType & fet,
358 : const ElemType inf_elem_type,
359 : const unsigned int i,
360 : const unsigned int j,
361 : const Point & p);
362 :
363 : /**
364 : * Generalized version of \p shape(), takes an \p Elem *. The \p data
365 : * contains both input and output parameters. For frequency domain
366 : * simulations, the complex-valued shape is returned. In time domain
367 : * both the computed shape, and the phase is returned.
368 : *
369 : * \note The phase (proportional to the distance of the \p Point \p
370 : * data.p from the envelope) is actually a measure how far into the
371 : * future the results are.
372 : */
373 : static void compute_data(const FEType & fe_t,
374 : const Elem * inf_elem,
375 : FEComputeData & data);
376 :
377 : /**
378 : * \returns The number of shape functions associated with
379 : * a finite element of type \p t and approximation order \p o.
380 : */
381 0 : static unsigned int n_shape_functions (const FEType & fet,
382 : const Elem * inf_elem)
383 0 : { return n_dofs(fet, inf_elem); }
384 :
385 : /**
386 : * \returns The number of shape functions associated with this
387 : * infinite element. Currently, we have \p o_radial+1 modes in
388 : * radial direction, and \code FE<Dim-1,T>::n_dofs(...) \endcode
389 : * in the base.
390 : */
391 : static unsigned int n_dofs(const FEType & fet,
392 : const Elem * inf_elem);
393 :
394 : #ifdef LIBMESH_ENABLE_DEPRECATED
395 : /**
396 : * \returns The number of dofs at infinite element node \p n
397 : * (not dof!) for an element of type \p t and order \p o.
398 : */
399 : static unsigned int n_dofs_at_node(const FEType & fet,
400 : const ElemType inf_elem_type,
401 : const unsigned int n);
402 : #endif // LIBMESH_ENABLE_DEPRECATED
403 :
404 : /**
405 : * \returns The number of dofs at infinite element node \p n
406 : * (not dof!) for an element of type \p t and order \p o.
407 : */
408 : static unsigned int n_dofs_at_node(const FEType & fet,
409 : const Elem * inf_elem,
410 : const unsigned int n);
411 :
412 : #ifdef LIBMESH_ENABLE_DEPRECATED
413 : /**
414 : * \returns The number of dofs interior to the element,
415 : * not associated with any interior nodes.
416 : *
417 : * \deprecated Call the version of this function that takes an Elem*
418 : * instead for consistency with other FEInterface::n_dofs() methods.
419 : */
420 : static unsigned int n_dofs_per_elem(const FEType & fet,
421 : const ElemType inf_elem_type);
422 : #endif // LIBMESH_ENABLE_DEPRECATED
423 :
424 : /**
425 : * \returns The number of dofs interior to the element,
426 : * not associated with any interior nodes.
427 : */
428 : static unsigned int n_dofs_per_elem(const FEType & fet,
429 : const Elem * inf_elem);
430 :
431 : /**
432 : * \returns The continuity of the element.
433 : */
434 0 : virtual FEContinuity get_continuity() const override
435 0 : { return C_ZERO; } // FIXME - is this true??
436 :
437 : /**
438 : * \returns \p true if the element's higher order shape functions are
439 : * hierarchic
440 : */
441 0 : virtual bool is_hierarchic() const override
442 0 : { return false; } // FIXME - Inf FEs don't handle p elevation yet
443 :
444 : /**
445 : * Usually, this method would build the nodal soln from the
446 : * element soln. But infinite elements require additional
447 : * simulation-specific data to compute physically correct
448 : * results. Use \p compute_data() to compute results. For
449 : * compatibility an empty vector is returned.
450 : */
451 : static void nodal_soln(const FEType & fet,
452 : const Elem * elem,
453 : const std::vector<Number> & elem_soln,
454 : std::vector<Number> & nodal_soln);
455 :
456 :
457 0 : static Point map (const Elem * inf_elem,
458 : const Point & reference_point)
459 : {
460 : // libmesh_deprecated(); // soon
461 0 : return InfFEMap::map(Dim, inf_elem, reference_point);
462 : }
463 :
464 :
465 0 : static Point inverse_map (const Elem * elem,
466 : const Point & p,
467 : const Real tolerance = TOLERANCE,
468 : const bool secure = true)
469 : {
470 : // libmesh_deprecated(); // soon
471 0 : return InfFEMap::inverse_map(Dim, elem, p, tolerance, secure);
472 : }
473 :
474 :
475 0 : static void inverse_map (const Elem * elem,
476 : const std::vector<Point> & physical_points,
477 : std::vector<Point> & reference_points,
478 : const Real tolerance = TOLERANCE,
479 : const bool secure = true)
480 : {
481 : // libmesh_deprecated(); // soon
482 0 : return InfFEMap::inverse_map(Dim, elem, physical_points,
483 0 : reference_points, tolerance, secure);
484 : }
485 :
486 :
487 : // The workhorses of InfFE. These are often used during matrix assembly.
488 :
489 : /**
490 : * This is at the core of this class. Use this for each
491 : * new element in the mesh. Reinitializes all the physical
492 : * element-dependent data based on the current element
493 : * \p elem.
494 : * \note pts need to be in reference space coordinates, not physical ones.
495 : */
496 : virtual void reinit (const Elem * elem,
497 : const std::vector<Point> * const pts = nullptr,
498 : const std::vector<Real> * const weights = nullptr) override;
499 :
500 : /**
501 : * Reinitializes all the physical
502 : * element-dependent data based on the \p side of an infinite
503 : * element.
504 : */
505 : virtual void reinit (const Elem * inf_elem,
506 : const unsigned int s,
507 : const Real tolerance = TOLERANCE,
508 : const std::vector<Point> * const pts = nullptr,
509 : const std::vector<Real> * const weights = nullptr) override;
510 :
511 : /**
512 : * Not implemented yet. Reinitializes all the physical
513 : * element-dependent data based on the \p edge of an infinite
514 : * element.
515 : */
516 : virtual void edge_reinit (const Elem * elem,
517 : const unsigned int edge,
518 : const Real tolerance = TOLERANCE,
519 : const std::vector<Point> * const pts = nullptr,
520 : const std::vector<Real> * const weights = nullptr) override;
521 :
522 : /**
523 : * Computes the reference space quadrature points on the side of
524 : * an element based on the side quadrature points.
525 : */
526 0 : virtual void side_map (const Elem * /* elem */,
527 : const Elem * /* side */,
528 : const unsigned int /* s */,
529 : const std::vector<Point> & /* reference_side_points */,
530 : std::vector<Point> & /* reference_points */) override
531 : {
532 0 : libmesh_not_implemented();
533 : }
534 :
535 : /**
536 : * The use of quadrature rules with the \p InfFE class is somewhat
537 : * different from the approach of the \p FE class. While the
538 : * \p FE class requires an appropriately initialized quadrature
539 : * rule object, and simply uses it, the \p InfFE class requires only
540 : * the quadrature rule object of the current \p FE class.
541 : * From this \p QBase *, it determines the necessary data,
542 : * and builds two appropriate quadrature classes, one for radial,
543 : * and another for base integration, using the convenient
544 : * \p QBase::build() method.
545 : */
546 : virtual void attach_quadrature_rule (QBase * q) override;
547 :
548 : /**
549 : * \returns The number of shape functions associated with
550 : * this infinite element.
551 : */
552 0 : virtual unsigned int n_shape_functions () const override
553 0 : { return _n_total_approx_sf; }
554 :
555 : /**
556 : * \returns The total number of quadrature points. Call this
557 : * to get an upper bound for the \p for loop in your simulation
558 : * for matrix assembly of the current element.
559 : */
560 2592 : virtual unsigned int n_quadrature_points () const override
561 2592 : { libmesh_assert(radial_qrule); return this->_n_total_qp; }
562 :
563 : /**
564 : * \returns the \p xyz spatial locations of the quadrature
565 : * points on the element.
566 : */
567 15 : virtual const std::vector<Point> & get_xyz () const override
568 6 : { libmesh_assert(!calculations_started || calculate_xyz);
569 15 : calculate_xyz = true; return xyz; }
570 :
571 : /**
572 : * \returns the Jacobian times quadrature weight.
573 : * Due to the divergence with increasing radial distance, this quantity
574 : * is numerically unstable.
575 : * Thus, it is safer to use \p get_JxWxdecay_sq() instead!
576 : */
577 0 : virtual const std::vector<Real> & get_JxW () const override
578 : {
579 0 : libmesh_assert(!calculations_started || calculate_jxw);
580 0 : calculate_jxw = true;
581 0 : return this->JxW;
582 : }
583 :
584 : /**
585 : * \returns Jacobian times quadrature weight times square of the
586 : * decaying function \f$ decay= r^{-\frac{dim+1}{2}}\f$
587 : *
588 : * This function is the variant of \p get_JxW() for \p InfFE.
589 : * Since J diverges there, a respectize decay-function must be
590 : * applied to obtain well-defined quantities.
591 : */
592 2597 : virtual const std::vector<Real> & get_JxWxdecay_sq () const override
593 866 : { libmesh_assert(!calculations_started || calculate_map_scaled);
594 2597 : calculate_map_scaled = true; return this->JxWxdecay;}
595 :
596 :
597 : /**
598 : * \returns The shape function \p phi weighted by r/decay
599 : * where \f$ decay = r^{-\frac{dim+1}{2}} \f$
600 : *
601 : * To compensate for the decay function applied to the Jacobian (see \p get_JxWxdecay_sq),
602 : * the wave function \p phi should be divided by this function.
603 : *
604 : * The factor r must be compensated for by the Sobolev \p weight.
605 : * (i.e. by using \p get_Sobolev_weightxR_sq())
606 : **/
607 2597 : virtual const std::vector<std::vector<OutputShape>> & get_phi_over_decayxR () const override
608 866 : { libmesh_assert(!calculations_started || calculate_phi_scaled);
609 2597 : calculate_phi_scaled = true; return phixr; }
610 :
611 :
612 : /**
613 : * \returns the gradient of the shape function (see \p get_dphi()),
614 : * but in case of \p InfFE, weighted with r/decay.
615 : * See \p get_phi_over_decayxR() for details.
616 : */
617 2597 : virtual const std::vector<std::vector<OutputGradient>> & get_dphi_over_decayxR () const override
618 866 : { libmesh_assert(!calculations_started || calculate_dphi_scaled);
619 2597 : calculate_dphi_scaled = true; return dphixr; }
620 :
621 :
622 : /**
623 : * \returns the gradient of the shape function (see \p get_dphi()),
624 : * but in case of \p InfFE, weighted with 1/decay.
625 : *
626 : * In contrast to the shape function, its gradient stays finite
627 : * when divided by the decay function.
628 : */
629 0 : virtual const std::vector<std::vector<OutputGradient>> & get_dphi_over_decay () const override
630 0 : { libmesh_assert(!calculations_started || calculate_dphi_scaled);
631 0 : calculate_dphi_scaled = true; return dphixr_sq; }
632 :
633 :
634 : /**
635 : * \returns The element tangents in xi-direction at the quadrature
636 : * points.
637 : */
638 0 : virtual const std::vector<RealGradient> & get_dxyzdxi() const override
639 0 : { calculate_map = true; libmesh_not_implemented();}
640 :
641 :
642 : /**
643 : * \returns The element tangents in eta-direction at the quadrature
644 : * points.
645 : */
646 0 : virtual const std::vector<RealGradient> & get_dxyzdeta() const override
647 0 : { calculate_map = true; libmesh_not_implemented();}
648 :
649 :
650 : /**
651 : * \returns The element tangents in zeta-direction at the quadrature
652 : * points.
653 : */
654 0 : virtual const std::vector<RealGradient> & get_dxyzdzeta() const override
655 0 : { calculate_map = true; libmesh_not_implemented();}
656 :
657 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
658 : /**
659 : * \returns The second partial derivatives in xi.
660 : */
661 0 : virtual const std::vector<RealGradient> & get_d2xyzdxi2() const override
662 0 : { calculate_map = true; libmesh_not_implemented();}
663 : /**
664 : * \returns The second partial derivatives in eta.
665 : */
666 0 : virtual const std::vector<RealGradient> & get_d2xyzdeta2() const override
667 0 : { calculate_map = true; libmesh_not_implemented();}
668 : /**
669 : * \returns The second partial derivatives in zeta.
670 : */
671 0 : virtual const std::vector<RealGradient> & get_d2xyzdzeta2() const override
672 0 : { calculate_map = true; libmesh_not_implemented();}
673 : /**
674 : * \returns The second partial derivatives in xi-eta.
675 : */
676 0 : virtual const std::vector<RealGradient> & get_d2xyzdxideta() const override
677 0 : { calculate_map = true; libmesh_not_implemented();}
678 : /**
679 : * \returns The second partial derivatives in xi-zeta.
680 : */
681 0 : virtual const std::vector<RealGradient> & get_d2xyzdxidzeta() const override
682 0 : { calculate_map = true; libmesh_not_implemented();}
683 : /**
684 : * \returns The second partial derivatives in eta-zeta.
685 : */
686 0 : virtual const std::vector<RealGradient> & get_d2xyzdetadzeta() const override
687 0 : { calculate_map = true; libmesh_not_implemented();}
688 : #endif
689 :
690 : /**
691 : * \returns The dxi/dx entry in the transformation
692 : * matrix from physical to local coordinates.
693 : */
694 5 : virtual const std::vector<Real> & get_dxidx() const override
695 2 : {libmesh_assert(!calculations_started || calculate_map);
696 5 : calculate_map = true; return dxidx_map;}
697 :
698 :
699 : /**
700 : * \returns The dxi/dy entry in the transformation
701 : * matrix from physical to local coordinates.
702 : */
703 5 : virtual const std::vector<Real> & get_dxidy() const override
704 2 : {libmesh_assert(!calculations_started || calculate_map);
705 5 : calculate_map = true; return dxidy_map;}
706 :
707 :
708 : /**
709 : * \returns The dxi/dz entry in the transformation
710 : * matrix from physical to local coordinates.
711 : */
712 5 : virtual const std::vector<Real> & get_dxidz() const override
713 2 : { libmesh_assert(!calculations_started || calculate_map);
714 5 : calculate_map = true; return dxidz_map;}
715 :
716 :
717 : /**
718 : * \returns The deta/dx entry in the transformation
719 : * matrix from physical to local coordinates.
720 : */
721 5 : virtual const std::vector<Real> & get_detadx() const override
722 2 : { libmesh_assert(!calculations_started || calculate_map);
723 5 : calculate_map = true; return detadx_map;}
724 :
725 :
726 : /**
727 : * \returns The deta/dy entry in the transformation
728 : * matrix from physical to local coordinates.
729 : */
730 5 : virtual const std::vector<Real> & get_detady() const override
731 2 : { libmesh_assert(!calculations_started || calculate_map);
732 5 : calculate_map = true; return detady_map;}
733 :
734 :
735 : /**
736 : * \returns The deta/dx entry in the transformation
737 : * matrix from physical to local coordinates.
738 : */
739 5 : virtual const std::vector<Real> & get_detadz() const override
740 2 : { libmesh_assert(!calculations_started || calculate_map);
741 5 : calculate_map = true; return detadz_map;}
742 :
743 :
744 : /**
745 : * \returns The dzeta/dx entry in the transformation
746 : * matrix from physical to local coordinates.
747 : */
748 10 : virtual const std::vector<Real> & get_dzetadx() const override
749 4 : { libmesh_assert(!calculations_started || calculate_map);
750 10 : calculate_map = true; return dzetadx_map;}
751 :
752 :
753 : /**
754 : * \returns The dzeta/dy entry in the transformation
755 : * matrix from physical to local coordinates.
756 : */
757 10 : virtual const std::vector<Real> & get_dzetady() const override
758 4 : { libmesh_assert(!calculations_started || calculate_map);
759 10 : calculate_map = true; return dzetady_map;}
760 :
761 :
762 : /**
763 : * \returns The dzeta/dz entry in the transformation
764 : * matrix from physical to local coordinates.
765 : */
766 10 : virtual const std::vector<Real> & get_dzetadz() const override
767 4 : { libmesh_assert(!calculations_started || calculate_map);
768 10 : calculate_map = true; return dzetadz_map;}
769 :
770 :
771 : /**
772 : * \returns The multiplicative weight at each quadrature point.
773 : * This weight is used for certain infinite element weak
774 : * formulations, so that weighted Sobolev spaces are
775 : * used for the trial function space. This renders the
776 : * variational form easily computable.
777 : */
778 15 : virtual const std::vector<Real> & get_Sobolev_weight() const override
779 6 : { libmesh_assert(!calculations_started || calculate_phi);
780 15 : calculate_phi = true; return weight; }
781 :
782 :
783 : /**
784 : * \returns The first global derivative of the multiplicative
785 : * weight at each quadrature point. See \p get_Sobolev_weight()
786 : * for details. In case of \p FE initialized to all zero.
787 : */
788 10 : virtual const std::vector<RealGradient> & get_Sobolev_dweight() const override
789 4 : {libmesh_assert(!calculations_started || calculate_dphi);
790 10 : calculate_dphi = true; return dweight; }
791 :
792 :
793 :
794 : /**
795 : * \returns The tangent vectors for face integration.
796 : */
797 5 : virtual const std::vector<std::vector<Point>> & get_tangents() const override
798 2 : { libmesh_assert(!calculations_started || calculate_map);
799 5 : calculate_map = true; return tangents; }
800 :
801 : /**
802 : * \returns The outward pointing normal vectors for face integration.
803 : */
804 5 : virtual const std::vector<Point> & get_normals() const override
805 2 : { libmesh_assert(!calculations_started || calculate_map);
806 5 : calculate_map = true; return normals; }
807 :
808 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
809 : /**
810 : * \returns The curvatures for use in face integration.
811 : */
812 0 : virtual const std::vector<Real> & get_curvatures() const override
813 0 : { calculate_map = true; libmesh_not_implemented();}
814 : #endif
815 :
816 : /**
817 : * \returns The multiplicative weight (see \p get_Sobolev_weight())
818 : * but weighted with the radial coordinate square.
819 : *
820 : */
821 2597 : virtual const std::vector<Real> & get_Sobolev_weightxR_sq() const override
822 866 : { libmesh_assert(!calculations_started || calculate_phi_scaled);
823 2597 : calculate_phi_scaled = true; return weightxr_sq; }
824 :
825 :
826 : /**
827 : * \returns The first global derivative of the multiplicative
828 : * weight (see \p get_Sobolev_weight())
829 : * but weighted with the radial coordinate square.
830 : *
831 : */
832 2597 : virtual const std::vector<RealGradient> & get_Sobolev_dweightxR_sq() const override
833 866 : {libmesh_assert(!calculations_started || calculate_dphi_scaled);
834 2597 : calculate_dphi_scaled = true; return dweightxr_sq; }
835 :
836 :
837 : #ifdef LIBMESH_ENABLE_AMR
838 :
839 : /**
840 : * Computes the constraint matrix contributions (for
841 : * non-conforming adapted meshes) corresponding to
842 : * variable number \p var_number, adapted to infinite elements.
843 : */
844 : static void inf_compute_constraints (DofConstraints & constraints,
845 : DofMap & dof_map,
846 : const unsigned int variable_number,
847 : const Elem * child_elem);
848 : #endif
849 :
850 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
851 : static void inf_compute_node_constraints (NodeConstraints & constraints,
852 : const Elem * elem);
853 : #endif
854 :
855 :
856 : protected:
857 :
858 : // static members used by the workhorses
859 :
860 : /**
861 : * \returns The value of the \f$ i^{th} \f$ polynomial evaluated
862 : * at \p v. This method provides the approximation
863 : * in radial direction for the overall shape functions,
864 : * which is defined in \p InfFE::shape().
865 : * This method is allowed to be static, since it is independent
866 : * of dimension and base_family. It is templated, though,
867 : * w.r.t. to radial \p FEFamily.
868 : *
869 : * \returns The value of the \f$ i^{th} \f$ mapping shape function
870 : * in radial direction evaluated at \p v when T_radial ==
871 : * INFINITE_MAP. Currently, only one specific mapping shape is
872 : * used. Namely the one by Marques JMMC, Owen DRJ: Infinite
873 : * elements in quasi-static materially nonlinear problems, Computers
874 : * and Structures, 1984.
875 : */
876 : static Real eval(Real v,
877 : Order o_radial,
878 : unsigned int i);
879 :
880 : /**
881 : * \returns The value of the first derivative of the
882 : * \f$ i^{th} \f$ polynomial at coordinate \p v.
883 : * See \p eval for details.
884 : */
885 : static Real eval_deriv(Real v,
886 : Order o_radial,
887 : unsigned int i);
888 :
889 :
890 :
891 : // Non-static members used by the workhorses
892 :
893 : /**
894 : * Updates the protected member \p base_elem to the appropriate base element
895 : * for the given \p inf_elem.
896 : */
897 : void update_base_elem (const Elem * inf_elem);
898 :
899 : /**
900 : * Do not use this derived member in \p InfFE<Dim,T_radial,T_map>.
901 : */
902 0 : virtual void init_base_shape_functions(const std::vector<Point> &,
903 : const Elem *) override
904 0 : { libmesh_not_implemented(); }
905 :
906 : /**
907 : * Determine which values are to be calculated, for both the FE
908 : * itself and for the FEMap.
909 : */
910 : virtual void determine_calculations() override;
911 :
912 : /**
913 : * Some of the member data only depend on the radial part of the
914 : * infinite element. The parts that only change when the radial
915 : * order changes, are initialized here.
916 : */
917 : void init_radial_shape_functions(const Elem * inf_elem,
918 : const std::vector<Point> * radial_pts = nullptr);
919 :
920 : /**
921 : * Initialize all the data fields like \p weight, \p mode,
922 : * \p phi, \p dphidxi, \p dphideta, \p dphidzeta, etc.
923 : * for the current element. This method prepares the data
924 : * related to the base part, and some of the combined fields.
925 : */
926 : void init_shape_functions(const std::vector<Point> & radial_qp,
927 : const std::vector<Point> & base_qp,
928 : const Elem * inf_elem);
929 :
930 : /**
931 : * Initialize all the data fields like \p weight,
932 : * \p phi, etc for the side \p s.
933 : */
934 : void init_face_shape_functions (const std::vector<Point> &,
935 : const Elem * inf_side);
936 :
937 : /**
938 : * After having updated the jacobian and the transformation
939 : * from local to global coordinates in FEAbstract::compute_map(),
940 : * the first derivatives of the shape functions are
941 : * transformed to global coordinates, giving \p dphi,
942 : * \p dphidx/y/z, \p dphasedx/y/z, \p dweight. This method
943 : * should barely be re-defined in derived classes, but
944 : * still should be usable for children. Therefore, keep
945 : * it protected.
946 : */
947 : void compute_shape_functions(const Elem * inf_elem,
948 : const std::vector<Point> & base_qp,
949 : const std::vector<Point> & radial_qp);
950 :
951 : void compute_face_functions();
952 :
953 : /**
954 : * Use \p compute_shape_functions(const Elem*, const std::vector<Point> &, const std::vector<Point> &)
955 : * instead.
956 : */
957 0 : virtual void compute_shape_functions(const Elem *, const std::vector<Point> & ) override
958 : {
959 : //FIXME: it seems this function cannot be left out because
960 : // it is pure virtual in \p FEBase
961 0 : libmesh_not_implemented();
962 : }
963 :
964 : /**
965 : * Are we calculating scaled mapping functions?
966 : */
967 : mutable bool calculate_map_scaled;
968 :
969 : /**
970 : * Are we calculating scaled shape functions?
971 : */
972 : mutable bool calculate_phi_scaled;
973 :
974 : /**
975 : * Are we calculating scaled shape function gradients?
976 : */
977 : mutable bool calculate_dphi_scaled;
978 :
979 :
980 : /**
981 : * Are we calculating the positions of quadrature points?
982 : */
983 : mutable bool calculate_xyz;
984 :
985 :
986 : /**
987 : * Are we calculating the unscaled jacobian?
988 : * We avoid it if not requested explicitly; this has the worst stability.
989 : */
990 : mutable bool calculate_jxw;
991 :
992 : // Miscellaneous static members
993 :
994 : /**
995 : * Computes the indices in the base \p base_node and in radial
996 : * direction \p radial_node (either 0 or 1) associated to the
997 : * node \p outer_node_index of an infinite element of type
998 : * \p inf_elem_type.
999 : */
1000 : static void compute_node_indices (const ElemType inf_elem_type,
1001 : const unsigned int outer_node_index,
1002 : unsigned int & base_node,
1003 : unsigned int & radial_node);
1004 :
1005 : /**
1006 : * Does the same as \p compute_node_indices(), but stores
1007 : * the maps for the current element type. Provided the
1008 : * infinite element type changes seldom, this is probably
1009 : * faster than using \p compute_node_indices () alone.
1010 : * This is possible since the number of nodes is not likely
1011 : * to change.
1012 : */
1013 : static void compute_node_indices_fast (const ElemType inf_elem_type,
1014 : const unsigned int outer_node_index,
1015 : unsigned int & base_node,
1016 : unsigned int & radial_node);
1017 :
1018 : #ifdef LIBMESH_ENABLE_DEPRECATED
1019 : /*
1020 : * \deprecated Call the version of this function that takes an Elem * instead.
1021 : */
1022 : static void compute_shape_indices (const FEType & fet,
1023 : const ElemType inf_elem_type,
1024 : const unsigned int i,
1025 : unsigned int & base_shape,
1026 : unsigned int & radial_shape);
1027 : #endif // LIBMESH_ENABLE_DEPRECATED
1028 :
1029 : /**
1030 : * Computes the indices of shape functions in the base \p base_shape and
1031 : * in radial direction \p radial_shape (0 in the base, \f$ \ge 1 \f$ further
1032 : * out) associated to the shape with global index \p i of an infinite element
1033 : * \p inf_elem.
1034 : */
1035 : static void compute_shape_indices (const FEType & fet,
1036 : const Elem * inf_elem,
1037 : const unsigned int i,
1038 : unsigned int & base_shape,
1039 : unsigned int & radial_shape);
1040 :
1041 : /**
1042 : * Physical quadrature points.
1043 : * Usually, this is obtained from the \p FEMap class,
1044 : * but here FEMap does not know enough to compute it.
1045 : */
1046 : std::vector<Point> xyz;
1047 :
1048 : std::vector<Real> weightxr_sq;
1049 : /**
1050 : * the additional radial weight \f$ 1/{r^2} \f$ in local coordinates,
1051 : * over all quadrature points. The weight does not vary in base
1052 : * direction. However, for uniform access to the data fields from the
1053 : * outside, this data field is expanded to all quadrature points.
1054 : */
1055 : std::vector<Real> dweightdv;
1056 :
1057 : std::vector<RealGradient> dweightxr_sq;
1058 :
1059 : /**
1060 : * the radial decay \f$ 1/r \f$ in local coordinates. Needed when
1061 : * setting up the overall shape functions.
1062 : *
1063 : * \note It is this decay which ensures that the Sommerfeld
1064 : * radiation condition is satisfied in advance.
1065 : */
1066 : std::vector<Real> som;
1067 : /**
1068 : * the first local derivative of the radial decay \f$ 1/r \f$ in local
1069 : * coordinates. Needed when setting up the overall shape functions.
1070 : */
1071 : std::vector<Real> dsomdv;
1072 :
1073 : /**
1074 : * the radial approximation shapes in local coordinates
1075 : * Needed when setting up the overall shape functions.
1076 : */
1077 : std::vector<std::vector<Real>> mode;
1078 :
1079 : /**
1080 : * the first local derivative of the radial approximation shapes.
1081 : * Needed when setting up the overall shape functions.
1082 : */
1083 : std::vector<std::vector<Real>> dmodedv;
1084 :
1085 : // mapping of reference element to physical element
1086 : // These vectors usually belong to \p this->fe_map
1087 : // but for infinite elements, \p FEMap cannot
1088 : // compute them.
1089 : std::vector<Real> dxidx_map;
1090 : std::vector<Real> dxidy_map;
1091 : std::vector<Real> dxidz_map;
1092 : std::vector<Real> detadx_map;
1093 : std::vector<Real> detady_map;
1094 : std::vector<Real> detadz_map;
1095 : std::vector<Real> dzetadx_map;
1096 : std::vector<Real> dzetady_map;
1097 : std::vector<Real> dzetadz_map;
1098 :
1099 : // scaled mapping: Similar to the above
1100 : // vectors, but scaled by radial (infinite) coordinate.
1101 : std::vector<Real> dxidx_map_scaled;
1102 : std::vector<Real> dxidy_map_scaled;
1103 : std::vector<Real> dxidz_map_scaled;
1104 : std::vector<Real> detadx_map_scaled;
1105 : std::vector<Real> detady_map_scaled;
1106 : std::vector<Real> detadz_map_scaled;
1107 : std::vector<Real> dzetadx_map_scaled;
1108 : std::vector<Real> dzetady_map_scaled;
1109 : std::vector<Real> dzetadz_map_scaled;
1110 :
1111 :
1112 : // respectively weighted shape functions.
1113 : //FIXME: these names are correct only in 3D
1114 : // but avoid long and clumbsy names...
1115 : std::vector<std::vector<Real>> phixr;
1116 : std::vector<std::vector<RealGradient>> dphixr;
1117 : std::vector<std::vector<RealGradient>> dphixr_sq;
1118 :
1119 : std::vector<Real> JxWxdecay;
1120 : std::vector<Real> JxW;
1121 :
1122 : std::vector<Point> normals;
1123 : std::vector<std::vector<Point>> tangents;
1124 :
1125 : // numbering scheme maps
1126 :
1127 : /**
1128 : * The internal structure of the \p InfFE
1129 : * -- tensor product of base element times radial
1130 : * nodes -- has to be determined from the node numbering
1131 : * of the current infinite element. This vector
1132 : * maps the infinite \p Elem node number to the
1133 : * radial node (either 0 or 1).
1134 : */
1135 : std::vector<unsigned int> _radial_node_index;
1136 :
1137 : /**
1138 : * The internal structure of the \p InfFE
1139 : * -- tensor product of base element times radial
1140 : * nodes -- has to be determined from the node numbering
1141 : * of the current element. This vector
1142 : * maps the infinite \p Elem node number to the
1143 : * associated node in the base element.
1144 : */
1145 : std::vector<unsigned int> _base_node_index;
1146 :
1147 : /**
1148 : * The internal structure of the \p InfFE
1149 : * -- tensor product of base element shapes times radial
1150 : * shapes -- has to be determined from the dof numbering
1151 : * scheme of the current infinite element. This vector
1152 : * maps the infinite \p Elem dof index to the radial
1153 : * \p InfFE shape index (\p 0..radial_order+1 ).
1154 : */
1155 : std::vector<unsigned int> _radial_shape_index;
1156 :
1157 : /**
1158 : * The internal structure of the \p InfFE
1159 : * -- tensor product of base element shapes times radial
1160 : * shapes -- has to be determined from the dof numbering
1161 : * scheme of the current infinite element. This vector
1162 : * maps the infinite \p Elem dof index to the associated
1163 : * dof in the base \p FE.
1164 : */
1165 : std::vector<unsigned int> _base_shape_index;
1166 :
1167 : // some more protected members
1168 :
1169 : /**
1170 : * The number of total approximation shape functions for
1171 : * the current configuration
1172 : */
1173 : unsigned int _n_total_approx_sf;
1174 :
1175 : /**
1176 : * this vector contains the combined integration weights, so
1177 : * that \p FEAbstract::compute_map() can still be used
1178 : */
1179 : std::vector<Real> _total_qrule_weights;
1180 :
1181 : /**
1182 : * The quadrature rule for the base element associated
1183 : * with the current infinite element
1184 : */
1185 : std::unique_ptr<QBase> base_qrule;
1186 :
1187 : /**
1188 : * The quadrature rule for the base element associated
1189 : * with the current infinite element
1190 : */
1191 : std::unique_ptr<QBase> radial_qrule;
1192 :
1193 : /**
1194 : * The "base" (aka non-infinite) element associated with the current
1195 : * infinite element. We treat is as const since the InfFE should not
1196 : * have to modify the geometric Elem in order to do its calculations.
1197 : */
1198 : std::unique_ptr<const Elem> base_elem;
1199 :
1200 : /**
1201 : * Have a \p FE<Dim-1,T_base> handy for base approximation.
1202 : * Since this one is created using the \p FEBase::build() method,
1203 : * the \p InfFE class is not required to be templated w.r.t.
1204 : * to the base approximation shape.
1205 : */
1206 : std::unique_ptr<FEBase> base_fe;
1207 :
1208 : /**
1209 : * This \p FEType stores the characteristics for which the data
1210 : * structures \p phi, \p phi_map etc are currently initialized.
1211 : * This avoids re-initializing the radial part.
1212 : *
1213 : * \note Currently only \p order may change, both the FE families
1214 : * and \p base_order must remain constant.
1215 : */
1216 : FEType current_fe_type;
1217 :
1218 :
1219 : private:
1220 :
1221 : /**
1222 : * \returns \p false, currently not required.
1223 : */
1224 : virtual bool shapes_need_reinit() const override;
1225 :
1226 : /**
1227 : * When \p compute_node_indices_fast() is used, this static
1228 : * variable remembers the element type for which the
1229 : * static variables in \p compute_node_indices_fast()
1230 : * are currently set. Using a class member for the
1231 : * element type helps initializing it to a default value.
1232 : */
1233 : static ElemType _compute_node_indices_fast_current_elem_type;
1234 :
1235 :
1236 : #ifdef DEBUG
1237 :
1238 : /**
1239 : * static members that are used to issue warning messages only once.
1240 : */
1241 : static bool _warned_for_nodal_soln;
1242 : static bool _warned_for_shape;
1243 : static bool _warned_for_dshape;
1244 :
1245 : #endif
1246 :
1247 : /**
1248 : * Make all \p InfFE<Dim,T_radial,T_map> classes
1249 : * friends of each other, so that the protected
1250 : * \p eval() may be accessed.
1251 : */
1252 : template <unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
1253 : friend class InfFE;
1254 :
1255 : friend class InfFEMap;
1256 : };
1257 :
1258 :
1259 :
1260 : // InfFERadial class inline members
1261 : // FIXME: decay in 3D is a/r
1262 : // thus, this function fixes the mapping v<->r explicitly.
1263 : // This is consistent with the current InfFEMap, which, however, is
1264 : // written such that more general mappings can be implemented.
1265 : inline
1266 1896 : Real InfFERadial::decay(const unsigned int dim, const Real v)
1267 : {
1268 1896 : switch (dim)
1269 : {
1270 1896 : case 3:
1271 1896 : return (1.-v)/2.;
1272 :
1273 0 : case 2:
1274 : // according to P. Bettess "Infinite Elements",
1275 : // the 2D case is rather tricky:
1276 : // - the sqrt() requires special integration rules
1277 : // if not both trial- and test function are involved
1278 : // - the analytic solutions contain not only the sqrt, but
1279 : // also a polynomial with rather many terms, so convergence
1280 : // might be bad.
1281 0 : return sqrt((1.-v)/2.);
1282 :
1283 0 : case 1:
1284 0 : return 1.;
1285 :
1286 0 : default:
1287 0 : libmesh_error_msg("Invalid dim = " << dim);
1288 : }
1289 : }
1290 :
1291 : inline
1292 560 : Real InfFERadial::decay_deriv (const unsigned int dim, const Real v)
1293 : {
1294 560 : switch (dim)
1295 : {
1296 224 : case 3:
1297 224 : return -0.5;
1298 :
1299 0 : case 2:
1300 : // according to P. Bettess "Infinite Elements",
1301 : // the 2D case is rather tricky:
1302 : // - the sqrt() requires special integration rules
1303 : // if not both trial- and test function are involved
1304 : // - the analytic solutions contain not only the sqrt, but
1305 : // also a polynomial with rather many terms, so convergence
1306 : // might be bad.
1307 :
1308 0 : libmesh_assert (-1.-1.e-5 <= v && v < 1.);
1309 0 : return -0.25/sqrt(1.-v);
1310 :
1311 0 : case 1:
1312 0 : return 0.;
1313 :
1314 0 : default:
1315 0 : libmesh_error_msg("Invalid dim = " << dim);
1316 : }
1317 : }
1318 :
1319 : } // namespace libMesh
1320 :
1321 :
1322 : #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1323 :
1324 :
1325 : #endif // LIBMESH_INF_FE_H
|
