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_FE_H
21 : #define LIBMESH_FE_H
22 :
23 : // Local includes
24 : #include "libmesh/fe_base.h"
25 : #include "libmesh/int_range.h"
26 : #include "libmesh/libmesh.h"
27 :
28 : // C++ includes
29 : #include <cmath>
30 : #include <cstddef>
31 :
32 : namespace libMesh
33 : {
34 :
35 : // forward declarations
36 : class DofConstraints;
37 : class DofMap;
38 : class QGauss;
39 :
40 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
41 :
42 : template <unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
43 : class InfFE;
44 :
45 : #endif
46 :
47 :
48 : /**
49 : * Most finite element types in libMesh are scalar-valued
50 : */
51 : template <FEFamily T>
52 : struct FEOutputType
53 : {
54 : typedef Real type;
55 : };
56 :
57 :
58 : /**
59 : * Specialize for non-scalar-valued elements
60 : */
61 : template<>
62 : struct FEOutputType<LAGRANGE_VEC>
63 : {
64 : typedef RealVectorValue type;
65 : };
66 :
67 : template<>
68 : struct FEOutputType<L2_LAGRANGE_VEC>
69 : {
70 : typedef RealVectorValue type;
71 : };
72 :
73 : template<>
74 : struct FEOutputType<HIERARCHIC_VEC>
75 : {
76 : typedef RealVectorValue type;
77 : };
78 :
79 : template<>
80 : struct FEOutputType<L2_HIERARCHIC_VEC>
81 : {
82 : typedef RealVectorValue type;
83 : };
84 :
85 : template<>
86 : struct FEOutputType<NEDELEC_ONE>
87 : {
88 : typedef RealVectorValue type;
89 : };
90 :
91 : template<>
92 : struct FEOutputType<MONOMIAL_VEC>
93 : {
94 : typedef RealVectorValue type;
95 : };
96 :
97 : template<>
98 : struct FEOutputType<RAVIART_THOMAS>
99 : {
100 : typedef RealVectorValue type;
101 : };
102 :
103 : template<>
104 : struct FEOutputType<L2_RAVIART_THOMAS>
105 : {
106 : typedef RealVectorValue type;
107 : };
108 :
109 :
110 : /**
111 : * A specific instantiation of the \p FEBase class. This
112 : * class is templated, and specific template instantiations
113 : * will result in different Finite Element families. Full specialization
114 : * of the template for specific dimensions(\p Dim) and families
115 : * (\p T) provide support for specific finite element types.
116 : * The use of templates allows for compile-time optimization,
117 : * however it requires that the specific finite element family
118 : * and dimension is also known at compile time. If this is
119 : * too restricting for your application you can use the
120 : * \p FEBase::build() member to create abstract (but still optimized)
121 : * finite elements.
122 : *
123 : * \author Benjamin S. Kirk
124 : * \date 2002-2007
125 : * \brief Template class which generates the different FE families and orders.
126 : */
127 : template <unsigned int Dim, FEFamily T>
128 : class FE : public FEGenericBase<typename FEOutputType<T>::type>
129 : {
130 : public:
131 :
132 : /**
133 : * Constructor.
134 : */
135 : explicit
136 : FE(const FEType & fet);
137 :
138 : typedef typename
139 : FEGenericBase<typename FEOutputType<T>::type>::OutputShape
140 : OutputShape;
141 :
142 : /**
143 : * \returns The value of the \f$ i^{th} \f$ shape function at
144 : * point \p p. This method allows you to specify the dimension,
145 : * element type, and order directly. This allows the method to
146 : * be static.
147 : *
148 : * On a p-refined element, \p o should be the total order of the element.
149 : */
150 : static OutputShape shape(const ElemType t,
151 : const Order o,
152 : const unsigned int i,
153 : const Point & p);
154 :
155 : /**
156 : * \returns The value of the \f$ i^{th} \f$ shape function at
157 : * point \p p. This method allows you to specify the dimension,
158 : * element type, and order directly. This allows the method to
159 : * be static.
160 : *
161 : * On a p-refined element, \p o should be the base order of the
162 : * element if \p add_p_level is left \p true, or can be the base
163 : * order of the element if \p add_p_level is set to \p false.
164 : */
165 : static OutputShape shape(const Elem * elem,
166 : const Order o,
167 : const unsigned int i,
168 : const Point & p,
169 : const bool add_p_level = true);
170 :
171 : /**
172 : * \returns The value of the \f$ i^{th} \f$ shape function at
173 : * point \p p. This method allows you to specify the dimension and
174 : * element type directly. The order is given by the FEType.
175 : * This allows the method to be static.
176 : *
177 : * On a p-refined element, \p o should be the base order of the
178 : * element if \p add_p_level is left \p true, or can be the base
179 : * order of the element if \p add_p_level is set to \p false.
180 : */
181 : static OutputShape shape(const FEType fet,
182 : const Elem * elem,
183 : const unsigned int i,
184 : const Point & p,
185 : const bool add_p_level = true);
186 :
187 : /**
188 : * Fills \p v with the values of the \f$ i^{th} \f$
189 : * shape function, evaluated at all points p. You must specify
190 : * element order directly. \p v should already be the appropriate
191 : * size.
192 : *
193 : * On a p-refined element, \p o should be the base order of the
194 : * element if \p add_p_level is left \p true, or can be the base
195 : * order of the element if \p add_p_level is set to \p false.
196 : */
197 : static void shapes(const Elem * elem,
198 : const Order o,
199 : const unsigned int i,
200 : const std::vector<Point> & p,
201 : std::vector<OutputShape> & v,
202 : const bool add_p_level = true);
203 :
204 :
205 : /**
206 : * Fills \p v[i][qp] with the values of the \f$ i^{th} \f$
207 : * shape functions, evaluated at all points in p. You must specify
208 : * element order directly. \p v should already be the appropriate
209 : * size.
210 : *
211 : * On a p-refined element, \p o should be the base order of the
212 : * element if \p add_p_level is left \p true, or can be the base
213 : * order of the element if \p add_p_level is set to \p false.
214 : */
215 : static void all_shapes(const Elem * elem,
216 : const Order o,
217 : const std::vector<Point> & p,
218 : std::vector<std::vector<OutputShape> > & v,
219 : const bool add_p_level = true);
220 :
221 : /**
222 : * \returns The \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
223 : * shape function at point \p p. This method allows you to
224 : * specify the dimension, element type, and order directly.
225 : *
226 : * On a p-refined element, \p o should be the total order of the element.
227 : */
228 : static OutputShape shape_deriv(const ElemType t,
229 : const Order o,
230 : const unsigned int i,
231 : const unsigned int j,
232 : const Point & p);
233 :
234 : /**
235 : * \returns The \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
236 : * shape function. You must specify element type, and order directly.
237 : *
238 : * On a p-refined element, \p o should be the base order of the
239 : * element if \p add_p_level is left \p true, or can be the base
240 : * order of the element if \p add_p_level is set to \p false.
241 : */
242 : static OutputShape shape_deriv(const Elem * elem,
243 : const Order o,
244 : const unsigned int i,
245 : const unsigned int j,
246 : const Point & p,
247 : const bool add_p_level = true);
248 :
249 : /**
250 : * \returns The \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
251 : * shape function. You must specify element type, and order (via
252 : * FEType) directly.
253 : *
254 : * On a p-refined element, \p o should be the base order of the
255 : * element if \p add_p_level is left \p true, or can be the base
256 : * order of the element if \p add_p_level is set to \p false.
257 : */
258 : static OutputShape shape_deriv(const FEType fet,
259 : const Elem * elem,
260 : const unsigned int i,
261 : const unsigned int j,
262 : const Point & p,
263 : const bool add_p_level = true);
264 :
265 : /**
266 : * Fills \p v with the \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
267 : * shape function, evaluated at all points p. You must specify
268 : * element order directly. \p v should already be the appropriate
269 : * size.
270 : *
271 : * On a p-refined element, \p o should be the base order of the
272 : * element if \p add_p_level is left \p true, or can be the base
273 : * order of the element if \p add_p_level is set to \p false.
274 : */
275 : static void shape_derivs(const Elem * elem,
276 : const Order o,
277 : const unsigned int i,
278 : const unsigned int j,
279 : const std::vector<Point> & p,
280 : std::vector<OutputShape> & v,
281 : const bool add_p_level = true);
282 :
283 : /**
284 : * Fills \p comps with dphidxi (and in higher dimensions, eta/zeta)
285 : * derivative component values for all shape functions, evaluated at
286 : * all points in p. You must specify element order directly.
287 : * Output component arrays in \p comps should already be the
288 : * appropriate size.
289 : *
290 : * On a p-refined element, \p o should be the base order of the
291 : * element if \p add_p_level is left \p true, or can be the base
292 : * order of the element if \p add_p_level is set to \p false.
293 : */
294 : static void all_shape_derivs(const Elem * elem,
295 : const Order o,
296 : const std::vector<Point> & p,
297 : std::vector<std::vector<OutputShape>> * comps[3],
298 : const bool add_p_level = true);
299 :
300 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
301 : /**
302 : * \returns The second \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
303 : * shape function at the point \p p.
304 : *
305 : * \note Cross-derivatives are indexed according to:
306 : * j = 0 ==> d^2 phi / dxi^2
307 : * j = 1 ==> d^2 phi / dxi deta
308 : * j = 2 ==> d^2 phi / deta^2
309 : * j = 3 ==> d^2 phi / dxi dzeta
310 : * j = 4 ==> d^2 phi / deta dzeta
311 : * j = 5 ==> d^2 phi / dzeta^2
312 : *
313 : * \note Computing second derivatives is not currently supported for
314 : * all element types: \f$ C^1 \f$ (Clough, Hermite and Subdivision),
315 : * Lagrange, Hierarchic, L2_Hierarchic, and Monomial are supported.
316 : * All other element types return an error when asked for second
317 : * derivatives.
318 : *
319 : * On a p-refined element, \p o should be the total order of the element.
320 : */
321 : static OutputShape shape_second_deriv(const ElemType t,
322 : const Order o,
323 : const unsigned int i,
324 : const unsigned int j,
325 : const Point & p);
326 :
327 : /**
328 : * \returns The second \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
329 : * shape function at the point \p p.
330 : *
331 : * \note Cross-derivatives are indexed according to:
332 : * j = 0 ==> d^2 phi / dxi^2
333 : * j = 1 ==> d^2 phi / dxi deta
334 : * j = 2 ==> d^2 phi / deta^2
335 : * j = 3 ==> d^2 phi / dxi dzeta
336 : * j = 4 ==> d^2 phi / deta dzeta
337 : * j = 5 ==> d^2 phi / dzeta^2
338 : *
339 : * \note Computing second derivatives is not currently supported for
340 : * all element types: \f$ C^1 \f$ (Clough, Hermite and Subdivision),
341 : * Lagrange, Hierarchic, L2_Hierarchic, and Monomial are supported.
342 : * All other element types return an error when asked for second
343 : * derivatives.
344 : *
345 : * On a p-refined element, \p o should be the base order of the
346 : * element if \p add_p_level is left \p true, or can be the base
347 : * order of the element if \p add_p_level is set to \p false.
348 : */
349 : static OutputShape shape_second_deriv(const Elem * elem,
350 : const Order o,
351 : const unsigned int i,
352 : const unsigned int j,
353 : const Point & p,
354 : const bool add_p_level = true);
355 :
356 : /**
357 : * \returns The second \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
358 : * shape function at the point \p p.
359 : *
360 : * \note Cross-derivatives are indexed according to:
361 : * j = 0 ==> d^2 phi / dxi^2
362 : * j = 1 ==> d^2 phi / dxi deta
363 : * j = 2 ==> d^2 phi / deta^2
364 : * j = 3 ==> d^2 phi / dxi dzeta
365 : * j = 4 ==> d^2 phi / deta dzeta
366 : * j = 5 ==> d^2 phi / dzeta^2
367 : *
368 : * \note Computing second derivatives is not currently supported for
369 : * all element types: \f$ C^1 \f$ (Clough, Hermite and Subdivision),
370 : * Lagrange, Hierarchic, L2_Hierarchic, and Monomial are supported.
371 : * All other element types return an error when asked for second
372 : * derivatives.
373 : *
374 : * On a p-refined element, \p o should be the total order of the element.
375 : */
376 : static OutputShape shape_second_deriv(const FEType fet,
377 : const Elem * elem,
378 : const unsigned int i,
379 : const unsigned int j,
380 : const Point & p,
381 : const bool add_p_level = true);
382 :
383 : #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
384 : /**
385 : * Build the nodal soln from the element soln.
386 : * This is the solution that will be plotted.
387 : *
388 : * On a p-refined element, \p o should be the base order of the element.
389 : */
390 : static void nodal_soln(const Elem * elem, const Order o,
391 : const std::vector<Number> & elem_soln,
392 : std::vector<Number> & nodal_soln,
393 : bool add_p_level = true,
394 : const unsigned vdim = 1);
395 :
396 : /**
397 : * Build the nodal soln on one side from the (full) element soln.
398 : * This is the solution that will be plotted on side-elements.
399 : *
400 : * On a p-refined element, \p o should be the base order of the element.
401 : */
402 : static void side_nodal_soln(const Elem * elem, const Order o,
403 : const unsigned int side,
404 : const std::vector<Number> & elem_soln,
405 : std::vector<Number> & nodal_soln_on_side,
406 : bool add_p_level = true,
407 : const unsigned vdim = 1);
408 :
409 : /**
410 : * \returns The number of shape functions associated with
411 : * this finite element.
412 : */
413 : virtual unsigned int n_shape_functions () const override;
414 :
415 : /**
416 : * \returns The number of shape functions associated with
417 : * a finite element of type \p t and approximation order \p o.
418 : *
419 : * On a p-refined element, \p o should be the total order of the element.
420 : *
421 : * This method does not support all finite element types; e.g. for an
422 : * arbitrary polygon or polyhedron type the number of shape
423 : * functions may depend on an individual element and not just its
424 : * type.
425 : */
426 0 : static unsigned int n_shape_functions (const ElemType t,
427 : const Order o)
428 1388 : { return FE<Dim,T>::n_dofs (t,o); }
429 :
430 : /**
431 : * \returns The number of shape functions associated with this
432 : * finite element.
433 : *
434 : * On a p-refined element, \p o should be the total order of the element.
435 : *
436 : * This method does not support all finite element types; e.g. for an
437 : * arbitrary polygon or polyhedron type the number of shape
438 : * functions may depend on an individual element and not just its
439 : * type.
440 : */
441 : static unsigned int n_dofs(const ElemType t,
442 : const Order o);
443 :
444 : /**
445 : * \returns The number of shape functions associated with this
446 : * finite element.
447 : *
448 : * On a p-refined element, \p o should be the total order of the element.
449 : *
450 : * \p e should only be a null pointer if using a FE family like
451 : * SCALAR that has degrees of freedom independent of any element.
452 : */
453 : static unsigned int n_dofs(const Elem * e,
454 : const Order o);
455 :
456 : /**
457 : * \returns The number of dofs at node \p n for a finite element
458 : * of type \p t and order \p o.
459 : *
460 : * On a p-refined element, \p o should be the total order of the element.
461 : *
462 : * This method does not support all finite element types; e.g. for an
463 : * arbitrary polygon or polyhedron type the meaning of a node index
464 : * \p n may depend on an individual element and not just its type.
465 : */
466 : static unsigned int n_dofs_at_node(const ElemType t,
467 : const Order o,
468 : const unsigned int n);
469 :
470 : /**
471 : * \returns The number of dofs at node \p n for a finite element
472 : * of type \p t and order \p o.
473 : *
474 : * On a p-refined element, \p o should be the total order of the element.
475 : */
476 : static unsigned int n_dofs_at_node(const Elem & e,
477 : const Order o,
478 : const unsigned int n);
479 :
480 : /**
481 : * \returns The number of dofs interior to the element,
482 : * not associated with any interior nodes.
483 : *
484 : * On a p-refined element, \p o should be the total order of the element.
485 : *
486 : * This method may not support all finite element types, e.g. higher
487 : * order polygons or polyhedra may differ from element to element.
488 : */
489 : static unsigned int n_dofs_per_elem(const ElemType t,
490 : const Order o);
491 :
492 : /**
493 : * \returns The number of dofs interior to the element,
494 : * not associated with any interior nodes.
495 : *
496 : * On a p-refined element, \p o should be the total order of the element.
497 : */
498 : static unsigned int n_dofs_per_elem(const Elem & e,
499 : const Order o);
500 :
501 : /**
502 : * \returns The continuity level of the finite element.
503 : */
504 : virtual FEContinuity get_continuity() const override;
505 :
506 : /**
507 : * \returns \p true if the finite element's higher order shape functions are
508 : * hierarchic
509 : */
510 : virtual bool is_hierarchic() const override;
511 :
512 : /**
513 : * Fills the vector di with the local degree of freedom indices
514 : * associated with side \p s of element \p elem
515 : *
516 : * On a p-refined element, \p o should be the base order of the element.
517 : */
518 : static void dofs_on_side(const Elem * const elem,
519 : const Order o,
520 : unsigned int s,
521 : std::vector<unsigned int> & di,
522 : bool add_p_level=true);
523 : /**
524 : * Fills the vector di with the local degree of freedom indices
525 : * associated with edge \p e of element \p elem
526 : *
527 : * On a p-refined element, \p o should be the base order of the element.
528 : */
529 : static void dofs_on_edge(const Elem * const elem,
530 : const Order o,
531 : unsigned int e,
532 : std::vector<unsigned int> & di,
533 : bool add_p_level=true);
534 :
535 0 : static Point inverse_map (const Elem * elem,
536 : const Point & p,
537 : const Real tolerance = TOLERANCE,
538 : const bool secure = true)
539 : {
540 : // libmesh_deprecated(); // soon
541 0 : return FEMap::inverse_map(Dim, elem, p, tolerance, secure, secure);
542 : }
543 :
544 0 : static void inverse_map (const Elem * elem,
545 : const std::vector<Point> & physical_points,
546 : std::vector<Point> & reference_points,
547 : const Real tolerance = TOLERANCE,
548 : const bool secure = true)
549 : {
550 : // libmesh_deprecated(); // soon
551 0 : FEMap::inverse_map(Dim, elem, physical_points, reference_points,
552 : tolerance, secure, secure);
553 0 : }
554 :
555 : /**
556 : * This is at the core of this class. Use this for each
557 : * new element in the mesh. Reinitializes all the physical
558 : * element-dependent data based on the current element
559 : * \p elem. By default the shape functions and associated
560 : * data are computed at the quadrature points specified
561 : * by the quadrature rule \p qrule, but may be any points
562 : * specified on the reference element specified in the optional
563 : * argument \p pts.
564 : */
565 : virtual void reinit (const Elem * elem,
566 : const std::vector<Point> * const pts = nullptr,
567 : const std::vector<Real> * const weights = nullptr) override;
568 :
569 : /**
570 : * This re-computes the dual shape function coefficients.
571 : * The dual shape coefficients are utilized when calculating dual shape functions.
572 : */
573 : virtual void reinit_dual_shape_coeffs (const Elem * elem,
574 : const std::vector<Point> & pts,
575 : const std::vector<Real> & JxW) override;
576 :
577 : /**
578 : * This computes the default dual shape function coefficients.
579 : * The dual shape coefficients are utilized when calculating dual shape functions.
580 : */
581 : virtual void reinit_default_dual_shape_coeffs (const Elem * elem) override;
582 :
583 : /**
584 : * Reinitializes all the physical element-dependent data based on
585 : * the \p side of \p face. The \p tolerance parameter is passed to
586 : * the involved call to \p inverse_map(). By default the shape
587 : * functions and associated data are computed at the quadrature
588 : * points specified by the quadrature rule \p qrule, but may be any
589 : * points specified on the reference \em side element specified in
590 : * the optional argument \p pts.
591 : */
592 : virtual void reinit (const Elem * elem,
593 : const unsigned int side,
594 : const Real tolerance = TOLERANCE,
595 : const std::vector<Point> * const pts = nullptr,
596 : const std::vector<Real> * const weights = nullptr) override;
597 :
598 : /**
599 : * Reinitializes all the physical element-dependent data based on
600 : * the \p edge. The \p tolerance parameter is passed to the
601 : * involved call to \p inverse_map(). By default the shape
602 : * functions and associated data are computed at the quadrature
603 : * points specified by the quadrature rule \p qrule, but may be any
604 : * points specified on the reference \em side element specified in
605 : * the optional argument \p pts.
606 : */
607 : virtual void edge_reinit (const Elem * elem,
608 : const unsigned int edge,
609 : const Real tolerance = TOLERANCE,
610 : const std::vector<Point> * const pts = nullptr,
611 : const std::vector<Real> * const weights = nullptr) override;
612 :
613 : /**
614 : * Computes the reference space quadrature points on the side of
615 : * an element based on the side quadrature points.
616 : */
617 : virtual void side_map (const Elem * elem,
618 : const Elem * side,
619 : const unsigned int s,
620 : const std::vector<Point> & reference_side_points,
621 : std::vector<Point> & reference_points) override;
622 :
623 : /**
624 : * Computes the reference space quadrature points on the side of
625 : * an element based on the edge quadrature points.
626 : */
627 : virtual void edge_map (const Elem * elem,
628 : const Elem * edge,
629 : const unsigned int e,
630 : const std::vector<Point> & reference_edge_points,
631 : std::vector<Point> & reference_points);
632 :
633 : /**
634 : * Provides the class with the quadrature rule, which provides the
635 : * locations (on a reference element) where the shape functions are
636 : * to be calculated.
637 : */
638 : virtual void attach_quadrature_rule (QBase * q) override;
639 :
640 : #ifdef LIBMESH_ENABLE_AMR
641 : /**
642 : * Computes the constraint matrix contributions (for
643 : * non-conforming adapted meshes) corresponding to
644 : * variable number \p var_number, using element-specific
645 : * optimizations if possible.
646 : */
647 : static void compute_constraints (DofConstraints & constraints,
648 : DofMap & dof_map,
649 : const unsigned int variable_number,
650 : const Elem * elem);
651 : #endif // #ifdef LIBMESH_ENABLE_AMR
652 :
653 : /**
654 : * \returns \p true when the shape functions (for
655 : * this \p FEFamily) depend on the particular
656 : * element, and therefore needs to be re-initialized
657 : * for each new element. \p false otherwise.
658 : */
659 : virtual bool shapes_need_reinit() const override;
660 :
661 0 : static Point map (const Elem * elem,
662 : const Point & reference_point)
663 : {
664 : // libmesh_deprecated(); // soon
665 0 : return FEMap::map(Dim, elem, reference_point);
666 : }
667 :
668 0 : static Point map_xi (const Elem * elem,
669 : const Point & reference_point)
670 : {
671 : // libmesh_deprecated(); // soon
672 0 : return FEMap::map_deriv(Dim, elem, 0, reference_point);
673 : }
674 :
675 0 : static Point map_eta (const Elem * elem,
676 : const Point & reference_point)
677 : {
678 : // libmesh_deprecated(); // soon
679 0 : return FEMap::map_deriv(Dim, elem, 1, reference_point);
680 : }
681 :
682 0 : static Point map_zeta (const Elem * elem,
683 : const Point & reference_point)
684 : {
685 : // libmesh_deprecated(); // soon
686 0 : return FEMap::map_deriv(Dim, elem, 2, reference_point);
687 : }
688 :
689 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
690 : /**
691 : * make InfFE classes friends, so that these may access
692 : * the private \p map, map_xyz methods
693 : */
694 : template <unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
695 : friend class InfFE;
696 : #endif
697 :
698 : protected:
699 :
700 : /**
701 : * Update the various member data fields \p phi,
702 : * \p dphidxi, \p dphideta, \p dphidzeta, etc.
703 : * for the current element. These data will be computed
704 : * at the points \p qp, which are generally (but need not be)
705 : * the quadrature points.
706 : */
707 : virtual void init_shape_functions(const std::vector<Point> & qp,
708 : const Elem * e);
709 :
710 : /**
711 : * A default implementation for all_shape_derivs
712 : */
713 : static void default_all_shape_derivs (const Elem * elem,
714 : const Order o,
715 : const std::vector<Point> & p,
716 : std::vector<std::vector<OutputShape>> * comps[3],
717 : const bool add_p_level = true);
718 :
719 :
720 : /**
721 : * Init \p dual_phi and potentially \p dual_dphi, \p dual_d2phi
722 : */
723 : void init_dual_shape_functions(unsigned int n_shapes, unsigned int n_qp);
724 :
725 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
726 :
727 : /**
728 : * Initialize the data fields for the base of an
729 : * an infinite element.
730 : */
731 : virtual void init_base_shape_functions(const std::vector<Point> & qp,
732 : const Elem * e) override;
733 :
734 : #endif
735 :
736 : /**
737 : * A default implementation for shapes
738 : */
739 1024326898 : static void default_shapes (const Elem * elem,
740 : const Order o,
741 : const unsigned int i,
742 : const std::vector<Point> & p,
743 : std::vector<OutputShape> & v,
744 : const bool add_p_level = true)
745 : {
746 89247114 : libmesh_assert_equal_to(p.size(), v.size());
747 5349271670 : for (auto vi : index_range(v))
748 4693649391 : v[vi] = FE<Dim,T>::shape (elem, o, i, p[vi], add_p_level);
749 1024326898 : }
750 :
751 : /**
752 : * A default implementation for all_shapes
753 : */
754 135826699 : static void default_all_shapes (const Elem * elem,
755 : const Order o,
756 : const std::vector<Point> & p,
757 : std::vector<std::vector<OutputShape>> & v,
758 : const bool add_p_level = true)
759 : {
760 1154887429 : for (auto i : index_range(v))
761 : {
762 88795200 : libmesh_assert_equal_to ( p.size(), v[i].size() );
763 1107015494 : FE<Dim,T>::shapes (elem, o, i, p, v[i], add_p_level);
764 : }
765 135826699 : }
766 :
767 : /**
768 : * A default implementation for shape_derivs
769 : */
770 1411502576 : static void default_shape_derivs (const Elem * elem,
771 : const Order o,
772 : const unsigned int i,
773 : const unsigned int j,
774 : const std::vector<Point> & p,
775 : std::vector<OutputShape> & v,
776 : const bool add_p_level = true)
777 : {
778 120287288 : libmesh_assert_equal_to(p.size(), v.size());
779 6918713954 : for (auto vi : index_range(v))
780 5951857093 : v[vi] = FE<Dim,T>::shape_deriv (elem, o, i, j, p[vi], add_p_level);
781 1411502576 : }
782 :
783 : /**
784 : * A default implementation for side_nodal_soln
785 : */
786 : static void default_side_nodal_soln(const Elem * elem, const Order o,
787 : const unsigned int side,
788 : const std::vector<Number> & elem_soln,
789 : std::vector<Number> & nodal_soln_on_side,
790 : bool add_p_level = true,
791 : const unsigned vdim = 1);
792 :
793 : /**
794 : * Vectors holding the node locations, edge and
795 : * face orientations of the last element we cached.
796 : */
797 : std::vector<Point> cached_nodes;
798 : std::vector<bool> cached_edges, cached_faces;
799 :
800 : /**
801 : * Repopulate the element cache with the node locations,
802 : * edge and face orientations of the element \p elem.
803 : */
804 : void cache(const Elem * elem);
805 :
806 : /**
807 : * Check if the node locations, edge and face orientations
808 : * held in the element cache match those of element \p elem.
809 : */
810 : bool matches_cache(const Elem * elem);
811 :
812 : /**
813 : * The last side and last edge we did a reinit on
814 : */
815 : ElemType last_side;
816 :
817 : ElemType last_edge;
818 : };
819 :
820 :
821 :
822 : /**
823 : * Clough-Tocher finite elements. Still templated on the dimension,
824 : * \p Dim.
825 : *
826 : * \author Roy Stogner
827 : * \date 2004
828 : */
829 : template <unsigned int Dim>
830 : class FEClough : public FE<Dim,CLOUGH>
831 : {
832 : public:
833 :
834 : /**
835 : * Constructor. Creates a hierarchic finite element
836 : * to be used in dimension \p Dim.
837 : */
838 : explicit
839 : FEClough(const FEType & fet) :
840 : FE<Dim,CLOUGH> (fet)
841 : {}
842 : };
843 :
844 :
845 :
846 : /**
847 : * Hermite finite elements. Still templated on the dimension,
848 : * \p Dim.
849 : *
850 : * \author Roy Stogner
851 : * \date 2005
852 : */
853 : template <unsigned int Dim>
854 : class FEHermite : public FE<Dim,HERMITE>
855 : {
856 : public:
857 :
858 : /**
859 : * Constructor. Creates a hierarchic finite element
860 : * to be used in dimension \p Dim.
861 : */
862 : explicit
863 : FEHermite(const FEType & fet) :
864 : FE<Dim,HERMITE> (fet)
865 : {}
866 :
867 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
868 : /**
869 : * 1D hermite functions on unit interval
870 : */
871 : static Real hermite_raw_shape_second_deriv(const unsigned int basis_num,
872 : const Real xi);
873 : #endif
874 : static Real hermite_raw_shape_deriv(const unsigned int basis_num,
875 : const Real xi);
876 : static Real hermite_raw_shape(const unsigned int basis_num,
877 : const Real xi);
878 : };
879 :
880 :
881 :
882 : /**
883 : * Subdivision finite elements.
884 : *
885 : * Template specialization prototypes are needed for calling from
886 : * inside FESubdivision::init_shape_functions
887 : */
888 : template <>
889 : Real FE<2,SUBDIVISION>::shape(const Elem * elem,
890 : const Order order,
891 : const unsigned int i,
892 : const Point & p,
893 : const bool add_p_level);
894 :
895 : template <>
896 : Real FE<2,SUBDIVISION>::shape_deriv(const Elem * elem,
897 : const Order order,
898 : const unsigned int i,
899 : const unsigned int j,
900 : const Point & p,
901 : const bool add_p_level);
902 :
903 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
904 : template <>
905 : Real FE<2,SUBDIVISION>::shape_second_deriv(const Elem * elem,
906 : const Order order,
907 : const unsigned int i,
908 : const unsigned int j,
909 : const Point & p,
910 : const bool add_p_level);
911 :
912 : #endif
913 :
914 : class FESubdivision : public FE<2,SUBDIVISION>
915 : {
916 : public:
917 :
918 : /**
919 : * Constructor. Creates a subdivision surface finite element.
920 : * Currently only supported for two-dimensional meshes in
921 : * three-dimensional space.
922 : */
923 : FESubdivision(const FEType & fet);
924 :
925 : /**
926 : * This is at the core of this class. Use this for each new
927 : * non-ghosted element in the mesh. Reinitializes all the physical
928 : * element-dependent data based on the current element
929 : * \p elem. By default the shape functions and associated
930 : * data are computed at the quadrature points specified
931 : * by the quadrature rule \p qrule, but may be any points
932 : * specified on the reference element specified in the optional
933 : * argument \p pts.
934 : */
935 : virtual void reinit (const Elem * elem,
936 : const std::vector<Point> * const pts = nullptr,
937 : const std::vector<Real> * const weights = nullptr) override;
938 :
939 : /**
940 : * This prevents some compilers being confused by partially
941 : * overriding this virtual function.
942 : */
943 0 : virtual void reinit (const Elem *,
944 : const unsigned int,
945 : const Real = TOLERANCE,
946 : const std::vector<Point> * const = nullptr,
947 : const std::vector<Real> * const = nullptr) override
948 0 : { libmesh_not_implemented(); }
949 :
950 : /**
951 : * Provides the class with the quadrature rule, which provides the
952 : * locations (on a reference element) where the shape functions are
953 : * to be calculated.
954 : */
955 : virtual void attach_quadrature_rule (QBase * q) override;
956 :
957 : /**
958 : * Update the various member data fields \p phi,
959 : * \p dphidxi, \p dphideta, \p dphidzeta, etc.
960 : * for the current element. These data will be computed
961 : * at the points \p qp, which are generally (but need not be)
962 : * the quadrature points.
963 : */
964 : virtual void init_shape_functions(const std::vector<Point> & qp,
965 : const Elem * elem) override;
966 :
967 : /**
968 : * \returns The value of the \f$ i^{th} \f$ of the 12 quartic
969 : * box splines interpolating a regular Loop subdivision
970 : * element, evaluated at the barycentric coordinates \p v,
971 : * \p w.
972 : */
973 : static Real regular_shape(const unsigned int i,
974 : const Real v,
975 : const Real w);
976 :
977 : /**
978 : * \returns The \f$ j^{th} \f$ derivative of the \f$ i^{th}
979 : * \f$ of the 12 quartic box splines interpolating a regular
980 : * Loop subdivision element, evaluated at the barycentric
981 : * coordinates \p v, \p w.
982 : */
983 : static Real regular_shape_deriv(const unsigned int i,
984 : const unsigned int j,
985 : const Real v,
986 : const Real w);
987 :
988 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
989 : /**
990 : * \returns The second \f$ j^{th} \f$ derivative of the
991 : * \f$ i^{th} \f$ of the 12 quartic box splines interpolating
992 : * a regular Loop subdivision element, evaluated at the
993 : * barycentric coordinates \p v, \p w.
994 : */
995 : static Real regular_shape_second_deriv(const unsigned int i,
996 : const unsigned int j,
997 : const Real v,
998 : const Real w);
999 :
1000 :
1001 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVE
1002 : /**
1003 : * Fills the vector \p weights with the weight coefficients
1004 : * of the Loop subdivision mask for evaluating the limit surface
1005 : * at a node explicitly. The size of \p weights will be
1006 : * 1 + \p valence, where \p valence is the number of neighbor
1007 : * nodes of the node where the limit surface is to be
1008 : * evaluated. The weight for the node itself is the first
1009 : * element of \p weights.
1010 : */
1011 : static void loop_subdivision_mask(std::vector<Real> & weights,
1012 : const unsigned int valence);
1013 :
1014 :
1015 : /**
1016 : * Builds the subdivision matrix \p A for the Loop scheme. The
1017 : * size depends on the element's \p valence.
1018 : */
1019 : static void init_subdivision_matrix(DenseMatrix<Real> & A,
1020 : unsigned int valence);
1021 : };
1022 :
1023 :
1024 :
1025 : /**
1026 : * Hierarchic finite elements. Still templated on the dimension,
1027 : * \p Dim.
1028 : *
1029 : * \author Benjamin S. Kirk
1030 : * \date 2002-2007
1031 : */
1032 : template <unsigned int Dim>
1033 : class FEHierarchic : public FE<Dim,HIERARCHIC>
1034 : {
1035 : public:
1036 :
1037 : /**
1038 : * Constructor. Creates a hierarchic finite element
1039 : * to be used in dimension \p Dim.
1040 : */
1041 : explicit
1042 : FEHierarchic(const FEType & fet) :
1043 : FE<Dim,HIERARCHIC> (fet)
1044 : {}
1045 : };
1046 :
1047 :
1048 :
1049 : /**
1050 : * Discontinuous Hierarchic finite elements. Still templated on the dimension,
1051 : * \p Dim.
1052 : *
1053 : * \author Truman E. Ellis
1054 : * \date 2011
1055 : */
1056 : template <unsigned int Dim>
1057 : class FEL2Hierarchic : public FE<Dim,L2_HIERARCHIC>
1058 : {
1059 : public:
1060 :
1061 : /**
1062 : * Constructor. Creates a hierarchic finite element
1063 : * to be used in dimension \p Dim.
1064 : */
1065 : explicit
1066 : FEL2Hierarchic(const FEType & fet) :
1067 : FE<Dim,L2_HIERARCHIC> (fet)
1068 : {}
1069 : };
1070 :
1071 :
1072 :
1073 : /**
1074 : * Lagrange finite elements. Still templated on the dimension,
1075 : * \p Dim.
1076 : *
1077 : * \author Benjamin S. Kirk
1078 : * \date 2002-2007
1079 : */
1080 : template <unsigned int Dim>
1081 : class FELagrange : public FE<Dim,LAGRANGE>
1082 : {
1083 : public:
1084 :
1085 : /**
1086 : * Constructor. Creates a Lagrange finite element
1087 : * to be used in dimension \p Dim.
1088 : */
1089 : explicit
1090 : FELagrange(const FEType & fet) :
1091 : FE<Dim,LAGRANGE> (fet)
1092 : {}
1093 : };
1094 :
1095 :
1096 : /**
1097 : * Discontinuous Lagrange finite elements.
1098 : */
1099 : template <unsigned int Dim>
1100 : class FEL2Lagrange : public FE<Dim,L2_LAGRANGE>
1101 : {
1102 : public:
1103 :
1104 : /**
1105 : * Constructor. Creates a discontinuous Lagrange finite element
1106 : * to be used in dimension \p Dim.
1107 : */
1108 : explicit
1109 : FEL2Lagrange(const FEType & fet) :
1110 : FE<Dim,L2_LAGRANGE> (fet)
1111 : {}
1112 : };
1113 :
1114 :
1115 : /**
1116 : * Monomial finite elements. Still templated on the dimension,
1117 : * \p Dim.
1118 : *
1119 : * \author Benjamin S. Kirk
1120 : * \date 2002-2007
1121 : */
1122 : template <unsigned int Dim>
1123 : class FEMonomial : public FE<Dim,MONOMIAL>
1124 : {
1125 : public:
1126 :
1127 : /**
1128 : * Constructor. Creates a monomial finite element
1129 : * to be used in dimension \p Dim.
1130 : */
1131 : explicit
1132 : FEMonomial(const FEType & fet) :
1133 : FE<Dim,MONOMIAL> (fet)
1134 : {}
1135 : };
1136 :
1137 :
1138 : /**
1139 : * The FEScalar class is used for working with SCALAR variables.
1140 : */
1141 : template <unsigned int Dim>
1142 : class FEScalar : public FE<Dim,SCALAR>
1143 : {
1144 : public:
1145 :
1146 : /**
1147 : * Constructor. Creates a SCALAR finite element
1148 : * which simply represents one or more
1149 : * extra DOFs coupled to all other DOFs in
1150 : * the system.
1151 : */
1152 : explicit
1153 16424 : FEScalar(const FEType & fet) :
1154 16424 : FE<Dim,SCALAR> (fet)
1155 464 : {}
1156 : };
1157 :
1158 :
1159 : /**
1160 : * XYZ finite elements. These require specialization
1161 : * because the shape functions are defined in terms of
1162 : * physical XYZ coordinates rather than local coordinates.
1163 : *
1164 : * \author Benjamin S. Kirk
1165 : * \date 2002-2007
1166 : */
1167 : template <unsigned int Dim>
1168 : class FEXYZ : public FE<Dim,XYZ>
1169 : {
1170 : public:
1171 :
1172 : /**
1173 : * Constructor. Creates a monomial finite element
1174 : * to be used in dimension \p Dim.
1175 : */
1176 : explicit
1177 990669 : FEXYZ(const FEType & fet) :
1178 990669 : FE<Dim,XYZ> (fet)
1179 34698 : {}
1180 :
1181 : /**
1182 : * Explicitly call base class method. This prevents some
1183 : * compilers being confused by partially overriding this virtual function.
1184 : * \note: pts need to be in reference space coordinates, not physical ones.
1185 : */
1186 1446281 : virtual void reinit (const Elem * elem,
1187 : const std::vector<Point> * const pts = nullptr,
1188 : const std::vector<Real> * const weights = nullptr) override
1189 1446281 : { FE<Dim,XYZ>::reinit (elem, pts, weights); }
1190 :
1191 : /**
1192 : * Reinitializes all the physical element-dependent data based on
1193 : * the \p side of \p face.
1194 : */
1195 : virtual void reinit (const Elem * elem,
1196 : const unsigned int side,
1197 : const Real tolerance = TOLERANCE,
1198 : const std::vector<Point> * const pts = nullptr,
1199 : const std::vector<Real> * const weights = nullptr) override;
1200 :
1201 :
1202 : protected:
1203 :
1204 : /**
1205 : * Update the various member data fields \p phi,
1206 : * \p dphidxi, \p dphideta, \p dphidzeta, etc.
1207 : * for the current element. These data will be computed
1208 : * at the points \p qp, which are generally (but need not be)
1209 : * the quadrature points.
1210 : */
1211 : virtual void init_shape_functions(const std::vector<Point> & qp,
1212 : const Elem * e) override;
1213 :
1214 : /**
1215 : * After having updated the jacobian and the transformation
1216 : * from local to global coordinates in \p FEAbstract::compute_map(),
1217 : * the first derivatives of the shape functions are
1218 : * transformed to global coordinates, giving \p dphi,
1219 : * \p dphidx, \p dphidy, and \p dphidz. This method
1220 : * should rarely be re-defined in derived classes, but
1221 : * still should be usable for children. Therefore, keep
1222 : * it protected.
1223 : */
1224 : virtual void compute_shape_functions(const Elem * elem, const std::vector<Point> & qp) override;
1225 :
1226 : /**
1227 : * Compute the map & shape functions for this face.
1228 : */
1229 : void compute_face_values (const Elem * elem,
1230 : const Elem * side,
1231 : const std::vector<Real> & weights);
1232 : };
1233 :
1234 :
1235 :
1236 : /**
1237 : * FELagrangeVec objects are used for working with vector-valued
1238 : * finite elements
1239 : *
1240 : * \author Paul T. Bauman
1241 : * \date 2013
1242 : */
1243 : template <unsigned int Dim>
1244 : class FELagrangeVec : public FE<Dim,LAGRANGE_VEC>
1245 : {
1246 : public:
1247 :
1248 : /**
1249 : * Constructor. Creates a vector Lagrange finite element
1250 : * to be used in dimension \p Dim.
1251 : */
1252 : explicit
1253 38890 : FELagrangeVec(const FEType & fet) :
1254 38890 : FE<Dim,LAGRANGE_VEC> (fet)
1255 3088 : {}
1256 : };
1257 :
1258 :
1259 : /**
1260 : * FEL2LagrangeVec objects are used for working with vector-valued
1261 : * finite elements
1262 : *
1263 : * \author Alexander Lindsay
1264 : * \date 2023
1265 : */
1266 : template <unsigned int Dim>
1267 : class FEL2LagrangeVec : public FE<Dim,L2_LAGRANGE_VEC>
1268 : {
1269 : public:
1270 :
1271 : /**
1272 : * Constructor. Creates a vector Lagrange finite element
1273 : * to be used in dimension \p Dim.
1274 : */
1275 : explicit
1276 4670 : FEL2LagrangeVec(const FEType & fet) :
1277 4670 : FE<Dim,L2_LAGRANGE_VEC> (fet)
1278 132 : {}
1279 : };
1280 :
1281 :
1282 :
1283 : /**
1284 : * FEHierarchicVec objects are used for working with vector-valued
1285 : * high-order finite elements
1286 : *
1287 : * \author Roy H. Stogner
1288 : * \date 2023
1289 : */
1290 : template <unsigned int Dim>
1291 : class FEHierarchicVec : public FE<Dim,HIERARCHIC_VEC>
1292 : {
1293 : public:
1294 :
1295 : /**
1296 : * Constructor. Creates a vector Hierarchic finite element
1297 : * to be used in dimension \p Dim.
1298 : */
1299 : explicit
1300 16572 : FEHierarchicVec(const FEType & fet) :
1301 16572 : FE<Dim,HIERARCHIC_VEC> (fet)
1302 1334 : {}
1303 : };
1304 :
1305 :
1306 : /**
1307 : * FEHierarchicVec objects are used for working with vector-valued
1308 : * high-order piecewise-continuous finite elements
1309 : *
1310 : * \author Roy H. Stogner
1311 : * \date 2023
1312 : */
1313 : template <unsigned int Dim>
1314 : class FEL2HierarchicVec : public FE<Dim,L2_HIERARCHIC_VEC>
1315 : {
1316 : public:
1317 :
1318 : /**
1319 : * Constructor. Creates a vector Hierarchic finite element
1320 : * to be used in dimension \p Dim.
1321 : */
1322 : explicit
1323 3550 : FEL2HierarchicVec(const FEType & fet) :
1324 3550 : FE<Dim,L2_HIERARCHIC_VEC> (fet)
1325 100 : {}
1326 : };
1327 :
1328 :
1329 : /**
1330 : * FENedelecOne objects are used for working with vector-valued
1331 : * Nedelec finite elements of the first kind.
1332 : *
1333 : * \author Paul T. Bauman
1334 : * \date 2013
1335 : */
1336 : template <unsigned int Dim>
1337 : class FENedelecOne : public FE<Dim,NEDELEC_ONE>
1338 : {
1339 : public:
1340 : /**
1341 : * Constructor. Creates a Nedelec finite element of the first kind
1342 : * to be used in dimension \p Dim.
1343 : */
1344 : explicit
1345 103851 : FENedelecOne(const FEType & fet) :
1346 103851 : FE<Dim,NEDELEC_ONE> (fet)
1347 4660 : {}
1348 : };
1349 :
1350 : /**
1351 : * FEMonomialVec objects are used for working with vector-valued
1352 : * discontinuous finite elements
1353 : *
1354 : * \author Alex D. Lindsay
1355 : * \date 2019
1356 : */
1357 : template <unsigned int Dim>
1358 : class FEMonomialVec : public FE<Dim,MONOMIAL_VEC>
1359 : {
1360 : public:
1361 :
1362 : /**
1363 : * Constructor. Creates a vector Monomial finite element
1364 : * to be used in dimension \p Dim.
1365 : */
1366 : explicit
1367 1470 : FEMonomialVec(const FEType & fet) :
1368 1470 : FE<Dim,MONOMIAL_VEC> (fet)
1369 42 : {}
1370 : };
1371 :
1372 : /**
1373 : * FERaviartThomas objects are used for working with vector-valued
1374 : * Raviart-Thomas finite elements.
1375 : *
1376 : * \author Nuno Nobre & Karthikeyan Chockalingam
1377 : * \date 2023
1378 : */
1379 : template <unsigned int Dim>
1380 : class FERaviartThomas : public FE<Dim,RAVIART_THOMAS>
1381 : {
1382 : public:
1383 :
1384 : /**
1385 : * Constructor. Creates a Raviart-Thomas finite element
1386 : * to be used in dimension \p Dim.
1387 : */
1388 : explicit
1389 438825 : FERaviartThomas(const FEType & fet) :
1390 438825 : FE<Dim,RAVIART_THOMAS> (fet)
1391 27522 : {}
1392 : };
1393 :
1394 : /**
1395 : * FEL2RaviartThomas objects are used for working with vector-valued
1396 : * discontinuous Raviart-Thomas finite elements, e.g. when constructing
1397 : * hybridized methods
1398 : *
1399 : * \author Alex D. Lindsay
1400 : * \date 2023
1401 : */
1402 : template <unsigned int Dim>
1403 : class FEL2RaviartThomas : public FE<Dim,L2_RAVIART_THOMAS>
1404 : {
1405 : public:
1406 :
1407 : /**
1408 : * Constructor. Creates a Raviart-Thomas finite element
1409 : * to be used in dimension \p Dim.
1410 : */
1411 : explicit
1412 8080 : FEL2RaviartThomas(const FEType & fet) :
1413 8080 : FE<Dim,L2_RAVIART_THOMAS> (fet)
1414 228 : {}
1415 : };
1416 :
1417 : /**
1418 : * Provide Typedefs for various element types.
1419 : */
1420 : namespace FiniteElements
1421 : {
1422 : /**
1423 : * Convenient definition for a 2D
1424 : * Clough-Tocher finite element.
1425 : */
1426 : typedef FEClough<2> FEClough2D;
1427 :
1428 : /**
1429 : * Convenient definition for a 1D
1430 : * Hierarchic finite element.
1431 : */
1432 : typedef FE<1,HIERARCHIC> FEHierarchic1D;
1433 :
1434 : /**
1435 : * Convenient definition for a 2D
1436 : * Hierarchic finite element.
1437 : */
1438 : typedef FE<2,HIERARCHIC> FEHierarchic2D;
1439 :
1440 : /**
1441 : * Convenient definition for a 3D
1442 : * Hierarchic finite element.
1443 : */
1444 : typedef FE<3,HIERARCHIC> FEHierarchic3D;
1445 :
1446 :
1447 : /**
1448 : * Convenient definition for a 1D
1449 : * Discontinuous Hierarchic finite element.
1450 : */
1451 : typedef FE<1,L2_HIERARCHIC> FEL2Hierarchic1D;
1452 :
1453 : /**
1454 : * Convenient definition for a 2D
1455 : * Discontinuous Hierarchic finite element.
1456 : */
1457 : typedef FE<2,L2_HIERARCHIC> FEL2Hierarchic2D;
1458 :
1459 : /**
1460 : * Convenient definition for a 3D
1461 : * Discontinuous Hierarchic finite element.
1462 : */
1463 : typedef FE<3,L2_HIERARCHIC> FEL2Hierarchic3D;
1464 :
1465 :
1466 : /**
1467 : * Convenient definition for a 1D
1468 : * Lagrange finite element.
1469 : */
1470 : typedef FE<1,LAGRANGE> FELagrange1D;
1471 :
1472 : /**
1473 : * Convenient definition for a 2D
1474 : * Lagrange finite element.
1475 : */
1476 : typedef FE<2,LAGRANGE> FELagrange2D;
1477 :
1478 : /**
1479 : * Convenient definition for a 3D
1480 : * Lagrange finite element.
1481 : */
1482 : typedef FE<3,LAGRANGE> FELagrange3D;
1483 :
1484 :
1485 : /**
1486 : * Convenient definition for a 1D
1487 : * Discontinuous Lagrange finite element.
1488 : */
1489 : typedef FE<1,L2_LAGRANGE> FEL2Lagrange1D;
1490 :
1491 : /**
1492 : * Convenient definition for a 2D
1493 : * Discontinuous Lagrange finite element.
1494 : */
1495 : typedef FE<2,L2_LAGRANGE> FEL2Lagrange2D;
1496 :
1497 : /**
1498 : * Convenient definition for a 3D
1499 : * Discontinuous Lagrange finite element.
1500 : */
1501 : typedef FE<3,L2_LAGRANGE> FEL2Lagrange3D;
1502 :
1503 :
1504 : /**
1505 : * Convenient definition for a 1D
1506 : * Monomial finite element.
1507 : */
1508 : typedef FE<1,MONOMIAL> FEMonomial1D;
1509 :
1510 : /**
1511 : * Convenient definition for a 2D
1512 : * Monomial finite element.
1513 : */
1514 : typedef FE<2,MONOMIAL> FEMonomial2D;
1515 :
1516 : /**
1517 : * Convenient definition for a 3D
1518 : * Monomial finite element.
1519 : */
1520 : typedef FE<3,MONOMIAL> FEMonomial3D;
1521 :
1522 : }
1523 :
1524 : /**
1525 : * Helper functions for finite differenced derivatives in cases where
1526 : * analytical calculations haven't been done yet.
1527 : */
1528 : template <typename OutputShape>
1529 : OutputShape fe_fdm_deriv(const Elem * elem,
1530 : const Order order,
1531 : const unsigned int i,
1532 : const unsigned int j,
1533 : const Point & p,
1534 : const bool add_p_level,
1535 : OutputShape(*shape_func)
1536 : (const Elem *, const Order,
1537 : const unsigned int, const Point &,
1538 : const bool));
1539 :
1540 : template <typename OutputShape>
1541 : OutputShape fe_fdm_deriv(const ElemType type,
1542 : const Order order,
1543 : const unsigned int i,
1544 : const unsigned int j,
1545 : const Point & p,
1546 : OutputShape(*shape_func)
1547 : (const ElemType, const Order,
1548 : const unsigned int, const Point &));
1549 :
1550 : template <typename OutputShape>
1551 : OutputShape fe_fdm_deriv(const ElemType type,
1552 : const Order order,
1553 : const Elem * elem,
1554 : const unsigned int i,
1555 : const unsigned int j,
1556 : const Point & p,
1557 : OutputShape(*shape_func)
1558 : (const ElemType type, const Order,
1559 : const Elem *, const unsigned int,
1560 : const Point &));
1561 :
1562 : /**
1563 : * The scaling that gives the \p i'th one-dimensional HIERARCHIC bubble function unit \f$H^1\f$
1564 : * seminorm on the reference interval, for \p i greater than one.
1565 : *
1566 : * The bubbles are \f$\xi^i - 1\f$ for even \p i and \f$\xi^i - \xi\f$ for odd \p i, up to this
1567 : * scaling. Integrating the square of their derivatives over \f$[-1,1]\f$ gives a seminorm of
1568 : * \f$\sqrt{2/(2i-1)}/(i-1)!\f$ in the even case and \f$(i-1)\sqrt{2/(2i-1)}/i!\f$ in the odd one, so
1569 : * the reciprocal of the seminorm cancels a factorial and what remains grows only as \f$\sqrt{i}\f$.
1570 : *
1571 : * The normalization matters because a shape function carrying \f$1/i!\f$ instead shrinks factorially
1572 : * with its order, and an operator assembled from such a basis inherits that spread on its diagonal.
1573 : * At order eight in two dimensions the smallest diagonal entry falls below the roundoff of the
1574 : * largest, which leaves the discretization numerically singular in double precision however it is
1575 : * solved.
1576 : *
1577 : * Only the bubbles are scaled. The two vertex functions are interpolatory, so their coefficients are
1578 : * values of the finite element solution at the vertices, and scaling them would change what a nodal
1579 : * boundary condition or a nodal output of a HIERARCHIC variable means.
1580 : */
1581 3782958227 : inline Real fe_hierarchic_bubble_scaling(const unsigned int i)
1582 : {
1583 312230987 : libmesh_assert_greater(i, 1);
1584 :
1585 : // An even bubble differentiates to xi^(i-1)/(i-1)!, whose square integrates to 2/(2i-1) over the
1586 : // interval. The linear term an odd bubble carries turns the i^2 of that calculation into (i-1)^2.
1587 3782958227 : const Real denominator = (i % 2) ? Real(i) - 1. : Real(i);
1588 :
1589 3782958227 : return std::sqrt((2. * Real(i) - 1.) / 2.) / denominator;
1590 : }
1591 :
1592 :
1593 : template <typename OutputShape>
1594 : OutputShape
1595 : fe_fdm_second_deriv(const Elem * elem,
1596 : const Order order,
1597 : const unsigned int i,
1598 : const unsigned int j,
1599 : const Point & p,
1600 : const bool add_p_level,
1601 : OutputShape(*deriv_func)
1602 : (const Elem *, const Order,
1603 : const unsigned int, const unsigned int,
1604 : const Point &, const bool));
1605 :
1606 : template <typename OutputShape>
1607 : OutputShape fe_fdm_second_deriv(const ElemType type,
1608 : const Order order,
1609 : const unsigned int i,
1610 : const unsigned int j,
1611 : const Point & p,
1612 : OutputShape(*deriv_func)
1613 : (const ElemType, const Order,
1614 : const unsigned int,
1615 : const unsigned int,
1616 : const Point &));
1617 :
1618 : template <typename OutputShape>
1619 : OutputShape fe_fdm_second_deriv(const ElemType type,
1620 : const Order order,
1621 : const Elem * elem,
1622 : const unsigned int i,
1623 : const unsigned int j,
1624 : const Point & p,
1625 : OutputShape(*deriv_func)
1626 : (const ElemType, const Order,
1627 : const Elem *,
1628 : const unsigned int,
1629 : const unsigned int,
1630 : const Point &));
1631 :
1632 : /**
1633 : * Helper functions for Lagrange-based basis functions.
1634 : */
1635 : void lagrange_nodal_soln(const Elem * elem,
1636 : const Order order,
1637 : const std::vector<Number> & elem_soln,
1638 : std::vector<Number> & nodal_soln,
1639 : bool add_p_level = true);
1640 :
1641 : /**
1642 : * Helper functions for Discontinuous-Pn type basis functions.
1643 : */
1644 : unsigned int monomial_n_dofs(const ElemType t, const Order o);
1645 :
1646 : unsigned int monomial_n_dofs(const Elem * e, const Order o);
1647 :
1648 : /**
1649 : * Helper functions for rational basis functions.
1650 : */
1651 : // shapes[i][j] is shape function phi_i at point p[j]
1652 : void rational_fe_weighted_shapes(const Elem * elem,
1653 : const FEType underlying_fe_type,
1654 : std::vector<std::vector<Real>> & shapes,
1655 : const std::vector<Point> & p,
1656 : const bool add_p_level);
1657 :
1658 : // shapes[i][q] is shape function phi_i at point p[q]
1659 : // derivs[j][i][q] is dphi_i/dxi_j at p[q]
1660 : void rational_fe_weighted_shapes_derivs(const Elem * elem,
1661 : const FEType fe_type,
1662 : std::vector<std::vector<Real>> & shapes,
1663 : std::vector<std::vector<std::vector<Real>>> & derivs,
1664 : const std::vector<Point> & p,
1665 : const bool add_p_level);
1666 :
1667 : Real rational_fe_shape(const Elem & elem,
1668 : const FEType underlying_fe_type,
1669 : const unsigned int i,
1670 : const Point & p,
1671 : const bool add_p_level);
1672 :
1673 : Real rational_fe_shape_deriv(const Elem & elem,
1674 : const FEType underlying_fe_type,
1675 : const unsigned int i,
1676 : const unsigned int j,
1677 : const Point & p,
1678 : const bool add_p_level);
1679 :
1680 : Real rational_fe_shape_second_deriv(const Elem & elem,
1681 : const FEType underlying_fe_type,
1682 : const unsigned int i,
1683 : const unsigned int j,
1684 : const Point & p,
1685 : const bool add_p_level);
1686 :
1687 : void rational_all_shapes (const Elem & elem,
1688 : const FEType underlying_fe_type,
1689 : const std::vector<Point> & p,
1690 : std::vector<std::vector<Real>> & v,
1691 : const bool add_p_level);
1692 :
1693 : template <typename OutputShape>
1694 : void rational_all_shape_derivs (const Elem & elem,
1695 : const FEType underlying_fe_type,
1696 : const std::vector<Point> & p,
1697 : std::vector<std::vector<OutputShape>> * comps[3],
1698 : const bool add_p_level);
1699 :
1700 : } // namespace libMesh
1701 :
1702 :
1703 : // Full specialization of all n_dofs type functions, for every
1704 : // dimension, with both original ElemType and new Elem signatures
1705 : #define LIBMESH_DEFAULT_NDOFS(MyType) \
1706 : template <> unsigned int FE<0,MyType>::n_dofs(const ElemType t, const Order o) { return MyType##_n_dofs(t, o); } \
1707 : template <> unsigned int FE<1,MyType>::n_dofs(const ElemType t, const Order o) { return MyType##_n_dofs(t, o); } \
1708 : template <> unsigned int FE<2,MyType>::n_dofs(const ElemType t, const Order o) { return MyType##_n_dofs(t, o); } \
1709 : template <> unsigned int FE<3,MyType>::n_dofs(const ElemType t, const Order o) { return MyType##_n_dofs(t, o); } \
1710 : \
1711 : template <> unsigned int FE<0,MyType>::n_dofs(const Elem * e, const Order o) { return MyType##_n_dofs(e, o); } \
1712 : template <> unsigned int FE<1,MyType>::n_dofs(const Elem * e, const Order o) { return MyType##_n_dofs(e, o); } \
1713 : template <> unsigned int FE<2,MyType>::n_dofs(const Elem * e, const Order o) { return MyType##_n_dofs(e, o); } \
1714 : template <> unsigned int FE<3,MyType>::n_dofs(const Elem * e, const Order o) { return MyType##_n_dofs(e, o); } \
1715 : \
1716 : template <> unsigned int FE<0,MyType>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(t, o, n); } \
1717 : template <> unsigned int FE<1,MyType>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(t, o, n); } \
1718 : template <> unsigned int FE<2,MyType>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(t, o, n); } \
1719 : template <> unsigned int FE<3,MyType>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(t, o, n); } \
1720 : \
1721 : template <> unsigned int FE<0,MyType>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(e, o, n); } \
1722 : template <> unsigned int FE<1,MyType>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(e, o, n); } \
1723 : template <> unsigned int FE<2,MyType>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(e, o, n); } \
1724 : template <> unsigned int FE<3,MyType>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(e, o, n); } \
1725 : \
1726 : template <> unsigned int FE<0,MyType>::n_dofs_per_elem(const ElemType t, const Order o) { return MyType##_n_dofs_per_elem(t, o); } \
1727 : template <> unsigned int FE<1,MyType>::n_dofs_per_elem(const ElemType t, const Order o) { return MyType##_n_dofs_per_elem(t, o); } \
1728 : template <> unsigned int FE<2,MyType>::n_dofs_per_elem(const ElemType t, const Order o) { return MyType##_n_dofs_per_elem(t, o); } \
1729 : template <> unsigned int FE<3,MyType>::n_dofs_per_elem(const ElemType t, const Order o) { return MyType##_n_dofs_per_elem(t, o); } \
1730 : \
1731 : template <> unsigned int FE<0,MyType>::n_dofs_per_elem(const Elem & e, const Order o) { return MyType##_n_dofs_per_elem(e, o); } \
1732 : template <> unsigned int FE<1,MyType>::n_dofs_per_elem(const Elem & e, const Order o) { return MyType##_n_dofs_per_elem(e, o); } \
1733 : template <> unsigned int FE<2,MyType>::n_dofs_per_elem(const Elem & e, const Order o) { return MyType##_n_dofs_per_elem(e, o); } \
1734 : template <> unsigned int FE<3,MyType>::n_dofs_per_elem(const Elem & e, const Order o) { return MyType##_n_dofs_per_elem(e, o); }
1735 :
1736 :
1737 : #define LIBMESH_DEFAULT_VEC_NDOFS(MyType) \
1738 : template <> unsigned int FE<0,MyType##_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,MyType>::n_dofs(t, o); } \
1739 : template <> unsigned int FE<1,MyType##_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,MyType>::n_dofs(t, o); } \
1740 : template <> unsigned int FE<2,MyType##_VEC>::n_dofs(const ElemType t, const Order o) { return 2*FE<2,MyType>::n_dofs(t, o); } \
1741 : template <> unsigned int FE<3,MyType##_VEC>::n_dofs(const ElemType t, const Order o) { return 3*FE<3,MyType>::n_dofs(t, o); } \
1742 : \
1743 : template <> unsigned int FE<0,MyType##_VEC>::n_dofs(const Elem * e, const Order o) { return FE<0,MyType>::n_dofs(e, o); } \
1744 : template <> unsigned int FE<1,MyType##_VEC>::n_dofs(const Elem * e, const Order o) { return FE<1,MyType>::n_dofs(e, o); } \
1745 : template <> unsigned int FE<2,MyType##_VEC>::n_dofs(const Elem * e, const Order o) { return 2*FE<2,MyType>::n_dofs(e, o); } \
1746 : template <> unsigned int FE<3,MyType##_VEC>::n_dofs(const Elem * e, const Order o) { return 3*FE<3,MyType>::n_dofs(e, o); } \
1747 : \
1748 : template <> unsigned int FE<0,MyType##_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<0,MyType>::n_dofs_at_node(t, o, n); } \
1749 : template <> unsigned int FE<1,MyType##_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<1,MyType>::n_dofs_at_node(t, o, n); } \
1750 : template <> unsigned int FE<2,MyType##_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 2*FE<2,MyType>::n_dofs_at_node(t, o, n); } \
1751 : template <> unsigned int FE<3,MyType##_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 3*FE<3,MyType>::n_dofs_at_node(t, o, n); } \
1752 : \
1753 : template <> unsigned int FE<0,MyType##_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<0,MyType>::n_dofs_at_node(e.type(), o, n); } \
1754 : template <> unsigned int FE<1,MyType##_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<1,MyType>::n_dofs_at_node(e.type(), o, n); } \
1755 : template <> unsigned int FE<2,MyType##_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return 2*FE<2,MyType>::n_dofs_at_node(e.type(), o, n); } \
1756 : template <> unsigned int FE<3,MyType##_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return 3*FE<3,MyType>::n_dofs_at_node(e.type(), o, n); } \
1757 : \
1758 : template <> unsigned int FE<0,MyType##_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<0,MyType>::n_dofs_per_elem(t, o); } \
1759 : template <> unsigned int FE<1,MyType##_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<1,MyType>::n_dofs_per_elem(t, o); } \
1760 : template <> unsigned int FE<2,MyType##_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return 2*FE<2,MyType>::n_dofs_per_elem(t, o); } \
1761 : template <> unsigned int FE<3,MyType##_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return 3*FE<3,MyType>::n_dofs_per_elem(t, o); } \
1762 : \
1763 : template <> unsigned int FE<0,MyType##_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<0,MyType>::n_dofs_per_elem(e.type(), o); } \
1764 : template <> unsigned int FE<1,MyType##_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<1,MyType>::n_dofs_per_elem(e.type(), o); } \
1765 : template <> unsigned int FE<2,MyType##_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return 2*FE<2,MyType>::n_dofs_per_elem(e.type(), o); } \
1766 : template <> unsigned int FE<3,MyType##_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return 3*FE<3,MyType>::n_dofs_per_elem(e.type(), o); }
1767 :
1768 :
1769 : #define LIBMESH_DEFAULT_VECTORIZED_FE(MyDim, MyType) \
1770 : template<> \
1771 : void FE<MyDim,MyType>::all_shapes \
1772 : (const Elem * elem, \
1773 : const Order o, \
1774 : const std::vector<Point> & p, \
1775 : std::vector<std::vector<OutputShape>> & v, \
1776 : const bool add_p_level) \
1777 : { \
1778 : FE<MyDim,MyType>::default_all_shapes \
1779 : (elem,o,p,v,add_p_level); \
1780 : } \
1781 : \
1782 : template<> \
1783 : void FE<MyDim,MyType>::shapes \
1784 : (const Elem * elem, \
1785 : const Order o, \
1786 : const unsigned int i, \
1787 : const std::vector<Point> & p, \
1788 : std::vector<OutputShape> & v, \
1789 : const bool add_p_level) \
1790 : { \
1791 : FE<MyDim,MyType>::default_shapes \
1792 : (elem,o,i,p,v,add_p_level); \
1793 : } \
1794 : \
1795 : template<> \
1796 : void FE<MyDim,MyType>::shape_derivs \
1797 : (const Elem * elem, \
1798 : const Order o, \
1799 : const unsigned int i, \
1800 : const unsigned int j, \
1801 : const std::vector<Point> & p, \
1802 : std::vector<OutputShape> & v, \
1803 : const bool add_p_level) \
1804 : { \
1805 : FE<MyDim,MyType>::default_shape_derivs \
1806 : (elem,o,i,j,p,v,add_p_level); \
1807 : } \
1808 : \
1809 : template<> \
1810 : void FE<MyDim,MyType>::all_shape_derivs \
1811 : (const Elem * elem, \
1812 : const Order o, \
1813 : const std::vector<Point> & p, \
1814 : std::vector<std::vector<OutputShape>> * comps[3], \
1815 : const bool add_p_level) \
1816 : { \
1817 : FE<MyDim,MyType>::default_all_shape_derivs \
1818 : (elem,o,p,comps,add_p_level); \
1819 : }
1820 :
1821 :
1822 : #endif // LIBMESH_FE_H
|
