{'fraction_non_alphanumeric': 0.04186341811832864, 'fraction_numerical': 0.022538788143327645, 'mean_word_length': 4.282330748270803, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Discontinuous Galerkin methods for nonvariational problems

We extend the finite element method introduced by Lakkis and Pryer [2011] to approximate the solution of second order elliptic problems in nonvariational form to incorporate the discontinuous Galerkin (DG) framework. This is done by viewing the NVFEM as a mixed method whereby the finite element Hessian is an auxiliary variable in the formulation. Representing the finite element Hessian in a discontinuous setting yields a linear system of the same size and having the same sparsity pattern of the compact DG methods for variational elliptic problems. Furthermore, the system matrix is very easy to assemble, Thus this approach greatly reduces the computational complexity of the discretisation compared to the continuous approach. We conduct a stability and consistency analysis making use of the unified framework set out in Arnold et. al. [2001]. We also give an apriori analysis of the method. The analysis applies to any consistent representation of the finite element Hessian, thus is applicable to the previous works making use of continuous Galerkin approximations.


Introduction
Nonvariational partial differential equations (PDEs) are those which are given in the form where X : Y = trace(X Y ) is the Frobenious inner product between matrices. If the matrix A is differentiable then there is an equivalence between this problem and its variational sibling Rewriting in this form is sometimes undesirable. For example, if the coefficient matrix A has near singular derivatives the problem will become advection dominated and possibly unstable for conforming finite element methods. There is a wealth of material on the treatment of advection dominated problems [EG04,ESW05,c.f.]. If A is not differentiable then the problem has no variational structure. In this case standard finite element methods cannot be applied.
In a previous work [LP11] a finite element method for the approximation of the nonvariational problem (1.1) was introduced. This involved the introduction of a finite element Hessian represented in the same finite element space as the solution (modulo boundary conditions). The applications of the discrete representation of a Hessian of a piecewise function are becoming broader, for example, it can be used to drive anisotropic adaptive algorithms [AV02, VMD + 07], as a notion of discrete convexity [AM09] and in the design of finite element methods for nonlinear fourth order problems [Pry12]. We are particularly interested in nonvariational PDEs due to their relation to general fully nonlinear PDEs which are of significant current research. There have been finite element methods presented for this general class of problem for example in [Böh08] the author presents a C 1 finite element method shows stability and consistency (hence convergence) of the scheme which requires a high degree of smoothness on the exact solution. In [FN09b,FN09a] the authors give a method in which they approximate the general second order fully nonlinear PDE by a sequence of fourth order quasilinear PDEs. This is reminiscent of the vanishing viscosity method introduced for classically studying first order fully nonlinear PDEs. Efficiency of any method used to approximate a problem such as this is key. Each of the methods are computationally costly due to their reliance on C 1 finite elements [Böh08,FN09b] or mixed methods [FN09a].
In [BS91] a generic framework was set up to prove convergence of numerical approximations to the solutions of fully nonlinear PDEs. This involved constructing monotone sequences of approximations which are typically applied to finite difference approximations of the nonlinear problem [Obe06,c.f.]. The assumption of consistency made in the [BS91] framework are incompatible with finite element methods, however, an extremely important observation made in [JS11] is that the consistency condition may be weakened to incorporate the finite element case using a localisation argument (in the case of isotropic diffusion).
In this contribution we present a method for the discontinuous approximation of the linear nonvariational problem (1.1). We also present convergence analysis for a certain subclass of the nonvariational problems, those which are coercive. This allows us to use variational techniques to analyse the problem. We prove optimal convergence rates for the finite element solution in broken Sobolev norms. Note that the results presented here are immediately applicable to the method derived for the continuous case given in [LP11].
The algebraic formulation of the continuous approximation of the nonvariational problem requires the solution of large sparse (d + 1) 2 × N 2 linear system [LP11, Lem 3.3], where d is the dimension of the problem and N the number of degrees of freedom. Equivalently, using a Schür complement argument, this can be reduced to an N 2 full linear system. The reason that this system is full is due to the global nature of the L 2 (Ω) projection operator into a continuous finite element space. The motivation for extending the nonvariational finite element method into the discontinuous setting is the massive gain in computational efficiency over the continuous case. Indeed, due to the local representation of the projection operators in these discontinuous spaces we are able to make massive computational savings, in that the system matrix will become sparse and is the same size as that of a standard discontinuous Galerkin stiffness matrix.
To test the method numerically we make use of the finite element package Dune [BBD + 08a, BBD + 08b]. In this work we are interested in the asymptotic behaviour of the discontinuous approximation. In a subsequent work we will study the computational gains using the discontinuous framework presented over the continuous one given in [LP11], as well as exploit the powerful parallelisation capabilities of the package.
The rest of the paper is set out as follows: In §2 we formally introduce the model problem and give a brief review of known classical facts about nonvariational PDEs. In §3 we examine the discretisation of the nonvariational method in the discontinuous Galerkin framework, making use of the unified framework set out in [ABCM02] to derive a very general formulation of the finite element Hessian represented as a discontinuous object. We present some examples and examine the natural question of what happens when we try to eliminate the finite element Hessian from the formulation. In §4 we look at the consistency and stability of the finite element Hessian and present our main analytical results of convergence. Finally, in §5 we detail a summary of extensive numerical experiments aimed at examining convergence and computational speed of the method presented.

Problem formulation
In this section we formulate the model problem, fix notation and give some basic assumptions. In addition we review the existence and uniqueness of the nonvariational problems. We begin by introducing the Lebesgue spaces α i and derivatives D α are understood in a weak sense. We pay particular attention to the cases k = 1, 2 and The model problem in strong form is: where the data f ∈ L 2 (Ω) is prescribed and L is a general linear, second order, uniformly elliptic partial differential operator. Let A ∈ L ∞ (Ω) d×d , we then define We assume that A is uniformly positive definite, i.e., there exists a γ > 0 such that for all x and we call γ the ellipticity constant.
Nonvariational PDEs are not as well studied as their variational brethren from a numerical analysis viewpoint. For the benefit of the reader we present a concise account of known results for strong solutions of this class of problem.
2.1. Definition (strong solution). A strong solution of (1.1) is a function u ∈ H 2 (Ω) ∩ H 1 0 (Ω), that is a twice weakly differentiable function, which satisfies the problem almost everywhere.
2.2. Theorem (existence and regularity of a strong solution of (1.1) [GT83, Thm 9.15]). Let Ω ⊂ R d be a C 1,1 domain. Suppose also that A ∈ C 0 (Ω) d×d and f ∈ L 2 (Ω) such that the problem is uniformly elliptic. Then (2.9) has a unique strong solution. There also exists a constant independent of u such that 2.3. Remark (less regular solutions). Note that the theory of viscosity solutions has been developed for non classical solutions of (2.6) if the problem data does not satisfy the regularity assumed above see [GT83].
2.4. Assumption (regularity of A). From hereon in we will assume that the problem data is sufficiently smooth such that solutions exist and belong to at least H 2 (Ω) ∩ H 1 0 (Ω). 2.5. Remark (regularity of Ω). Theorem 2.2 specifies that Ω must be a C 1,1 domain. We will be approximating such a domain with one which is only C 0,1 (i.e., a polyhedral one). We thus assume that the model problem admits a unique strong solution even when Ω is only C 0,1 . To circumvent this assumption curved finite elements could be used to fit the boundary exactly [Ber89]. For simplicity we will not present this case here, although we believe our analysis can be extended to this case.

Discretisation
Let T be a conforming, shape regular triangulation of Ω, namely, T is a finite family of sets such that (1) K ∈ T implies K is an open simplex (segment for d = 1, triangle for d = 2, tetrahedron for d = 3), (2) for any K, J ∈ T we have that K ∩ J is a full subsimplex (i.e., it is either ∅, a vertex, an edge, a face, or the whole of K and J) of both K and J and (3) K∈T K = Ω. We use the convention where h : Ω → R denotes the meshsize function of T , i.e., where h K is the diameter of K. We let E be the skeleton (set of common interfaces) of the triangulation T and say e ∈ E if e is on the interior of Ω and e ∈ ∂Ω if e lies on the boundary ∂Ω. The assumptions on the tessellation made here are typical in the finite element analysis. For the presentation of the method and its analysis, some assumption could be relaxed (e.g. the form of the elements or the assumption on a conforming triangulation) but this would lead to an unnecessary increase in the complexity of the presentation. Let P k (T ) denote the space of piecewise polynomials of degree k over the triangulation T ,i.e., and introduce the finite element spaces to be the usual spaces of discontinuous piecewise polynomial functions which are compactly and non compactly supported respectively.
3.1. Remark (generalised Hessian). Assume a function v ∈ H 2 (Ω), let n : ∂Ω → R d be the outward pointing normal of Ω then the Hessian D 2 v of v, satisfies the following identity: If v ∈ H 1 (Ω) (3.5) is still well defined in view of duality, in this case we set where the last term is understood as a pairing between H −1/2 (Ω) and H 1/2 (Ω).
3.2. Definition (broken Sobolev spaces, trace spaces). We introduce the broken Sobolev space We also make use of functions defined in these broken spaces restricted to the skeleton of the triangulation. This requires an appropriate trace space 3.3. Definition (jumps, averages and tensor jumps). We define average, jump and tensor jump operators for arbitrary scalar Note that on the boundary of the domain ∂Ω the jump and average operators are defined as We will often use the following Proposition which we state in full for clarity but whose proof is merely using the identities in Definition 3.3.
3.4. Proposition (elementwise integration). For a generic vector valued function p and scalar valued function φ we have An equivalent tensor formulation of (3.15)-(3.16) is In addition for matrix valued V we have that 3.5. Construction of an appropriate discrete Hessian. We now use the framework set out in [ABCM02] to construct a general notion of discrete Hessian. We first give a definition using a flux formulation: 3.6. Definition (generalised finite element Hessian: flux formulation). Let u ∈ H 2 (T ) andÛ : H 1 (T ) → T (E ∪ ∂Ω) be a linear form andp : d a bilinear form representing approximations to u and ∇u over the skeleton of the triangulation. Then we define the generalized finite element Hessian We now present the primal formulation for the generalized finite element Hessian: 3.7. Theorem (generalised finite element Hessian: primal form). Let u ∈ H 2 (T ) and letÛ andp be defined as in Definition 3.6, then the generalised finite element (3.23) Proof Note that in view of Definition 3.3 for generic vector fields q ∈ W and v ∈ V we have the following identity Then summing (3.21) over K ∈ T and making use of the identity (3.24) we see (3.25) Using the same argument for (3.22) (3.26) Note that, again making use of (3.24) we have for each q ∈ H 1 (T ) d and v ∈ H 1 (T ) that Taking v = u in (3.27) and substituting into (3.22) we see (3.28) Now choosing q = ∇ h Φ and substituting (3.28) into (3.21) we arrive at the fully generalised finite element Hessian given by (3.23).
3.8. Remark (consistent representations of the gradient operator). If one were interested in consistent representations of other derivatives, for example the gradient operator, one would need to modify the proof of Theorem 3.7. Examples of consistent gradient representations can be found in [ABCM02]. See also [BO09,DPE10,BE08]. Using this methodology it should be possible to construct an entire hierarchy of derivatives. 3.9.
Example. An example of a dG formulation for the approximation to the Hessian, D 2 u, can be derived by taking the fluxes in the following waŷ (3.31) 3.10. The discontinuous nonvariational finite element method. We are now in a position to state the numerical method for the approximation of (1.1). We look to find u h ∈ where the penalisation parameter σ > 0 is to be chosen sufficiently large to guarantee coercivity.
Using the L 2 projection operator P V : it is possible to elliminate the finite element Hessian from the bilinear form for sufficiently smooth A: 3.11. Lemma (elimination of the finite element Hessian in a general setting). If ó d×d and the fluxes are chosen as in Example 3.9 then (3.36) Proof. This follows from the following identity (3.37) 3.12. Remark. The solution of the problem in this form is nontrivial due to the global L 2 (Ω) projection appearing in the formulation. However, in the discontinuous setting the global L 2 (Ω) projection is in fact computable locally. We may actually exploit this fact to optimise our schemes efficiency. We will discuss this further in the sequel.
The nonvariational finite element method thus coincides with the classical (symmetric) interior penalty method for the Laplacian [DD76].
3.14. Remark (relation to standard dG methods). It is not difficult to prove that choosing to numerical fluxes in the same way as presented in [ABCM02, Table 3.2] results in the same correlation to the dG methods summarised in the aforementioned paper for the case that A is constant. For brevity we will not prove this here.
Note that when A is not constant we have that the nonvariational finite element method does not coincide with its standard variational finite element counterpart. There is an extra stability property which allows the method to successfully cope with advection dominated problems [LP11, §4.2] which is illustrated by the result of Lemma 3.11.
We conclude this section with a proof consistency of the method and then show that Galerkin orthogonality holds.
3.15. Lemma (consistency). Let u ∈ H 2 (T ) and assume that the numerical fluxes are chosen in a consistent fashion in the sense of [ABCM02, §3.1], that is, Proof Applying Proposition 3.4 to the first term in the definition of 3.17. Remark. If A is piecewise constant then since H[u] = P V D 2 u we have J(Ψ) = 0 and we recover the usual Galerkin orthogonality A h (u h − u, Ψ) = 0. We will show in the next section that in general the error functional J is of higher order for smooth enough u.

Coercivity, continuity and convergence
In this section we examine the coercivity, continuity and convergence of the method. We will focus on the fluxes given in Example 3.9 to simplify the presentation. Furthermore we make the following additional assumption on the problem data.

Assumption (coercivity of the nonvariational problem). For the rest of this section we will assume that the nonvariational operator is coercive, that is
ó d×d and that div(DA) ≤ 0.

4.2.
Remark (variational nature of the coercive problem). Under Assumption 4.1 the problem can be written variationally. The solution to the nonvariational problem is the minimiser to the (degenerate) second order variational problem: Find u ∈ H 2 (Ω) ∩ H 1 0 (Ω) such that

Definition (H 1 (T ), H 2 (T ) and H −1 (T ) norms). We introduce the broken H 1 (T ) and H 2 (T ) norms as
These are equivalent to their continuous equivalent norms for functions in V D .

4.4.
Proposition (projection approximation in V D ). Let P V : L 2 (Ω) → V D be the L 2 (Ω) orthogonal projection operator defined by (3.35). Using standard approximation arguments we have that In particular, let A h denote the L 2 orthogonal projection of A into the space of piecewise constant functions, then we have . 4.5. Theorem (stability of H [Pry12, Theorem 4.10]). Let H be defined as in Example 3.9 then the dG Hessian is stable in the sense that Consequently we have We now state the following technical Lemmata.

Lemma (upper bound on
and v h ∈ V then it holds that (4.10) Proof Let A h denote the L 2 orthogonal projection of A into the space of piecewise constant functions. Then adding and subtracting appropriate terms we see Using the approximation properties of the L 2 (Ω) projection given in Proposition 4.4 we see that (4.12) Now using the properties of A h (4.7) and inverse inequalities we have (4.13) as required.

Lemma (upper bound on some skeletal terms).
Let v h ∈ V D , B ∈ W 1 ∞ (Ω) d×d and B h ∈ V d×d D be the L 2 orthogonal projection of B, then in view of trace and inverse inequalities we have the following: (4.14) Proof For brevity we prove only the first inequality, the second and third follow similar arguments. In view of the definition of the average operator (3.9) it follows that Now by a trace inequality we see that Using the stability of the L 2 (Ω) projection operator in H 1 (Ω) we have (4.17) as required.
4.8. Theorem (discrete continuity and coercivity). Let the conditions in Assumption 4.1 hold. Suppose also that DA is sufficiently small such that where γ is the ellipticity constant, C 1 is a constant appearing in Proposition 4.4 with k = 0, C 2 , C 4 and C 5 are the constants appearing in Lemmata 4.6 and 4.7 and > 0 is some parameter. In addition assume σ the penalisation term is sufficiently large, specifically where A h is some piecewise constant approximation to A. Then there exist positive constants C B and C C such that We postpone the proof of this theorem to the end of this section and first prove the error estimates for our discrete solution.
The results of Theorem 4.8 allow us to invoke Strang's second Lemma.

Corollary (Strang [EG04, c.f.]).
There exists a C > 0 such that 4.10. Lemma (discrete negative norm convergence). Let A ∈ L ∞ (Ω) and u ∈ H k+3 (Ω). Then we have that there exists a constant C > 0 such that Proof We have, in view of Cauchy-Schwarz inequality, that for Ψ ∈ Proof The proof of (4.26) is immediate from applying Proposition 4.4 to Corollary 4.9 with v h = P V u and noting that the bound for the consistency error nothing but the result of Lemma 4.10, concluding the proof.
To conclude this section we prove Theorem 4.8. (4.27) Now making use of the definition of H from Example 3.9 we see (4.28) Adding and subtracting appropriate terms we have that which rewriting variationally gives (4.30) Note that and hence we see that upon substituting (4.31) into (4.30) that (4.32) We proceed by applying Cauchy-Schwartz componentwise to (4.32) and estimating θ by 1 arriving at (4.33) In view of Lemma 4.6 and the Poincaré inequality we have (4.34) For the skeletal terms we apply the result of Lemma 4.7 which upon substituting these into (4.34) we see that as required.
For coercivity we use the equality given in (4.32) with v h = u h to find (4.36) We proceed by bounding each term individually. By the ellipticity of the problem we have that (4.37) By the coercivity of the problem we have By the Cauchy-Schwartz inequality and making use of Lemma 4.6 (4.39) We combine the fourth and fifth terms and let A h denote the L 2 orthogonal projection of A onto the space of piecewise constant functions. Upon adding and subtracting appropriate terms (4.40) Using the identities we have that (4.43) Using Cauchy-Schwartz we see (4.44) Making use of the various bounds from Lemma 4.7 we have (4.45) We now apply a Cauchy inequality and use the approximation properties of A h from Proposition 4.4 to find for any > 0 that (4.46) The other terms are bounded similarly in that and (4.48) Note that the final two terms are already in their desired form since (4.49) Collecting the bounds from (4.46)-(4.50) shows (4.51) The final term in (4.36) is given by (4.52) Finally, collecting the bounds from (4.37), (4.38), (4.39), (4.51) and (4.52) shows (4.53) Coercivity of the discrete bilinear form follows using the assumption in Theorem 4.8, by choosing sufficiently small and the penalisation parameter σ sufficiently large for small enough h.

4.13.
Remark (the coercivity bound). We note that the coercivity bound relies on the term DA not becoming too large, as specified in Theorem 4.8. If it is we view this as an advection dominated problem. Our numerical experiments suggest that there is no condition on the size of this term. If the coefficient matrix A is divergence free, i.e., DA = 0 then the bound simplifies considerably. For example, in the case that A is constant we regain the same theoretical results as for the method given in Example 3.13.

Numerical experiments
In this section we detail numerical experiments carried out in the finite element package Dune-Fem [DKNO10] which is based on the Dune software framework [BBD + 08a, BBD + 08b]. The code used to test the method will be made freely available within the Dune-Fem-Howto in a future release.
We present some benchmark problems designed such that the exact solution is known. In each of the experiments the domain Ω = [0, 1] 2 and we consider the coefficient matrix to be In each of the numerical experiments we make use a stabilised conjugate gradient solver taken from the Dune-Istl module [BB07] preconditioned with an incomplete LU factorisation. We choose the penalty parameter σ = 20. 5.1. Test 1 : a coercive operator. In this test we take the components of A such that the operator is coercive, fitting into the analytical framework presented in §4. With x = (x 1 , x 2 ), we set We choose the problem data such that the exact solution is given by u(x) = sin (πx 1 ) sin (πx 2 ) and approximate this using the formulation (3.32). In Tables 1a-1b we present the results for the cases k = 1, 2, numerically demonstrating that the analytical rates of convergence are achieved in the dG energy norm, moreover, optimal convergence is achieved in L 2 (Ω). We choose the exact solution as in 5.1 and conduct the same tests. Tables 2a-2b detail the results. Note that this is not a coercive operator and as such, does not fit into the analytical framework presented in §4, we do however still achieve optimal convergence in · and · dG,1 .  5.3. Test 3 : irregular solutions. In this test we consider the case the exact solution does not satisfy the regularity requirements presented in the analytical framework of §4, i.e., u ∈ H k+3 (Ω). In addition we consider the case that u ∈ H 2 (Ω), demonstrating the method converges even for viscocity solutions of the problem. We consider the coercive operator from §5.1 and choose the problem data such that Note that this function is H 2 (Ω) but not H 3 (Ω). We also take the problem data such that This function is H 1 (Ω) but not H 2 (Ω). The results are given in Tables 3a-3b.
In the case u is given by (5.7) the scheme converges with optimal rate in the · dG,1 norm even if the solution is not in H 3 . The convergence in the L 2 is more erratic, but we observe the same behavior testing the standard IP FEM taking A to be the identity.
In the case u is given by (5.8) the convergence rates are suboptimal since the solution is not H 2 . Table 3. 5.3 -Test 3. We present errors and convergence rates of the approximation given by solving (3.32). In both cases we consider k = 1.

Conclusions and outlook
In this work we have extended the framework from [LP11] for linear nonvariational problems to incorporate discontinuous approximations.
We have shown the method presented (and subsequently that of the continuous case from [LP11]) is well posed and converges optimally under coercivity assumptions on the coefficient matrix A.
In the numerical experiments we note the the method is well posed and converges optimally even for A which do not satisfy the coercivity assumptions or u which do not satisfy the regularity needed in the analytical framework. This motivates that another analytical approach needs to be developed. This approach can not be variational in nature as such will be completely non standard. This is the topic of ongoing research.