 # Trefftz Finite Elements on Curvilinear Polygons

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 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.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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. [10, 2, 23, 11, 12, 5, 38, 8, 13, 4, 6, 9]), 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. [25, 49, 48, 78, 66, 67, 79, 50, 80, 82, 81]), which have looked more toward the older Trefftz methods for motivation. A similar strategy has been followed in our previous work , 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. [41, 36, 65, 40, 57, 69] and the references in ), 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. [18, 19, 17, 32, 30, 31, 20]) and Weak Galerkin (WG) schemes (cf. [75, 76, 74, 62, 61, 60, 77]), as well as the recent adaptations of the discontinuous Petrov-Galerkin method (cf. ), in this regard.

The present work builds upon  for second-order, linear, elliptic boundary value problems posed on possibly curved domains : Find such that

 (1) ∫ΩA∇u⋅∇v+(b⋅∇u+cu)vdx=∫Ωfvdx+∫∂ΩNgvds for all v∈H ,

where is some appropriate subspace of incorporating homogeneous Dirichlet boundary conditions, and standard assumptions on the data 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. [68, 54, 14]), involve (local) mappings of standard mesh cells to fit curved boundaries, and these methods remain popular today. More recently, isogeometric analysis (cf. ), 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 [22, 20]. 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 three, all of which are very recent. In , 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,  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 , i.e. the polynomials of degree smaller or equal  are locally not contained in the local approximation space of order . This introduces additional difficulties in the study of approximation properties. Our earlier work  discusses the natural incorporation of Dirichlet boundary data on straight or curved portions of the boundary, allowing for elements having curved edges along portions of the boundary where Dirichlet data is prescribed. In that work, we described two potential approaches for handling defining elements having curved interior edges, and the present work realizes that goal both in practice, and by providing supporting interpolation theory. Our choice of discrete spaces is such that .

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 , 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 ) 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 in our local spaces provides the expected approximation power, and the presence of other (possibly singular) functions in is not detrimental. In , 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 , 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—some of which have only curved edges. 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 [24, 34], let be a connected subset of , with non-empty interior and compact closure, whose Lipschitz boundary, , is a simple closed contour consisting of a finite union of smooth arcs, see Figure 1. We will refer to as a mesh cell, the arcs as edges, and will implicitly assume that adjacent edges meet at an (interior) angle strictly between and , i.e. has no slits or cusps. We allow adjacent edges to meet at a straight angle. The vertices of are those points where two adjacent edges meet. Given an integer and a mesh cell , we define the space to be the polynomials of (total) degree at most on , with for , and the space to be continuous functions on whose trace on each edge is the trace of a function from (equivalently, from ) on , and we denote by this edge trace space. In , we refer to this definition of as its Type 2 version; the Type 1 version consists of functions on that are polynomials with respect to a natural parameter, such as arc length, in a parametrization of . 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 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 be a bounded domain with Lipschitz boundary. Given a partition of a , we define by

 (2) Vp(T)={v∈C(¯¯¯¯Ω):v|K∈Vp(K) for all K∈T} ,

where we define the space as follows,

 (3) v∈Vp(K) if and only if Δv∈Pp−2(K) in K and v|∂K∈Pp(∂K) .

The space clearly contains , but it typically contains other functions as well. A natural decomposition of is , where

 (4) v∈VKp(K) if and only if Δv∈Pp−2(K) in K and v=0 on ∂K , (5) v∈V∂Kp(K) if and only if Δv=0%inK and v|∂K∈Pp(∂K) .

The dimension of is

 (6) dimVp(K)=dimVKp(K)+dimV∂Kp(K)=dimPp−2(K)+dimPp(∂K)=(p2)+dimPp(∂K) .

The dimension of depends on the number and nature of the edges of . If is a straight edge, , but if is not a straight edge, the dimension of can be as high as , as it is when . The dimension of more generally is given in the following proposition, a proof of which may be found in [58, Theorem 7.1], for example.

###### Proposition 2.1.

Suppose that is an irreducible polynomial of degree , and that all points satisfy . It holds that . If does not lie on a real algebraic curve in the plane, then .

###### Example 2.2.

We consider the dimensions of the spaces , , for each of the three mesh cells in Figure 1. We have and . The continuity of functions in implies that

 (7) dimV∂Kp(K)=dimPp(∂K)=∑e⊂∂KdimPp(e)−# edges .

This formula holds for arbitrary . For , we have

 dimV∂K1(K) =(# straight edges)+2(# curved edges) , dimV∂K2(K) =2(# straight edges)+4(# curved conic edges)+5(# curved non-conic edges) .

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

Shuriken Half-Washer Two-Edge Circle
0+(0)+2(4)=8 0+(5)+2(1)=7 0+(0)+2(2)=4
1+2(0)+4(0)+5(4)=21 1+2(5)+4(1)+5(0)=15 1+2(0)+4(2)+5(0)=9

A basis for implicitly defines a basis for . In Section 4, we describe how we form the associated local finite element linear systems over , 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 be a point in , and be a multi-index. In , the authors provide an explicit formula for a polynomial satisfying . A basis of , , is given by , where

 (8) ΔψKα=0 in K,ψKα=−qα on ∂K .

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

 (9) ℓj(x)=1−(x−zj)⋅R(zj−1−zj+1)(√3/2)h2 ,

where we understand the subscripts modulo (i.e. and ), and

 R=(01−10),h=|z1−z0|,z2=z1+z02−R⋅√3(z1−z0)2 .

Any basis for yields a spanning set for 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, 16, 71]). For example, a hierarchical basis for is

 (10) {ℓ0,ℓ1,ℓ2,4ℓ1ℓ2,4ℓ0ℓ2,4ℓ0ℓ1,3√32ℓ1ℓ2(ℓ1−ℓ2),3√32ℓ0ℓ2(ℓ0−ℓ2),3√32ℓ0ℓ1(ℓ0−ℓ1),27ℓ0ℓ1ℓ2} ,

Here, we have chosen the scaling on each function so that its maximum value on the triangle, in magnitude, is . In any hierarchical basis for , the only functions that do not vanish at both and are and . A simple consequence of this fact is that

###### Proposition 2.3.

For any edge , a hierarchical basis of contains both and .

As stated in Proposition 2.1, if lies on an algebraic curve of order , we know the dimension of . However, it may be undesirable to make this determination in practice. Regardless, we need a practical method for paring down a spanning set for to a basis. Let , and suppose that is a hierarchical spanning set of , as described above. The functions are listed in increasing order of degree. The Gram matrix may be used to determine the remaining basis functions (in addition to ) for . We recall that (cf. [51, Theorem 7.2.10]). A basis for consisting of some subset of these functions may be determined using a rank-revealing Cholesky decomposition of (cf. [47, 46, 45]). 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 , and hence of .

###### Proposition 2.4.

Let , be the Gram matrix associated with an inner-product

and a list of vectors

. Let be a rank-revealing Cholesky decomposition, where is a permutation matrix, and , with the matrix having strictly positive entries. Then is a basis for , where is the permutation on defined by , and are the standard coordinate vectors.

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

###### Algorithm 2.5.

Let , , be the Gram matrix associated with an inner-product and a list of vectors . Upon termination of the following algorithm, , the indices of a basis of :

while end

Here, is the column of the current .

In practice, one replaces the condition with for some suitably small tolerance . Some speed-up of this basic algorithm may be achieved by exploiting the fact that previous reduction reduction steps, , 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 in our determination of a basis for .

###### Example 2.6.

For any edge , a hierarchical spanning set for is given by (10) where we have restricted the domains of these functions to . Let be parameterized by , , so is part of the hyperbola . We know in advance that , and will be part of our basis for , so we must select five of the remaining eight functions,

 {ℓ2,4ℓ1ℓ2,4ℓ0ℓ2,4ℓ0ℓ1,3√32ℓ1ℓ2(ℓ1−ℓ2),3√32ℓ0ℓ2(ℓ0−ℓ2),3√32ℓ0ℓ1(ℓ0−ℓ1),27ℓ0ℓ1ℓ2} ,

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 is given by

 {ℓ0,ℓ1,4ℓ0ℓ2,4ℓ0ℓ1,3√32ℓ1ℓ2(ℓ1−ℓ2),3√32ℓ0ℓ1(ℓ0−ℓ1),27ℓ0ℓ1ℓ2} .

These basis functions are plotted, as functions of the parameter , in Figure 2, together with the edge and associated triangle used to define the barycentric coordinates . Figure 2. At left, the edge e and associated triangle (dashed) for Example 2.6. At right, plots of a basis for P3(e) with respect to a parametrization, x=x(t), of e.

As a matter of interest, we note that, when the reduction algorithm was used, with as before, on the entire spanning set (10), a different basis was obtained,

 {ℓ0,ℓ1,4ℓ1ℓ2,4ℓ0ℓ2,4ℓ0ℓ1,3√32ℓ1ℓ2(ℓ1−ℓ2),3√32ℓ0ℓ2(ℓ0−ℓ2)} .
###### Remark 2.7.

We describe an alternative to the barycentric coordinates (9) associated with that acts more like a local cartesian coordinate system. Using the same notation , and , we define

 (11) ~ℓ0(x)=(z1−x)⋅(z1−z0)h2,~ℓ1(x)=(x−z0)⋅(z1−z0)h2,~ℓ2(x)=(z1−x)⋅R(z1−z0)h2 .

Straight-forward manipulations reveal that

 ℓ0=~ℓ0−~ℓ2√3,ℓ1=~ℓ1−~ℓ2√3,ℓ2=~ℓ2√3/2 ,

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

Having now properly defined , the discrete version of (1) is to find such that

 (12) ∫ΩA∇^u⋅∇v+(b⋅∇u+c^u)vdx=∫Ωfvdx+∫∂Ωgvds for all v∈Vp(T)∩H .

The intersection, , ensures that we respect any homogeneous Dirichlet boundary conditions inherent in . 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 is controlled by interpolation error , where is some appropriately defined interpolant of . 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

 (13) IK:W(K)={v∈C(¯¯¯¯¯K)∩H1(K):Δv∈L2(K)}→Vp(K) ,

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

 (14) I:W={v∈C(¯¯¯¯Ω)∩H1(Ω):Δv∈L2(Ω)}→Vp(T)

by .

Our definition of is motived by the decomposition . We begin with a related decomposition of as , where

 (15) {ΔvK=Δv in KvK=0 on ∂K,{Δv∂K=0 in Kv∂K=v on ∂K .

We define by an analogous decomposition , where and are given by

 (16)

In order to complete this definition, we must define and . We define by

 (17) ∫K(Δv−qK)ϕdx=0 for all ϕ∈Pp−2(K) .

We define by defining it on each edge of . For a non-trivial open subset , we use the inner-product

 (18) (ϕ,ψ)H1/2(Γ)=∫Γϕψds+∫Γ∫Γ(ϕ(x)−ϕ(y))(ψ(x)−ψ(y))|x−y|2ds(x)ds(y) ,

with as the associated norm. Below, we take to be either the entire boundary, , or a single edge, . Fix an edge of , having endpoints , and let . We define by the conditions

 (19) qe(z0)=v(z0),qe(z1)=v(z1),(v−qe,w)H1/2(e)=0 for all w∈Pp,0(e) .

Finally, is defined by .

In several places below, it will be convenient to use the following basic result. Suppose that , and in and on . Then , so . For example, we have

 (20) |v−IKv|2H1(K)=|vK−IKKv|2H1(K)+|v∂K−I∂KKv|2H1(K) ,

and we consider both contributions to the interpolation error in turn. In the proofs below, we use as a constant that may vary from one appearance to the next. Throughout, denotes the diameter of . We first consider .

###### Proposition 3.1.

Suppose that and for some . There is a scale-invariant constant for which

 |vK−IKKv|H1(K)
###### Proof.

It holds that

 |vK−IKKv|2H1(K) =−∫K(Δv−qK)(vK−IKKv)dx≤∥Δv−qK∥L2(K)∥vK−IKKv∥L2(K) .

Since , the Poincaré-Friedrichs Inequality, for , ensures that

 (21) |vK−IKKv|H1(K)

The estimate follows from this by applying the Bramble-Hilbert Lemma. The norm result follows from this by applying the Poincaré-Friedrichs Inequality again. ∎

###### Remark 3.2.

If is convex, can be replaced by in (21) (cf. ). Furthermore, for convex , the dependence on of the constant 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 depends on the shape of have been established [73, 33, 28, 27].

For our analysis of the following result will be useful.

###### Proposition 3.3.

If and in , there is a scale-invariant constant for which

 (22) ∥v∥L∞(K)≤cinf{hK|w|H2(K)+|w|H1(K)+h−1K∥w∥L2(K):w∈H2(K) and w=v on ∂K} .
###### Proof.

For , we have by a Sobolev imbedding result. A standard scaling argument then yields

 ∥w∥L∞(K)≤c(hK|w|H2(K)+|w|H1(K)+h−1K∥w∥L2(K)) ,

where is scale-invariant. Now (22) follows from the fact that harmonic functions on attain their extrema on , so, if have the same Dirichlet trace on , and is harmonic on , then . ∎

###### Remark 3.4.

Since we are working in , is continuously imbedded in for any . Therefore, Proposition 3.3 is readily generalized to such spaces, with the obvious bound

 (23) ∥v∥L∞(K)≤cinf{hsK|w|H1+s(K)+|w|H1(K)+h−1K∥w∥L2(K):w∈H2(K) and w=v on ∂K} .

Typical assumptions on the domain and the data for the problem guarantee that for some (cf.  [44, 43, 83]).

We now consider the term . Let be an edge of , with endpoints . We begin with a further decomposition of , namely , where and . This induces a natural decomposition of , , where satisfies at each vertex of , and vanishes at the vertices.

###### Proposition 3.5.

Suppose that for some . There is a scale-invariant constant for which

 |v∂K−I∂KKv|H1(K) ≤chpK|v|Hp+1(K),∥v∂K−I∂KKv∥L2(K)≤chp+1K|v|Hp+1(K) .
###### Proof.

We decompose as , where

It follows that .

We denote the set of vertices of by , and the set of edges of by . For , we define as follows: if is not adjacent to , then vanishes on , and if is adjacent to , then on , where for one of the endpoints of . It follows that , so

 |w1|H1(K)≤∥v∂K∥L∞(K)∑z∈V(K)|ℓz|H1(K)≤c∥v∂K∥L∞(K)≤c(hK|v|H2(K)+|v|H1(K)+h−1K∥v∥L2(K)) .

A similar argument shows that .

From (19), we see that , so , for each edge . Now,

 |w0|2H1(K)=∫∂K(∂w0/∂n)q∂K0ds≤c∥∂w0/∂n∥H−1/2(∂K)∥q∂K0∥H1/2(∂K)≤c|w0|H1(K)∥q∂K0∥H1/2(∂K) .

Here we have used the standard trace inequality , where is scale-invariant. Since is harmonic, . From this it follows that

 |w0|2H1(K)≤c∥q∂K0∥2H1/2(∂K)≤c∑e∈E(K)∥q∂K0∥2H1/2(e)≤c∑e∈E(K)∥v−qe,1∥2H1/2(e)≤c∥v−q∂K1∥2H1/2(∂K) .

The second inequality holds because vanishes at the vertices, see Remark 3.6. At this stage, .

Another standard trace inequality ensures that . Combining this with our estimates above, we obtain

 |w0|H1(K) ≤c(|v−w1|H1(K)+h−1K∥v−w1∥L2(K)) ≤c(hK|v|H2(K)+|v|H1(K)+h−1K∥v∥L2(K)) ,

and it follows that

 (24) |I∂KKv|H1(K)≤c(hK|v|H2(K)+|v|H1(K)+h−1K∥v∥L2(K)) .

A standard inverse inequality yields the obvious analogue in ,

 (25) ∥I∂KKv∥L2(K)≤c(