A general approach to transforming finite elements

TL;DR

This paper introduces a general method for transforming finite element bases, enabling the use of reference elements for a broader class of finite elements that previously could not be handled in this way.

Contribution

It develops a unified approach to map bases for complex finite elements like Hermite and Argyris, expanding the applicability of reference element techniques.

Findings

Provides a new framework for finite element basis transformation

Enables efficient implementation of complex finite elements

Broadens the class of finite elements compatible with reference element methods

Abstract

The use of a reference element on which a finite element basis is constructed once and mapped to each cell in a mesh greatly expedites the structure and efficiency of finite element codes. However, many famous finite elements such as Hermite, Morley, Argyris, and Bell, do not possess the kind of equivalence needed to work with a reference element in the standard way. This paper gives a generalizated approach to mapping bases for such finite elements by means of studying relationships between the finite element nodes under push-forward.

Equations

$$\begin{array}{rl}X& \equiv {\left(P\right)}^{\nu },\\ X^{†}& \equiv {\left(C_b^k(K{)}^{\prime }\right)}^{\nu }.\end{array}$$

$$\begin{array}{rl}\hat{X}& \equiv {\left(\hat{P}\right)}^{\nu },\\ {\hat{X}}^{†}& \equiv {\left(C_b^k(\hat{K}{)}^{\prime }\right)}^{\nu }.\end{array}$$

$$F^{\ast }\left(\hat{f}\right)=\hat{f}\circ F$$

$$F_{\ast }(n)=n\circ F^{\ast }$$

$$\begin{array}{rl}{F}_{\ast }(N)\in {\hat{X}}^{†},{\left(F_{\ast }(N)\right)}_i& =F_{\ast }(n_i),\\ F^{\ast }(\hat{\Phi })\in X,{\left(F^{\ast }(\hat{\Phi })\right)}_i& =F^{\ast }(\widehat{{\varphi }_i}).\end{array}$$

$${\left(N(\Phi )\right)}_{ij}=n_i({\varphi }_j).$$

$$N(M\Phi )=N(\Phi )M^T.$$

$${\left(N(M\Phi )\right)}_{ij}=n_i\left({\left(M\Phi \right)}_j\right)=n_i\left(\sum _{k=1}^{\nu }{M}_{jk}{\varphi }_k\right)=\sum _{k=1}^{\nu }{n}_i\left({\varphi }_k\right)=\sum _{k=1}^{\nu }{\left(N(\Phi )\right)}_{ik}{M}_{jk}.$$

$$N(F^{\ast }(\hat{\Phi }))=F_{\ast }(N)(\hat{\Phi })$$

$$\mathcal{I}(f)=\sum _{i=1}^{\nu }{n}_i(f){\psi }_i.$$

$${\delta }_{\mathbf{x}}(p)=p(\mathbf{x}).$$

$${\delta }_{\mathbf{x}}^{\mathbf{s}}(p)={\mathbf{s}}^T\nabla p(\mathbf{x}).$$

$${\nabla }_{\mathbf{x}}={\left[\begin{array}{cc}{\delta }_{\mathbf{x}}^{\mathbf{x}}& {\delta }_{\mathbf{x}}^{\mathbf{y}}\end{array}\right]}^T$$

$${\nabla }_{\mathbf{x}}^{\mathbf{nt}}={\left[\begin{array}{cc}{\delta }_{\mathbf{x}}^{\mathbf{n}}& {\delta }_{\mathbf{x}}^{\mathbf{t}}\end{array}\right]}^T$$

$${△}_{\mathbf{v}}={\left[\begin{array}{ccc}{\delta }_{\mathbf{x}}^{\mathbf{xx}}& {\delta }_{\mathbf{x}}^{\mathbf{xy}}& {\delta }_{\mathbf{x}}^{\mathbf{yy}}\end{array}\right]}^T$$

$$\Psi =M{F}^{\ast }(\hat{\Psi }),$$

$$B_{ij}\equiv F_{\ast }(n_i)({\hat{\psi }}_j)=n_i(F^{\ast }({\hat{\psi }}_j))$$

$$\hat{\pi }\hat{N}=V{F}_{\ast }(\pi N)$$

$$M=V^T.$$

$$I=N(\Psi )=N(M{F}^{\ast }(\hat{\Psi }))=N(F^{\ast }(\hat{\Psi }))M^T=B{M}^T.$$

$$M=B^{-T}.$$

$$I=\left(V{F}_{\ast }(N)\right)(\hat{\Psi })=V{F}_{\ast }(N)(\hat{\Psi })=VB,$$

$$N(F^{\ast }(\hat{\Psi }))=F_{\ast }(N)(\hat{\Psi })=\hat{N}(\hat{\Psi })=I$$

$${\hat{\Psi }}_{iq}={\hat{\psi }}_i({\hat{\xi }}_q)$$

$${\psi }_i({\xi }_q)=\sum _{k=1}^{\nu }{M}_{i,k}{\hat{\Psi }}_{k,q}$$

$$D{\hat{\Psi }}_{i,q,:}=\hat{\nabla }{\hat{\psi }}_i(\hat{\xi }{)}_q.$$

$$D{\Psi }_{i,q,:}^{\prime }:=\sum _{k=1}^{\nu }{M}_{i,k}D{\hat{\Psi }}_{k,q,:},$$

$$D{\Psi }_{i,q,:}:=J^{T}D{\Psi }_{i,q,:}^{\prime }.$$

$$N={\left[\begin{array}{ccccccc}{\delta }_{{\mathbf{v}}_1}& {\nabla }_{{\mathbf{v}}_1}^T& {\delta }_{{\mathbf{v}}_2}& {\nabla }_{{\mathbf{v}}_2}^T& {\delta }_{{\mathbf{v}}_3}& {\nabla }_{{\mathbf{v}}_3}^T& {\delta }_{{\mathbf{v}}_4}\end{array}\right]}^T,$$

$$\nabla (\hat{\psi }\circ F)=J^T\hat{\nabla }\hat{\psi }\circ F.$$

Full text

1 Introduction

At the heart of any finite element implementation lies the evaluation of basis functions and their derivatives on each cell in a mesh. These values are used to compute local integral contributions to stiffness matrices and load vectors, which are assembled into a sparse matrix and then passed on to an algebraic solver. While it is fairly easy to parametrize local integration routines over basis functions, one must also provide an implementation of those basis functions. Frequently, finite element codes use a reference element, on which a set of basis functions is constructed once and mapped via coordinate change to each cell in a mesh. Alternately, many finite element bases can be expressed in terms of barycentric coordinates, in which case one must simply convert between the physical and barycentric coordinates on each cell in order evaluate basis functions. Although we refer the reader to recent results on Bernstein polynomials [1, 22] for interesting algorithms in the latter case, the prevelance of the reference element paradigm in modern high-level finite element software [4, 6, 24, 25, 30, 31] we shall restrict ourselves to the former.

The development of FIAT [21] has had a significant impact on finite element software, especially through its adoption in high-level software projects such as FEniCS [24] and Firedrake [31]. FIAT provides tools to describe and construct reference bases for arbitrary-order instances of many common and unusual finite elements. Composed with a domain-specific language for variational problems like UFL [2] and a form compiler mapping UFL into efficient code for element integrals [18, 23, 26] gives a powerful, user-friendly tool chain.

However, any code based on the reference element paradigm operates under the assumption that finite elements satisfy a certain kind of equivalence. Essentially, one must have a pull-back operation that puts basis functions on each cell into one-to-one correspondence with the reference basis functions. Hence, the original form of ffc [23] used only (arbitrary order) Lagrange finite elements, although this was generalized to $H(\mathrm{div})$ and $H(\mathrm{curl})$ elements using Piola transforms in [32]. Current technology captures the full simplicial discrete de Rham complex and certain other elements, but many famous elements are not included. Although it is possible to construct reference elements in FIAT or some other way, current form compilers or other high-level libraries do not provide correct code for mapping them.

Elements such as Hermite [11], Argyris [3], Morley [28], and Bell [5], shown alongside the Lagrange element in Figure 1, do not satisfy the proper equivalence properties to give a simple relationship between the reference basis and nodal basis on a general cell. Typically, implementations of such elements require special-purpose code for constructing the basis functions separately on each element, which can cost nearly as much in terms of work and storage as building the element stiffness matrix itself. It also requires a different internal workflow in the code. Although Domínguez and Sayas [29] give a technique for mapping bases for the Argyris element and a separate computer implementation is available (https://github.com/VT-ICAM/ArgyrisPack) and Jardin [19] gives a per-element construction technique for the Bell element, these represents the exception rather than the rule. The literature contains no general approach for constructing and mapping finite element bases in the absence of affine equivalence or a suitable generalization thereof.

In this paper we provide such a general theory for transforming finite elements that supplements the theory on which FIAT is based for constructing those elements. Our focus is on the case of scalar-valued elements in affine spaces, although we indicate how the techniques generalize on both counts. We begin the rest of the paper by recalling definitions in § 2. The bulk of the paper occurs in § 3, where we show how to map finite element bases under affine equivalence, affine-interpolation equivalence, and when neither holds. We also sketch briefly how the theory is adapted to the case of more general pullbacks such as non-affine coordinate mappings or Piola transforms. All the theory in § 3 assumes that the natural pull-back operation (i.e. composition with coordinate change) exactly preserves the function spaces between reference and physical space. However, in certain notable cases such as the Bell element, this condition fails to hold. In § 4, we give a more general theory with application to the Bell element. Finally, in § 5, we present some numerical results using these elements.

2 Definitions and preliminaries

Througout, we let $C_{b}^{k}(\Omega)$ denote the space of functions with continuous and bounded derivatives up to and including order $k$ over $\Omega$, and $C_{b}^{k}(\Omega)^{\prime}$ its topological dual.

Definition 2.1**.**

A finite element is a triple $(K,P,N)$ such that

• •

$K\subset\mathbb{R}^{d}$* is a bounded domain.*

• •

$P\subset C^{k}_{b}(K)$* for some integer $k\geq 0$ is a finite-dimensional function space.*

• •

$N=\{n_{i}\}_{i=1}^{\nu}\subset C^{k}_{b}(K)^{\prime}$* is a collection of linearly independent functionals whose actions restricted to $P$ form a basis for $P^{\prime}$.*

The nodes in $N$ are taken as objects in the full infinite-dimensional dual, although sometimes we will only require their restrictions to members of $P$. For any $n\in C^{k}_{b}(K)^{\prime}$, define $\pi n\in P^{\prime}$ by restriction. That is, define $\pi n(p)=n(p)$ for any $p\in P$.

Further, with a slight abuse in notation, we will let $N=\begin{bmatrix}n_{1}&n_{2}&\dots&n_{\nu}\end{bmatrix}^{T}$ denote a functional on $P^{\nu}$, or equivalently, a vector of $\nu$ members of the dual space.

As shorthand, we define these spaces consisting of vectors of functions or functionals by

[TABLE]

We can “vectorize” the restriction operator $\pi$, so that for any $N\in X^{\dagger}$, $\pi N\in(P^{\nu})^{\prime}$ has $(\pi N)_{i}=\pi(n_{i})$.

Galerkin methods work in terms of a basis for the approximating space, and these are typically built out of local bases for each element:

Definition 2.2**.**

Let $(K,P,N)$ be a finite element with $\dim P=\nu$. The nodal basis for $P$ is the set $\{\psi_{i}\}_{i=1}^{\nu}$ such that $n_{i}(\psi_{j})=\delta_{i,j}$ for each $1\leq i,j\leq\nu$.

The nodal basis also can be written as $X\ni\Psi=\begin{bmatrix}\psi_{1}&\psi_{2}&\dots&\psi_{\nu}\end{bmatrix}$.

Traditionally, finite element codes construct the nodal basis for a reference finite element $\left(\hat{K},\hat{P},\hat{N}\right)$ and then map it into the basis for $\left(K,P,N\right)$ for each $K$ in the mesh. Let $F:K\rightarrow\hat{K}$ be the geometric mapping, as in Figure 2. We let $J$ denote the Jacobian matrix of this transformation.

Similarly to (1), we define the vector spaces relative to the reference cell:

[TABLE]

As with $\pi$, we define $\hat{\pi}\hat{n}$ as the restriction of $\hat{n}$ to $\hat{P}$, and can vectorize it over $\hat{X}^{\dagger}$ accordingly.

This geometric mapping induces a mapping between spaces of functions over $K$ and $\hat{K}$ as well as between the dual spaces. These are called the pull-back, and push-forward operations, respectively:

Definition 2.3**.**

The pull-back operation mapping $C^{k}_{b}(\hat{K})\rightarrow C^{k}_{b}(K)$ is defined by

[TABLE]

for each $\hat{f}\in C^{k}_{b}(\hat{K})$.

Definition 2.4**.**

The push-forward operation mapping the dual space $C^{k}_{b}(K)^{\prime}$ into $C^{k}_{b}(\hat{K})^{\prime}$ is defined by

[TABLE]

for each $n\in C^{k}_{b}(K)^{\prime}$.

It is easy to verify that the pull-back and push-forward are linear operations preserving the vector space operations. Moreover, they are invertible iff $F$ itself is. Therefore, we have

Proposition 2.1**.**

Given finite elements $(K,P,N)$ and $(\hat{K},\hat{P},\hat{N})$ such that $F(K)=\hat{K}$ and $F^{*}(\hat{P})=P$, $F^{*}:\hat{P}\rightarrow P$ and $F_{*}:P^{\prime}\rightarrow\hat{P}^{\prime}$ are isomorphisms.

The pull-back and push-forward operations are also defined over the vector spaces $X$, $X^{\dagger}$, $\hat{X}$, and $\hat{X}^{\dagger}$. If $N$ is a vector of functionals and $\Phi$ a vector of functions, then the vector push-forward and pull-back are, respectively

[TABLE]

It will also be useful to consider vectors of functionals acting on vectors of functions. We define this to produce a matrix as follows. If $N=\begin{bmatrix}n_{1}&n_{2}&\dots&n_{k}\end{bmatrix}^{T}$ is a collection of functionals and $\Phi=\begin{bmatrix}\phi_{1}&\phi_{2}&\dots&\phi_{\ell}\end{bmatrix}^{T}$ a collection of functions, then we define the (outer) product $N(\Phi)$ to be the $k\times\ell$ matrix

[TABLE]

For example, if $N$ is the vector of nodes of a finite element and $\Psi$ contains the nodal basis functions, then the Kronecker delta property is expressed as $N(\Psi)=I.$

If $M$ is a matrix of numbers of appropriate shape and $\Phi\in X$ members of a function space $P$, then $M\Phi$ is just defined by $(M\Phi)_{i}=\sum_{j=1}^{\nu}M_{ij}\Phi_{j},$ according to the usual rule for matrix-vector multiplication.

Lemma 2.1**.**

Let $N\in X^{\dagger}$ and $\Phi\in X$ and $M\in\mathbb{R}^{\nu\times\nu}$. Then

[TABLE]

Proof.

The proof is a simple calculation:

[TABLE]

∎

The relationship between pull-back and push-forward also leads to the vectorized relation

Lemma 2.2**.**

Let $N\in X^{\dagger}$ and $\hat{\Phi}\in\hat{X}$. Then

[TABLE]

Definition 2.5**.**

Let $(K,P,N)$ and $(\hat{K},\hat{P},\hat{N})$ be finite elements and $F$ an affine mapping on $K$. Then $(K,P,N)$ and $(\hat{K},\hat{P},\hat{N})$ are affine equivalent if

• •

$F(K)=\hat{K}$,

• •

The pullback maps $F^{*}(\hat{P})=P$ (in the sense of equality of vector spaces),

• •

$F_{*}(N)=\hat{N}$* (in the sense of equality of finite sets).*

Definition 2.6**.**

Let $(K,P,N)$ be a finite element of class $C^{k}$ and $\Psi\in X$ its nodal basis. The nodal interpolant $\mathcal{I}_{N}:C_{b}^{k}(K)\rightarrow P$ is defined by

[TABLE]

This interpolant plays a fundamental role in establishing approximation properties of finite elements via the Bramble-Hilbert Lemma [7, 14]. The homogeneity arguments in fact go through for the following generalized notion of element equivalence:

Definition 2.7**.**

Two finite elements $(K,P,N)$ and $(K,P,\tilde{N})$ are interpolation equivalent if $\mathcal{I}_{N}=\mathcal{I}_{\tilde{N}}$.

Definition 2.8**.**

If $(K,P,\tilde{N})$ is affine equivalent to $(\hat{K},\hat{P},\hat{N})$ and interpolation equivalent to $(K,P,N)$, then $(K,P,N)$ and $(\hat{K},\hat{P},\hat{N})$ are affine-interpolation equivalent.

Brenner and Scott [8] give the following result, of which we shall make use:

Proposition 2.2**.**

Finite elements $(K,P,N)$ and $(K,P,\tilde{N})$ are interpolation equivalent iff the spans of $N$ and $\tilde{N}$, (viewed as subsets of $C_{b}^{k}(K)^{\prime}$), are equal.

For Lagrange and certain other finite elements, one simply has that $F^{*}(\hat{\Psi})=\Psi$, which allows for the traditional use of reference elements used in FEniCS, Firedrake, and countless other codes. However, for many other elements this is not the case. It is our goal in this paper to give a general approach that expresses $\Psi$ as a linear transformation $M$ applied to $F^{*}(\hat{\Psi})$.

Before proceeding, we note that approximation theory for Argyris and other families without affine-interpolation equivalence can proceed by means of establishing the almost-affine property [10]. Such proofs can involve embedding the inequivalent element family into an equivalent one with the requisite approximation properties. For example, the Argyris element is proved almost-affine by comparison to the “type (5)” quintic Hermite element. Although we see definite computational consequences of affine-equivalence, affine-interpolation equivalence, and neither among our element families, we our approach to transforming inequivalent families does not make use of any almost-affine properties.

3 Transformation theory when $F^{*}(\hat{P})=P$

For now, we assume that the pull-back operation (3) appropriately converts the reference element function space into the physical function space and discuss the construction of nodal bases based on relationships between the reference nodes $\hat{N}$ and the pushed-forward physical nodes $F_{*}(N)$.

We focus on the simplicial case, although generalizations do not have a major effect, as we note later. Throughout, we will use following convention, developed in [32] for handling facet orientation in mixed methods but also useful in order higher-order Lagrange degrees of freedom. Since our examples are triangles (2-simplices), it is not necessary to expand on the entire convention. Given a triangle with vertices $\left(\mathbf{v}_{1},\mathbf{v}_{2},\mathbf{v}_{3}\right)$, we define edge $\gamma_{i}$ of the triangle to connect the vertices other than $\mathbf{v}_{i}$. The (unit) tangent vector $\mathbf{t}_{i}=\begin{bmatrix}t^{\mathbf{x}}_{i}&t^{\mathbf{y}}_{i}\end{bmatrix}^{T}$, points in the direction from the lower- to the higher-numbered vertex. When triangles share an edge, then, they agree on its orientation. The normal to an edge is defined by rotating the tangent by applying the matrix $R=\begin{bmatrix}0&1\\ -1&0\end{bmatrix}$ so that $\mathbf{n}_{i}=R\mathbf{t}_{i}=\begin{bmatrix}n^{\mathbf{x}}_{i}&n^{\mathbf{y}}_{i}\end{bmatrix}^{T}$ We also let $\mathbf{e}_{i}$ denote the midpoint of $\gamma_{i}$.

Now, we fix some notation for describing nodes. First, we define $\delta_{\mathbf{x}}$ acting on any continuous function by pointwise evaluation. That is:

[TABLE]

We let $\delta^{\mathbf{s}}_{\mathbf{x}}$ denote the directional derivative in direction $\mathbf{s}$ at a point $\mathbf{x}$, so that

[TABLE]

We use repeated superscripts to indicate higher-order derivatives, so that $\delta^{\mathbf{x}\mathbf{x}}_{\mathbf{x}}$ defines the second directional derivative along the $x$-axis at point $\mathbf{x}$.

It will also be convenient to use block notation, with a single symbol representing two or items. For example, the gradient notation

[TABLE]

gives the pair of functionals evaluating the Cartesian derivatives at a point $\mathbf{x}$. To denote a gradient in a different basis, we append the directions as superscripts so that

[TABLE]

contains the normal and tangential derivatives at a point $\mathbf{x}$.

Similarly, we let

[TABLE]

denote the vector of three functionals evaluating the unique (supposing sufficient smoothness) second partials at $\mathbf{x}$.

Let $\Psi=\{\psi_{i}\}_{i=1}^{\nu}$ be the nodal basis for a finite element $(K,P,N)$ and $\hat{\Psi}=\{\hat{\psi}_{i}\}_{i=1}^{\nu}$ that for a reference element $\left(\hat{K},\hat{P},\hat{N}\right)$. We also assume that $F(K)=\hat{K}$ and $F^{*}(\hat{P})=P$. Because the pull-back is invertible, it maps linearly independent sets to linearly independent sets. So, $F^{*}(\hat{\Psi})$ must also be a basis for $P$. There exists an invertible $\nu\times\nu$ matrix $M$ such that

[TABLE]

or equivalently, that each nodal basis function is some linear combination of the pull-backs of the reference nodal basis functions.

Our theory for transforming the basis functions (i.e. computing the matrix $M$) will work via duality – relating the matrix $M$ to how the nodes, or at least their restrictions to the finite-dimensional spaces, push forward.

It will be useful to define as an intermediate $\nu\times\nu$ matrix $B=F_{*}(N)(\hat{\Psi})$. Recall from (6) that its entries for $1\leq i,j\leq\nu$ are

[TABLE]

This matrix, having nodes only applied to members of $P$ is indifferent to restrictions and so $B=F_{*}(\pi N)(\hat{\Psi})$ as well.

Because of Proposition 2.1 and finite-dimensionality, the the nodal sets $\hat{\pi}\hat{N}$ and $F_{*}(\pi N)$ are both bases for $\hat{P}^{\prime}$, and so there exists an invertible $\nu\times\nu$ matrix $V$ such that

[TABLE]

Frequently, it may be easier to express the pushed-forward nodes as a linear combination of the reference nodes. In this case, one obtains the matrix $V^{-1}$. At any rate, the matrices $V$ and $M$ are closely related.

Theorem 3.1**.**

For finite elements $(K,P,N)$ and $(\hat{K},\hat{P},\hat{N})$ with $F(K)=\hat{K}$ and $F_{*}(\hat{P})=P$, the matrices in (12) and (14) satisfy

[TABLE]

Proof.

We proceed by relating both matrices to $B$ defined in (13) via the Kronecker property of nodal bases. First, we have

[TABLE]

so that

[TABLE]

Similarly,

[TABLE]

so that $V=B^{-1}$ and the result follows. ∎

That is, to relate the pullback of the reference element basis functions to any element’s basis functions, it is sufficient to determine the relationship between the nodes.

3.1 Affine equivalence: The Lagrange element

When elements form affine-equivalent families, the matrix $M$ has a particularly simple form.

Theorem 3.2**.**

If $(K,P,N)$ and $(\hat{K},\hat{P},\hat{N})$ are affine-equivalent finite elements then the transformation matrix $M$ is the identity.

Proof.

Suppose the two elements are affine-equivalent, so that $F_{*}(N)=\hat{N}$. Then, a direct calculation gives

[TABLE]

so that $M=I$. ∎

The Lagrange elements are the most widely used finite elements and form the prototypical affine-equivalent family [8]. For a simplex $K$ in dimension $d$ and integer $r\geq 1$, one defines $P=P_{r}(K)$ to be the space of polynomials over $K$ of total degree no greater than $r$, which has dimension $\binom{r+d}{d}$. The nodes are taken to be pointwise evaluation at a lattice of $\binom{r+d}{d}$ points. Classically, these are taken to be regular and equispaced, although options with superior interpolation and conditioning properties for large $r$ are also known [17]. One must ensure that nodal locations are chosen at the boundary to enable $C^{0}$ continuity between adjacent elements. A cubic Lagrange triangle ($r=3$ and $d=2$) is shown earlier in Figure 1(a).

The practical effect of Theorem 3.2 is that the reference element paradigm “works.” That is, a computer code contains a routine to evaluate the nodal basis $\hat{\Psi}$ and its derivatives for a reference element $(\hat{K},\hat{P},\hat{N})$. Then, this routine is called at a set of quadrature points in $\hat{K}$. One obtains values of the nodal basis at quadrature points on each cell $K$ by pull-back, so no additional work is required. To obtain the gradients of each basis function at each quadrature point, one simply multiplies each basis gradient at each point by $J^{T}$.

On the other hand, when $M\neq I$, the usage of tabulated reference values is more complex. Given a table

[TABLE]

of the reference basis at the reference quadrature points, one finds the nodal basis for $(K,P,N)$ by constructing $M$ for that element and then computing the matrix-vector product $M\hat{\varPsi}$ so that

[TABLE]

Mapping gradients from the reference element requires both multiplication by $M$ as well as application of $J^{T}$ by the chain rule. We define $D\hat{\varPsi}\in\mathbb{R}^{\nu\times|\xi|\times 2}$ by

[TABLE]

Then, the basis gradients requires contraction with $M$

[TABLE]

followed by the chain rule

[TABLE]

In fact, the application of $M$ and $J^{T}$ can be performed in either order. Note that applying $M$ requires an $\nu\times\nu$ matrix-vector multiplication and in principle couples all basis functions together, while applying $J^{T}$ works pointwise on each basis function separately. When $M$ is quite sparse, one expects this to be a small additional cost compared to the other required arithmetic. We present further details for this in the case of Hermite elements, to which we now turn.

3.2 The Hermite element: affine-interpolation equivalence

The Hermite triangle [11], show in Figure 1(b) is based cubic polynomials, although higher-order instances can also be defined [8]. In contrast to the Lagrange element, its node set includes function values and derivatives at the nodes, as well as an interior function value. The resulting finite element spaces have $C^{0}$ continuity with $C^{1}$ continuity at vertices. They provide a classic example of elements that are not affine equivalent but instead give affine-interpolation equivalent families.

We will let $(K,P,N)$ be a cubic Hermite triangle, specifying the gradient at each vertex in terms of the Cartesian derivatives – see Figure 3(b). Let $\{\mathbf{v}_{i}\}_{i=1}^{3}$ be the three vertices of $K$ and $\mathbf{v}_{4}$ its barycenter. We order the nodes $N$ by

[TABLE]

using block notation.

Now, we fix the reference element $(\hat{K},\hat{P},\hat{N})$ with $\hat{K}$ as the unit right triangle and express the gradient by the derivatives in the direction of the reference Cartesian coordinates, as in Figure 3(a). Let $\{\hat{\mathbf{v}}_{i}\}_{i=1}^{3}$ be the three vertices of $\hat{K}$ and $\hat{\mathbf{v}}_{4}$ its barycenter. We define $\hat{N}$ analogously to $N$.

Consider the relationship between the nodal basis functions $\Psi$ and the pulled-back $F^{*}(\hat{\Psi})$. For any $\hat{\psi}\in\hat{P}$, the chain rule leads to

[TABLE]

Now, suppose that $\hat{\psi}$ is a nodal basis function corresponding to evaluation at a vertex or the barycenter, so that $\delta_{\hat{\mathbf{v}}_{i}}\hat{\psi}=1$ for some $1\leq i\leq 4$, with the remaining reference nodes vanishing on $\hat{\psi}$. We compute that

[TABLE]

while $\delta_{\mathbf{v}_{j}}F^{*}(\hat{\psi})=0$ for $1\leq j\leq 4$ with $j\neq i$. Also, since the reference gradient of $\hat{\psi}$ vanishes at each vertex, (23) implies that the physical gradient of $F^{*}(\hat{\psi})$ must also vanish at each vertex. So, pulling back $\hat{\psi}$ gives the corresponding nodal basis function for $(K,P,N)$.

The situation changes for the derivative basis functions. Now take $\hat{\psi}$ to be the basis function with unit-valued derivative in, say, the $\hat{\mathbf{x}}$ direction at vertex $\hat{\mathbf{v}}_{i}$ and other degrees of freedom vanishing. Since it vanishes at each vertex and the barycenter of $\hat{K}$, $F^{*}(\hat{\psi})$ will vanish at each vertex and the barycenter of $K$. The reference gradient of $\hat{\psi}$ vanishes at the vertices other than $i$, so the physical gradient of its pullback must also vanish at the corresponding vertices of $K$. However, (23) shows that $\nabla(\hat{\psi}\circ F)$ will typically not yield $\begin{bmatrix}1&0\end{bmatrix}^{T}$ at $\mathbf{v}_{i}$. Consequently, the pull-backs of the reference derivative basis functions do not produce the physical basis functions.

Equivalently, we may express this failure in terms of the nodes – pushing forward $N$ does not yield $\hat{N}$. We demonstrate this pictorially in Figure 4, showing the images of the derivative nodes under push-forward do not correspond to the reference derivative nodes. Taking this view allows us to address the issue using Theorem 15.

This discussion using the chain rule can be summarized by the matrix-valued equation

[TABLE]

noting that the second, fourth, and sixth rows and columns of this matrix are blocks of two, and each “[math]” is taken to be the zero matrix of appropriate size. This is exactly the inverse of $V$ from Theorem 15.

In this case, the transformation $V$ is quite local – that is, only the push-forward of nodes at a given point are used to construct the reference nodes at the image of that point. This seems to be generally true for interpolation-equivalent elements, although functionals with broader support (e.g. integral moments over the cell or a facet thereof) would require a slight adaptation. We will see presently for Morley and Argyris elements that the transformation neeed not be block diagonal for elements without interpolation equivalence. At any rate, the following elementary observation from linear algebra suggests the sparsity of $V$:

Proposition 3.1**.**

Let $W$ be a vector space with sets of vectors $W_{1}=\{w^{1}_{i}\}_{i=1}^{m}\subset W$ and $W_{2}=\{w^{2}_{i}\}_{i=1}^{n}$. Suppose that $\mathrm{span}W_{1}\subset\mathrm{spanW_{2}}$ so that there exists a matrix $A\in\mathbb{R}^{m\times n}$ such that $w^{1}_{i}=\sum_{k=1}^{n}A_{ik}w^{2}_{k}$. If we further have that some $w^{1}_{i}\in\mathrm{span}\{w^{2}_{j}\}_{j\in\mathcal{J}}$ for some $\mathcal{J}\subset[1,n]$, then $A_{ij}=0$ for all $j\notin\mathcal{J}$.

Our theory applies equally to the general family of Hermite triangles of degree $k\geq 3$. In those cases, the nodes consist of gradients at vertices together with point-wise values at appropriate places. All higher-order cases generate $C^{0}$ families of elements with $C^{1}$-continuity at vertices. The $V$ matrix remains analogous to the cubic case, with $J^{-T}$ on the diagonal in three places corresponding to the vertex derivative nodes. No major differences appear for the tetradral Hermite elements, either.

As we saw earlier, Hermite and other elements for which $M\neq I$ incur an additional cost in mapping from the reference element, as one must compute basis function values and gradients via (18) and (21). The key driver of this additional cost is the application of $M$. Since $M$ is very sparse for Hermite elements – just 12 nonzeros counting the 1’s on the diagonal – evaluating (18) requires just $12$ operations per column, so a 10-point quadrature rule requires 120 operations. Evaluating (20) requires twice this, or 240 operations. Applying $J^{T}$ in (21) is required whether Hermite or Lagrange elements are used. It requires $4\times 10$ times the number of quadrature points used – so a 10-point rule would require 400 operations. Hence, the chain rule costs more than the application of $M$ in this situation. On the other hand, building an element stiffness matrix requires a double loop over these 10 basis functions nested with a loop over the, say, 10 quadrature points. Hence, the loop body requires 1000 iterations, and with even a handful of operations will easily dominate the additional cost of multiplying by $M$.

3.3 The Morley and Argyris elements

The construction of $C^{1}$ finite elements, required for problems such as plate bending or the Cahn-Hilliard equations, is a long-standing difficulty. Although it is possible to work around this requirement by rewriting the fourth-order problem as a lower order system or by using $C^{0}$ elements in conjunction with variational form penalizing the jumps in derivatives [15, 33], this doesn’t actually give a $C^{1}$ solution.

The quadratic Morley triangle [28], shown in Figure 1(c), finds application in plate-bending problems and also provides a relatively simple motivation for and application of the theory developed here. The six degrees of freedom, vertex values and the normal derivatives on each edge midpoint, lead to an assembled finite element space that is neither $C^{0}$ nor $C^{1}$, but it is still suitable as a convergent nonconforming approximation for fourth-order problems.

The quintic Argyris triangle [3], shown in Figure 1(d), with its 21 degrees, gives a proper $C^{1}$ finite element. Hence it can be used generically for fourth-order problems as well as second-order problems for which a continuously differentiable solution is desired. The Argyris elements use the values, gradients, and second derivatives at each triangle vertex plus the normal derivatives at edge midpoints as the twenty-one degrees of freedom.

It has been suggested that the Bell element [5] represents a simpler $C^{1}$ element than the Argyris element, on the account that it has fewer degrees of freedom. Shown in Figure 1(e), we see that the edge normal derivatives have been removed from the Argyris element. However, this comes with a (smaller but) more complicated function space. Rather than full quintic polynomials, the Bell element uses quintic polynomials that have normal derivatives on each edge of only third degree. This constraint on the polynomial space turns out to complicate the transformation of Bell elements compared to Hermite or even Argyris. For the rest of this section, we focus on Morley and Argyris, returning to Bell later.

It can readily be seen that, like the Hermite element, the standard affine mapping will not preserve nodal bases. Unlike the Hermite element, however, the Morley and Argyris elements do not form affine-interpolation equivalent families – the spans of the nodes are not preserved under push-forward thanks to the edge normal derivatives – see Figure 5. As the Morley and Aryris nodal sets do not contain a full gradient at edge midpoints, the technique used for Hermite elements cannot be directly applied.

To work around this, we introduce the following idea:

Definition 3.1**.**

Let $(K,P,N)$ and $(\hat{K},\hat{P},\hat{N})$ be finite elements of class $C^{k}$ with affine mapping $F:K\rightarrow\hat{K}$ and associated pull-back and push-forward $F^{*}$ and $F_{*}$. Suppose also that $F^{*}(\hat{P})=P$. Let $N^{c}=\left\{n^{c}_{i}\right\}_{i=1}^{\mu}\subset C_{b}^{k}(K)^{\prime}$ and $\hat{N}^{c}=\left\{\hat{n}^{n}_{i}\right\}_{i=1}^{\mu}\subset C^{k}(\hat{K})^{\prime}$ be such that

• •

$N\subset N^{c}$* (taken as sets rather than vectors),*

• •

$\hat{N}\subset\hat{N}^{c}$* (again as sets),*

• •

$\mathrm{span}(F_{*}(N^{c}))=\mathrm{span}(\hat{N}^{c})$* in $C^{k}(\hat{K})^{\prime}$.*

Then $N^{c}$ and $\hat{N}^{c}$ form a compatible nodal completion of $N$ and $\hat{N}$.

Example 3.1**.**

Let $(K,P,N)$ and $(\hat{K},\hat{P},\hat{N})$ be the Morley triangle and reference triangle. Take $N^{c}$ to contain all the nodes of $N$ together with the tangential derivatives at the midpoint of each edge of $K$ and similarly for $\hat{N}^{c}$. In this case, $\mu=9$. Then, both $N^{c}$ and $\hat{N}^{c}$ contain complete gradients at each edge midpoint and function values at each vertex. The push-forward of $N^{c}$ has the same span as $\hat{N}^{c}$ and so $N^{c}$ and $\hat{N}^{c}$ form a compatible nodal completion of $N$ and $\hat{N}$. This is shown pictorially in Figure 6.

A similar completion – supplementing the nodes with tangential derivatives at edge midpoints – exists for the Argyris nodes and reference nodes [29].

Now, since the spans of $\hat{N}^{c}$ and $F_{*}(N^{c})$ agree (even in $C_{b}^{k}(\hat{K})^{\prime}$), there exists a $\mu\times\mu$ matrix $V^{c}$, typically block diagonal, such that

[TABLE]

Let $E\in\mathbb{R}^{\nu\times\mu}$ be the Boolean matrix with $E_{ij}=1$ iff $\hat{n}_{i}=\hat{n}_{j}^{c}$ so that

[TABLE]

and it is clear that

[TABLE]

That is, the reference nodes are linear combinations of the pushed-forward nodes and the extended nodes, but we must have the linear combination in terms of the pushed-forward nodes alone.

Recall that building the nodal basis only requires the action of the nodes on the polynomial space. Because $\mu>\nu$, the set of nodes $\pi N^{c}$ must be linearly dependent. So, we seek a matrix $D\in\mathbb{R}^{\mu\times\nu}$ such that

[TABLE]

Since $F_{*}$ is an isomorphism, such a $D$ also gives

[TABLE]

Rows $i$ of the matrix $D$ such that $n^{c}_{i}=n_{j}$ for some $j$ will just have $D_{ik}=\delta_{kj}$ for $1\leq k\leq\nu$. The remaining rows must be constructed somehow via an interpolation argument, although the details will vary by element.

This discussion suggests a three-stage process, each encoded by matrix multiplication, for converting the push-forwards of the physical nodes to the reference nodes, hence giving a factored form of $V$ in (14). Before working examples, we summarize this in the following theorem:

Theorem 3.3**.**

Let $(K,P,N)$ and $(\hat{K},\hat{P},\hat{N})$ be finite elements with affine mapping $F:K\rightarrow\hat{K}$ and suppose that $F^{*}(\hat{P})=P$. Let $N^{c}$ and $\hat{N}^{c}$ be a compatible nodal completion of $N$ and $\hat{N}$. Then given matrices $E\in\mathbb{R}^{\nu\times\mu}$ from (26), $V^{c}\in\mathbb{R}^{\mu\times\mu}$ from (25) and $D\in\mathbb{R}^{\mu\times\nu}$ from (28) that builds the (restrictions of) the extended nodes out of the given physical nodes, the nodal transformation matrix $V$ satisfies

[TABLE]

This gives a general outline for mapping finite elements, and we illustrate now by turning to the Morley element.

3.3.1 The Morley element

Following our earlier notation for the geometry and nodes, we order the nodes of a Morley triangle by

[TABLE]

Nodes $N^{C}$ will also include tangential derivatives at the edge midpoint. We put

[TABLE]

Again, this is a block vector the last three entries each consist of two values. We give the same ordering of reference element nodes $\hat{N}$ and $\hat{N}^{c}$.

The matrix $E$ simply extracts the members of $N^{C}$ that are also in $N$, so with $\eta=\begin{bmatrix}1&0\end{bmatrix}$, we have the block matrix

[TABLE]

Because the gradient nodes in $N^{c}$ use normal and tangential coordinates, $V^{c}$ will be slightly more more complicated than $V$ for the Hermite element. For local edge $\gamma_{i}$, we define the (orthogonal) matrix

[TABLE]

with the normal and tangent vector in the rows. Similarly, we let

[TABLE]

contain the unit normal and tangent to edge $\hat{\gamma}_{i}$ of the reference cell $\hat{K}$. It is clear that

[TABLE]

so, defining

[TABLE]

we have that

[TABLE]

Now, we turn to the matrix $D\in\mathbb{R}^{9\times 6}$, writing members of $\pi N^{c}$ in terms of $\pi N$ alone. The challenge is to express the tangential derivative nodes in terms of the remaining six nodes – vertex values and normal derivatives. In fact, only the vertex values are needed. Along any edge, any member of $P$ is just a univariate quadratic polynomial, and so the tangential derivative is linear. Linear functions attain their average value over an interval at its midpoint. But the average value of the derivative over the edge is just the difference between vertex values divided by the edge length. The matrix $D$ must be

[TABLE]

We can also arrive at this formulation of $D$ in another way, that sets up the discussion used for Argyris and later Bell elements. Consider the following univariate result:

Proposition 3.2**.**

Let $p(x)$ any quadratic polynomial on $[-1,1]$. Then

[TABLE]

Proof.

Write $p(x)=a+bx+cx^{2}$. Then $p^{\prime}(x)=b+2cx$ so that $p^{\prime}(0)=b$. Also note that $p(1)=a+b+c$ and $p(-1)=a-b+c$. Wanting to write $p^{\prime}(0)=d_{1}p(1)+d_{-1}p(-1)$ for constants $d_{1}$ and $d_{-1}$ leads to a $2\times 2$ linear system, which is readily solved to give $d_{1}=-d_{-1}=\tfrac{1}{2}$. ∎

Then, by a change of variables, this rule can be mapped to $\left[-\tfrac{\ell}{2},\tfrac{\ell}{2}\right]$ so that

[TABLE]

Finally, one can apply this rule on the edge of a triangle running from $\mathbf{v}_{a}$ to $\mathbf{v}_{b}$ to find that

[TABLE]

It is interesting to explicitly compute the product $V=EV^{C}D$, as giving a single formula rather than product of matrices is more useful in practice. Multiplying through gives:

[TABLE]

From the definition of $B^{i}$, it is possibly to explicitly calculate its entries in terms of the those of the Jacobian and the normal and tangent vectors for $K$ and $\hat{K}$. Only the first row of each $B^{i}$ is needed

[TABLE]

We can also recall that the normal and tangent vectors are related by $n^{\mathbf{x}}=t^{\mathbf{y}}$ and $n^{\mathbf{y}}=-t^{\mathbf{x}}$ to express these entries purely in terms of either the normal or tangent vectors. Each entry of the Jacobian and normal and tangent vectors of $K$ and $\hat{K}$ enter into the transformation.

In this form, $V$ has 12 nonzero entries, although the formation of those entries, which depend on normal and tangent vectors and the Jacobian, from the vertex coordinates requires an additional amount of arithmetic. The Jacobian will typically be computed anyway in a typical code, and the cost of working with $M=V^{T}$ will again be subdominant to the nested loops over basis functions and quadrature points required to form element matrices, much like Hermite.

3.3.2 The Argyris element

Because it is higher degree than Morley and contains second derivatives among the nodes, the Argyris transformation is more involved. However, it is a prime motivating example and also demonstrates that the general theory here reproduces the specific technique in [29]. The classical Argyris element has $P$ as polynomials of degree 5 over a triangle $K$, a 21-dimensional space. The 21 associated nodes $N$ are selected as the point values, gradients, and all three unique second derivatives at the vertices together with the normal derivatives evaluated at edge midpoints. These nodal choices lead to a proper $C^{1}$ element, and $C^{2}$ continuity is obtained at vertices.

Since the Argyris elements do not form an affine-interpolation equivalent family, we will need to embed the physical nodes into a larger set. Much as with Morley elements, the edge normal derivatives will be augmented by the tangential derivatives.

With this notation, $N$ is a vector of 21 functionals and $N^{C}$ a vector of 24 functions written as

[TABLE]

with corresponding ordering of reference nodes $\hat{N}$ and $\hat{N}^{c}$. The $21\times 24$ matrix $E$ just selects out the items in $N^{C}$ that are also in $N$, so that

[TABLE]

The matrix $V^{C}$ relating the push-forward of the extended nodes to the extended reference nodes is block diagonal and similar to our earlier examples. We use (23) to map the vertex gradient nodes as in the Hermite case. Mapping the three unique second derivatives by the chain rule requires the matrix:

[TABLE]

The edge midpoint nodes transform by $B$ just as in (35), so that the $V^{C}$ is

[TABLE]

Constructing $D$, like for Morley, is slightly more delicate. The additional nodes acting on quintic polynomials – tangential derivatives at edge midpoints – must be written in terms of the remaining nodes. The first aspect of this involves a univariate interpolation-theoretic question. On the biunit interval $[-1,1]$, we seek a rule of the form

[TABLE]

that is exact when $f$ is a quintic polynomial. The coefficients may be determined to by writing a $6\times 6$ linear system asserting correctness on the monomial basis. The answer, given in [29], is that

Proposition 3.3**.**

Any quintic polynomial $p$ defined on $[-1,1]$ satisfies

[TABLE]

This can be mapped to the interval $[-\tfrac{\ell}{2},\tfrac{\ell}{2}]$ by a change of variables:

[TABLE]

Now, we can use this to compute the tangential derivative at an edge midpoint, expanding the tangential first and second derivatives in terms of the Cartesian derivatives. If $\mathbf{v}_{a}$ and $\mathbf{v}_{b}$ are the beginning and ending vertex of edge $\gamma_{i}$ with midpoint $\mathbf{e}_{i}$ and length $\ell_{i}$, we write the tangential derivative acting on quintics as

[TABLE]

For each edge $\gamma_{i}$, define the vector $\mathbf{\tau}_{i}$ by

[TABLE]

The end result is that

[TABLE]

If this transformation is kept in factored form, $D$ contains 57 nonzero entries and $V^{c}$ contains 54 nonzero entries. $E$ is just a Boolean matrix and its application requires copies. So, application of $M$ requires no more than $111$ floating-point operations, besides the cost of forming the entries themselves. While this is about ten times the cost of the Hermite transformation, it is for about twice the number of basis functions and still well-amortized over the cost of integration loops. Additionally, one can multiply out the product $EV^{c}D$ symbolically and find only 81 nonzero entries, which reduces the cost of multiplication accordingly.

3.4 Generalizations

3.4.1 Non-affine mappings

Non-affine geometric transformations, whether for simplicial or other element shapes, present no major complications to the theory. In this case, $K$ and $\hat{K}$ are related by a non-affine map, and $P$ is taken to be the image of $\hat{P}$ under pull-back

[TABLE]

although this space need not consist of polynomials for non-affine $F$. At any rate, one may define Hermite elements on curvilinear cells [10, 13]. In this case, the Jacobian matrix varies spatially so that each instance of $J^{T}$ in (24) must be replaced by the particular value of $J^{T}$ at each vertex.

3.4.2 Generalized pullbacks

Many vector-valued finite element spaces make use of pull-backs other than composition with affine maps. For example, the Raviart-Thomas and Nédélec elements use contravariant and covariant Piola maps, respectively. Because these preserve either normal or tangential components, one can put the nodal basis functions of a given element $(K,P,N)$ and reference element $(\hat{K},\hat{P},\hat{N})$ into one-to-one correspondence by means of the Piola transform, a fact used heavily in [32] possible. It would be straightforward to give a generalization of affine equivalence to equivalence under an arbitrary pull-back $F^{*}$, with push-forward defined in terms of $F^{*}$. In this case, the major structure of § 3.1 would be unchanged.

However, not all $H(\mathrm{div})$ elements form equivalent families under the contravariant Piola transform. For example, Mardal, Tai, and Winther [27] give an element that can be paired with discontinuous polynomials to give uniform inf-sup stability on a scale of spaces between $H(\mathrm{div})$ and $(H^{1})^{2}$, although it is $H^{1}$-nonconforming. The degrees of freedom include constant and linear moments of normal components on edges, which are preserved under Piola mapping. However, the nodes also include the constant moments of the tangential component on edges, which are not preserved under Piola transform. One could push-forward both the normal and tangential constant moments, then express them as a linear combination of the normal and tangential moments on the reference cell in a manner like (24). One could see the Mardal–Tai–Winther element as satisfying a kind of “Piola-interpolation equivalence” and readily adapt the techniques for Hermite elements,

3.5 A further note on computation

We have commented on the added cost of multiplying the set of basis functions by $M$ during local integration. It also also possible to apply the transformation in a different way that perhaps more fully leverages pre-existing computer routines. With this approach, $M$ can also be included in local matrix assembly by means means of a congruence transform acting on the “wrong” element matrix as follows.

Given a finite element $(K,P,N)$ with nodal basis $\Psi=\{\psi_{i}\}_{i=1}^{\nu}$ and bilinear form $a_{K}(\cdot,\cdot)$ over the domain $K$, we want to compute the matrix

[TABLE]

Suppose that a computer routine existed for evaluating $A^{K}$ via a reference mapping for affine-equivalent elements. That is, given the mapping $F:\hat{K}\rightarrow K$, this routine maps all integration to the reference domain $\hat{K}$ assuming that the integrand over $K$ is just the affine pull-back of something on $\hat{K}$. Consider the following computation:

[TABLE]

Now, this is just expressed in terms of the affine pullback of reference-element integrands and so could use the hypothesized computer routine. We then have

[TABLE]

or, more compactly,

[TABLE]

where $\tilde{A}^{K}$ is the matrix one would obtain by using the pull-back of the reference element nodal basis functions instead of the actual nodal basis for $(K,P,N)$. Hence, rather than applying $M$ invasively at each quadrature point, one may use existing code for local integration and pre- and post-multiply the resulting matrix by the basis transformation. In the case of Hermite, for example, applying $M$ to a vector costs 12 operations, so applying $M$ to all 10 columns of $\tilde{A}^{K}$ costs 120 operations, plus another 120 for the transpose. This adds 240 extra operations to the cost of building $\tilde{A}^{K}$, or just 2.4 extra FLOPs per entry of the matrix.

One may also apply this idea in a “matrix-free” context. Given a routine for applying $\tilde{A}^{K}$ to a vector, one may simply apply $M^{T}$ to the input vector, apply $\tilde{A}^{K}$ to the result, and post-multiply by $M$. Hence, one has the cost of muliplying by $\tilde{A}^{K}$ plus the cost of applying $M$ and its transpose to a single vector. In the case of Hermite, one has the cost of computing the “wrong” local matrix-vector product via an existing kernel plus 24 additional operations.

Finally, we comment on evaluating discrete functions over elements requiring such transforms. Discrete function evaluation is frequently required in matrix-free computation, nonlinear residual evaluation, and in bilinear form evaluation when a coefficient is expressed in a finite element space. Suppose one has on a local element $K$ a function expressed by

[TABLE]

where $c\in\mathbb{R}^{\nu}$ is the vector of coefficients and $\{\psi_{j}\}$ is the nodal basis for $(K,P,N)$. In terms of pulled-back reference basis functions, $u$ is given by

[TABLE]

which can also be written as

[TABLE]

Just as one can build element matrices by means of the “wrong” basis functions and a patch-up operation, one can also evaluate functions by transforming the coefficients and then using the standard pullback of the reference basis functions. Such observations may make incorporating nonstandard element transformations into existing code more practical.

4 What if $P\neq F^{*}(\hat{P})$?

The theory so far has been predicated on $F^{*}$ providing an isomorphism between the reference and physical function spaces. In certain cases, however, this fails. Our main motivation here is to transform the Bell element, a near-relative of the quintic Argyris element. In this case, one takes $P$ to be the subspace of $P_{5}$ that has cubic normal derivatives on edges rather than the typical quartic values. This reduction of $P$ by three dimensions is accompanied by removing the three edge normal derivatives at midpoints from $N$. In general, however, the pull-back $F^{*}(\hat{P})$ does not coincide with $P$. Instead of cubic normal derivatives on edges, $F^{*}(\hat{P})$ has reduced degree in some other direction corresponding to the image of the normal under affine mapping. The theory developed earlier can be extended somewhat to resolve this situation.

4.1 General theory: extending the finite element

Abstractly, one may view the Bell element or other spaces built by constraint as the intersection of the null spaces of a collection of functionals acting on some larger space as follows. Let $(K,P,N)$ be a finite element. Suppose that $P\subset\tilde{P}$ and that $\{\lambda_{i}\}_{i=1}^{\kappa}\subset\left(C^{k}_{b}\right)^{\prime}$ are linearly independent functionals that when acting on $\tilde{P}$ satisfy

[TABLE]

The following result is not difficult to prove:

Proposition 4.1**.**

Let $(K,P,N)$ be a finite element with $\cap_{i=1}^{\kappa}\mathrm{null}(\lambda_{i})=P\subset\tilde{P}$ as per (53). Similarly, let Let $(\hat{K},\hat{P},\hat{N})$ be a reference element with $\cap_{i=1}^{\kappa}\mathrm{null}(\hat{\lambda}_{i})=\hat{P}\subset\tilde{\hat{P}}$. Suppose that $\tilde{P}=F^{*}(\tilde{\hat{P}}).$ Then $P=F^{*}(\hat{P})$ iff

[TABLE]

In the case of the Bell element, the span condition (54) fails and so that the function space is not preserved under affine mapping. Consequently, the theory of the previous section predicated on this preservation does not directly apply. Instead, we proceed by making the following observation.

Proposition 4.2**.**

Let $(K,P,N)$ be a finite element with $P\subset\tilde{P}$ satisfying $P=\cap_{i=1}^{\kappa}\mathrm{null}(\lambda_{i})$ for linearly independent functionals $\{\lambda_{i}\}_{i=1}^{\kappa}$. Define

[TABLE]

to include the nodes of $N$ together with $L=\begin{bmatrix}\lambda_{1}&\lambda_{2}&\dots&\lambda_{\kappa}\end{bmatrix}^{T}$. Then $(K,\tilde{P},\tilde{N})$ is a finite element.

Proof.

Since we have a finite-dimensional function space, it remains to show that $\tilde{N}$ is linearly independent and hence spans $\tilde{P}^{\prime}$. Consider a linear combination in $\tilde{P}^{\prime}$

[TABLE]

Apply this linear combination to any $p\in P$ to find

[TABLE]

since $\lambda_{i}(p)=0$ for $p\in P$. Because $(K,P,N)$ is a finite element, the $n_{i}$ are linearly independent in $P^{\prime}$ so $c_{i}=0$ for $1\leq i\leq\nu$. Applying the same linear combination to any $\in\tilde{P}\backslash P$ then gives that $d_{i}=0$ since the constraint functionals are also linearly independent. ∎

Given a nodal basis $(K,\tilde{P},\tilde{N})$, it is easy to obtain one for $(K,P,N)$.

Proposition 4.3**.**

Let $(K,P,N)$, $\{\lambda_{i}\}_{i=1}^{\kappa}$, and $(K,\tilde{P},\tilde{N})$ be as in Proposition 4.2. Order the nodes in $\tilde{N}$ by $\tilde{N}=\begin{bmatrix}N\\ L\end{bmatrix}$ with $L_{i}=\lambda_{i}$ for $1\leq i\leq\kappa$. Let $\{\tilde{\psi}_{i}\}_{i=1}^{\nu+\kappa}$ be the nodal basis for $(K,\tilde{P},\tilde{N})$. Then $\{\tilde{\psi}_{i}\}_{i=1}^{\nu}$ is the nodal basis for $(K,P,N)$.

Proof.

Clearly, $n_{i}(\tilde{\psi}_{j})=\delta_{ij}$ for $1\leq i,j\leq\nu$ by the ordering of the nodes in $\tilde{N}$. Moreover, $\{\tilde{\psi}_{i}\}_{i=1}^{\nu}\subset P$ because $\lambda_{i}(\tilde{\psi}_{j})=0$ for each $1\leq i\leq\kappa$. ∎

4.2 The Bell element

So, we can obtain a nodal basis for the Bell element or others with similarly constrained function spaces by mapping the nodal basis for a slightly larger finite element and extracting a subset of the basis functions. Let $(K,P,N)$ and $(\hat{K},\hat{P},\hat{N})$ be the Bell elements over $K$ and reference cell $\hat{K}$.

Recall that the Legendre polynomial of degree $n$ is orthogonal to polynomials of degree $n-1$ or less. Let $\mathcal{L}^{n}$ be the Legendre polynomial of degree $n$ mapped from the biunit interval to edge $\gamma_{i}$ of $K$. Define a functional

[TABLE]

For any $p\in P_{5}(K)$, its normal derivative on edge $i$ is cubic iff $\lambda_{i}(p)=0$. So, the constraint functionals are given in $L=\begin{bmatrix}\lambda_{1}&\lambda_{2}&\lambda_{3}\end{bmatrix}^{T}$ and $\tilde{N}=\begin{bmatrix}N\\ L\end{bmatrix}$ as in Proposition 4.2. We define

[TABLE]

and hence $(\hat{K},\hat{P},\hat{N})$ as well as $\hat{L}$ and $\tilde{\hat{N}}$ in a similar way.

$P$ and $\hat{P}$ are the constrained spaces – quintic polynomials with cubic normal derivatives on edges, while $\tilde{P}$ and $\tilde{\hat{P}}$ are the spaces of full quintic polynomials over $K$ and $\hat{K}$, respectively. We must construct a nodal basis for $(\hat{K},\tilde{\hat{P}},\tilde{\hat{N}})$, map it to a nodal basis for $(K,\tilde{P},\tilde{N})$ by the techniques in Section 3, and then take the subset of basis functions corresponding to the Bell basis.

This is accomplished by specifying a compatible nodal extension of $\tilde{N}$ and $\tilde{\hat{N}}$ by including the edge moments of tangential derivatives against $\mathcal{L}^{4}$ with those of $\tilde{N}$ and $\tilde{\hat{N}}$. We define

[TABLE]

We must specify the $E$, $V^{c}$, and $D$ matrices for this extended set of finite element nodes. We focus first on $D$, needing to compute each $\lambda_{i}^{\prime}$ in terms of the remaining functionals. As with Morley and Argyris, we begin with univariate results.

The following is readily confirmed, for example, by noting the right-hand side is a quintic polynomial and computing values and first and second derivatives at $\pm 1$:

Proposition 4.4**.**

Let $p$ be any quintic polynomial on $[-1,1]$. Then

[TABLE]

The formula (58) can be differentiated and then integrated against $\mathcal{L}^{4}$ to show that

[TABLE]

Then, this can be mapped to a general interval $[\tfrac{-\ell}{2},\tfrac{\ell}{2}]$ by a simple change of variables:

[TABLE]

Now, we can use this to express the functionals $\lambda_{i}^{\prime}$ from (57) as linear combinations of the Bell nodes:

Proposition 4.5**.**

Let $K$ be a triangle and $\mathbf{v}_{a}$ and $\mathbf{v}_{b}$ are the beginning and ending vertex of edge $\gamma_{i}$ with length $\ell_{i}$. Let $p$ be any bivariate quintic polynomial over $K$ and $\lambda_{i}^{\prime}$ defined in (57). Then the restriction of $\lambda_{i}^{\prime}$ to bivariate quintic polynomials satisfies

[TABLE]

and hence

[TABLE]

Now, $V^{c}$ is quite similar to that for the Argyris element. There is a slight difference in the handling the edge nodes, for we have an integral moment instead of a point value and must account for the edge length accordingly. By converting between normal/tangent and Cartesian coordinates via the matrix $G_{i}$ and mapping to the reference element, we find that for any $p$,

[TABLE]

This calculation shows that $V^{C}$ for the Bell element is identical to (43) for Argyris, except with a geometric scaling of the $B$ matrices.

The extraction matrix $E$ for the extended Bell elements consisting of full quintics now is identical to that for Argyris. Then, when evaluating basis functions, one multiplies the affinely mapped set of basis values by $V^{T}$ and then takes only the first 18 entries to obtain the local Bell basis.

4.3 A remark on the Brezzi-Douglas-Fortin-Marini element

In [21], we describe a two-part process for computing the triangular Brezzi-Douglas-Fortin-Marini (BDFM) element [16], an $H(\mathrm{div})$ conforming finite element based on polynomials of degree $k$ with normal components constrained to have degree $k-1$. This is a reduction of the Brezzi-Douglas-Marini element [9] somewhat as Bell is of Argyris. However, as both elements form Piola-equivalent families, the transformation techniques developed here are not needed.

Like the Bell element, one can define constraint functionals (integral moments of normal components against the degree $k$ Legendre polynomial) for BDFM. In [21], we formed a basis for the intersection of the null spaces of these functionals by means of a singular value decomposition. A nodal basis for the BDFM space then followed by building and inverting a generalized Vandermonde matrix on the basis for this constrained space.

In light of Propositions 4.2 and 4.3, however, this process was rather inefficient. Instead, we could have merely extended the BDFM nodes by the constraint functionals, building and inverting a single Vandermonde-like matrix. If one takes the BDM edge degrees of freedom as moments of normal components against Legendre polynomials up to degree $(k-1)$ instead of pointwise normal values, then one can even build a basis for BDM that includes a a basis for BDFM as a proper subset.

5 Numerical results

Incorporation of these techniques into high-level software tools such as Firedrake is the subject of ongoing investigation. In the meantime, we provide some basic examples written in Python, with sparse matrix assemble and solvers using petsc4py [12].

5.1 Scaling degrees of freedom

Before considering the accuracy of the $L^{2}$ projection, achieved via the global mass matrix, we comment on the conditioning of the mass and other matrices when both derivative and point value degrees of freedom appear. The Hermite element is illustrative of the situation.

On a cell of typical diameter $h$, consider a basis function corresponding to the point value at a given vertex. Since the vertex basis function has a size of $\mathcal{O}(1)$ on a triangle of size $\mathcal{O}(h^{2})$, its $L^{2}$ norm should be $\mathcal{O}(h)$. Now, consider a basis function corresponding to a vertex derivative. Its derivative is now $\mathcal{O}(1)$ on the cell, so that the $H^{1}$ seminorm is $\mathcal{O}(h)$. Inverse inequalities suggest that the $L^{2}$ norm could then be as large as $\mathcal{O}(1)$. That is, the different kinds of nodes introduce multiple scales of basis function sizes under transformation, which manifests in ill-conditioning. Where one expects a mass matrix to have an $\mathcal{O}(1)$ condition number, one now obtains an $\mathcal{O}(h^{-2})$ condition number. This is observed even on a unit square mesh, in Figure 7. All condition numbers are computed by converting the PETSc mass matrix to a dense matrix and using LAPACK via scipy [20]

However, there is a simple solution. For the Hermite element, one can scale the derivative degrees of freedom locally by an “effective $h$”. All cells sharing a given vertex must agree on that $h$, which could be the average cell diameter among cells sharing a vertex. Scaling the nodes/basis functions (which amounts to multiplying $V$ on the right by a diagonal matrix with 1’s or $h$’s) removes the scale separation among basis functions and leads again to an $\mathcal{O}(1)$ condition number for mass matrices, also seen in Figure 7. From here, we will assume that all degrees of freedom are appropriately scaled to give $\mathcal{O}(1)$ conditioning for the mass matrix.

5.2 Accuracy of $L^{2}$ projection

Now, we demonstrate that optimal-order accuracy is obtained by performing $L^{2}$ projection of smooth functions into the Lagrange, Hermite, Morley, Argyris, and Bell finite element spaces. In each case we use an $N\times N$ mesh divided into right triangles. Defining $u(x,y)=\sin(\pi x)\sin(2\pi y)$ on $[0,1]^{2}$, we seek $u_{h}$ such that

[TABLE]

for each $v_{h}\in V_{h}$, where $V_{h}$ is one of the the finite element spaces. Predicted asymptotic convergence rates – third for Morley, fourth for Hermite and Lagrange, fifth for Bell, and sixth for Argyris, are observed in Figure 8.

Note that the Hermite and Lagrange elements have the same order of approximation, but the Lagrange element delivers a slightly lower error. This is to be expected, as the space spanned by cubic Hermite triangles is a proper subset of that spanned by Lagrange.

5.3 The Laplace operator

As a simple second-order elliptic operator, we consider the Dirichlet problem for the Laplace operator on the unit square $\Omega$:

[TABLE]

equipped with homogeneous Dirichlet boundary conditions $u=0$ on $\partial\Omega$.

We divide $\Omega$ into an $N\times N$ mesh of triangles and let $V_{h}$ be one of the Lagrange, Hermite, Argyris, or Bell finite element spaces, all of which are $H^{1}$-conforming, over this mesh. The Morley element is not a suitable $H^{1}$ nonconforming element, so we do not use it here. We then seek $u_{h}\in V_{h}$ such that

[TABLE]

for all $v_{h}\in V_{h}$.

Enforcing strong boundary conditions on elements with derivative degrees of freedom is delicate in general. However, with grid-aligned boundaries, it is less difficult. To force a function to be zero on a given boundary segment, we simply require the vertex values and all derivatives tangent to the edge vanish. This amounts to setting the $x$-derivatives on the top and bottom edges of the box and $y$-derivative on the left and right for Hermite, Argyris, and Bell elements. Dirichlet conditions for Lagrange are enforced in the standard way.

By the method of manufactured solutions, we select $f(x,y)=8\pi^{2}\sin(2\pi x)\sin(2\pi y)$ so that $u(x,y)=\sin(2\pi x)\sin(2\pi y)$. In Figure 9, we show the $L^{2}$ error in the computed solution for both element families. As the mesh is refined, both curves approach the expected order of convergence – fourth for Hermite and Lagrange, fifth for Bell, and sixth for Argyris. Again, the error for Lagrange is slightly smaller than for Hermite, albeit with more global degrees of freedom.

5.4 The clamped plate problem

We now turn to a fourth-order problem for which the Argyris and Bell elements provide conforming $H^{2}$ discretizations and Morley a suitable nonconforming one. Following [8], we take the bilinear form defined on $H^{2}(\Omega)$ to be

[TABLE]

where $0<\nu<1$ yields a coercive bilinear form for any closed subspace of $H^{2}$ that does not contain nontrivial linear polynomials. We fix $\nu=0.5$.

Then, we consider the variational problem

[TABLE]

posed over suitable subspaces of $H^{2}$. It is known [8] that solutions of (68) that lie in $H^{4}(\Omega)$ satisfy the biharmonic equation $\Delta^{2}u=f$ in an $L^{2}$ sense.

We consider the clamped plate problem, in which both the function value and outward normal derivative are set to vanish, which removes nontrivial linear polynomials from the space. Again, we use the method of manufactured solutions on the unit square to select $f(x,y)$ such that $u(x,y)=\left(x(1-x)y(1-y)\right)^{2}$, which satifies clamped boundary conditions. We solve this problem with Argyris and Bell elements, and then also use the nonconforming Morley element in the bilinear form. Again, expected orders of convergence are observed in Figure 10.

6 Conclusions

Many users have wondered why FEniCS, Firedrake, and most other high-level finite element tools lack the full array of triangular elements, including Argyris and Hermite. One answer is that fundamental mathematical aspects of mapping such elements have remained relatively poorly understood. This work demonstrates the challenges involved with mapping such elements from a reference cell, but also proposes a general paradigm for overcoming those challenges by embedding the nodes into a larger set that transforms more cleanly and using interpolation techniques to relate the additional nodes back to original ones. In the future, we hope to incorporate these techniques in FInAT (https://github.com/FInAT/FInAT), a successor project to FIAT that produces abstract syntax for finite element evaluation rather than flat tables of numerical values. TSFC [18] already relies on FInAT to enable sum-factorization of tensor-product bases. If FInAT can provide rules for evaluating the matrix $M$ in terms of local geometry on a per-finite element basis, then TSFC and other form compilers should be able to seamlessly (from the end-users’ perspective) generate code for many new kinds of finite elements.