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