Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 80725 : static unsigned int n_dofs (const Order o_radial)
114 234839 : { 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 8540 : static unsigned int n_dofs_per_elem (const Order o_radial)
137 22806 : { 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 2208 : ~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 : #ifdef LIBMESH_ENABLE_DEPRECATED
386 : /*
387 : * \deprecated Call the version of this function that takes a
388 : * pointer-to-Elem instead.
389 : */
390 0 : static unsigned int n_shape_functions (const FEType & fet,
391 : const ElemType t)
392 0 : { return n_dofs(fet, t); }
393 :
394 : /**
395 : * \deprecated Call the version of this function that takes an Elem*
396 : * instead for consistency with other FEInterface::n_dofs() methods.
397 : */
398 : static unsigned int n_dofs(const FEType & fet,
399 : const ElemType inf_elem_type);
400 : #endif // LIBMESH_ENABLE_DEPRECATED
401 :
402 : /**
403 : * \returns The number of shape functions associated with this
404 : * infinite element. Currently, we have \p o_radial+1 modes in
405 : * radial direction, and \code FE<Dim-1,T>::n_dofs(...) \endcode
406 : * in the base.
407 : */
408 : static unsigned int n_dofs(const FEType & fet,
409 : const Elem * inf_elem);
410 :
411 : #ifdef LIBMESH_ENABLE_DEPRECATED
412 : /**
413 : * \returns The number of dofs at infinite element node \p n
414 : * (not dof!) for an element of type \p t and order \p o.
415 : */
416 : static unsigned int n_dofs_at_node(const FEType & fet,
417 : const ElemType inf_elem_type,
418 : const unsigned int n);
419 : #endif // LIBMESH_ENABLE_DEPRECATED
420 :
421 : /**
422 : * \returns The number of dofs at infinite element node \p n
423 : * (not dof!) for an element of type \p t and order \p o.
424 : */
425 : static unsigned int n_dofs_at_node(const FEType & fet,
426 : const Elem * inf_elem,
427 : const unsigned int n);
428 :
429 : #ifdef LIBMESH_ENABLE_DEPRECATED
430 : /**
431 : * \returns The number of dofs interior to the element,
432 : * not associated with any interior nodes.
433 : *
434 : * \deprecated Call the version of this function that takes an Elem*
435 : * instead for consistency with other FEInterface::n_dofs() methods.
436 : */
437 : static unsigned int n_dofs_per_elem(const FEType & fet,
438 : const ElemType inf_elem_type);
439 : #endif // LIBMESH_ENABLE_DEPRECATED
440 :
441 : /**
442 : * \returns The number of dofs interior to the element,
443 : * not associated with any interior nodes.
444 : */
445 : static unsigned int n_dofs_per_elem(const FEType & fet,
446 : const Elem * inf_elem);
447 :
448 : /**
449 : * \returns The continuity of the element.
450 : */
451 0 : virtual FEContinuity get_continuity() const override
452 0 : { return C_ZERO; } // FIXME - is this true??
453 :
454 : /**
455 : * \returns \p true if the element's higher order shape functions are
456 : * hierarchic
457 : */
458 0 : virtual bool is_hierarchic() const override
459 0 : { return false; } // FIXME - Inf FEs don't handle p elevation yet
460 :
461 : /**
462 : * Usually, this method would build the nodal soln from the
463 : * element soln. But infinite elements require additional
464 : * simulation-specific data to compute physically correct
465 : * results. Use \p compute_data() to compute results. For
466 : * compatibility an empty vector is returned.
467 : */
468 : static void nodal_soln(const FEType & fet,
469 : const Elem * elem,
470 : const std::vector<Number> & elem_soln,
471 : std::vector<Number> & nodal_soln);
472 :
473 :
474 0 : static Point map (const Elem * inf_elem,
475 : const Point & reference_point)
476 : {
477 : // libmesh_deprecated(); // soon
478 0 : return InfFEMap::map(Dim, inf_elem, reference_point);
479 : }
480 :
481 :
482 0 : static Point inverse_map (const Elem * elem,
483 : const Point & p,
484 : const Real tolerance = TOLERANCE,
485 : const bool secure = true)
486 : {
487 : // libmesh_deprecated(); // soon
488 0 : return InfFEMap::inverse_map(Dim, elem, p, tolerance, secure);
489 : }
490 :
491 :
492 0 : static void inverse_map (const Elem * elem,
493 : const std::vector<Point> & physical_points,
494 : std::vector<Point> & reference_points,
495 : const Real tolerance = TOLERANCE,
496 : const bool secure = true)
497 : {
498 : // libmesh_deprecated(); // soon
499 0 : return InfFEMap::inverse_map(Dim, elem, physical_points,
500 0 : reference_points, tolerance, secure);
501 : }
502 :
503 :
504 : // The workhorses of InfFE. These are often used during matrix assembly.
505 :
506 : /**
507 : * This is at the core of this class. Use this for each
508 : * new element in the mesh. Reinitializes all the physical
509 : * element-dependent data based on the current element
510 : * \p elem.
511 : * \note pts need to be in reference space coordinates, not physical ones.
512 : */
513 : virtual void reinit (const Elem * elem,
514 : const std::vector<Point> * const pts = nullptr,
515 : const std::vector<Real> * const weights = nullptr) override;
516 :
517 : /**
518 : * Reinitializes all the physical
519 : * element-dependent data based on the \p side of an infinite
520 : * element.
521 : */
522 : virtual void reinit (const Elem * inf_elem,
523 : const unsigned int s,
524 : const Real tolerance = TOLERANCE,
525 : const std::vector<Point> * const pts = nullptr,
526 : const std::vector<Real> * const weights = nullptr) override;
527 :
528 : /**
529 : * Not implemented yet. Reinitializes all the physical
530 : * element-dependent data based on the \p edge of an infinite
531 : * element.
532 : */
533 : virtual void edge_reinit (const Elem * elem,
534 : const unsigned int edge,
535 : const Real tolerance = TOLERANCE,
536 : const std::vector<Point> * const pts = nullptr,
537 : const std::vector<Real> * const weights = nullptr) override;
538 :
539 : /**
540 : * Computes the reference space quadrature points on the side of
541 : * an element based on the side quadrature points.
542 : */
543 0 : virtual void side_map (const Elem * /* elem */,
544 : const Elem * /* side */,
545 : const unsigned int /* s */,
546 : const std::vector<Point> & /* reference_side_points */,
547 : std::vector<Point> & /* reference_points */) override
548 : {
549 0 : libmesh_not_implemented();
550 : }
551 :
552 : /**
553 : * The use of quadrature rules with the \p InfFE class is somewhat
554 : * different from the approach of the \p FE class. While the
555 : * \p FE class requires an appropriately initialized quadrature
556 : * rule object, and simply uses it, the \p InfFE class requires only
557 : * the quadrature rule object of the current \p FE class.
558 : * From this \p QBase *, it determines the necessary data,
559 : * and builds two appropriate quadrature classes, one for radial,
560 : * and another for base integration, using the convenient
561 : * \p QBase::build() method.
562 : */
563 : virtual void attach_quadrature_rule (QBase * q) override;
564 :
565 : /**
566 : * \returns The number of shape functions associated with
567 : * this infinite element.
568 : */
569 0 : virtual unsigned int n_shape_functions () const override
570 0 : { return _n_total_approx_sf; }
571 :
572 : /**
573 : * \returns The total number of quadrature points. Call this
574 : * to get an upper bound for the \p for loop in your simulation
575 : * for matrix assembly of the current element.
576 : */
577 2592 : virtual unsigned int n_quadrature_points () const override
578 2592 : { libmesh_assert(radial_qrule); return this->_n_total_qp; }
579 :
580 : /**
581 : * \returns the \p xyz spatial locations of the quadrature
582 : * points on the element.
583 : */
584 15 : virtual const std::vector<Point> & get_xyz () const override
585 6 : { libmesh_assert(!calculations_started || calculate_xyz);
586 15 : calculate_xyz = true; return xyz; }
587 :
588 : /**
589 : * \returns the Jacobian times quadrature weight.
590 : * Due to the divergence with increasing radial distance, this quantity
591 : * is numerically unstable.
592 : * Thus, it is safer to use \p get_JxWxdecay_sq() instead!
593 : */
594 0 : virtual const std::vector<Real> & get_JxW () const override
595 : {
596 0 : libmesh_assert(!calculations_started || calculate_jxw);
597 0 : calculate_jxw = true;
598 0 : return this->JxW;
599 : }
600 :
601 : /**
602 : * \returns Jacobian times quadrature weight times square of the
603 : * decaying function \f$ decay= r^{-\frac{dim+1}{2}}\f$
604 : *
605 : * This function is the variant of \p get_JxW() for \p InfFE.
606 : * Since J diverges there, a respectize decay-function must be
607 : * applied to obtain well-defined quantities.
608 : */
609 2597 : virtual const std::vector<Real> & get_JxWxdecay_sq () const override
610 866 : { libmesh_assert(!calculations_started || calculate_map_scaled);
611 2597 : calculate_map_scaled = true; return this->JxWxdecay;}
612 :
613 :
614 : /**
615 : * \returns The shape function \p phi weighted by r/decay
616 : * where \f$ decay = r^{-\frac{dim+1}{2}} \f$
617 : *
618 : * To compensate for the decay function applied to the Jacobian (see \p get_JxWxdecay_sq),
619 : * the wave function \p phi should be divided by this function.
620 : *
621 : * The factor r must be compensated for by the Sobolev \p weight.
622 : * (i.e. by using \p get_Sobolev_weightxR_sq())
623 : **/
624 2597 : virtual const std::vector<std::vector<OutputShape>> & get_phi_over_decayxR () const override
625 866 : { libmesh_assert(!calculations_started || calculate_phi_scaled);
626 2597 : calculate_phi_scaled = true; return phixr; }
627 :
628 :
629 : /**
630 : * \returns the gradient of the shape function (see \p get_dphi()),
631 : * but in case of \p InfFE, weighted with r/decay.
632 : * See \p get_phi_over_decayxR() for details.
633 : */
634 2597 : virtual const std::vector<std::vector<OutputGradient>> & get_dphi_over_decayxR () const override
635 866 : { libmesh_assert(!calculations_started || calculate_dphi_scaled);
636 2597 : calculate_dphi_scaled = true; return dphixr; }
637 :
638 :
639 : /**
640 : * \returns the gradient of the shape function (see \p get_dphi()),
641 : * but in case of \p InfFE, weighted with 1/decay.
642 : *
643 : * In contrast to the shape function, its gradient stays finite
644 : * when divided by the decay function.
645 : */
646 0 : virtual const std::vector<std::vector<OutputGradient>> & get_dphi_over_decay () const override
647 0 : { libmesh_assert(!calculations_started || calculate_dphi_scaled);
648 0 : calculate_dphi_scaled = true; return dphixr_sq; }
649 :
650 :
651 : /**
652 : * \returns The element tangents in xi-direction at the quadrature
653 : * points.
654 : */
655 0 : virtual const std::vector<RealGradient> & get_dxyzdxi() const override
656 0 : { calculate_map = true; libmesh_not_implemented();}
657 :
658 :
659 : /**
660 : * \returns The element tangents in eta-direction at the quadrature
661 : * points.
662 : */
663 0 : virtual const std::vector<RealGradient> & get_dxyzdeta() const override
664 0 : { calculate_map = true; libmesh_not_implemented();}
665 :
666 :
667 : /**
668 : * \returns The element tangents in zeta-direction at the quadrature
669 : * points.
670 : */
671 0 : virtual const std::vector<RealGradient> & get_dxyzdzeta() const override
672 0 : { calculate_map = true; libmesh_not_implemented();}
673 :
674 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
675 : /**
676 : * \returns The second partial derivatives in xi.
677 : */
678 0 : virtual const std::vector<RealGradient> & get_d2xyzdxi2() const override
679 0 : { calculate_map = true; libmesh_not_implemented();}
680 : /**
681 : * \returns The second partial derivatives in eta.
682 : */
683 0 : virtual const std::vector<RealGradient> & get_d2xyzdeta2() const override
684 0 : { calculate_map = true; libmesh_not_implemented();}
685 : /**
686 : * \returns The second partial derivatives in zeta.
687 : */
688 0 : virtual const std::vector<RealGradient> & get_d2xyzdzeta2() const override
689 0 : { calculate_map = true; libmesh_not_implemented();}
690 : /**
691 : * \returns The second partial derivatives in xi-eta.
692 : */
693 0 : virtual const std::vector<RealGradient> & get_d2xyzdxideta() const override
694 0 : { calculate_map = true; libmesh_not_implemented();}
695 : /**
696 : * \returns The second partial derivatives in xi-zeta.
697 : */
698 0 : virtual const std::vector<RealGradient> & get_d2xyzdxidzeta() const override
699 0 : { calculate_map = true; libmesh_not_implemented();}
700 : /**
701 : * \returns The second partial derivatives in eta-zeta.
702 : */
703 0 : virtual const std::vector<RealGradient> & get_d2xyzdetadzeta() const override
704 0 : { calculate_map = true; libmesh_not_implemented();}
705 : #endif
706 :
707 : /**
708 : * \returns The dxi/dx entry in the transformation
709 : * matrix from physical to local coordinates.
710 : */
711 5 : virtual const std::vector<Real> & get_dxidx() const override
712 2 : {libmesh_assert(!calculations_started || calculate_map);
713 5 : calculate_map = true; return dxidx_map;}
714 :
715 :
716 : /**
717 : * \returns The dxi/dy entry in the transformation
718 : * matrix from physical to local coordinates.
719 : */
720 5 : virtual const std::vector<Real> & get_dxidy() const override
721 2 : {libmesh_assert(!calculations_started || calculate_map);
722 5 : calculate_map = true; return dxidy_map;}
723 :
724 :
725 : /**
726 : * \returns The dxi/dz entry in the transformation
727 : * matrix from physical to local coordinates.
728 : */
729 5 : virtual const std::vector<Real> & get_dxidz() const override
730 2 : { libmesh_assert(!calculations_started || calculate_map);
731 5 : calculate_map = true; return dxidz_map;}
732 :
733 :
734 : /**
735 : * \returns The deta/dx entry in the transformation
736 : * matrix from physical to local coordinates.
737 : */
738 5 : virtual const std::vector<Real> & get_detadx() const override
739 2 : { libmesh_assert(!calculations_started || calculate_map);
740 5 : calculate_map = true; return detadx_map;}
741 :
742 :
743 : /**
744 : * \returns The deta/dy entry in the transformation
745 : * matrix from physical to local coordinates.
746 : */
747 5 : virtual const std::vector<Real> & get_detady() const override
748 2 : { libmesh_assert(!calculations_started || calculate_map);
749 5 : calculate_map = true; return detady_map;}
750 :
751 :
752 : /**
753 : * \returns The deta/dx entry in the transformation
754 : * matrix from physical to local coordinates.
755 : */
756 5 : virtual const std::vector<Real> & get_detadz() const override
757 2 : { libmesh_assert(!calculations_started || calculate_map);
758 5 : calculate_map = true; return detadz_map;}
759 :
760 :
761 : /**
762 : * \returns The dzeta/dx entry in the transformation
763 : * matrix from physical to local coordinates.
764 : */
765 10 : virtual const std::vector<Real> & get_dzetadx() const override
766 4 : { libmesh_assert(!calculations_started || calculate_map);
767 10 : calculate_map = true; return dzetadx_map;}
768 :
769 :
770 : /**
771 : * \returns The dzeta/dy entry in the transformation
772 : * matrix from physical to local coordinates.
773 : */
774 10 : virtual const std::vector<Real> & get_dzetady() const override
775 4 : { libmesh_assert(!calculations_started || calculate_map);
776 10 : calculate_map = true; return dzetady_map;}
777 :
778 :
779 : /**
780 : * \returns The dzeta/dz entry in the transformation
781 : * matrix from physical to local coordinates.
782 : */
783 10 : virtual const std::vector<Real> & get_dzetadz() const override
784 4 : { libmesh_assert(!calculations_started || calculate_map);
785 10 : calculate_map = true; return dzetadz_map;}
786 :
787 :
788 : /**
789 : * \returns The multiplicative weight at each quadrature point.
790 : * This weight is used for certain infinite element weak
791 : * formulations, so that weighted Sobolev spaces are
792 : * used for the trial function space. This renders the
793 : * variational form easily computable.
794 : */
795 15 : virtual const std::vector<Real> & get_Sobolev_weight() const override
796 6 : { libmesh_assert(!calculations_started || calculate_phi);
797 15 : calculate_phi = true; return weight; }
798 :
799 :
800 : /**
801 : * \returns The first global derivative of the multiplicative
802 : * weight at each quadrature point. See \p get_Sobolev_weight()
803 : * for details. In case of \p FE initialized to all zero.
804 : */
805 10 : virtual const std::vector<RealGradient> & get_Sobolev_dweight() const override
806 4 : {libmesh_assert(!calculations_started || calculate_dphi);
807 10 : calculate_dphi = true; return dweight; }
808 :
809 :
810 :
811 : /**
812 : * \returns The tangent vectors for face integration.
813 : */
814 5 : virtual const std::vector<std::vector<Point>> & get_tangents() const override
815 2 : { libmesh_assert(!calculations_started || calculate_map);
816 5 : calculate_map = true; return tangents; }
817 :
818 : /**
819 : * \returns The outward pointing normal vectors for face integration.
820 : */
821 5 : virtual const std::vector<Point> & get_normals() const override
822 2 : { libmesh_assert(!calculations_started || calculate_map);
823 5 : calculate_map = true; return normals; }
824 :
825 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
826 : /**
827 : * \returns The curvatures for use in face integration.
828 : */
829 0 : virtual const std::vector<Real> & get_curvatures() const override
830 0 : { calculate_map = true; libmesh_not_implemented();}
831 : #endif
832 :
833 : /**
834 : * \returns The multiplicative weight (see \p get_Sobolev_weight())
835 : * but weighted with the radial coordinate square.
836 : *
837 : */
838 2597 : virtual const std::vector<Real> & get_Sobolev_weightxR_sq() const override
839 866 : { libmesh_assert(!calculations_started || calculate_phi_scaled);
840 2597 : calculate_phi_scaled = true; return weightxr_sq; }
841 :
842 :
843 : /**
844 : * \returns The first global derivative of the multiplicative
845 : * weight (see \p get_Sobolev_weight())
846 : * but weighted with the radial coordinate square.
847 : *
848 : */
849 2597 : virtual const std::vector<RealGradient> & get_Sobolev_dweightxR_sq() const override
850 866 : {libmesh_assert(!calculations_started || calculate_dphi_scaled);
851 2597 : calculate_dphi_scaled = true; return dweightxr_sq; }
852 :
853 :
854 : #ifdef LIBMESH_ENABLE_AMR
855 :
856 : /**
857 : * Computes the constraint matrix contributions (for
858 : * non-conforming adapted meshes) corresponding to
859 : * variable number \p var_number, adapted to infinite elements.
860 : */
861 : static void inf_compute_constraints (DofConstraints & constraints,
862 : DofMap & dof_map,
863 : const unsigned int variable_number,
864 : const Elem * child_elem);
865 : #endif
866 :
867 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
868 : static void inf_compute_node_constraints (NodeConstraints & constraints,
869 : const Elem * elem);
870 : #endif
871 :
872 :
873 : protected:
874 :
875 : // static members used by the workhorses
876 :
877 : /**
878 : * \returns The value of the \f$ i^{th} \f$ polynomial evaluated
879 : * at \p v. This method provides the approximation
880 : * in radial direction for the overall shape functions,
881 : * which is defined in \p InfFE::shape().
882 : * This method is allowed to be static, since it is independent
883 : * of dimension and base_family. It is templated, though,
884 : * w.r.t. to radial \p FEFamily.
885 : *
886 : * \returns The value of the \f$ i^{th} \f$ mapping shape function
887 : * in radial direction evaluated at \p v when T_radial ==
888 : * INFINITE_MAP. Currently, only one specific mapping shape is
889 : * used. Namely the one by Marques JMMC, Owen DRJ: Infinite
890 : * elements in quasi-static materially nonlinear problems, Computers
891 : * and Structures, 1984.
892 : */
893 : static Real eval(Real v,
894 : Order o_radial,
895 : unsigned int i);
896 :
897 : /**
898 : * \returns The value of the first derivative of the
899 : * \f$ i^{th} \f$ polynomial at coordinate \p v.
900 : * See \p eval for details.
901 : */
902 : static Real eval_deriv(Real v,
903 : Order o_radial,
904 : unsigned int i);
905 :
906 :
907 :
908 : // Non-static members used by the workhorses
909 :
910 : /**
911 : * Updates the protected member \p base_elem to the appropriate base element
912 : * for the given \p inf_elem.
913 : */
914 : void update_base_elem (const Elem * inf_elem);
915 :
916 : /**
917 : * Do not use this derived member in \p InfFE<Dim,T_radial,T_map>.
918 : */
919 0 : virtual void init_base_shape_functions(const std::vector<Point> &,
920 : const Elem *) override
921 0 : { libmesh_not_implemented(); }
922 :
923 : /**
924 : * Determine which values are to be calculated, for both the FE
925 : * itself and for the FEMap.
926 : */
927 : virtual void determine_calculations() override;
928 :
929 : /**
930 : * Some of the member data only depend on the radial part of the
931 : * infinite element. The parts that only change when the radial
932 : * order changes, are initialized here.
933 : */
934 : void init_radial_shape_functions(const Elem * inf_elem,
935 : const std::vector<Point> * radial_pts = nullptr);
936 :
937 : /**
938 : * Initialize all the data fields like \p weight, \p mode,
939 : * \p phi, \p dphidxi, \p dphideta, \p dphidzeta, etc.
940 : * for the current element. This method prepares the data
941 : * related to the base part, and some of the combined fields.
942 : */
943 : void init_shape_functions(const std::vector<Point> & radial_qp,
944 : const std::vector<Point> & base_qp,
945 : const Elem * inf_elem);
946 :
947 : /**
948 : * Initialize all the data fields like \p weight,
949 : * \p phi, etc for the side \p s.
950 : */
951 : void init_face_shape_functions (const std::vector<Point> &,
952 : const Elem * inf_side);
953 :
954 : /**
955 : * After having updated the jacobian and the transformation
956 : * from local to global coordinates in FEAbstract::compute_map(),
957 : * the first derivatives of the shape functions are
958 : * transformed to global coordinates, giving \p dphi,
959 : * \p dphidx/y/z, \p dphasedx/y/z, \p dweight. This method
960 : * should barely be re-defined in derived classes, but
961 : * still should be usable for children. Therefore, keep
962 : * it protected.
963 : */
964 : void compute_shape_functions(const Elem * inf_elem,
965 : const std::vector<Point> & base_qp,
966 : const std::vector<Point> & radial_qp);
967 :
968 : void compute_face_functions();
969 :
970 : /**
971 : * Use \p compute_shape_functions(const Elem*, const std::vector<Point> &, const std::vector<Point> &)
972 : * instead.
973 : */
974 0 : virtual void compute_shape_functions(const Elem *, const std::vector<Point> & ) override
975 : {
976 : //FIXME: it seems this function cannot be left out because
977 : // it is pure virtual in \p FEBase
978 0 : libmesh_not_implemented();
979 : }
980 :
981 : /**
982 : * Are we calculating scaled mapping functions?
983 : */
984 : mutable bool calculate_map_scaled;
985 :
986 : /**
987 : * Are we calculating scaled shape functions?
988 : */
989 : mutable bool calculate_phi_scaled;
990 :
991 : /**
992 : * Are we calculating scaled shape function gradients?
993 : */
994 : mutable bool calculate_dphi_scaled;
995 :
996 :
997 : /**
998 : * Are we calculating the positions of quadrature points?
999 : */
1000 : mutable bool calculate_xyz;
1001 :
1002 :
1003 : /**
1004 : * Are we calculating the unscaled jacobian?
1005 : * We avoid it if not requested explicitly; this has the worst stability.
1006 : */
1007 : mutable bool calculate_jxw;
1008 :
1009 : // Miscellaneous static members
1010 :
1011 : /**
1012 : * Computes the indices in the base \p base_node and in radial
1013 : * direction \p radial_node (either 0 or 1) associated to the
1014 : * node \p outer_node_index of an infinite element of type
1015 : * \p inf_elem_type.
1016 : */
1017 : static void compute_node_indices (const ElemType inf_elem_type,
1018 : const unsigned int outer_node_index,
1019 : unsigned int & base_node,
1020 : unsigned int & radial_node);
1021 :
1022 : /**
1023 : * Does the same as \p compute_node_indices(), but stores
1024 : * the maps for the current element type. Provided the
1025 : * infinite element type changes seldom, this is probably
1026 : * faster than using \p compute_node_indices () alone.
1027 : * This is possible since the number of nodes is not likely
1028 : * to change.
1029 : */
1030 : static void compute_node_indices_fast (const ElemType inf_elem_type,
1031 : const unsigned int outer_node_index,
1032 : unsigned int & base_node,
1033 : unsigned int & radial_node);
1034 :
1035 : #ifdef LIBMESH_ENABLE_DEPRECATED
1036 : /*
1037 : * \deprecated Call the version of this function that takes an Elem * instead.
1038 : */
1039 : static void compute_shape_indices (const FEType & fet,
1040 : const ElemType inf_elem_type,
1041 : const unsigned int i,
1042 : unsigned int & base_shape,
1043 : unsigned int & radial_shape);
1044 : #endif // LIBMESH_ENABLE_DEPRECATED
1045 :
1046 : /**
1047 : * Computes the indices of shape functions in the base \p base_shape and
1048 : * in radial direction \p radial_shape (0 in the base, \f$ \ge 1 \f$ further
1049 : * out) associated to the shape with global index \p i of an infinite element
1050 : * \p inf_elem.
1051 : */
1052 : static void compute_shape_indices (const FEType & fet,
1053 : const Elem * inf_elem,
1054 : const unsigned int i,
1055 : unsigned int & base_shape,
1056 : unsigned int & radial_shape);
1057 :
1058 : /**
1059 : * Physical quadrature points.
1060 : * Usually, this is obtained from the \p FEMap class,
1061 : * but here FEMap does not know enough to compute it.
1062 : */
1063 : std::vector<Point> xyz;
1064 :
1065 : std::vector<Real> weightxr_sq;
1066 : /**
1067 : * the additional radial weight \f$ 1/{r^2} \f$ in local coordinates,
1068 : * over all quadrature points. The weight does not vary in base
1069 : * direction. However, for uniform access to the data fields from the
1070 : * outside, this data field is expanded to all quadrature points.
1071 : */
1072 : std::vector<Real> dweightdv;
1073 :
1074 : std::vector<RealGradient> dweightxr_sq;
1075 :
1076 : /**
1077 : * the radial decay \f$ 1/r \f$ in local coordinates. Needed when
1078 : * setting up the overall shape functions.
1079 : *
1080 : * \note It is this decay which ensures that the Sommerfeld
1081 : * radiation condition is satisfied in advance.
1082 : */
1083 : std::vector<Real> som;
1084 : /**
1085 : * the first local derivative of the radial decay \f$ 1/r \f$ in local
1086 : * coordinates. Needed when setting up the overall shape functions.
1087 : */
1088 : std::vector<Real> dsomdv;
1089 :
1090 : /**
1091 : * the radial approximation shapes in local coordinates
1092 : * Needed when setting up the overall shape functions.
1093 : */
1094 : std::vector<std::vector<Real>> mode;
1095 :
1096 : /**
1097 : * the first local derivative of the radial approximation shapes.
1098 : * Needed when setting up the overall shape functions.
1099 : */
1100 : std::vector<std::vector<Real>> dmodedv;
1101 :
1102 : // mapping of reference element to physical element
1103 : // These vectors usually belong to \p this->fe_map
1104 : // but for infinite elements, \p FEMap cannot
1105 : // compute them.
1106 : std::vector<Real> dxidx_map;
1107 : std::vector<Real> dxidy_map;
1108 : std::vector<Real> dxidz_map;
1109 : std::vector<Real> detadx_map;
1110 : std::vector<Real> detady_map;
1111 : std::vector<Real> detadz_map;
1112 : std::vector<Real> dzetadx_map;
1113 : std::vector<Real> dzetady_map;
1114 : std::vector<Real> dzetadz_map;
1115 :
1116 : // scaled mapping: Similar to the above
1117 : // vectors, but scaled by radial (infinite) coordinate.
1118 : std::vector<Real> dxidx_map_scaled;
1119 : std::vector<Real> dxidy_map_scaled;
1120 : std::vector<Real> dxidz_map_scaled;
1121 : std::vector<Real> detadx_map_scaled;
1122 : std::vector<Real> detady_map_scaled;
1123 : std::vector<Real> detadz_map_scaled;
1124 : std::vector<Real> dzetadx_map_scaled;
1125 : std::vector<Real> dzetady_map_scaled;
1126 : std::vector<Real> dzetadz_map_scaled;
1127 :
1128 :
1129 : // respectively weighted shape functions.
1130 : //FIXME: these names are correct only in 3D
1131 : // but avoid long and clumbsy names...
1132 : std::vector<std::vector<Real>> phixr;
1133 : std::vector<std::vector<RealGradient>> dphixr;
1134 : std::vector<std::vector<RealGradient>> dphixr_sq;
1135 :
1136 : std::vector<Real> JxWxdecay;
1137 : std::vector<Real> JxW;
1138 :
1139 : std::vector<Point> normals;
1140 : std::vector<std::vector<Point>> tangents;
1141 :
1142 : // numbering scheme maps
1143 :
1144 : /**
1145 : * The internal structure of the \p InfFE
1146 : * -- tensor product of base element times radial
1147 : * nodes -- has to be determined from the node numbering
1148 : * of the current infinite element. This vector
1149 : * maps the infinite \p Elem node number to the
1150 : * radial node (either 0 or 1).
1151 : */
1152 : std::vector<unsigned int> _radial_node_index;
1153 :
1154 : /**
1155 : * The internal structure of the \p InfFE
1156 : * -- tensor product of base element times radial
1157 : * nodes -- has to be determined from the node numbering
1158 : * of the current element. This vector
1159 : * maps the infinite \p Elem node number to the
1160 : * associated node in the base element.
1161 : */
1162 : std::vector<unsigned int> _base_node_index;
1163 :
1164 : /**
1165 : * The internal structure of the \p InfFE
1166 : * -- tensor product of base element shapes times radial
1167 : * shapes -- has to be determined from the dof numbering
1168 : * scheme of the current infinite element. This vector
1169 : * maps the infinite \p Elem dof index to the radial
1170 : * \p InfFE shape index (\p 0..radial_order+1 ).
1171 : */
1172 : std::vector<unsigned int> _radial_shape_index;
1173 :
1174 : /**
1175 : * The internal structure of the \p InfFE
1176 : * -- tensor product of base element shapes times radial
1177 : * shapes -- has to be determined from the dof numbering
1178 : * scheme of the current infinite element. This vector
1179 : * maps the infinite \p Elem dof index to the associated
1180 : * dof in the base \p FE.
1181 : */
1182 : std::vector<unsigned int> _base_shape_index;
1183 :
1184 : // some more protected members
1185 :
1186 : /**
1187 : * The number of total approximation shape functions for
1188 : * the current configuration
1189 : */
1190 : unsigned int _n_total_approx_sf;
1191 :
1192 : /**
1193 : * this vector contains the combined integration weights, so
1194 : * that \p FEAbstract::compute_map() can still be used
1195 : */
1196 : std::vector<Real> _total_qrule_weights;
1197 :
1198 : /**
1199 : * The quadrature rule for the base element associated
1200 : * with the current infinite element
1201 : */
1202 : std::unique_ptr<QBase> base_qrule;
1203 :
1204 : /**
1205 : * The quadrature rule for the base element associated
1206 : * with the current infinite element
1207 : */
1208 : std::unique_ptr<QBase> radial_qrule;
1209 :
1210 : /**
1211 : * The "base" (aka non-infinite) element associated with the current
1212 : * infinite element. We treat is as const since the InfFE should not
1213 : * have to modify the geometric Elem in order to do its calculations.
1214 : */
1215 : std::unique_ptr<const Elem> base_elem;
1216 :
1217 : /**
1218 : * Have a \p FE<Dim-1,T_base> handy for base approximation.
1219 : * Since this one is created using the \p FEBase::build() method,
1220 : * the \p InfFE class is not required to be templated w.r.t.
1221 : * to the base approximation shape.
1222 : */
1223 : std::unique_ptr<FEBase> base_fe;
1224 :
1225 : /**
1226 : * This \p FEType stores the characteristics for which the data
1227 : * structures \p phi, \p phi_map etc are currently initialized.
1228 : * This avoids re-initializing the radial part.
1229 : *
1230 : * \note Currently only \p order may change, both the FE families
1231 : * and \p base_order must remain constant.
1232 : */
1233 : FEType current_fe_type;
1234 :
1235 :
1236 : private:
1237 :
1238 : /**
1239 : * \returns \p false, currently not required.
1240 : */
1241 : virtual bool shapes_need_reinit() const override;
1242 :
1243 : /**
1244 : * When \p compute_node_indices_fast() is used, this static
1245 : * variable remembers the element type for which the
1246 : * static variables in \p compute_node_indices_fast()
1247 : * are currently set. Using a class member for the
1248 : * element type helps initializing it to a default value.
1249 : */
1250 : static ElemType _compute_node_indices_fast_current_elem_type;
1251 :
1252 :
1253 : #ifdef DEBUG
1254 :
1255 : /**
1256 : * static members that are used to issue warning messages only once.
1257 : */
1258 : static bool _warned_for_nodal_soln;
1259 : static bool _warned_for_shape;
1260 : static bool _warned_for_dshape;
1261 :
1262 : #endif
1263 :
1264 : /**
1265 : * Make all \p InfFE<Dim,T_radial,T_map> classes
1266 : * friends of each other, so that the protected
1267 : * \p eval() may be accessed.
1268 : */
1269 : template <unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
1270 : friend class InfFE;
1271 :
1272 : friend class InfFEMap;
1273 : };
1274 :
1275 :
1276 :
1277 : // InfFERadial class inline members
1278 : // FIXME: decay in 3D is a/r
1279 : // thus, this function fixes the mapping v<->r explicitly.
1280 : // This is consistent with the current InfFEMap, which, however, is
1281 : // written such that more general mappings can be implemented.
1282 : inline
1283 1896 : Real InfFERadial::decay(const unsigned int dim, const Real v)
1284 : {
1285 1896 : switch (dim)
1286 : {
1287 1896 : case 3:
1288 1896 : return (1.-v)/2.;
1289 :
1290 0 : case 2:
1291 : // according to P. Bettess "Infinite Elements",
1292 : // the 2D case is rather tricky:
1293 : // - the sqrt() requires special integration rules
1294 : // if not both trial- and test function are involved
1295 : // - the analytic solutions contain not only the sqrt, but
1296 : // also a polynomial with rather many terms, so convergence
1297 : // might be bad.
1298 0 : return sqrt((1.-v)/2.);
1299 :
1300 0 : case 1:
1301 0 : return 1.;
1302 :
1303 0 : default:
1304 0 : libmesh_error_msg("Invalid dim = " << dim);
1305 : }
1306 : }
1307 :
1308 : inline
1309 560 : Real InfFERadial::decay_deriv (const unsigned int dim, const Real v)
1310 : {
1311 560 : switch (dim)
1312 : {
1313 224 : case 3:
1314 224 : return -0.5;
1315 :
1316 0 : case 2:
1317 : // according to P. Bettess "Infinite Elements",
1318 : // the 2D case is rather tricky:
1319 : // - the sqrt() requires special integration rules
1320 : // if not both trial- and test function are involved
1321 : // - the analytic solutions contain not only the sqrt, but
1322 : // also a polynomial with rather many terms, so convergence
1323 : // might be bad.
1324 :
1325 0 : libmesh_assert (-1.-1.e-5 <= v && v < 1.);
1326 0 : return -0.25/sqrt(1.-v);
1327 :
1328 0 : case 1:
1329 0 : return 0.;
1330 :
1331 0 : default:
1332 0 : libmesh_error_msg("Invalid dim = " << dim);
1333 : }
1334 : }
1335 :
1336 : } // namespace libMesh
1337 :
1338 :
1339 : #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1340 :
1341 :
1342 : #endif // LIBMESH_INF_FE_H
|
