MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
fe_base.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#ifndef MFEM_FE_BASE
13#define MFEM_FE_BASE
14
15#include "../intrules.hpp"
16#include "../geom.hpp"
17#include "../doftrans.hpp"
19
20#include <map>
21#include <memory>
22#include <unordered_map>
23#include <utility>
24
25namespace mfem
26{
27
28/// Possible basis types. Note that not all elements can use all BasisType(s).
30{
31public:
32 enum
33 {
34 Invalid = -1,
35 GaussLegendre = 0, ///< Open type
36 GaussLobatto = 1, ///< Closed type
37 Positive = 2, ///< Bernstein polynomials
38 OpenUniform = 3, ///< Nodes: x_i = (i+1)/(n+1), i=0,...,n-1
39 ClosedUniform = 4, ///< Nodes: x_i = i/(n-1), i=0,...,n-1
40 OpenHalfUniform = 5, ///< Nodes: x_i = (i+1/2)/n, i=0,...,n-1
41 Serendipity = 6, ///< Serendipity basis (squares / cubes)
42 ClosedGL = 7, ///< Closed GaussLegendre
43 IntegratedGLL = 8, ///< Integrated GLL indicator functions
44 NumBasisTypes = 9 /**< Keep track of maximum types to prevent
45 hard-coding */
46 };
47 /** @brief If the input does not represents a valid BasisType, abort with an
48 error; otherwise return the input. */
49 static int Check(int b_type)
50 {
51 MFEM_VERIFY(0 <= b_type && b_type < NumBasisTypes,
52 "unknown BasisType: " << b_type);
53 return b_type;
54 }
55 /** @brief If the input does not represents a valid nodal BasisType, abort
56 with an error; otherwise return the input. */
57 static int CheckNodal(int b_type)
58 {
59 MFEM_VERIFY(Check(b_type) != Positive && b_type != IntegratedGLL,
60 "invalid nodal BasisType: " << Name(b_type));
61 return b_type;
62 }
63 /** @brief Get the corresponding Quadrature1D constant, when that makes
64 sense; otherwise return Quadrature1D::Invalid. */
65 static int GetQuadrature1D(int b_type)
66 {
67 switch (b_type)
68 {
71 case Positive: return Quadrature1D::ClosedUniform; // <-----
78 }
80 }
81 /// Return the nodal BasisType corresponding to the Quadrature1D type.
82 static int GetNodalBasis(int qpt_type)
83 {
84 switch (qpt_type)
85 {
92 }
93 return Invalid;
94 }
95 /// Check and convert a BasisType constant to a string identifier.
96 static const char *Name(int b_type)
97 {
98 static const char *name[] =
99 {
100 "Gauss-Legendre", "Gauss-Lobatto", "Positive (Bernstein)",
101 "Open uniform", "Closed uniform", "Open half uniform",
102 "Serendipity", "Closed Gauss-Legendre",
103 "Integrated Gauss-Lobatto indicator"
104 };
105 return name[Check(b_type)];
106 }
107 /// Check and convert a BasisType constant to a char basis identifier.
108 static char GetChar(int b_type)
109 {
110 static const char ident[]
111 = { 'g', 'G', 'P', 'u', 'U', 'o', 'S', 'c', 'i' };
112 return ident[Check(b_type)];
113 }
114 /// Convert char basis identifier to a BasisType constant.
115 static int GetType(char b_ident)
116 {
117 switch (b_ident)
118 {
119 case 'g': return GaussLegendre;
120 case 'G': return GaussLobatto;
121 case 's': return GaussLobatto;
122 case 'P': return Positive;
123 case 'u': return OpenUniform;
124 case 'U': return ClosedUniform;
125 case 'o': return OpenHalfUniform;
126 case 'S': return Serendipity;
127 case 'c': return ClosedGL;
128 case 'i': return IntegratedGLL;
129 }
130 MFEM_ABORT("unknown BasisType identifier");
131 return -1;
132 }
133};
134
135/** @brief Structure representing the matrices/tensors needed to evaluate (in
136 reference space) the values, gradients, divergences, or curls of a
137 FiniteElement at the quadrature points of a given IntegrationRule. */
138/** Objects of this type are typically created and owned by the respective
139 FiniteElement object. */
141{
142public:
143 /// The FiniteElement that created and owns this object.
144 /** This pointer is not owned. */
145 const class FiniteElement *FE;
146
147 /** @brief IntegrationRule that defines the quadrature points at which the
148 basis functions of the #FE are evaluated. */
149 /** This pointer is not owned. */
151
152 /// Type of data stored in the arrays #B, #Bt, #G, and #Gt.
153 enum Mode
154 {
155 /** @brief Full multidimensional representation which does not use tensor
156 product structure. The ordering of the degrees of freedom is as
157 defined by #FE */
159
160 /** @brief Tensor product representation using 1D matrices/tensors with
161 dimensions using 1D number of quadrature points and degrees of
162 freedom. */
163 /** When representing a vector-valued FiniteElement, two DofToQuad objects
164 are used to describe the "closed" and "open" 1D basis functions. */
166
167 /** @brief Full multidimensional representation which does not use tensor
168 product structure. The ordering of the degrees of freedom is the
169 same as TENSOR, but the sizes of B and G are the same as FULL.*/
171 };
172
173 /// Describes the contents of the #B, #Bt, #G, and #Gt arrays, see #Mode.
175
176 /** @brief Number of degrees of freedom = number of basis functions. When
177 #mode is TENSOR, this is the 1D number. */
178 int ndof;
179
180 /** @brief Number of quadrature points. When #mode is TENSOR, this is the 1D
181 number. */
182 int nqpt;
183
184 /// Basis functions evaluated at quadrature points.
185 /** The storage layout is column-major with dimensions:
186 - #nqpt x #ndof, for scalar elements, or
187 - #nqpt x dim x #ndof, for vector elements,
188
189 where
190
191 - dim = dimension of the finite element reference space when #mode is
192 FULL, and dim = 1 when #mode is TENSOR. */
194
195 /// Transpose of #B.
196 /** The storage layout is column-major with dimensions:
197 - #ndof x #nqpt, for scalar elements, or
198 - #ndof x #nqpt x dim, for vector elements. */
200
201 /** @brief Gradients/divergences/curls of basis functions evaluated at
202 quadrature points. */
203 /** The storage layout is column-major with dimensions:
204 - #nqpt x dim x #ndof, for scalar elements, or
205 - #nqpt x #ndof, for H(div) vector elements, or
206 - #nqpt x cdim x #ndof, for H(curl) vector elements,
207
208 where
209
210 - dim = dimension of the finite element reference space when #mode is
211 FULL, and 1 when #mode is TENSOR,
212 - cdim = 1/1/3 in 1D/2D/3D, respectively, when #mode is FULL, and cdim =
213 1 when #mode is TENSOR. */
215
216 /// Transpose of #G.
217 /** The storage layout is column-major with dimensions:
218 - #ndof x #nqpt x dim, for scalar elements, or
219 - #ndof x #nqpt, for H(div) vector elements, or
220 - #ndof x #nqpt x cdim, for H(curl) vector elements. */
222};
223
224/// Describes the function space on each element
226{
227public:
228 enum
229 {
230 Pk, ///< Polynomials of order k
231 Qk, ///< Tensor products of polynomials of order k
232 rQk, ///< Refined tensor products of polynomials of order k
233 Uk ///< Rational polynomials of order k
234 };
235};
236
237class ElementTransformation;
238class Coefficient;
239class VectorCoefficient;
240class MatrixCoefficient;
241
242/// Abstract class for all finite elements.
244{
245protected:
246 int dim; ///< Dimension of reference space
247 int vdim; ///< Vector dimension of vector-valued basis functions
248 int cdim; ///< Dimension of curl for vector-valued basis functions
249 Geometry::Type geom_type; ///< Geometry::Type of the reference element
252 mutable
253 int dof, ///< Number of degrees of freedom
254 order; ///< Order/degree of the shape functions
255 mutable int orders[Geometry::MaxDim]; ///< Anisotropic orders
257#ifndef MFEM_THREAD_SAFE
258 mutable DenseMatrix vshape; // Dof x Dim
259#endif
260 /// Container for all DofToQuad objects created by the FiniteElement.
261 /** Multiple DofToQuad objects may be needed when different quadrature rules
262 or different DofToQuad::Mode are used. */
264
265public:
266 /// Enumeration for range_type and deriv_range_type
268
269 /** @brief Enumeration for MapType: defines how reference functions are
270 mapped to physical space.
271
272 A reference function $ \hat u(\hat x) $ can be mapped to a function
273 $ u(x) $ on a general physical element in following ways:
274 - $ x = T(\hat x) $ is the image of the reference point $ \hat x $
275 - $ J = J(\hat x) $ is the Jacobian matrix of the transformation T
276 - $ w = w(\hat x) = det(J) $ is the transformation weight factor for square J
277 - $ w = w(\hat x) = det(J^t J)^{1/2} $ is the transformation weight factor in general
278 */
280 {
281 UNKNOWN_MAP_TYPE = -1, /**< Used to distinguish an unset MapType variable
282 from the known values below. */
283 VALUE, /**< For scalar fields; preserves point values
284 $ u(x) = \hat u(\hat x) $ @anchor map_type_value */
285 INTEGRAL, /**< For scalar fields; preserves volume integrals
286 $ u(x) = (1/w) \hat u(\hat x) $ */
287 H_DIV, /**< For vector fields; preserves surface integrals of the
288 normal component $ u(x) = (J/w) \hat u(\hat x) $ */
289 H_CURL /**< For vector fields; preserves line integrals of the
290 tangential component
291 $ u(x) = J^{-t} \hat u(\hat x) $ (square J),
292 $ u(x) = J(J^t J)^{-1} \hat u(\hat x) $ (general J) */
293 };
294
295 /** @brief Enumeration for DerivType: defines which derivative method
296 is implemented.
297
298 Each FiniteElement class implements up to one type of derivative. The
299 value returned by GetDerivType() indicates which derivative method is
300 implemented.
301 */
303 {
304 NONE, ///< No derivatives implemented
305 GRAD, ///< Implements CalcDShape methods
306 DIV, ///< Implements CalcDivShape methods
307 CURL ///< Implements CalcCurlShape methods
308 };
309
310 /** @brief Construct FiniteElement with given
311 @param D Reference space dimension
312 @param G Geometry type (of type Geometry::Type)
313 @param Do Number of degrees of freedom in the FiniteElement
314 @param O Order/degree of the FiniteElement
315 @param F FunctionSpace type of the FiniteElement
316 */
317 FiniteElement(int D, Geometry::Type G, int Do, int O,
318 int F = FunctionSpace::Pk);
319
320 /// Returns the reference space dimension for the finite element.
321 int GetDim() const { return dim; }
322
323 /** @brief Returns the vector dimension for vector-valued finite elements,
324 which is also the dimension of the interpolation operation. */
325 int GetRangeDim() const { return vdim; }
326
327 /// Returns the dimension of the curl for vector-valued finite elements.
328 int GetCurlDim() const { return cdim; }
329
330 /// Returns the Geometry::Type of the reference element.
332
333 /// Returns the number of degrees of freedom in the finite element.
334 int GetDof() const { return dof; }
335
336 /** @brief Returns the order of the finite element. In the case of
337 anisotropic orders, returns the maximum order. */
338 int GetOrder() const { return order; }
339
340 /** @brief Returns true if the FiniteElement basis *may be using* different
341 orders/degrees in different spatial directions. */
342 bool HasAnisotropicOrders() const { return orders[0] != -1; }
343
344 /// Returns an array containing the anisotropic orders/degrees.
345 const int *GetAnisotropicOrders() const { return orders; }
346
347 /// Returns the type of FunctionSpace on the element.
348 int Space() const { return func_space; }
349
350 /// Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
351 int GetRangeType() const { return range_type; }
352
353 /** @brief Returns the FiniteElement::RangeType of the element derivative, either
354 SCALAR or VECTOR. */
355 int GetDerivRangeType() const { return deriv_range_type; }
356
357 /** @brief Returns the FiniteElement::MapType of the element describing how reference
358 functions are mapped to physical space, one of {VALUE, INTEGRAL
359 H_DIV, H_CURL}. */
360 int GetMapType() const { return map_type; }
361
362 /** @brief Returns the FiniteElement::DerivType of the element describing the
363 spatial derivative method implemented, one of {NONE, GRAD,
364 DIV, CURL}. */
365 int GetDerivType() const { return deriv_type; }
366
367 /** @brief Returns the FiniteElement::DerivType of the element describing how
368 reference function derivatives are mapped to physical space, one of {VALUE,
369 INTEGRAL, H_DIV, H_CURL}. */
370 int GetDerivMapType() const { return deriv_map_type; }
371
372 /** @brief Evaluate the values of all shape functions of a scalar finite
373 element in reference space at the given point @a ip. */
374 /** The size (#dof) of the result Vector @a shape must be set in advance. */
375 virtual void CalcShape(const IntegrationPoint &ip,
376 Vector &shape) const = 0;
377
378 /** @brief Evaluate the values of all shape functions of a scalar finite
379 element in physical space at the point described by @a Trans. */
380 /** The size (#dof) of the result Vector @a shape must be set in advance. */
381 void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const;
382
383 /** @brief Evaluate the gradients of all shape functions of a scalar finite
384 element in reference space at the given point @a ip. */
385 /** Each row of the result DenseMatrix @a dshape contains the derivatives of
386 one shape function. The size (#dof x #dim) of @a dshape must be set in
387 advance. */
388 virtual void CalcDShape(const IntegrationPoint &ip,
389 DenseMatrix &dshape) const = 0;
390
391 /** @brief Evaluate the gradients of all shape functions of a scalar finite
392 element in physical space at the point described by @a Trans. */
393 /** Each row of the result DenseMatrix @a dshape contains the derivatives of
394 one shape function. The size (#dof x SDim) of @a dshape must be set in
395 advance, where SDim >= #dim is the physical space dimension as described
396 by @a Trans. */
397 void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const;
398
399 /// Get a const reference to the nodes of the element
400 const IntegrationRule & GetNodes() const { return Nodes; }
401
402 /** @brief Evaluate the Hessians of all shape functions of a scalar finite
403 element in reference space at the given point @a ip. */
404 /** Each row of the result DenseMatrix @a Hessian contains upper triangular
405 part of the Hessian of one shape function.
406 The order in 2D is {u_xx, u_xy, u_yy}.
407 The size (#dof x (#dim (#dim+1)/2) of @a Hessian must be set in advance.*/
408 virtual void CalcHessian(const IntegrationPoint &ip,
409 DenseMatrix &Hessian) const;
410
411 /** @brief Evaluate the Hessian of all shape functions of a scalar finite
412 element in physical space at the given point @a ip. */
413 /** The size (#dof, #dim*(#dim+1)/2) of @a Hessian must be set in advance. */
415 DenseMatrix& Hessian) const;
416
417 /** @brief Evaluate the Laplacian of all shape functions of a scalar finite
418 element in physical space at the given point @a ip. */
419 /** The size (#dof) of @a Laplacian must be set in advance. */
421 Vector& Laplacian) const;
422
423 /** @brief Evaluate the Laplacian of all shape functions of a scalar finite
424 element in physical space at the given point @a ip. */
425 /** The size (#dof) of @a Laplacian must be set in advance. */
427 Vector& Laplacian) const;
428
429 /** @brief Evaluate the values of all shape functions of a *vector* finite
430 element in reference space at the given point @a ip. */
431 /** Each row of the result DenseMatrix @a shape contains the components of
432 one vector shape function. The size (#dof x #dim) of @a shape must be set
433 in advance. */
434 virtual void CalcVShape(const IntegrationPoint &ip,
435 DenseMatrix &shape) const;
436
437 /** @brief Evaluate the values of all shape functions of a *vector* finite
438 element in physical space at the point described by @a Trans. */
439 /** Each row of the result DenseMatrix @a shape contains the components of
440 one vector shape function. The size (#dof x SDim) of @a shape must be set
441 in advance, where SDim >= #dim is the physical space dimension as
442 described by @a Trans. */
443 virtual void CalcVShape(ElementTransformation &Trans,
444 DenseMatrix &shape) const;
445
446 /// Equivalent to the CalcVShape() method with the same arguments.
448 { CalcVShape(Trans, shape); }
449
450 /** @brief Evaluate the divergence of all shape functions of a *vector*
451 finite element in reference space at the given point @a ip. */
452 /** The size (#dof) of the result Vector @a divshape must be set in advance.
453 */
454 virtual void CalcDivShape(const IntegrationPoint &ip,
455 Vector &divshape) const;
456
457 /** @brief Evaluate the divergence of all shape functions of a *vector*
458 finite element in physical space at the point described by @a Trans. */
459 /** The size (#dof) of the result Vector @a divshape must be set in advance.
460 */
461 void CalcPhysDivShape(ElementTransformation &Trans, Vector &divshape) const;
462
463 /** @brief Evaluate the curl of all shape functions of a *vector* finite
464 element in reference space at the given point @a ip. */
465 /** Each row of the result DenseMatrix @a curl_shape contains the components
466 of the curl of one vector shape function. The size (#dof x CDim) of
467 @a curl_shape must be set in advance, where CDim = 3 for #dim = 3 and
468 CDim = 1 for #dim = 2. */
469 virtual void CalcCurlShape(const IntegrationPoint &ip,
470 DenseMatrix &curl_shape) const;
471
472 /** @brief Evaluate the curl of all shape functions of a *vector* finite
473 element in physical space at the point described by @a Trans. */
474 /** Each row of the result DenseMatrix @a curl_shape contains the components
475 of the curl of one vector shape function. The size (#dof x CDim) of
476 @a curl_shape must be set in advance, where CDim = 3 for #dim = 3 and
477 CDim = 1 for #dim = 2. */
478 virtual void CalcPhysCurlShape(ElementTransformation &Trans,
479 DenseMatrix &curl_shape) const;
480
481 /** @brief Get the dofs associated with the given @a face.
482 @a *dofs is set to an internal array of the local dofc on the
483 face, while *ndofs is set to the number of dofs on that face.
484 */
485 virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const;
486
487 /** @brief Return the local interpolation matrix @a I (Dof x Dof) where the
488 fine element is the image of the base geometry under the given
489 transformation. */
491 DenseMatrix &I) const;
492
493 /** @brief Return a local restriction matrix @a R (Dof x Dof) mapping fine
494 dofs to coarse dofs.
495
496 The fine element is the image of the base geometry under the given
497 transformation, @a Trans.
498
499 The assumption in this method is that a subset of the coarse dofs can be
500 expressed only in terms of the dofs of the given fine element.
501
502 Rows in @a R corresponding to coarse dofs that cannot be expressed in
503 terms of the fine dofs will be marked as invalid by setting the first
504 entry (column 0) in the row to infinity().
505
506 This method assumes that the dimensions of @a R are set before it is
507 called. */
508 virtual void GetLocalRestriction(ElementTransformation &Trans,
509 DenseMatrix &R) const;
510
511 /** @brief Return interpolation matrix, @a I, which maps dofs from a coarse
512 element, @a fe, to the fine dofs on @a this finite element. */
513 /** @a Trans represents the mapping from the reference element of @a this
514 element into a subset of the reference space of the element @a fe, thus
515 allowing the "coarse" FiniteElement to be different from the "fine"
516 FiniteElement as when h-refinement is combined with p-refinement or
517 p-derefinement. It is assumed that both finite elements use the same
518 FiniteElement::MapType. */
519 virtual void GetTransferMatrix(const FiniteElement &fe,
521 DenseMatrix &I) const;
522
523 /** @brief Given a coefficient and a transformation, compute its projection
524 (approximation) in the local finite dimensional space in terms
525 of the degrees of freedom. */
526 /** The approximation used to project is usually local interpolation of
527 degrees of freedom. The derived class could use other methods not
528 implemented yet, e.g. local L2 projection. */
529 virtual void Project(Coefficient &coeff,
530 ElementTransformation &Trans, Vector &dofs) const;
531
532 /** @brief Given a vector coefficient and a transformation, compute its
533 projection (approximation) in the local finite dimensional space
534 in terms of the degrees of freedom. (VectorFiniteElements) */
535 /** The approximation used to project is usually local interpolation of
536 degrees of freedom. The derived class could use other methods not
537 implemented yet, e.g. local L2 projection. */
538 virtual void Project(VectorCoefficient &vc,
539 ElementTransformation &Trans, Vector &dofs) const;
540
541 /** @brief Given a vector of values at the finite element nodes and a
542 transformation, compute its projection (approximation) in the local
543 finite dimensional space in terms of the degrees of freedom. Valid for
544 VectorFiniteElements. */
545 virtual void ProjectFromNodes(Vector &vc, ElementTransformation &Trans,
546 Vector &dofs) const;
547
548 /** @brief Given a matrix coefficient and a transformation, compute an
549 approximation ("projection") in the local finite dimensional space in
550 terms of the degrees of freedom. For VectorFiniteElements, the rows of
551 the coefficient are projected in the vector space. */
552 virtual void ProjectMatrixCoefficient(
553 MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const;
554
555 /** @brief Project a delta function centered on the given @a vertex in
556 the local finite dimensional space represented by the @a dofs. */
557 virtual void ProjectDelta(int vertex, Vector &dofs) const;
558
559 /** @brief Compute the embedding/projection matrix from the given
560 FiniteElement onto 'this' FiniteElement. The ElementTransformation is
561 included to support cases when the projection depends on it. */
562 virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
563 DenseMatrix &I) const;
564
565 /** @brief Compute the discrete gradient matrix from the given FiniteElement
566 onto 'this' FiniteElement. The ElementTransformation is included to
567 support cases when the matrix depends on it. */
568 virtual void ProjectGrad(const FiniteElement &fe,
570 DenseMatrix &grad) const;
571
572 /** @brief Compute the discrete curl matrix from the given FiniteElement onto
573 'this' FiniteElement. The ElementTransformation is included to support
574 cases when the matrix depends on it. */
575 virtual void ProjectCurl(const FiniteElement &fe,
577 DenseMatrix &curl) const;
578
579 /** @brief Compute the discrete divergence matrix from the given
580 FiniteElement onto 'this' FiniteElement. The ElementTransformation is
581 included to support cases when the matrix depends on it. */
582 virtual void ProjectDiv(const FiniteElement &fe,
584 DenseMatrix &div) const;
585
586 /** @brief Return a DofToQuad structure corresponding to the given
587 IntegrationRule using the given DofToQuad::Mode. */
588 /** See the documentation for DofToQuad for more details. */
589 virtual const DofToQuad &GetDofToQuad(const IntegrationRule &ir,
590 DofToQuad::Mode mode) const;
591
592
593 /** @brief Return the mapping from lexicographic face DOFs to lexicographic
594 element DOFs for the given local face @a face_id. */
595 /** Given the @a ith DOF (lexicographically ordered) on the face referenced
596 by @a face_id, face_map[i] gives the corresponding index of the DOF in
597 the element (also lexicographically ordered).
598
599 @note For L2 spaces, this is only well-defined for "closed" bases such as
600 the Gauss-Lobatto or Bernstein (positive) bases.
601
602 @warning GetFaceMap() is currently only implemented for tensor-product
603 (quadrilateral and hexahedral) elements. Its functionality may change
604 when simplex elements are supported in the future. */
605 virtual void GetFaceMap(const int face_id, Array<int> &face_map) const;
606
607 /** @brief Return a DoF transformation object for this particular type of
608 basis.
609 */
611 { return NULL; }
612
613 /// Deconstruct the FiniteElement
614 virtual ~FiniteElement();
615
616 /** @brief Return true if the BasisType of @a b_type is closed
617 (has Quadrature1D points on the boundary). */
618 static bool IsClosedType(int b_type)
619 {
620 const int q_type = BasisType::GetQuadrature1D(b_type);
621 return ((q_type != Quadrature1D::Invalid) &&
623 }
624
625 /** @brief Return true if the BasisType of @a b_type is open
626 (doesn't have Quadrature1D points on the boundary). */
627 static bool IsOpenType(int b_type)
628 {
629 const int q_type = BasisType::GetQuadrature1D(b_type);
630 return ((q_type != Quadrature1D::Invalid) &&
632 }
633
634 /** @brief Ensure that the BasisType of @a b_type is closed
635 (has Quadrature1D points on the boundary). */
636 static int VerifyClosed(int b_type)
637 {
638 MFEM_VERIFY(IsClosedType(b_type),
639 "invalid closed basis type: " << b_type);
640 return b_type;
641 }
642
643 /** @brief Ensure that the BasisType of @a b_type is open
644 (doesn't have Quadrature1D points on the boundary). */
645 static int VerifyOpen(int b_type)
646 {
647 MFEM_VERIFY(IsOpenType(b_type), "invalid open basis type: " << b_type);
648 return b_type;
649 }
650
651 /** @brief Ensure that the BasisType of @a b_type nodal
652 (satisfies the interpolation property). */
653 static int VerifyNodal(int b_type)
654 {
655 return BasisType::CheckNodal(b_type);
656 }
657};
658
659/** @brief Class for finite elements with basis functions
660 that return scalar values. */
662{
663protected:
665 {
666 MFEM_VERIFY(fe.GetRangeType() == SCALAR,
667 "'fe' must be a ScalarFiniteElement");
668 return static_cast<const ScalarFiniteElement &>(fe);
669 }
670
671public:
672 /** @brief Construct ScalarFiniteElement with given
673 @param D Reference space dimension
674 @param G Geometry type (of type Geometry::Type)
675 @param Do Number of degrees of freedom in the FiniteElement
676 @param O Order/degree of the FiniteElement
677 @param F FunctionSpace type of the FiniteElement
678 */
679 ScalarFiniteElement(int D, Geometry::Type G, int Do, int O,
680 int F = FunctionSpace::Pk)
681 : FiniteElement(D, G, Do, O, F)
683
684 /** @brief Set the FiniteElement::MapType of the element to either VALUE or
685 INTEGRAL. Also sets the FiniteElement::DerivType to GRAD if the
686 FiniteElement::MapType is VALUE. */
687 virtual void SetMapType(int M)
688 {
689 MFEM_VERIFY(M == VALUE || M == INTEGRAL, "unknown MapType");
690 map_type = M;
691 deriv_type = (M == VALUE) ? GRAD : NONE;
692 }
693
694 /** @brief Get the matrix @a I that defines nodal interpolation
695 @a between this element and the refined element @a fine_fe. */
697 DenseMatrix &I,
698 const ScalarFiniteElement &fine_fe) const;
699
700 /** @brief Get matrix @a I "Interpolation" defined through local
701 L2-projection in the space defined by the @a fine_fe. */
702 /** If the "fine" elements cannot represent all basis functions of the
703 "coarse" element, then boundary values from different sub-elements are
704 generally different. */
706 DenseMatrix &I,
707 const ScalarFiniteElement &fine_fe) const;
708
709 /** @brief Get restriction matrix @a R defined through local L2-projection
710 in the space defined by the @a coarse_fe. */
711 /** If the "fine" elements cannot represent all basis functions of the
712 "coarse" element, then boundary values from different sub-elements are
713 generally different. */
715 DenseMatrix &R,
716 const ScalarFiniteElement &coarse_fe) const;
717};
718
719/// Class for standard nodal finite elements.
721{
722private:
723 /// Create and cache the LEXICOGRAPHIC_FULL DofToQuad maps.
724 void CreateLexicographicFullMap(const IntegrationRule &ir) const;
725protected:
727 void ProjectCurl_2D(const FiniteElement &fe,
729 DenseMatrix &curl) const;
730
731public:
732 /** @brief Construct NodalFiniteElement with given
733 @param D Reference space dimension
734 @param G Geometry type (of type Geometry::Type)
735 @param Do Number of degrees of freedom in the FiniteElement
736 @param O Order/degree of the FiniteElement
737 @param F FunctionSpace type of the FiniteElement
738 */
739 NodalFiniteElement(int D, Geometry::Type G, int Do, int O,
740 int F = FunctionSpace::Pk)
741 : ScalarFiniteElement(D, G, Do, O, F) { }
742
743 const DofToQuad &GetDofToQuad(const IntegrationRule &ir,
744 DofToQuad::Mode mode) const override;
745
747 DenseMatrix &I) const override
748 { NodalLocalInterpolation(Trans, I, *this); }
749
751 DenseMatrix &R) const override;
752
755 DenseMatrix &I) const override
756 { CheckScalarFE(fe).NodalLocalInterpolation(Trans, I, *this); }
757
758 void Project(Coefficient &coeff,
759 ElementTransformation &Trans, Vector &dofs) const override;
760
761 void Project(VectorCoefficient &vc,
762 ElementTransformation &Trans, Vector &dofs) const override;
763
764 // (mc.height x mc.width) @ DOFs -> (Dof x mc.width x mc.height) in dofs
766 MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const override;
767
768 void Project(const FiniteElement &fe, ElementTransformation &Trans,
769 DenseMatrix &I) const override;
770
771 void ProjectGrad(const FiniteElement &fe,
773 DenseMatrix &grad) const override;
774
775 void ProjectDiv(const FiniteElement &fe,
777 DenseMatrix &div) const override;
778
779 /** @brief Get an Array<int> that maps lexicographically ordered indices to
780 the indices of the respective nodes/dofs/basis functions.
781
782 Lexicographic ordering of nodes is defined in terms of reference-space
783 coordinates (x,y,z). Lexicographically ordered nodes are listed first in
784 order of increasing x-coordinate, and then in order of increasing
785 y-coordinate, and finally in order of increasing z-coordinate.
786
787 For example, the six nodes of a quadratic triangle are lexicographically
788 ordered as follows:
789
790 5
791 |\
792 3 4
793 | \
794 0-1-2
795
796 The resulting array may be empty if the DOFs are already ordered
797 lexicographically, or if the finite element does not support creating
798 this permutation. The array returned is the same as the array given by
799 TensorBasisElement::GetDofMap, but it is also available for non-tensor
800 elements. */
802
803 /// Given a lexicographically ordered Vector @a dofs, containing @a ncomp
804 /// components of the size of the scalar FiniteElement, reorder its entries
805 /// into native (H1) ordering.
806 /// The function assumes that GetLexicographicOrdering() is not empty.
807 void ReorderLexToNative(int ncomp, Vector &dofs) const;
808};
809
810/** @brief Intermediate class for finite elements whose basis functions return
811 vector values. */
813{
814 // Hide the scalar functions CalcShape and CalcDShape.
815private:
816 /// Overrides the scalar CalcShape function to print an error.
817 void CalcShape(const IntegrationPoint &ip,
818 Vector &shape) const override;
819
820 /// Overrides the scalar CalcDShape function to print an error.
821 void CalcDShape(const IntegrationPoint &ip,
822 DenseMatrix &dshape) const override;
823
824protected:
826#ifndef MFEM_THREAD_SAFE
829#endif
830 void SetDerivMembers();
831
833 DenseMatrix &shape) const;
834
836 DenseMatrix &shape) const;
837
838 /** @brief Project a vector coefficient onto the RT basis functions
839 @param nk Face normal vectors for this element type
840 @param d2n Offset into nk for each degree of freedom
841 @param vc Vector coefficient to be projected
842 @param Trans Transformation from reference to physical coordinates
843 @param dofs Expansion coefficients for the approximation of vc
844 */
845 void Project_RT(const real_t *nk, const Array<int> &d2n,
847 Vector &dofs) const;
848
849 /// Projects the vector of values given at FE nodes to RT space
850 /** Project vector values onto the RT basis functions
851 @param nk Face normal vectors for this element type
852 @param d2n Offset into nk for each degree of freedom
853 @param vc Vector values at each interpolation point
854 @param Trans Transformation from reference to physical coordinates
855 @param dofs Expansion coefficients for the approximation of vc
856 */
857 void Project_RT(const real_t *nk, const Array<int> &d2n,
858 Vector &vc, ElementTransformation &Trans,
859 Vector &dofs) const;
860
861 /// Project the rows of the matrix coefficient in an RT space
863 const real_t *nk, const Array<int> &d2n,
864 MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const;
865
866 /** @brief Project vector-valued basis functions onto the RT basis functions
867 @param nk Face normal vectors for this element type
868 @param d2n Offset into nk for each degree of freedom
869 @param fe Vector-valued finite element basis
870 @param Trans Transformation from reference to physical coordinates
871 @param I Expansion coefficients for the approximation of each basis
872 function
873
874 Note: If the FiniteElement, fe, is scalar-valued the projection will
875 assume that a FiniteElementSpace is being used to define a vector
876 field using the scalar basis functions for each component of the
877 vector field.
878 */
879 void Project_RT(const real_t *nk, const Array<int> &d2n,
880 const FiniteElement &fe, ElementTransformation &Trans,
881 DenseMatrix &I) const;
882
883 // rotated gradient in 2D
884 void ProjectGrad_RT(const real_t *nk, const Array<int> &d2n,
885 const FiniteElement &fe, ElementTransformation &Trans,
886 DenseMatrix &grad) const;
887
888 // Compute the curl as a discrete operator from ND FE (fe) to ND FE (this).
889 // The natural FE for the range is RT, so this is an approximation.
890 void ProjectCurl_ND(const real_t *tk, const Array<int> &d2t,
891 const FiniteElement &fe, ElementTransformation &Trans,
892 DenseMatrix &curl) const;
893
894 void ProjectCurl_RT(const real_t *nk, const Array<int> &d2n,
895 const FiniteElement &fe, ElementTransformation &Trans,
896 DenseMatrix &curl) const;
897
898 /** @brief Project a vector coefficient onto the ND basis functions
899 @param tk Edge tangent vectors for this element type
900 @param d2t Offset into tk for each degree of freedom
901 @param vc Vector coefficient to be projected
902 @param Trans Transformation from reference to physical coordinates
903 @param dofs Expansion coefficients for the approximation of vc
904 */
905 void Project_ND(const real_t *tk, const Array<int> &d2t,
907 Vector &dofs) const;
908
909 /// Projects the vector of values given at FE nodes to ND space
910 /** Project vector values onto the ND basis functions
911 @param tk Edge tangent vectors for this element type
912 @param d2t Offset into tk for each degree of freedom
913 @param vc Vector values at each interpolation point
914 @param Trans Transformation from reference to physical coordinates
915 @param dofs Expansion coefficients for the approximation of vc
916 */
917 void Project_ND(const real_t *tk, const Array<int> &d2t,
918 Vector &vc, ElementTransformation &Trans,
919 Vector &dofs) const;
920
921 /// Project the rows of the matrix coefficient in an ND space
923 const real_t *tk, const Array<int> &d2t,
924 MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const;
925
926 /** @brief Project vector-valued basis functions onto the ND basis functions
927 @param tk Edge tangent vectors for this element type
928 @param d2t Offset into tk for each degree of freedom
929 @param fe Vector-valued finite element basis
930 @param Trans Transformation from reference to physical coordinates
931 @param I Expansion coefficients for the approximation of each basis
932 function
933
934 Note: If the FiniteElement, fe, is scalar-valued the projection will
935 assume that a FiniteElementSpace is being used to define a vector
936 field using the scalar basis functions for each component of the
937 vector field.
938 */
939 void Project_ND(const real_t *tk, const Array<int> &d2t,
940 const FiniteElement &fe, ElementTransformation &Trans,
941 DenseMatrix &I) const;
942
943 void ProjectGrad_ND(const real_t *tk, const Array<int> &d2t,
944 const FiniteElement &fe, ElementTransformation &Trans,
945 DenseMatrix &grad) const;
946
949 DenseMatrix &I) const;
950
952 const real_t *nk, const Array<int> &d2n,
954 DenseMatrix &I) const;
955
958 DenseMatrix &I) const;
959
961 const real_t *tk, const Array<int> &d2t,
963 DenseMatrix &I) const;
964
965 void LocalRestriction_RT(const real_t *nk, const Array<int> &d2n,
967 DenseMatrix &R) const;
968
969 void LocalRestriction_ND(const real_t *tk, const Array<int> &d2t,
971 DenseMatrix &R) const;
972
974 {
975 if (fe.GetRangeType() != VECTOR)
976 { mfem_error("'fe' must be a VectorFiniteElement"); }
977 return static_cast<const VectorFiniteElement &>(fe);
978 }
979
980public:
981 VectorFiniteElement(int D, Geometry::Type G, int Do, int O, int M,
982 int F = FunctionSpace::Pk);
983};
984
985/// @brief Class for computing 1D special polynomials and their associated basis
986/// functions
988{
989public:
990 /// One-dimensional basis evaluation type
992 {
993 ChangeOfBasis = 0, ///< Use change of basis, O(p^2) Evals
994 Barycentric = 1, ///< Use barycentric Lagrangian interpolation, O(p) Evals
995 Positive = 2, ///< Fast evaluation of Bernstein polynomials
996 Integrated = 3, ///< Integrated indicator functions (cf. Gerritsma)
997 NumEvalTypes = 4 ///< Keep count of the number of eval types
998 };
999
1000 /// @brief Class for evaluating 1D nodal, positive (Bernstein), or integrated
1001 /// (Gerritsma) bases.
1002 class Basis
1003 {
1004 private:
1005 EvalType etype; ///< Determines how the basis functions should be evaluated.
1007 mutable Vector x, w;
1008 /// The following data members are used for "integrated basis type", which
1009 /// is defined in terms of nodal basis of one degree higher.
1010 ///@{
1011 mutable Vector u_aux, d_aux, d2_aux;
1012 ///@}
1013 /// @brief An auxiliary nodal basis used to evaluate the integrated basis.
1014 /// This member variable is NULL whenever etype != Integrated.
1015 Basis *auxiliary_basis;
1016 /// Should the integrated basis functions be scaled? See ScaleIntegrated.
1017 bool scale_integrated;
1018
1019 public:
1020 /// Create a nodal or positive (Bernstein) basis of degree @a p
1021 Basis(const int p, const real_t *nodes, EvalType etype = Barycentric);
1022 /// Evaluate the basis functions at point @a x in [0,1]
1023 void Eval(const real_t x, Vector &u) const;
1024 /// @brief Evaluate the basis functions and their derivatives at point @a
1025 /// x in [0,1]
1026 void Eval(const real_t x, Vector &u, Vector &d) const;
1027 /// @brief Evaluate the basis functions and their first two derivatives at
1028 /// point @a x in [0,1]
1029 void Eval(const real_t x, Vector &u, Vector &d, Vector &d2) const;
1030 /// @brief Evaluate the "integrated" basis type using pre-computed closed
1031 /// basis derivatives.
1032 ///
1033 /// This basis is given by the negative partial sum of the corresponding
1034 /// closed basis derivatives. The closed basis derivatives are given by @a
1035 /// d, and the result is stored in @a i.
1036 void EvalIntegrated(const Vector &d, Vector &i) const;
1037 /// @brief Set whether the "integrated" basis should be scaled by the
1038 /// subcell sizes. Has no effect for non-integrated bases.
1039 ///
1040 /// Generally, this should be true for mfem::FiniteElement::MapType VALUE
1041 /// and false for all other map types. If this option is enabled, the
1042 /// basis functions will be scaled by the widths of the subintervals, so
1043 /// that the basis functions represent mean values. Otherwise, the basis
1044 /// functions represent integrated values.
1045 void ScaleIntegrated(bool scale_integrated_);
1046 /// Returns true if the basis is "integrated", false otherwise.
1047 bool IsIntegratedType() const { return etype == Integrated; }
1048 ~Basis();
1049 };
1050
1051private:
1052 /// key: (btype, p), value: underlying storage Array
1053 typedef std::unordered_map<std::pair<int, int>,
1054 std::unique_ptr<Basis>, PairHasher>
1055 BasisMap;
1056 /// key: (btype, p), value: underlying storage Array
1057 typedef std::unordered_map<std::pair<int, int>,
1058 std::unique_ptr<Array<real_t>>, PairHasher>
1059 PointsMap;
1060
1061 MemoryType h_mt;
1062 PointsMap points_container;
1063 BasisMap bases_container;
1064
1065 static Array2D<int> binom;
1066
1067 static void CalcMono(const int p, const real_t x, real_t *u);
1068 static void CalcMono(const int p, const real_t x, real_t *u, real_t *d);
1069
1070 static void CalcChebyshev(const int p, const real_t x, real_t *u);
1071 static void CalcChebyshev(const int p, const real_t x, real_t *u, real_t *d);
1072 static void CalcChebyshev(const int p, const real_t x, real_t *u, real_t *d,
1073 real_t *dd);
1074
1075 QuadratureFunctions1D quad_func;
1076
1077public:
1079
1080 /** @brief Get a pointer to an array containing the binomial coefficients "p
1081 choose k" for k=0,...,p for the given p. */
1082 static const int *Binom(const int p);
1083
1084 /** @brief Get the coordinates of the points of the given BasisType,
1085 @a btype.
1086
1087 @param[in] p The polynomial degree; the number of points is `p+1`.
1088 @param[in] btype The BasisType.
1089
1090 @return A pointer to an array containing the `p+1` coordinates of the
1091 points. Returns NULL if the BasisType has no associated set of
1092 points. */
1093 const Array<real_t>* GetPointsArray(const int p, const int btype);
1094
1095 /** @brief Get the coordinates of the points of the given BasisType,
1096 @a btype.
1097
1098 @param[in] p The polynomial degree; the number of points is `p+1`.
1099 @param[in] btype The BasisType.
1100 @param[in] on_device true if the requested pointer should be accessible
1101 from the device.
1102
1103 @return A pointer to an array containing the `p+1` coordinates of the
1104 points. Returns NULL if the BasisType has no associated set of
1105 points. */
1106 const real_t *GetPoints(const int p, const int btype,
1107 bool on_device = false)
1108 {
1109 return GetPointsArray(p, btype)->Read(on_device);
1110 }
1111
1112 /// Get coordinates of an open (GaussLegendre) set of points if degree @a p
1113 const real_t *OpenPoints(const int p,
1114 const int btype = BasisType::GaussLegendre,
1115 bool on_device = false)
1116 {
1117 return GetPoints(p, btype, on_device);
1118 }
1119
1120 /// Get coordinates of a closed (GaussLegendre) set of points if degree @a p
1121 const real_t *ClosedPoints(const int p,
1122 const int btype = BasisType::GaussLobatto,
1123 bool on_device = false)
1124 {
1125 return GetPoints(p, btype, on_device);
1126 }
1127
1128 /** @brief Get a Poly_1D::Basis object of the given degree and BasisType,
1129 @a btype.
1130
1131 @param[in] p The polynomial degree of the basis.
1132 @param[in] btype The BasisType.
1133
1134 @return A reference to an object of type Poly_1D::Basis that represents
1135 the requested basis type. */
1136 Basis &GetBasis(const int p, const int btype);
1137
1138 /** @brief Evaluate the values of a hierarchical 1D basis at point x
1139 hierarchical = k-th basis function is degree k polynomial */
1140 static void CalcBasis(const int p, const real_t x, real_t *u)
1141 // { CalcMono(p, x, u); }
1142 // Bernstein basis is not hierarchical --> does not work for triangles
1143 // and tetrahedra
1144 // { CalcBernstein(p, x, u); }
1145 // { CalcLegendre(p, x, u); }
1146 { CalcChebyshev(p, x, u); }
1147
1148 /** @brief Evaluate the values of a hierarchical 1D basis at point x
1149 hierarchical = k-th basis function is degree k polynomial */
1150 static void CalcBasis(const int p, const real_t x, Vector &u)
1151 { CalcBasis(p, x, u.GetData()); }
1152
1153 /// Evaluate the values and derivatives of a hierarchical 1D basis at point @a x
1154 static void CalcBasis(const int p, const real_t x, real_t *u, real_t *d)
1155 // { CalcMono(p, x, u, d); }
1156 // { CalcBernstein(p, x, u, d); }
1157 // { CalcLegendre(p, x, u, d); }
1158 { CalcChebyshev(p, x, u, d); }
1159
1160 /** @brief Evaluate the values and derivatives of a hierarchical 1D basis at
1161 point @a x. */
1162 static void CalcBasis(const int p, const real_t x, Vector &u, Vector &d)
1163 { CalcBasis(p, x, u.GetData(), d.GetData()); }
1164
1165 /// Evaluate the values, derivatives and second derivatives of a hierarchical 1D basis at point x
1166 static void CalcBasis(const int p, const real_t x, real_t *u, real_t *d,
1167 real_t *dd)
1168 // { CalcMono(p, x, u, d); }
1169 // { CalcBernstein(p, x, u, d); }
1170 // { CalcLegendre(p, x, u, d); }
1171 { CalcChebyshev(p, x, u, d, dd); }
1172
1173 /** @brief Evaluate the values, derivatives and second derivatives of a
1174 hierarchical 1D basis at point @a x. */
1175 static void CalcBasis(const int p, const real_t x, Vector &u, Vector &d,
1176 Vector &dd)
1177 { CalcBasis(p, x, u.GetData(), d.GetData(), dd.GetData()); }
1178
1179 /// Evaluate a representation of a Delta function at point x
1180 static real_t CalcDelta(const int p, const real_t x)
1181 { return pow(x, (real_t) p); }
1182
1183 /** @brief Compute the points for the Chebyshev polynomials of order @a p
1184 and place them in the already allocated @a x array. */
1185 static void ChebyshevPoints(const int p, real_t *x);
1186
1187 /** @brief Compute the @a p terms in the expansion of the binomial (x + y)^p
1188 and store them in the already allocated @a u array. */
1189 static void CalcBinomTerms(const int p, const real_t x, const real_t y,
1190 real_t *u);
1191 /** @brief Compute the terms in the expansion of the binomial (x + y)^p and
1192 their derivatives with respect to x assuming that dy/dx = -1. Store the
1193 results in the already allocated @a u and @a d arrays.*/
1194 static void CalcBinomTerms(const int p, const real_t x, const real_t y,
1195 real_t *u, real_t *d);
1196 /** @brief Compute the derivatives (w.r.t. x) of the terms in the expansion
1197 of the binomial (x + y)^p assuming that dy/dx = -1. Store the results
1198 in the already allocated @a d array.*/
1199 static void CalcDBinomTerms(const int p, const real_t x, const real_t y,
1200 real_t *d);
1201 /** @brief Compute the derivatives (w.r.t. x) of the terms in the expansion
1202 of the binomial (x + y)^p. Store the results in the already allocated
1203 @a d array.*/
1204 static void CalcDxBinomTerms(const int p, const real_t x, const real_t y,
1205 real_t *d);
1206 /** @brief Compute the derivatives (w.r.t. y) of the terms in the expansion
1207 of the binomial (x + y)^p. Store the results in the already allocated
1208 @a d array.*/
1209 static void CalcDyBinomTerms(const int p, const real_t x, const real_t y,
1210 real_t *d);
1211
1212 /** @brief Compute the values of the Bernstein basis functions of order
1213 @a p at coordinate @a x and store the results in the already allocated
1214 @a u array. */
1215 static void CalcBernstein(const int p, const real_t x, real_t *u)
1216 { CalcBinomTerms(p, x, 1. - x, u); }
1217
1218 /** @brief Compute the values of the Bernstein basis functions of order
1219 @a p at coordinate @a x and store the results in the already allocated
1220 @a u array. */
1221 static void CalcBernstein(const int p, const real_t x, Vector &u)
1222 { CalcBernstein(p, x, u.GetData()); }
1223
1224 /** @brief Compute the values and derivatives of the Bernstein basis functions
1225 of order @a p at coordinate @a x and store the results in the already allocated
1226 @a u and @a d arrays. */
1227 static void CalcBernstein(const int p, const real_t x, real_t *u, real_t *d)
1228 { CalcBinomTerms(p, x, 1. - x, u, d); }
1229
1230 /** @brief Compute the values and derivatives of the Bernstein basis
1231 functions of order @a p at coordinate @a x and store the results in the
1232 already allocated @a u and @a d arrays. */
1233 static void CalcBernstein(const int p, const real_t x, Vector &u, Vector &d)
1234 { CalcBernstein(p, x, u.GetData(), d.GetData()); }
1235
1236 static void CalcLegendre(const int p, const real_t x, real_t *u);
1237 static void CalcLegendre(const int p, const real_t x, real_t *u, real_t *d);
1238
1239 ~Poly_1D() = default;
1240};
1241
1242extern MFEM_EXPORT Poly_1D poly1d;
1243
1244/// An element defined as an ND tensor product of 1D elements on a segment,
1245/// square, or cube
1247{
1248protected:
1253
1254public:
1256 {
1259 Sr_DOF_MAP = 2, // Sr = Serendipity
1260 };
1261
1262 TensorBasisElement(const int dims, const int p, const int btype,
1263 const DofMapType dmtype);
1264
1265 int GetBasisType() const { return b_type; }
1266
1267 const Poly_1D::Basis &GetBasis1D() const { return basis1d; }
1268
1269 /** @brief Get an Array<int> that maps lexicographically ordered indices to
1270 the indices of the respective nodes/dofs/basis functions. If the dofs are
1271 ordered lexicographically, i.e. the mapping is identity, the returned
1272 Array will be empty. */
1273 const Array<int> &GetDofMap() const { return dof_map; }
1274
1276 {
1277 switch (dim)
1278 {
1279 case 1: return Geometry::SEGMENT;
1280 case 2: return Geometry::SQUARE;
1281 case 3: return Geometry::CUBE;
1282 default:
1283 MFEM_ABORT("invalid dimension: " << dim);
1284 return Geometry::INVALID;
1285 }
1286 }
1287
1288 /// Return @a base raised to the power @a dim.
1289 static int Pow(int base, int dim)
1290 {
1291 switch (dim)
1292 {
1293 case 1: return base;
1294 case 2: return base*base;
1295 case 3: return base*base*base;
1296 default: MFEM_ABORT("invalid dimension: " << dim); return -1;
1297 }
1298 }
1299
1300 static const DofToQuad &GetTensorDofToQuad(
1301 const FiniteElement &fe, const IntegrationRule &ir,
1302 DofToQuad::Mode mode, const Poly_1D::Basis &basis, bool closed,
1303 Array<DofToQuad*> &dof2quad_array);
1304};
1305
1307 public TensorBasisElement
1308{
1309public:
1310 NodalTensorFiniteElement(const int dims, const int p, const int btype,
1311 const DofMapType dmtype);
1312
1313 const DofToQuad &GetDofToQuad(const IntegrationRule &ir,
1314 DofToQuad::Mode mode) const override;
1315
1316 void SetMapType(const int map_type_) override;
1317
1319 ElementTransformation &Trans,
1320 DenseMatrix &I) const override
1321 {
1323 {
1324 CheckScalarFE(fe).ScalarLocalInterpolation(Trans, I, *this);
1325 }
1326 else
1327 {
1329 }
1330 }
1331
1332 void GetFaceMap(const int face_id, Array<int> &face_map) const override;
1333};
1334
1336 public TensorBasisElement
1337{
1338private:
1339 mutable Array<DofToQuad*> dof2quad_array_open;
1340
1341protected:
1343
1344public:
1345 VectorTensorFiniteElement(const int dims, const int d, const int p,
1346 const int cbtype, const int obtype,
1347 const int M, const DofMapType dmtype);
1348
1349 // For 1D elements: there is only an "open basis", no "closed basis"
1350 VectorTensorFiniteElement(const int dims, const int d, const int p,
1351 const int obtype, const int M,
1352 const DofMapType dmtype);
1353
1355 DofToQuad::Mode mode) const override
1356 {
1357 return (mode == DofToQuad::TENSOR) ?
1358 GetTensorDofToQuad(*this, ir, mode, basis1d, true, dof2quad_array) :
1360 }
1361
1363 DofToQuad::Mode mode) const
1364 {
1365 MFEM_VERIFY(mode == DofToQuad::TENSOR, "invalid mode requested");
1366 return GetTensorDofToQuad(*this, ir, mode, obasis1d, false,
1367 dof2quad_array_open);
1368 }
1369
1371};
1372
1373void InvertLinearTrans(ElementTransformation &trans,
1374 const IntegrationPoint &pt, Vector &x);
1375
1376} // namespace mfem
1377
1378#endif
Dynamic 2D array using row-major layout.
Definition array.hpp:392
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition array.hpp:337
Possible basis types. Note that not all elements can use all BasisType(s).
Definition fe_base.hpp:30
static int GetQuadrature1D(int b_type)
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
Definition fe_base.hpp:65
static int Check(int b_type)
If the input does not represents a valid BasisType, abort with an error; otherwise return the input.
Definition fe_base.hpp:49
static char GetChar(int b_type)
Check and convert a BasisType constant to a char basis identifier.
Definition fe_base.hpp:108
@ ClosedGL
Closed GaussLegendre.
Definition fe_base.hpp:42
@ OpenHalfUniform
Nodes: x_i = (i+1/2)/n, i=0,...,n-1.
Definition fe_base.hpp:40
@ Serendipity
Serendipity basis (squares / cubes)
Definition fe_base.hpp:41
@ GaussLobatto
Closed type.
Definition fe_base.hpp:36
@ GaussLegendre
Open type.
Definition fe_base.hpp:35
@ Positive
Bernstein polynomials.
Definition fe_base.hpp:37
@ OpenUniform
Nodes: x_i = (i+1)/(n+1), i=0,...,n-1.
Definition fe_base.hpp:38
@ IntegratedGLL
Integrated GLL indicator functions.
Definition fe_base.hpp:43
@ ClosedUniform
Nodes: x_i = i/(n-1), i=0,...,n-1.
Definition fe_base.hpp:39
static int CheckNodal(int b_type)
If the input does not represents a valid nodal BasisType, abort with an error; otherwise return the i...
Definition fe_base.hpp:57
static const char * Name(int b_type)
Check and convert a BasisType constant to a string identifier.
Definition fe_base.hpp:96
static int GetType(char b_ident)
Convert char basis identifier to a BasisType constant.
Definition fe_base.hpp:115
static int GetNodalBasis(int qpt_type)
Return the nodal BasisType corresponding to the Quadrature1D type.
Definition fe_base.hpp:82
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
Structure representing the matrices/tensors needed to evaluate (in reference space) the values,...
Definition fe_base.hpp:141
Mode mode
Describes the contents of the B, Bt, G, and Gt arrays, see Mode.
Definition fe_base.hpp:174
const IntegrationRule * IntRule
IntegrationRule that defines the quadrature points at which the basis functions of the FE are evaluat...
Definition fe_base.hpp:150
Array< real_t > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition fe_base.hpp:214
Mode
Type of data stored in the arrays B, Bt, G, and Gt.
Definition fe_base.hpp:154
@ FULL
Full multidimensional representation which does not use tensor product structure. The ordering of the...
Definition fe_base.hpp:158
@ LEXICOGRAPHIC_FULL
Full multidimensional representation which does not use tensor product structure. The ordering of the...
Definition fe_base.hpp:170
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:165
Array< real_t > B
Basis functions evaluated at quadrature points.
Definition fe_base.hpp:193
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:178
Array< real_t > Gt
Transpose of G.
Definition fe_base.hpp:221
const class FiniteElement * FE
The FiniteElement that created and owns this object.
Definition fe_base.hpp:145
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:182
Array< real_t > Bt
Transpose of B.
Definition fe_base.hpp:199
Abstract class for all finite elements.
Definition fe_base.hpp:244
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_base.cpp:40
static int VerifyNodal(int b_type)
Ensure that the BasisType of b_type nodal (satisfies the interpolation property).
Definition fe_base.hpp:653
int dof
Number of degrees of freedom.
Definition fe_base.hpp:253
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_base.cpp:101
virtual void ProjectMatrixCoefficient(MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Given a matrix coefficient and a transformation, compute an approximation ("projection") in the local...
Definition fe_base.cpp:144
virtual ~FiniteElement()
Deconstruct the FiniteElement.
Definition fe_base.cpp:513
int GetDerivMapType() const
Returns the FiniteElement::DerivType of the element describing how reference function derivatives are...
Definition fe_base.hpp:370
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
Compute the discrete divergence matrix from the given FiniteElement onto 'this' FiniteElement....
Definition fe_base.cpp:175
virtual void ProjectFromNodes(Vector &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector of values at the finite element nodes and a transformation, compute its projection (ap...
Definition fe_base.cpp:138
virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const
Get the dofs associated with the given face. *dofs is set to an internal array of the local dofc on t...
Definition fe_base.cpp:96
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition fe_base.cpp:365
RangeType
Enumeration for range_type and deriv_range_type.
Definition fe_base.hpp:267
static bool IsOpenType(int b_type)
Return true if the BasisType of b_type is open (doesn't have Quadrature1D points on the boundary).
Definition fe_base.hpp:627
int GetRangeDim() const
Returns the vector dimension for vector-valued finite elements, which is also the dimension of the in...
Definition fe_base.hpp:325
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition fe_base.hpp:338
void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in physical space at the poi...
Definition fe_base.cpp:192
int GetDerivType() const
Returns the FiniteElement::DerivType of the element describing the spatial derivative method implemen...
Definition fe_base.hpp:365
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_base.cpp:65
virtual void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Compute the discrete gradient matrix from the given FiniteElement onto 'this' FiniteElement....
Definition fe_base.cpp:161
IntegrationRule Nodes
Definition fe_base.hpp:256
virtual void GetFaceMap(const int face_id, Array< int > &face_map) const
Return the mapping from lexicographic face DOFs to lexicographic element DOFs for the given local fac...
Definition fe_base.cpp:507
virtual const StatelessDofTransformation * GetDofTransformation() const
Return a DoF transformation object for this particular type of basis.
Definition fe_base.hpp:610
int GetDim() const
Returns the reference space dimension for the finite element.
Definition fe_base.hpp:321
static int VerifyClosed(int b_type)
Ensure that the BasisType of b_type is closed (has Quadrature1D points on the boundary).
Definition fe_base.hpp:636
MapType
Enumeration for MapType: defines how reference functions are mapped to physical space.
Definition fe_base.hpp:280
int vdim
Vector dimension of vector-valued basis functions.
Definition fe_base.hpp:247
void CalcPhysHessian(ElementTransformation &Trans, DenseMatrix &Hessian) const
Evaluate the Hessian of all shape functions of a scalar finite element in physical space at the given...
Definition fe_base.cpp:288
FiniteElement(int D, Geometry::Type G, int Do, int O, int F=FunctionSpace::Pk)
Construct FiniteElement with given.
Definition fe_base.cpp:23
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_base.cpp:150
int orders[Geometry::MaxDim]
Anisotropic orders.
Definition fe_base.hpp:255
virtual void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
Definition fe_base.cpp:119
static bool IsClosedType(int b_type)
Return true if the BasisType of b_type is closed (has Quadrature1D points on the boundary).
Definition fe_base.hpp:618
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition fe_base.hpp:360
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition fe_base.hpp:351
virtual void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs.
Definition fe_base.cpp:113
const int * GetAnisotropicOrders() const
Returns an array containing the anisotropic orders/degrees.
Definition fe_base.hpp:345
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:400
static int VerifyOpen(int b_type)
Ensure that the BasisType of b_type is open (doesn't have Quadrature1D points on the boundary).
Definition fe_base.hpp:645
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
bool HasAnisotropicOrders() const
Returns true if the FiniteElement basis may be using different orders/degrees in different spatial di...
Definition fe_base.hpp:342
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition fe_base.hpp:331
int cdim
Dimension of curl for vector-valued basis functions.
Definition fe_base.hpp:248
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
Definition fe_base.cpp:52
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
Definition fe_base.cpp:107
int Space() const
Returns the type of FunctionSpace on the element.
Definition fe_base.hpp:348
DerivType
Enumeration for DerivType: defines which derivative method is implemented.
Definition fe_base.hpp:303
@ DIV
Implements CalcDivShape methods.
Definition fe_base.hpp:306
@ NONE
No derivatives implemented.
Definition fe_base.hpp:304
@ CURL
Implements CalcCurlShape methods.
Definition fe_base.hpp:307
@ GRAD
Implements CalcDShape methods.
Definition fe_base.hpp:305
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition fe_base.cpp:126
Geometry::Type geom_type
Geometry::Type of the reference element.
Definition fe_base.hpp:249
void CalcPhysDivShape(ElementTransformation &Trans, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in physical space at the po...
Definition fe_base.cpp:58
Array< DofToQuad * > dof2quad_array
Container for all DofToQuad objects created by the FiniteElement.
Definition fe_base.hpp:263
void CalcPhysLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Evaluate the Laplacian of all shape functions of a scalar finite element in physical space at the giv...
Definition fe_base.cpp:203
void CalcPhysVShape(ElementTransformation &Trans, DenseMatrix &shape) const
Equivalent to the CalcVShape() method with the same arguments.
Definition fe_base.hpp:447
DenseMatrix vshape
Definition fe_base.hpp:258
int GetCurlDim() const
Returns the dimension of the curl for vector-valued finite elements.
Definition fe_base.hpp:328
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void ProjectCurl(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Compute the discrete curl matrix from the given FiniteElement onto 'this' FiniteElement....
Definition fe_base.cpp:168
virtual void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
Definition fe_base.cpp:71
int GetDerivRangeType() const
Returns the FiniteElement::RangeType of the element derivative, either SCALAR or VECTOR.
Definition fe_base.hpp:355
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:334
int order
Order/degree of the shape functions.
Definition fe_base.hpp:254
void CalcPhysLinLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Evaluate the Laplacian of all shape functions of a scalar finite element in physical space at the giv...
Definition fe_base.cpp:244
int dim
Dimension of reference space.
Definition fe_base.hpp:246
void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in physical space at the point ...
Definition fe_base.cpp:182
Describes the function space on each element.
Definition fe_base.hpp:226
@ Pk
Polynomials of order k.
Definition fe_base.hpp:230
@ Qk
Tensor products of polynomials of order k.
Definition fe_base.hpp:231
@ Uk
Rational polynomials of order k.
Definition fe_base.hpp:233
@ rQk
Refined tensor products of polynomials of order k.
Definition fe_base.hpp:232
static const int MaxDim
Definition geom.hpp:47
Class for integration point with weight.
Definition intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
Base class for Matrix Coefficients that optionally depend on time and space.
Class for standard nodal finite elements.
Definition fe_base.hpp:721
NodalFiniteElement(int D, Geometry::Type G, int Do, int O, int F=FunctionSpace::Pk)
Construct NodalFiniteElement with given.
Definition fe_base.hpp:739
void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const override
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition fe_base.cpp:796
const Array< int > & GetLexicographicOrdering() const
Get an Array<int> that maps lexicographically ordered indices to the indices of the respective nodes/...
Definition fe_base.hpp:801
void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const override
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
Definition fe_base.hpp:753
void ProjectCurl_2D(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition fe_base.cpp:724
void ReorderLexToNative(int ncomp, Vector &dofs) const
Definition fe_base.cpp:976
void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const override
Compute the discrete divergence matrix from the given FiniteElement onto 'this' FiniteElement....
Definition fe_base.cpp:945
void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const override
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
Definition fe_base.hpp:746
void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const override
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs.
Definition fe_base.cpp:764
void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const override
Compute the discrete gradient matrix from the given FiniteElement onto 'this' FiniteElement....
Definition fe_base.cpp:916
void ProjectMatrixCoefficient(MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const override
Given a matrix coefficient and a transformation, compute an approximation ("projection") in the local...
Definition fe_base.cpp:835
Array< int > lex_ordering
Definition fe_base.hpp:726
const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const override
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition fe_base.cpp:703
void GetFaceMap(const int face_id, Array< int > &face_map) const override
Return the mapping from lexicographic face DOFs to lexicographic element DOFs for the given local fac...
Definition fe_base.cpp:2686
NodalTensorFiniteElement(const int dims, const int p, const int btype, const DofMapType dmtype)
Definition fe_base.cpp:2648
const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const override
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition fe_base.cpp:2672
void SetMapType(const int map_type_) override
Set the FiniteElement::MapType of the element to either VALUE or INTEGRAL. Also sets the FiniteElemen...
Definition fe_base.cpp:2659
void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const override
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
Definition fe_base.hpp:1318
Class for evaluating 1D nodal, positive (Bernstein), or integrated (Gerritsma) bases.
Definition fe_base.hpp:1003
void Eval(const real_t x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition fe_base.cpp:1777
bool IsIntegratedType() const
Returns true if the basis is "integrated", false otherwise.
Definition fe_base.hpp:1047
void ScaleIntegrated(bool scale_integrated_)
Set whether the "integrated" basis should be scaled by the subcell sizes. Has no effect for non-integ...
Definition fe_base.cpp:2034
void EvalIntegrated(const Vector &d, Vector &i) const
Evaluate the "integrated" basis type using pre-computed closed basis derivatives.
Definition fe_base.cpp:2009
Basis(const int p, const real_t *nodes, EvalType etype=Barycentric)
Create a nodal or positive (Bernstein) basis of degree p.
Definition fe_base.cpp:1713
Class for computing 1D special polynomials and their associated basis functions.
Definition fe_base.hpp:988
static void CalcBasis(const int p, const real_t x, Vector &u, Vector &d, Vector &dd)
Evaluate the values, derivatives and second derivatives of a hierarchical 1D basis at point x.
Definition fe_base.hpp:1175
static void CalcDBinomTerms(const int p, const real_t x, const real_t y, real_t *d)
Compute the derivatives (w.r.t. x) of the terms in the expansion of the binomial (x + y)^p assuming t...
Definition fe_base.cpp:2157
static void CalcBernstein(const int p, const real_t x, real_t *u, real_t *d)
Compute the values and derivatives of the Bernstein basis functions of order p at coordinate x and st...
Definition fe_base.hpp:1227
const real_t * GetPoints(const int p, const int btype, bool on_device=false)
Get the coordinates of the points of the given BasisType, btype.
Definition fe_base.hpp:1106
static void CalcBasis(const int p, const real_t x, real_t *u, real_t *d, real_t *dd)
Evaluate the values, derivatives and second derivatives of a hierarchical 1D basis at point x.
Definition fe_base.hpp:1166
static void CalcBernstein(const int p, const real_t x, Vector &u)
Compute the values of the Bernstein basis functions of order p at coordinate x and store the results ...
Definition fe_base.hpp:1221
static const int * Binom(const int p)
Get a pointer to an array containing the binomial coefficients "pchoose k" for k=0,...
Definition fe_base.cpp:2044
Basis & GetBasis(const int p, const int btype)
Get a Poly_1D::Basis object of the given degree and BasisType, btype.
Definition fe_base.cpp:2372
const Array< real_t > * GetPointsArray(const int p, const int btype)
Get the coordinates of the points of the given BasisType, btype.
Definition fe_base.cpp:2344
static void CalcDyBinomTerms(const int p, const real_t x, const real_t y, real_t *d)
Compute the derivatives (w.r.t. y) of the terms in the expansion of the binomial (x + y)^p....
Definition fe_base.cpp:2216
static real_t CalcDelta(const int p, const real_t x)
Evaluate a representation of a Delta function at point x.
Definition fe_base.hpp:1180
static void CalcLegendre(const int p, const real_t x, real_t *u)
Definition fe_base.cpp:2245
static void CalcBernstein(const int p, const real_t x, real_t *u)
Compute the values of the Bernstein basis functions of order p at coordinate x and store the results ...
Definition fe_base.hpp:1215
EvalType
One-dimensional basis evaluation type.
Definition fe_base.hpp:992
@ ChangeOfBasis
Use change of basis, O(p^2) Evals.
Definition fe_base.hpp:993
@ Integrated
Integrated indicator functions (cf. Gerritsma)
Definition fe_base.hpp:996
@ Positive
Fast evaluation of Bernstein polynomials.
Definition fe_base.hpp:995
@ NumEvalTypes
Keep count of the number of eval types.
Definition fe_base.hpp:997
@ Barycentric
Use barycentric Lagrangian interpolation, O(p) Evals.
Definition fe_base.hpp:994
static void CalcBasis(const int p, const real_t x, real_t *u, real_t *d)
Evaluate the values and derivatives of a hierarchical 1D basis at point x.
Definition fe_base.hpp:1154
static void ChebyshevPoints(const int p, real_t *x)
Compute the points for the Chebyshev polynomials of order p and place them in the already allocated x...
Definition fe_base.cpp:2061
const real_t * ClosedPoints(const int p, const int btype=BasisType::GaussLobatto, bool on_device=false)
Get coordinates of a closed (GaussLegendre) set of points if degree p.
Definition fe_base.hpp:1121
const real_t * OpenPoints(const int p, const int btype=BasisType::GaussLegendre, bool on_device=false)
Get coordinates of an open (GaussLegendre) set of points if degree p.
Definition fe_base.hpp:1113
~Poly_1D()=default
static void CalcBasis(const int p, const real_t x, Vector &u, Vector &d)
Evaluate the values and derivatives of a hierarchical 1D basis at point x.
Definition fe_base.hpp:1162
static void CalcBernstein(const int p, const real_t x, Vector &u, Vector &d)
Compute the values and derivatives of the Bernstein basis functions of order p at coordinate x and st...
Definition fe_base.hpp:1233
static void CalcBasis(const int p, const real_t x, real_t *u)
Evaluate the values of a hierarchical 1D basis at point x hierarchical = k-th basis function is degre...
Definition fe_base.hpp:1140
static void CalcBinomTerms(const int p, const real_t x, const real_t y, real_t *u)
Compute the p terms in the expansion of the binomial (x + y)^p and store them in the already allocate...
Definition fe_base.cpp:2093
static void CalcBasis(const int p, const real_t x, Vector &u)
Evaluate the values of a hierarchical 1D basis at point x hierarchical = k-th basis function is degre...
Definition fe_base.hpp:1150
static void CalcDxBinomTerms(const int p, const real_t x, const real_t y, real_t *d)
Compute the derivatives (w.r.t. x) of the terms in the expansion of the binomial (x + y)^p....
Definition fe_base.cpp:2187
static int CheckOpen(int type)
If the Quadrature1D type is not open return Invalid; otherwise return type.
Definition intrules.cpp:941
@ ClosedUniform
aka closed Newton-Cotes
Definition intrules.hpp:408
@ ClosedGL
aka closed Gauss Legendre
Definition intrules.hpp:410
@ OpenHalfUniform
aka "open half" Newton-Cotes
Definition intrules.hpp:409
@ OpenUniform
aka open Newton-Cotes
Definition intrules.hpp:407
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
Definition intrules.cpp:928
A Class that defines 1-D numerical quadrature rules on [0,1].
Definition intrules.hpp:376
Class for finite elements with basis functions that return scalar values.
Definition fe_base.hpp:662
ScalarFiniteElement(int D, Geometry::Type G, int Do, int O, int F=FunctionSpace::Pk)
Construct ScalarFiniteElement with given.
Definition fe_base.hpp:679
void NodalLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
Get the matrix I that defines nodal interpolation between this element and the refined element fine_f...
Definition fe_base.cpp:522
virtual void SetMapType(int M)
Set the FiniteElement::MapType of the element to either VALUE or INTEGRAL. Also sets the FiniteElemen...
Definition fe_base.hpp:687
void ScalarLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
Get matrix I "Interpolation" defined through local L2-projection in the space defined by the fine_fe.
Definition fe_base.cpp:561
void ScalarLocalL2Restriction(ElementTransformation &Trans, DenseMatrix &R, const ScalarFiniteElement &coarse_fe) const
Get restriction matrix R defined through local L2-projection in the space defined by the coarse_fe.
Definition fe_base.cpp:602
static const ScalarFiniteElement & CheckScalarFE(const FiniteElement &fe)
Definition fe_base.hpp:664
const Array< int > & GetDofMap() const
Get an Array<int> that maps lexicographically ordered indices to the indices of the respective nodes/...
Definition fe_base.hpp:1273
const Poly_1D::Basis & GetBasis1D() const
Definition fe_base.hpp:1267
static Geometry::Type GetTensorProductGeometry(int dim)
Definition fe_base.hpp:1275
Poly_1D::Basis & basis1d
Definition fe_base.hpp:1251
static const DofToQuad & GetTensorDofToQuad(const FiniteElement &fe, const IntegrationRule &ir, DofToQuad::Mode mode, const Poly_1D::Basis &basis, bool closed, Array< DofToQuad * > &dof2quad_array)
Definition fe_base.cpp:2599
TensorBasisElement(const int dims, const int p, const int btype, const DofMapType dmtype)
Definition fe_base.cpp:2399
static int Pow(int base, int dim)
Return base raised to the power dim.
Definition fe_base.hpp:1289
Base class for vector Coefficients that optionally depend on time and space.
Intermediate class for finite elements whose basis functions return vector values.
Definition fe_base.hpp:813
void ProjectGrad_RT(const real_t *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Definition fe_base.cpp:1222
void LocalRestriction_RT(const real_t *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &R) const
Definition fe_base.cpp:1628
void Project_ND(const real_t *tk, const Array< int > &d2t, VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Project a vector coefficient onto the ND basis functions.
Definition fe_base.cpp:1306
void Project_RT(const real_t *nk, const Array< int > &d2n, VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Project a vector coefficient onto the RT basis functions.
Definition fe_base.cpp:1082
void LocalInterpolation_RT(const VectorFiniteElement &cfe, const real_t *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &I) const
Definition fe_base.cpp:1502
void ProjectMatrixCoefficient_ND(const real_t *tk, const Array< int > &d2t, MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Project the rows of the matrix coefficient in an ND space.
Definition fe_base.cpp:1335
void LocalRestriction_ND(const real_t *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &R) const
Definition fe_base.cpp:1671
void LocalL2Projection_RT(const VectorFiniteElement &cfe, ElementTransformation &Trans, DenseMatrix &I) const
Definition fe_base.cpp:1455
void ProjectGrad_ND(const real_t *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Definition fe_base.cpp:1434
void CalcVShape_ND(ElementTransformation &Trans, DenseMatrix &shape) const
Definition fe_base.cpp:1071
void CalcVShape_RT(ElementTransformation &Trans, DenseMatrix &shape) const
Definition fe_base.cpp:1059
VectorFiniteElement(int D, Geometry::Type G, int Do, int O, int M, int F=FunctionSpace::Pk)
Definition fe_base.cpp:993
void LocalInterpolation_ND(const VectorFiniteElement &cfe, const real_t *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &I) const
Definition fe_base.cpp:1589
void ProjectCurl_RT(const real_t *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition fe_base.cpp:1287
static const VectorFiniteElement & CheckVectorFE(const FiniteElement &fe)
Definition fe_base.hpp:973
void LocalL2Projection_ND(const VectorFiniteElement &cfe, ElementTransformation &Trans, DenseMatrix &I) const
Definition fe_base.cpp:1543
void ProjectCurl_ND(const real_t *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition fe_base.cpp:1249
void ProjectMatrixCoefficient_RT(const real_t *nk, const Array< int > &d2n, MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Project the rows of the matrix coefficient in an RT space.
Definition fe_base.cpp:1119
const DofToQuad & GetDofToQuadOpen(const IntegrationRule &ir, DofToQuad::Mode mode) const
Definition fe_base.hpp:1362
VectorTensorFiniteElement(const int dims, const int d, const int p, const int cbtype, const int obtype, const int M, const DofMapType dmtype)
Definition fe_base.cpp:2692
const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const override
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition fe_base.hpp:1354
Vector data type.
Definition vector.hpp:82
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:235
int dim
Definition ex24.cpp:53
void trans(const Vector &u, Vector &x)
Definition ex27.cpp:412
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void mfem_error(const char *msg)
Definition error.cpp:154
float real_t
Definition config.hpp:43
MemoryType
Memory types supported by MFEM.
@ HOST
Host memory; using new[] and delete[].
Poly_1D poly1d
Definition fe.cpp:28
void InvertLinearTrans(ElementTransformation &trans, const IntegrationPoint &pt, Vector &x)
Definition fe_base.cpp:748
real_t p(const Vector &x, real_t t)
Helper class for hashing std::pair. Usable in place of std::hash<std::pair<T,U>>
Definition hash.hpp:495
std::array< int, NCMesh::MaxFaceNodes > nodes