Trefftz Discontinuous Galerkin discretization for the Stokes problem

We introduce a new discretization based on the Trefftz-DG method for solving the Stokes equations. Discrete solutions of a corresponding method fulfill the Stokes equation pointwise within each element and yield element-wise divergence-free solutions. Compared to standard DG methods, a strong reduction of the degrees of freedom is achieved, especially for higher order polynomial degrees. In addition, in contrast to many other Trefftz-DG methods, our approach allows to easily incorporate inhomogeneous right hand sides (driving forces) by using the concept of the embedded Trefftz-DG method. On top of a detailed a priori error analysis, we further compare our approach to standard discontinuous Galerkin Stokes discretizations and present numerical examples.


Introduction
Trefftz methods, originating from the work of E. Trefftz [1], use approximation spaces that lie in the kernel of the target differential operator.Compared to standard approximation spaces, e.g.polynomial spaces, the corresponding Trefftz approximation spaces are constructed in such a way that they still have similar approximation properties as the original space, but, due to the kernel property, have a significantly reduced number of degrees of freedom.Trefftz methods provide a natural way to reduce degrees of freedom in discontinuous Galerkin (DG) methods by replacing the element-wise basis functions of the DG method by Trefftz basis functions.
So far, the Trefftz-DG methodology has been restricted to a small set of particular PDEs due to the need to explicitly construct basis functions for the corresponding PDE-dependent Trefftz-DG spaces and the restrictions mentioned above.Recently, a simple way to circumvent the explicit construction of Trefftz spaces for the Trefftz-DG method has been introduced with the embedded Trefftz-DG method in [14].This approach allows an easy construction of Trefftz-DG spaces for a broad class of PDEs, e.g.problems with differential operators that have not been considered in the Trefftz-DG literature before or PDEs with varying coefficients.Further, the approach allows for a generic construction of particular solutions to handle inhomogeneous source terms.
Trefftz methods for the Stokes equations that can be found in the literature are spectral-type methods.Explicit Trefftz basis functions have been constructed in [15,16].In [17] a so-called 'quasi-Trefftz spectral-method' is presented, that solves an eigenvalue problem on an encompassing domain to construct Trefftz-like functions and can also treat the Stokes problem with a source term.In [18] the Stokes problem is considered in two dimensions and singular solutions are introduced into a collocation Trefftz methods to deal with corner singularities.
While there is, to the best of our knowledge, no Trefftz-DG method for the Stokes equations in the literature, standard DG discretizations for the Stokes equations have been well-established for decades.As one way to reduce the computational costs of DG methods, hybridization, see [19], has become very popular leading to the class of hybrid DG methods; with several hybrid DG methods developed also for the Stokes equations.However, in the search for efficient discretizations based on DG formulations, Trefftz-DG and hybrid DG methods are incompatible.By incompatible we don't mean that a combination of both methods is impossible, but that a combination of both does not give any significant improvement over one of the two basic methods.That being said, instead of reducing the number of basis functions, non-polynomial Trefftz type functions have been used successfully in a hybrid setting to enrich the element-wise basis in [20].Since we focus on Trefftz-DG methods in this work we skip the large amount of literature on Stokes discretizations based on hybrid DG formulations in this introduction.Note, however, that several (hybrid) DG formulations are discussed in Section 5.3.

Main contributions and structure of the article
In this work, we extend the (embedded) Trefftz-DG methodology to the Stokes problem.We introduce a DG method with local basis functions that are polynomials and solve the Stokes problem pointwise.This extends the approach introduced in [21,22] where a DG method using element-wise divergence-free polynomials are presented for the Stokes and Navier-Stokes equations.In view of the Trefftz methodology the approach in [21,22] can be seen as a Trefftz-DG method with respect to the mass conservation equation, while the Trefftz-DG discretization considered in this paper includes also the momentum equation in the construction of the Trefftz-DG space.
The analysis of the discretization reveals that higher-order pressure functions are locally determined by the velocity and only the piecewise constant pressures appear as Lagrange multipliers in the Stokes saddle point problem.We prove the stability of the saddle point problem and derive a priori error bounds in the energy as well as the L 2 -norm of the velocity.In contrast to most other works on Trefftz-DG methods we allow inhomogeneous right-hand side source terms in the method and its analysis.While the method for treating inhomogeneous right-hand sides have already been introduced in [14], to the best of our knowledge, this is the first time a rigorous error analysis is provided for this generic approach in dealing with inhomogeneities.
The article proceeds as follows: We start with some preliminaries in Section 2 and then introduce a DG and the corresponding Trefftz-DG method for the Stokes problem in Section 3. The analysis is carried out in Section 4. Numerical results are presented in Section 5. Final comments are given in Section 6.

Preliminaries
Consider an open bounded Lipschitz domain Ω ⊂ R d with d = 2, 3.The Stokes equations determine a velocity u and a pressure p such that where f, g are external body forces and ν > 0 is the dynamic viscosity.For the ease of presentation we only consider homogeneous Dirichlet boundary conditions.The method (and its analysis) can be extended to more general boundary conditions using standard DG techniques, as demonstrated in the numerical examples.The weak formulation of the problem (2.1) is then given by: Find (u, p) where L 2 0 (Ω) is the space of square-integrable functions with a zero mean value.The (Trefftz-) DG formulation introduced in this work is set on a sequence of shape regular simplicial triangulations T h of a polygonal domain Ω.We denote by F h the set of facets in the mesh T h (edges/faces in 2d/3d, respectively) and define F i h as the subset of interior facets.With respect to a triangulation T h we denote by h the piecewise constant field of local mesh sizes with h| T = h T = diam(T ) for T ∈ T h .With abuse of notation, when h is evaluated on facets F ∈ F h we set h| F = h F = diam(F ), and also denote h = max T ∈T h h T as the maximum mesh size when h appears as a scalar quantity.
We denote by P k (S) the space of polynomials up to degree k on an domain S and set P k = P k (R d ).Further, for k < 0 we set P k = {0}.With P k (T h ) and P k (F h ) we denote the broken, i.e. element-or facet-wise polynomial space, such that for instance for v ∈ P k (T h ) it holds v| T ∈ P k (T ), for all T ∈ T h .We use an according notation for other broken function spaces as well, e.g.H 2 (T h ) denotes the space of functions with a local element-wise H 2 regularity.Due to its frequent appearance, we abbreviate the broken spaces on the mesh T h by setting P k = P k (T h ) (and similarly for other broken function spaces) when the mesh is clear from the context.
The L 2 (S) inner product over a domain S is denoted by (•, •) S .Specifically, we introduce the following notation for certain inner products: (•, •) T h sum of element-wise L 2 -inner products over all mesh elements, (•, •) ∂T h sum of facet-wise L 2 -inner products over the boundary of mesh elements, (•, •) F h sum of facet-wise L 2 -inner products over the facets of the mesh.
The unit outer normal is denoted by n.Correspondingly we use ∥ • ∥ S for the L 2 -norm on a domain or set of elements/facets S and define By Π k S : L 2 (S) → P k (S) we denote the L 2 projection into the scalar-valued polynomial space of order k on S. With abuse of notation we use the same symbol for vectorial L 2 -projections and denote by Π k the element-wise L 2 -projection on P k (T h ) or [P k (T h )] d , respectively.
Let T and T ′ be two neighboring elements sharing a common facet F ∈ F i h where T and T ′ are uniquely ordered.On F the functions u T and u T ′ denote the two limits of a discrete function from the different sides of the element interface.The vectors n T and n T ′ are the unit outer normals to T and T ′ .For F ∈ F i h we then define the jump and the mean value by and while on facets on the domain boundary we set Finally, we use the notation A ≃ B when there are constants c, C > 0 independent of h and ν such that A ≤ CB and B ≤ cA.Similarly, we also use the symbols ≲ and ≳.

Trefftz-DG for Stokes
In this section we derive the Trefftz-DG formulation for the Stokes problem.For that, we first introduce a standard DG method and then include the Trefftz spaces in the formulation.In Section 3.4 we show a possible way to implement the Trefftz spaces and construct local particular solutions to treat the inhomogeneous problem using the embedded Trefftz-DG method, see [14].

The underlying DG Stokes discretization
As a basis for the following, we consider the established symmetric interior penalty DG discretization of the Stokes problem (2.1), cf.e.g.[23,Section 6.1.5]: with the bilinear forms where the interior penalty parameter α = O(k 2 ) is chosen sufficiently large and we used the notation ∂ n w := ∇w • n.
In this work, we want to introduce a Trefftz-DG method by reducing the finite element space to a subspace, the Trefftz-DG space.For the underlying DG space, we introduce the notation For simplicity, we omit the superscript k and write X k h = X h and X k h (T ) = X h (T ).Further, we introduce the bilinear form K h on the product space X h by (3.2) Then (3.1) also reads as: Find (u h , p h ) ∈ X h such that

The Trefftz-DG Stokes discretization
The main idea of the Trefftz-DG method is to select a proper (affinely shifted) lowerdimensional subspace of X h that allows to reduce the computational costs of the underlying DG discretization without harming the approximation quality too much.
To this end, we pick the affine subspace of X h that fulfills the Stokes equations (2.1) pointwise inside each element up to data approximation: Again, when the relation to the mesh T h is clear from the context, we will abbreviate ).The definition of a local Trefftz space T k f,g (T ) is correspondingly based on X k h (T ).Note that we allow for the case k = 1 for which the first constraint −∆u h + ∇p h = Π k−2 f is automatically fulfilled as u h is piecewise linear and p h is piecewise constant while the r.h.s.projection maps to zero.Only the constraint − div u h = Π 0 g remains non-trivial in this case.
To work with linear spaces we decompose T k f,g into a (non-unique) particular solution to the Trefftz constraints (u par h , p par h ) ∈ T k f,g and a linear space 0,0 (T h ), typically denoted as the Trefftz-DG space.Again we omit the superscript k and simply write T k f,g = T f,g and T k = T = T(T h ) = T 0,0 (T h ).We prove below in Section 3.3 that a particular solution always exists, cf.Lemma 1, and can be constructed in a straight-forward manner, cf.Section 3.4.
This yields the Trefftz problem: Equivalently, to highlight the homogenization, we can also write: Find (u hom h , p hom h ) ∈ T(T h ) so that for all (v h , q h ) ∈ T(T h ) there holds

Trefftz constraints and the dimension of T
In this subsection we want to discuss some properties of the Trefftz spaces T f,g and T.
First, as we will typically have g = 0, we note that the velocities in T f,0 (and T) will be divergence-free on each element with the additional constraint of −∆u h + ∇ p h = Π k−2 f , making T f,0 (and T) a subspace of the solenoidal vector fields considered in [21,24].Next, we prove that a particular solution to the Trefftz constraints always exist.Furthermore, we state the dimension of the local Trefftz space.Lemma 1.The pointwise Stokes operator L : Proof.We give the proof for d = 3.The proof for d = 2 follows similar lines.To prove surjectivity, first note that div[P k ] d = P k−1 and hence it remains only to show that We prove this in two steps.
1. First we note that [ , where for the first equality we have used the surjectivity of the Laplace operator, see Corollary 12 in the appendix.Now using that curl The other direction, curl d is obvious and hence equality holds.We used the notation for the sum of two spaces.Note that this is not an orthogonal decomposition, and that the spaces considered here have a non-trivial intersection.In total this gives that for every where we used that to every u ∈ [P k ] d with div u = 0 we can find a w ∈ [P k+1 ] d such that u = curl w (and vice versa).Remark 1.We note that in the definition of the space X k h (T h ) -which the Trefftz space T k (T h ) is a subspace of -we factored out the (globally) constant pressure (R).This however does not affect the local spaces X h (T ) and T(T ) which are defined without consideration of the global constant pressure.
Remark 2. The formula derived in (3.7) shows the number of degrees of freedom per element (ndofs) for the Trefftz-DG space.We observe that the ndofs are exactly reduced by the number of scalar constraints that are imposed.The Stokes equations can be seen as the momentum equation formulated on the subspace of divergence-free functions where the pressure takes the role of the corresponding Lagrange multiplier.The costs of these additional unknowns are directly removed by the corresponding divergence-free constraint so that the ndofs remaining coincide with those of harmonic vector polynomials up to degree k.Note however that the remaining system is still a saddle-point problem, cf. also Remark 3. Let us now take a look at the dimension reduction.While the dimension of X h grows cubic and quadratic with respect to the polynomial degree k for d = 3 and d = 2, respectively, the dimension of T grows only quadratic and linear for d = 3 and d = 2.In Table 1 we present concrete numbers for k ∈ {1, .., 6} to give an idea of the reduction.This dimension reduction gives the same complexity in k that is also obtained in hybrid DG methods after static condensation.The pressure scaling is different between the first seven and the last three basis functions.Note, that these basis functions are obtained from the generic approach of the embedded Trefftz-DG method [14], cf.Section 3.4, and hence do not offer a complete insight into an available structure of the Trefftz space such as a clean decomposition into lower and higher order basis functions.Nevertheless, we observe that velocity and pressure functions are coupled for most basis functions except for the three lowest order Trefftz basis functions: The first two basis functions (in the upper row) have a zero pressure and are constant and linear, respectively, and divergence-free and the last (in the lower row) basis function corresponds to a zero velocity and a constant pressure.
Remark 3 (Lowest order subspaces).Due to the Trefftz constraint, basis functions of the Trefftz-DG space can no longer be separated into velocity and pressure functions in general.See Fig. 1 for an example set of basis functions for k = 2.However, for the lowest-order polynomial degrees, some degrees of freedom can be associated only with velocities and others only with pressure functions.The piecewise constant pressure functions P 0 (T h ) are not seen by the Trefftz constraints as the pressures only appear as gradients, i.e. these degrees of freedom remain in the Trefftz space.From the velocity space, the velocity-pressure coupling is removed (from the space) for all linear velocity functions, i.e. [P 1 (T h )] d , as here the Laplacian vanishes.The element-wise divergence constraint removes exactly one dof per element from this space (div[P 1 (T h )] d = P 0 (T h )).Especially, piecewise constant functions [P 0 (T h )] d remain in the Trefftz-DG space.Finally, let us note that the problem formulated in the Trefftz-DG space is still of saddle-point form with the Lagrange multiplier space being only the space of piecewise constant functions.

Implementation via embedded Trefftz method
So far we characterized the Trefftz-DG space only implicitly through the Trefftz constraints.In principle a basis for the homogeneous Trefftz space T 0,0 could be constructed, see e.g.[15,16].An alternative approach, based on the implicit characterization, is the use of the embedded Trefftz-DG method, presented in [14].The method allows to construct an embedding T : T(T h ) → X h (T h ) based on the Trefftz constraint equations.Then, all essential operations can be done by exploiting the underlying DG space.In this approach it is also straightforward to find a generic element-wise particular solution (u par h , p par h ) ∈ X h (T h ) required for the homogenization of the Trefftz problem.In the following, we shortly recap the procedure of the embedded Trefftz method for the Stokes setting, but refer to [14] for more details.
Let us define the matrix and vector associated to the discrete linear system (3.1) for i, j = 1, . . ., N and (ϕ i , ψ i ) ∈ X h (T h ) a set of basis functions with N the number of degrees of freedom.The Trefftz embedding can then be represented by the kernel of the matrix ).The Trefftz embedding matrix T then allows to characterize a Trefftz basis function as a linear combination of polynomial DG basis functions.This allows the reduction of the size of global and local finite element matrices reducing the overall costs associated with the solution of the arising linear systems.

A-priori error analysis
In this section we derive a-priori error bounds for the Trefftz-DG method.To this end we define the norms For the overall problem we define the norm In addition, as common in the analysis of DG methods, we introduce the stronger norms We note that on T(T h ) the norms ∥(•, •)∥ T and ∥(•, •)∥ T, * are equivalent with constants independent of h and ν.
For the analysis we will repeatedly rely on discrete trace and inverse inequalities, see e.g.[23,, and optimality of the Π 0 -projection on mesh elements and on faces, see e.

Saddle-point structure and space decomposition
In this section, we introduce some structures to analyse and exploit the saddle-point structure of the variational problem.We first note that the space {0}×P 0 (T h )/ R, with P 0 (T h ) the space of piecewise constant pressures form a subspace of T(T h ).Accordingly, we can introduce the decomposition into two subspaces: one that contains all velocity functions and the high order (≥ 1) pressure functions and a second subspace only consisting of the zero velocity and low-order pressures We recall T = T(T h ), H = H(T h ) and L = L(T h ).Corresponding projection operators for L and H are denoted by Π L and Π H = id −Π L , respectively, so that for (u Note that the decomposition effectively operates on the pressure only and is hence L 2 -orthogonal on the pressure.We can then introduce norms for the subspaces H and L by These norms are natural in the sense that Next, we observe that the higher-order pressures in H are controlled by the velocity due to the Trefftz constraint while on L norms are completely determined by the lowest order pressure contributions.This norm control is highlighted in the following lemma for general functions in T. Proof.From the Trefftz constraint we have element-wise ∇p h = ν∆u h for (u h , p h ) ∈ T. Exploiting this, together with standard local inverse inequalities, yields Similarly we also have ∥(u h , p h )∥ 2 L = ν −1 ∥p h ∥ 2 0,h which follows from u h = 0 and ∇ p h = 0 for a function (u h , p h ) ∈ L. Now, splitting (u h , p h ) ∈ T into its components from L and H and exploiting the previous statements we have For the stability analysis, we will make use of the special saddle-point structure, that develops from the subspace splitting T = H ⊕ L, to show that the bilinear form K h (•, •) satisfies a discrete inf-sup condition.To this end we can split the bilinear form K h (•, •) into its contributions that we obtain by restricting to the corresponding subspaces L and H.This gives where we used that K h (Π L (u h , p h ), Π L (v h , q h )) = 0 and Π L (u h , p h ) = (0, Π 0 p h ) so that the velocity contributions of Π L (u h , p h ) vanish.

Continuity and kernel-coercivity
Lemma 3. The bilinear form K h (•, •) is continuous, i.e. we have for all (u, p), Proof.The bound for a h (u, v) ≤ ∥u∥ 1,h, * ∥v∥ 1,h, * follows by standard arguments, see e.g.[23,Lemma 4.16].Where as for the mixed terms we compute Using optimality of the Π 0 -projection and a discrete trace inequality we get and by similar reasoning for the volume term ∥p∥ 0 ≤ ∥h∇p∥ 0 + ∥Π 0 p∥ 0 .Finally, we use [26, Remark 1.1], with the fact that Π 0 p ∈ P 0 (T h )/R, to obtain ∥Π 0 p∥ 0 ≤ ∥h , finishing the proof.Next, we show coercivity of K h | H×H ("the upper left block").Lemma 4. Let the stabilisation parameter α > 0 be sufficiently large, then for (u h , p h ) ∈ H(T h ) there holds for a constant independent of the mesh size h and ν.
Proof.Using decomposition (4.3) we first have For the first term a h (u h , u h ) it is well-known that there holds coercivity on P k (T h ) with respect to the discrete norm ν 1 2 ∥•∥ 1,h for α > 0 sufficiently large, see for example [27] or [23,Lemma 4.12].We denote by c α a sufficiently lower bound such that for α ≥ c α coercivity with corresponding coercivity constant c a is guaranteed.
Thus it remains to bound the latter term.Let ε > 0, then Cauchy-Schwarz and Young's inequality (2ab Using inverse inequalities for polynomials and recalling that p h = (id −Π 0 )p h (as (u h , p h ) ∈ H) we have with the element-wise Trefftz constraint and hence, for α sufficiently large we obtain the claim with the norm equivalence of Lemma 2.
For the LBB-stability we use a stability result of an auxiliary (low order) discrete problem: Given ph Here, the bilinear form a ν=1 h (•, •) is to be understood as the bilinear form that is obtained from a h (•, •) when setting ν to 1, so that the auxiliary problem is independent of ν.In this auxiliary problem, r h and r h are the Lagrange multipliers for enforcing the pointwise divergence-constraint div w h = 0 and the normal jump condition [[w h • n]] = ph , respectively.For ph = 0 the solution w h is normal-continuous, i.e.H(div)conforming.More precisely we would have with BDM 1 the Brezzi-Douglas-Marini space, cf.[28, Section 2.3.1].For general ph we instead have w h ∈ BDM −1 (T h ) = [P 1 (T h )] d with BDM −1 (T h ) the Brezzi-Douglas-Marini element with broken H(div)-continuity.Lemma 5.The velocity solution w h ∈ [P 1 (T h )] d of the auxiliary problem (4.6) is element-wise divergence-free div w h = 0 and further there holds Proof.Very closely related statements can be found in the literature of hybridized mixed and H(div)-conforming DG and hybrid DG methods, cf.e.g.[29,30].For completeness, we give a proof here.We analyse the twofold saddle-point problem (4.6) with the bilinear forms A similar analysis has been carried out in [31,Lemma 4.3.6].From element-wise H(div)-interpolation, cf.e.g.[28, Proposition 2.5.1], and standard scaling arguments we have that [[w h • n]] can be matched with scalar data on facets in the following sense: For all s h ∈ P 0 (F h ), and there holds sup The kernel of e h (•, •) are H(div)-conforming functions in BDM 1 (T h ).For this subspace, we have the following well-known inf-sup result, cf.[30,32]: Hence, we have inf-sup stability of d h (•, •) on ker e h (i.e.functions with normal continuity), inf-sup stability of e h (•, •) and coercivity of a h (•, •) on ker d h ∩ ker e h (i.e.divergence free and normal continuous functions).Here we use a slight abuse of notation associating the kernel of the bilinear form with the kernel of the corresponding operator.This implies stability of the twofold saddle-point problem, cf.e.g.[33] and thus the stability estimate in the claim.
We can now state the LBB-stability result for the Trefftz-DG problem: Lemma 6.There holds inf for a constant c b > 0 independent of h and ν.
Proof.The first equivalence immediately follows by Lemma 2. Now let (v h , q h ) ∈ L be arbitrary, i.e. v h = 0 and q h ∈ P 0 .We define ph = −[[q h ]] ∈ P 0 (F i h ) on F i h and will construct a suitable velocity field u h that allows to control ph .For this we use the element-wise divergence-free space ) there holds on each element that ∆v h = 0 (since v h is linear) and div v h = 0, the tuple (v h , 0) ∈ H, thus BDM −1 0 × {0} ⊂ H.This gives inf Now choose u h = w h with w h being the solution of the auxiliary problem (4.6) with data (used for the right hand side) ph .This gives and thus by the stability estimates of Lemma 6 for u h = w h we get which concludes the proof.

Inf-Sup-Stability
Theorem 7.For (v h , q h ) ∈ T(T h ) there holds sup for constants c T , c * T independent of h and ν and hence the Trefftz-DG problem (3.5) admits a unique solution that continuously depends on the data.
Proof.Together with Lemma 4 and Lemma 6 and the continuity (4.4) we have shown that K h (•, •) is inf-sup-stable on T and thus well-posed, see e.g.[23,Lemma 1.30].

Quasi-best approximation and Aubin-Nitsche
In the following, we consider a best approximation result in the Trefftz space and discuss an Aubin-Nitsche-like result for the L 2 -error of the velocity.Lemma 8. Let (u, p) ∈ [H 2 (T h )] d × H 1 (T h )/ R be the solution of the Stokes problem (2.2) and (u h , p h ) ∈ T f,g (T h ) be the discrete solution to (3.5).Then, there holds where we made use of the consistency of the formulation, and continuity (4.4).Finally, the claim follows with an application of the triangle inequality.
To show error estimates in the L 2 -norm we use a duality argument in the style of Aubin-Nitsche.This requires additional regularity for the solution of the Stokes problem, see e.g.[34] for regularity results of the Stokes problem mentioned below.Theorem 9. Assume the domain boundary to be sufficiently smooth or Ω to be convex so that L 2 -H 2 -regularity holds and let (u, p) Due to the consistency of the primal problem, we can add the following expression which adds up to zero for any (w h , r h ) ∈ T: With w h = Π BDM,1 w, the BDM 1 -interpolant of w, we have div w h = Π 0 div w = 0 and hence with r h = Π 0 r there holds (w h , r h ) ∈ T 1 ⊂ T. Now applying standard interpolation results for the BDM-interpolator as well as exploiting the L 2 -H 2 -regularity yields the bound We can then conclude with continuity (4.4) to obtain Dividing by ∥u − u h ∥ Ω then yields the claim.

Approximation and a-priori error bounds
We conclude the stability analysis with the following optimal error estimate.
Hence (v h , q h ) ∈ T f,g .Now we bound: The first part, I, can be directly bounded by the right-hand side of (4.12) due to the approximation properties of the averaged Taylor polynomial, cf.[35, Lemma 4.3.8].
We hence turn our attention to II.As the terms in II are discrete functions (piecewise polynomials), we can apply inverse inequalities to obtain Then we have by the properties of the averaged Taylor polynomial, cf.[35, Lemma 4.3.8],[36].Now, using (in order), Π C v ∈ C 0 (Ω), a trace inequality, triangle inequalities, the interpolation properties of T k and Π C , and (local) H 1 -continuity of T k and Π C we have Similarly we also have We note that (u − u ′ , p − p ′ ) solves the Stokes problem for the data (f Exploiting linearity and stability of the continuous problem (see [23 what concludes the proof. R be the solution of the Stokes problem and (u h , p h ) ∈ T f,g (T h ) be the discrete solution to (3.5).Then, there holds Proof.This is a direct consequence of the previous Lemma and Theorem 9.

Exact solution
For the first numerical example, we consider the domain Ω = (0, 1) d and the solution for d = 2, 3, respectively.We then compute the numerical solution with the given right hand side f = −ν∆u + ∇p (using above exact solution), g = 0 and homogeneous Dirichlet boundary data.Further, for simplicity, we fix ν = 1.The penalty parameter is chosen as α = 20.
In Figure 2 we compare the L 2 -error of the Trefftz-DG method (3.5) (marked with T k ) to the solution of the standard DG method (3.3) (marked with P k ) for orders k = 2, 3, 4. The error of the velocity and pressure approximation with respect to the exact solution is measured in the L 2 -norm.As expected, see (4.12) and (4.11), we observe that the solution of the Trefftz-DG method converges with the same (optimal) order as the solution of the standard DG method for both the velocity and the pressure error.

Moffatt eddies
In [40] Moffatt presents an example on a wedge that produces an infinite amount of eddies that differ greatly in magnitude, making it a challenging numerical example.This benchmark has also been considered in [41].We consider a triangular domain, shown in Figure 3 Example in 2 dimensions (3.5) the following boundary terms

Example in 3 dimensions
where u D is the Dirichlet boundary data.We impose an inflow on the part of the boundary given by the line y = 0 with a parabolic velocity profile and a no-slip condition (homogeneous Dirichlet boundary) is imposed on the remaining boundary, i.e. the boundary data u D is given by The domain Ω is given by a wedge with a sharp angle of approximately 2α = 36.87• on the corner opposite of the boundary with the inflow.The velocity solution comprises an endless series of swirls, with each subsequent swirl being approximately 400 times less intense than the one preceding it.The series of swirls converges to (0, −2) Additionally, the pressure field exhibits two point singularities at (−1, 0) and (1, 0).In Figure 3 we show the numerical results on a mesh with 28 elements for k = 10.We see that the Trefftz-DG method is able to capture the sharp corner and the eddies, resolving four to five eddies, encompassing a scale range of 10 13 .

Comparison of computational effort to different Stokes discretizations
In this section, we want to compare the Trefftz-DG method with other popular discretizations for the Stokes problem with respect to the number of unknowns and their sparsity patterns that finally determine the computational costs to a large extent.For a comparison, we only consider simplicial meshes and mostly restrict to different variants of DG methods.We focus on a comparison of the total number of degrees of freedom (ndof), the ndof that remains after elimination of all interior unknowns (i.e. after static condensation) denoting them as coupling ndof (ncdof), and the resulting sparsity pattern of the methods (after static condensation) in terms of the non-zero entries in the resulting sparse system matrices (nnze).We include several hybrid Discontinuous Galerkin (HDG) methods, where the velocity couplings across element interfaces stemming from DG-terms are avoided by introducing additional facet unknowns for the velocities.As an exception, we also include the Taylor-Hood method (of arbitrarily high order) which has continuous velocity and pressure approximations.All methods share the same order of convergence in a discrete H 1 -type norm for the velocity and an L 2 -type norm for the pressure.The list of methods presented here is by no means complete, for additional material on HDG methods we refer to e.g.[42][43][44], and for further exploration, we refer to the relevant literature, including the books [23,25,28] and the references therein.
First, we consider discontinuous Galerkin methods, with basis functions that are discontinuous across inter-element boundaries, including the method presented in this work.The methods can be considered the closest relatives to the Trefftz-DG method as the inter-element regularity of the solution is solely enforced in the bilinear form.This makes them very flexible in the mesh choice, allowing easy for e.g.polygonal meshes.We consider the following:

• Standard DG:
As a representative of a "standard" Discontinuous Galerkin method we consider the interior penalty discretization from Section 3.1, cf.especially (3.1), respectively from [23, Section 6.1.5]with u h ∈ [P k ] d and p h ∈ P k−1 , cf.Section 3.1.

• Trefftz-DG:
The Trefftz-DG method, presented in this work, uses the same variational formulation as the standard DG method, but with a different finite element space so that (u h , p h ) ∈ T k with T k as in (3.4).The method only has unknowns on the elements and static condensation cannot be applied.• Solenoidal DG: In [21,22] the velocity space is reduced to element-wise divergence-free polynomial functions.For the enforcement of normal-continuity, a facet variable for the pressure is then introduced leading to exactly divergence-free discrete solutions.Velocities

• H(div)-DG:
By considering discrete velocities from H(div)-conforming finite element spaces the pressure facet variable can be avoided while still obtaining exactly solenoidal solutions (after computation).This has been proposed in [30] with u h ∈ BDM k and p h ∈ P k−1 .• High order divergence-free (hodf) H(div)-DG: Exploiting the a-priori knowledge that velocity solutions in H(div)-DG are exactly divergence-free allows to reduce the velocity basis functions a-priorily to those with a piecewise constant divergence, u h ∈ {v ∈ [P k ] d | div v ∈ P 0 }, cf.[46].Note that the basis functions in BDM k generating a divergence in P 0 are not local and remain in the system.As pressure couplings across element boundaries vanish in an H(div)-conforming setting we can correspondingly remove the pressure variables of higher order so that w.r.t. the pressure only p h ∈ P 0 remains for the solution of the global linear system.Higher order pressures can again be reconstructed elementby-element in a simple post-processing and are not considered in this comparison.Compared to the previous H(div)-DG method this step can be considered as a version of static condensation.

• H(div)-HDG:
In the Rhebergen-Wells-HDG method the normal component is made continuous essentially twice, through a pressure facet variable and as one component of the weak continuity stemming from the viscosity term.This redundancy can be removed.
In [47] an H(div)-conforming space for the velocity is considered, u h ∈ BDM k with p h ∈ P k−1 and a tangential vector function on the facets u ∂T h with div BDM ℓ ∂T h = P 0 , cf. [46].Then, the interior normal-bubbles (which determine the higher-order divergence of the velocity) can be statically condensed alongside the higher-order pressures.Only the facet variables (normal dof of the BDM k ∂T h space and the tangential vector functions in P k τ (F h )) and the piecewise constant pressure functions in P 0 (T h ) remain in the global linear system.
• Projected jumps modification (pj) of the H(div)-HDG method: An improvement of the H(div)-HDG method is obtained when reducing the polynomial degree of the tangential facet variables by one degree using a projected jumps modification of the hybrid interior penalty formulation, cf.[47].• Highest order discontinuous facet modification (hodc) of the H(div)-HDG method: A second improvement of the H(div)-HDG method is obtained when normal continuity is relaxed by one degree leading to a special velocity space Then, the facet dofs of BDM ⋆,k coincide with those of BDM k−1 ∂T h and together with P k−1 τ (F h ) they suffice to couple velocity functions on neighboring elements, cf.[48,49].
Finally, we also include an H 1 -conforming method, here we consider: • Taylor-Hood: The Taylor-Hood method is an H 1 -conforming Galerkin method with u h ∈ [P k ∩ C 0 (Ω)] d and p h ∈ P k−1 ∩ C 0 (Ω) using bilinear and linear form as in the continuous weak formulation.By P ℓ = P ℓ • ⊕P ℓ ∂T h we denote the decomposition of the polynomial space P ℓ into interior bubbles P ℓ • = {v ∈ P ℓ | v| ∂T h = 0} and interface polynomials P ℓ ∂T h .For the Taylor-Hood element, the interior bubbles can be eliminated by static condensation and only the interface polynomials remain in the reduced linear system.
In Table 2 we summarized the considered methods in a table with an emphasis on the used discretization spaces and those that remain after static condensation.In Fig. 4 we illustrate the different degrees of freedom for some of the methods in two dimensions and k = 4. Tables B1 to B6 we show the corresponding results for all methods mentioned above.We observe that the Trefftz-DG method brings a significant improvement over the Standard DG method, especially in terms of non-zero entries in the matrix.Compared to the Hybrid DG methods it is competitive in terms of ndof and ncdof and only slightly worse in terms of nnze when compared to the optimized HDG variants.It even compares quite well to the popular Taylor-Hood method for higher orders.

Conclusion & Outlook
In this paper we introduced a new Stokes discretization based on local basis functions that are Trefftz, i.e. that fulfill the Stokes equations pointwise (w.r.t. to an approximated r.h.s.).This leads to a strong reduction of unknowns compared to standard DG methods.To construct the corresponding basis we use the embedded Trefftz-DG method which allows us to deal with inhomogeneous forces and sources.The method is analyzed and a priori error bounds are derived.To the best of our knowledge, this is the first Trefftz-DG discretization for the Stokes problem.Crucial and new components in the analysis are a splitting of the pressure unknowns into higher order and lower order parts and the analysis for inhomogeneous forces and sources.Especially the latter will also be useful for the analysis of (embedded) Trefftz-DG methods for other PDEs.
A topic that we have not addressed in this work is pressure robustness.Although discrete solutions of the presented scheme will be pointwise divergence-free, discrete solutions will in general not be normal-continuous or pressure robust.We leave improvements of the proposed scheme in that regard for future research.

Figure 1
Figure 1Example basis functions of a Stokes-Trefftz space for k = 2 on a triangle.The coloring corresponds to the pressure value while the arrows indicate the velocity.The pressure scaling is different between the first seven and the last three basis functions.Note, that these basis functions are obtained from the generic approach of the embedded Trefftz-DG method[14], cf.Section 3.4, and hence do not offer a complete insight into an available structure of the Trefftz space such as a clean decomposition into lower and higher order basis functions.Nevertheless, we observe that velocity and pressure functions are coupled for most basis functions except for the three lowest order Trefftz basis functions: The first two basis functions (in the upper row) have a zero pressure and are constant and linear, respectively, and divergence-free and the last (in the lower row) basis function corresponds to a zero velocity and a constant pressure.
be the solution of the Stokes problem (2.2) and (u h , p h ) ∈ T f,g (T h ) be the discrete solution to (3.5).Then, there holds ∥u − u h ∥ Ω ≲ h∥(u h − u, p h − p)∥ T, * .(4.11) Proof.We apply the Aubin-Nitsche trick and pose the auxiliary adjoint Stokes problem L(w, r) = (u − u h , 0) on Ω.With u − u h ∈ [L 2 (Ω)] d and the assumed L 2 -H 2 -regularity, we have (w, r) ∈ [H 2 (Ω)] d × H 1 (Ω).Now, from the (adjoint) consistency of the bilinear forms a h (•, •) and b denote the elementwise (potentially vectorial) averaged Taylor polynomial (averaged on a proper inner ball of each element) of degree k, cf.[35, Section 4.1-4.3].Then, we set v h = T k u ′ and q h = T k−1 p ′ and have on each element on the right, with non-homogeneous Dirichlet boundary conditions on one side, and zero source term, f = (0, 0).To implement the non-homogeneous Dirichlet boundary conditions we add to the right hand side of ourTrefftz-

Figure 2
Figure 2 Numerical results for the DG and Trefftz-DG method for the two dimensional example on the top row and three dimensional example on the bottom.The gray (solid and dashed) lines indicate the expected convergence rates.

Figure 3
Figure3 Moffatt eddies with k = 10 on the left, including a zoom on the bottom eddies on the bottom left.On the right we show the computational mesh.

Figure 5
Figure 5 Comparison of ndof, ncdof and nnze for a selection of methods in Section 5.3 for the 2D and 3D Stokes problem on the displayed mesh.

Table 1
Dimensions of the local finite element spaces X h (T ) and T(T ).
where a α=cα h (u h , v h ) is the bilinear form with interior penalty stabilization parameter c α for which there holds a α=cα h (u h , u h ) ≥ c a ν∥u h ∥ 2 1,h .As u h is pointwise divergence-free the volume contribution in b h (•, •) vanishes and we have [45]pressures p h ∈ P k (F h ) are assumed to appear in the global linear system, i.e. no static condensation is applied.Volume pressure unknowns do not appear in the global linear system (neither as ndof nor ncdof), but can be reconstructed in element-by-element post-processing.•Rhebergen-Wells-HDG:In[45]Rhebergenand Wells introduce d + 1 scalar facet variables of degree k to enforce normal continuity and the weak continuity of the velocity vector.All element interior unknowns for velocity and pressure can then be eliminated by static condensation.
Next, we will consider H(div)-conforming methods, i.e. methods based on a discrete velocity space that is normal-continuous across element interfaces.We recall that BDM k = [P k ] d ∩ H(div; Ω) is the Brezzi-Douglas-Marini space, see e.g.[28, Section 2.3.1].