Unfitted Trefftz discontinuous Galerkin methods for elliptic boundary value problems

We propose a new geometrically unfitted finite element method based on discontinuous Trefftz ansatz spaces. Trefftz methods allow for a reduction in the number of degrees of freedom in discontinuous Galerkin methods, thereby, the costs for solving arising linear systems significantly. This work shows that they are also an excellent way to reduce the number of degrees of freedom in an unfitted setting. We present a unified analysis of a class of geometrically unfitted discontinuous Galerkin methods with different stabilisation mechanisms to deal with small cuts between the geometry and the mesh. We cover stability and derive a-priori error bounds, including errors arising from geometry approximation for the class of discretisations for a model Poisson problem in a unified manner. The analysis covers Trefftz and full polynomial ansatz spaces, alike. Numerical examples validate the theoretical findings and demonstrate the potential of the approach.


Introduction
In the last two decades, unfitted finite element methods became popular as an alternative to more traditional body-fitted methods to solve partial differential equations on complex geometries numerically.The idea is to separate the computational mesh from the geometry description to remove the burden of mesh generation, mesh adaptation and remeshing when dealing with complex (and possibly time-dependent) geometries.This is accomplished by embedding the domain of interest into an unfitted background mesh.In the context of finite element methods, unfitted discretisations go under different names such as CutFEM [9], extended FEM (X-FEM) [4,44], Finite Cell [46] and several more.In many cases, discontinuous Galerkin (DG) methods are attractive on unfitted meshes as they are for fitted meshes, e.g., for convection-dominated convection-diffusion problems, because of their flexibility in regards to the polynomial basis function as exemplified by Trefftz and ℎ-methods or computational aspects.In the spirit of CutFEM and X-FEM, there are several unfitted discretisations based on a discontinuous Galerkin formulation, cf.[3,10,19,21,30,42,43].
A challenge for unfitted methods arises from the fact that the boundary can cut through mesh elements arbitrarily.Trimming parts of the mesh outside the domain of interest may lead to ill-shaped elements.Two mechanisms have proven fruitful in making such methods robust with respect to ill-shaped elements: ghost penalties and element aggregation.
Element aggregation (AG) joins mesh elements to ensure that the support of basis functions does not degenerate for bad-cut configurations in the mesh.This has been done among others in [19,30] for DG methods and in [11] for unfitted Hybrid-High-Order methods.In [1] and the proceeding works of Verdugo and Badia, the idea has also been generalised to continuous finite elements.The method is also referred to as cell merging or cell agglomeration.
Ghost penalty (GP) introduces an additional stabilisation term, the ghost penalty stabilisation [6,9], that introduces a volumetric coupling, i.e., a coupling that involves all unknowns of two adjacent elements in the vicinity of shape-irregular cuts.No mesh elements or basis functions are changed for this approach.In order to reduce the number of coupled elements, the stabilisation term can be applied within patches only.These patches are built by the same machinery used by element aggregation, see [2].We refer to this method as patch-wise ghost penalty (GP); it also appears under the names of weak ghost penalty or weak aggregated elements.We also mention the work [17] for an approach where ghost penalties are avoided.The ill-cut problem is resolved here by removing basis functions whose supports have small intersections with the computational domain.
Discontinuous Galerkin methods come with increased degrees of freedom (dofs) when compared to their continuous counterpart.When solving linear systems for corresponding discretisations, the duplication of degrees of freedom affects the efficiency of numerical methods even more.There are two well-known remedies in the literature for the body-fitted case.
The first is to use Hybrid Discontinuous Galerkin (HDG) methods [15], where additional degrees of freedom are introduced on element interfaces.These additional unknowns allow (in many cases) the elimination of interior (volumetric) DG unknowns by a Schur complement strategy (known as static condensation).The remaining degrees of freedom in the global linear system are then significantly reduced, especially in the case of higherorder discretisations.A difficulty with the extension of HDG methods from the fitted to the unfitted case is the robustness with respect to shape-irregular cuts.Applying a ghost penalty stabilisation is not possible as the couplings between direct element neighbours introduced by the stabilisation terms contradict the decoupling exploited in the hybridisation.Cell merging strategies are possible but require the handling of polygonal meshes with corresponding facet functions as it is done in the unfitted Hybrid High Order (HHO) method [7,8].In [22], neither ghost penalties nor cell merging needs to be used by changing the background mesh to avoid ill-shaped cuts (in two space dimensions).
The second remedy to overcome the computational costs of DG methods is the class of Trefftz DG methods.In Trefftz methods, originating from [54], the ansatz space is constructed to lie in the kernel of the differential operator of the PDE at hand.Compared with a corresponding DG method, the same accuracy can be achieved for these methods at significantly reduced costs.Trefftz-DG methods for the Laplace problem are analysed in [26,39,40].Indeed, the complexity reduction is comparable to that of the HDG method, cf.[37].In contrast to the HDG mechanism, the mechanism to reduce the computational complexity does not interfere with either the ghost penalty stabilisation or the cell merging strategy.

Main contributions and outline of the paper
In this work, we will consider the combination of unfitted discontinuous Galerkin formulations with Trefftz DG finite element spaces on the example of the Poisson problem.Our objective is to develop the tools for analyzing and advancing the method beyond this particular model problem.To the best of our knowledge, this is the first occasion of a geometrically unfitted Trefftz DG method in the literature.
Our analysis covers the unfitted DG method with two choices for the basis functions: either the full polynomial space or the Trefftz space.These choices of basis functions are then combined with either the cell-merging, the ghost penalty, or the patch-wise ghost penalty stabilisation -cf.also Figure 3 below for a comparison of stabilisation strategies -to arrive at robust unfitted discretisations, which are then analysed.The analysis Figure 1.Overview on different discretisation settings treated in this manuscript.Twelve different settings arise from the possible choices in the basis discretisation (DG vs. Trefftz DG), stability measure for small cuts (ghost penalty in two versions vs. aggregation) and the assumption on the geometry handling (exact vs. approximated).covers a higher-order a-priori error analysis, which we first present for exact geometries and then extend to the case with geometry approximation errors.A summary of these methods presented here is given in Figure 1.
The unfitted DG method with ghost penalty has already been proposed and analysed in [21], except for the analysis of the geometrical errors.The element aggregation and patch-wise ghost penalty follow from the works [2] and [1] on the (more general) aggregated FEM.Our unified analysis is possible, as both the Trefftz and aggregated finite element spaces are subsets of the standard DG space.The patch-wise ghost penalty consists of a minor modification of the ghost-penalty term.
Novel contributions of the work are the -description of unfitted DG and unfitted Trefftz DG methods with three different stabilisation mechanisms, both with and without geometry approximation errors, -unified analysis leading to a priori error estimates for these unfitted (Trefftz) DG methods, including geometry errors, and -numerical experiments for these methods and the discussion of implementational aspects.
The remainder of this paper is structured as follows: In Section 1.2, we present the model problem under consideration in this paper.In Section 2, we recap the unfitted DG method under the assumption of exact geometry handling, as covered in previous literature, and introduce the Trefftz DG method under the same assumption on the geometry.In Section 3, we consider the different approaches to deal with the problem of small cuts, namely element aggregation in Section 3.1 and ghost-penalty stabilisation in Section 3.2.Section 4.1 presents the unified error analysis for the considered methods.The analysis extends the work on unfitted DG methods in [21].A crucial role for this extension is the special choice for the interpolant of smooth functions.Section 4.2 then covers the error analysis of the unfitted DG and TDG methods, including geometry approximation errors inherent in unfitted finite element methods.We then discuss several variants and implementational aspects of the covered methods in Section 5, including embedding the Trefftz and aggregated spaces into standard discontinuous finite element spaces.We present numerical examples of the methods in Section 6.Here we include examples with and without geometry approximation errors.Finally, in Appendix A, we give some proofs for completeness that consist only of minor adaptations of proofs available in the literature.Here, we used the inner product (, ) Ω := ∫︀ Ω d.We note that for the analysis below, we will homogenise the problem with respect to the volume forcing, meaning that we compute a particular solution of (1.1a) and then discretise an homogeneous problem, see (1.1 hom ) below.This is standard for Trefftz methods.We will also discuss the numerical realisation of this homogenisation in Section 5 and present a numerical example with  ̸ = 0.

Unfitted DG and Trefftz DG methods with exact geometry
In this section, we first recap the unfitted DG method assuming exact geometries, i.e. neglecting errors arising from inaccurate integration over cut elements.Based on this, we will then introduce the unfitted Trefftz DG method.

Preliminaries: Geometry, mesh and cut elements
We start by introducing some notation and assumptions: Let Ω be a background domain, sufficiently large such that Ω ⊆ Ω and let ̃︀  ℎ = { } be a division of Ω into non-overlapping shape-regular elements.The local mesh size of a mesh element  ∈ ̃︀  ℎ is defined as ℎ  = diam( ) := sup 1,2∈ ‖ 1 −  2 ‖ 2 .We allow for quite general polyhedral meshes including possibly curved elements, but require that the mesh conforms to the following assumption or its subsequent relaxed version.
(c) The element boundary  ′ is the union of a finite (yet, arbitrarily large) number of closed  1 surfaces.
Here and in what follows, we use the notation   if there exists a constant  > 0, independent of the mesh size and mesh-interface cut position, such that  ≤ .Similarly, we use if  ≥ , and  ≃  if both   and   holds.These assumptions are essentially based on Assumptions 4.1 and 4.3 from [13] to guarantee that the trace and inverse estimates cited from this work are valid here.For further details on these mesh assumptions, we refer to Section 4, Figure 2 and Figure 3 from [13].For a simpler but more restrictive mesh assumption for polytopic meshes under which appropriate trace estimates are available, we also refer to [12].
Let us note that (possibly curved) simplicial, hexahedral or quadrilateral meshes that are shape-regular in the usual sense fulfil Assumption 2.1.In the following, we want to enlarge the class of admissible meshes to those arising from merging a (uniformly bounded) number of (neighbouring) elements from meshes fulfilling Assumption 2. Remark 2.4.The first part of Assumption 2.1 (even in its relaxed version of Assumption 2.2) guarantees a shape regularity property sufficient for an inverse estimate and the interpolation via Taylor polynomials below, as needed for the analysis of the Trefftz method.The second assumption is essentially a bound on the curvature of the elements, although this assumption is a weak restriction compared with Assumption 2.1(a).By construction, starting from a mesh fulfilling Assumption 2.1 and applying a cell merging strategy (where always only a uniformly bounded number of neighbouring elements are merged) results in a mesh fulfilling Assumption 2.2.Obviously, further merging of a resulting mesh only fulfilling Assumption 2.2 will still yield a mesh fulfilling Assumption 2.2.However, these merging steps will decrease the shape regularity bound by a multiplicative constant (∼ −1  ).
Of specific interest will be the following parts of the background mesh, which we call the active mesh and the cut mesh, i.e. the parts of the background mesh that contribute to the covering of Ω and the part that is intersected by the domain boundary Γ, respectively: We collect the domain of the active mesh as   := Int ⋃︀  ∈ ℎ  .We further introduce sets of facets needed for the unfitted DG method.For the communication between all direct neighbour elements in a set of elements  = { } we collect and denote With abuse of notation we denote by ℎ the global mesh size ℎ = max  ∈ ℎ ℎ  when used as a scalar, as well as the piecewise constant field on Ω, ℎ : Ω → R with ℎ|  = ℎ  ,  ∈  ℎ or as the piecewise constant field on the skeleton, ℎ : Note that due to shape regularity ℎ  ≤ ℎ  ℎ  for  ⊂  .

Unfitted DG methods with exact geometry
Starting point and first method is the unfitted DG discretisation as in [21].The discrete function spaces are given as the discontinuous polynomials of order  on  ℎ : where P  ( ) is the space of polynomials up to degree  on  .
As usual with interior penalty DG methods, we penalise jumps across facets  ∈ ℱ ℎ .For this, we need the following average and jump operations: and   denotes some fixed facet normal to  .Next, in preparation of the discrete variational formulation, we introduce the bilinear form  ℎ as the corresponding linear form ℒ ℎ for the right-hand side as and a ghost penalty stabilisation form  ℎ (•, •), which we will specify in Section 3.2.
In terms of these definitions, we can then introduce the discrete problem as follows: Find  ℎ ∈ P  ( ℎ ) such that (DG)

Unfitted Trefftz DG methods with exact geometry
One major drawback of discontinuous Galerkin methods is the high computational cost related to the additional degrees of freedom due to the discontinuities.A popular remedy for this is to consider hybridised discontinuous Galerkin (HDG) methods.In HDG methods, additional unknowns are introduced at the facets of elements in order to remove direct couplings between neighbouring element unknowns.Thereby the size of the linear systems to be solved can be reduced by static condensation significantly, especially for higher order.This approach does not fit well with unfitted finite element methods in general, as stabilisation techniques, such as the ghost penalty method, rely on direct couplings between neighbouring elements and can not be hybridised efficiently.Only if ghost penalty methods can be circumvented, e.g. by the cell aggregation techniques discussed above, a hybridised approach can be applied as in the unfitted HHO methods, see, e.g., [11].
An alternative approach to reduce the size of the system to be solved is to consider Trefftz DG methods.Here the essential idea is to use a DG space consisting only of polynomials1 which fulfil the homogeneous PDE problem on every element.This leads to a complexity reduction of the number of degrees of freedom as well as the global couplings, which is similar to the HDG method, see the discussion in [37].The construction of the Trefftz DG method does not interfere with the usual coupling pattern of DG methods and can hence be combined with the ghost penalty method straightforwardly.
Let us now introduce a Trefftz version of the previous unfitted DG method simply by replacing the discrete function space to discontinuous polynomials of order , which satisfy the Trefftz condition of the Laplace problem: In order to make sense of this subspace of the previous DG space with respect to the inhomogeneous Poisson equation, we assume for now that an element-wise smooth particular solution   , with Δ  =  on each element (possibly discontinuous across element interfaces), is given.
Several approaches exist to homogenise the Laplace problem to apply Trefftz methods; see [37,47,55,56].The approaches [47,55,56] focus on collocation-based methods.The embedded Trefftz method, presented in [37], is DG based and gives an easy way to homogenise the system, which we review in Section 5.
In terms of these definitions, we can then introduce the discrete problem as follows: Find  T ∈ T  ( ℎ ) such that ∀ ℎ ∈ T  ( ℎ ) there holds (TDG) Remark 2.5.The bilinear form  ℎ (•, •) with Trefftz test and trial spaces was analysed in [26] for the fitted case.
Using integration by parts again on the first term in  ℎ (•, •), one obtains an 'ultra-weak formulation', where the volume term vanishes due to the properties of Trefftz test functions, then  ℎ (•, •) only requires integration over facets and no additional terms are needed.In the discrete setting, this provides an equivalent formulation and can bring considerable savings for the assembly of the linear system.For the Helmholtz problem, such an ultra-weak formulation has been studied in [14,27].

Stabilisation techniques
In this section, we consider different approaches to deal with stability in the presence of shape-irregular cut configurations.These approaches are element aggregation and different ghost penalty stabilisations.In the remainder of this work, we will cover all variants in a mostly unified manner.

Element aggregation
We can avoid the presence of shape irregular cut elements by cell-merging strategies.This has been done, among others, in [19,29,30] for DG methods and in [11] for unfitted Hybrid-High-Order methods.In [1] and the proceeding works of Verdugo and Badia, the idea has also been generalised to continuous finite elements.
We will use the clustering strategies which group certain elements  ∈  ℎ together in patches   ℎ = { }, as, e.g., presented in [1].These patches are directly used in a cell-merging strategy resulting in an aggregated element  = Int( ⋃︀  ∈  ℎ  ) per patch.However, they will also prove useful for the strategy based on ghost penalties; see also [2].
We denote  ℎ = {} as the set of aggregated elements obtained from the aggregation of two or more elements, and   ℎ = {  ℎ ,  ∈  ℎ } as the set of non-trivially aggregated elements.To every aggregated element we define the diameter ℎ  = diam().Finally,  ag ℎ denotes the active mesh after aggregation.Note that an aggregated element, i.e. an element in  ag ℎ , may originate only from a single element in the interior of the domain.In the next paragraph, we summarise the crucial properties of these aggregated elements.
The purpose of the non-trivial patches is to group every shape-irregular cut element together with at least one shape-regular root element, cf. Figure 2 for a sketch.Elements that are not adjacent to cut elements are not affected and form trivial patches.The number of elements in each patch is uniformly bounded so that the set of aggregated elements  ag ℎ itself fulfils Assumption 2.2.For the construction of the patches we formulate the following assumption.The fulfilment of this assumption is the (achieved) goal of the strategies in [1].Assumption 3.1.We assume that from any cut cell ℎ , and (c)  ≤  max , where  max is a fixed integer.Remark 3.2 (Shape-regular vs. shape-irregular elements).We note that the construction of the patches in the previous assumption is based on a distinction between cut and uncut elements rather than shape-regular and shape-irregular cut elements.This could be generalised further, e.g., based on | ∩ Ω|/| |,  ∈  ℎ .As it is standard in the CutFEM literature and avoids the introduction of additional notational burden, we stay with the simpler choice for ease of presentation.
Furthermore, the intersection of the aggregated element with Ω is shape regular in the sense that there is a constant  > 0 uniform in  ∈  ℎ so that | ∩ Ω| ≥ /||.
Proof.The first statement is stated and proven in Lemma 2.2 from [1].The second follows directly from quasi uniformity.
We note that the lemma implies that for all  ∈  ℎ we have ℎ  ≃ ℎ  for all  ∈   ℎ .After (proper) cell merging, an aggregated mesh is guaranteed to be cut-shape-regular in the following sense: Definition 3.4.An active mesh  ℎ that fulfils Assumption 2.2 is denoted as cut-shape-regular if there is a constant  > 0 (independent of ℎ) so that min  ∈ ℎ | ∩ Ω|/| | > .
Lemma 3.5.For a cut-shape-regular mesh  ℎ , it holds for  ℎ ∈ P  ( ℎ ) that Proof.This is an immediate consequence of Definition 3.4, and the constants in (3.1)only depend on the cut-shape-regularity constant.

Ghost penalty and patch-wise ghost penalty
In this section, we assume that we do not have a cut-shape-regular mesh and introduce a stabilisation (in two variants).In this case, we will still use the patches from the previous section to define regions where the stabilisation needs to act but not to change the mesh.
Using (2.1), we define the set of facets ℱ ℎ (  ℎ ).This is the set of all interior facets of a patch   ℎ ∈   ℎ and is only needed for the variant (GP).
For the ghost-penalty stabilisation, we require a set of facets that connects (possibly indirectly over several elements and facets) every shape-irregular cut element with a shape-regular root element.We denote this set as ℱ gp ℎ .As a minimal choice for stability, we take the set of all interior facets of all patches In this case, the connection from shape irregular to shape regular element is dealt with within each patch separately.A larger set, more often used in the literature, takes all facets between cut and uncut elements For the ghost penalty method, we set ℱ gp ℎ = ℱ gp⋆ ℎ ; for the patch-wise ghost penalty method, we set ℱ gp ℎ = ℱ gp,min ℎ .In the analysis, we will need ℱ gp ℎ ⊃ ℱ gp,min ℎ to prove stability and ℱ gp ℎ ⊂ ℱ gp⋆ ℎ to prove approximation properties.The global ghost penalty method is the most common approach in the literature, while the patchwise ghost penalty method has been considered in [2,31] and is sometimes referred to as the weak aggregation approach.
Different realisations of the ghost penalty stabilisation method are possible.These have essentially the same properties and decompose into facet contributions: where  > 0 is a corresponding stabilisation parameter.Two possible and popular choices for  ℎ, are: Here, Π  ℎ denotes the element-wise L 2 projection onto P  ( ℎ ),   = Int( 1 ∪  2 ) denotes the element aggregation to a facet  ∈ ℱ gp ℎ ,  =  1 ∩  2 .The facet (volumetric) patch jump    of a polynomial  ℎ ∈ P  ( ℎ ) is defined as where ℰ  denotes the canonical extension of a polynomial from  to Ω, i.e.
To keep the discussion simple, we only use the direct version  ℎ, =  2 ℎ, introduced in [48] in the following.Essentially, the ghost penalty stabilisation ensures control of finite element functions on cut elements by borrowing it from interior neighbours.For a more detailed discussion on ghost penalties and different realisations, we refer to [21,35].The main required property of the ghost penalty operator is that a stabilised version of Lemma 3.5 holds.Lemma 3.6.Under Assumption 3.1, we have for  ℎ ∈ P  ( ℎ ) that Proof.We refer to Lemma 5.2 from [35].

Summary of approaches
We finalise this section by collecting the three different methods of stabilisation under consideration: The ghost penalty method uses a stabilisation term, given in (3.2), acting on ℱ gp ℎ = ℱ gp⋆ ℎ .The element clusters are then used solely in the analysis for interpolation.In this case, the shape regularity of the aggregated elements (and the nestedness of finite element spaces) can be exploited to obtain optimal approximation results.

(GP)
For the patch-wise ghost penalty method, the same stabilisation term (3.2) is used over the minimal set of faces.The element clusters are used to reduce the regions where stabilisation is applied, as stability for bad-cut elements can be supported from within each cluster by choosing For the element aggregation method, mesh elements are merged, resulting in an active mesh  ℎ =  ag ℎ that is cut-shape-regular in the sense of Definition 3.4.Hence, no ghost penalty-like stabilisation is needed, and we set  ℎ (•, •) ≡ 0.

(AG)
We note that one could also think about a situation where one directly starts with a cut-shape-regular mesh  ℎ and then considers an unstabilised discretisation.In the analysis below, we will take this viewpoint for the case (ag), i.e., we simply assume to be given a cut-shape-regular mesh  ℎ for which no stabilisation is required.Given a particular mesh, we illustrate the choices available in Figure 3 and the resulting consequences for our mesh notation.

Error analysis
In this section, we analyse the unfitted DG and Trefftz DG methods under the assumption of exact geometry handling.The analysis for the unfitted DG method has already been treated in [21].We repeat the analysis with slight generalisations concerning mesh assumptions for inverse inequalities, the aggregated DG formulation and a patch-wise ghost penalty formulation.In particular, we generalise the approximation result by a special interpolation operator in Section 4.1.3.The particular interpolation operator allows us to conveniently extend the analysis further for the unfitted Trefftz DG method.

Error analysis of the unfitted DG and Trefftz DG methods with exact geometry
In this section, we will analyse the unfitted DG as well as the unfitted Trefftz DG method introduced above under the assumption that we are given a smooth particular solution   ∈ H ℓ (Ω), ℓ ≥ 2 of (1.1a).We homogenise the problem a priori, i.e. instead of , the solution to (1.1), we look for , the solution to the following variant: Hence, no additional homogenisation is needed in the numerical methods.The corresponding DG and Trefftz DG method (both with right-hand side ℒ ℎ (0,  hom ;  ℎ )) will be denoted as (DG hom ) and (TDG hom ).A numerical solution for the more generic case that no suitable   ∈ H ℓ (Ω) is known will be discussed in Section 5.
In the remainder of this work, we make the following assumptions.Firstly, we assume that the computational mesh  ℎ fulfils Assumption 2.2.Second, either the stabilisation bilinear form  ℎ (•, •) is chosen as described in Section 3.2 and Assumption 3.1 is valid, or  ℎ is already cut-shape-regular so that  ℎ (•, •) = 0.In the latter case, we set  ag ℎ =  ℎ , i.e. all patches considered below are trivial patches only consisting of one element in  ℎ .

Trace inequalities and coercivity
We gather several trace estimates which hold due to Assumption 2.2.Lemma 4.1 (Continuous trace inequality).Let  ∈  ℎ .For all  ∈ H 1 ( ), we have Proof.For meshes that fulfil Assumption 2.1, the proof can be found in Lemma 4.7 from [13].Let  ∈  ℎ with subelement { ′ } as in Assumption 2.2.We can then divide the boundary of  into the subelement contribution and apply the trace inequality for the subelement  ′ , .Lemma 4.2 (Discrete trace inequality).Let  ∈  ℎ .For all  ℎ ∈ P  ( ), we have with a constant  > 0, independent of the local mesh size ℎ  and the order .
Proof.The claim follows from Lemma 4.1 and an inverse inequality applied on the subelements for all  ℎ ∈ P  ( ), and a constant  > 0 independent of ℎ and .See Corollary 4.24 from [13] for details on the inverse inequality applied on  ′ under our mesh assumptions.
Proof.Similar proofs to similar statements are given e.g. in [7,16,23].For the framework at hand, the proof of Lemma 5.2 from [41] can be used with only slight adaptations in the setting.For completeness, we included the proof in Lemma A.1.
In conjunction with the problem (DG), we introduce the following (semi-)norms for the analysis below if  > 0 is sufficiently large, independent of ℎ and .Finally, the problem (DG) admits a unique solution.
Proof.The ghost-penalty stabilised case is covered by Proposition 2.6 from [21] for shape-regular (but not necessarily cut-shape-regular) simplicial meshes.We briefly repeat the arguments for this setting.The continuity estimates follow immediately from the Cauchy-Schwarz inequality.For the coercivity, we have with the Cauchy-Schwarz and weighted Young's inequality We then observe with   =   • ∇|   ℎ and the discrete trace inequality in Lemma 4.2 that with a constant  1 > 0 independent of the local mesh size ℎ  and order .With a cut version of the discrete trace inequality, obtained by combining Lemma 4.3 and (4.4), we also have with a constant  2 > 0 independent of the mesh size (field) ℎ and order .In the cut-shape-regular case, we can use Lemma 3.5 to bound the right-hand side of (4.9) and (4.10) by a norm on Ω, and in the ghost penalty stabilised case, we can use Lemma 3.6 to the same effect.Consequently, we have where we recall that in the cut shape regular case the ghost penalty terms are empty, i.e.  ℎ (•, Combining this with the above estimate, choosing  ≤ 1 2 2 and subsequently  ≥ 4 2 , we have In the following we will always assume that  is sufficiently large so that Lemma 4.4 holds.Proof.Coercivity and continuity as stated in Lemma 4.4 are directly inherited on the subspace T  ( ℎ ) ⊂ P  ( ℎ ), which follows per definition of the DG and Trefftz spaces in (2.2) and (2.3).

Approximation
For the approximation results, we will use averaged Taylor polynomials on adapted domains as the interpolation operator, cf.Section 4.1 from [5].We repeat its crucial properties and formulate them in a suitable setting for the following analysis.Lemma 4.7.Let  be a domain with Lipschitz boundary and assume that there are two balls   and   with   ⊂  ⊂   and associated diameters ℎ  and ℎ  so that ℎ  /ℎ  =   1.We define the interpolation operator    : L 1 (  ) →   (  ) as the operator realizing the averaged Taylor polynomial of degree  by averaging over   , see [5].
(1) For  ∈ W +1, (  ) and  ≥ 1 it holds that with a constant only depending on   ,  and .The second property of commuting interpolation and differentiation is crucial in the context of the Trefftz DG subspace.Note that in (4.12b)  and    are effectively only evaluated on the small ball   .It specifically implies that for  ∈ H  ( ) with Δ = 0 (on   ) and  ≥ 2 there holds Δ    =  −2  Δ = 0 on each element (for  < 2 there trivially holds Δ    = 0).In other words, the averaged Taylor polynomial of a harmonic function is also harmonic.
As a consequence of Lemma 2.3 and the fact that the aggregated elements in  ag ℎ again fulfil Assumption 2.2 we make the following observation in preparation of the subsequent lemma.
where the constant depends on the maximum number of subelements  ′ in each element,   , cf.Assumption 2.2, and the maximum number of elements in a patch  max , cf., Lemma 3.3.
Proof.We first recall that in the case where  ℎ is cut-shape-regular we have  ag ℎ =  ℎ .The first inequality is obvious due to T  ( ag ℎ ) ⊂ T  ( ℎ ).For the second, we will find an interpolator  ℎ ∈ T  ( ag ℎ ) with the corresponding bound.Before we discuss the construction in more detail we want to stress that the interpolator will be constructed w.r.t. the finite element space on  ag ℎ which is smaller than the one used in the discretisation only in the case where  ℎ is not cut-shape-regular (otherwise  ag ℎ =  ℎ holds).First, we bound all different norm contributions by H  (T)-semi-norms on elements  in the active mesh  ℎ .To this end, we recall the trace inequalities (4.2) and (4.5).Applying these estimates to all norm contributions of |||•|||  ℎ to  ∈ H  ( ℎ ) yields For the construction of the interpolant  ℎ =   ℎ , we apply the averaged Taylor polynomial from Lemma 4.7 patch-wise for all  ∈  ag ℎ so that   ℎ |  =    .Making use of Corollary 4.8 the averaged Taylor polynomial   ℎ |  is based only on |  with a ball   ⊂  ∩ Ω.By construction  ℎ depends only on  in Ω and hence with (4.12b) is harmonic, i.e.  ℎ ∈ T  ( ag ℎ ) ⊂ T  ( ℎ ).We cannot directly plug in the approximation error bound for , as  is not defined on   .We hence make use of a linear extension operator ℰ : H  (Ω) → H  (  ), for  ≥ 0, for which it holds ℰ| Ω =  and ‖ℰ‖ H  (  ) ‖‖ H  (Ω) , see for example Section VI.3 from [52].For the interpolant of the extension ℰ on the entire aggregated element , we choose the same averaged Taylor polynomial, so that   ℎ ℰ =   ℎ  ∈ T  ( ag ℎ ).We note that we constructed the averaged Taylor polynomial especially so that   ℎ only depend on values in Ω.We set  ℎ =   ℎ (ℰ).Applying (4.12a) gives for  ∈ H +1 (Ω) with Δ = 0 and each  ∈  ag ℎ the estimate for ,  ≤ .Hence, we have (with a finite overlap argument for the domains   ,  ∈  ag ℎ ) )︁ .
For the ghost-penalty semi-norm, let us consider two aligned elements  1 ,  2 and an aggregated element from a single facet patch   = Int( 1 ∪  2 ) and    = { 1 ,  2 }.We denote with Π   the L 2 projection onto the polynomial space over the domain   , and by Π   the element-wise projection onto the broken polynomial space over the elements    .Let us denote  ℎ, = ℰ  Π   ℎ ∈ P  (  ),  = 1, 2. Then for any  ∈ P  (  ) where the second last step follows using the shape regularity of the mesh, by which we can bound the norm of the discrete function on  1 by the norm on  2 and vice versa as the arguments are polynomials (not only element-wise) on the aggregated element.Let  = Π   (ℰ) and recall The operators   ℎ , Π   have the usual optimal approximation bounds.We hence obtain We note that in the case of the patch-wise ghost-penalty operator, i.e., ℱ gp ℎ = ℱ gp,min Proof.Follows from T  ( ag ℎ ) ⊂ T  ( ℎ ) ⊂ P  ( ℎ ) and Lemma 4.9.

.16)
Proof.Due to the elliptic regularity assumption, we have that for the auxiliary problem , it follows that ∇ •   = 0 and  = 0 on all  ∈ ℱ ℎ ∩ Ω. Due to the symmetry and consistency of the symmetric interior penalty form  ℎ (•, •), we have that Furthermore, we have the perturbed Galerkin-orthogonality Now let  1 ℎ be the interpolation operator by average Taylor polynomials onto P 1 ( ag ℎ ) = T 1 ( ag ℎ ) ⊂  ℎ .We note that all piecewise linear functions are naturally harmonic.Then by (4.17) and Lemma 4.4, we have where the penultimate estimate follows from Lemma 4.7 and the consistency of the ghost-penalty semi-norm shown in Lemma 4.9.The claim then follows from Corollary 4.11.

Unfitted DG and Trefftz DG methods with geometry approximation
Unfitted finite element discretisations pose the additional challenge of geometry approximation, i.e., the demand for an accurate representation of a geometry not described through the computational mesh and the need for robust and accurate numerical integration over cut elements.Several techniques to achieve these are known in the literature; see for example, [20,28,32,43,45,49].In our numerical examples below, we will consider an approach based on piecewise linear reference configuration and a (small) local mesh deformation [32].In this section, we introduce the methods with respect to a discrete approximated geometry and carry out an error analysis based on a Strang-type lemma.Similar techniques have been applied in [33,36] and the works by Deckelnick, Elliott, Ranner, e.g.[16,18].For ease of presentation we restrict to the case of quasi-uniform meshes, i.e. ℎ  ≃ ℎ for all  ∈  ℎ .

Geometry approximation
In the remainder we assume that a geometry approximation Ω ℎ of higher order is given, i.e., where  is the geometry order of approximation and we assume that integrals on Ω ℎ can be computed accurately.We further assume that there is a mapping Φ ℎ : Ω ℎ → Ω that allows to map the approximated domain onto the exact domain.This mapping is assumed to be a piecewise smooth bijection, and fulfils Φ ℎ (Γ ℎ ) = Γ and In this setting, we introduce adjusted versions of the previous bi-and linear forms and discrete norms of the DG discretisations.To this end, we effectively replace Ω by its approximation Ω ℎ yielding slightly modified discrete regions Let us stress that all previous dependencies on Ω, e.g., the selection of active elements  ℎ , are now replaced by dependencies on the discrete domain Ω ℎ .As in Section 4.1, we restrict ourselves to the homogeneous case  = 0.For the boundary data given by  hom ∈ H 1,∞ (Γ), we assume that there exists a sufficiently smooth extension into a domain Γ  ⊃ Γ∪Γ ℎ .We refer to the extension by  hom , abusing the notation.The plain-DG discretisation with geometry errors then reads as: Find  ℎ ∈ P  ( ℎ ), such that (geoDG hom ) Similarly, the Trefftz discretisation reads: Find  T ∈ T  ( ℎ ), such that

A Priori error bounds
For the error analysis, we introduce the auxiliary bilinear form with respect to the exact geometry.For We split the bilinear form  ℎ into a part containing the symmetry and boundary control terms and another containing the remainder, i.e. the bulk, consistency and DG terms We note that it is not necessary to split the linear form  ℎ due to  = 0. We have the following Strang type lemma for the DG method: Lemma 4.13.Let  ∈ H 2 (Ω) be the solution to (1.1) with data  hom ∈ H 1 2 + (Ω ∪ Ω ℎ ),  > 0. Furthermore, let either  ℎ = P  ( ℎ ) and  ℎ ∈  ℎ be the solution to (geoDG hom ), or let  ℎ = T  ( ℎ ) and  ℎ ∈  ℎ be the solution to (geoTDG hom ).Then Proof.The proof follows the lines of [33,36].Let us denote ũ =  ∘  ℎ .With the triangle inequality we have for arbitrary With  ℎ :=  ℎ −  ℎ and Lemma 4.4, we see Dividing by ||| ℎ |||  ℎ and observing that both Due to integration by parts, we have (, For the final term, we observe Proposition 4.14 (Unfitted DG error estimate.).
is extended sufficiently smooth into the discrete domain.For the solution  ℎ ∈ P  to (geoDG hom ) there holds , and () := Proof.We need to estimate the terms on the right-hand side of the Strang estimate in (4.19).With a suitably continuous extension operator ℰ : H +1 (Ω) → H +1 (  ) (for  ≥ 3) as above or ℰ : W 3,∞ (Ω) → W 3,∞ (  ) (for  ≥ 2), we set   := ℰ.We then have For the first term, we have The proof of this bound follows along the lines of Lemma 12 from [33].The additional interior-penalty facet contributions not treated in Lemma 12 from [33] have the same structure as the L 2 (Γ ℎ )-expressions treated in Lemma 12 from [33].We note that the cases  ≤ 2 and  ≥ 3 are distinguished as the proof requires  ∈ W 2,∞ (Ω) which is implied by  ℓ (Ω)-regularity only for ℓ ≥ 4. For the approximation part in (4.22), we use the unfitted projection based on averaged Taylor/polynomials as above.With  ℎ =   ℎ ℰ =   ℎ  per construction, we have as above For the ghost-penalty semi-norm, we further have from the consistency error estimate (4.13) from above, that For consistency error contributions in the Strang estimate (4.19), we have for  hom ∈ H 1,∞ (Γ) and assuming the data extension is bounded on for  ℎ ∈  ℎ .We provide the proofs in Appendix A.2. Combining these estimates proves the claim.
is extended sufficiently smooth into the discrete domain.For the solution  T ∈ T  to (geoTDG hom ) there holds Proof.The proof is an immediate consequence of the proof of Proposition 4.14, by the observation that due to Δ = 0 and the choice of our interpolation operator, we have in (4.23) that  ℎ =   ℎ  hom ∈ T  .

Implementational aspects
This section discusses some implementational aspects of the Trefftz and aggregated DG methods.In particular, we discuss how these approaches can be implemented in a general DG code without having to implement the basis functions for the Trefftz space or the DG space on general patches.

Embedded Trefftz method
In the previous section, we did not discuss the construction of T  .One way is to set up the basis of harmonic polynomials; see, e.g., [25].A more flexible choice is to construct T  through an embedding in the corresponding DG space P  .This has recently been introduced as the embedded Trefftz DG method [37].One of the advantages of using the embedded Trefftz DG method in this setting is that it allows the easy construction of a particular solution so that we can also deal with the inhomogeneous problem.We sketch the approach in this section but note that there are no specific adjustments for the unfitted setting that needs to be taken compared to the embedding procedure introduced in [37].
As the Trefftz space is the kernel of Δ in P  ( ℎ ), we are equivalently looking for a basis of ker(W) to characterise the Trefftz space.As W is block diagonal, with the blocks corresponding to the elements of the active mesh, we construct W element wise.
For  ∈  ℎ , let W  ∈ R   ×  be the block in W ∈ R  × corresponding to the element  with   = dim(P  ( )).The dimension of the kernel of W  is   = dim(ker(Δ)) =   −   with   = dim(range(Δ)) (e.g. in 2D   = ( + 2)( + 1)/2 − ( − 1)/2 = 2 + 1).We can compute the kernel of W  and collect the set of orthogonal basis vectors in the matrix T  ∈ R   ×  .This can be done numerically, for example using a QR-decomposition or singular value decomposition (SVD), see also Figure 4.These blocks are then collected into a block matrix T ∈ R  × , with which we have the characterisation ker(W) = T • R  .
With the help of T, we can now naturally define the Trefftz Galerkin isomorphism as  T : R  → T  ( ℎ ),  T (x) = (Tx).Having assembled B, l and T for the standard DG method, the system corresponding to (TDG) reads as follows: Find u T (=  −1 T ( T )), such that Note that the Trefftz system with matrices B T = T  BT and l T = T  l can also be assembled in an elementby-element fashion avoiding the setup of the (larger) matrices and vectors B, l.
The embedded Trefftz approach allows the construction of a particular solution   ∈ P  ( ℎ ) in a generic way.To this end, we compute element-wise (w  )  = (, (e  )) = (, Δ  ) and define (u  )  = W †  w  .Here W †  denotes the pseudoinverse of the matrix W  , which can be obtained using the QR decomposition or SVD of the matrix W  , which may have already been computed when numerically computing the kernel of W  .Then a particular solution is given by   =  T (u  ), and the Trefftz solution to the corresponding homogenised problem corresponds to the solution u T to Let us summarise that the key idea of the embedded Trefftz approach is to exploit the characterisation of the Trefftz DG space as the kernel of differential operator on a standard DG space on a linear level.A similar approach can also be applied to implement the aggregated finite element method as we discuss next.

Embedded aggregated FEM
In the aggregated finite element method, degrees of freedom on problematic elements, i.e. cut elements or ill-shaped cut elements, are marked as dependent degrees of freedom which are determined from degrees of freedom on uncut elements using an extrapolation.In [1], this interpolation is done geometrically based on a nodal representation of finite element degrees of freedom.Here, we want to discuss a different approach in the virtue of the embedded Trefftz DG approach.Since P  ( ag ℎ ) ⊂ P  ( ℎ ), we would like to characterise P  ( ag ℎ ) as the kernel of an operator in P  ( ℎ ).In fact, we can identify the patch-wise jump operator, i.e. the facet-patch-local version of the ghost-penalty operator, as a corresponding suitable operator.
As before, the matrix W is block diagonal, with each non-trivial block corresponding to one patch and the zero blocks corresponding to all trivial patches.As before, we can therefore compute the blocks W  independently, compute the kernel of dimension   = dim(P  ()) numerically and collect the orthogonal basis vectors of ker(W  ) in small matrices T  .Let  = dim(P  ( ag ℎ )) = dim(P  ( ℎ ∖   ℎ )) +   | ℎ |.We then build a block diagonal matrix T ∈ R  × from these blocks, together with an identity block corresponding to all trivial patches.
With the help of T we can now naturally define the aggregated Galerkin isomorphism as  ag : R  → T  ( ℎ ),  ag (x) = (Tx).
With the notation as for the embedded Trefftz method, the embedded aggregate finite element method corresponds to finding u  such that T  BTu  = T  l.
Note that here the ghost penalty is part of B, but is used to determine T.
The generic nature of this formulation is an advantage over the extrapolation based approach in [1].The extrapolation in [1] is based on a geometrical extrapolation exploiting a nodal representation of degrees of freedom.In contrast, our generic formulation allows for an implementation of the aggregation, which is independent of the choice of basis functions and can easily be generalised to curved elements.

Embedded aggregated Trefftz DG methods
Both previously discussed approaches to embed a special DG finite element space into the standard finite element space can be combined.We simply combine the previous components to obtain an embedding for the aggregated Trefftz DG method into the standard finite element space.Defining (W)  =  ℎ (  ,   ) + ⟨Δ  , Δ  ⟩ 0,ℎ , and determining the kernel of the corresponding operator in P  ( ℎ ) yields the space of patch-wise harmonic polynomials, i.e. the aggregated Trefftz space.The procedure can be executed -as in the previous sectionpatch by patch and results in a linear system of the form with patch-wise block diagonal T. Note that here the ghost penalty is again not part of B, but is used to determine T. The results can be seen in Figure 6.We again observe optimal convergence in the L 2 -norm for all considered polynomial orders after some initial pre-asymptotic behaviour.In contrast to the previous examples, we see that the error constant favours the Trefftz methods rather than the DG method.Concerning performance, we again see that the TDG method is significantly faster than the DG method.The exact compute times are again presented in Appendix B. As we can now apply the zero Dirichlet boundary condition on the discrete interface, the geometry approximation error is present in this example.To observe higher-order convergence, we, therefore, apply the strategy of isoparametric unfitted finite elements [32,33,36].Here, a cheaply computable and small (magnitude ℎ 2 ) mesh deformation Θ ℎ is applied on the background mesh in a way such that a piecewise linear geometry approximation Ω lin is mapped close to the exact geometry (in a higher order way), Ω ℎ = Θ ℎ (Ω lin ) ≈ Ω. Thereby, the problem of numerical integration can be reformulated as the much simpler problem of numerical integration on a piecewise linear geometry.This enables robust higher-order convergence of the quadrature problem.At a first glance, changing the mesh may seem to contradict the unfitted finite element paradigm.However, we recall that the difficulty that is to be circumvented in unfitted finite element methods is the initial meshing problem (or the The exact boundary condition is applied on the approximated boundary and no higher-order geometry approximation is applied.Right: Compute times for specific steps of the methods with ℎ = 2 −3 ,  = 4 and using 1, 2, 4 and 12 parallel threads respectively (Timings computed on an Intel Xenon E5-2687W v4 Processor @3.0 GHz).remeshing for moving domain problems) which is a non-local problem determining the mesh topology.In the isoparametric unfitted finite element approach, however, the necessary mesh alteration, the mesh deformation, is only local and does not change the mesh topology, i.e. the advantages of unfitted finite elements in terms of the meshing problem are not touched.We refer to [32] for further details.In contrast to previous works with the isoparametric unfitted finite element methods, we exploit the flexibility of discontinuous Galerkin methods 5and keep the coordinate system on the curved elements in world coordinates, see Figure 7.

Example 3: inhomogeneous problem and higher-order geometry approximation
We have seen in the previous two examples that the Trefftz and embedded Trefftz methods give indistinguishable results.While the standard Trefftz method does not yield optimal approximation rates when applied to the inhomogeneous problem directly, the embedded Trefftz method immediately provides a particular solution to homogenise the problem with, as described in Section 5.1.Therefore we only show results for the embedded Trefftz method here.We again only use the global ghost-penalty choice.
The results can be seen in Figure 8.Here we observe some pre-asymptotic behaviour on the three coarsest meshes.We attribute this to the fact that the width of the ring-shaped domain is of the same order of magnitude as the coarsest mesh size.As a result, the geometry and mesh deformation cannot be approximated properly on the coarser meshes.In particular, this becomes worse with increasing order of the deformation.However, we see optimal-order convergence for all orders and both methods once the mesh is sufficiently fine.We also see the slightly superior error constant for the DG method in the higher-order cases.
Figure 7. Left: straight triangle and straight cut configuration, Right: curved element with higher order geometry approximation and coordinate system with respect to world coordinates Figure 8. Example 3: L 2 -error convergence for the unfitted embedded Trefftz and DG methods over a series of meshes with different polynomial orders.The boundary is approximated to higher-order using a deformed mesh.

Example 4: species dissolution from a circle subject to a convection-diffusion equation
This manuscript focuses on the equation as a model PDE to introduce concepts for stable, reliable and arbitrarily high-order accurate unfitted (Trefftz) DG methods and their analysis.In this last example, we demonstrate -without error analysis -that the same concepts can be applied to more general applications.Here, we consider a convection-dominated convection-diffusion equation on a square background domain ̃︀ Ω = (−1, 1) 2 with a circular obstacle Ω cir =   ((0, 0)) of radius  = 1/4, i.e., Ω = ̃︀ Ω ∖ Ω cir .A divergence-free flow field that is tangential on the obstacle and essentially describes the advective transport from the left boundary ( = −1) to the right boundary ( = 1) is given by  = (1 +  2 ( 2 −  2 )/ 4 , −2 2 / 4 )  and we set  w = 2 ≃ ‖w‖ ∞ as a constant for the velocity magnitude.The convection-diffusion equation reads: Find  : Ω → R such that −Δ +  • ∇ = 0 in Ω, (6.1a) = 1 on Ω cir .(6.1d) We define Γ := Γ  ∪ Ω cir as the part of the boundary Ω where Dirichlet boundary conditions are prescribed (so that the notation for the Nitsche formulation in ( ℎ ) fits this setting).We consider  = 10 −3 , i.e. a strongly convection-dominated configuration with Pclet number Pe = w(2)  = 1000.This problem is challenging for standard  1 -conforming Galerkin discretisations due to a strong (parabolic) boundary layer that forms near the obstacle, and proper convection stabilisation becomes necessary if the boundary layer is not resolved, cf., e.g., [34].In the context of unfitted DG methods considered here, we can use a comparably simple upwind discretisation for the convection discretisation.
The discrete DG problem reads: Find  ℎ ∈ P  ( ℎ ) such that where  ℎ (, ) is the upwind DG bilinear form: with the upwind numerical flux ̂︂    = w • n  lim →0 + (x − w) on interior facets and outflow boundary facets and ̂︂ = 0 on all inflow boundaries, i.e. boundaries with w • n Ω < 0. For simplicity, we only consider the global ghost penalty stabilisation and adjust the ghost penalty scaling to consider the contribution from the convection and choose  =  0 ( + ℎ w ) with  0 = 0.001.
To define a proper Trefftz method, we can no longer use harmonic polynomials but have to use a more generic construction of Trefftz basis functions that are at least approximately in the kernel of ℒ = −Δ + w • ∇.To this end, we apply the idea of weak Trefftz methods, implemented through an embedding into the DG space as introduced in detail in [37] and define T  ( ℎ ) := { ℎ ∈ P  ( ℎ ) | Π  (−Δ ℎ + w • ∇ ℎ ) = 0} where Π  is the  2 projection into P −2 ( ℎ ).
Note that this is indeed a generalisation of the previously used Trefftz space as we recover the space of harmonic polynomials for w = 0 and  = 1.Further, the thusly defined Trefftz DG space has the same dimension as the space of harmonic polynomials, i.e. the same computational advantages over DG methods as for the Laplace problem.The discrete Trefftz DG problem then reads: Find  ℎ ∈ T  ( ℎ ) such that  ℎ ( ℎ ,  ℎ ) +  ℎ ( ℎ ,  ℎ ) +  ℎ ( ℎ ,  ℎ ) =  ℎ (, ;  ℎ ) ∀ ℎ ∈ T  ( ℎ ). (CD-TDG) Figure 9 shows the results for the unfitted DG and the unfitted Trefftz DG method of three decreasingly finer meshes (with  = 4) where the meshes are curvilinear to obtain higher order geometry approximation (as in Section 6.3).We observe that strong oscillations are avoided by the upwind (and the ghost penalty) stabilisation for both discretisations.Both schemes can capture the characteristics of the problem even in these under-resolved situations with mesh Pclet numbers Pe ℎ, = w  ℎ  ∈ {50, 25, 12.5}.Corresponding to the characterisation of the weak Trefftz DG basis functions, we observe a slightly different pattern in the small oscillations around 0 for the Trefftz DG than for the DG method.While the high oscillations for the DG method do not show an explicit direction, the small oscillations of the Trefftz DG solution align with the flow field.However, the magnitude of the oscillations for both methods are similar and decay under refinement, as one would expect.
We have presented these results -without going further into the details -to illustrate that the methodology presented in this manuscript has considerable potential for a broader class of unfitted PDE discretisations beyond the Laplace equation.Proof.We apply a change of coordinates, the corresponding mapping is denoted by Φ  .We recall that we assumed Γ ∈ C 2 .Let   be a point in  Γ =  ∩ Γ and   be the tangential plane to   .First, we apply a translation with −  and scaling with ℎ  to shrink  to a domain of size  (1).Second, we apply a rotation, denoted by the orthogonal matrix   , so that the tangential plane after transformation, i.e.Φ  (  ) aligns with the  1 - 2 plane (or  1 axis in 2D).We then have  = Φ  () = ℎ −1    •(−  ) and denote T = {Φ  () |  ∈  }.For sufficiently fine mesh sizes the unit outer normal of T Γ = Φ  ( Γ ) is close to   , the th unit vector, especially there holds (,   ) ≥ 1 −   ℎ ≥  ∈ (0, 1) where   only depends on the maximum curvature of  Γ which is uniformly bounded.Furthermore we have that either  ∩ Ω or  ∖ Ω is shape regular (and hence allows for the application of Lem.4.1).We denote this part as  * and T * = Φ  ( * ).Hence, we have for each v ∈ H 1 ( T ) that there holds  A.2. Proof of estimates (4.24) Proof of (4.24a).See also [36,Appendix A.4].For ease of notation let ̃︀  =  ∘  ℎ and ̃︀  ℎ =  ℎ ∘  −1 ℎ .We deal with the individual contribution in (4.24a) individually and begin with the volume terms.First we see using the chain-rule )︁ • ∇ ℎ d The final bound follows from (4.18) and the norm equivalence ‖̃︀ ‖ Ω ℎ ∼ = ‖‖ Ω which can be proven along the same lines of argument.
For the symmetric interior penalty consistency term, let us consider single facet  ∈ ℱ ℎ .For simplicity, we identify  with  ∩ Ω.Furthermore, we denote F = Φ ℎ ( ),   as the unit normal vector on  and n as the unit normal vector on F .From [36], we then have

Figure 2 .
Figure 2. Sketches of the different sets of elements and facets for the unfitted methods.Left: Global ghost-penalty stabilisation elements and facets.Right: Aggregated elements (elements gathered in a patch have the same colour) and corresponding inner-patch facets corresponding to patch-wise ghost-penalty stabilisation.

Figure 3 .
Figure 3. Stabilisation choices dependent on the input mesh and resulting mesh and facet sets.

Figure 4 .
Figure 4. Sketch of extraction of the kernel of the local matrix W  based on a local SV decomposition.

Figure 5 .
Figure 5. Example 1: Left: L 2 -error convergence for the unfitted Trefftz, embedded Trefftz and DG methods over a series of unstructured meshes with different polynomial orders.The exact boundary condition is applied on the approximated boundary and no higher-order geometry approximation is applied.Right: Compute times for specific steps of the methods with ℎ = 0.1,  = 5 and using 1, 2, 4 and 12 parallel threads respectively (Timings computed on an Intel Xenon E5-2687W v4 Processor @3.0 GHz).

Figure 6 .
Figure 6.Example 2. Left: L 2 -error convergence for the unfitted Trefftz, embedded Trefftz and DG methods over a series of structured tetrahedral meshes with different polynomial orders.The exact boundary condition is applied on the approximated boundary and no higher-order geometry approximation is applied.Right: Compute times for specific steps of the methods with ℎ = 2 −3 ,  = 4 and using 1, 2, 4 and 12 parallel threads respectively (Timings computed on an Intel Xenon E5-2687W v4 Processor @3.0 GHz).

Figure 9 .
Figure 9. Example 4: Top row: Unfitted DG solution to the convection-diffusion problem on curved meshes with ℎ = 0.2, 0.1 and 0.05 for polynomial degree  = 4, resulting in 3375, 13530 and 53565 dofs respectively.Bottom row: Unfitted weak Trefftz DG solution with the same discretisation parameters and on the identical meshes resulting in 2025, 8118 and 32139 dofs, respectively.The discrete colour scale emphasises small oscillations around zero.

Figure A. 1 .
Figure A.1.Sketch of configuration in proof of Lemma 4.3.
=0 with diam( ′ ) ≤  diam(  ),  = 0, . . .,   ′ , where   ′ and  are uniformly bounded, satisfying (i) There exists a sub-element  ′  of  ′ with  planar facets meeting at a vertex  0  ∈  ′ , such that  ′  is star-shaped with respect to  0  and ℎ  ′ Assumption 2.1.We assume that for every  ′ ∈ ̃︀  ℎ there holds (a) There are two balls   ′ ⊂  ′ ⊂   ′ , such that  ′ is star shaped with respect to the ball   ′ and diam(  ′ )/ diam(  ′ ) 1. (b) The element boundary can be divided into mutually exclusive subsets {  } 1.A corresponding relaxed version of Assumption 2.1 is the following.Assumption 2.2.For every  ∈ ̃︀  ℎ , we assume that  is a Lipschitz domain and that it can be decomposed in   non-overlapping elements { ′ }, with   uniformly bounded and for each element  ′ the assumptions in Assumption 2.1 holds.Lemma 2.3.Assumption 2.1 directly implies Assumption 2.2 with   = 1 and { ′ } = { } for every  ∈ ︀  ℎ .Further, it immediately follows that for every sub-element  ′ of an element  ∈ ̃︀  ℎ with ̃︀  ℎ fulfilling Assumption 2.2 there are two balls   ′ ⊂  ′ ⊂  ⊂   such that  ′ is star shaped w.r.t. the ball   ′ and diam(  )/ diam(  ′ )   .