Trefftz Finite Elements on Curvilinear Polygons

TL;DR

This paper introduces a Trefftz finite element method for meshes with curvilinear polygons, utilizing integral equations for basis functions, and demonstrates optimal convergence even with singularities and curved edges.

Contribution

It develops a novel Trefftz finite element approach on curved polygonal meshes, including new definitions of polynomial edge functions and an interpolation operator with proven optimal convergence.

Findings

Finite element solutions achieve optimal order convergence.

Method works effectively on meshes with curved edges.

Exploits singular functions for improved convergence without adaptive refinement.

Abstract

We present a Trefftz-type finite element method on meshes consisting of curvilinear polygons. Local basis functions are computed using integral equation techniques that allow for the efficient and accurate evaluation of quantities needed in the formation of local stiffness matrices. To define our local finite element spaces in the presence of curved edges, we must also properly define what it means for a function defined on a curved edge to be "polynomial" of a given degree on that edge. We consider two natural choices, before settling on the one that yields the inclusion of complete polynomial spaces in our local finite element spaces, and discuss how to work with these edge polynomial spaces in practice. An interpolation operator is introduced for the resulting finite elements, and we prove that it provides optimal order convergence for interpolation error under reasonable assumptions. We provide a description of the integral equation approach used for the examples in this paper, which was recently developed precisely with these applications in mind. A few numerical examples illustrate this optimal order convergence of the finite element solution on some families of meshes in which every element has at least one curved edge. We also demonstrate that it is possible to exploit the approximation power of locally singular functions that may exist in our finite element spaces in order to achieve optimal order convergence without the typical adaptive refinement toward singular points.

Tables

Table 1. Table 1. Discretization errors, | u − u ^ | 1 subscript 𝑢 ^ 𝑢 1 |u-\hat{u}|_{1} , and error ratios for the Shuriken, Pegboard, and Jigsaw meshes. For the Shuriken meshes, both Type 1 and Type 2 elements are used. For the other meshes, only Type 2 elements are used.

	Shuriken, Type 1		Shuriken, Type 2		Pegboard		Jigsaw
$r$	${\left|u-\hat{u}\right|}_{H^1\left(\Omega \right)}$	ratio	${\left|u-\hat{u}\right|}_{H^1\left(\Omega \right)}$	ratio	${\left|u-\hat{u}\right|}_{H^1\left(\Omega \right)}$	ratio	${\left|u-\hat{u}\right|}_{H^1\left(\Omega \right)}$	ratio
4	7.882e-02		5.629e-02		3.470e-02		3.209e-02
8	8.200e-02	0.961	2.847e-02	1.977	1.726e-02	2.010	1.538e-02	2.086
16	8.813e-02	0.930	1.429e-02	1.993	8.537e-03	2.022	7.559e-03	2.035
32	9.232e-02	0.955	7.150e-03	1.998	4.233e-03	2.017	3.754e-03	2.014
64	9.470e-02	0.975	3.576e-03	1.999	2.106e-03	2.010	1.871e-03	2.006
128	9.600e-02	0.987	1.788e-03	2.000	1.050e-03	2.006	9.327e-04	2.006
256	9.662e-02	0.993	8.937e-04	2.001	5.239e-04	2.004	4.628e-04	2.015
512	9.694e-02	0.997	4.463e-04	2.003	2.611e-04	2.007	2.247e-04	2.060

Table 2. Table 2. Discretization errors, | u − u ^ | H 1 ​ ( Ω ) subscript 𝑢 ^ 𝑢 superscript 𝐻 1 Ω |u-\hat{u}|_{H^{1}(\Omega)} , and error ratios for the sequence of meshes having one L-shaped cell of fixed size and increasingly fine square cells. Also included are the spectral condition numbers of the global (sparse) stiffness matrix A 𝐴 A and the small (dense) stiffness matrix A L subscript 𝐴 𝐿 A_{L} associated with the basis functions on the L-shaped cell.

$r$	${\left|u-\hat{u}\right|}_{H^1\left(\Omega \right)}$	ratio	${\kappa }_2\left(A\right)$	${\kappa }_2\left(A_L\right)$
1	1.3629e-01		1.2908e+01	5.7834e+00
2	6.7610e-02	2.0158	3.7679e+01	7.5699e+00
4	3.3734e-02	2.0042	1.2527e+02	8.9390e+00
8	1.6855e-02	2.0014	4.6339e+02	9.8197e+00
16	8.4305e-03	1.9993	1.7919e+03	1.0377e+01
32	4.2289e-03	1.9935	7.0585e+03	1.0728e+01
64	2.1468e-03	1.9699	2.8029e+04	1.2450e+01

Equations

$${\int }_{\Omega }A\nabla u\cdot \nabla v+(\mathbf{b}\cdot \nabla u+cu)vdx={\int }_{\Omega }fvdx+{\int }_{\partial {\Omega }_N}gvds\text{mbox}forallv\in \mathcal{H},$$

$$V_p(\mathcal{T})=\{v\in C(\overline{\Omega }):v_{{|}_K}\in V_p(K)\text{mbox}forallK\in \mathcal{T}\},$$

$$v\in V_p(K)\text{mbox}ifandonlyif\Delta v\in {\mathbb{P}}_{p-2}(K)\text{mbox}inK\text{mbox}andv{|}_{\partial K}\in {\mathbb{P}}_p(\partial K).$$

$$v\in V_p^K(K)\text{mbox}ifandonlyif\Delta v\in {\mathbb{P}}_{p-2}(K)\text{mbox}inK\text{mbox}andv=0\text{mbox}on\partial K,$$

$$v\in V_p^{\partial K}(K)\text{mbox}ifandonlyif\Delta v=0\text{mbox}inK\text{mbox}andv{|}_{\partial K}\in {\mathbb{P}}_p(\partial K).$$

$$\dim V_p(K)=\dim V_p^K(K)+\dim V_p^{\partial K}(K)=\dim {\mathbb{P}}_{p-2}(K)+\dim {\mathbb{P}}_p(\partial K)=\left(\genfrac{}{}{0ex}{}{p}{2}\right)+\dim {\mathbb{P}}_p(\partial K).$$

$$\dim V_p^{\partial K}(K)=\dim {\mathbb{P}}_p(\partial K)=\sum _{e\subset \partial K}\dim {\mathbb{P}}_p(e)-\#\text{mbox}edges.$$

$$\dim V_1^{\partial K}(K)$$

$$\dim V_2^{\partial K}(K)$$

$$\Delta {\psi }_{\alpha }^K=0\text{mbox}inK,{\psi }_{\alpha }^K=-q_{\alpha }\text{mbox}on\partial K.$$

$${\ell }_j(x)=1-\frac{(x-z_j)\cdot R(z_{j-1}-z_{j+1})}{(\sqrt{3}/2)h^2},$$

$$R=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right),h=|z_1-z_0|,z_2=\frac{z_1+z_0}{2}-R\cdot \frac{\sqrt{3}(z_1-z_0)}{2}.$$

$$\left\{{\ell }_0,{\ell }_1,{\ell }_2,4{\ell }_1{\ell }_2,4{\ell }_0{\ell }_2,4{\ell }_0{\ell }_1,\frac{3\sqrt{3}}{2}{\ell }_1{\ell }_2({\ell }_1-{\ell }_2),\frac{3\sqrt{3}}{2}{\ell }_0{\ell }_2({\ell }_0-{\ell }_2),\frac{3\sqrt{3}}{2}{\ell }_0{\ell }_1({\ell }_0-{\ell }_1),27{\ell }_0{\ell }_1{\ell }_2\right\},$$

$$\left\{{\ell }_2,4{\ell }_1{\ell }_2,4{\ell }_0{\ell }_2,4{\ell }_0{\ell }_1,\frac{3\sqrt{3}}{2}{\ell }_1{\ell }_2({\ell }_1-{\ell }_2),\frac{3\sqrt{3}}{2}{\ell }_0{\ell }_2({\ell }_0-{\ell }_2),\frac{3\sqrt{3}}{2}{\ell }_0{\ell }_1({\ell }_0-{\ell }_1),27{\ell }_0{\ell }_1{\ell }_2\right\},$$

$$\left\{{\ell }_0,{\ell }_1,4{\ell }_0{\ell }_2,4{\ell }_0{\ell }_1,\frac{3\sqrt{3}}{2}{\ell }_1{\ell }_2({\ell }_1-{\ell }_2),\frac{3\sqrt{3}}{2}{\ell }_0{\ell }_1({\ell }_0-{\ell }_1),27{\ell }_0{\ell }_1{\ell }_2\right\}.$$

$$\left\{{\ell }_0,{\ell }_1,4{\ell }_1{\ell }_2,4{\ell }_0{\ell }_2,4{\ell }_0{\ell }_1,\frac{3\sqrt{3}}{2}{\ell }_1{\ell }_2({\ell }_1-{\ell }_2),\frac{3\sqrt{3}}{2}{\ell }_0{\ell }_2({\ell }_0-{\ell }_2)\right\}.$$

$$\left\{{\ell }_0,{\ell }_1,{\ell }_2,\frac{3\sqrt{3}}{2}{\ell }_1{\ell }_2({\ell }_1-{\ell }_2),\frac{3\sqrt{3}}{2}{\ell }_0{\ell }_2({\ell }_0-{\ell }_2),\frac{3\sqrt{3}}{2}{\ell }_0{\ell }_1({\ell }_0-{\ell }_1),27{\ell }_0{\ell }_1{\ell }_2\right\},$$

$${\tilde{\ell }}_0(x)=\frac{(z_1-x)\cdot (z_1-z_0)}{h^2},{\tilde{\ell }}_1(x)=\frac{(x-z_0)\cdot (z_1-z_0)}{h^2},{\tilde{\ell }}_2(x)=\frac{(z_1-x)\cdot R(z_1-z_0)}{h^2}.$$

$${\ell }_0={\tilde{\ell }}_0-\frac{{\tilde{\ell }}_2}{\sqrt{3}},{\ell }_1={\tilde{\ell }}_1-\frac{{\tilde{\ell }}_2}{\sqrt{3}},{\ell }_2=\frac{{\tilde{\ell }}_2}{\sqrt{3}/2},$$

$${\int }_{\Omega }A\nabla \hat{u}\cdot \nabla v+(\mathbf{b}\cdot \nabla u+c\hat{u})vdx={\int }_{\Omega }fvdx+{\int }_{\partial {\Omega }_N}gvds\text{mbox}forallv\in V_p(\mathcal{T})\cap \mathcal{H}.$$

$${\mathcal{I}}_K:W(K)=\{v\in C(\overline{K})\cap H^1(K):\Delta v\in L^2(K)\}\to V_p(K),$$

$$\mathcal{I}:W=\{v\in C(\overline{\Omega })\cap H^1(\Omega ):\Delta v\in L^2(\Omega )\}\to V_p(\mathcal{T})$$

$$\{\begin{array}{ll}\Delta v^K=\Delta v& \text{mbox}inK\\ v^K=0& \text{mbox}on\partial K\end{array},\{\begin{array}{ll}\Delta v^{\partial K}=0& \text{mbox}inK\\ v^{\partial K}=v& \text{mbox}on\partial K\end{array}.$$

$$\{\begin{array}{ll}\Delta ({\mathcal{I}}_K^{K}v)=q^K& \text{mbox}inK\\ {\mathcal{I}}_K^{K}v=0& \text{mbox}on\partial K\end{array},\{\begin{array}{ll}\Delta ({\mathcal{I}}_K^{\partial K}v)=0& \text{mbox}inK\\ {\mathcal{I}}_K^{\partial K}v=q^{\partial K}& \text{mbox}on\partial K\end{array}.$$

$${\int }_K(\Delta v-q^K)\varphi dx=0\text{mbox}forall\varphi \in {\mathbb{P}}_{p-2}(K).$$

$$(\varphi ,\psi {)}_{H^{1/2}(\Gamma )}={\int }_{\Gamma }\varphi \psi ds+{\int }_{\Gamma }{\int }_{\Gamma }\frac{(\varphi (x)-\varphi (y))(\psi (x)-\psi (y))}{|x-y{|}^2}ds(x)ds(y),$$

$$q_e(z_0)=v(z_0),q_e(z_1)=v(z_1),(v-q_e,w{)}_{H^{1/2}(e)}=0\text{mbox}forallw\in {\mathbb{P}}_{p,0}(e).$$

$$|v-{\mathcal{I}}_{K}v{|}_{H^1(K)}^2=|v^K-{\mathcal{I}}_K^{K}v{|}_{H^1(K)}^2+|v^{\partial K}-{\mathcal{I}}_K^{\partial K}v{|}_{H^1(K)}^2,$$

$$|v^K-{\mathcal{I}}_K^{K}v{|}_{H^1(K)}$$

$$|v^K-{\mathcal{I}}_K^{K}v{|}_{H^1(K)}^2$$

$$|v^K-{\mathcal{I}}_K^{K}v{|}_{H^1(K)}$$

Full text

This work was partially supported by the NSF Grants DMS-1522471 and DMS-1624776. Numerical studies were facilitated by the Portland Institute for Computational Sciences.

Key words. finite element methods, curvilinear polygons, interpolation error analysis, Trefftz methods

AMS subject classifications. 65N30, 65R10, 65N12, 65N15

1. Introduction

Polygonal and polyhedral meshes in finite element analysis for the numerical treatment of boundary value problems have attracted a lot of interest during the last few years due to their enormous flexibility. They resolve the paradigm of a small class of element shapes (e.g. triangles, quadrilaterals, tetrahedra, etc.) in finite element methods (FEM) and therefore open the possibility for very problem adapted mesh handling. This comes with an easy realization of local mesh refinement, coarsening and adaptation near singularities and interfaces. In particular, the notion of such general meshes naturally deal with “hanging nodes”—allowing two edges of a polygon to meet at a straight angle removes the notion of hanging nodes altogether. Virtual Element Methods (VEM) (cf. [11, 2, 25, 12, 13, 5, 40, 9, 14, 4, 6, 10]), which have drawn inspiration from mimetic finite difference schemes, constitute one active line of research in this direction. Another involves Boundary Element-Based Finite Element Methods (BEM-FEM) (cf. [27, 51, 50, 81, 68, 69, 82, 52, 83, 86, 84, 85]), which have looked more toward the older Trefftz methods for motivation. A similar strategy has been followed in our previous work [3], where a Nyström approximation is applied for the treatment of local boundary integral equations instead of a boundary element method. The gained insights and flexibilities in that work build the basis of the development in this paper. A third line of research involves generalized barycentric coordinates (cf. [43, 38, 67, 42, 59, 72] and the references in [39]), that mimic certain key properties of standard barycentric coordinates over general element shapes. The before mentioned approaches yield globally-conforming discretizations, which is challenging on general meshes. However, there has also been significant interest in various non-conforming methods for polyhedral meshes. We mention Compatible Discrete Operator (CDO), Hybrid High-Order (HHO) schemes (cf. [19, 20, 18, 34, 32, 33, 21]) and Weak Galerkin (WG) schemes (cf. [78, 79, 77, 64, 63, 62, 80]), as well as the recent adaptations of the discontinuous Petrov-Galerkin method (cf. [75]), in this regard.

The present work builds upon [3] for second-order, linear, elliptic boundary value problems posed on possibly curved domains $\Omega\subset\mathbb{R}^{2}$: Find $u\in\mathcal{H}$ such that

[TABLE]

where $\mathcal{H}$ is some appropriate subspace of $H^{1}(\Omega)$ incorporating homogeneous Dirichlet boundary conditions, and standard assumptions on the data $A,\mathbf{b},c,f,g$ ensure that the problem is well-posed. Although polygonal meshes are quite flexible and have been studied intensively in recent years, there are relatively few results in this direction that allow for curved elements in the spirit of polygonal meshes, despite their natural appeal in fitting curved domain boundaries and interfaces. Early efforts at treating curved boundaries in the finite element context, such as isoparametric elements (cf. [71, 56, 15]), involve (local) mappings of standard mesh cells to fit curved boundaries, and these methods remain popular today. More recently, isogeometric analysis (cf. [28]), which integrates the use of splines both for modeling complex (curved) geometries and in constructing finite elements on the resulting meshes. This remains an active area of research. Two recent contributions employing non-conforming methods over curved polygonal elements are described in [24, 21]. In terms of conforming methods for treating curved boundaries that are in the same vein as the conforming polygonal methods mentioned in the first paragraph, we mention four, all of which are very recent. In [16], the curved boundary of the domain is approximated by polygonal elements with straight edges and a stabilization is constructed such that optimal rates of convergence are retained for high order methods. In contrast, [10] gives a first study of VEM with polygonal elements having curved edges in 2D for the treatment of curved boundaries and interfaces, but the construction results in $\mathbb{P}_{p}(K)\not\subset V_{p}(K)$, i.e. the polynomials of degree smaller or equal $p$ are locally not contained in the local approximation space of order $p$. This introduces additional difficulties in the study of approximation properties. More recently, these authors work with a richer finite element space of functions for which $\mathbb{P}_{p}(K)\subset V_{p}(K)$, see [8]. This richer space is referred to as “Type 2 elements” in our previous contribution [3], which considered the natural incorporation of Dirichlet data on curved (or straight) portions. Our present work provides a practical realization, as well as supporting interpolation theory, for Type 2 elements on very general planar meshes consisting of curvilinear polygons. As both our work and [8] must address many of the same theoretical and practical concerns, it is unsurprising that there are strong similarities between the approaches, and comparisons between them will be of interest as both are developed further.

The paper is organized as follows: In Section 2 we describe local and global finite element spaces allowing for mesh cells that are fairly general curvilinear polygons. As is done in VEM and BEM-FEM, as well as our previous Trefftz-Nyström contribution [3], the local spaces are defined in terms of Poisson problems with polynomial data on the mesh cells. For curved edges, we discuss two natural choices (those suggested in [3]) for what it means to have polynomial boundary data on curved edges. The one that we believe is the more appropriate of the two requires further explanation concerning how these edge polynomial spaces and their bases can be constructed in practice, and the bulk of Section 2 is devoted to doing so. Having defined the local and global spaces, Section 3 provides an interpolation operator, and establishes that interpolation in these spaces is at least as good as interpolation by polynomials in more standard (e.g. triangular, quadrilateral) meshes, as well as interpolation in straight-edged polygonal meshes. In brief, it is established that the inclusion of $\mathbb{P}_{p}(K)$ in our local spaces $V_{p}(K)$ provides the expected approximation power, and the presence of other (possibly singular) functions in $V_{p}(K)$ is not detrimental. In [3], an example illustrated that such locally singular functions can even be beneficial for approximation, and we develop that argument further in the final example of Section 5. Section 4 provides a description of the integral equation approach we use for computing the information about our basis functions that is needed in forming finite element stiffness matrices. This approach, which was developed with our present application in mind, is discussed in detail in our previous work [65], so we here provide a broader description of the approach and some of its key features. Finally, Section 5 provides several examples illustrating the convergence of the finite element solution on different families of meshes whose elements each have at least one curved edge, including a comparison of convergence and conditioning for Type 1 and Type 2 elements on families of meshes whose curved edges are very close to being straight. As mentioned above, in the final example of this section we both argue and demonstrate that it is possible to exploit the approximation power of locally singular functions that may exist in our finite element spaces in order to achieve optimal order convergence without the typical adaptive refinement toward singular points.

2. Local and Global Spaces

Following [26, 36], let $K$ be a connected subset of $\mathbb{R}^{2}$, with non-empty interior and compact closure, whose Lipschitz boundary, $\partial K$, is a simple closed contour consisting of a finite union of smooth arcs, see Figure 1. We will refer to $K$ as a mesh cell, the arcs as edges, and will implicitly assume that adjacent edges meet at an (interior) angle strictly between [math] and $2\pi$, i.e. $K$ has no slits or cusps. We allow adjacent edges to meet at a straight angle. The vertices of $K$ are those points where two adjacent edges meet. Given an integer $p$ and a mesh cell $K$, we define the space $\mathbb{P}_{p}(K)$ to be the polynomials of (total) degree at most $p$ on $K$, with $\mathbb{P}_{p}(K)=\{0\}$ for $p<0$, and the space $\mathbb{P}_{p}(\partial K)$ to be continuous functions on $\partial K$ whose trace on each edge $e$ is the trace of a function from $\mathbb{P}_{p}(K)$ (equivalently, from $\mathbb{P}_{p}(\mathbb{R}^{2})$) on $e$, and we denote by $\mathbb{P}_{p}(e)$ this edge trace space. In [3], we refer to this definition of $\mathbb{P}_{p}(e)$ as its Type 2 version; the Type 1 version consists of functions on $e$ that are polynomials with respect to a natural parameter, such as arc length, in a parametrization of $e$. In order to avoid unnecessary complications in our description, we will assume that no edge is a closed contour, i.e. each edge has two distinct endpoints, and that $K$ is simply-connected. This is not a necessary constraint in practice, but allowing for even more general elements, such as those having no vertices, or those that are not simply-connected (i.e. have holes) requires using different integral equation techniques. We briefly highlight this issue in Section 4.

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded domain with Lipschitz boundary. Given a partition $\mathcal{T}=\{K\}$ of $\Omega$, we define $V_{p}(\mathcal{T})$ by

[TABLE]

where we define the space $V_{p}(K)$ as follows,

[TABLE]

The space $V_{p}(K)$ clearly contains $\mathbb{P}_{p}(K)$, but it typically contains other functions as well. A natural decomposition of $V_{p}(K)$ is $V_{p}(K)=V_{p}^{K}(K)\oplus V_{p}^{\partial K}(K)$, where

[TABLE]

The dimension of $V_{p}(K)$ is

[TABLE]

The dimension of $\mathbb{P}_{p}(\partial K)$ depends on the number and nature of the edges of $K$. If $e$ is a straight edge, $\dim\mathbb{P}_{p}(e)=p+1$, but if $e$ is not a straight edge, the dimension of $\mathbb{P}_{p}(e)$ can be as high as $\binom{p+2}{2}$, as it is when $p=1$. The dimension of $\mathbb{P}_{p}(e)$ more generally is given in the following proposition, a proof of which may be found in [60, Theorem 7.1], for example.

Proposition 2.1**.**

Suppose that $f_{m}\in\mathbb{P}_{m}(\mathbb{R}^{2})$ is an irreducible polynomial of degree $m$, and that all points $x\in e$ satisfy $f_{m}(x)=0$. It holds that $\dim\mathbb{P}_{p}(e)=\binom{p+2}{2}-\binom{p-m+2}{2}$. If $e$ does not lie on a real algebraic curve in the plane, then $\dim\mathbb{P}_{p}(e)=\binom{p+2}{2}$.

Example 2.2*.*

We consider the dimensions of the spaces $V_{p}(K)=V_{p}^{K}(K)\oplus V_{p}^{\partial K}(K)$, $p=1,2$, for each of the three mesh cells in Figure 1. We have $\dim V_{1}^{K}(K)=0$ and $\dim V_{2}^{K}(K)=1$. The continuity of functions in $\mathbb{P}_{p}(\partial K)$ implies that

[TABLE]

This formula holds for arbitrary $p$. For $p=1,2$, we have

[TABLE]

For the Half-Washer and Two-Edge Circle, the curved edges are circular arcs. For the Shuriken, the curved edges are not segments of curved conic sections (ellipses, parabolas, hyperbolas). Therefore, the dimensions of these spaces are

[TABLE]

A basis for $\mathbb{P}_{p-2}(K)$ implicitly defines a basis for $V_{p}^{K}(K)$. In Section 4, we describe how we form the associated local finite element linear systems over $V_{p}(K)$, using integral equations to get the relevant information about our basis functions. At this stage, we merely state that it is convenient to compute harmonic functions in this context. To this end, let $z$ be a point in $K$, and $\alpha=(\alpha_{1},\alpha_{2})$ be a multi-index. In [54], the authors provide an explicit formula for a polynomial $q_{\alpha}\in\mathbb{P}_{p}(K)$ satisfying $\Delta q_{\alpha}=(x-z)^{\alpha}$; see also the beginning of Section 4. A basis of $V_{p}^{K}(K)$, $\{\phi_{\alpha}^{K}\in\mathbb{P}_{p}(K):\,|\alpha|\leq p\}$, is given by $\phi_{\alpha}^{K}=\psi_{\alpha}^{K}+q_{\alpha}$, where

[TABLE]

Similarly, a basis of $\mathbb{P}_{p}(\partial K)$ naturally leads to a basis of $V^{\partial K}_{p}(K)$. Given an edge $e$ in the mesh, we describe an approach for obtaining a basis of $\mathbb{P}_{p}(e)$ that is independent of the mesh cell(s) of which it is an edge. Let $e$ have vertices $z_{0},z_{1}$. We choose a third point $z_{2}$ such that $z_{0},z_{1},z_{2}$ are the vertices of an equilateral triangle (see Figure 2)—note that $z_{2}$ typically has nothing to do with the underlying mesh $\mathcal{T}$. Given a global numbering of the vertices of the mesh, this can be done in a consistent way by choosing $z_{2}$ such that a counter-clockwise traversal of the boundary of the triangle is consistent with traversing the edge $e$ from its smaller to its larger vertex numbers. Let $\ell_{0},\ell_{1},\ell_{2}\in\mathbb{P}_{1}(\mathbb{R}^{2})$ be the three barycentric coordinates associated with these vertices. Formulas for these three functions are given by

[TABLE]

where we understand the subscripts modulo $3$ (i.e. $z_{-1}=z_{2}$ and $z_{3}=z_{0}$), and

[TABLE]

Any basis for $\mathbb{P}_{p}(\mathbb{R}^{2})$ yields a spanning set for $\mathbb{P}_{p}(e)$ by restriction, and such a basis may be expressed in terms of linear combinations of products of the barycentric coordinates. We will consider hierarchical bases expressed in this way (cf. [1, 17, 74]). For example, a hierarchical basis for $\mathbb{P}_{3}(\mathbb{R}^{2})$ is

[TABLE]

Here, we have chosen the scaling on each function so that its maximum value on the triangle, in magnitude, is $1$. In any hierarchical basis for $\mathbb{P}_{p}(\mathbb{R}^{2})$, the only functions that do not vanish at both $z_{0}$ and $z_{1}$ are $\ell_{0}$ and $\ell_{1}$. A simple consequence of this fact is that

Proposition 2.3**.**

For any edge $e$, a hierarchical basis of $\mathbb{P}_{p}(e)$ contains both $\ell_{0}$ and $\ell_{1}$.

As stated in Proposition 2.1, if $e$ lies on an algebraic curve of order $m$, we know the dimension of $\mathbb{P}_{p}(e)$. However, it may be undesirable to make this determination in practice. Regardless, we need a practical method for paring down a spanning set for $\mathbb{P}_{p}(e)$ to a basis. Let $N=\binom{p+2}{2}$, and suppose that $\{\ell_{0},\ell_{1},b_{1},\ldots,b_{N-2}\}$ is a hierarchical spanning set of $\mathbb{P}_{p}(e)$, as described above. The functions are listed in increasing order of degree. The Gram matrix $m_{ij}=\int_{e}b_{i}b_{j}\,ds$ may be used to determine the remaining basis functions (in addition to $\ell_{0},\ell_{1}$) for $\mathbb{P}_{p}(e)$. We recall that $\mathrm{rank}(M)=\dim\mathrm{span}\{b_{1},\ldots,b_{N-2}\}$ (cf. [53, Theorem 7.2.10]). A basis for $\mathrm{span}\{b_{1},\ldots,b_{N-2}\}$ consisting of some subset of these functions may be determined using a rank-revealing Cholesky decomposition of $M$ (cf. [49, 48, 47]). We state a slightly more general version of this result in the following proposition, and then provide a simple algorithm for selecting a basis of $\mathrm{span}\{b_{1},\ldots,b_{N-2}\}$, and hence of $\mathbb{P}_{p}(e)$.

Proposition 2.4**.**

Let $m_{ij}=\langle b_{j},b_{i}\rangle$, $1\leq i,j\leq n$ be the Gram matrix associated with an inner-product $\langle\cdot,\cdot\rangle$ and a list of vectors $(b_{1},\dots,b_{n})$. Let $P^{T}MP=R^{T}R$ be a rank-revealing Cholesky decomposition, where $P$ is a permutation matrix, and $R=\begin{pmatrix}R_{11}&R_{12}\end{pmatrix}$, with the $r\times r$ matrix $R_{11}$ having strictly positive entries. Then $\{b_{p(j)}:\,1\leq j\leq r\}$ is a basis for $\mathrm{span}\{b_{1},\ldots,b_{n}\}$, where $p$ is the permutation on $\{1,\ldots,n\}$ defined by $P\mathbf{e}_{j}=\mathbf{e}_{p(j)}$, and $\{\mathbf{e}_{1},\ldots,\mathbf{e}_{n}\}$ are the standard coordinate vectors.

The following algorithm is essentially Gaussian elimination with complete pivoting for positive semi-definite matrices, where the pivoting is done in place.

Algorithm 2.5**.**

Let $m_{ij}=\langle b_{j},b_{i}\rangle$, $1\leq i,j\leq n$, be the Gram matrix associated with an inner-product $\langle\cdot,\cdot\rangle$ and a list of vectors $(b_{1},\dots,b_{n})$. Upon termination of the following algorithm, $\mathrm{index}=\{p(j):\,1\leq j\leq r\}$, the indices of a basis $\{b_{p(j)}:\,1\leq j\leq r\}$ of $\mathrm{span}\{b_{1},\ldots,b_{n}\}$:

• $\mathrm{index}=\{\}$

• $k=\mathrm{arg\,max}\{m_{jj}:\,1\leq j\leq n\}$

• while $m_{kk}>0$

• $\mathrm{index}=\mathrm{index}\cup\{k\}$

• $M=M-m_{kk}^{-1}\mathbf{m}_{k}\mathbf{m}_{k}^{T}$

• $k=\mathrm{arg\,max}\{m_{jj}:\,1\leq j\leq n\}$

• end

Here, $\mathbf{m}_{k}$ is the $k^{th}$ column of the current $M$.

In practice, one replaces the condition $m_{kk}>0$ with $m_{kk}>\tau$ for some suitably small tolerance $\tau>0$. Some speed-up of this basic algorithm may be achieved by exploiting the fact that previous reduction reduction steps, $M=M-m_{kk}^{-1}\mathbf{m}_{k}\mathbf{m}_{k}^{T}$, have zeroed out the rows and columns in the index set, so these are no longer needed for further reductions. In the following example, we use $\tau=10^{-12}$ in our determination of a basis for $\mathbb{P}_{3}(e)$.

Example 2.6*.*

For any edge $e$, a hierarchical spanning set for $\mathbb{P}_{3}(e)$ is given by (10) where we have restricted the domains of these functions to $e$. Let $e$ be parameterized by $x(t)=(\cosh t,(\sinh t)/2)$, $0\leq t\leq 1$, so $e$ is part of the hyperbola $x^{2}-4y^{2}=1$. We know in advance that $\dim\mathbb{P}_{3}(e)=10-3=7$, and $\ell_{0},\ell_{1}$ will be part of our basis for $\mathbb{P}_{3}(e)$, so we must select five of the remaining eight functions,

[TABLE]

to complete our basis. Taking the functions in this order, and forming the associated Gram matrix, we determine that the indices are (given in the order computed): 4, 7, 8, 3, 5; the knowledge that we only needed five functions was not used in this computation. Therefore, our basis for $\mathbb{P}_{3}(e)$ is given by

[TABLE]

These basis functions are plotted, as functions of the parameter $t$, in Figure 2, together with the edge $e$ and associated triangle used to define the barycentric coordinates $\ell_{0},\ell_{1},\ell_{2}$.

As a matter of interest, we note that, when the reduction algorithm was used, with $\tau=10^{-12}$ as before, on the entire spanning set (10), a different basis was obtained,

[TABLE]

Remark 2.7*.*

A natural variant of Algorithm 2.5 that may be used if the diagonal entries of $M$ are all non-zero is to diagonally rescale its entries, $m_{ij}\longleftarrow m_{ij}/\sqrt{m_{ii}m_{jj}}$, before beginning the elimination loop. If we do this for Example 2.6, the resulting basis for $\mathbb{P}_{3}(e)$ is

[TABLE]

regardless of whether we use the entire spanning set (10), or remove $\{\ell_{0},\ell_{1}\}$, in constructing $M$.

Remark 2.8*.*

We describe an alternative to the barycentric coordinates (9) associated with $e$ that acts more like a local cartesian coordinate system. Using the same notation $z_{0},z_{1}$, $h$ and $R$, we define

[TABLE]

Straight-forward manipulations reveal that

[TABLE]

so it is simple to translate between coordinate systems if desired.

Having now properly defined $V_{p}(\mathcal{T})\subset H^{1}(\Omega)$, the discrete version of (1) is to find $\hat{u}\in V_{p}(\mathcal{T})\cap\mathcal{H}$ such that

[TABLE]

The intersection, $V_{p}(\mathcal{T})\cap\mathcal{H}$, ensures that we respect any homogeneous Dirichlet boundary conditions inherent in $\mathcal{H}$. For standard finite elements, as well as those defined on more general polygonal meshes, common assumptions on the data ensure that the finite element error $\|u-\hat{u}\|_{H^{1}(\Omega)}$ is controlled by interpolation error $\|u-\mathcal{I}u\|_{H^{1}(\Omega)}$, where $\mathcal{I}u\in V_{p}(\mathcal{T})$ is some appropriately defined interpolant of $u$. In the next section, we define a projection-based interpolation operator appropriate for our setting, and prove that it yields the desired approximation properties.

3. Interpolation

In this section, we describe a local interpolation scheme

[TABLE]

and establish local error estimates under stronger regularity assumptions. By construction, the local interpolation operator will define a global interpolation operator

[TABLE]

by $(\mathcal{I}v){|_{K}}=\mathcal{I}_{K}v$.

Our definition of $\mathcal{I}_{K}$ is motivated by the decomposition $V_{p}(K)=V_{p}^{K}(K)\oplus V_{p}^{\partial K}(K)$. We begin with a related decomposition of $v$ as $v=v^{K}+v^{\partial K}$, where

[TABLE]

We define $\mathcal{I}_{K}$ by an analogous decomposition $\mathcal{I}_{K}v=\mathcal{I}_{K}^{K}v+\mathcal{I}_{K}^{\partial K}v$, where $\mathcal{I}_{K}^{K}v\in V_{p}^{K}(K)$ and $\mathcal{I}_{K}^{\partial K}v\in V_{p}^{\partial K}(K)$ are given by

[TABLE]

In order to complete this definition, we must define $q^{K}\in\mathbb{P}_{p-2}(K)$ and $q^{\partial K}\in\mathbb{P}_{p}(\partial K)$. We define $q^{K}$ by

[TABLE]

We define $q^{\partial K}$ by defining it on each edge of $\partial K$. For a non-trivial open subset $\Gamma\subset\partial K$, we use the inner-product

[TABLE]

with $\|\cdot\|_{H^{1/2}(\Gamma)}$ as the associated norm. Below, we take $\Gamma$ to be either the entire boundary, $\partial K$, or a single edge, $e$. Fix an edge $e$ of $\partial K$, having endpoints $z_{0},z_{1}$, and let $\mathbb{P}_{p,0}(e)=\{w\in\mathbb{P}_{p}(e):\,w(z_{0})=w(z_{1})=0\}$. We define $q_{e}\in\mathbb{P}_{m}(e)$ by the conditions

[TABLE]

Finally, $q^{\partial K}$ is defined by $(q^{\partial K}){|_{e}}=q_{e}$.

In several places below, it will be convenient to use the following basic result. Suppose that $\psi,\phi\in H^{1}(K)$, and $\Delta\psi=0$ in $K$ and $\phi=0$ on $\partial K$. Then $\int_{K}\nabla\psi\cdot\nabla\phi\,dx=0$, so $|\psi+\phi|_{H^{1}(K)}^{2}=|\psi|_{H^{1}(K)}^{2}+|\phi|_{H^{1}(K)}^{2}$. For example, we have

[TABLE]

and we consider both contributions to the interpolation error in turn. In the proofs below, we use $c$ as a constant that may vary from one appearance to the next. Throughout, $h_{K}$ denotes the diameter of $K$. We first consider $v^{K}-\mathcal{I}_{K}^{K}v$.

Proposition 3.1**.**

Suppose that $v\in H^{1}(K)$ and $\Delta v\in H^{p-1}(K)$ for some $p\geq 1$. There is a scale-invariant constant $c=c(p,K)>0$ for which

[TABLE]

Proof.

It holds that

[TABLE]

Since $v^{K}-\mathcal{I}_{K}^{K}v\in H^{1}_{0}(K)$, the Poincaré-Friedrichs Inequality, $\|w\|_{L^{2}(K)}\leq h_{K}|w|_{K}$ for $w\in H^{1}_{0}(K)$, ensures that

[TABLE]

The estimate $|v^{K}-\mathcal{I}_{K}^{K}v|_{H^{1}(K)}\leq ch_{K}^{p}|\Delta v|_{H^{p-1}(K)}$ follows from this by applying the Bramble-Hilbert Lemma. The $L^{2}(K)$ norm result follows from this by applying the Poincaré-Friedrichs Inequality again. ∎

Remark 3.2*.*

If $K$ is convex, $h_{K}$ can be replaced by $h_{K}/\pi$ in (21) (cf. [66]). Furthermore, for convex $K$, the dependence on $K$ of the constant $c(p,K)$ coming from the Bramble-Hilbert Lemma in Proposition 3.1 can be removed, and for non-convex domains that are star-shaped with respect to a point, ball, or more general subdomain, various estimates of how $c(p,K)$ depends on the shape of $K$ have been established [76, 35, 30, 29].

For our analysis of $v^{\partial K}-\mathcal{I}_{K}^{\partial K}v$ the following result will be useful.

Proposition 3.3**.**

If $v\in H^{2}(K)$ and $\Delta v=0$ in $K$, there is a scale-invariant constant $c=c(K)$ for which

[TABLE]

Proof.

For $w\in H^{2}(K)$, we have $\|w\|_{L^{\infty}(K)}\leq c\|w\|_{H^{2}(K)}$ by a Sobolev embedding result. A standard scaling argument then yields

[TABLE]

where $c=c(K)$ is scale-invariant. Now (22) follows from the fact that harmonic functions on $K$ attain their extrema on $\partial K$, so, if $v,w\in H^{2}(K)$ have the same Dirichlet trace on $\partial K$, and $v$ is harmonic on $K$, then $\|v\|_{L^{\infty}(K)}\leq\|w\|_{L^{\infty}(K)}$. ∎

Remark 3.4*.*

Since we are working in $\mathbb{R}^{2}$, $H^{1+s}(K)$ is continuously imbedded in $C(\overline{K})$ for any $s\in(0,1)$. Therefore, Proposition 3.3 is readily generalized to such spaces, with the obvious bound

[TABLE]

Typical assumptions on the domain $\Omega$ and the data for the problem guarantee that $u\in H^{1+s}(\Omega)$ for some $s>0$ (cf.  [46, 45, 87]).

We now consider the term $|v^{\partial K}-\mathcal{I}_{K}^{\partial K}v|_{1,K}$. Let $e$ be an edge of $\partial K$, with endpoints $z_{0},z_{1}$. We begin with a further decomposition of $q_{e}$, namely $q_{e}=q_{e,1}+q_{e,0}$, where $q_{e,1}=v(z_{0})\ell_{0}+v(z_{1})\ell_{1}\in\mathbb{P}_{1}(e)$ and $q_{e,0}=q_{e}-q_{e,1}\in\mathbb{P}_{p,0}(e)$. This induces a natural decomposition of $q^{\partial K}$, $q^{\partial K}=q^{\partial K}_{1}+q^{\partial K}_{0}$, where $q^{\partial K}_{1}\in\mathbb{P}_{1}(\partial K)$ satisfies $q^{\partial K}_{1}(z)=v(z)$ at each vertex $z$ of $K$, and $q^{\partial K}_{0}$ vanishes at the vertices.

Proposition 3.5**.**

Suppose that $v\in H^{p+1}(K)$ for some $p\geq 1$. There is a scale-invariant constant $c=c(p,K)$ for which

[TABLE]

Proof.

We decompose $\mathcal{I}^{\partial K}_{K}v$ as $\mathcal{I}^{\partial K}_{K}v=w_{1}+w_{0}$, where

[TABLE]

It follows that $|\mathcal{I}^{\partial K}_{K}v|_{H^{1}(K)}\leq|w_{1}|_{H^{1}(K)}+|w_{0}|_{H^{1}(K)}$.

We denote the set of vertices of $K$ by $\mathcal{V}(K)$, and the set of edges of $K$ by $\mathcal{E}(K)$. For $z\in\mathcal{V}(K)$, we define $\ell_{z}\in\mathbb{P}_{1}(\partial K)$ as follows: if $e$ is not adjacent to $z$, then $\ell_{z}$ vanishes on $e$, and if $e$ is adjacent to $z$, then $\ell_{z}=\ell_{j}$ on $e$, where $z=z_{j}$ for one of the endpoints $z_{0},z_{1}$ of $e$. Let $\phi_{z}$ be the harmonic function on $K$ whose Dirichlet trace on $\partial K$ is $\ell_{z}$. It follows that $w_{1}=\sum_{z\in\mathcal{V}(K)}v(z)\phi_{z}$, so

[TABLE]

where we have used (22) in the final inequality. A similar argument shows that $\|w_{1}\|_{L^{2}(K)}\leq c(h_{K}^{2}|v|_{H^{2}(K)}+h_{K}|v|_{H^{1}(K)}+\|v\|_{L^{2}(K)})$.

From (19), we see that $(q_{e,0},q_{e,0})_{H^{1/2}(e)}=(v-q_{e,1},q_{e,0})_{H^{1/2}(e)}$, so $\|q_{e,0}\|_{H^{1/2}(e)}\leq\|v-q_{e,1}\|_{H^{1/2}(e)}$, for each edge $e$. Now,

[TABLE]

Here we have used applied the trace inequality $\|\partial w_{0}/\partial n\|_{H^{-1/2}(\partial K)}\leq c(|w_{0}|_{H^{1}(K)}+\|\Delta w_{0}\|_{L^{2}(K)})=c|w_{0}|_{H^{1}(K)}$ (cf. [70, Theorem A.33]), where $c=c(K)$ is scale-invariant. From this it follows that

[TABLE]

The second inequality holds because $q_{0}$ vanishes at the vertices, see Remark 3.6. At this stage, $c=c(p,K)$.

Another standard trace inequality ensures that $\|v-q_{1}^{\partial K}\|_{H^{1/2}(\partial K)}\leq c\left(|v-w_{1}|_{H^{1}(K)}+h_{K}^{-1}\|v-w_{1}\|_{L^{2}(K)}\right)$. Combining this with our estimates above, we obtain

[TABLE]

and it follows, by applying the estimates for $|w_{0}|_{H^{1}(K)}$ and $|w_{1}|_{H^{1}(K)}$, that

[TABLE]

A standard inverse inequality, and the fact that our interpolation scheme preserves constants, yields the obvious analogue in $L^{2}(K)$,

[TABLE]

Now, let $\phi\in V_{p}(K)$ and decompose it as $\phi=\phi^{K}+\phi^{\partial K}$, with $\phi^{K}\in V_{p}^{K}(K)$ and $\phi^{\partial K}\in V_{p}^{\partial K}(K)$. We have $|v-\phi|_{H^{1}(K)}^{2}=|v^{K}-\phi^{K}|_{H^{1}(K)}^{2}+|v^{\partial K}-\phi^{\partial K}|_{H^{1}(K)}^{2}$, so $|v^{\partial K}-\phi^{\partial K}|_{H^{1}(K)}\leq|v-\phi|_{H^{1}(K)}$. Noting that $\mathcal{I}^{\partial K}_{K}\phi=\phi^{\partial K}$, and applying (24) to $v-\phi$, we see that

[TABLE]

Since $\mathbb{P}_{p}(K)\subset V_{p}(K)$, the Bramble-Hilbert Lemma now implies that $|v^{\partial K}-\mathcal{I}^{\partial K}_{K}v|_{H^{1}(K)}\leq ch_{K}^{p}|v|_{H^{p+1}(K)}$, as claimed. The result for the $L^{2}(K)$ norm follows the same pattern, but we briefly lay out the argument anyway. It holds that

[TABLE]

It remains to estimate $\|v^{\partial K}-\phi^{\partial K}\|_{L^{2}(K)}$, for which we have

[TABLE]

Combining this with our previous estimate yields,

[TABLE]

and Bramble-Hilbert Lemma completes the argument. ∎

Remark 3.6*.*

The claim that $\|q_{0}^{\partial K}\|_{H^{1/2}(\partial K)}^{2}\leq c\sum_{e\in\mathcal{E}(K)}\|q_{0}^{\partial K}\|_{H^{1/2}(e)}^{2}$ in the proof of Proposition 3.5 requires further comment. Superficially, this holds because both quantities are (squares of) norms on the finite dimensional vector space $\mathbb{P}_{p,0}(\partial K)=\{w\in P_{p}(\partial K):\,w(z)=0\mbox{ for all }z\in\mathcal{V}(K)\}$. Although such an argument allows for the dependence of $c$ on $\dim\mathbb{P}_{p,0}(\partial K)$ (hence on $p$), we want to ensure that $c$ is scale-invariant. For that, we look a little closer at the norms. We have $\|w\|_{H^{1/2}(\Gamma)}^{2}=\|w\|_{L^{2}(\Gamma)}^{2}+|w|_{H^{1/2}(\Gamma)}^{2}$, where

[TABLE]

As suggested by the notation, $|\cdot|_{H^{1/2}(\Gamma)}$ is generally a semi-norm, with constant functions as its kernel. However, for $w\in\mathbb{P}_{p,0}(\partial K)$, both $|w|_{H^{1/2}(\partial K)}$ and $\left(\sum_{e\in\mathcal{E}(K)}|w|_{H^{1/2}(e)}^{2}\right)^{1/2}$ are norms, so there is a constant $c$ such that $|w|_{H^{1/2}(\partial K)}\leq c\left(\sum_{e\in\mathcal{E}(K)}|w|_{H^{1/2}(e)}^{2}\right)^{1/2}$. Since both norms in this inequality are scale invariant, so is $c$. Since $\|w\|_{L^{2}(\partial K)}^{2}=\sum_{e\in\mathcal{E}(K)}\|w\|_{L^{2}(e)}^{2}$, we have the result that was claimed.

We briefly mention two earlier contributions that have considered some of the same issues that we do here concerning working with the $H^{1/2}$ norm on all versus individual parts of the boundary of a mesh cell or a polyhedral subdomain , but in the context of standard finite element meshes. The first is [22], and it concerns domain decomposition-type preconditioners for linear solvers. Though we were unable to use the results of Section 3 in that paper related to localization of the $H^{1/2}$ norm our context, they provide the first discussion and treatment of this issue of which we are aware in the finite element literature. The second contribution is [31], in which the authors set forth projection-based interpolation schemes that are conforming in $H^{1}$, $H(\mathrm{curl})$ and $H(\mathrm{div})$ spaces. Our interpolation scheme is also projection based, but because their results were for standard element shapes, they could not be readily applied in our context.

Remark 3.7*.*

The proof of Proposition 3.5 revealed that

[TABLE]

In fact, by the same reasoning as discussed in Remark 3.4, we have the expected versions for fractional order spaces as well, for $s\in(0,1]$,

[TABLE]

Combining Propositions 3.1 and 3.5, we obtain our key interpolation error result,

Theorem 3.8**.**

Suppose that $v\in H^{p+1}(K)$ for some $p\geq 1$. There is a scale-invariant constant $c=c(p,K)$ for which

[TABLE]

Once a proper notion of “shape regularity” is determined for families of meshes $\{\mathcal{T}_{h}\}$ consisting of curvilinear polygons, a result such as $|v-\mathcal{I}v|_{H^{1}(\Omega)}\leq ch^{p}|v|_{H^{p+1}(\Omega)}$, where $c=c(p)$ and $h=\max\{h_{K}:\,K\in\mathcal{T}_{h}\}$, follows immediately. A meaningful analysis of how the constant $c(p,K)$ in Theorem 3.8 depends on $p$ and the geometric features of $K$ is beyond the scope of the present work. One might expect measures such as a “chunkiness parameter” (a natural generalization of aspect ratio, cf. [23, Defintion 4.2.16]), the number of edges, the curvature of edges, and the length of edges with respect to the element diameter, to play an important role in determining the dependence of $c=c(K)$ on element geometry. Indeed, the number of edges of $K$ clearly arises in the proof of Proposition 3.5, when we bound $|w_{1}|_{H^{1}(K)}$ using a sum of seminorms of functions associated the vertices; see also Remark 3.6, in which a sum over edges is used. In contrast, the special analysis given for the L-shaped elements of Example 5.3, which have fixed size but increasing number of edges as the mesh is “refined”, provides an example in which neither the number of edges nor their relation to the diameter of the element have any bearing on the associated interpolation constant. In each of the other examples in Section 5, the maximal curvature of edges grows without bound as the diameters of the elements shrink, and this has no apparent negative effect on the convergence of the discretization error, which suggests that edge curvature may not ultimately play such an important role in interpolation error analysis either. Additionally, some of the families of meshes in Section 5 consist entirely of elements that are not star-shaped with respect to any ball, in which case discussion of a chunkiness parameter is either meaningless, or would have to take on a different form if it were to be applicable at all. In summary, a thorough analysis of how local interpolation error depends (or does not depend) on geometric features of elements is needed. Further extensions of our interpolation error analysis of interest include:

1. (a)

Estimates that directly involve both the element diameter $h_{K}$ and a local “polynomial degree” $p_{K}$, in the manner of standard $hp$-finite element analysis. 2. (b)

Estimates that exploit the fact that $V_{p}(K)$ is a richer space than $\mathbb{P}_{p}(K)$, often containing singular functions that may allow similar convergence results for interpolation under weaker regularity assumptions on $v$, as suggested by Remarks 3.4 and 3.7.

We plan to pursue these extensions in subsequent work.

4. Computing with Curved Trefftz Finite Elements

We recall that the functions that we wish to compute in $V_{p}(K)$ satisfy one of two types of equations:

[TABLE]

where $f\in\mathbb{P}_{p-2}(K)$ and $g\in\mathbb{P}_{p}(\partial K)$. The first type of equation is readily converted to the second type as follows. Given $f\in\mathbb{P}_{p-2}(K)$, one can explicitly construct a $\hat{f}\in\mathbb{P}_{p}(K)$ such that $\Delta\hat{f}=f$. With such a function in hand, the first type of problem is reduced to finding $\hat{v}$ satisfying $\Delta\hat{v}=0$ in $K$ and $\hat{v}=-\hat{f}$ on $\partial K$. Then $v=\hat{v}+\hat{f}$ satisfies the first problem. In [54, Theorem 2], the authors show that, if $p$ is a homogeneous polynomial of degree $j$, then the polynomial $q$ of degree $j+2$ given by

[TABLE]

where $[j/2]$ denotes the integer part of $j/2$, satisfies $\Delta q=p$. Having reduced either type of problem to the computation of a harmonic function with piecewise smooth boundary data, we may now employ any number of boundary integral equation techniques to compute such functions. One such technique is to use Boundary Element Methods for first-kind integral equations, as is done in BEM-FEM, to directly compute the outward normal derivative $\partial v/\partial n$; interior point values are computed from layer potentials, as needed, for quadrature approximation of the element stiffness matrix. The limitations in extending this kind of approach to curved element boundaries in a natural way was one of the reasons that we opted for Nyström discretizations in [3]. In that work, we employed second-kind integral equations, which do not directly yield $\partial v/\partial n$, but offered greater flexibility in other areas that offset this downside.

Before describing the approach we use in the current work, we recall the types of integrals we must compute in order to form the finite element stiffness matrix. They include integrals of the following forms,

[TABLE]

where $v,w\in V_{p}(K)$. It is sometimes advantageous to use integration-by-parts on the first of these integrals, yielding integrals of the forms

[TABLE]

The benefits of such an approach become clear when $A$ is a constant scalar on $K$, in which case the two integrals above simplify to

[TABLE]

We note that $\Delta v\in\mathbb{P}_{p-2}(K)$ and $w\in\mathbb{P}_{p}(\partial K)$. Further simplifications occur when $v\in V_{p}^{\partial K}(K)$ or $w\in V_{p}^{K}(K)$—at least one of these two integrals vanishes. We see then that, in quadrature approximations of these kinds of integrals, we should have access to function values and derivatives (up to second partials) of functions in $V_{p}(K)$ in the interior of $K$, and normal derivatives of such functions on $\partial K$—function values and tangential derivatives of such functions on $\partial K$ are straightforward.

With these goals in mind, in [65] we developed an approach that delivers each of these quantities efficiently and with very high accuracy, while performing all computations on the boundary $\partial K$. The method is based on the fact that, on simply-connected domains $K\subset\mathbb{R}^{2}$, for each harmonic function $u$, there is a family of harmonic conjugates that differ from each other only by additive constants. We recall that $v$ is a harmonic conjugate of $u$ on $K$ when $\Delta v=0$ in $K$ and $R\nabla v=\nabla u$ in $K$, where the matrix $R$ rotates vectors clockwise by $\pi/2$. Such a pair of harmonic functions satisfy the Cauchy-Reimann equations, and thus can be taken as the real and imaginary parts of a complex analytic function in $K$. More precisely, making the natural identification between $z=x_{1}+\mathfrak{i}x_{2}\in\mathbb{C}$ and $x=(x_{1},x_{2})\in K$, the function defined by $w(z)=u(x)+\mathfrak{i}v(x)$ is analytic in $K$. Given both $u$ and $v$ on the boundary $\partial K$, the value $w$ and its derivatives at points inside $K$ can be obtained via Cauchy’s integral formula,

[TABLE]

and the desired $k$th partial derivatives of $u$ (or $v$) can be extracted from the real and imaginary parts of $w^{(k)}(z)$. For example, $w^{\prime}(z)=u_{x_{1}}(x)-\mathfrak{i}u_{x_{2}}(x)$, where $u_{x_{j}}$ denotes the partial derivative of $u$ in its $j$th argument. For the second partials, we have $w^{\prime\prime}(z)=u_{x_{1}x_{1}}(x)-\mathfrak{i}u_{x_{1}x_{2}}(x)$, with $u_{x_{2}x_{1}}(x)=u_{x_{1}x_{2}}(x)$ and $u_{x_{2}x_{2}}(x)=-u_{x_{1}x_{1}}(x)$. Furthermore, the orthogonality of $\nabla u$ and $\nabla v$ in $K$ ensures the following relationship between the normal and tangential derivatives of $u$ and $v$ on $\partial K$.

[TABLE]

where $\partial v/\partial t$ denotes the tangential derivative along $\partial K$ in the counter-clockwise direction. With these relationships, we see that it is possible to compute the normal derivative of one harmonic function as the (much more convenient) tangential derivative of a harmonic conjugate. Therefore, this general approach allows us to compute all of the quantities of interest related to our harmonic function while performing all computations on the boundary $\partial K$, as claimed.

Given the piecewise smooth boundary Dirichlet data $g$ of a harmonic function $u$, any harmonic conjugate $v$ satisfies the complementary Neumann problem

[TABLE]

As stated earlier, solutions of (30) are only unique up to additive constants, and we fix a particular member by specifying that $\int_{\partial K}v\,ds=0$. The trace of $v$ on $\partial K$ is computed as the solution of the following second-kind integral equation,

[TABLE]

where $\Phi(x,y)=-(2\pi)^{-1}\ln|x-y|$ is the fundamental solution of Laplace’s equation. The addition of $1$ to the integral kernel $\partial\Phi(x,y)/\partial n(y)$ above, ensures that (31) is well-posed by enforcing that $\int_{\partial K}v\,ds=0$. Since $K$ will typically have corners, the integral equation must be modified in their vicinity if we are to understand the equation pointwise. More specifically, if $x_{c}\in\partial K$ is a corner point, and $x\in\partial K$ is not a corner point, then

[TABLE]

The case of multiple corners is handled similarly. In the present work, as in [65], we solve (31) via a Nyström discretization. Key to the practical success of this approach is the choice of quadrature schemes that are well-suited for the types of singularities present in the integrands. The interested reader may find the details in that paper [65]. Having computed $v$ on $\partial K$, we now have access to the quantities of interest for $u$ as described above.

Remark 4.1* (Multiply connected mesh cells).*

If $K$ is not simply connected, the existence of harmonic conjugate pairs is not guaranteed, so the approach described above cannot be used. In such cases, one can employ different integral equation techniques to efficiently and accurately compute the quantities necessary for assembling local stiffness matrices. We mention the contribution [44] (and references therein) in this regard. The discussion in that work assumes smooth boundaries (at least $C^{2}$), so their approach must be modified in order to handle mesh cells having corners. We intend to pursue this in subsequent work.

5. Numerical Experiments

The experiments in this section illustrate the linear convergence rate indicated by Theorem 3.8 ($p=1$) for $|u-\hat{u}|_{H^{1}(\Omega)}$ on simple model problems which nonetheless illustrate the theoretical claims are achieved in practical computations. In the first set of experiments, three increasingly complex families of meshes are used for the same problem on the unit square. In the second set of experiments, we explore the effects of “nearly straight” edges on convergence and conditioning for Type 1 and Type 2 elements. For the final set of experiments, we consider a problem for which the exact solution is known to have a singularity due to a non-convex corner of the domain, and use it to illustrate the approximation power of locally singular functions in our finite element spaces.

Example 5.1* (Three Curved Mesh Families).*

Let $\Omega=(0,1)\times(0,1)$, and suppose that $u$ satisfies

[TABLE]

Series representations of $u$ and $|u|_{H^{1}(\Omega)}^{2}$ are

[TABLE]

We approximate $u$ by the finite element solution $\hat{u}\in V_{1}(\mathcal{T})$ on three different families of meshes, each indexed by a mesh parameter $r$ that is inversely proportional to the characteristic diameter of its mesh cells, see Figure 3 for the case $r=16$ of each, together with the corresponding computed finite element solution $\hat{u}\in V_{1}(\mathcal{T})$. Since the cells are of uniform size, $r$ is the number of cells touching each edge of $\partial\Omega$.

We refer to the first family of meshes as the Shuriken meshes, because it consists of shuriken elements, as seen in Figure 1), which are naturally modified at the boundary to properly fit it. There are three types of elements in this case, the corner elements, edge elements and interior elements. The second family is called the Pegboard meshes, and it consists of two types of elements, the half-washers and two-edge circles, as seen in Figure 1). The third family is called the Jigsaw meshes, and it has four different types of elements: corner pieces, two different types of edge pieces, and interior pieces. Of the nine different types of elements that appear in each of these families, only the two-edge circles are convex. In fact, none of the other types of elements are even star-shaped. Furthermore, each of the jigsaw elements have at least two non-convex corners, which implies that the local space $V_{1}(K)$ for such an element will include functions that are singular, i.e. not in $H^{2}(K)$.

The discretization error $|u-\hat{u}|_{H^{1}(\Omega)}$ satisfies $|u-\hat{u}|_{H^{1}(\Omega)}^{2}=|u|_{H^{1}(\Omega)}^{2}-|\hat{u}|_{H^{1}(\Omega)}^{2}$, making it straight-forward to compute once $\hat{u}\in V_{1}(\mathcal{T})$ has been computed. These errors, and ratios of consecutive errors, are given in Table 1 for each of the families, demonstrating the expected linear convergence. As an interesting comparison, we also include the errors and ratios for the discretizations that would arise had we chosen the Type 1 definition of $\mathbb{P}_{1}(e)$ for the shuriken elements. This choice leads to local spaces having dimension $4$ for mesh cells not touching the boundary, in contrast to the dimension $8$ local spaces using the Type 2 definition of $\mathbb{P}_{1}(e)$. We recall that $\mathbb{P}_{1}(K)\not\subset V_{1}(K)$ for Type 1 elements (unless $K$ is a straight-edge polygon), and we see that there is no convergence at all in $H^{1}$ in this case! In all three cases, optimal order convergence, $|u-\hat{u}|_{H^{1}(\Omega)}=\mathcal{O}(h^{2})$, is obtained when Type 2 elements are used.

Example 5.2* (Perturbed Triangle Mesh).*

Let $\Omega$ and $u$ be as in the previous example. We again approximate $u$ by its finite element solution $\hat{u}\in V_{1}(\mathcal{T})$, where the mesh $\mathcal{T}$ consists perturbed triangular elements, as shown in Figure 4, each of which has one curved edge. Reference elements, one convex and the other non-convex, are obtained by splitting the unit square $(0,1)\times(0,1)$ using a circular arc through the vertices $(1,0)$ and $(0,1)$ whose center is $(-a,-a)$, for some $a>0$. The radius of curvature of this curved edge is $\sqrt{a^{2}+(1+a)^{2}}$, and it approaches a straight line as $a$ increases. More specifically, the maximum distance between a point on the curved edge and the closest point to it on the line between $(1,0)$ and $(0,1)$ is $(\sqrt{(2a+1)^{2}+1}-(2a+1))/\sqrt{2}$, which behaves like $1/(4a\sqrt{2})$ as $a\to\infty$. For the corresponding finite element meshes, these reference element pairs are scaled so that their straight edges have length $1/r$.

In Figure 5, we report both the errors $|u-\hat{u}|_{H^{1}(\Omega)}$ and the spectral condition numbers $\kappa_{2}(A)$ of the stiffness matrices for Type 1 and Type 2 elements as the meshes is refined, for several values of $a$. For both types of elements, $A$ is diagonally rescaled, $a_{ij}\longleftarrow a_{ij}/\sqrt{a_{ii}a_{jj}}$, before computing the condition numbers. As expected, the Type 2 elements exhibit optimal order convergence throughout the refinements. For Type 1 elements, the convergence curves improve as the $a$ is increased, in the sense that they stay roughly parallel to their Type 2 counterparts through more levels of refinement, but the convergence curves for Type 1 elements eventually level off, indicating a threshold beyond which the error does not decrease. The condition number plots for the Type 1 and Type 2 elements provide a complementary comparison, for which the Type 2 elements yield condition numbers that are eventually close to, and grow at the same rate as, those of the Type 1 elements, but may be significantly larger than their Type 1 counterparts for coarser meshes when the curved edges are nearly straight. The condition numbers for Type 1 elements grow like $r^{2}$ as the mesh is refined (i.e. like $\dim(V)$), which is accordance with standard linear (bilinear) elements on triangular (rectangular) meshes. We observe, based on computations done for $a=1/4,1,4,\ldots,4096$, that the condition numbers on the coarsest meshes ($r=4$) for Type 2 elements appear to grow quadratically in $a$. We also observe an apparent correlation between when the convergence curves for Type 1 elements tend to level off and when the condition numbers for Type 2 elements transition from a relatively flat phase to behaving like their Type 1 counterparts. It is a topic of future investigation to better understand how element shapes and “polynomial orders” $p$ affect convergence and conditioning for both types of elements.

Example 5.3* (L-Shaped Domain).*

For our final set of experiments, we again consider the problem

[TABLE]

but on the (rotated) L-shaped domain $\Omega=(-1,1)\times(-1,1)\setminus[0,1]\times[-1,0]$ (see Figure 6). This solution, though not known explicitly, is known to have a singularity at the origin, behaving asymptotically like $|x|^{2/3}$ near the origin (cf. [87, 46, 55]). In standard finite element computations the efficient approximation of such a singular solution would be achieved by targeted refinement of mesh cells toward the singular point that is either guided by local error indicators (computed a posteriori) or by specific knowledge of the local singular behavior to determine an a priori mesh grading strategy. A head-to-head empirical comparison of these two types of refinement strategies is provided in [58]. Others have sought to address the issue of singularities by augmenting standard polynomial finite elements with (local) enrichment functions having the types of singularities expected of the solution based on a priori knowledge (cf. [37], [61] XFEM, [73] GFEM).

Our approach for this problem is different. We use the fact that, if a mesh cell $K$ has a non-convex corner, then the local space $V_{1}(K)$ automatically contains functions having the correct type of singularity for that geometry. Remark 3.7 suggests improved approximation power for interpolation of functions having the same kind of singularities, and optimal order convergence was demonstrated in [3] for interpolation error in $L^{2}(\Omega)$ of a harmonic function having an $|x|^{2/3}$-type singularity on precisely the kinds of meshes shown in Figure 6. That work did not, however, consider interpolation error or discretization error in $H^{1}(\Omega)$ for such a problem. As before, we use a parameter $r$ to describe the meshes in this family. The $r$th mesh in this family, $\mathcal{T}_{r}$, consists of one L-shaped element, $K_{L}=(-1/3,1/3)\times(-1/3,1/3)\setminus[0,1/3]\times[1/3,0]$, and $24r^{2}$ congruent squares of size $(3r)^{-1}\times(3r)^{-1}$, see Figure 6. We note that there are $r$ squares touching each of the short edges of $\partial K\setminus\partial\Omega$, and $\partial K_{L}$ has $6r+2$ vertices. Although none of the edges in these meshes are curved, the optimal convergence rates enabled by the single L-shaped element illustrates how a result like (26) might be used to prove what is empirically observed. For this example, we provide such an analysis.

As with Example 5.1, we use a highly accurate approximation of $|u|_{H^{1}(\Omega)}^{2}$, together with the identity $|u-\hat{u}|_{H^{1}(\Omega)}^{2}=|u|_{H^{1}(\Omega)}^{2}-|\hat{u}|_{H^{1}(\Omega)}^{2}$. In the previous example, we used a Fourier expansion to obtain our approximation of $|u|_{H^{1}(\Omega)}^{2}$. Here, we use the techniques developed in [65]. Letting $w=-|x|^{2}/4$, and recognizing that $-\Delta w=1$, we have $u=v+w$, where $\Delta v=0$ in $\Omega$ and $v=-w$ on $\partial\Omega$, and it follows that

[TABLE]

The integral $\int_{\partial\Omega}\frac{\partial v}{\partial n}\,v\,ds$ is approximated using the techniques from [65].

For the square elements $K$ in $\mathcal{T}_{r}$, $V_{1}(K)$ consists of the standard bilinear finite elements. The element stiffness matrices for these elements remain that same (up to symmetric permutation) for all meshes,

[TABLE]

The local space $V_{1}(K_{L})$ for single L-shaped element $K_{L}$ changes from mesh to mesh, so its element stiffness matrix $A_{L}$ must be recomputed on each mesh. The number of rows/columns of $A_{L}$ on $\mathcal{T}_{r}$ is $6r+2$.

Revisiting the interpolation identity (20) in our present context, we have

[TABLE]

because $V_{1}^{K}(K)=\{0\}$. Since $-\Delta u^{K}=1$ in $K$ and $u^{K}=0$ on $\partial K$, we have

[TABLE]

For all of the square elements $K$, which shrink as $r$ increases, the term $|u^{K}|_{H^{1}(K)}$ is not problematic. However, the L-shaped element $K=K_{L}$ does not shrink as $r$ increases so the estimate $|u^{K}|_{H^{1}(K)}\leq|K|^{1/2}h_{K}$ gives no guarantee of convergence at all, much less at the optimal rate. In fact, if we approximate $u$ by $\hat{u}\in V_{1}(\mathcal{T}_{r})$ for this family of meshes, we do not get convergence in $H^{1}(\Omega)$!

This issue is simple to fix, however, and the remedy we now describe is suggestive of a more general principle that we aim to explore in detail in subsequent work. Because $-\Delta u^{K}=-\Delta u=1$ on $K=K_{L}$, we include the interior bubble function $\phi\in V_{2}^{K}(K_{L})$ satisfying $-\Delta\phi=1$ in $K_{L}$ and $\phi=0$ on $\partial K_{L}$. The necessary quantities associated with $\phi$ can be computed in the same manner as described (33) and its paragraph. We take $\widetilde{V}_{1}(K_{L})=\mathrm{span}(V_{1}(K_{L})\cup\{\phi\})$ and $\widetilde{V}_{1}(\mathcal{T}_{r})=\mathrm{span}(V_{1}(\mathcal{T}_{r})\cup\{\phi\})$. We recall that $\int_{\Omega}\nabla\psi\cdot\nabla\phi\,dx=\int_{K_{L}}\nabla\psi\cdot\nabla\phi\,dx=0$ for all $\psi\in V_{1}(\mathcal{T}_{r})$, so adding this function does not increase the cost of assembling and solving the necessary linear system.

With this enrichment by $\phi$, we have, on $K=K_{L}$,

[TABLE]

Since $u^{\partial K}-\mathcal{I}_{K}^{\partial K}u$ is harmonic on $K_{L}$, Dirchlet’s principle ensures that, on $K=K_{L}$,

[TABLE]

This can be estimated by a technique similar in spirit to that given in [41, Theorem 4.1]. The argument involves creating a (fictitious) sub-triangulation of $K_{L}$, taking $\phi$ to be the piecewise linear interpolant of $u^{\partial K}$ (or $u$) on this sub-triangulation, and using standard interpolation error estimates. However, unlike the estimate in [41, Theorem 4.1], which assumes $H^{2}$ regularity of $u^{\partial K}$, we use geometrically graded sub-triangulations, as pictured in Figure 7, and use estimates from [7, 57] to deduce that $|u^{\partial K}-\mathcal{I}_{K}^{\partial K}u|_{H^{1}(K)}\leq C_{K}r^{-1}$ on $K=K_{L}$, where $C_{K}$ depends only on the mesh grading parameter and the norm of $u^{\partial K}$ (or $u$) in an appropriate weighted Sobolev space. Combining this with the simple estimates from Theorem 3.8 for the square elements, and we see that $|u-\mathcal{I}u|_{H^{1}(\Omega)}\leq Cr^{-1}$, which is the rate of convergence observed in Table 2. We emphasize that sub-triangulations of $K_{L}$ are purely for the purpose of this interpolation error estimate, and are not used at any point in the actual computations.

The number of rows and columns (and non-zeros) of the global stiffness matrix $A$ grows quadratically with $r$, and the number of rows and columns of the (dense) submatrix $A_{L}$ corresponding to $V_{1}(K_{L})$ grows linearly with $r$. We comment briefly on the condition numbers of these matrices, which are reported in Table 2. As with standard (low-order) finite elements, the condition number of $A$ grows linearly with $\dim(\widetilde{V}_{1}(\mathcal{T}_{r}))$, or quadratically with $r$. In contrast, the growth of the condition number of $A_{L}$ seems to be leveling off—it is certainly not growing quadratically, or even linearly, with $r$. These issues of conditioning will be explored further in subsequent work.

In [3, Example 4.4] we compared interpolation errors in $L^{2}(\Omega)$ for the harmonic function $u=r^{2/3}\sin(2(\theta-\pi/2)/3)$ on $\Omega=(-1,1)\times(-1,1)\setminus[0,1]\times[0,1]$, a rotated version of the $\Omega$ used here, on three different families of meshes. As noted above, one of these families of meshes was the one used here, and it led to optimal order convergence. The two other families yielded sub-optimal convergence at a theoretically predicted rate. One of these families of meshes consisted solely of congruent squares, and the other family was the same except right near the corner, where it had a single small L-shaped cell obtained by merging three of these squares. Although the local space on this L-shaped element could approximate the singular function at the optimal rate, the neighboring square elements, which got increasingly closer to the singularity on finer meshes, could not, so the overall convergence was spoiled. This motivates our choice to keep the L-shaped element of fixed size, as we have here. The particular size of this element is not crucial to the overall asymptotic behavior of convergence. In fact, we could have chosen $K_{L}=\Omega$ in this case, and just solved the problem using integral equation techniques, as we did above for (33). The point of using the kinds of meshes that we did here is to demonstrate that they can offer a feasible alternative to more traditional refinement techniques.

6. Conclusions

We have provided analysis and a practical low-order realization of a novel finite element method on meshes consisting of quite general curvilinear polygons. Allowing for such curved elements introduces both theoretical and computational challenges, including the proper definition and treatment of polynomial spaces defined on curves, determining an appropriate interpolation operator and obtaining meaningful error estimates, and efficiently computing with the implicitly-defined basis functions. Concerning the first of these challenges, we described and demonstrated simple methods for constructing a spanning set for a polynomial space on an edge, and then pairing it down to a basis. Concerning the second challenge, we proved local interpolation estimates in $L^{2}$ and $H^{1}$ for a projection-based scheme, showing that interpolation in these spaces is at least as good as interpolation in standard polynomial spaces on typical element shapes (e.g. triangles and quadrilaterals). The optimal order convergence of finite element approximations of a function having an unbounded gradient without employing small cells near the singularity, as well as the analysis provided for that specific example, indicates an even richer approximation theory that will be explored in subsequent work. In terms of practical computations, we described a boundary integral approach that was very recently developed with precisely these applications in mind. The numerical examples illustrated our convergence results on several families of meshes whose mesh cells are far from being simple perturbations of straight-edged polygons. We also numerically compared the approximation power of two types of harmonic bases on a mesh consisting of triangles with a single perturbed edge.

As highlighted at the end of Section 3, a better understanding of how geometric features of elements and choice of $p$ affect constants appearing in the interpolation analysis for $V_{p}(K)$ is needed. Additionally, an accounting of errors made in the approximation of quantities required for the formation of element stiffness matrices, and those subsequently arising from quadratures, should be taken into account as part of a more complete analysis of the method. Since the experiments in this work really only involved harmonic basis functions, all quadratures were performed on the boundaries of elements, but higher-order elements will require volumetric quadratures as well, and the development of efficient and robust volumetric quadratures in our context is another topic for further investigation. The method should also be tested on PDEs modeling more complex phenomena, and problems in which there are curved (and moving) interfaces between materials are of particular interest in this regard. Extending this approach to 3D problems in which general curved cells are permitted presents both theoretical and practical/computational challenges beyond the obvious analogues discussed above, and we aim to address them in future work. Among these is a definition of $V_{p}(K)$ that leads to a conforming space $V_{p}(\mathcal{T})$ without making the dimension of $V_{p}(K)$ much larger than is necessary to achieve optimal approximation properties. A second challenge is the efficient and accurate evaluation of quantities that are needed to form the local finite element linear systems; the approach outlined in Section 4 is inherently 2D, but there are integral equation approaches that may prove beneficial in our setting.