A High-Order Lower-Triangular Pseudo-Mass Matrix for Explicit Time Advancement of hp Triangular Finite Element Methods

TL;DR

This paper introduces a high-order lower-triangular pseudo-mass matrix for explicit time advancement in hp triangular finite element methods, enabling efficient computations without losing spatial accuracy.

Contribution

It proposes a novel pseudo-mass matrix approach for triangular elements that maintains accuracy while improving computational efficiency in explicit time-stepping.

Findings

Pseudo-mass matrix enables efficient explicit time-stepping.

Maintains exact projection of functions in polynomial spaces.

Applicable to high-order triangular finite element methods.

Abstract

Explicit time advancement for continuous finite elements requires the inversion of a global mass matrix. For spectral element simulations on quadrilaterals and hexahedra, there is an accurate approximate mass matrix which is diagonal, making it computationally efficient for explicit simulations. In this article it is shown that for the standard space of polynomials used with triangular elements, denoted $\mathcal{T}(p)$ where $p$ is the degree of the space, there is no diagonal approximate mass matrix that permits accurate solutions. Accuracy is defined as giving an exact projection of functions in $\mathcal{T}(p-1)$. In light of this, a lower-triangular pseudo-mass matrix method is introduced and demonstrated for the space $\mathcal{T}(3)$. The pseudo-mass matrix and accompanying high-order basis allow for computationally efficient time-stepping techniques without sacrificing the accuracy of the spatial approximation for unstructured triangular meshes.

Figures

Tables

Table 1. Table 1: Errors and convergence rates of the p = 3 𝑝 3 p=3 lower triangular inversion.

$k$	1	2	3	4	5	6
$h$	${\left(\frac{1}{2}\right)}^2$	${\left(\frac{1}{2}\right)}^3$	${\left(\frac{1}{2}\right)}^4$	${\left(\frac{1}{2}\right)}^5$	${\left(\frac{1}{2}\right)}^6$	${\left(\frac{1}{2}\right)}^7$
error	$4.0e-3$	$5.3e-4$	$6.8e-5$	$8.6e-6$	$1.1e-7$	$1.3e-8$
log-log slope		$2.6$	$2.8$	$2.9$	$3.0$	$3.0$

Equations

$$\frac{\partial u(\mathbf{x},t)}{\partial t}=\mathcal{L}(u)\text{ for }\mathbf{x}\in \Omega \text{ and }0<t<T_{final},$$

$${\Omega }^h=\bigcup _{1\le k\le N_{el}}{\Omega }_k$$

$${\mathcal{J}}_k({\Omega }_{ref})={\Omega }_k.$$

$${\Omega }_{ref}=\left\{(r,s)|-1\le r,s\text{ and }r+s\le 0\right\}.$$

$$\mathcal{T}(p)=\text{span}\left(\left\{\left(r^n{s}^m\right)|0\le m,n\text{ and }m+n\le p\right\}\right)$$

$$\text{dim}\left(\mathcal{T}(p)\right)=\frac{1}{2}(p+1)(p+2).$$

$$V_p^h=V_p\bigcap H^1(\Omega ).$$

$${\int }_{\Omega }{v}^h\frac{\partial u^h(\mathbf{x},t)}{\partial t}d\Omega ={\int }_{\Omega }{v}^h\mathcal{L}(u^h(\mathbf{x},t))d\Omega \forall v^h\in V_p^h.$$

$$u^h(\mathbf{x},t)=\sum u_i^h(t){\varphi }_i(\mathbf{x})={\mathbf{U}}^T(t)\vec{\varphi },$$

$${\mathbf{V}}^T{\int }_{\Omega }\vec{\varphi }{\vec{\varphi }}^{T}d\Omega \frac{\partial \mathbf{U}(t)}{\partial t}={\mathbf{V}}^T{\int }_{\Omega }\vec{\varphi }\mathcal{L}\left({\vec{\varphi }}^T\mathbf{U}(t)\right)d\Omega$$

$$M={\int }_{\Omega }\vec{\varphi }{\vec{\varphi }}^{T}d\Omega .$$

$$M\frac{\partial }{\partial t}\mathbf{U}(t)=R(\mathbf{U}(t)).$$

$$\frac{d}{dt}\mathbf{U}(t)=M^{-1}R(\mathbf{U}(t)).$$

$$\xi =-1+2\left(\frac{1+r}{1-s}\right)\text{ and }\eta =s.$$

$$\begin{array}[]{c}\phi_{A}(\xi,\eta)\,=\,\left(\frac{1+\eta}{2}\right),\\ \phi_{B}(\xi,\eta)\,=\,\left(\frac{1-\xi}{2}\right)\left(\frac{1-\eta}{2}\right),\\ \phi_{C}(\xi,\eta)\,=\,\left(\frac{1+\xi}{2}\right)\left(\frac{1-\eta}{2}\right).\end{array}$$

$$\begin{array}[]{c c}\phi_{{\cal A},m}\,=\,\left(\frac{1+\xi}{2}\right)\left(\frac{1-\xi}{2}\right)P_{m-1}^{2,2}(\xi)\left(\frac{1-\eta}{2}\right)^{m+1}&\text{ for }1\leq m\leq p-1,\\ \phi_{{\cal B},m}\,=\,\left(\frac{1+\xi}{2}\right)\left(\frac{1-\eta}{2}\right)\left(\frac{1+\eta}{2}\right)P_{m-1}^{2,2}(\eta)&\text{ for }1\leq m\leq p-1,\\ \phi_{{\cal C},m}\,=\,(-1)^{m-1}\left(\frac{1-\xi}{2}\right)\left(\frac{1-\eta}{2}\right)\left(\frac{1+\eta}{2}\right)P_{m-1}^{2,2}(\eta)&\text{ for }1\leq m\leq p-1.\end{array}$$

$${\varphi }_{\mathcal{I},m,n}=\left(\frac{1-\xi }{2}\right)\left(\frac{1+\xi }{2}\right)P_m^{2,2}(\xi ){\left(\frac{1-\eta }{2}\right)}^{m+2}\left(\frac{1+\eta }{2}\right)P_n^{2m+5,2}(\eta )$$

$$j(m,n)=\frac{1}{2}(m+n)(m+n+1)+n+1.$$

$$G:=\left\{g_1,g_2,\dots ,g_N\right\}.$$

$$g_n=\sum _{k=1,k\ne n}^N{a}_k{g}_k$$

$$\hat{G}:=\left\{f{g}_1,f{g}_2,\dots ,f{g}_N\right\}.$$

$$f{g}_n=f\left(\sum _{k=1,k\ne n}^N{a}_k{g}_k\right)=\sum _{k=1,k\ne n}^N{a}_{k}f{g}_k.$$

$${\vec{u}}_{\psi }=M_{\psi }^{-1}{\int }_{{\Omega }_{ref}}f\vec{\psi }d\Omega .$$

$${\vec{u}}_{\psi }^T\vec{\psi }=f\text{ for all }f\in \text{span}\left(\vec{\psi }\right).$$

$${\vec{u}}_{\psi }={\underline{M}}^{-1}\left({\int }_{\Omega }\vec{\psi }fd\Omega \right)$$

$${\vec{u}}_{\psi }^T\vec{\psi }=f\text{ for all }f\in \mathcal{T}(k)\subseteq \mathcal{T}(p).$$

$$D\,=\,\left[\begin{array}[]{cccc}D_{1}&&&\\ &D_{2}&&\\ &&\ddots&\\ &&&D_{\text{dim}\left(\mathcal{T}(p)\right)}\end{array}\right].$$

$${\vec{u}}_{\varphi }^T\vec{\varphi }=f.$$

$$\vec{\psi }=T\vec{\varphi }.$$

$${\psi }_A={\varphi }_A+\sum _{i=1}^{p-1}{A}_{\mathcal{B},i}{\varphi }_{\mathcal{B},i}+\sum _{i=1}^{p-1}{A}_{\mathcal{C},i}{\varphi }_{\mathcal{C},i}+\sum _{j=1}^{\frac{1}{2}(p-1)(p-2)}{A}_{\mathcal{I},j}{\varphi }_{\mathcal{I},j},$$

Full text

{AMS}

65M60, 65D30, 65M20, 65N35, 65N40

1 Introduction

The two-dimensional spectral element method (SEM) discussed in [19] is a continuous Galerkin finite element method (FEM) for quadrilateral meshes that utilizes a high degree polynomial basis. It is often coupled with explicit finite differences in time for the simulation of unsteady problems. Explicit simulations are computationally efficient for many problems and also easy to implement in parallel [26, 20]. For FEM, explicit time advancement schemes typically require the inversion of a global mass matrix at each time step, which can be computationally expensive [1]. This expense is avoided in the SEM because it uses Gauss-Lobatto (GL) integration, which provides the numerical integration scheme, the nodes (the GL points), and the basis functions (Lagrange interpolants over the GL points) for the method [22]. This approach constitutes a nodal collocation FEM because each basis function is only nonzero at one GL integration point. When the SEM mass matrix is evaluated using GL numerical integration, it becomes diagonal, which eliminates the expense of the mass matrix inversion. Note that with exact integration, the mass matrix is not diagonal. The diagonality is a consequence of the numerical integration; there is no polynomial basis with continuity properties appropriate for continuous FEM that is exactly orthogonal (i.e. a basis with a diagonal mass matrix). The SEM approximate mass matrix, which we call a pseudo-mass matrix, has the property that when used to project a function into the SEM approximation space it is exact if the function to be projected is a polynomial in the SEM space of one less degree than the approximation space [16].

Although the SEM has become a popular numerical technique for simulating challenging problems in complicated domains [12, 5], the need for high-order methods on unstructured meshes with robust adaptivity has motivated the development of a triangular finite element method (TFEM) that can compare to the SEM [23, 2, 18]. A comparable TFEM would allow an implementation of numerical methods for differential equations that possesses high-order spatial accuracy, low-memory usage, and efficient parallelizability. Such a method would lend itself well to structural mechanics [19], direct numerical simulations of computational fluid dynamics [19, 16, 24, 6], atmospheric modeling [14, 11], etc. [25, 28, 21, 22]. This has driven research into high-order continuous Galerkin TFEM possessing diagonal mass matrices that maintain a high level of accuracy [15, 4, 27, 8, 17]. In [16], Helenbrook showed that there is no nodal $C^{0}$ TFEM comparable to the SEM; specifically, he showed that there is no GL integration rule for triangles that can be used to create a $p^{th}$ degree nodal basis and is also accurate to order $2p-1$, which is the case on quadrilaterals. An artifact of the proof was the derivation a unique set of modal vertex functions that provide an accurate diagonal vertex block of a pseudo-mass matrix. The existence of these modal vertex functions inspire this investigation into an entire modal basis with a diagonal pseudo-mass matrix.

To investigate whether this basis and diagonal pseudo-mass matrix pair exists, the paper begins with a derivation of an explicit $C^{0}$ TFEM for an arbitrary basis and introduces the continuity requirements for a high-order basis. The modified Dubiner basis [9] is presented as a popular choice that satisfies these requirements. The concept of a pseudo-mass matrix and its accuracy requirements are then mathematically defined. A change of basis is introduced that can be used to map the Dubiner basis to any other high-order basis suitable for $C^{0}$ TFEM thus allowing the accuracy requirements to be defined for all bases. Then, these requirements are used to show that on triangles there is no basis and diagonal pseudo-mass matrix pair that satisfies the accuracy constraints.

To provide an alternative approach, we relax the diagonality requirement and instead look for a pseudo-mass matrix that can be inverted with only local operations. To this end, a new lower-triangular pseudo-mass matrix for $p=3$ is introduced that provides a desirable alternative to the full mass matrix approach by avoiding the need for a global mass matrix. This new method not only serves as a demonstration of concept, but is both a viable higher-order option and a segue into future work toward arbitrarily high-order continuous methods for triangles.

2 A High-Order Continuous Explicit TFEM

This section provides the formulation of a high-order $C^{0}$ TFEM for a transient problem, thereby introducing the necessary concepts and definitions for the following sections. The primary concern of an explicit time advancement method is the handling of the transient term. Of particular interest in this paper are transient partial differential equations (PDEs) that can be separated into a spatial operator ${\cal L}$ and a temporal derivative $\frac{\partial}{\partial t}$:

[TABLE]

in which $u$ is a scalar function, $\Omega\subset\mathbb{R}^{2}$ is a closed and bounded spatial domain, and $T_{final}$ is the end time of the simulation. As our analysis is entirely about the treatment of the unsteady term, the boundary conditions and specific form of ${\cal L}$ are unimportant.

The TFEM is based on a $N_{el}$ element partition $\Omega^{h}$ of the domain $\Omega$,

[TABLE]

in which $\Omega_{k}$ for all $k$, $1\leq k\leq N_{el}$ are affine transformations, $\left\{\mathcal{J}_{k}\right\}_{k=1}^{N_{el}}$, of a reference triangular element

[TABLE]

The typical triangular reference element is given by

[TABLE]

Locally defined polynomial approximation spaces, denoted $\mathcal{T}_{k}(p)$ for $1\leq k\leq N_{el}$, are given by mapping

[TABLE]

to the partition elements $\Omega_{k}$. The dimension of $\mathcal{T}(p)$ is

[TABLE]

In two dimensional PDEs with homogeneous boundary conditions, a continuous finite element method seeks a piecewise polynomial solution restricted to $H^{1}(\Omega)$, a subset of $C^{0}(\Omega)$ by the Sobolev Embedding Theorem. Let $V_{p}$ be the finite element space over $\Omega^{h}$ defined by the local polynomial spaces $\bigcup\limits_{k=1}^{N_{el}}\mathcal{T}_{k}(p)$. Define $V^{h}_{p}$ as

[TABLE]

then the finite element problem is to find $u^{h}({\bf x},t)\in V^{h}_{p}$, such that

[TABLE]

The space $V^{h}_{p}$ is finite dimensional due to the dimensionality of the local approximation spaces, and therefore $u^{h}$ may be expressed as a linear combination of global basis functions, $\vec{\phi}$, by

[TABLE]

and $v^{h}$ allows a similar form except with ${\bf U}^{T}$ replaced by ${\bf V}^{T}$. Eq. (8) can then be written as

[TABLE]

The global mass matrix, denoted $M$, is defined in the classical way as

[TABLE]

The result is a Galerkin TFEM

[TABLE]

where $R({\bf U}(t))$ is the discrete form of the spatial operator and the ${\bf V}^{T}$ term has been removed, as equation (10) must be true for all ${\bf V}^{T}$.

The model problem (1) has thus become a first order differential equation in time:

[TABLE]

In an explicit method the approximate solution at each time step is determined by a matrix multiplication of the inverse mass matrix $M^{-1}$. To be computationally efficient, it is therefore necessary to have a mass matrix that is trivially inverted, i.e. a diagonal mass matrix.

2.1 High Order Bases for Continuous Methods

Following along with the work of Dubiner in [9], we present necessary conditions for the direct enforcement of a globally $C^{0}$ approximation space. Specifically, the reference element $\Omega_{ref}$ is further decomposed into the vertices $V_{A}$, $V_{B}$, $V_{C}$, and the edges $E_{\mathcal{A}}$, $E_{\mathcal{B}}$, and $E_{\mathcal{C}}$ as shown in Fig. 1.

The concept of vertex, edge, and interior spaces is used to categorize local basis functions: A vertex function is a function $f:\Omega_{ref}\rightarrow\mathbb{R}$ with non-zero evaluation, typically $1$, at only one of the vertices and zero evaluation along the opposing edge. An edge function is a function $f:\Omega_{ref}\rightarrow\mathbb{R}$ with non-zero evaluation along only one edge and zero evaluation along the other two edges. And, an interior function is a non-zero function $f:\Omega_{ref}\rightarrow\mathbb{R}$ such that the evaluation of $f$ is zero along all edges (and at all vertices). For a basis for the space ${\cal T}(p)$ to be useful for $C^{0}$ finite elements, there must be one vertex function for each vertex, $(p-1)$ edge functions for each edge, and $(p-1)(p-2)/2$ interior functions.

In [9], Dubiner introduced a basis, $\vec{\phi}$ for $\mathcal{T}(p)$ that is commonly used for $C^{0}$ TFEM. This basis is achieved by mapping $\Omega_{ref}$ to the reference quadrilateral element, $\{(\xi,\eta)|-1\leq\xi,\eta\leq 1\}$, by

[TABLE]

The reference triangle can then be represented in terms of a warped coordinate system $(\xi,\eta)$, see Fig. 2.

In terms of the warped coordinates $(\xi,\eta)$, the modified Dubiner basis functions are originally defined in [9], then slightly modified in [16]. The three vertex functions are given by

[TABLE]

For $p\geq 2$ the basis function definitions rely on the classical Jacobi Polynomials $P^{\alpha,\beta}_{n}(z)$, these polynomials are orthogonal over the interval $(-1,1)$ with respect to the weighting function $(1-z)^{\alpha}(1+z)^{\beta}$. The $(p-1)$ edge functions for each edge are

[TABLE]

The interior functions for $0\leq m\leq(p-3)$ and $0\leq n\leq(p-3-m)$ are defined

[TABLE]

and given one-dimensional indices using

[TABLE]

An important property of this basis, and linearly independent sets of polynomials in general, regarding subsets and their linear independence is given in the following lemma which states specifically that the linear independence of a set of polynomials is not affected by the removal of common factors. This property will later be used to show the invertibility of a matrix with elements defined by the $L_{2}$ inner product.

Lemma 2.1**.**

Any finite dimensional, linearly independent set of polynomials possessing a common factor $f$, $\{fg_{i}\}$ defines the linearly independent set, $\{g_{i}\}$.

Proof 2.2**.**

To prove this consider the contrapositive: The linear dependence of a set of polynomials is maintained after element-wise multiplication of any function. Let $G$ be a set $N$ of linearly dependent functions over some space $\Omega$,

[TABLE]

Then there exists some $n$, $1\leq n\leq N$, for which $g_{n}\in G$ is a linear combination of the other functions:

[TABLE]

Let $f$ be any function defined on $\Omega$. Define the set of functions $\hat{G}$,

[TABLE]

It is then clear that

[TABLE]

*Therefore, $\hat{G}$ is linearly dependent as well. *

As an example application of this lemma, the set of edge $E_{{\cal B}}$ basis polynomials for $\mathcal{T}(p)$ all possess common factors of $\phi_{A}$ and $\phi_{C}$, the vertex $V_{A}$ and $V_{C}$ functions. The factor $\phi_{A}$ ensures that the functions $\phi_{{\cal B},i}$ for $1\leq i\leq p-1$ are all zero along $E_{{\cal A}}$, and the factor $\phi_{C}$ forces a zero evaluation along $E_{\cal C}$. The functions $\phi_{{\cal B},i}$ are polynomials of degree varying between $2$ and $p$, whereas $\phi_{A}$ and $\phi_{B}$ are first degree polynomials. By factoring $\phi_{A}$ and $\phi_{C}$ out of the edge functions, $\phi_{\cal B}/\left(\phi_{A}\phi_{C}\right)$, the resulting set of polynomials remains linearly independent with degree ranging between [math] and $p-2$. In fact, the resulting set is a linearly independent subset of $\mathcal{T}(p-2)$.

3 Diagonal pseudo-mass matrix

The $L_{2}$ projection of a function $f$ in $L_{2}(\Omega_{ref})$ onto $\mathcal{T}(p)$ is the unique function in $\mathcal{T}(p)$ that minimizes the error measured in the $L_{2}$ norm. In finite element methods the $L_{2}$ projection is evaluated by first identifying the coefficient vector of the projection and then taking the inner product of this coefficient vector with the basis. For any basis $\vec{\psi}$ of $\mathcal{T}(p)$ over $\Omega_{ref}$ there exists an invertible mass matrix $M_{\psi}$ defined by (11). The coefficient vector of the $L_{2}$ projection is

[TABLE]

The unique $L_{2}$ projection is then given by $\vec{u}_{\psi}^{T}\vec{\psi}$, and for all $f\in\mathcal{T}(p)$ the $L_{2}$ projection of $f$ is identically $f$:

[TABLE]

There is no basis that is both orthogonal and divisible into vertex, edge, and interior functions. This implies that there is no basis for $C^{0}$ FEM with a diagonal mass matrix. There are however pseudo-mass matrix approaches; a pseudo-mass matrix is simply a non-singular real-valued matrix $\underline{M}\in\mathbb{R}^{\text{dim}(\mathcal{T}(p))\times\text{dim}(\mathcal{T}(p))}$ that acts in place of the mass matrix. Mass-lumping techniques, which replace the mass matrix with a row-summed diagonal matrix, clearly constitute pseudo-mass matrix approaches [3, 8, 10, 13, 7]. However, these approaches often suffer from poor accuracy.

SEM is a pseudo-mass approach as the approximate mass matrix is cleverly generated by the GL integration scheme at the expense of only one polynomial degree of accuracy. In hopes of matching this success, we seek a high-order method that uses an accurate diagonal pseudo-mass matrix. The $p^{th}$ degree SEM pseudo-mass matrix is capable of exactly projecting all polynomials of one less polynomial degree than the basis, namely polynomials of degree $(p-1)$. We describe the accuracy of the pseudo-mass matrix by the degree of the subspace for which the projections are exact. For a FEM using a $p$ degree basis, the pseudo-mass matrix $\underline{M}$ is said to be $k$-exact for $1\leq k\leq p$ when the coefficient vector

[TABLE]

is exact for all $f$ of polynomial degree up to $k$; for TFEM

[TABLE]

The SEM uses a $(p-1)$-exact pseudo-mass matrix and it would be ideal for the TFEM to have a basis $\vec{\psi}$ for which there exists a $(p-1)$-exact diagonal pseudo-mass matrix, $D$:

[TABLE]

3.1 An Arbitrary Change of Basis for Continuous Methods

If a $(p-1)$-exact diagonal pseudo-mass matrix exists, the basis for which it functions is unknown. Therefore, it is beneficial to identify all sets of bases appropriate for the $C^{0}$ TFEM. Recall the modified Dubiner basis $\vec{\phi}$, (15) – (17). As $\vec{\phi}$ is a basis for $\mathcal{T}(p)$, for any function $f\in\mathcal{T}(p)$, by (20) $\vec{u}_{\phi}$ defined via (19) satisfies

[TABLE]

Furthermore by the equivalence of bases, any basis $\vec{\psi}$ may be written as a change of basis $T\in\mathbb{R}^{\text{dim}(\mathcal{T}(p))\times\text{dim}(\mathcal{T}(p))}$ of the modified Dubiner basis $\vec{\phi}$:

[TABLE]

Any basis $\vec{\psi}$ need only be defined by an invertible $T$ coupled with $\vec{\phi}$. However, an arbitrary invertible $T$ does not guarantee that the basis $\vec{\psi}$ satisfies the continuity constraints. Fortunately, the modified Dubiner basis is a continuous basis. Since the modified Dubiner basis satisfies the continuity constraints, the mapping $T$ must only enforce the preservation of these constraints in the new basis. The continuity constraints given earlier, §2.1, state the need of vertex, edge, and interior basis functions.

Let the vertex $V_{A}$, vertex $V_{B}$, and vertex $V_{C}$ basis functions be represented by the subscripts $A$, $B$, and $C$ respectively. For example the vertex $V_{A}$ function of $\vec{\phi}$ is written $\phi_{A}$. There will be $(p-1)$ functions along each edge, and these will be denoted by the subscripts ${\cal A}_{i}$, ${\cal B}_{i}$, and ${\cal C}_{i}$ for $1\leq i\leq(p-1)$. The functions $\phi_{{\cal A}_{i}}$, $\phi_{{\cal B}_{i}}$, and $\phi_{{\cal C}_{i}}$ will be non-zero on the edges $E_{\cal A}$, $E_{\cal B}$, and $E_{\cal C}$ respectively (shown in Fig. 1). Finally, ${\cal I}_{i}$ will be used to specify the interior functions, organized by $i=j(m,n)$, (18).

The vertex functions of $\vec{\phi}$ must remain vertex functions under transformation by $T$; the first three rows of $T$ define the three vertex functions of the new basis. Recall that any vertex $V_{A}$ function, $\psi_{A}$ must be zero along the edge $E_{\cal A}$. Therefore $\psi_{A}$ exists in the subspace of $\mathcal{T}(p)$ spanned by basis functions of $\vec{\phi}$ that are zero on $E_{\cal A}$: $\phi_{A}$, $\phi_{{\cal B},i}$, $\phi_{{\cal C},i}$, and $\phi_{{\cal I},i}$. Any new vertex $V_{A}$ function $\psi_{A}$ is then a linear combination of these functions:

[TABLE]

in which $\psi_{A}$ has been normalized by $\phi_{A}$ and the coefficients $A_{{\cal B},i}$, $A_{{\cal C},i}$, and $A_{{\cal I},j}$ represent the additions of edge $E_{\cal B}$, edge $E_{\cal C}$, and interior functions to $\phi_{A}$ respectively. The first three rows of $T$ are then given by

[TABLE]

for which $A_{\cal B}$, $A_{\cal C}$, $B_{\cal A}$, $B_{\cal C}$, $C_{\cal A}$, and $C_{\cal B}$ are all scalar valued row vectors of length $(p-1)$, and the interior augmentations are row vectors of length $(p-1)(p-2)/2$.

The next $3(p-1)$ rows of $T$ are used to map the edge functions of $\vec{\phi}$. The modified Dubiner basis edge functions are organized first by edge and then by degree such that basis functions $4$ to $4+(p-1)$ are edge $E_{\cal A}$ functions, followed by the edge $E_{\cal B}$ functions, then $E_{\cal C}$ functions. Note that this organization of basis functions is used in this section to place emphasis on the structure of individual edge subspaces of $\mathcal{T}(p)$, it will be altered in a later section. By definition, an edge function exists in a space spanned by functions that are non-zero along only one edge. The same-edge edge functions and the interior functions span this subspace. Similarly, only the interior functions satisfy the interior function requirements. Therefore, the most general transformation from the modified Dubiner basis to any other basis adhering to the continuity constraints is given by

[TABLE]

in which ${\cal A}_{\cal A}$, ${\cal B}_{\cal B}$, ${\cal C}_{\cal C}$, and ${\cal A}_{\cal I}$, ${\cal B}_{\cal I}$, and ${\cal C}_{\cal I}$ denote the additions of same-edge and interior functions. The blocks ${\cal A}_{\cal A}$, ${\cal B}_{\cal B}$, and ${\cal C}_{\cal C}$ are all $(p-1)\times(p-1)$ non-singular matrices. The dimensions of ${\cal A}_{\cal I}$, ${\cal B}_{\cal I}$, and ${\cal C}_{\cal I}$ are all $(p-1)\times(p-1)(p-2)/2$; and the $(p-1)(p-2)/2\times(p-1)(p-2)/2$ non-singular block ${\cal I}$ is used for the interior change of basis. As all diagonal blocks of $T$ are non-singular, $T$ is invertible. Thus, any basis $\vec{\psi}$ satisfying the continuity constraints can be represented by $\vec{\phi}$ via $T\vec{\phi}$ with $T$ restricted to the structure of (24).

3.2 No $(p-1)$-exact Diagonal TFEM

From [16], there are unique vertex functions for $p\geq 1$ that can be coupled with a $(p-1)$-exact diagonal pseudo-mass matrix. Therefore the first three rows of the transformation matrix exist. It remains to be determined if edge functions and interior functions can be found. Assume that the entire $(p-1)$-exact pseudo-mass matrix $D$ and basis $\vec{\psi}$ exist. So, for any function $f\in\mathcal{T}(p-1)$ there exists a unique $\vec{u}_{\psi}$ defined by

[TABLE]

such that $\vec{u}_{\psi}^{T}\vec{\psi}\,=\,f$. This differs from (19) because $D$ is not the mass matrix, and it is not exact for functions in $\mathcal{T}(p)/\mathcal{T}(p-1)$.

As the three vertex functions have been identified by Helenbrook [16], we look more closely to the first edge function $\psi_{{\cal A},1}$. Associated with the fourth row of (25), this edge function is defined in (22) by the following linear combination of modified Dubiner basis functions:

[TABLE]

By assumption, the fourth row of (25) should hold for all functions $f\in\mathcal{T}(p-1)$. These constraints will be identified by means of the modified Dubiner basis.

Let $f\in\mathcal{T}(p-1)$ and let $\vec{\psi}$ be a basis for which a $(p-1)$-exact diagonal pseudo-mass matrix exists. Note that because $f\in\mathcal{T}(p-1)\subset\mathcal{T}(p)$, there is a unique coefficient vector defined in (19) for which

[TABLE]

From (22), $\vec{\psi}=T\vec{\phi}$, and $\vec{u}_{\psi}^{T}\left(T\vec{\phi}\right)\,=\,\vec{u}_{\phi}^{T}\vec{\phi},$ or

[TABLE]

As the coefficient vectors are unique, $T^{T}\vec{u}_{\psi}\,=\,\vec{u}_{\phi}.$ The coefficient vector for $\vec{\psi}$ is then

[TABLE]

The inverse of $T$ can be found by block inversion of (24):

[TABLE]

with

[TABLE]

[TABLE]

and

[TABLE]

The constraints will consist of the $L_{2}$ projections of the basis functions of $\mathcal{T}(p-1)$. As the modified Dubiner basis is hierarchical, $f$ in (25) can be replaced by $\phi_{i}$ for all indices $i$ such that $\phi_{i}\in\mathcal{T}(p-1)$. This in turn implies that $\vec{u}_{\phi}$ is zero at every entry except for the $i^{th}$ entry which is one. Substituting (26) and (27) on the right and left hand side of (25) respectively, the fourth row of (25) can be written as

[TABLE]

Equation (29) is a system of equations and may be more easily solved by isolating the unknowns from ${\cal A}_{\cal A}$ and ${\cal A}_{\cal I}$ on the right hand side of the equation. Denote these unknowns as the vector $\vec{T}^{*}_{4}$ defined as

[TABLE]

Using (30), (29) is

[TABLE]

for all $i$ such that $\phi_{i}\in\mathcal{T}(p-1)$. In §3.2.1 we formally state, and reduce, a list of indices $i$ of constraining functions. Then in §3.2.2 we prove that this system results in $\vec{T}^{*}_{4}$ being identically zero which contradicts the fact that $T$ must be invertible.

3.2.1 Defining the Constraining Indices

The constraints are all of the functions of one less polynomial degree, $(p-1)$. Recall that the indices are currently arranged as vertices, $[A,B,C]$, edge ${\cal A}$, ${\cal B}$, and ${\cal C}$; then finally the interior functions, ${\cal I}$. The edge and interior functions are organized by increasing order. For $p\geq 2$, which is the case of interest here, vertex functions are always part of the constraints, which implies that the indices $[1,2,3]$, representing $[A,B,C]$, are included. For $p>2$ all edge functions of $\vec{\phi}$ exist in $\mathcal{T}(p-1)$ except for the last function on each side. Thus the indices $\left[4\!:\!(p+1),(p+3)\!:\!2p,(2p+2)\!:\!(3p-1)\right]$ are included. Lastly, for $p>3$ all interior functions of $\mathcal{T}(p-1)$ are also a component of $\vec{\phi}$. The indexing system, (18), for the interior functions was created such that as $p$ increases, functions are added to the end of the indexing system, thus the interior constraints are the first $(p-2)(p-3)/2$ interior basis functions. This results in the following list of indices

[TABLE]

The inverse of a matrix transposed is the transpose of the inverse of that matrix; rewriting $T^{-T}_{4,i}$ as $T^{-1}_{i,4}$ in (31) gives

[TABLE]

The desired basis $\vec{\psi}$ is unknown, and therefore $T$ is unknown as well. However, a complete knowledge of the fourth column of $T^{-1}$ is unnecessary in the derivation of a contradiction. From (28), $D_{4,4}T^{-1}_{i,4}$ has only $(p+1)$ non-zero entries, namely one from the term $-B_{\cal A}{\cal A}_{\cal A}^{-1}$, denoted as $T^{-1}_{2,4}$, one from the term $-C_{\cal A}{\cal A}_{\cal A}^{-1}$ denoted as $T^{-1}_{3,4}$, and $p-1$ from the term ${\cal A}_{\cal A}^{-1}$ denoted as $T^{-1}_{4:(p+2),4}$. The rows of $T^{-1}_{i,4}$ associated with the first vertex function, the edge two and three functions, and the interior functions are zero. The subset of indices that result in a zero value for $D_{4,4}T^{-1}_{i,4}$ will be referred to as the zero constraints of the fourth row for $p$, denoted as $\mu_{0}$:

[TABLE]

Then for $i\in\mu_{0}$,

[TABLE]

3.2.2 No Accurate Diagonal Scheme Exists

If there exist sufficient non-degenerate constraints, an invertible linear system defining the elements of the fourth row of $T$ and $D_{4,4}$ is defined. Adding up the number of entries in $\mu_{0}$ determines that

[TABLE]

Similarly, $\vec{T}^{*}_{4}$ as defined by (30) has $(p-1)$ degrees of freedom from the first row of ${\cal A}_{\cal A}$ and $(p-1)(p-2)/2$ degrees of freedom from ${\cal A}_{\cal I}$ for a total of

[TABLE]

degrees of freedom.

Theorem 3.1**.**

The system defined by (35) for $i\in\mu_{0}$ is both square and non-singular.

Proof 3.2**.**

Let $\vec{\phi}_{\mu_{0}}$ denote the set of basis functions from $\vec{\phi}$ that force $D_{4,4}T^{-1}_{i,4}=0$, i.e. $\phi_{i}$ for $i\in\mu_{0}$. Let $\vec{\phi}_{*}$ denote the set of basis functions of $\vec{\phi}$ associated with $\vec{T}^{*}_{4}$ as shown in (35). As discussed above, both $\vec{\phi}_{\mu_{0}}$ and $\vec{\phi}_{*}$ have $\frac{1}{2}p(p-1)$ elements. The elements of $\vec{\phi}_{\mu_{0}}$ share a common factor $\phi_{A}$ because all of the functions are zero along edge ${\cal A}$. This can be seen by referring to (15), (16), and (17). Similarly, every element of $\vec{\phi}_{*}$ shares the common factors $\phi_{B}$ and $\phi_{C}$ because the function is required to be zero along edges ${\cal B}$ and ${\cal C}$. This can also be verified by direct examination of the basis. Define $\hat{\phi}_{\mu_{0}}$ and $\hat{\phi}_{*}$ such that

[TABLE]

and

[TABLE]

The elements of $\vec{\phi}_{\mu_{0}}$ range in degree from $1$ through $p-1$. A polynomial of degree $1$ was factored from this. This leaves us with $\frac{1}{2}p(p-1)$ polynomials of degree [math] through $p-2$, or consequently $\text{dim}\left(\mathcal{T}(p-2)\right)$. By Corollary 2.1, these are linearly independent in $\mathcal{T}(p-2)$. So, the elements of $\hat{\phi}_{\mu_{0}}$ form a basis for $\mathcal{T}(p-2)$. Similarly, the elements of $\vec{\phi}_{*}$ range in degree from $2$ through $p$. Out of $\vec{\phi}_{*}$ two vertex functions were factored, leaving the elements of $\hat{\phi}_{*}$ ranging in degree [math] through $p-2$. Therefore, the elements of $\hat{\phi}_{*}$ form a basis for $\mathcal{T}(p-2)$ as well.

By the equivalence of basis functions, there exists an invertible transformation, not necessarily of similar shape to (24), but invertible none-the-less, that maps between $\hat{\phi}_{\mu_{0}}$ and $\hat{\phi}_{*}$. Call the transformation $H$. Therefore

[TABLE]

The matrix that is required to identify the fourth row of $T$ for the constraints defined by $\vec{\phi}_{\mu_{0}}$ is

[TABLE]

This is rewritten using the transformation of bases $H$:

[TABLE]

As the transformation $H$ has only constant elements, it may be factored out of the resulting integral

[TABLE]

To see that this matrix is non-singular first notice that as $H$ is a change of basis, it is non-singular. Therefore, as the product of invertible matrices is necessarily invertible, it suffices to show that the matrix defined by the integral, which we define as $\hat{M}$, is symmetric positive definite, i.e. invertible.

The symmetry is a byproduct of the symmetry of multiplication:

[TABLE]

To see that $\hat{M}$ is positive definite, consider an arbitrary vector $\vec{x}$, then

[TABLE]

We rewrite this as

[TABLE]

Therefore

[TABLE]

But, $\phi_{A}$, $\phi_{B}$, $\phi_{C}$, and the sum squared are all greater than or equal to zero on $\Omega_{ref}$, so

[TABLE]

*Furthermore, because the product $\phi_{A}\phi_{B}\phi_{C}$ is nonzero on the entirety of the interior, and the sum squared term is non-negative, $\vec{x}^{T}\hat{M}\vec{x}$ can only be zero if $\vec{x}\,=\,0$. With that, it has been shown that $\hat{M}$ is symmetric positive definite and invertible. And, as the product of two invertible matrices is invertible, $H\hat{M}$ is necessarily non-singular. *

Given the invertibility of $H\hat{M}$, we may now prove that no $(p-1)$-exact diagonal $C^{0}$ TFEM exist.

Theorem 3.3**.**

There are no continuous triangular $(p-1)$-exact diagonal pseudo-mass matrix methods for $\mathcal{T}(p)$.

Proof 3.4**.**

*To begin, assume that there is a basis $\vec{\psi}$ for $\mathcal{T}(p)$ for which there exists an associated $(p-1)$-exact diagonal pseudo-mass matrix. This basis $\vec{\psi}$ is a transformation of the modified Dubiner basis $\vec{\phi}$ under a change of basis $T$ described in (24). In particular, the first edge function of $\vec{\psi}$, $\psi_{4}$ is defined by $\vec{T}^{*}_{4}$ in (29). But Thm. 3.1 shows that $\vec{T}^{*}_{4}$ is zero, which contradicts the invertibility of $T$. *

4 A Lower-Triangular Method

In this section, we define a $(p-1)$-exact, lower-triangular pseudo-mass matrix approach for $p=3$. The intent of this section is not to develop a general method for arbitrary $p$, but rather to demonstrate that further work could lead to a $(p-1)$-exact triangular finite element method that does not require the inversion of a full mass matrix. That is, this section demonstrates hope for arbitrarily high-order explicit methods by extending the tools created in the contradiction of a $(p-1)$-exact diagonal pseudo-mass matrix to define other less restrictive approaches.

In [16], Helenbrook found a unique set of vertex functions for $\mathcal{T}(p)$ whose coefficients could be determined exactly by a diagonal operation when $f\in V_{p-1}^{h}$. Since there is no diagonal approach to determine edge function coefficients, we relax the diagonality constraint of the pseudo-mass matrix in order to define a higher-order method that is appropriate for explicit time-stepping. We rearrange the basis functions and consequent degrees of freedom to emphasize the desired lower-triangular structure. Let $L$ be the single element pseudo-mass matrix for $p=3$ defined as

[TABLE]

The degrees of freedom in (40) are organized by vertex functions $\psi_{A}$, $\psi_{B}$, $\psi_{C}$, edge functions $\psi_{{\cal A},1}$, $\psi_{{\cal B},1}$, $\psi_{{\cal C},1}$, $\psi_{{\cal A},2}$, $\psi_{{\cal B},2}$, $\psi_{{\cal C},2}$, and then the interior function $\psi_{{\cal I},1}$. The first three rows of $L$ correspond to the diagonal equations to determine the vertex coefficients where $v$ is the diagonal entry. The next three rows determine the first edge function coefficients independently for each edge, but allow a dependence on the known vertex coefficients. $e_{i}$ denotes the diagonal term for the $i^{th}$ edge function and the edge-vertex couplings are denoted $e_{i,o}$, $e_{i,-}$, and $e_{i,+}$ where the $o$, $-$, or $+$ indicates the opposite, clockwise, or counterclockwise vertex from the edge, respectively. For the second function on each edge, coupling to the known first edge coefficients are allowed. These couplings are denoted $e_{2,1-}$ and $e_{2,1+}$ in which the pair $2,1$ denotes the coupling from function 2 to function 1 and the $-$ or $+$ describe the relative edge location of the function 1 coefficient (clockwise or counterclockwise from the function 2 coefficient edge). For the one interior function at $p=3$, the subscripted $i$ entries are the couplings to the known vertex and edge coefficients as indicated by the subscript.

The pseudo-mass matrix is intended to work with a basis $\vec{\psi}$, which is expressed as a change of basis, $T$, of $\vec{\phi}$, the modified Dubiner basis. We include an additional constraint that the change of basis of the edge functions (${\cal A}_{\cal{A}},{\cal B}_{\cal{B}}$ and ${\cal C}_{\cal{C}}$ in (24)) must be upper triangular. This implies that on any edge, only the addition of higher edge functions to lower order functions is permitted. This restriction eases the analysis of the systems defined in this section. In order to work along with $L$, the organization of $T$ used in this section differs from that which was used previously in the paper. Following the ordering of $L$, the change of basis for $p=3$ is defined as

[TABLE]

4.1 Defining $T$ and $L$ for $p=3$

The matrices $L$ and $T$ are unknown. Being that the desired method is $(p-1)$-exact, the pseudo-mass matrix $L$ will define a coefficient vector $\vec{u}_{\psi}$ by

[TABLE]

that satisfies $\vec{u}_{\psi}^{T}\vec{\psi}=f$ for any $f\in\mathcal{T}(p-1)$.

Following the same arguments that led to (27), $\vec{u}_{\psi}=T^{-T}\vec{u}_{\phi}$, where $\vec{u}_{\phi}$ is again the coefficients of the representation of $f$ in the modified Dubiner basis. Thus,

[TABLE]

Letting $\mu$ denote the set of indices of the Dubiner functions of degree $\leq 2$, i.e. $(p-1)$-exact, then similar to (33), we have

[TABLE]

in which

[TABLE]

and $T^{-1}_{\mu,:}$ can be derived from (28) as

[TABLE]

To determine a vertex function we start with the first row of $L$. The only non-zero entry of the first row of $L$ is $v=L_{1,1}$, so just information in the first row of $T^{-T}$ is necessary. Multiplying by the first row of $L$ for each $i\in\mu$ gives

[TABLE]

in which $\delta(\cdot,\cdot)$ denotes the Kronecker delta function. This corresponds to 6 linear equations in 6 unknowns – $v$ and the 5 coefficients in the first row of $T$. These equations can be inverted which determines the transformation to define the vertex function $\psi_{A}$. A similar procedure can be used to find the other two vertex functions. This defines the first three rows of $L$ and $T$.

The first function on edge ${\cal A}$ is determined using row 4 of $L$. For rows 4-6 of $L$, only the first 6 column entries are non-zero. These are multiplied by the first 6 rows of $T^{-T}$ for each $i\in\mu$ (shown in (45)). The entries in the first 6 rows of this matrix are all known as they have been determined by the equations for the vertex functions. Thus, for row 4 of $L$, the left hand side only involves linear combinations of the unknowns $e_{1,o}$, $e_{1,-}$, $e_{1,+}$ and $e_{1}$. The right hand side, $T_{4,:}M_{:,i}$, involves the two unknowns ${{\cal A}}_{{\cal A},1,2}$ and ${{\cal A}}_{{\cal I},1,1}$. This again is a linear system of 6 equations in 6 unknowns. A similar process can be followed using rows 5 and 6 to determine the first edge function on each edge.

Using row 7 of $L$ a similar process is followed to determine the next edge function on edge ${\cal A}$. Now, the first 9 rows of (45) are all known and the left hand side becomes linear equations in the variables $e_{2,o}$, $e_{2,-}$, $e_{2,+}$, $e_{2,1+}$, $e_{2,1-}$, and $e_{2}$. The right hand side only involves the unknown ${{\cal A}}_{{\cal I},2,1}$. There are 6 equations in 7 unknowns, so the solution is not unique. To obtain a unique solutions, somewhat arbitrarily, $e_{2}$ was chosen to be 1. Repeating the process for rows 8 and 9 of $L$ completely determines $T$.

The last row of $L$ gives 6 equations in 10 unknowns which is the dimension of the space. One solution is to let the last row of $L$ be equal to the last row of the exact mass matrix produced by $\vec{\psi}$ which is now determined. The equations corresponding to the last row of $L$ will then be satisfied for all $f\in\mathcal{T}(p)$.

The lower triangular pseudo-mass matrix determined by this process is

[TABLE]

and the corresponding change of basis $T$ defining the basis $\vec{\psi}$, with which $L$ functions is

[TABLE]

Figure 3 shows the basis functions determined by $T$. The vertex functions are localized near their corresponding vertex similar to a nodal basis. The edge functions remain modal along the edge but become localized near the edge in the interior of the element.

4.2 A Brief Numerical Test

This section presents a rudimentary numerical test of the $(p-1)$-exact lower-triangular method given by $L$ and $T$ in (47) and (48) for $p=3$. The method is implemented as in (42) via the two matrices $L$ and $T$ with the modified Dubiner basis functions. The method was used on a series of refined meshes to project a function outside $\mathcal{T}(3)$, and the resulting $L_{2}$ errors of the projection were determined. The convergence rate of the projection is then compared to the expected rate of $3$.

To implement the lower triangular approach, (42) is written as

[TABLE]

In this case, $T$ is the global change of basis which is defined using the local $T$ matrices. The global $L$ is assembled from element $L$ matrices in the normal FEM assembly process. The term to the right of $L^{-1}$ is assembled element by element by integrating with respect to the Dubiner basis on each element, multiplying by the element $T$ matrix and then performing an FEM assembly. $L$ is inverted by first inverting all of the vertex rows, which are diagonal, then doing all of the first edge functions in the mesh, then the second functions, and finally the interior functions. Each of these operations only involves bringing the the known lower-diagonal terms of $L$ to the right hand side and then dividing by the diagonal. This then determines the coefficients $\vec{u}_{\psi}$. Multiplying by $T^{T}$, which again involves only local operations, determines $\vec{u}_{\phi}$.

The spatial convergence of the method is tested using the function

[TABLE]

over a sequence of structured meshes with reducing mesh size $h$ as seen in Fig. 4.

The approximation to $f(x,y)$ is $u^{h}(x,y)$ as defined in (49). The $L_{2}$ error is defined in the typical way as

[TABLE]

This integral is numerically computed over each element of $\Omega^{h}$ by means of GL integration for each $h$ refinement, and the results are summarized in Table 1. The step size $h$ ranges from $2^{-2}$ through $2^{-7}$. The errors are measured in the typical $L_{2}$-norm, (51). The errors and step sizes are used to approximate the spatial convergence rates by $\frac{log(error_{k}/error_{k-1})}{log(h_{k}/h_{k-1})}$. The convergence rate approaches $3$ as the meshes are refined, which is what is expected for a method that is exact for quadratic functions.

5 Conclusion

We have proven that there is no $(p-1)$-exact, diagonal pseudo-mass matrix associated with any basis $\vec{\psi}$ of $\mathcal{T}(p)$ (Thm. 3.3). However, by adapting the techniques used for this proof, we have developed a $(p-1)$-exact $C^{0}$ TFEM for $p=3$ that is appropriate for explicit time advancement as it avoids the need to perform a global inversion of the mass matrix. Similar to the spectral element method, this lower triangular inversion attains a $p^{th}$-order spatial convergence rate while avoiding the need to work with a large mass matrix. In our future work, we will explore whether this method extends to arbitrarily high-order.