HDGlab: An Open-Source Implementation of the Hybridisable Discontinuous Galerkin Method in MATLAB

This paper presents HDGlab, an open source MATLAB implementation of the hybridisable discontinuous Galerkin (HDG) method. The main goal is to provide a detailed description of both the HDG method for elliptic problems and its implementation available in HDGlab. Ultimately, this is expected to make this relatively new advanced discretisation method more accessible to the computational engineering community. HDGlab presents some features not available in other implementations of the HDG method that can be found in the free domain. First, it implements high-order polynomial shape functions up to degree nine, with both equally-spaced and Fekete nodal distributions. Second, it supports curved isoparametric simplicial elements in two and three dimensions. Third, it supports non-uniform degree polynomial approximations and it provides a flexible structure to devise degree adaptivity strategies. Finally, an interface with the open-source high-order mesh generator Gmsh is provided to facilitate its application to practical engineering problems.


Introduction
In recent years, hybrid discretisation methods have received increasing attention by the applied mathematics and computational engineering community. The main interest in these methodologies is due to their reduced computational cost with respect to classical discontinuous Galerkin (DG) methods, see [135,163,170,267], from which they inherit appealing stability and convergence properties as well as the flexibility to devise highorder, non-uniform degree and adaptive discretisations and the capability to efficiently exploit parallel computing architectures [42,91,108,158,226].
The purpose of the present contribution is two-fold: to present a review on the state-ofthe-art of hybrid discretisation methods including both fundamental and applied contributions; to provide an educational implementation of the hybridisable discontinuous Galerkin (HDG) method in MATLAB, the so-called HDGlab library, and describe its structure, capabilities and functioning.
HDGlab is an open-source library released under GNU GPL licence and designed for rapid prototyping and testing. It supports simplicial meshes and it provides a seamless 2D and 3D implementation with vectorised loops on the integration points. In addition, HDGlab presents four specific features, currently not available in existing open-source HDG implementations in MATLAB: 1. Availability of high-order polynomial shape functions up to degree 9, with both equally-spaced and Fekete nodal distributions.
3. Support of non-uniform degree polynomial approximations and flexibility to devise degree adaptivity strategy.

4.
Interface with the open-source high-order mesh generator Gmsh.
The remainder of this paper is organised as follows. First, a review of the state-ofthe-art on hybrid discretisation methods is presented in section 2. The formulation of the HDG method for the Poisson and Stokes problems is briefly recalled in sections 3 and 4, respectively. Section 5 provides a description of the structure of the HDGlab library and the url of the repository available under GNU GPL licence. The data structures for the storage of the mesh information, the reference element and the reference face are presented in section 6. Section 7 is devoted to the preprocessing operations, whereas the core of the HDGlab solver for the scalar Poisson equation is described in section 8. Its extension to vectorial problems involving incompressible Stokes flows is discussed in section 9. The visualisation library is introduced in section 10. Section 11 is devoted to numerical examples, in 2D and 3D, validating the optimal convergence properties of the HDG method and showing the potentialities of the HDGlab implementation. Finally, section 12 summarises the capabilities of the presented library and three appendices provide implementation details for the Poisson (Appendix A) and Stokes (Appendix B) solvers and for the interface with the mesh generator Gmsh (Appendix C).

Literature review
The common idea of all hybrid discretisation methods stems from the seminal works of Guyan on static condensation of primal formulations [157] and of Fraeijs de Veubeke on hybridisation of mixed formulations [138] of the finite element method. In the context of element-by-element discontinuous approximations, these techniques allow to remedy the drawback of node duplication in DG methods by considering only the unknowns on the mesh faces (edges in 2D) as globally-coupled degrees of freedom. More precisely, the unknowns in each element are expressed as a function of the degrees of freedom on the element faces by solving a local boundary value problem with purely Dirichlet data, whereas appropriate transmission conditions are imposed to guarantee the interelement continuity of the solution and the fluxes, see [75].
Three families of hybrid numerical schemes lay within this description, namely, (i) hybrid/hybridised DG, (ii) hybridisable DG, henceforth referred to as HDG, and (iii) hybrid high order (HHO) methods. Stemming from classical DG primal formulations, the hybrid or hybridised DG method reduces the number of globally coupled degrees of freedom by performing static condensation [124][125][126]. In addition, improved efficiency can be achieved using polynomial spaces of degree p+1 and p for the primal and hybrid variables, respectively and resorting to the reduced stabilisation approach [206,207]. The hybridisable DG method, henceforth named HDG, is derived from the mixed formulation of the local DG method [77,87,99] with hybridisation. The main advantage of HDG with respect to other hybridised DG methods relies in the introduction of a mixed variable approximating the gradient of the primal unknown [88,89]. This approach is of special interest in the context of engineering problems where quantities of interest often depend on the flux of the solution or on the stress. Finally, HHO bridges the two approaches above by utilising a primal formulation and introducing a local reconstruction operator for the gradient of the solution and an appropriate stabilisation term in the static condensation problem [109,110]. It is worth noting that many hybrid discretisation schemes can be interpreted in a unique framework as HDG-type methods via appropriate definitions of the stabilisation term, see e.g. [68,69] for the staggered DG method and [106] for HHO.
Unified presentations of hybrid discretisation techniques and their relationship with other known numerical methods are available in [37,89,114]. Interested readers are also referred to the review papers [75,149] and to the recent monograph [112]. In the following subsections, an overview of the contributions on hybrid discretisation methods according to the authors' vision is presented.

From linear to nonlinear scalar equations
Second-order scalar elliptic problems have been extensively studied using HDG [89], HHO [110,115] and the hybridised DG method [206], whereas their extension to linear convectiondiffusion problems is discussed in [79,113,124,200]. Cases of higher-order partial differential equations (PDEs) are presented in [78] and [62] for HDG discretisations of biharmonic and third-order equations, respectively, whereas an HHO approximation of the Cahn-Hilliard equation is proposed in [53]. In addition, time-fractional diffusion problems are discussed in [76,197].
More recently, there has been growing interest towards the analysis and simulation of quasilinear and semilinear problems, including the quasilinear p-Laplace operator [95,216] and the semilinear Grad-Shafranov equation [233,234]. To reduce the computational cost of semilinear problems, the interpolatory HDG method was recently devised introducing an interpolation procedure for the efficient and accurate approximation of nonlinear terms [54,100].
Concerning nonlinear problems, HDG discretisations were proposed for nonlinear convection diffusion [201] and nonlinear Schrödinger [46] equations, whereas an HHO formulation of the nonlinear Leray-Lions equation is presented in [111]. Recent applications involving HDG approximations of nonlinear scalar equations focus on the optoelectronic simulation of photovoltaic solar cells. This problem couples a high-order HDG method for the drift-diffusion electronic model in the semiconductor layer of the solar cells with an efficient approximation of the time-harmonic Maxwell's equations [21,56].

Incompressible flows
In the context of incompressible flows, HDG formulations of the Stokes equations were devised and analysed in [90,92,98,199]. The corresponding analysis of the HDG method for Oseen flow is presented in [49]. In [92], it was observed that the HDG method based on Cauchy stress tensor formulation experiences suboptimal convergence of the mixed variable and loss of superconvergence of the postprocessed velocity, when low-order polynomial approximations are considered. The M -decomposition approach [84] remedies this issue by appropriately enriching the discrete local spaces of approximation. An alternative strategy imposing the symmetry of the mixed variable pointwise via Voigt notation is discussed in [148] along with a postprocessing procedure to handle translational and rotational rigid body modes. Divergence-conforming HDG [94], hybridised DG [125] and embedded-hybridised DG (EHDG) [225] discretisations were also studied for the incompressible Stokes equations, whereas a pressure-robust HHO method for viscosity-dependent Stokes flows is proposed in [116].
HDG formulations for the nonlinear incompressible Navier-Stokes equations using equal order and different order of polynomial approximations for the primal, mixed and hybrid variables are described in [50,204] and [215], respectively. The former approach is also employed in [152] to devise a degree adaptive strategy relying on the local superconvergence of the postprocessed velocity. Stemming from the work in [117], different HHO formulations of the incompressible Navier-Stokes equations were proposed, incorporating a skew-symmetric form of the convection term [31] and a globally divergence-free velocity approximation to achieve robustness in presence of large irrotational body forces [45].

Compressible flows and gas kinetics equations
Hybrid formulations for inviscid Euler and laminar compressible Navier-Stokes equations are proposed in [209] in the context of HDG and in [205] for the embedded DG (EDG) method. Extension to viscous turbulent compressible flows using RANS equations with Spalart-Allmaras turbulence model is presented in [193], whereas a large-eddy simulation framework is introduced in [133]. In addition, an entropy stable space-time discretisation was proposed for the compressible Navier-Stokes equations using an HDG approach in space and a discontinuous approximation in time [266]. More recently, special attention was dedicated to the development of positivity-preserving Riemann solvers in the context of hybridised DG methods [263]. For a complete review on HDG methods for compressible flows, interested readers are referred to [263], whereas the application to gas kinetics modelled by means of the linearised Bhatnagar-Gross-Krook equation is discussed in [254].

Plasma physics and magnetohydrodynamics
Computational physics community is showing increasing interest towards the application of hybrid discretisation methods to the simulation of magnetic plasma physics. Promising preliminary results concerning the HDG approximation of the Grad-Shafranov equation in axisymmetric confinement devices modelling fusion reactors are described in [233,234]. In the context of magnetohydrodynamics (MHD), an HDG method for steady-state linearised incompressible MHD equations is proposed in [180]. Approximation strategies for the unsteady compressible MHD equations using HDG, EDG and the interior embedded DG (IEDG) methods with DIRK time integrators are explored in [74].

Shallow water equations
The shallow water equations have been extensively studied in the context of hybrid DG methods, starting from the linearised shallow water system in [38] to the nonlinear Kortewegde Vries equation in [63,231]. In both the above mentioned works, time integration is performed implicitly using a backward Euler method. Extension to high-order backward differentiation formulas is discussed in [128] in the context of the Benjamin-Bona-Mahony equation. To reduce the computational cost of fully-implicit procedures, in [228] an operator splitting is applied to the Green-Naghdi equation and the nonlinear hyperbolic subproblem is solved using an explicit approach, whereas the implicit time integrator is only applied to the linear dispersive subproblem. A similar idea is presented in [169] to devise an implicit-explicit (IMEX) HDG-DG scheme in which the linear part of the problem is solved using a hybridised DG method and a singly diagonally implicit Runge-Kutta (SDIRK) scheme and the nonlinear one is approximated by means of an explicit Runge-Kutta (RK) DG discretisation. A detailed comparison of explicit and implicit approaches to the nonlinear shallow water equations is provided in [229].

Wave propagation phenomena
The benefits of high-order methods in the simulation of wave propagation prompted extensive research on hybrid discretisation methods in the fields of electromagnetics, elastodynamics and acoustics. A detailed review on HDG and EDG approaches for these problems is available in [132].
Starting from the work in [203], research on time-harmonic Maxwell's equations tackled the analysis and development of HDG formulations [186], including methods suitable for simulations at large wave numbers [189] and Schwarz-type domain decomposition (DD) strategies designed for HDG [18,185]. Recent applications of HDG to time-harmonic Maxwell's equations focus on wave propagation in heterogeneous media modelling photovoltaic cells [41], coupling with nonlocal hydrodynamic Drude and generalised nonlocal optical response models [184] and with hydrodynamic models for metals [257][258][259]271] to simulate plasmonic nanostructures. In the context of time-domain Maxwell's equations, HDG methods are presented and analysed in [55,59,122], whereas implicit hybridised DG discretisations are proposed in [67].
In the framework of elastodynamics, HDG with DIRK time integrators were introduced in [202], whereas in [256] an HDG spectral element method (HDG-SEM) is utilised to simulate wave propagation in coupled elastic-acoustic media. In the frequency-domain, HDG methods for elastodynamics are analysed and presented in [28,164].
The first HDG solver for acoustics, introduced in [202], relied on a fully-implicit approach based on DIRK time integrators. Since then, explicit HDG formulations utilising strong stability-preserving RK (SSPRK) and explicit RK integrators were proposed in [252], wheras an explicit arbitrary derivative (ADER) approach is discussed in [236]. In addition, a comparison of implicit and explicit HDG schemes for acoustic wave propagation is performed in [173]. More recently, an HDG-based cut finite element startegy with local time stepping was presented in [237]. It is worth recalling that devising a conservative numerical scheme is a critical aspect for the accurate simulation of acoustic wave propagation. To correct the dissipative nature of the method analysed in [93], an energy-conservative HDG formulation with a two-step Stormer-Numerov time-marching is proposed in [86]. Moreover, symplectic [232] and multisymplectic [190] HDG schemes preserving the Hamiltonian structure of the PDEs under analysis are developed to achieve energy conservation.
Among the applications of HDG to wave propagation phenomena, it is also worth mentioning the degree adaptive approximation of the mild slope equation to perform harbour simulations [152] and the cardiac electrophysiology simulations of the monodomain model [159,227].

Linear and nonlinear elasticity
In linear elasticity, the imposition of the symmetry of the stress tensor using HDG methods based on mixed formulations has been extensively studied in the literature. Indeed, the first formulations introduced in [141,251] experienced suboptimal convergence of the mixed variable and a loss of superconvergence of the postprocessed displacement field. To remedy this issue, a formulation considering a weakly symmetric stress tensor was presented in [97]. The strong imposition of the symmetry can be achieved via several strategies: in [214], different degrees of polynomial approximation are considered for the primal and hybrid variables; the M -decomposition framework [82,83,85] is applied to the linear elastic problem [81] to enrich the discrete spaces of approximation utilised in the local problem; an alternative formulation imposing the symmetry of the mixed variable pointwise via Voigt notation is proposed for high-order and the lowest-order HDG discretisations in [245] and [244], respectively. In the context of the high-order discretisation, a novel postprocessing strategy accounting for rigid translation and rotation is also devised. It is worth noting that hybrid methods relying on primal formulations do not suffer from these issues, see e.g. HHO [109].
Timoshenko beams are discussed in [47,48], whereas the case of Kirchhoff plates is considered in [162] using HDG and in [27] using HHO methods.
In the context of nonlinear elasticity, the first hybrid discretisation formulation was presented in [167]. In this work, it was observed that the method may not converge to the exact solution if the interelement jumps are not appropriately penalised and a detailed numerical study on the choice of the HDG stabilisation is discussed in [96]. More recently, a locking-free HDG formulation for nonlinear elasticity of thin structures subject to large deformations was proposed [255]. In addition, HHO discretisations of hyperelastic materials in small and finite deformations were presented in [34] and [15], respectively. HHO discretisations for problems involving plastic and elastoplatic simulations are discussed in [16,17], whereas contact phenomena are addressed in [66].

Interface problems and immersed discretisations
The first attempt to solve interface problems using hybrid discretisation techniques was proposed in [165] using a body-fitted mesh. In this context, a superparametric HDG formulation was considered to limit the geometric error due to the polygonal approximation of curved interfaces.
second-order elliptic problems in [39,40]. In the framework of HDG, Poisson interface problems are treated in [121] by means of an unfitted method introducing appropriately defined ansatz functions in the vicinity of the interface. An alternative approach to handle curved interfaces is proposed in [217] where a fictitious domain strategy is developed coupling a mesh of planar faces and a transferring function for the imposition of the transmission conditions on the fictitious subdomain. Inspired by the cut finite element method, in [237], a high-order HDG strategy employing a level-set function to describe the immersed interfaces and a cell agglomeration procedure is described for the wave equation. Similarly, the extended HDG (X-HDG) method introduces a framework in which the HDG local problem is modified only in the elements cut by the interface. In this context, cut instabilities are handled by displacing the mesh nodes responsible for the bad cuts [154][155][156]. Finally, an HDG-based phase-field model for brittle fracture was recently proposed in [194].

High-order and exact geometry representations
Geometry representation plays a crucial role in the capability of high-order methods to achieve optimal accuracy. In the context of HDG, high-order isoparametric approaches in presence of curved meshes are utilised in many references, see e.g. [148,193], whereas this technique is addressed for HHO in [29]. An alternative approach relying on meshes with planar faces and the extension to a fictitious subdomain is discussed in [101,102,250] for several linear problems and was recently extended to the semilinear Grad-Shafranov equation [233,234]. It is worth noting that all the techniques mentioned above introduce geometric errors due to the polynomial approximation of the boundaries. In order to exploit the exact CAD representation of the boundaries, the NURBS-enhanced finite element method (NEFEM) [241,242] is employed in [149,239,247] to devise HDG formulations with exact geometry for Stokes, linear elastic and electrostatics problems, respectively.

Lowest-order hybrid discretisations
Hybrid discretisation methods have been traditionally developed in the context of highorder approximations. Nonetheless, it is well-known that lowest-order discretisations, e.g. the finite volume (FV) method, are more robust than high-order techniques. In this framework, a new class of lowest-order hybrid discretisations was developed, with unknowns approximated by means of constant functions on the mesh faces. The recently proposed face-centred finite volume (FCFV) for Poisson, Stokes [243] and linear elasticity [244] can be interpreted as an HDG method of degree zero. Variants of this approach achieving optimal second-order convergence of the primal variable are discussed in [150,262]. Stemming from HHO, lowest-order nonconforming discretisations are proposed in [32] for linear elasticity and in [72] for elliptic obstacle problems. As their high-order counterparts, the above mentioned methodologies allow the use of generic polygonal and polyhedral elements and provide a workaround to the sensitivity issues of FV methods to mesh distortion and 9 stretching [118,119].

Iterative solvers and preconditioning
Although hybrid discretisation methods are responsible for a substantial reduction of degrees of freedom with respect to classical DG methods, their applicability to realistic problems of engineering interest still rely on the development of efficient solution strategies for large scale systems.
In [26], a DD strategy based on restricted additive Schwarz methods is proposed for hybridised DG approximations, whereas an optimised Schwarz DD approach suitable to handle the many-subdomain case is discussed in [144]. Starting from [80], several works also explored the capabilities of multigrid solvers for HDG formulations, including hierarchical scale separation [238], geometric multigrid [265], nested geometric multigrid on manycore processors [129], p-multigrid in the context of second-order elliptic problems [174] and compressible Navier-Stokes flows [139] and GPU-accelerated p-multigrid for linear elasticity [127]. Finally, iterative algorithms inspired by the Gauss-Seidel method were proposed in [196] and tested on massively parallel architectures up to 16,384 cores. A block symmetric Gauss-Seidel type preconditioner was also introduced in [224], whereas a multilevel solver coupled with a block-Jacobi fine scale solver is proposed in [195].

A posteriori error estimates and adaptivity
The quality of hybrid discretisation methods has been assessed in several works by means of a posteriori estimates of the error in the primal, mixed and hybrid variables, as well as in quantities of interest.
Starting from the seminal works [104,105] establishing reliability and efficiency of error estimates for the HDG approximations of second-order elliptic equations, a posteriori estimates were developed for steady and unsteady scalar convection-diffusion problems in [57] and [182], respectively, and for the vectorial case of incompressible Oseen [23] and Brinkman [22] flows. In addition, constant-free computable a posteriori error estimates are devised in [19] for second-order elliptic problems using an equilibrated fluxes approach, whereas residual-based estimates are established for Maxwell's equations in [58].
In the context of adaptivity, on the one hand, the analysis of HDG approximations based on non-uniform polynomial degrees [60,61] and the superconvergence property of the postprocessed solution [77,89] prompted the development of degree adaptive procedures based on superparametric HDG methods [152,153] and on isoparametric HDG-NEFEM approaches [149,239,247]. Degree adaptivity is also applied to the simulation of cardiac electrophysiology in [159]. On the other hand, mesh adaptivity procedures to capture localised abrupt changes in the solution were devised in [234] and [194] for the Grad-Shafranov equation and the phase-field model for brittle fracture, respectively. Octree-based mesh refinement is performed in [230] for anisotropic inhomogeneous diffusion problems. Mesh adaptivity driven by local error indicators is also employed in the context of second-order FCFV approximations [150,262]. Concerning the error in quantities of interest, an adjointbased method allowing to achieve superconvergent approximations of linear functionals is described in [103] and goal-oriented mesh adaptation strategies are proposed in [136,268].

Coupling HDG with other numerical methods
The accuracy of high-order HDG approximations has been recently exploited to develop efficient algorithms coupling different numerical methodologies in different regions of the computational domain.
In [172], a strategy coupling HDG and a vertex-centred finite volume method is proposed to simulate transient inviscid flows using coarse meshes designed for steady-state problems. In addition, different couplings of HDG and continuous Galerkin (CG) discretisations were explored in the literature. A strategy inspired by a non-overlapping DD method is presented in [208] in the context of incompressible Navier-Stokes flows coupled with conjugate heat transfer phenomena. An alternative minimally-intrusive coupling based on a Nitsche's formulation of the CG method was first introduced in [175] for linear elastic problems involving nearly incompressible materials and was extended to FSI problems with weakly compressible flows in [176].

HDG-based reduced order models
In recent years, the accuracy of the HDG method and its flexibility to devise high-order adaptive discretisation have been also employed to devise high-fidelity reduced and surrogate models. In [260,261], a reduced order model to accelerate the Monte-Carlo simulation of stochastic elliptic PDEs is constructed coupling a high-order HDG method with a reduced basis and empirical interpolation approach. The combination of an HDG solver for time-harmonic Maxwell's equations and a proper orthogonal decomposition (POD) strategy to design parametrised plasmonic nanogap structures is proposed in [259]. An HGD-POD reduced order model (ROM) is also discussed in [249] for the fast simulation of the unsteady heat equation. More recently, an a priori ROM based on HDG and the proper generalised decomposition was proposed to simulate Stokes flows in geometrically parametrised domains [147,240].

HDG formulation of the Poisson equation
In this section, the formulation of the HDG method for the Poisson equation is briefly recalled. Special attention is devoted to the identification of the building blocks of the numerical scheme whose implementation will be detailed in section 8. Interested readers are referred to [89] for a complete theoretical introduction to the HDG method for Poisson equation and to [246] for a tutorial on its derivation.
Let Ω ⊂ R nsd be an open bounded domain in n sd spatial dimensions such that its boundary is ∂Ω=Γ D ∪ Γ N and Γ D ∩ Γ N =∅. The strong form of the Poisson equation is where the unknown u represents the solution field, κ denotes the material parameter (e.g. conductivity in a thermal problem) and s is a volumetric source term. On the boundary, Dirichlet, u D , and Neumann, g, data prescribe the values of the unknown and its flux on Γ D and Γ N , respectively. The vector n denotes the outward unit normal vector to the boundary. Following the HDG rationale [89,92,200,201,204,246], a mixed variable q=− √ κ∇u is introduced and problem (1) is rewritten as a system of first-order equations element-byelement, that is,

HDG local and global problems: strong form
in Ω e , e=1, . . . , n el , in Ω e , e=1, . . . , n el , where the jump operator · is defined as being i and j the evaluations of the quantity in two neighbouring elements Ω i and Ω j sharing a given interface [191]. The last two conditions in (2), known as transmission conditions, enforce the continuity of the solution and of its normal flux across the internal mesh skeleton Γ.
The HDG algorithm solves equation (2) in two stages. First, an independent hybrid variableû is introduced to represent the trace of the solution on ∂Ω e \ Γ D and the primal and mixed variables (u e , q e ) in each element Ω e , e=1, . . . , n el are expressed as functions of the unknownû, namely in Ω e , e=1, . . . , n el , Remark 1. Equation (3) represents the n el HDG local problems. This stage corresponds to the hybridisation of the mixed problem, see [138], and is equivalent to the static condensation procedure in classical continuous Galerkin methods [157].
Second, the hybrid variable is computed by solving the HDG global problem, which accounts for the transmission conditions on the mesh skeleton Γ and the Neumann boundary condition on Γ N , that is, Remark 2. The first condition is automatically fulfilled owing to the Dirichlet boundary condition u e =û on ∂Ω e \ Γ D imposed in the local problem and to the uniqueness of the hybrid variable on each mesh face (respectively, edge in 2D).
The solution (u e , q e ) in each element Ω e , e=1, . . . , n el is thus efficiently retrieved by solving n el independent problems, see equation (3), element-by-element.

HDG local and global problems: weak form
Following the rationale introduced in [246], the discrete functional spaces V h (Ω):={v ∈ L 2 (Ω) : v| Ωe ∈ P p (Ω e )∀Ω e , e=1, . . . , n el }, (5a) are defined for the approximation of the element-based and face-based variables, respectively. In (5), P p (Ω e ) and P p (Γ i ) stand for the spaces of polynomial functions of complete degree at most p in Ω e and on Γ i , respectively.
Similarly, the weak form of the HDG global problem is: for allv ∈ V h (Γ ∪ Γ N ).

HDG local and global problems: discrete form
An isoparametric formulation is considered for the primal, mixed and hybrid variables in the discrete spaces (5), that is, where u i , q i andû i are the nodal values of the unknowns, N i andN i are the polynomial shape functions of degree p defined in a reference element and on a reference face, respectively and n en and n fn denote the number of nodes per element and per face, respectively.
Hence, for each element Ω e , e = 1, . . . , n el , the local problem (6) leads to the linear system of equations from which the following solution is computed with the matrices and the vectors The corresponding discretisation of the global problem (7) leads to By plugging the local elemental solution (10a) into (11), the HDG problem Kû =f involving only the globally-coupled degrees of freedom is obtained, where the matrix and the right-hand side vector are obtained by assembling the elemental contributions The expressions of the matrices and vectors introduced in this section are detailed in appendix A.
The HDG postprocess procedure allows to compute a superconvergent approximation u of the primal variable by solving an independent local problem in each element, namely with the constraint (u , 1) Ωe = (u e , 1) Ωe on the mean value of the solution in the element.
For each element Ω e , e = 1, . . . , n el , the weak form of the postprocess procedure is: for all v ∈ V h (Ω e ).
Using an isoparametric approximation for the functions in the space V h (Ω), the HDG local postprocess gives rise to the linear system where the saddle-point structure of the problem follows from the imposition of the constraint (16b) via the Lagrange multiplier λ and I : V h → V h and I nsd : denote the interpolation operators from the spaces of polynomial functions of degree p to the ones of degree p+1 for scalar and vector-valued functions.
The expressions of the matrices and vectors introduced in this section are detailed in appendix A.

HDG formulation of the Stokes equations
This section presents the formulation of the HDG method for the Stokes equations, extending the framework previously introduced for the Poisson equation. For the sake of simplicity, the present work focuses on the velocity pressure formulation of the Stokes equations. For the Cauchy stress tensor formulation, a tutorial to devise an HDG method based on equal order approximation for all the variables and pointwise symmetric mixed variable is presented in [151].
The open bounded domain Ω ⊂ R nsd is characterised now by a boundary partitioned in three portions disjoint by pairs such that ∂Ω=Γ D ∪ Γ N ∪ Γ S , where Dirichlet, Neumann and slip conditions are imposed. The strong form of the Stokes equations is where the pair (u, p) denotes the unknown velocity and pressure fields, ν > 0 is the kinematic viscosity, n is the outward unit normal to the boundary, s is the vector of the body forces and u D and g represent the imposed velocity and pseudo-traction on the Dirichlet and Neumann boundaries, respectively. On the slip boundary, matrices D and E are defined as D:=[n, βt 1 , . . . , βt nsd−1 ] and E:=[αn, t 1 , . . . , t nsd−1 ], the tangential vectors t j , j=1, . . . , n sd −1 being such that {n, t 1 , . . . , t nsd−1 } form an orthonormal system of vectors. Two scalars, α and β, represent the penetration and friction coefficient, respectively. For α, β → 0, the case of a perfectly slip condition is retrieved [151].
Remark 4. The divergence-free condition in equation (18) induces the following compatibility condition on the velocity field Remark 5. In case of a purely Dirichlet boundary value problem, that is ∂Ω=Γ D , an additional constraint needs to be introduced to retrieve uniqueness of the pressure field. A common approach relies on imposing the mean value of the pressure in the domain [120,218] 1 or on the boundary of the domain [88,98,148,199], namely

HDG local and global problems: strong form
Following the rationale presented in section 3, equation (18) is rewritten element-byelement as a system of first-order equations by introducing a mixed variable L= − √ ν∇u, for e=1, . . . , n el .
First, the HDG algorithm performes the hybridisation step by expressing (u e , p e , L e ) in each element Ω e as functions of the unknown trace of the velocity u on the element faces via the HDG local problem for e=1, . . . , n el . Note that equation (23a) is a purely Dirichlet boundary value problem. Hence, following remark 5, the constraint is introduced, where ρ e is an independent variable representing the mean value of the pressure on the boundary ∂Ω e . It is worth noting that the variable ρ was not present in the HDG approximation of the Poisson equation and its treatment in the HDGlab code will be detailed in section 9.
The HDG global problem thus accounts for the transmission and non-Dirichlet boundary conditions, that is where, following remark 2, the first condition is automatically fulfilled. In addition, the constraint in remark 4 is rewritten element-by-element in terms of the hybrid variable u leading to for e=1, . . . , n el .

HDG local and global problems: weak form
The corresponding weak form of the HDG local problem (23) is:  (24)

HDG local and global problems: discrete form
The discretisation of the local problem (25) leads to the following linear system of equations where the constraint (23b) is imposed using the Lagrange multiplier ζ. It is worth noting that the Lagrange multiplier is required to guarantee that equation (27) is well-posed and the computed pressure field is unique but is not utilised in the solution of the global HDG problem. The resulting local elemental solution is given by with the matrix and vectors defined as The discrete form of the global problem (26) reads as By inserting the solution (28a) into (29), the following system involving the globallycoupled unknownsû and ρ is obtained where the matrices and vectors have the form (31c) The expressions of the matrices and vectors introduced above are detailed in appendix B.
Remark 6. Differently from the Poisson case, the global problem (30) features a saddlepoint structure, as classical in the context of incompressible flows [120]. The proof of the symmetry of the HDG matrix in equation (30) can be devised following the rationale described in [151].

HDG local postprocess
The HDG postprocess procedure presented in section 3.4 for the Poisson equation can be extended straightforwardly to the case of the vectorial variable u.
Remark 7. The postprocessing procedure utilised here is inspired by the work of Stenberg [253] and was exploited in the framework of HDG to obtain an improved approximation of the velocity field via the solution of an additional inexpensive element-by-element problem [199,247]. Nonetheless, in the context of incompressible flows, it is often of interest retrieving an H(div)-conforming and exactly divergence-free approximation of the velocity field. For this purpose, alternative postprocessing strategies inspired by the Brezzi-Douglas-Marini (BDM) projection operator [36] were proposed [90,92]. It is worth noting that the above mentioned procedures are suitable only for the velocity-pressure formulation of the incompressible flow equations. In case a formulation based on the Cauchy stress tensor is considered, an additional constraint is required to handle the rigid rotational modes as discussed in [148,151,245].
A superconvergent velocity field u is thus obtained by solving an independent local problem in each element, namely with the constraint (u , 1) Ωe = (u e , 1) Ωe on the mean value of the velocity in the element.
Hence, for each element Ω e , e=1, . . . , n el , the weak form of the postprocess procedure is: Using an isoparametric approximation for the functions in the space [V h (Ω)] nsd , the HDG local postprocess for the Stokes equations leads to where λ is the vector of Lagrange multipliers imposing the constraint (34b) on the elemental mean value and I nsd×nsd : ] nsd×nsd denotes the interpolation operator for tensorvalued functions from the space of polynomials of degree p to the one of degree p+1.
The expressions of the matrices and vectors introduced above are detailed in appendix B.

The HDGlab repository
The implementation of the HDG solver for the Poisson and Stokes equations has been made available as an open-source software, released under the terms of the GNU General Public License version 3.0 or any later version (https://www.gnu.org/licenses) and is freely available from the repository: https://git.lacan.upc.edu/hybridLab/HDGlab.
The structure of the repository, illustrated in figure 1, is described in this section.
The directory dat contains the data files. This includes two and three dimensional meshes in the directories meshes2D and meshes3D respectively, the pre-computed reference elements in two and three dimensions in the directory refElem and the data structure required to postprocess high-order HDG solutions in the directory postprocess.
The directory Poisson contains the HDG solver for the Poisson problem as described in section 3. The directory hdgPoisson contains the core HDG library for the Poisson equation. The directory testsPoisson is used to organise the functions that describe the setup of the problems to be solved, including the definition of the boundary conditions, source term and, if known, the analytical solution. The directory resPoisson is where the results are saved after an execution of the Poisson solver.
The directory Stokes contains the HDG solver for the Stokes problem as described in section 4. The structure of this directory follows the same rationale as the one corresponding to the Poisson solver.
The directory common contains a set of functions that are common to both the Poisson and the Stokes solvers.
Finally, the directory importMesh contains a library that is provided to import a mesh generated with the open source software Gmsh [146] in HDGlab. The core routines to convert a mesh from .msh to .mat format are located in the directory GMSH. The directory examples contains some test cases including .geo and .msh files, whereas the output of the imported mesh is stored in the directory meshFiles. This library is described in detail in Appendix C.

Data structures
Three data structures are used to manage the mesh, the reference element and the face information required to compute the elemental contributions of the global HDG problem. These three variables are assumed to be an input of the HDG library and a detailed description is provided in this section. To make the developed software more accessible variables of type struct are used in this work rather than class types. 24

Mesh
The variable mesh contains the following information: • nsd: Number of spatial dimensions.
• X: Array of dimension nOfNodes×nsd containing the nodal coordinates of the mesh.
• indexT: Array of dimension nOfElements×2 containing the connectivity indices of the mesh. The first column contains the first node of the element and the second column contains the last node of the element.
• pElem: Array of dimension 1×nOfElements containing the degree of approximation of each element.
• matElem: Array of dimension 1×nOfElements containing the material flag for each element.
• nOfIntFaces: Number of interior faces (i.e. faces not on the boundary).
• intFaces: Array of dimension nOfIntFaces×5. The first two columns contain the first element sharing this face and its local face number. The next two columns contain the second element sharing this face and its local face number. The last column contains the local node number of the second face that matches with the first local node of the first face.
• nOfExtFaces: Number of exterior faces (i.e. faces on the boundary).
• extFaces: Array of dimension nOfIntFaces×4. The first two columns contains the element sharing this face and its local face number. The third column contains the boundary condition flag and the last column contains the flag of the boundary curve or surface.
The information stored in the field intFaces is characteristic of a DG formulation, where integrals on interior faces need to be computed. This is in contrast with a standard CG formulation, where only the field extFaces is needed to impose the boundary conditions. The last column of intFaces is needed to account for the different orientation of an interior face as seen from the element on the left and on the right of a face. It is worth noting that in two dimensions the information could be omitted because the local node  number of the second face that matches with the first local node of the first face is always equal to two, as illustrated in figure 2. However, in three dimensions this information is required as there are three possible rotations of the local face nodes that do not alter the orientation of the element/face, as depicted in figure 3.
To illustrate the mesh data structure, a coarse mesh of the domain Ω = [0, 1] 2 , with four triangular elements is considered. Figure 4 represents the mesh with the global numbering of the mesh elements, the element nodes, the mesh faces and the face nodes.
An overview of the data contained in the mesh data structure for the example of figure 4 is shown in figure 5, with the details of some of the arrays depicted in figure 6.
It is worth noting that the connectivity of an element is obtained by using the function shown in figure 7. This means that the mesh nodes in the array X are duplicated. This functionality is meant to help the handling of meshes with a non-uniform degree of approximation and to provide a seamless way to partition the mesh in case future users would like to parallelise the code. If desired, a different array can be introduced for the connectivity information and the user only needs to redefine the function getElemConnectivity.    Figure 5: Overview of the data contained in the mesh data structure for the example of figure 4.

Reference element
As usual in an isoparametric finite element context, the information related to the approximation and numerical integration in an element is stored by means of a reference element, with local coordinates, ξ = (ξ 1 , . . . , ξ nsd ). To easily handle meshes with a non-uniform degree of approximation, the variable refElem is considered an array of dimension 1 × p max , where p max is the maximum degree of approximation used in all the elements. For each component of the refElem array, the following information is stored: • nsd: Number of spatial dimensions.
• p: Degree of approximation.
• nOfVertices: Number of element vertices.   Figure 7: Function to retrieve the connectivity of an element.
• coordinates: Array of dimension nOfNodes × nsd containing the local nodal coordinates.
• face: Array of dimension 1 × nOfFaces. Each position is a structure containing the following information about an element face: nodes: Array containing the element local number of the nodes in the current face.
-nodesPerm: Array containing nsd + 1 permutations of the face nodes in the field nodes, using the element local number. Each row provides the permutation as required by the field intFaces in the mesh data structure.
-nodesPermHDG: Array containing nsd + 1 permutations of the face nodes in the field nodes, using the face local number. Each row provides the permutation as required by the field intFaces in the mesh data structure.
• shapeFunctions: Array of dimension nOfGauss × nOfNodes × (nsd + 1). When the third index is equal to 1, the array contains the value of the shape functions, for all the nodes at all integration points. When the third index is equal to r > 1, the array contains the value of the derivative of the shape functions in the ξ k direction, for all the nodes at all integration points.
• shapeFunctionsNodesPPp1: Array of dimension nOfNodes × nOfNodes, where nOfNodes denotes the number of nodes of an element with degree of approximation p+1. It contains the value of the shape functions of the current element at the nodes of an element with one extra degree of approximation. This array is only used to perform the HDG postprocess described in sections 3.4 and 4.4.
• shapeFunctionsGaussPPp1: Array of dimension nOfGauss × nOfNodes, where nOfGauss denotes the number of integration points of an element with degree of approximation p+1. It contains the value of the shape functions of the current element and its derivatives at the integration points of an element with one extra degree of approximation. This array is only used to perform the HDG postprocess described in sections 3.4 and 4.4.
To illustrate the refElem data structure, figure 8 depicts the reference quadratic triangular element. An overview of the data contained in the refElem data structure for the example of figure 8 is shown in figure 9.
Similarly, figure 10 shows the reference quadratic tetrahedral element and figure 11 depicts and overview of the data contained in the refElem data structure.
It is worth noting that the code provided utilises nodal basis functions for the polynomial approximation. However, it is straightforward for future users to change the basis used for the approximation by simply changing the reference element information. With minimum effort it is also possible to incorporate other element types.      Figure 11: Overview of the data contained in the refElem data structure for the reference tetrahedral quadratic element of figure 10.

Reference face
The information related to the approximation and numerical integration on a face is stored by means of a reference face, with local coordinates, η = (η 1 , . . . , η nsd−1 ). To easily handle meshes with a non-uniform degree of approximation, the variable refFace is considered an array of dimension p max × p max , where p max is the maximum degree of approximation used in all the elements. For each diagonal component of the refFace array, the information stored is a subset of the information stored in the refElem data structure. As an example, figure 12 offers an overview of the data contained in the diagonal term of refFace corresponding to  Figure 12: Overview of the data contained in the refFace data structure for the reference face of a quadratic element in two dimensions.
a quadratic face of a triangular element.
The upper triangular portion of the refFace data structure contains the information associated with the field shapeFunctions. This is only required when a mesh with a nonuniform degree of approximation is employed. The position (r, s) of the array refFace contains the shape functions of degree r evaluated at the integration points of a face with degree of approximation s. This is required when computing the integrals in an interior face where the degree of approximation of the elements sharing this face is different. It is worth noting that only the entries in the upper triangle of the array refFace contain relevant information because the degree of approximation used for the hybrid variable in a given face is defined as the maximum between the degree of approximation of the two elements sharing the face.

Preprocess
This section describes the preprocessing stage of the HDG solver. The implementation can be found in the function hdgPreprocess, which produces as an output an updated version of the data structures mesh and hdg.
The data structure mesh is taken as an input, containing the fields described in section 6.1, and a new field, called indexTf is added. This field contains an array of dimension nOfFaces × 2 featuring the connectivity indices of the mesh skeleton, where nOfFaces is the total number of mesh faces (i.e. interior faces plus exterior faces). The first column contains the first degree of freedom of a face and the second column contains the last degree of freedom of the face. Figure 13 shows the data contained in indexTf, after the preprocessing is performed, for the mesh of figure 4 and for a Poisson problem (i.e. scalar unknown). The same information for a Stokes problem is shown in figure 14.
For the Poisson problem, each node has one degree of freedom associated as the hybrid variable contains an approximation of the trace of the solution, which is a scalar quantity. In contrast, for the Stokes problem the global vector of unknowns contains an approximation of the trace of the velocity and an approximation of the mean pressure. Therefore, in two    dimensions each face contains 2(p + 1) degrees of freedom for the velocity and one extra degree of freedom per element is introduced for the mean pressure.
The hdg data structure is also an input of the function hdgPreprocess, containing two user defined parameters, namely: • tau: Value of the constant stabilisation parameter.
• problem: String containing the name of the problem to be solved. The code provided contains two model problems, namely 'Poisson' and 'Stokes'.
The structure is updated in the preprocess stage with the following information: • faceInfo: Structure of dimension 1 × nOfElements. For each element of the array, the following information provides a link between the element and face information of the mesh: -local2global: Array of dimension 1×nOfElementFaces, containing the global face number of the faces of the current element.
-localNumFlux: Array of dimension 1 × nOfElementFaces, containing a flag for the numerical flux function associated with the faces of the current element.
For an interior face, a value of zero is set. For a boundary face, the number corresponds to the boundary condition to be imposed and it is linked to the third column of the array extFaces of the mesh data structure described in section 6.1.  Figure 15: Data contained in the hdg data structure for the mesh of figure 4.

Reference element 2D
permutations: Array of dimension 1 × nOfElementFaces, containing a flag for the permutation required to ensure that the ordering of the nodes in the global face matches the ordering of the face nodes in the current element.
-pHat: Array of dimension 1 × nOfElementFaces, containing the degree of approximation used for the hybrid variable in the corresponding faces of the current element.
• nDOFglobal: Number of global degrees of freedom.
• vDOFtoSolve: Number of global degrees of freedom associated with nodes not on a Dirichlet boundary.
In the case of the Stokes equations, three additional fields are introduced in the hdg data structure during the preprocess routine: • pureDirichlet: Boolean variable identifying whether the problem has purely Dirichlet boundary conditions.
• columnsGlobalSystem: Number of columns in the global system, corresponding to the number of unknowns given by the hybrid velocity and the mean pressure.
• rowsGlobalSystem: Number of rows in the global system, including the constraint for the uniqueness of pressure. Note that the value of this variable will differ from columnsGlobalSystem only in the case of purely Dirichlet boundary value problems.
The data contained in the hdg data structure, after the preprocess stage, is shown in figure 15 for the Poisson problem on the mesh of figure 4. The data contained in the field faceInfo for the first two elements is also depicted in figure 16.
Two additional simple data structures are defined at the preprocess stage, namely problemParams and ctt. The structure problemParams contains problem-specific information. The current implementation uses this data structure to carry the following information: • nOfMat: Number of materials in the domain.
• charLength: A characteristic length, used to define the value of the stabilisation parameter of the HDG formulation.  Figure 16: Data contained in the hdg.faceInfo for the first two elements of the mesh of figure 4.
• example: An integer that points to the number of a user-defined example.
In addition, problemParams stores the information on the material parameters. For the Poisson problem: • conductivity: Array of dimension 1×nOfMat that contains the conductivity of each material present in the domain.
For the Stokes problem: • viscosity: A scalar value representing the viscosity coefficient of the fluid.
• alphaSlip: A scalar value describing the penetration coefficient for the slip boundary condition.
• betaSlip: A scalar value describing the friction coefficient for the slip boundary condition.
Finally, the structure ctt contains the following information: • iBC Interior: A flag for the numerical flux function to be used on an interior face.
• iBC Dirichlet: A flag for the numerical flux function to be used on an exterior face where a Dirichlet boundary condition is imposed.
• iBC Neumann: A flag for the numerical flux function to be used on an exterior face where a Neumann boundary condition is imposed.
• iBC Slip: A flag for the numerical flux function to be used on an exterior face where a slip boundary condition is imposed (only supported for the Stokes problem).
• nOfComponents: Number of components of the primal variable.
The flags used to distinguish the type of face and numerical flux to be considered can be specified by the user and they are linked to the third column of the array extFaces of the mesh data structure described in section 6.1.

The HDGlab Poisson solver
A code workflow diagram of the HDGlab Poisson code is shown in figure 17. This section focuses on the core part of the code that involves the assembly and solution of the global system of equations, the element-by-element solution of the local problems and the local postprocess to obtain a superconvergent solution.

Global problem
The sponding to the computation of the element integrals and the face integrals respectively. An extract of this function, showing the computation of the elemental matrices A qq and A uq and the elemental vector f u , is displayed in figure 18. It is worth noting that the loop on integration points is vectorised by using the binary singleton expansion function bsxfun.
In a similar fashion, the second part of the function hdg Poisson ElementalMatrices performs a loop on the faces of the current element and computes the face integrals that lead to the matrices A uu , A uû , A qû and Aûû. This computation distinguishes between interior and exterior faces and, for the exterior faces, utilises the flag in hdg.faceInfo.localNumFlux to enforce the correct boundary condition. For a Dirichlet face, the vector f u is updated and the vector f q is computed. For a Neumann face, the vector fû is computed.
One of the distinctive parts of the HDG formulation is found in the loop on faces, where a vector called indexFlip is used to ensure that the ordering of the degrees of freedom corresponding to the hybrid variable, as seen from the current element, matches the global ordering of the degrees of freedom of the vector of unknowns of the global HDG system. After all the elemental matrices and vectors are computed, the elemental contributions to the global system are prepared to be assembled, namely K e andf e , as shown in the extract of figure 19. It is worth noting that at this stage the matrices Z uû and Z qû and the vectors z f u and z f q , defined in equation (10), are stored in the data structure local in order to be used during the solution of the element-by-element local problems. For large problems, the user might choose to save these matrices to disk before solving the global system of equations. the global vectors for the primal and mixed variables. This is managed by the function hdgElemToFaceIndex, which is designed to work for variables with any number of components. It is also worth noting that this function accounts for the possibility to have a non-uniform degree of approximation.

Local postprocess
As discussed in section 3.4, once the primal and mixed variables are computed, it is possible to perform a local, element-by-element, postprocess procedure to obtain a more accurate, superconvergent, approximation of the solution. This is implemented in the function hdg Poisson LocalPostprocess, shown in figure 21.
A key aspect in the function hdg Poisson LocalPostprocess is the computation of the elemental matrices and vectors in equation (17), which is implemented in the subroutine hdg Poisson LocalPostprocessElemMat.
It is worth emphasising that the implementation assumes that no extra geometric information is known at this stage. Therefore, to compute the high-order nodal distribution of degree p+1, the nodal distribution in the reference element is mapped to the physical space by using the isoparametric mapping of degree p. This implies that, for curved elements, a subparametric formulation is considered. This formulation can lead to a suboptimal rate of convergence for the postprocessed solution as demonstrated in [239,247], where a NURBS-enhanced implementation was proposed.  Figure 20: Function hdg Poisson LocalProblem to solve the element-by-element local problems.

The HDGlab Stokes solver
In this section, the HDGlab solver for the Stokes equations is presented. It is worth noting that the code features the same structure introduced in the previous section for the Poisson case. Hence, only the differences with respect to the Poisson solver will be detailed.

A vector-valued problem
The HDG global system of equations is assembled and solved in the hdg Stokes GlobalSystem function. More precisely, the element and face integrals are computed by the function hdg Stokes ElementalMatrices for each element.

Slip boundary conditions
In the second part of hdg Stokes ElementalMatrices, the face integrals are computed within a loop on the element faces. More precisely, the matrices A Lû , A uu , A uû and A pû are computed for the local problem, whereas the matrix Aûû is computed for the global problem. In addition, on the Neumann faces the vector fû is computed, whereas on the Dirichlet faces the vectors f L and f p are computed and the vector f u is updated.
Remark 8. In case of Dirichlet and Neumann boundary conditions, the remaining matrices involved in the global problem are such that In case slip boundary conditions are also considered, the properties in equation (36) no longer hold. In this context, the matrices Aû L , Aû u and Aû p are computed in the function hdg Stokes ElementalMatrices. Figure 23 displays the initialisation of the above elemental matrices which are then computed in the loop on faces when the flag in hdg.faceInfo.localNumFlux matches ctt.iBC Slip.
A specific treatment of the case in which slip boundary conditions are considered is also required for the definition of the elemental matrices of the global problem with appropriate ordering of the degrees of freedom of the hybrid variable using the indexFlip vector, see figure 24.

Additional constraint in the local problem
A major difference between the Poisson and Stokes case is the structure of the local problems (9) and (27). Despite both matrices are symmetric, the one arising from the discretisation of the Poisson equation is positive definite, whereas it is indefinite in the Stokes case. More precisely, the matrix in (27) features a saddle-point structure, as classical in the context of finite element approximations of incompressible flow problems [120]. In addition, since the HDG local problem is a purely Dirichlet boundary value problem, the constraint (23b) is introduced using a Lagrange multiplier ζ.
The structure of the symmetric indefinite matrix involved in the local problem is displayed on the left-hand side of figure 25. On the right-hand side, the blocks of the first and last columns are associated with the contribution ofû and ρ, respectively, whereas the second column is related to the independent term of the equation. The figure reports the hybridisation stage in which the elemental matrices and vectors defined in (28) are computed for each element. The output of this computation is stored in the data structure local to be successively utilised in the solution of the element-by-element local problems.  Figure 25: Extract of the hdg Stokes ElementalMatrices function that computes the matrices and vectors in equation (27).   Finally, hdg Stokes ElementalMatrices computes the elemental contribution to the matrix in the global problem (30) as reported in figure 26. It is worth noting that the variable hdg.pureDirichlet is utilised here to discriminate the construction of the global matrix of a purely Dirichlet boundary value problem. More precisely, besides the blocks K, G and G T , also the vector a ρρ (i.e. ArlExtra) arising from the imposition of the global constraint (20) is taken into account. An extract of the hdg Stokes ElementalMatrices function computing such a vector is displayed in figure 27. % Initialisation mat.i = zeros(1, hdg.rowsGlobalSystem); mat.j = zeros(1, hdg.columnsGlobalSystem); mat.Kij = zeros(1, hdg.nDOFglobal); f = zeros(hdg.rowsGlobalSystem, 1); 1 Figure 28: Extract of the hdg Stokes Globalystem function that initialises the structures required for the assembly of the global system. Remark 9. The assembly of the block matrix reported in figure 26 does not exploit the symmetry of the terms in equation (30) to its full potential. Indeed, the rationale of HDGlab being educational, the matrix G is constructed by inserting the solution (28a) of the local problem into the global problem (29), leading to

Assembly of the global system
As described in section 4, the global problem for the Stokes equations features a saddlepoint structure. Hence, the assembly of the global system presents major differences with respect to the Poisson case previously discussed. First, figure 28 reports the initialisation of the structures required to perform the assembly. It is worth recalling that the value of hdg.rowsGlobalSystem differs from the one of hdg.columnsGlobalSystem only if Dirichlet conditions are imposed on all boundaries. In this case, the additional constraint (20) is considered to guarantee the well-posedness of the problem.
The construction of the structures to perform the assembly of the global system is displayed in figure 29. According to the variable hdg.pureDirichlet, the above mentioned constraint is introduced as an extra line in the system.
Of course, the matrix arising from the operations displayed in figure 29 is rectangular. In order to retrieve a square matrix, an extra column is added to the matrix and the constraint is imposed via a Lagrange multiplier (Fig. 30).

Three unknowns in the local problem
The structure of the code for the local problem in the Stokes case replicates the one presented in figure 20  if hdg.pureDirichlet f(hdg.rowsGlobalSystem) = .
..  in the computation of three variables in each element, namely the pressure p, the velocity u and the gradient of velocity L. An extract of the function hdg Stokes LocalProblem is displayed in figure 31, focusing on the operations in equation (28).

Visualisation
The use of a high-order functional approximation means that non-standard techniques are required to produce a reliable representation of the solution in each element. With the increased popularity of high-order methods in recent years, different strategies to efficiently display high-order solutions have been proposed [188,198,220]. The HDGlab provides the capability to accurately display high-order solutions, including curved isoparametric elements.
Three data structures are employed in the visualisation. The data structure plotOpts is used to collect all the user defined options available for the visualisation. It contains the following information: • resolution: Takes value 1 for a faster but less accurate representation of high-order solutions and value 2 for a slower but more accurate representation of high-order solutions.
• fieldsWithMesh: Plots the solution whilst displaying the mesh.
• fieldsWithNodes: Plots the solution whilst displaying the nodes.
• componentsToPlot: A user-defined vector containing the components of the solution to be represented.
In addition, for the Stokes equation, two problem-specific options are available in the data structure plotOpts: • componentsU: A boolean variable allowing the predefined visualisation of all the components of the velocity vector. This functionality relies on the definition of the vector componentsToPlot.
• moduleU: A boolean variable allowing the visualisation of the module of the velocity.
The data structure postproc is provided for triangular and tetrahedral elements with equally-spaced and Fekete nodal sets in the directory dat/postprocess. In two dimensions, the data structure postproc contains the following information: • nOfNodesPlot: Number of nodes used to display the high-order solution in each element.
• nodesPlot: Array of dimension nOfNodesPlot × 2. Coordinates of the nodes, in the reference element, used to display the high-order solution. Figure 32: Submesh of the reference element used to provide an accurate representation of high-order solutions in each element. The left picture corresponds to resolution=1 and the right picture to resolution=2.
• nSubElemsOnePlot: Number of subelements used to display the high-order solution in each element.
• connecNodesPlot: Array of dimension nSubElemsOnePlot × 3 accounting for the connectivity of the submesh used to display the high-order solution in each element.
• nOfEdges: Number of edges of the element.
• edgeNodesSplit: Array of dimension 5r × nOfEdges, where r is the resolution selected by the user in the data structure plotOpts. The i-th column contains the local number of the list of nodesPlot that belong to the i-th edge.
• elem: Array of dimension 1×p max , where p max is the maximum degree of approximation used in all the elements. The component p of elem contains a field called N, of dimension nOfNodesPlot × p(p + 1)/2, that stores the value of the shape functions of order p at the positions given by nodesPlot. This information is used to interpolate the solution at the nodes of the submesh, providing a more accurate representation of the high-order solution.
• face: Array of dimension 1 × p max . The component p of elem contains a field called N, of dimension nOfNodesPlot × (p + 1), that stores the value of the shape functions of order p at the positions of an edge. This information is used to interpolate the solution at the edges of the submesh.
The submeshes used for a triangular element with resolution=1 and resolution=2 are displayed in figure 32. To illustrate the effect of the user-defined parameter resolution on the visualisation, figure 33 depicts the shape function associated to the fourth node of the reference quadratic triangular element, shown in figure 8, using resolution=1 and resolution=2.
In three dimensions the data structure postproc contains the following information: • Face: Structure that contains: Figure 33: Shape function of the fourth node of a quadratic triangular element using resolution=1 (left) and resolution=2 (right).
-nElemPlot: Number of subelements used to display the high-order solution in each element face.
-connecPlot: Array of dimension nElemPlot × 3 accounting for the connectivity of the submesh used to display the high-order solution in each element face.
-nNodesPlot: Number of nodes used to display the high-order solution on each element face. It is given by (5r + 1)(5r + 2)/2, where r is the resolution selected by the user in the data structure plotOpts.
-edgeNodesPlot: Array of dimension 1 × (15r + 1), where r is the resolution selected by the user in the data structure plotOpts. It contains the local number of the face nodes that belong to the edges of the face. coord: Array of dimension p(p+1)(p+2)/6 that contains the nodal distribution on the reference tetrahedral for an approximation of degree p.
• faceVertices: Array of dimension 4 × 3. The i-th row contains the local number of the vertices on the i-th face.
The third data structure utilised during the postprocess stage is called visual and it is built by the function buildSubmeshPostprocess2D in two dimensions and by the routine buildSubmeshPostprocess3D in three dimensions. This data structure contains the following information: • X: Physical coordinates of all the nodes of the submesh used to display the high-order solution.
• T: Connectivities of all the subelements used to display the high-order solution.
• Xnodes: Physical coordinates of all the vertices of the mesh. This field is only used if the user sets fieldsWithNodes=1 in the data structure plotOpts.
• edges: Structure containing the list of nodes of the submesh that form the high-order representation of the physical edges of the mesh.
In a separate function, the data structure is updated by adding the field U that contains the interpolated values of the solution on the submesh used to display the high-order solution. This action is performed by the function called interpolateSolutionPostprocess2D and interpolateSolutionPostprocess3D in two and three dimensions respectively.
It is worth noting that the function that creates the data structure visual is independent on the field to be represented and, therefore, it is only called once. Instead, the second function that updates the data structure visual with the field U depends upon the field to be represented. Therefore, several calls can be made to the function updating visual without the need to build the submesh again. It is also important to note that the function that updates the data structure visual with the field U accepts elemental and nodal fields.
Once all information is available in the data structure visual, the postprocessField2D function is used to plot the high-order solution.
In three dimensions there is an extra option available that consists of representing the solution only in a region of the computational domain. The user can set the value of a string, called conditionPlot, that specifies a region in the physical space. Before visual is computed, the function selectFacesToPlot3D computes the list of faces in the computational mesh that satisfy the condition given by conditionPlot. Then, the submesh and interpolation of the solution is only performed over the faces that satisfy the condition specified by the user.
Finally, the visualisation function also accounts for the need to represent the superconvergent solution obtained after the local postprocess described in sections 3.4 and 4.4. To simplify the implementation, the visualisation builds a mesh data structure where the degree of approximation in each element is the degree used for the computation plus one. With this information, the same functions used to display the high-order primal solution can be used to postprocess the higher order superconvergent solution. where is a characteristic length of the domain and c P a scaling factor selected equal to 1 [175]. Following [151], the stabilisation for the Stokes cases is defined as τ =c S ν/ , , the scaling factor being c S =3.

Optimal convergence properties
The optimal convergence properties of the proposed HDG implementation are presented for the Stokes flow, using test cases with analytical solution, in two and three dimensions. Uniform meshes of triangular and tetrahedral elements with Fekete nodal sets are utilised.
First, the two-dimensional Wang flow [264] in the unit square domain Ω=[0, 1] 2 is considered. The analytical velocity field is whereas the pressure field is uniformly zero in Ω. The coefficients a, b and λ in (38) are selected such that a=b=1 and λ=10 and the kinematic viscosity ν is set to 1. The source term s and the boundary conditions are computed starting from the analytical solution above. More precisely, a pseudo-traction g is applied on the bottom surface Γ N :={(x 1 , x 2 ) ∈ Ω | x 2 =0} and Dirichlet data u D are imposed on the remaining boundaries Γ D =∂Ω \ Γ N . Figure 34 displays the convergence history of the relative error, measured in the L 2 (Ω) norm, of the primal, mixed and postprocessed variables as a function of the characteristic mesh size. Optimal convergence of order p+1 is observed for velocity, u, pressure, p, and gradient of velocity, L, whereas superconvergence of order p+2 is achieved by the postprocessed velocity u . The following example involves a three-dimensional Stokes flow in the unit cube Ω=[0, 1] 3 , with the following manufactured solution where m=1 and n= 1 2 . The kinematic viscosity is set to ν=1, a Neumann datum is imposed on the boundary Γ N :={(x 1 , x 2 , x 3 ) ∈ Ω | x 3 =0}, whereas Dirichlet conditions are prescribed on Γ D =∂Ω \ Γ N .
The optimal convergence of order p+1 of the relative L 2 (Ω) error for velocity, pressure and gradient of velocity and the superconvergence of the postprocessed velocity are confirmed in the 3D case by the results in figure 35.

High-order curved meshes
The coaxial Couette flow [64] is considered to show the optimal convergence properties of the HDGlab solver using a high-order isoparametric approximation in a domain featuring curved boundaries.
This test consists of an incompressible viscous flow within two coaxial circular cylinders of infinite length and radius R int =1 and R ext =5, respectively. The computational domain is defined as a section of the 3D cylinders, that is, Ω={(x 1 , Of course, being a purely Dirichlet boundary value problem, the constraint on the mean value of pressure is introduced to enforce the field to be uniformly equal to 1 in the domain.
A set of high-order uniformly refined meshes with Fekete nodal distribution is constructed using the strategy described in [212]. Figure 36 displays the first level of mesh refinement featuring 128 triangular elements of polynomial degree 3 and the module of the computed velocity.
The convergence of the relative error, measured in the L 2 (Ω) norm, of the primal, mixed and postprocessed variables is reported in figure 37 as a function of the characteristic mesh size. Optimal convergence of the primal and mixed variables and superconvergence of the postprocessed variable is achieved also in presence of high-order curved meshes.

Non-uniform degree of approximation
In this section, the flexibility of HDGlab to devise a non-uniform polynomial degree approximation in the domain is presented. This case, inspired by the study on micromixers in [131], consists of the flow in a microchannel with five obstacles. The problem setup features a parabolic inlet velocity profile and homogeneous Dirichlet and Neumann conditions on the top/bottom walls and on the outlet, respectively. The channel has dimensions [0, 6.6] × [−0.5, 0.5] and the obstacles, attached to the top and bottom walls have thickness 0.2 and height 0.5. A mesh with local element size ranging between 0.08 and 0.19 is generated without any specific a priori refinement. It is worth noting that only two mesh elements are defined along the thickness of the obstacles.
To capture the complex flow features among the obstacles, a high-order non-uniform polynomial degree distribution is generated following the adaptivity strategy described in [247]. The resulting degree of approximation in each element is displayed in figure 38. High-order polynomials are employed in correspondance of the tip of the obstacles where localised flow features appear, whereas low-order approximations are utilised in region further away.
The module of the velocity on the described mesh is also reported in figure 38. It is worth noting that the mesh structure featuring the non-uniform polynomial approximation is provided as a datum for this test case. The corresponding simulation is performed seamlessly in HDGlab and no specific intervention is required to the user.   A high-order mesh featuring 1,036 tetrahedral elements is generated via the solid mechanics analogy described in [212,269]. Figure 39 displays the module of the velocity and the pressure field computed using an isoparametric approximation of degree 6.

Stokes flow past a sphere
It is worth noting that the computations in sections 11.1, 11.2, 11.3 and 11.4 are performed using a unique HDGlab solver for the Stokes equations, independently on the number of spatial dimensions of the problem under analysis. Indeed, HDGlab provides a seamless implementation of the HDG method, in which all relevant information is extracted from the mesh structure and the user is required to specify only the physical parameters and the boundary conditions to setup a test case.

Applications of the Poisson solver
The next example, taken from [187], shows the solution of an electrostatic problem governed by the Poisson equation in three dimensions. The domain of interest corresponds to the exterior of 11 conducting spheres and it is discretised with a mesh of 35,895 quadratic tetrahedral elements, as shown in figure 40. This figure shows one of the implemented capabilities of the postprocessing library to display the faces corresponding to the exterior and interior faces of the mesh separately, with different colour and transparencies in each case. The first plot in figure 40 is produced by selecting the faces corresponding to the far field boundary and using a transparency. In a second phase, the exterior faces corresponding to the conducting spheres are displayed with no transparency. It is worth noting that the far field boundary and the boundary corresponding to the conducting spheres could be distinguished using either the boundary condition flag or simply imposing a condition on the faces to be displayed. The second plot of figure 40 is also obtained in two stages.    First, the interior faces satisfying a condition corresponding to a positive x 2 coordinate are displayed with a transparency. Second, the exterior faces corresponding to the conducting spheres is displayed with no transparency. When representing the interior and exterior faces, a constant element field is used to assign different colours and aid the visualisation.
Dirichlet boundary conditions are considered in the whole boundary of the computational domain. A positive electrostatic potential of magnitude 5 is imposed on the central sphere and five of the surrounding spheres, whereas a negative potential of magnitude −5 is imposed on the remaining spheres. On the outer boundary a zero potential is imposed. Figure 41(a) shows the 11 conducting spheres coloured according to the boundary condition imposed. This figure is produced by using the boundary condition flag to select the first set and the second set of spheres separately. After solving the problem with HDGlab, not only the primal variable corresponding to the electrostatic potential is obtained but also its gradient, which corresponds to the electric field in this example. Figure 41(b) shows the magnitude of the electric field on the surface of the 11 spheres.
The following example involves the solution of a heat transfer problem in a three dimensional mechanical component. The domain is discretised with a mesh of 6,465 elements with p = 4, as illustrated in figure 42. The first plot in figure 42 includes the high-order nodal distribution and the second plot in figure 42 shows the possibilities offered by the postprocessing library provided within HDGlab. In fact, it shows an extra functionality that can be added to display the intersection curves of the underlying CAD geometry when the user have access to this data.
Dirichlet and Neumann boundary conditions are imposed on the blue and red portions of the boundary, respectively, as shown in figure 43(a). In the Dirichlet part of the boundary a temperature equal to 10 is imposed, whereas the Neumann boundary condition is homogeneous, imposing that part of the boundary is perfectly insulated. The temperature     field obtained after solving the problem with HDGlab is shown in figure 43(b).
The last example considers the computation of the potential flow past a generic unmanned aerial vehicle (UAV). The domain is discretised using 101,923 tetrahedra of degree p=2, as represented in figure 44. A Neumann boundary condition, corresponding to a unit velocity, is imposed on the inflow part of the boundary and a Dirichlet boundary condition corresponding to zero potential on the outflow. On the rest of the boundary a homogeneous Neumann boundary condition is enforced to represent a physical wall. After solving the problem with HDGlab, the flow potential and the velocity field are obtained. Figure 45 shows the magnitude of the velocity field on the surface of the UAV and the pressure field, computed using the Bernoulli equation.
This last example shows the applicability of the developed HDGlab to problems involving complex geometries, where the resulting global system has more than one million of equations. Despite the code was not developed with computational efficiency in mind, it still enables the interested users to solve relatively large problems and, with some additional improvements in terms of performance, it can be used for larger three dimensional problems.

Concluding remarks
An open-source Matlab implementation of the HDG method for elliptic problems has been presented, the so-called HDGlab. The code is capable of solving problems governed by 57 the Poisson and Stokes equations using high-order simplicial elements, including curved elements, by means of an isoparametric formulation.
HDGlab provides a suitable environment to those interested in the programming of hybrid methods in general and the HDG method in particular. The paper describes the HDG formulation employed for the Poisson and Stokes solvers and provides a very detailed description of the code, with particular emphasis in the data structures utilised. The presentation of the code involves the preprocess, computation and postprocess stages, and includes detailed explanations of the most relevant functions. A set of examples is presented to illustrate the use and the potential of HDGlab. The examples go from simple test cases with a known analytical solution, used to demonstrate the numerical properties of the HDG method, to more complex such as the computation of the potential flow around a UAV using high-order curved meshes.
HDGlab is available as an open-source software, released under the terms of the GNU General Public License version 3.0 or any later version (https://www.gnu.org/licenses) and is freely available from the repository: https://git.lacan.upc.edu/hybridLab/ HDGlab.

A Matrices and vectors for the Poisson solver
In this appendix, the expressions of the matrices and vectors appearing in the discrete form of the HDG local and global problems for the Poisson equation are presented.
For the sake of readability, introduce the compact forms of the shape functions and their derivatives, namely where n denotes the outward unit normal vector to a face and J is the Jacobian of the isoparametric mapping.
The solution of the HDG local problem (9)  N where n e fa denotes the number of faces of the element Ω e , ξ e g are the n e ip integration points defined on the reference element and ξ f g are the n f ip integration points of the reference face. The corresponding weights for the integration points are denoted by w e g and w f g , respectively. In addition, the indicator function χ f is defined as Similarly, for the HDG global problem (11)  N To describe the terms appearing in the HDG postprocess (17), the compact forms ..,n en denote the shape functions for the discretisation of u , J is the Jacobian of the corresponding isoparametric transformation and n en is the number of elemental nodes of the approximation.
The following matrices and vector are thus defined as N (ξ g )|J (ξ g )|w g , where ξ g and w g are the n ip integration points and the weights associated with the higher order approximation in the space V h .

B Matrices and vectors for the Stokes solver
In this appendix, the expressions of the matrices and vectors appearing in the discrete form of the HDG local and global problems for the Stokes equations are presented. It is worth recalling that in this problem, the primal, u, and hybrid,û, variables are n sd -dimensional vector-valued unknowns and the mixed variable L is a tensor-valued unknown of dimension n sd ×n sd .
In addition to the compact forms of the shape functions presented in equations (42) and (44), the following definitions are introduced for the vectorial problem where m sd :=n 2 sd . Moreover, for each component j=1, . . . , n sd , the compact forms N n j := N 1 n j N 2 n j . . . N nen n j T , (45f) N Similarly, the matrices and vectors for the global problem (29) are In addition, if a problem with purely Dirichlet boundary conditions is solved, the con- which is used to enforce the constraint using an appropriately defined Lagrange multiplier.
Finally, for the postprocess (32), the following matrices are required

C An interface with the mesh generator Gmsh
In this appendix, the interface between the high-order open-source mesh generator Gmsh [146] and HDGlab is described. Starting from the structure of the variable mesh introduced in section 6.1, the algorithm to convert a mesh file from the .msh to the .mat format of HDGlab is presented (Fig. 46).
Remark 10. The routines described in this appendix are designed starting from the .msh ASCII file format version 2.2.
First, it is worth recalling that HDGlab offers the feature of utilising either equallyspaced or Fekete points for the approximation. Gmsh provides mesh discretisation based on equally-spaced nodal sets, whence the variable optionNodes is set to 0, see section 6.1.
The function convertMSHtoMAT requires the definition of the following data concerning the mesh to be imported: • fileName: String that specifies the name of the .msh mesh file without extension.
The mesh files are archived in the directory example.
• nsd: Scalar variable defining the number of spatial dimensions (either 2 or 3).
• pDegree: Scalar variable defining the degree of the polynomial order of the mesh.
• outputPath: String that specifies the location where the imported mesh will be stored. The default location is the directory meshFiles.
Given the number of spatial dimensions and the polynomial degree defined above, the refElem and refFace data structures are constructed.

C.1 Reading the .msh file
The .msh file is read by the function scanMeshFileMSH and the following information is extracted: • X: Array containing the coordinates of the mesh nodes.
• T: Connectivity matrix featuring the identifiers of the nodes associated with each element.
• faces: Array containing the information on the boundary faces, namely the identifier of the element the face belongs to, the list of vertices, the flags of the boundary and of the type of boundary condition imposed.
• matElem: Array containing the flag of the region each element belongs to.  Figure 47: Extract of the scanMeshFileMSH function that constructs the permutation vector for the appropriate renumbering of the mesh nodes.
A detailed description of the structure of the .msh file as well as the notation utilised by the mesh generator to identify element and face types is available in the official Gmsh documentation [13]. It is worth noting that the numbering of the nodes in HDGlab differs from the one provided by Gmsh. Hence, an appropriate permutation is performed as reported in figure 47.

C.2 Retrieving the information on the faces
As described in section 6.1, the interior and exterior faces of the mesh are stored in the data structures intFaces and extFaces, respectively.
To construct the structure intFaces, the getInternalFaces function explores the mesh and, for each face of each element, it identifies the neighbouring element by comparing the list of face nodes, see figure 48. In order to optimise this operation, a list of potential neighbours is preliminarily stored by identifying the list of elements containing each node, as detailed by the extract in figure 49.
The function getBoundaryFaces, not reported here for brevity, is responsible for constructing the data structure of the exterior faces. More precisely, the structure extFaces is obtained by rearranging the information previously stored in the data structure faces, according to the rationale described in section 6.1.

C.3 Assemblying the mesh structure
The structure mesh is finally assembled by the buildMeshStruct function using the information obtained by the mesh generator, namely X, T and matElem, and the structures of the interior and exterior faces, intFaces and extFaces, previously generated. More precisely, the connectivity matrix is constructed via a loop on the mesh elements as reported in figure 50.
It is worth noting that the framework provided in HDGlab to import meshes generated using Gmsh can be easily extended to any mesh generator by introducing appropriate functions to construct the intFaces and extFaces data structures, as described above.  hdg.pureDirichlet f(hdg.rowsGlobalSystem) = . ..
edges to appear in the interior of the domain as well (Fig. 51a-b). In both cases, the nodes are equally-spaced.
In addition, the L 2 (Ω) error of the pressure and the gradient of velocity (Fig 51c-d) and the primal and postprocessed velocities (Fig 51e-f) is reported as a function of the characteristic mesh size. The results diplay the expected optimal convergence of order p+1 for pressure, velocity and gradient of velocity and the superconvergence of the postprocessed variable u , showing the capability of HDGlab to work using both equally-spaced and Fekete nodal distributions. It is worth noting that the results obtained using Fekete nodal sets (Fig. 37) provides up to one order of magnitude of extra accuracy for a given mesh size compared to the equally-spaced nodes in figure 51.