1 Introduction

The numerical solution of the wave equation in time involves a large sparse mass matrix. Mass lumping avoids the cost of its lower-upper decomposition, leaving a diagonal matrix that can be readily inverted. On the reference element, its entries are proportional to numerical quadrature weights. If applied to the regular higher-order polynomial elements, there is a zero weight for degree 2, or a negative one for degree 3. For higher degrees, the resulting spatial accuracy is too low for a second-order partial differential equation. The linear degree-1 element is the exception. For higher degrees, positive weights and sufficient spatial accuracy can be obtained by choosing the degree of the polynomials in the interior higher than on the element boundary. In this way, elements of degree 2 [1,2,3], 3 [4, 5], 4 [6], 5 [7], 6 [8], 7 and 8 [9, 10], and 9 [9] were found. The results for degrees 7 and 8 in [9] and [10] are identical, but for degree 9, the one in [10] has more nodes and degree 10 on the edges.

The generalization of this approach to tetrahedra led to an element of degree 2 [6] and two elements of degree 3 [6]. The number of additional nodes in 3D is rather large, compared to the standard tetrahedral element.

Recently, it was found that a sharper, less restrictive condition can be imposed to preserve spatial accuracy after mass lumping [11]. The resulting tetrahedral elements for degree 2 and 3 have less nodes and are more efficient than the older ones. A number of degree-4 elements were found as well.

Here, the same accuracy condition is used in 2D for the construction of triangular mass-lumped elements. In addition, the moment equations for exact quadrature are based on symmetric polynomials instead of bi-variate monomials modulo symmetry to ensure that the number of independent equations is the same as the number of polynomials. With the old accuracy criterion, this could be accomplished by choosing a suitable subset [7, 12]. This enables quicker scanning through potential nodes patterns and elements that result in a number equations that does not exceed the number of unknowns. Though not strictly necessary, the last condition increases the chance of finding a solution.

The approach comprises the following steps. Given the polynomial degree of the element and a higher interior degree, one or more symmetric node patterns have to be found to support the polynomials. Next, a set of moment equations is formed by requiring numerical quadrature to be exact for elements of the Cartesian product of the basis polynomials and polynomials of the element degree minus two. This condition is less restrictive than the classic one that requires exactness for the element degree plus interior degree minus two. The resulting system is linear in the quadrature weights with coefficients that are symmetric polynomials in the node parameters. To increase the chance of finding a solution, the number of unknowns, consisting in weights and node parameters, is required to not exceed the number of independent equations. Among the solutions of the system, if any, only the real ones with strictly positive weights and node parameters that do not move nodes outside their symmetry class or the triangle are acceptable. From multiple solutions that correspond to the same set of nodes and weights modulo symmetry, only one is selected. Because the number of solutions grows quickly with the degree and system size, numerical root finding is required for the higher degrees.

When multiple acceptable solutions are found, which one is the best? The answer depends on the problem. When the wave equation is solved, the performance is affected by the number of nodes, the maximum allowable size of the time step and the resulting error. For elements of a given degree, on a rectangular domain with a decent mesh and material properties piecewise constant per element, the error behaves as element size to a power that is one higher than the degree, but the error constant varies among the multiple solutions, as does the Courant-Friedrichs-Lewy (CFL) number [13] that determines the maximum size of the time step.

To assess the performance of the various element, the numerical errors for a homogeneous problem with a point source were measured as a function of element size as well as run time. A comparison in terms of efficiency to other higher-order discretisation schemes on triangles, such as discontinuous Galerkin methods [14,15,16], the spectral difference method [17] and summation-by-parts operators [18], has not been considered. The last method avoids problems with the unisolvency condition, which will be encountered here for some choices of nodes patterns.

2 Construction Method

2.1 Node Sets and Basis Functions

A point on the triangle with vertices \({\mathbf {x}}_A\), \({\mathbf {x}}_B\) and \({\mathbf {x}}_C\) is represented by \({\mathbf {x}}= x_0 {\mathbf {x}}_A + x_1 {\mathbf {x}}_B+x_2 {\mathbf {x}}_C\) with \(x_0=1-x_1-x_2\) and \(x_k\in [0,1]\) for \(k=0,1,2\). Given a single node with natural coordinates \((x_1,x_2)\), its symmetric counterparts on the triangle and itself are elements of the set \(\Pi (x_1,x_2)\), which can be obtained by taking the first two of each of the 6 permutations of \((x_0,x_1,x_2)\) and removing duplicates. Table 1 lists the 6 equivalence classes for symmetric sets of nodes on the triangle. The second column contains the symmetry class, following definition 2.3 of [19], which will be referred to by the index j of the first column in what follows. The third column describes the node position on the triangle in natural coordinates modulo symmetry, and the fourth column their number \(v_j\) based on symmetry.

Table 1 Classes for the triangle
Full size table

A node pattern or rule set is given by \(K=\{K_1,K_2,\ldots ,K_6\}\). The total number of nodes per triangle is \(n_\mathrm {p}=\sum _{j=1}^6 K_j v_j\). The number of node parameters is \(n_\mathrm {xpar}=K_3 + K_5 + 2 K_6\). Together with the quadrature weights, of which there are \(n_\mathrm {weight}=\sum _{j=1}^6 K_j\), the total number of parameters becomes \(n_\mathrm {par}=n_\mathrm {xpar}+n_\mathrm {weight}\).

The polynomial set of basis function of degree \(p\) is \(P_p=\{ x_1^i x_2^j \, \vert \, i \ge 0, \, j\ge 0, \, i+j\le p; \, i,j\in {\mathbb {Z}} \}\). Since this choice leads to non-positive quadrature weights already for degree 2, it can be augmented by higher-degree polynomials in the interior that vanish on the edges: \({\overline{P}}_{{p'}}=\{\eta \}\otimes P_{{p'}-3}=\{ \eta \,x_1^i x_2^j \,\vert \, i\ge 0,\, j\ge 0,\, i+j\le {p'}-3;\, i,j\in {\mathbb {Z}} \}\) with a bubble function \(\eta =x_0 x_1 x_2=(1-x_1-x_2)x_1 x_2\) that vanishes on the edges. Here, \({p'}\ge 3\) and \({p'}\ge p\).

Instead of \({\overline{P}}_{{p'}}=\{\eta \}\otimes P_{{p'}-3}\), the product of the bubble function \(\eta \) and some subset \({{{\tilde{P}}}}_{{p'}-3}\) of \(P_{{p'}-3}\) can be chosen, with \(P_{p-3}\subseteq {{{\tilde{P}}}}_{{p'}-3} \subseteq P_{{p'}-3}\), similar to the tetrahedral elements of degree 4 in [11]. This option will not be considered here, because it is less obvious how to automate the search for new elements, and is left for further research.

Conformity between elements is obtained if the restriction to the edges has degree \(p\). A necessary condition for unisolvence is that there are as many nodes as polynomial basis functions \(n_\mathrm {pol}\). Together, they define the node pattern for given degrees \(p\) and \({p'}\), as described by Theorem 1 in Appendix C.

Conjecture 3.1 in [10] states another requirement for unisolvence: the node pattern in the interior should be the same as that of the regular element of degree \({p'}\) with equidistant nodes in the natural coordinates. Lemma 1 in Appendix C proves its consistency with Theorem 1.

2.2 Quadrature

Lumping of the mass matrix amounts to numerical quadrature. The usual requirement is that polynomials up to degree \(q=p+{p'}-2\) should be integrated exactly for accuracy to be preserved [5, 20, e.g.]. The elements up to degree 9 mentioned in the introduction are all based on this criterion. A less restrictive condition was introduced in [11] and enabled the construction of tetrahedral elements with less nodes. Quadrature should be exact for polynomials in \(Q_{p,{p'}}=P_{p-2}\otimes U\), where \(U=P_{p}\oplus {\overline{P}}_{{p'}}\) contains the basis functions. If \(\psi (x_1,x_2)\) is an element of \(Q_{p,{p'}}\), the condition becomes

$$\begin{aligned} \int _{x_1=0}^1 {\mathrm {d}}x_1 \int _{x_2=0}^{1-x_1} {\mathrm {d}}x_2\ \psi (x_1,x_2) = \sum _{k=1}^{n_\mathrm {weight}} w_{k} \Psi _k, \quad \Psi _k=\sum _{({{{\tilde{x}}}}_1,{{{\tilde{x}}}}_2)\in \Pi ({{{\tilde{x}}}}_{k,1},~{{{\tilde{x}}}}_{k,2})} \psi ({{{\tilde{x}}}}_1,{{{\tilde{x}}}}_2). \end{aligned}$$
(1)

Here k runs over the nodes \(({{{\tilde{x}}}}_{k,1},{{{\tilde{x}}}}_{k,2})\) modulo symmetry. Since U is a subset of \(P_{{p'}}\), this potentially leads to elements with less nodes.

The requirement of exact quadrature results in a polynomial system of equations, linear in the weights \(w_k\), with coefficients that are rational numbers for the polynomials used here. It has a subsystem that corresponds to polynomials of \({\overline{P}}_{{p'}}\) and the interior nodes.

To avoid inconsistent systems, it is natural to require the number of equations \(n_\mathrm {eq}\) not to exceed the number of parameters \(n_\mathrm {par}\), although it may still happen that the overdetermined system has one or more solutions. Among the solutions, only the real ones with strictly positive weights \(w_k\) and nodes not outside the triangle are useful. In addition, degeneracies where a node moves out of its class into another should be avoided. A class 3 node should therefore have \(0<a<\tfrac{1}{2}\), a class 5 node \(0<a<\tfrac{1}{2}\) and \(a\ne \frac{1}{3}\), and a class 6 node \(a\ne b\), \(a>0\), \(b>0\), \(1-a-b>0\).

The number of equations \(n_\mathrm {eq}\) tends to be less than the number of bivariate monomials modulo symmetry in \(Q_{p,{p'}}=P_{p-2}\otimes U\), because some of them are integrated to zero. In the earlier approach, which required exactness up to degree \(q=p+{p'}-2\), only the monomials \(x_0^{i_0} x_1^{i_1} x_2^{i_2} \) with \(i_2=i_1\ge i_0\ge 0\) and \(\sum _{j=0}^2 i_j = i_0+2 i_1\le q\) had to be considered [7, 12]. In the current setting, it is more convenient to work with symmetric polynomials, c.f. [21].

A polynomial \(\Psi _k\) in equation (1) can be expressed in the elementary symmetric polynomials \(\{1,x_0+x_1+x_2,\xi ,\eta \}\), with \(\xi (x_1,x_2)=x_0 x_1 + x_0 x_2 + x_1 x_2\) and bubble function \(\eta (x_1,x_2)=x_0 x_1 x_2\). Because \(x_0+x_1+x_2=1\), the second polynomial drops out. The set of symmetric polynomials up to degree \(p\) is \(P^\mathrm {sym}_{p}= \{\xi ^i \eta ^j \,\vert \, i\ge 0,\, j\ge 0,\, 2i+3j\le p;\,\) \(i,j\in {\mathbb {Z}} \}\). For quadrature, the set \(Q_{p,{p'}}\) can be replaced by

$$\begin{aligned} {{{\overline{Q}}}}_{p,{p'}}=P^\mathrm {sym}_{2p-2}\cup \left[ \{\eta \}\otimes \left( P^\mathrm {sym}_{p+{p'}-5}\setminus P^\mathrm {sym}_{2p-5} \right) \right] . \end{aligned}$$
(2)

Equation (1) for an element \({{\overline{\psi }}}\) of \({{{\overline{Q}}}}_{p,{p'}}\) can be replaced by

$$\begin{aligned} \int _{x_1=0}^1 {\mathrm {d}}x_1 \int _{x_2=0}^{1-x_1} {\mathrm {d}}x_2\ {{\overline{\psi }}}(x_1,x_2) = \sum _{k=1}^{n_\mathrm {weight}} w_{k} v_{j(k)}{{\overline{\psi }}}({{{\tilde{x}}}}_{k,1},{{{\tilde{x}}}}_{k,2}), \end{aligned}$$
(3)

where j(k) is the class of node number k modulo symmetry, i.e., for each symmetric set of nodes, only one is included. The number of independent polynomials in \(P^\mathrm {sym}_{p}\) is \(n^\mathrm {sym}_{p}=1+\mathrm {floor}\left( {p}({p}+6)/12\right) \) [21]. Their number in \({{{\overline{Q}}}}_{p,{p'}}\) is the same as the resulting number of equations: \(n_\mathrm {eq}=n^\mathrm {sym}_{2p-2}+n^\mathrm {sym}_{p+{p'}-5}-n^\mathrm {sym}_{2p-5}\), assuming \({p'}>p\ge 2\).

The use of symmetric polynomials suggests that (a, 0) be better replaced by \(a=\tfrac{1}{2}(1-\sqrt{1-\alpha })\) and \(0\le \alpha \le 1\), as in [8]. A similar substitution for a class 6 node (a, b) is \({\alpha }=3\xi (a,b) =3[a(1-a)+b(1-b)-a b]\), \({\beta }=27\eta (a,b)=27 a b (1-a-b) \). Appendix D describes this parametrization and its inverse, which proved to be useful for the element of degree 5. It lowers the degrees in the system of polynomial equations but may be less attractive for numerical root finding with bounds.

2.3 Polynomial System

Given a polynomial degree \(p\), the values of \(K_{1,2,3}\) of the node pattern are prescribed by Theorem 1 in Appendix C. Because the vertices are always included, \(K_1=1\). Then, a loop over the interior degree \({p'}\) is started, defining the remaining \(K_{4,5,6}\) by Theorem 1, to find combinations with a number of free parameters \(n_\mathrm {par}\) not exceeding the number of independent equations \(n_\mathrm {eq}\). The search is limited to values of the total number of nodes \(n_\mathrm {p}\) not much larger than that for the known elements, constructed with the more restrictive accuracy criterion.

For each candidate node pattern and degrees \(p\), \({p'}\), equation (3) provides a system of equations of the form

$$\begin{aligned} {\mathbf {f}}({\mathbf {a}},{\mathbf {w}})={\mathbf {A}}({\mathbf {a}})\begin{pmatrix}{\mathbf {w}}\\ 1\end{pmatrix}={\mathbf {0}}, \end{aligned}$$
(4)

where the vector \({\mathbf {w}}\) contains the \(n_\mathrm {weight}\) integration weights and the vector \({\mathbf {a}}\) the \(n_\mathrm {xpar}\) node parameters. The matrix \({\mathbf {A}}({\mathbf {a}})\) with \(n_\mathrm {eq}\) rows and \(n_\mathrm {weight}+1\) columns contains polynomials in the node parameters up to degree \(p+{p'}-2\) and rational coefficients. A subsystem corresponding to the product of the bubble function and polynomials only contains the interior weights and node parameters.

The system may have no solution, which typically but not necessarily happens for \(n_\mathrm {eq}>n_\mathrm {xpar}\). It may be underdetermined and then either has no solution, if inconsistent, or infinitely many. This happens typically but not necessarily if \(n_\mathrm {eq}<n_\mathrm {xpar}\). If the number of complex solutions is finite, the system is called zero-dimensional. In that case, the number of solutions does not exceed the Bézout bound [22,  e.g.], which grows rapidly with the degree and size of the system. Among those solutions, only the real ones with positive weights and node parameters within the given bounds are useful, if they exist at all.

There are several packages for solving polynomial systems of equations but their unfavourable computational complexity limits their usefulness to small systems. For the results presented later on, the function Gröbnerbasis in Mathematica [23], an implementation of the Buchberger’s algorithm [24], was nevertheless used with a time limit because it returns relatively quickly if the system is inconsistent. For small systems, the function Reduce provided solution sets, from which the real-valued ones within the parameter bounds were selected, if present, and only one for each set of symmetric solutions was retained. For larger systems with \(p\ge 5\), an empirical approach was adopted with a hand-coded damped Newton-type algorithm using extended numerical precision in Mathematica, starting from randomly chosen initial values. Better results were obtained with starting values obtained with lsqnonlin in Matlab [25], a nonlinear least-squares algorithm with bounds using a trust-region method [26], and refining them with the Newton-type algorithm in Mathematica. The authors of [27, 28] take a similar approach with additional refinements. In some cases, this provided convergence with an extended precision of 64 or even 128 digits. In other cases, the Matlab solution proved to be false, due to its insufficient precision.

Table 2 Node patterns for the less restrictive accuracy criterion
Full size table

3 New Elements

Table 2 lists candidate node patterns for a small number of nodes \(n_\mathrm {p}\) and a number of parameters \(n_\mathrm {par}\) not too much larger than the number of equations \(n_\mathrm {eq}\).

For degrees 2 to 4, the elements obtained with the less restrictive accuracy criterion are the same as the known ones. New elements are found for degree 5 to 7.

The degree-5 element has the same node pattern and number of nodes as the known element but, with the less restrictive accuracy criterion, the number of equations is now one less than the number of parameters and infinitely many solutions exist. Some of the elements listed in Table 3 have a much better CFL number than the old element.

The value the CFL number in the table is a crude estimate computed on a single reference element with natural boundary conditions, as in [8]. Note that for the performance comparison later on, the power method is used.

Figure 1 depicts the node distributions. Versions ‘A’ to ‘E’ were found by numerical root finding. As it turns out, the node parameters and weights could actually be obtained as functions of a single parameter related to the class-6 node. The other parameter for that node is given by a root of a quartic equation. Two of the four roots provide acceptable solutions, each for a certain range of the one parameter. From the two parameters, the other unknowns follow. Appendix E describes the details. With this approach, the elements ‘F’ and ‘G’ in Table 3 were found.

Table 3 Elements of degree \(p=5\) and \({p'}=7\) for node pattern \(K=\{1,0,2,0,3,1\}\) and 30 nodes
Full size table
Fig. 1
figure 1

Elements of degree 5 with 30 nodes

Full size image
Table 4 As Table 3 but for elements of degree \(p=6\) and \({p'}=8\), node pattern \(K=\{1,1,2,0,3,2\}\) and 39 nodes.
Full size table

For degree 6, Table 4 lists two new elements for \({p'}=8\) and \(K=\{1,1,2,0,3,2\}\), one with a much larger estimated CFL than the other. Figure 2 displays their node patterns. They have 39 nodes, significantly less than the 46 for the elements found in [8]. Note that the related polynomial system of equations may have more acceptable solutions than the two found so far.

Table 5 contains a new degree-7 element with less nodes than the known element, but with a negative weight, making it unsuited for explicit time stepping with the wave equation. Figure 3 shows the distribution of nodes. Table 6 lists a new element for degree 7 with version name ‘old,2’, obtained with the old accuracy criterion. It has an unfavourable estimated CFL number compared to the 0.0288 for the degree-7 element in [9, 10], not included in the table.

For the same node pattern, the new accuracy criterion provided the elements ‘new,A’, ‘new,B’ and ‘new,C’, among others. Note that the elements obtained with the old accuracy criterion also obey the new criterion if the same node pattern is used, but with the latter the number of equations is less than the number of parameters. There are infinitely many solutions, similar to the elements of degree 5. Figure 4 displays the node distributions for various degree-7 elements.

The node pattern for these elements results in a polynomial system that has one equation less than the number of unknowns. Once an element is found with a Newton-type method, the Jacobian of the polynomial system has a nullspace characterized by a single vector. A small perturbation in its direction followed by a few Newton iterations will provide a new element. With this continuation approach, one can, for instance, minimize or maximize a parameter. Note that a larger value of the smallest weight, typically the first one that corresponds to the vertices, as suggested by [10], does not necessarily imply a larger CFL number, as can be seen in Table 3 when comparing versions ‘A’ and ‘B’, for instance.

Fig. 2
figure 2

Two new elements of degree 6 with 39 nodes.

Full size image
Table 5 As Table 3 but for elements of degree \(p=7\) and \({p'}=9\), node pattern \(K=\{1,0,3,1,3,3\}\), 49 nodes, but one negative weight at the centroid
Full size table
Table 6 As Table 3 but for elements of degree \(p=7\) and \({p'}=10\), node pattern \(K=\{1,0,3,0,4,4\}\) and 57 nodes.
Full size table

4 Numerical Tests

Accuracy and performance tests were carried out for the wave equation

$$\begin{aligned} \frac{1}{\rho c^2} \frac{\partial ^2 p}{\partial t^2} + \nabla \cdot \left( \frac{1}{\rho }\nabla p\right) = f, \end{aligned}$$
(5)

with solution \(p(t,{\mathbf {x}})\), wave speed \(c({\mathbf {x}})\), density \(\rho ({\mathbf {x}})\) and forcing function \(f(t,{\mathbf {x}})=w(t)\delta ({\mathbf {x}}-{\mathbf {x}}_s)\), representing a point source located at \({\mathbf {x}}_s\) with signature w(t). The same homogeneous test problem as in [8] is used, on the unit square for constant \(c({\mathbf {x}})=1\) and \(\rho ({\mathbf {x}})=1\) and with zero Dirichlet boundary conditions (not Neumann as erroneously mentioned in [8]). The source is located at the centre and has a signature \(w(t)=\left[ 4(t/T)\{1-(t/T)\}\right] ^{16}\) for \(0\le t\le T\) and zero otherwise, with a duration \(T=0.2\).

The spatial discretisation is reviewed in Appendix A. Geompack [30] was used to construct unstructured meshes. Assembly is performed on the fly in each time step, as in [6, 11, 31,32,33], although on modern hardware, a pre-assembled matrix might be more efficient in 2D, in particular for lower orders.

Fig. 3
figure 3

New element of degree 7 with 49 nodes but a negative weight

Full size image
Fig. 4
figure 4

Three new and two old elements of degree 7 with 57 nodes. Version ‘7,1 old’ (c) was found earlier [9, 10]. The newly found ‘7,2 old’ has an unfavourable CFL number

Full size image

Higher-order time stepping [34,35,36,37] is applied instead of Stork’s dispersion correction [38,39,40,41], which would be more efficient in practice. The time-stepping order is taken to be \(M_t=2\,\mathrm {ceil}((p+1)/2)\). Details, such as the increase of the CFL number for higher orders, can be found in [8,  e.g.]. If the signature is not smooth enough at the endpoints, at time 0 or T, the numerical error for high-order time stepping may become too large because a higher derivative of the source signature ceases to exist. A higher power will avoid that problem and is chosen here as 16 instead of the power 8 used [8].

Fig. 5
figure 5

Convergence under mesh refinement for degree \(p=5\). The legend contains the degree and the order for time stepping, the number of nodes, 30 in this case, and the version name. The left panel a shows the relative L\(_2\) error as a function of the square root of the number of degrees of freedom N. The grey line indicates the expected trend. The panel on the right b depicts errors as a function of the elapsed time for time stepping. Similar results for the maximum error, relative to the maximum of the solution, are shown in c and d

Full size image
Fig. 6
figure 6

As Fig. 5, but for degree 6

Full size image
Fig. 7
figure 7

As Fig. 5, but for degree 7. Versions ‘1’ from [9, 10] and ‘2’ were obtained with the old, more restrictive accuracy criterion, whereas ‘A’ to ’C’ were found with the new criterion

Full size image

For second-order times stepping, the time step is limited by \({\Delta t}_{\max }=\sqrt{2/\sigma _{\max }}\), where \(\sigma _{\max }\) is the largest eigenvalue of the spatial operator, the product of the inverse lumped mass matrix and the stiffness operator. The power method [42] was used to find an estimate for each run. The results turned out be slightly larger than the crude estimates listed in the tables.

The numerical solution at time \(t=1.25\) is compared to the exact solution on all the nodes. Their number equals the number of degrees of freedom N, which includes the zero Dirichlet boundary values. The relative L\(_2\) error is defined as \(\Vert u-u_{\mathrm {exact}}\Vert _2 /\Vert u_{\mathrm {exact}}\Vert _2\), with a norm given by

$$\begin{aligned} \Vert u \Vert _2^2 = \sum _{j=1}^{n_\mathrm {tri}} a_j^{} \sum _{k=1}^{n_\mathrm {p}} w_k^{} u_{j,k}^2, \end{aligned}$$
(6)

where j runs over all triangle and k over all nodes in each triangle with area \(a_j\), quadrature weights \(w_k\) and solution values \(u_{j,k}\). With a number of degrees of freedom \(N\), the expected convergence behaviour is \(N^{-(p+1)/2}\) for the current problem. Note that the use of these weights in equation (6) results in an additional error, which should be neglegible for \({p+{p'}-2}>p\) or \({p'}>2\).

The computational efficiency of an element is not only determined by the number of nodes, but also by its CFL number and the error constant. A crude indication of the efficiency is obtained by performing tests on the given problem and recording its run time. The single-threaded code was written in C and dates back to 1995. The wall clock time for the time stepping part of the code is by no means representative for what can be obtained on modern architectures, but is nevertheless used as a measure for the performance comparison.

Figure 5a shows the measured relative L\(_2\) error on a sequence of unstructured grids for three versions of the scheme, including the old one, as a function of the square root of number of degrees of freedom. The trend follows the expected behaviour, \(N^{-(p+1)/2}\), indicated by the grey line. Of the three version in this comparison, ‘B’ has the smallest error. Figure 5b plots the error against the elapsed time and provides a rough indication of the relative performance. Versions ‘B’ and ‘G’ appear to be slightly more efficient than the other, and both do better than the old scheme. Although version ‘G’ has a smaller CFL number than ‘F’, its better accuracy make its more efficient. Figure 5(c) and (d) show similar results measured for the maximum error, relative to the maximum absolute value of the solution. The maximum error tends to be less smooth than the L\(_2\) norm [43,  e.g.]. The performance for the new elements is similar, whereas the gain over the old element is small.

Figure 6a confirms the expected convergence behaviour for degree 6. Both new elements ‘A’ and ‘B’ have a smaller error than the known scheme with 46 nodes, version ‘D’, in [8]. In terms of efficiency, Fig. 6b shows that version ‘A’ outperforms the others. Version ‘B’ is less efficient than the old one. Although it has less nodes per element, its smaller CFL number make it less attractive. The results in the relative maximum norm, displayed in Fig. 6(c) and (d), show a similar behaviour.

Figure 7 display errors for two degree-7 elements, version ‘1’ of [9, 10] and ‘2’ obtained with the old, more restrictive accuracy criterion. Elements ‘A’,‘B’ and ‘C’ were found with the new criterion. Their errors are similar, but ‘C’ is more efficient.

5 Conclusions

The search for continuous mass-lumped triangular finite elements with a less restrictive accuracy criterion than the classic one has led to a number of new elements of degree 5, 6 and 7. These elements are enriched with polynomials of a higher degree in their interior to preserve accuracy after mass lumping. Given the polynomial degrees on the boundary of the triangle and in its interior, it is proven that the requirements of conformity and unisolvence define a unique node pattern.

The elements for lower degrees turned out to be same as the known ones. For degree 5 and 7, there are infinitely many elements with the same number of nodes as the known element. For degree 6, two new elements with less nodes were found.

Numerical tests on a homogeneous wave-equation problem with a point source at the centre show that some of the degree-5 elements have a somewhat better accuracy at a given cost than the classic element. The same is true for one of the new elements of degree 7. The gain is most pronounced for one of the two degree-6 elements.