Automatic symbolic computation for discontinuous Galerkin finite element methods

The implementation of discontinuous Galerkin finite element methods (DGFEMs) represents a very challenging computational task, particularly for systems of coupled nonlinear PDEs, including multiphysics problems, whose parameters may consist of power series or functionals of the solution variables. Thereby, the exploitation of symbolic algebra to express a given DGFEM approximation of a PDE problem within a high level language, whose syntax closely resembles the mathematical definition, is an invaluable tool. Indeed, this then facilitates the automatic assembly of the resulting system of (nonlinear) equations, as well as the computation of Fr\'echet derivative(s) of the DGFEM scheme, needed, for example, within a Newton-type solver. However, even exploiting symbolic algebra, the discretisation of coupled systems of PDEs can still be extremely verbose and hard to debug. Thereby, in this article we develop a further layer of abstraction by designing a class structure for the automatic computation of DGFEM formulations. This work has been implemented within the FEniCS package, based on exploiting the Unified Form Language. Numerical examples are presented which highlight the simplicity of implementation of DGFEMs for the numerical approximation of a range of PDE problems.

1. Introduction. The finite element method (FEM) represents an indispensable computational tool for the accurate, efficient, and rigorous numerical approximation of continuum models arising within a wide range of scientific and engineering application areas. Key reasons for the success of FEMs include their applicability to very general classes of partial differential equations (PDEs), simple treatment of complicated computational geometries and enforcement of boundary conditions, ease of adaptivity including both local mesh subdivision (h-refinement) and local polynomial enrichment (p-refinement), and, from a mathematical point of view, the availability of tools for their rigorous error analysis. However, when compared with their finite difference counterparts, FEMs are typically regarded as being complicated to implement. Indeed, assembly of the underlying matrix stemming from the FEM discretisation of a given (linear, for example) PDE problem typically involves mapping each element present in the computational mesh, defined in the global coordinate system, to a given reference or canonical element, where both the local FEM basis and corresponding quadrature, needed to approximate the underlying integral, is defined. In this manner, local elemental stiffness matrices and load vectors may be computed; these entries are then inserted into the global matrix and right-hand side vector, respectively, according to the elementwise local-to-global degree of freedom mapping, subject to the enforcement of inter-element continuity constraints and the imposition of boundary conditions. This general assembly strategy is almost universally employed within both open source and commercial software, e.g., FreeFem++ [23], DUNE [8], deal.II [4] and OpenFOAM [40], to name just a few. Thereby, in principle, assuming that a user can define both the element stiffness matrix and load vector, then a numerical approximation may be readily determined. However, most open source software packages do not provide easy-to-use user interfaces; indeed, typically the user must first understand the low level language in which the code is written, e.g., C++, Fortran, etc, and develop an understanding of the underlying datastructures and function/subroutine calls defined within their chosen package in order to be able to develop code specific to their own application. This can, of course, be a rather time-consuming exercise, and often requires one to use multiple software libraries, for example, when certain features are not available within the chosen package which the user needs to utilise.
In the nonlinear setting, assuming a Newton-type iteration is exploited to solve the underlying system of nonlinear equations stemming from the given FEM employed, the general strategy of assembly of the resulting linearised equations is similar to the linear case; in this setting the load vector represents the residual of the numerical scheme. However, in this case, the Fréchet derivative of the FEM must now be computed; we point out that, in the context of the numerical detection of bifurcation points, the number of derivatives that must be computed just to form the so-called extended system needed to accurately compute the bifurcation point, even before a Newton iteration is implemented, is dependent on the codimension of the singularity being sought, cf. [10,11]. In general, the evaluation and implementation of both the FEM residual and, moreover, the Fréchet derivative(s) of the scheme is a difficult and time consuming task, which is inevitably prone to human error. This is particularly the case when exploiting FEMs for the numerical approximation of systems of coupled nonlinear PDEs, including multiphysics problems, whose parameters may consist of power series or functionals of the solution variables, cf. [28].
Thus far, our discussion has primarily focussed on the application of conforming FEMs, whereby, excluding Neumann, or weakly imposed boundary conditions, only element contributions need to be evaluated. The exploitation of more general FEMs, and in particular discontinuous Galerkin FEMs (DGFEMs) requires the implementation of inter-element flux terms, which involves combinations of both innerand outer-traces of the FEM solution, relative to a given orientation of the element face. In recent years, DGFEMs have become an increasingly popular class of FEMs, most notably due to their local conservation properties, inherent numerical stability for convection-dominated diffusion problems, limited interelement communication, which is restricted only to neighbouring elements, and has important advantages for the implementation of boundary conditions and the parallel efficiency of the method, and finally the ease in which so-called hanging nodes can be treated, and the efficient implementation of hp-adaptivity. Indeed, tremendous progress has been made on both the analytical and computational aspects of DGFEMs; for a review of some of the key developments in the subject, we refer to the recent monographs [13,14,24,35]. For a historical review of DGFEMs, we refer to the articles [3,12], and the references cited therein.
The addition of inter-element face terms within DGFEMs further complicates the implementation of such schemes within standard FEM packages. Assuming for a moment that the underlying PDE problem is linear, we note that for a given interior face, shared by two neighbouring elements, there are, in general, four local face matrices, which stem from the different combinations of traces of the FEM solution from the two elements whose boundaries form the current face. Once these local face-wise matrices have been assembled, they can then be inserted into the global FEM matrix in an analogous manner to the treatment of the element stiffness matrix. On boundary faces, contributions to the load vector must also be computed. Again, in the nonlinear setting, the task of computing both the residual vector and Fréchet derivative(s) of the DGFEM scheme is a very technical and time consuming task.
The purpose of this article is to discuss the use of symbolic algebra to facilitate the assembly of FEM matrix problems, and in particular those arising from the application of DGFEMs, for the numerical approximation of general nonlinear systems of PDEs. The general approach is to develop a high level language syntax, which closely corresponds to the mathematical formulation of the underlying FEM. Thereby, through this layer of abstraction, the user needs to only specify the FEM residual in a concise and easy to read manner, whereby the evaluation of the Fréchet derivative(s) of the scheme are automatically computed symbolically and the resulting low level C++/Fortran code for element stiffness and face matrices and load/residual vectors are automatically generated by the so-called form compiler. In this way, any existing open source FEM package may be utilised, subject to the implementation of a suitable interface which directly calls the automatically generated snippets of code when assembling element and face matrices and load/residual vectors. Most importantly, we stress that once the user has selected a particular FEM package which is appropriate for their purposes, the interface to that software platform which links to the automatically generated code only needs to be written and debugged once; it may then subsequently be exploited to solve a potentially huge variety of PDE problems, with a plethora of FEM schemes. At this point it is pertinent to mention that some FEM packages do indeed include such a symbolic interface at the heart of their design; most notably we mention the excellent FEniCS package, cf. [1,31], for example. The Unified Form Language (UFL) component of FEniCS provides an easy to use python interface which allows for the automatic FEM numerical approximation of systems of PDEs in a user-friendly manner. Other form compilers include SyFi [2], and Manycore Form Compiler [32], for example.
However, even exploiting a package such as UFL, as powerful as it is, the definition of DGFEMs in this framework is still rather verbose, particularly for nonlinear systems of coupled PDEs. With this in mind, we present a further layer of abstraction for the automatic computation of DGFEM formulations employing symbolic algebra. Indeed, in this article DGFEM utility functions for general PDE operators are developed, which significantly simplifies the specification of the DGFEM discretisation of a given problem. This work was originally inspired by the need to numerically model the formation of a hydrogen plasma in a microwave power assisted chemical vapour deposition reactor employed for the manufacture of synthetic diamond; this includes constituent equations for the background gas mass average velocity, gas temperature, electromagnetic field energy and plasma density, cf. [28,36].
In [36] we originally developed easy-to-use DGFEM utility functions as part of our inhouse software package AptoPy [36], for application with our own FEM package AptoFEM [26]. AptoPy is written in Python and exploits the open source symbolic algebra package SymPy [37]. We stress that these choices are entirely user-dependent; indeed, in the past we have employed the symbolic algebra packages REDUCE [22], Mathematica, and Maple. In addition to AptoFEM, ENTWIFE has also been employed as the low level FEM package. The DGFEM utility functions written in Ap-toPy have been ported to UFL for use within the FEniCS package; full open source codes are available from https://bitbucket.org/nate-sime/dolfin _ dg. With this in mind, for simplicity of presentation, throughout this article, we only show snippets of UFL code in order to highlight how the DGFEM utility functions may be exploited in practice; the corresponding AptoPy syntax is quite similar, cf. [36].
The outline of this article is as follows. In Section 2 we briefly outline the DGFEM discretisation of general systems of nonlinear conservation laws. Then, in Section 3 we propose a computational framework for the automatic generation of DGFEM schemes within a simple unified setting. On the basis of this work, in Section 4 we provide some examples to illustrate the flexibility and ease within which systems of PDEs may be numerically approximated using the software developed in this article. Finally, in Section 5 we provide a summary of the work undertaken in this article, as well as outlining potential future developments.
2. Discontinuous Galerkin Finite Element Methods. As a representative PDE example, in this section we outline the DGFEM discretisation for the following system of conservation laws: . . , f v d (u, ∇u)) represent the convective and diffusive fluxes, respectively, which are assumed to be continuously differentiable, and f is a given source function. For simplicity of presentation, we assume that (2.1) may be supplemented with the boundary conditions: where Γ = Γ D ∪ Γ N , and Γ D and Γ N are two disjoint subsets, with Γ D nonempty and relatively open in Γ. We stress that more general boundary conditions can also be considered; for example, in the case when (2.1) represents the compressible Euler or Navier-Stokes equations, we refer to, for example, [17] and [20], respectively. For the purposes of discretisation, we rewrite (2.1) in the following equivalent form: Here, the matrices G kl (u) = ∂f v k (u, ∇u)/∂u x l , k, l = 1, . . . , d, (2.4) are the homogeneity tensors defined by f v k (u, ∇u) = G kl (u)∂u/∂x l , k = 1, . . . , d. We write the homogeneity tensor product For simplicity of presentation, we now proceed to discretise (2.1)-(2.3) based on employing the symmetric interior penalty (SIPG) formulation presented in [20]; however, we stress that other DGFEMs can easily be included within this general setting, cf. below. To this end, we partition Ω into a mesh T h = {κ} consisting of non-overlapping open element domains κ, such thatΩ = ∪ κ∈T hκ . For each κ ∈ T h , we denote by n κ the unit outward normal vector to the boundary ∂κ. We assume that each κ ∈ T h is an image of a fixed reference elementκ, that is, κ = σ κ (κ) for all κ ∈ T h . On the reference elementκ we define spaces of polynomials of degree ≥ 0 as follows: Thereby, with this notation, we define the following DGFEM finite element space An interior face of T h is defined as the (non-empty) (d−1)-dimensional interior of ∂κ + ∩ ∂κ − , where κ + and κ − are two adjacent elements of T h , not necessarily matching. A boundary face f of T h is defined as a (non-empty) (d − 1)-dimensional facet of κ, κ ∈ T h , where κ is a boundary element of T h , such that f ⊂ ∂κ ∩ Γ. We denote by Γ I the union of all interior faces of T h . Let κ + and κ − be two adjacent elements of T h , and x an arbitrary point on the interior face f = ∂κ + ∩ ∂κ − . Furthermore, let v and τ be vector-and matrix-valued functions, respectively, that are smooth inside each element κ ± . By (v ± , τ ± ), we denote the traces of (v, τ ) on f taken from within the interior of κ ± , respectively. Then, the averages of v and τ at x ∈ f are given by where n denotes the unit outward normal vector to Γ.
The interior penalty DGFEM discretisation of (2.1)-(2.3) is given by: The subscript h on the operator ∇ h is used to denote the discrete counterpart of ∇, defined elementwise. Here, H(·, ·, ·) denotes the (convective) numerical flux function; this may be chosen to be any two-point monotone Lipschitz function which is both consistent and conservative; typical choices include the (local) Lax-Friedrichs flux, the HLLE flux, the Roe flux, and the Vijayasundaram flux, cf. [30,39]. The penalisation term δ(·) arising in the DGFEM scheme (2.6) is given by where C IP is a (sufficiently large) positive constant, cf. [15]. Finally, we define the boundary terms present in the form N Γ (·, ·) by Here, the viscous boundary flux F v Γ and the corresponding homogeneity tensor G Γ are defined by The convective boundary flux H Γ is defined by Finally, the boundary function u Γ (u) is given according to the type of boundary condition imposed; in the current setting u Γ (u) = g D on Γ D and u Γ (u) = u on Γ N . For further details regarding the imposition of more general boundary conditions, we refer to [20], for example.

Computational Framework for DGFEMs.
In this section we present the general computational framework for the automatic generation of DGFEM (semi-) linear forms in a concise and easy to use manner. We begin by outlining the treatment of both convective and viscous components arising in conservation laws. In this setting, the DGFEM formulations can be constructed in a consistent manner. We exploit this by designing utility functions to automatically generate these symbolic DGFEM formulations. We then proceed by proposing a hierarchical framework for computing DGFEM formulations of PDE operators. This hierarchy takes advantage of the DGFEM utility functions, providing a modular framework for a suite of DGFEM formulations for various operators arising in PDE problems of engineering interest.

Automatic Treatment of Convective Terms.
In this section we discuss the automatic symbolic representation of the DGFEM discretisation of the term involving the convective flux function F c (·). To this end, we recall that the convective numerical flux function H(·, ·, ·), cf. (2.6), must be defined on the element interfaces. Defining a callable function F _ c which specifies F c (·) we construct the abstraction of the numerical flux function as shown in Listing 1. The methods interior and exterior will be called to construct the flux function on interior and exterior faces, respectively. The method setup is provided for the inheriting class to initialise any members prior to calls to interior and exterior. This design allows for the flux function to be different on the two types of faces present in the mesh.
For simplicity, here we consider two implementations of the ConvectiveFlux class, based on employing the local Lax-Friedrichs and HLLE fluxes, though we stress that other fluxes, such as the Vijayasundaram or Roe flux, for example, may also be employed. Thereby, on the boundary ∂κ of an element κ ∈ T h , we define the local-

Lax Friedrichs flux H LF and HLLE flux H HLLE by
respectively. Here, α is the local dissipation parameter, which is selected to be the maximum of the (absolute value) of the eigenvalues of the Jacobi matrix evaluated on the element face. More precisely, we have that , where λ max denotes the largest eigenvalue (in modulus) of B (·, n κ ). Additionally, where λ min denotes the smallest eigenvalue (in modulus) of B (·, n κ ). Given F _ c, which specifies F c (·), the numerical flux functions H LF (·, ·, ·) and H HLLE (·, ·, ·) can be automatically generated in order to yield the discretisation of the convective term present in the underlying PDE problem; we note that the constructors of both classes LocalLaxFriedrichs and HLLE require the symbolic representation of the eigenvalues of B. As an example, we consider the application of the DGFEM employing the local Lax-Friedrichs flux for the numerical approximation of the linear advection equation bu, and f is some given forcing function. Here, the dissipation parameter can be shown to be α| ∂κ = |b · n κ |. The UFL code required to generate the numerical flux function, together with the corresponding element boundary term arising in the DGFEM scheme for this problem is given in Listing 2. We note that the use of the class HLLE defining the HLLE numerical flux follows in an analogous manner. In many cases an analytical expression for the eigenvalues of the Jacobi matrix B(u, n κ ), κ ∈ T h , may be difficult to evaluate; thereby, packages such as SymPy [37] may be used to compute them symbolically. To this end, B(u, n κ ) can be assembled using symbolic differentiation; the eigenvalues of this matrix may then be computed symbolically by exploiting the Berkowitz algorithm [7]. In Listing 3 we show an example of this method applied to the convective component of the incompressible Navier-Stokes equations, i.e., F c (u) = u ⊗ u.

Automatic Treatment of Viscous Terms.
In this section we develop utility functions which automatically generate the DGFEM discretisation of secondorder PDE operators. To this end, we recall from Section 2 that the viscous component of the underlying PDE is given by , dg _ outer(v, n)) * dExt \ + inner(self.sig * hyper _ tensor _ product(G, dg _ outer(u -u _ gamma, n)), dg _ outer(v, n)) * dExt return residual Listing 4: The class DGFemViscousTerm provides the abstract interface for interior and exterior residual formulation using symbolic algebra. The class DGFemSIPG uses UFL to automatically formulate the interior and exterior terms of the SIPG formulation of the viscous component of (2.6).
here, we have set the right-hand in (3.3) to zero, for simplicity, in order to concentrate on the discretisation of the viscous term. The semilinear DGFEM formulation of (3.3), shown in equation (2.6), is encapsulated by implementations of the abstract class is a callable function, which defines the viscous flux function, u and v denote the trial and test functions, respectively, gamma ≡ γ = C IP 2 /h is the DGFEM penalisation coefficient, G ≡ G(u) is the homogeneity tensor, and n ≡ n is the face normal.
Recalling the definition of the homogeneity tensor (2.4) and the homogeneity tensor product (2.5), G(u) is automatically computed using the function homogeneity _ tensor (F _ v, u) and the function hyper _ tensor _ product(G, tau) computes the homogeneity Listing 5: Example UFL code for the DGFEM discretisation of the quasi-linear PDE (3.5) using the DGFemSIPG utility class.
tensor product. These functions may then be employed for the generation of the corresponding DGFEM formulation. On the basis of these two functions, we introduce the abstract class DGFemViscousTerm; this offers the following three methods for handling each of the boundary components present in the DGFEM scheme: 1. DGFemViscousTerm.interior _ residual(dS) automatically generates terms associated with the interior boundaries Γ I present in (2.6); 2. DGFemViscousTerm.exterior _ residual(u _ gamma, ds _ i) generates the terms associated with exterior boundary component ds _ i present in (2.7) with boundary condition u Γ (u) = u _ gamma; 3. Finally, DGFemViscousTerm.neumann _ residual(gN, ds _ i) generates any terms arising from Neumann boundary conditions with flux specification F v (u, ∇u)· n = g N . We demonstrate an example of the implementation of the DGFemViscousTerm with the SIPG method in Listing 4; for the sake of brevity, we have removed input and error checking. Further implementations such as the non-symmetric interior penalty (NIPG) and Baumann-Oden schemes are available in the classes DGFemNIPG and DGFemBO, respectively. The DGFEM formulation proposed by Bassi & Rebay [6] is challenging in the symbolic framework due to the requirement of solving elementwise problems for the local lifting operator. Such operations are not easily formulated in the UFL for use with DOLFIN, for example; the implementation of DGFEMs defined based on employing local lifting operators will be considered as part of our programme of future research. For a more detailed discussion of the formulation of lifting operators in DOLFIN, we refer to [34,Chapter 5].
An example of the UFL code required to generate the DGFEM semilinear residual discretisation of the quasi-linear second-order PDE problem is shown in Listing 5.

Automatic Generation of DGFEM Formulations.
Even with the utility functions outlined in Sections 3.1 and 3.2, the specification of the DGFEM scheme can be very verbose; this is particularly the case when discretising systems comprising many PDE variables with multiple boundary conditions. In this section we propose a hierarchical scheme for the management of DGFEM formulations of PDE operators bcs = [DGDirichletBC(bc _ 1, ds _ 1), DGDirichletBC(bc _ 2, ds _ 2), ...] vt = DGFemViscousTerm(F _ v, u, v, delta) ext = sum(vt.exterior _ residual(bc.get _ function(), bc.get _ boundary()) for bc in bcs) Listing 6: Example of the automatic generation of the exterior residual terms of a given DGFemViscousTerm for a list of Dirichlet boundary conditions. of increasing complexity.
3.3.1. Boundary Conditions. Firstly, we outline a simple framework for managing the boundary conditions enforced within a given PDE problem. We define the DGBC (discontinuous Galerkin boundary condition) abstract class from which the classes DGDirichletBC and DGNeumannBC inherit. These implementations simply serve to store the boundary condition and the boundary component over which the condition should be enforced. For example, applying a Dirichlet boundary condition as required by the quasi-linear PDE problem stated in (3.5) simply requires instantiation of DGDirichletBC(ds _ D, g _ D). Similarly, for the imposition of the Neumann boundary condition we construct DGNeumannBC(ds _ N, g _ N); note that Robin conditions can be implemented in an analogous manner. By generating a list of boundary conditions in this manner, we may easily formulate their imposition within the DGFEM scheme. As an example of a series of Dirichlet boundary conditions being automatically generated using DGFemViscousTerm, we refer to Listing 6.

Abstract DGFEM Formulation.
We encapsulate the abstraction of a DGFEM scheme in the class DGFemFormulation which prescribes one abstract method generate _ fem _ formulation. The DGFemFormulation constructor requires the mesh T h , the function space V m and the list of boundary conditions. This class will serve as the base class for the DGFEM formulation of all derived PDE operators. In this work we describe two direct children of DGFemFormulation: HyperbolicOperator and EllipticOperator. We stress that the flexibility of this design permits DGFEM formulations of other PDE operators, cf. Sections 4.6, 4.7, and 4.8 below.
The class HyperbolicOperator inherits DGFemFormulation and its purpose is to generate DGFEM formulations of the PDE operator ∇ · F c (·). The class implementation is shown in Listing 7; here, generate _ fem _ formulation is overridden to construct the DGFEM formulation of the provided convective flux operator on the interior and exterior faces, as well as on the elements in the mesh. The numerical flux function provided in the constructor H, which implements ConvectiveFlux, is used for the DGFEM formulation. To highlight to modularity of this design, consider an extension of the HyperbolicOperator for the time-dependent Burgers' equation Setting t = y, we may recast (3.6) in the following equivalent form the implementation of the class SpacetimeBurgersOperator (using LocalLaxFriedrichs by default) is depicted in Listing 8. Here, we note that the derived class must simply specify the form of the convective flux F c (u) and the flux function H(·, ·, ·); indeed, once these are defined, the automatic generation of the DGFEM formulation is then handled by the parent DGFemFormulation. Secondly, we introduce the class EllipticOperator, which inherits DGFemFormulation, and requires the specification of F v (·, ∇·) at instantiation. The overridden method generate _ fem _ formulation is then written to automatically generate the interior and exterior boundary, as well as the element, integration terms, implementing all of the concepts of the utility functions for elliptic operators in DGFemViscousTerm. The outline of the EllipticOperator class, which generates the SIPG formulation by default, is given in Listing 9. To highlight the modularity of EllipticOperator in

Hierarchy of PDE Operators.
A natural extension of this framework is a hierarchy of automatically generated DGFEM operators. As shown in Listing 8, the SpacetimeBurgersOperator inherits from HyperbolicOperator. Consider now, for example, the compressible Navier Stokes equations. Here the inheritance chain may begin with the HyperbolicOperator, from which a CompressibleEulerOperator would inherit, and in turn a CompressibleNavierStokesOperator would inherit Compress-ibleEulerOperator and additionally EllipticOperator for the viscosity terms. Further implementations of each member of this class hierarchy need not only be undertaken by inheritance. Operators deriving from models of large physical systems may more appropriately aggregate sub-operators as necessary. The derived classes must simply manage the function spaces and boundary conditions amongst the aggregated DGFemFormulation members. This framework significantly reduces the code required for subsequent development of DGFEM formulations of PDE operators of increasing complexity. By ensuring that each layer of the hierarchy is correctly verified and fully tested means that DGFEM formulations may be debugged in a very straightforward manner.

Examples.
In this section we present a series of examples of increasing complexity to highlight the ease in which each PDE problem may be discretised using a DGFEM formulation, based on employing the class hierarchy proposed and implemented within this article. We stress that the verbosity and complexity of specifying each DGFEM discretisation is vastly reduced within this modular framework. To this end, we consider the discretisation of a simple scalar advection-diffusion equation, the compressible Euler equations, the compressible Navier-Stokes equations posed in both conserved and entropy variables, the indefinite Maxwell problem, and a hyperelasticity problem. In each case, for simplicity of presentation, we employ the SIPG DGFEM discretisation of the second-order PDE operator and a local Lax-Friedrichs flux for the numerical approximation of the convective terms. For brevity, in some of the examples given below only code snippets will be shown, however, full versions of the corresponding python scripts are available from https://bitbucket.org/nate-sime/dolfin _ dg. the problem: find u such that The specification of the numerical fluxes F c and F v , when K = 1 + u, along with α(u, n) used with the local Lax-Friedrichs flux function within UFL syntax is provided in the code listing presented in Listing 12. Exploiting the software concepts outlined in Section 3 for the automatic computation of DGFEM formulations, the UFL code to construct and solve the resulting DGFEM approximation of (4.1), (4.2) is presented in Listing 13, where we have specified that b = (1, 1) , f = −4e 2(x−y) −2e x−y , and g D = e x−y ; thereby, the analytical solution to (4.1), (4.2) is given by u = e x−y . We note that on the final line of the code listing given in Listing 13, the call to DOLFIN's solve function automatically computes the Gâteaux derivative of the DGFEM semilinear form N (·, ·), cf. (2.6), employed for the numerical approximation of (4.1), (4.2) and utilises a Newton solver to evaluate the DGFEM solution; for further details, we refer to [2]. Finally, in Figure 4.1 we show the asymptotic behaviour of the underlying DGFEM on a sequence of uniformly refined triangular meshes with polynomial orders = 1, 2, 3, 4; as we expect, we observe that the L 2 (Ω) and H 1 (Ω) norms of the error tend to zero at the respective rates of O(h +1 ) and O(h ) as the mesh size h tends to zero. The complete code employed for this numerical example is provided in the file advection _ diffusion.py.

Example 2: Compressible Euler
where the convective flux is given by Here, ρ, u v , p, and E denote the density, velocity vector, pressure, and specific total energy, respectively. The equation of state of an ideal gas is given by where γ = c p /c v is the ratio of specific heat capacities at constant pressure, c p , and constant volume, c v ; for dry air, γ = 1.4. Finally, I denotes the d × d identity matrix. The automatic construction of the DGFEM formulation for the numerical approximation of (4.3) is encapsulated within the class CompressibleEulerOperator. We note that this class simply inherits the HyperbolicOperator class, with the specification of the Euler flux and dissipation parameter required for the definition of the Lax Friedrichs flux, cf. Listing 14 for the case when d = 2. We note that in this setting λ max = max{|u v |, |u v | ± c}, where c = γp/ρ is the speed of sound.
As an illustration of the exploitation of CompressibleEulerOperator, we consider the DGFEM approximation of Ringleb's flow problem for which an analytical solution may be obtained using the hodograph method, cf. [9]. The DGFEM discretisation of this test case has also been studied in [18]. This problem represents a transonic flow in a curved channel domain, where the flow is mainly subsonic, with a small supersonic  must be treated in a careful manner to ensure optimal convergence of the underlying DGFEM discretisation. Indeed, our computational experience suggests that a curved polynomial description of the boundary of order q ≥ + 1 is necessary to ensure that the L 2 (Ω) and H 1 (Ω) norms of the error tend to zero at the optimal rates of O(h +1 ) and O(h ) as the mesh size h tends to zero, for fixed , cf. Figure 4.3. The complete python code for this numerical example is provided ringleb _ example.py; a snippet of this code is depicted in Listing 15 in order to highlight the simplicity of the specification of the DGFEM for this example.

Example 3a: Compressible Navier-Stokes Equations.
In this example we consider the DGFEM discretisation of the compressible Navier-Stokes equations: find u = (ρ, ρu v , ρE) : where F c is defined as in (4.4) and the viscous flux is given by where T denotes the temperature and K is the thermal conductivity coefficient. The stress tensor is defined by where µ is the dynamic viscosity coefficient. Finally, we note that KT = µγ Pr E − 1 2 u 2 v , where Pr = 0.72 is the Prandtl number.
Given the UFL code in Listing 14, together with the specification of the viscous flux in Listing 16, in the case when d = 2, we have implemented the class CompressibleNavierStokesOperator class, which inherits both the CompressibleEuler -Operator and EllipticOperator classes. Indeed, the CompressibleEulerOperator component is treated in an identical manner as in the previous section, while the EllipticOperator component is employed with the UFL specification of the viscous flux. On the basis of these classes, the homogeneity tensor and the resulting symbolic DGFEM formulation can then be automatically generated. Moreover, the Gâteaux derivative of the DGFEM formulation is automatically computed by UFL for use within the Newton solver managed in DOLFIN by invoking the call to solve. As a simple test, we consider the example outlined in [20]; namely, we set Ω = (0, π) 2 and select f so that the analytical solution to (4.3) is given by  sin(2(x + y)) + 4 1 /5 sin(2(x + y)) + 4 1 /5 sin(2(x + y)) + 4 (sin(2(x + y)) + 4) where u v = (u v,1 , u v,2 ) . Furthermore, we set µ = 1, Γ = Γ D and Γ N = ∅. The snippet of UFL code required to solve this problem is given in Listing 17; the complete code is provided in compressible _ navierstokes _ square.py. The orders of convergence of u−u h L2(Ω) and u−u h H 1 (Ω) are reported in Figure 4.4 as the mesh is uniformly refined for = 1, 2, 3, 4; as in [20] we observe that u − u h L2(Ω) = O(h +1 ) and u − u h H 1 (Ω) = O(h ) as h tends to zero, for each fixed polynomial order .

Example 3b: Compressible Flow Around a NACA0012 Airfoil.
To demonstrate the application of the automatic formulation of DGFEMs applied to exterior flow problems, in this section we consider laminar flow around a NACA0012 airfoil whose geometry is defined by where the thickness fraction t = 0.12. Here, we set the angle of attack α = 2 • , Reynolds number Re = 5 × 10 3 , inlet flow Mach number M = 0.5, and impose an adiabatic no slip condition on the airfoil; cf. [20,19]. We note that the adiabatic no slip condition imposes a zero heat flux condition K∇T · n = 0 on the airfoil surface boundary. We have implemented this by defining a new boundary condition in the class DGAdiabaticWallBC inheriting DGBC. This boundary condition is then recognised by the CompressibleNavierStokesOperator implementation which automatically generates the DGFEM formulation accordingly.
In this example, we shall undertake goal-oriented adaptive mesh refinement based on employing a dual weighted residual (DWR) a posteriori error estimator; for details, we refer to [19]. In particular, we select the quantity of interest to be the drag coefficient where Γ airfoil denotes the surface of the airfoil, ψ drag = (cos(α), sin(α)) , C ∞ = code to generate these results is available in the file dg _ naca0012 _ 2d.py.

Example 3c: Entropy Variable
Formulation of the Compressible Navier-Stokes Equations. In this section, we consider the DGFEM approximation of the compressible Navier-Stokes equations written in terms of so-called (symmetrisation) entropy variables; for further details, we refer, for example, to [5,29]. For simplicity of presentation we only consider the two-dimensional case, though the extension to d = 3 follows in an analogous manner. Thereby, we introduce the change of variable u → V(u) given by and its inverse Here, e is the internal energy density and s is a nondimensional entropy. The implementation of this mapping and its inverse in the UFL symbolic algebra framework is depicted in Listing 18. On the basis of this transformation, the symmetrised formulation of the compressible Navier-Stokes equations is given by . On the basis of the compressible Euler and compressible Navier-Stokes examples presented in the previous two sections, the DGFEM formulation of problem (4.12) is straightforward in the symbolic algebra framework of the UFL. We use the convective and viscous fluxes shown in Listings 14 and 16 and implement the transformations (4.10) and (4.11) noting that ∇u(V) = ∂u(V) ∂V ∇V. Indeed, the UFL formulation of the convective and viscous fluxes, F c V and F v V , respectively, are presented in Listing 19. These constructs are then exploited within the DGFEM utility functions in the same manner as in the examples presented in Listings 15 and 17. As a simple test, we consider the numerical approximation of the compressible Navier-Stokes example presented in Section 4.3, based on employing the above entropy variable formulation. To this end, the underlying test problem is first transformed according to the change of variable u → V(u), cf. (4.10), whereby the DGFEM is computed; we then transform the numerical solution according to the inverse mapping (4.11) in order to evaluate the error in the underlying approximation. As in the previous example, we again ob-   4.6. Example 4: Maxwell Operator. In this penultimate example we demonstrate the extensibility of the symbolic DGFEM framework developed in this article by considering the automatic discretisation of the Maxwell problem: find u such that in Ω, (4.13) n × u = n × g D on Γ, (4.14) where, for simplicity, we set F m (u, ∇ × u) = ∇×u and k > 0 is the wave number; for solvability, we assume that k 2 is not a Maxwell eigenvalue. The class MaxwellOperator implements the symmetric interior penalty discretisation of (4.13), (4.14) outlined in [27]. Thereby, given the ufl representation of F m (u, ∇ × u) depicted in Listing 21, the corresponding code needed to compute the DGFEM approximation of (4.13), (4.14) can be written in the very compact form shown in Listing 22; see maxwell.py. As a simple test case, here we have set Ω = (−1, 1) 2 , k = 2, g D = (sin(kx), sin(ky)) and f = 0 such that the analytical solution is given by u = g D . Numerical experiments demonstrating the performance of the resulting scheme on a sequence of uniformly refined triangular meshes are presented in Figure 4.7.
4.7. Example 5a: Hyperelasticity. Our final two examples highlight the flexibility of the DGFEM framework proposed in this article for the discretisation of hyperelasticity problems; in this setting DGFEM schemes offer computational benefits in the nearly-incompressible regime, cf. [16,38]. Given a domain Ω 0 defining an elastic body's reference configuration with boundary ∂Ω 0 and outward pointing unit normal vector N, we seek the displacement vector at all points in the domain u(X) : Ω 0 → R 3 which determines the mapping from the reference X ∈ Ω 0 to the deformed x ∈ Ω configuration, such that x = X + u. The constitutive model demands mass balance, so we seek u such that where (∇ X ) i := ∂/∂X i is the gradient in the reference domain, P = ∂Ψ(F)/∂F is the first Piola-Kirchoff stress tensor, Ψ(·) is the strain-energy function, B is a body force, F = I + ∇ X u is the strain tensor, T is a traction force, and ∂Ω 0,D and ∂Ω 0,N denote the Dirichlet and Neumann boundaries, respectively. Here, we examine a neo-Hookean hyperelasticity model where the strain energy density function is given by where µ = E/ (2 (1 + ν)) and λ = Eν/ ((1 + ν) (1 − 2ν)) are the first and second Lamé parameters, respectively, given in terms of Young's modulus E and Poisson ratio ν. Clearly we observe that equation (4.15) can be rewritten with F v (u, ∇ X u) ≡ P, cf. Listing 23; thereby, the DGFEM framework outlined previously may be directly applied. We highlight that the flexibility of this approach naturally allows us to consider other potential strain energy models, such as Saint Venant-Kirchhoff, Ogden, and Mooney-Rivlin [25].
Convergence of the relative errors of the DGFEM approximation of this problem are shown in Figure 4.8. The code to generate these results is available in the file dg _ cantilever.py.
4.8. Example 5b: Near-Incompressible Hyperelastic Buckling. To demonstrate the capability of the DGFEM framework applied to larger three-dimensional problems we draw inspiration from a numerical example of an elastic body buckled in a compressive state, cf. Section 6.7 in [38]. Thereby, we set Ω 0 = (0, 1 /2) 2 × (0, 2), E = 1 × 10 8 , ν = 0.46, B = T = 0 and prescribe u| X3=0 = 0 and u| X3=2 = (0, 0.07, −0.5) . Using the proposed DGFEM framework the resulting deformed mesh configuration is shown in Figure 4.9. Here, we have employed both uniform mesh refinement, as well as adaptive refinement of the computational mesh based on employing a DWR a posteriori error indicator, where we select the internal strain energy Ω0 Ψ(u) dX to be the quantity of interest. As expected, the adaptive refinement algorithm selects regions to enrich the computational mesh near the edges of the Dirichlet boundary, where large changes in the stress occur, as well as regions under compression, which lead to large changes in the internal strain energy. The implementation of this example is provided in the file dg _ compression.py.

Concluding remarks.
In this article we have exploited the use of symbolic algebra for the automatic computation of DGFEMs for the numerical approximation of general systems of nonlinear PDEs. In particular, we have proposed and implemented a class structure in order to allow for the specification of a given DGFEM in a clear and concise manner. While the examples we have presented have primarily focused on flow problems, we stress that the generality of this approach allows for the treatment of a wide range of PDEs stemming from diverse application areas. Indeed, in the final two examples, we have included the DGFEM approximation of the indefinite Maxwell problem and a hyperelasticity problem, respectively. The exploitation of the software developed within this article allows the user to build up the DGFEM Degrees of freedom  discretisation of systems of multi-physics PDEs in a simple and concise manner; as an example, we have employed the software here for the DGFEM approximation of microwave power assisted chemical vapour deposition reactor employed for the manufacture of synthetic diamond; see [28,36]. Moreover, the code can easily be extended to other problems, and indeed other DGFEM schemes, so that users can tailor the software to their own applications.