HDGlab: An open-source implementation of the hybridisable discontinuous Galerkin method in MATLAB
HHDGlab : An open-source implementation of thehybridisable discontinuous Galerkin method inMATLAB
M. Giacomini , R. Sevilla * and A. Huerta
September 21, 2020
Abstract
This paper presents
HDGlab , an open source MATLAB implementation of thehybridisable discontinuous Galerkin (HDG) method. The main goal is to provide adetailed description of both the HDG method for elliptic problems and its implemen-tation available in
HDGlab . Ultimately, this is expected to make this relatively newadvanced discretisation method more accessible to the computational engineeringcommunity.
HDGlab presents some features not available in other implementationsof the HDG method that can be found in the free domain. First, it implementshigh-order polynomial shape functions up to degree nine, with both equally-spacedand Fekete nodal distributions. Second, it supports curved isoparametric simplicialelements in two and three dimensions. Third, it supports non-uniform degree polyno-mial approximations and it provides a flexible structure to devise degree adaptivitystrategies. Finally, an interface with the open-source high-order mesh generator
Gmsh is provided to facilitate its application to practical engineering problems.
Keywords: hybridizable discontinuous Galerkin, high-order, elliptic problems,MATLAB, open-source Laboratori de C`alcul Num`eric (LaC`aN), ETS de Ingenieros de Caminos, Canales y Puertos, UniversitatPolit`ecnica de Catalunya, Barcelona, Spain International Centre for Numerical Methods in Engineering (CIMNE), Barcelona, Spain. Zienkiewicz Centre for Computational Engineering, College of Engineering, Swansea University, Wales,UK* Corresponding author: Ruben Sevilla.
E-mail: [email protected] a r X i v : . [ c s . M S ] S e p Introduction
In recent years, hybrid discretisation methods have received increasing attention by theapplied mathematics and computational engineering community. The main interest inthese methodologies is due to their reduced computational cost with respect to classicaldiscontinuous Galerkin (DG) methods, see [135, 163, 170, 267], from which they inheritappealing stability and convergence properties as well as the flexibility to devise high-order, non-uniform degree and adaptive discretisations and the capability to efficientlyexploit parallel computing architectures [42, 91, 108, 158, 226].The purpose of the present contribution is two-fold: to present a review on the state-of-the-art of hybrid discretisation methods including both fundamental and applied contribu-tions; to provide an educational implementation of the hybridisable discontinuous Galerkin(HDG) method in MATLAB, the so-called
HDGlab library, and describe its structure, ca-pabilities and functioning.
HDGlab is an open-source library released under GNU GPLlicence and designed for rapid prototyping and testing. It supports simplicial meshes andit provides a seamless 2D and 3D implementation with vectorised loops on the integra-tion points. In addition,
HDGlab presents four specific features, currently not available inexisting open-source HDG implementations in MATLAB:1. Availability of high-order polynomial shape functions up to degree 9, with bothequally-spaced and Fekete nodal distributions.2. Support of curved isoparametric simplicial elements in 2D and 3D.3. Support of non-uniform degree polynomial approximations and flexibility to devisedegree 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-of-the-art on hybrid discretisation methods is presented in section 2. The formulation ofthe HDG method for the Poisson and Stokes problems is briefly recalled in sections 3and 4, respectively. Section 5 provides a description of the structure of the
HDGlab libraryand the url of the repository available under GNU GPL licence. The data structuresfor the storage of the mesh information, the reference element and the reference face arepresented in section 6. Section 7 is devoted to the preprocessing operations, whereasthe 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 insection 9. The visualisation library is introduced in section 10. Section 11 is devotedto numerical examples, in 2D and 3D, validating the optimal convergence properties ofthe HDG method and showing the potentialities of the
HDGlab implementation. Finally,section 12 summarises the capabilities of the presented library and three appendices provideimplementation details for the Poisson (Appendix A) and Stokes (Appendix B) solvers andfor the interface with the mesh generator
Gmsh (Appendix C).2
Literature review
The common idea of all hybrid discretisation methods stems from the seminal works ofGuyan on static condensation of primal formulations [157] and of Fraeijs de Veubeke onhybridisation of mixed formulations [138] of the finite element method. In the contextof element-by-element discontinuous approximations, these techniques allow to remedythe drawback of node duplication in DG methods by considering only the unknowns onthe mesh faces (edges in 2D) as globally-coupled degrees of freedom. More precisely, theunknowns in each element are expressed as a function of the degrees of freedom on theelement faces by solving a local boundary value problem with purely Dirichlet data, whereasappropriate transmission conditions are imposed to guarantee the interelement continuityof the solution and the fluxes, see [75].Three families of hybrid numerical schemes lay within this description, namely, (i) hy-brid/hybridised DG, (ii) hybridisable DG, henceforth referred to as HDG, and (iii) hybridhigh order (HHO) methods. Stemming from classical DG primal formulations, the hybridor hybridised DG method reduces the number of globally coupled degrees of freedom byperforming static condensation [124–126]. In addition, improved efficiency can be achievedusing polynomial spaces of degree p +1 and p for the primal and hybrid variables, respec-tively and resorting to the reduced stabilisation approach [206, 207]. The hybridisable DGmethod, henceforth named HDG, is derived from the mixed formulation of the local DGmethod [77, 87, 99] with hybridisation. The main advantage of HDG with respect to otherhybridised DG methods relies in the introduction of a mixed variable approximating thegradient of the primal unknown [88, 89]. This approach is of special interest in the contextof engineering problems where quantities of interest often depend on the flux of the solu-tion or on the stress. Finally, HHO bridges the two approaches above by utilising a primalformulation and introducing a local reconstruction operator for the gradient of the solu-tion 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 uniqueframework as HDG-type methods via appropriate definitions of the stabilisation term, seee.g. [68, 69] for the staggered DG method and [106] for HHO.Unified presentations of hybrid discretisation techniques and their relationship withother known numerical methods are available in [37, 89, 114]. Interested readers are alsoreferred to the review papers [75, 149] and to the recent monograph [112]. In the followingsubsections, an overview of the contributions on hybrid discretisation methods accordingto the authors’ vision is presented. 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 convection-diffusion problems is discussed in [79,113,124,200]. Cases of higher-order partial differential3quations (PDEs) are presented in [78] and [62] for HDG discretisations of biharmonic andthird-order equations, respectively, whereas an HHO approximation of the Cahn-Hilliardequation is proposed in [53]. In addition, time-fractional diffusion problems are discussedin [76, 197].More recently, there has been growing interest towards the analysis and simulation ofquasilinear and semilinear problems, including the quasilinear p -Laplace operator [95, 216]and the semilinear Grad-Shafranov equation [233, 234]. To reduce the computational costof semilinear problems, the interpolatory HDG method was recently devised introducing aninterpolation procedure for the efficient and accurate approximation of nonlinear terms [54,100].Concerning nonlinear problems, HDG discretisations were proposed for nonlinear con-vection diffusion [201] and nonlinear Schr¨odinger [46] equations, whereas an HHO formu-lation of the nonlinear Leray-Lions equation is presented in [111]. Recent applicationsinvolving HDG approximations of nonlinear scalar equations focus on the optoelectronicsimulation of photovoltaic solar cells. This problem couples a high-order HDG method forthe drift-diffusion electronic model in the semiconductor layer of the solar cells with anefficient approximation of the time-harmonic Maxwell’s equations [21, 56]. In the context of incompressible flows, HDG formulations of the Stokes equations weredevised and analysed in [90, 92, 98, 199]. The corresponding analysis of the HDG methodfor Oseen flow is presented in [49]. In [92], it was observed that the HDG method basedon Cauchy stress tensor formulation experiences suboptimal convergence of the mixedvariable and loss of superconvergence of the postprocessed velocity, when low-order poly-nomial approximations are considered. The M -decomposition approach [84] remedies thisissue by appropriately enriching the discrete local spaces of approximation. An alterna-tive strategy imposing the symmetry of the mixed variable pointwise via Voigt notationis discussed in [148] along with a postprocessing procedure to handle translational androtational rigid body modes. Divergence-conforming HDG [94], hybridised DG [125] andembedded-hybridised DG (EHDG) [225] discretisations were also studied for the incom-pressible Stokes equations, whereas a pressure-robust HHO method for viscosity-dependentStokes flows is proposed in [116].HDG formulations for the nonlinear incompressible Navier-Stokes equations using equalorder and different order of polynomial approximations for the primal, mixed and hybridvariables are described in [50, 204] and [215], respectively. The former approach is alsoemployed in [152] to devise a degree adaptive strategy relying on the local superconver-gence of the postprocessed velocity. Stemming from the work in [117], different HHOformulations of the incompressible Navier-Stokes equations were proposed, incorporatinga skew-symmetric form of the convection term [31] and a globally divergence-free veloc-ity approximation to achieve robustness in presence of large irrotational body forces [45].4oreover, special attention was devoted in recent years to the development of hybridisedDG schemes [126] with pointwise divergence-free velocity [171, 181, 223] and with relaxed H (div)-conformity [177, 178], as well as divergence-conforming hybrid DG discretisationsfor incompressible flows on surfaces [179]. It is worth noting that all the above mentionedreferences focus on viscous laminar flows and preliminary promising results on the incom-pressible Reynolds averaged Navier-Stokes (RANS) equations coupled with the SpalartAllmaras turbulence model were recently presented in [210].Besides classical approaches to steady and unsteady Navier-Stokes equations, HDG-based space-time formulations were studied for their ability to effectively handle movingand deforming domains. More precisely, stemming from the HDG formulation introducedin [222], H (div)-conforming hybridised DG [160] and EHDG [161] methods were proposed.Hybridised DG and HDG methods with arbitrary Lagrangian Eulerian (ALE) formula-tions were thus presented in [134, 140] and the resulting HDG-ALE framework was appliedto fluid-structure interaction (FSI) problems involving incompressible [248] and weakly-compressible flows [176].Among the applications of hybrid discretisation methods to incompressible flows, it isalso worth mentioning the recent attempts to simulate quasi-Newtonian fluids [145] andviscoplastic materials [44]. HDG simulations of immiscible incompressible two-phase flows in heterogeneous porousmedia were first proposed in [130] and coupled with high-order diagonally implicit Runge-Kutta (DIRK) time integrators in [107]. Moreover, in [168] a linear degenerate ellipticproblem modelling two-phase mixture is approximated using a hybridised DG approach.Darcy flow and two-phase flow simulations in highly heterogeneous media are performedin [270] via the so-called generalised multiscale HDG (GMsHDG) method which is con-nected to the mortar mixed finite element method described in [24]. GMsHDG was alsoemployed for multiscale simulations of elliptic PDEs in heterogeneous media [70, 123] andperforated domains [65] and of parabolic PDEs in heterogeneous media [192].An alternative to GMsHDG is the HHO frameowrk for highly oscillatory elliptic prob-lems introduced in [73]. Moreover, in the context of coupled problems involving porousmedia, HHO simulations of passive transport of a solute in a fractured medium are pre-sented in [52], whereas nonlinear poroelastic phenomena in a saturated porous mediumwith a slightly compressible fluid are described in [33, 35].Extensive research has been also devoted to coupled Stokes/Darcy and Brinkman mod-els. In [51], an EHDG formulation of the Stokes/Darcy system is described. Concerningthe Brinkman model, an analysis of its HDG approximation is presented in [22], its simula-tion in the context of heterogeneous media with high-contrast is discussed in [183] and an H (div)-conforming discretisation is proposed in [142]. In [30], an HHO formulation with5ivergence-conforming Darcy velocity and higher-order Stokes velocity is devised. Hybrid formulations for inviscid Euler and laminar compressible Navier-Stokes equationsare 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 withSpalart-Allmaras turbulence model is presented in [193], whereas a large-eddy simulationframework is introduced in [133]. In addition, an entropy stable space-time discretisationwas proposed for the compressible Navier-Stokes equations using an HDG approach inspace and a discontinuous approximation in time [266]. More recently, special attentionwas dedicated to the development of positivity-preserving Riemann solvers in the contextof hybridised DG methods [263]. For a complete review on HDG methods for compress-ible flows, interested readers are referred to [263], whereas the application to gas kineticsmodelled by means of the linearised Bhatnagar-Gross-Krook equation is discussed in [254].
Computational physics community is showing increasing interest towards the applicationof hybrid discretisation methods to the simulation of magnetic plasma physics. Promisingpreliminary results concerning the HDG approximation of the Grad-Shafranov equation inaxisymmetric confinement devices modelling fusion reactors are described in [233, 234]. Inthe context of magnetohydrodynamics (MHD), an HDG method for steady-state linearisedincompressible MHD equations is proposed in [180]. Approximation strategies for theunsteady compressible MHD equations using HDG, EDG and the interior embedded DG(IEDG) methods with DIRK time integrators are explored in [74].
The shallow water equations have been extensively studied in the context of hybrid DGmethods, starting from the linearised shallow water system in [38] to the nonlinear Korteweg-de Vries equation in [63, 231]. In both the above mentioned works, time integration isperformed implicitly using a backward Euler method. Extension to high-order backwarddifferentiation formulas is discussed in [128] in the context of the Benjamin-Bona-Mahonyequation. To reduce the computational cost of fully-implicit procedures, in [228] an op-erator splitting is applied to the Green-Naghdi equation and the nonlinear hyperbolicsubproblem is solved using an explicit approach, whereas the implicit time integrator isonly applied to the linear dispersive subproblem. A similar idea is presented in [169] todevise an implicit-explicit (IMEX) HDG-DG scheme in which the linear part of the prob-lem is solved using a hybridised DG method and a singly diagonally implicit Runge-Kutta6SDIRK) scheme and the nonlinear one is approximated by means of an explicit Runge-Kutta (RK) DG discretisation. A detailed comparison of explicit and implicit approachesto the nonlinear shallow water equations is provided in [229].
The benefits of high-order methods in the simulation of wave propagation prompted exten-sive research on hybrid discretisation methods in the fields of electromagnetics, elastody-namics and acoustics. A detailed review on HDG and EDG approaches for these problemsis available in [132].Starting from the work in [203], research on time-harmonic Maxwell’s equations tackledthe analysis and development of HDG formulations [186], including methods suitable forsimulations at large wave numbers [189] and Schwarz-type domain decomposition (DD)strategies designed for HDG [18, 185]. Recent applications of HDG to time-harmonicMaxwell’s equations focus on wave propagation in heterogeneous media modelling pho-tovoltaic cells [41], coupling with nonlocal hydrodynamic Drude and generalised nonlocaloptical response models [184] and with hydrodynamic models for metals [257–259, 271] tosimulate plasmonic nanostructures. In the context of time-domain Maxwell’s equations,HDG methods are presented and analysed in [55, 59, 122], whereas implicit hybridised DGdiscretisations are proposed in [67].In the framework of elastodynamics, HDG with DIRK time integrators were introducedin [202], whereas in [256] an HDG spectral element method (HDG-SEM) is utilised tosimulate 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 ap-proach based on DIRK time integrators. Since then, explicit HDG formulations utilis-ing strong stability-preserving RK (SSPRK) and explicit RK integrators were proposedin [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 prop-agation is performed in [173]. More recently, an HDG-based cut finite element startegywith local time stepping was presented in [237]. It is worth recalling that devising aconservative numerical scheme is a critical aspect for the accurate simulation of acousticwave propagation. To correct the dissipative nature of the method analysed in [93], anenergy-conservative HDG formulation with a two-step Stormer-Numerov time-marchingis proposed in [86]. Moreover, symplectic [232] and multisymplectic [190] HDG schemespreserving the Hamiltonian structure of the PDEs under analysis are developed to achieveenergy conservation.Among the applications of HDG to wave propagation phenomena, it is also worthmentioning the degree adaptive approximation of the mild slope equation to perform har-bour simulations [152] and the cardiac electrophysiology simulations of the monodomain7odel [159, 227].
In linear elasticity, the imposition of the symmetry of the stress tensor using HDG methodsbased on mixed formulations has been extensively studied in the literature. Indeed, thefirst formulations introduced in [141,251] experienced suboptimal convergence of the mixedvariable and a loss of superconvergence of the postprocessed displacement field. To rem-edy this issue, a formulation considering a weakly symmetric stress tensor was presentedin [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 andhybrid variables; the M -decomposition framework [82,83,85] is applied to the linear elasticproblem [81] to enrich the discrete spaces of approximation utilised in the local prob-lem; an alternative formulation imposing the symmetry of the mixed variable pointwisevia Voigt notation is proposed for high-order and the lowest-order HDG discretisationsin [245] and [244], respectively. In the context of the high-order discretisation, a novelpostprocessing strategy accounting for rigid translation and rotation is also devised. It isworth noting that hybrid methods relying on primal formulations do not suffer from theseissues, see e.g. HHO [109].Timoshenko beams are discussed in [47, 48], whereas the case of Kirchhoff plates isconsidered in [162] using HDG and in [27] using HHO methods.In the context of nonlinear elasticity, the first hybrid discretisation formulation waspresented in [167]. In this work, it was observed that the method may not converge tothe exact solution if the interelement jumps are not appropriately penalised and a detailednumerical 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 largedeformations was proposed [255]. In addition, HHO discretisations of hyperelastic mate-rials in small and finite deformations were presented in [34] and [15], respectively. HHOdiscretisations for problems involving plastic and elastoplatic simulations are discussedin [16, 17], whereas contact phenomena are addressed in [66]. The first attempt to solve interface problems using hybrid discretisation techniques wasproposed in [165] using a body-fitted mesh. In this context, a superparametric HDGformulation was considered to limit the geometric error due to the polygonal approximationof curved interfaces.Recently, immersed methods have received special attention, both in the context ofHHO and HDG formulations. More precisely, unfitted HHO methods relying on a cell ag-glomeration procedure to remedy small cut instabilities are analysed for scalar and vectorial8econd-order elliptic problems in [39,40]. In the framework of HDG, Poisson interface prob-lems are treated in [121] by means of an unfitted method introducing appropriately definedansatz functions in the vicinity of the interface. An alternative approach to handle curvedinterfaces is proposed in [217] where a fictitious domain strategy is developed coupling amesh of planar faces and a transferring function for the imposition of the transmissionconditions 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 inter-faces and a cell agglomeration procedure is described for the wave equation. Similarly, theextended HDG (X-HDG) method introduces a framework in which the HDG local problemis modified only in the elements cut by the interface. In this context, cut instabilities arehandled by displacing the mesh nodes responsible for the bad cuts [154–156]. Finally, anHDG-based phase-field model for brittle fracture was recently proposed in [194].
Geometry representation plays a crucial role in the capability of high-order methods toachieve optimal accuracy. In the context of HDG, high-order isoparametric approachesin presence of curved meshes are utilised in many references, see e.g. [148, 193], whereasthis technique is addressed for HHO in [29]. An alternative approach relying on mesheswith 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-Shafranovequation [233, 234]. It is worth noting that all the techniques mentioned above introducegeometric errors due to the polynomial approximation of the boundaries. In order to exploitthe exact CAD representation of the boundaries, the NURBS-enhanced finite elementmethod (NEFEM) [241, 242] is employed in [149, 239, 247] to devise HDG formulationswith exact geometry for Stokes, linear elastic and electrostatics problems, respectively.
Hybrid discretisation methods have been traditionally developed in the context of high-order 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 frame-work, a new class of lowest-order hybrid discretisations was developed, with unknownsapproximated by means of constant functions on the mesh faces. The recently proposedface-centred finite volume (FCFV) for Poisson, Stokes [243] and linear elasticity [244] canbe interpreted as an HDG method of degree zero. Variants of this approach achievingoptimal second-order convergence of the primal variable are discussed in [150, 262]. Stem-ming from HHO, lowest-order nonconforming discretisations are proposed in [32] for linearelasticity and in [72] for elliptic obstacle problems. As their high-order counterparts, theabove mentioned methodologies allow the use of generic polygonal and polyhedral elementsand provide a workaround to the sensitivity issues of FV methods to mesh distortion and9tretching [118, 119].
Although hybrid discretisation methods are responsible for a substantial reduction of de-grees of freedom with respect to classical DG methods, their applicability to realistic prob-lems of engineering interest still rely on the development of efficient solution strategies forlarge scale systems.In [26], a DD strategy based on restricted additive Schwarz methods is proposed forhybridised DG approximations, whereas an optimised Schwarz DD approach suitable tohandle the many-subdomain case is discussed in [144]. Starting from [80], several works alsoexplored the capabilities of multigrid solvers for HDG formulations, including hierarchicalscale separation [238], geometric multigrid [265], nested geometric multigrid on many-core 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 linearelasticity [127]. Finally, iterative algorithms inspired by the Gauss-Seidel method wereproposed in [196] and tested on massively parallel architectures up to 16,384 cores. Ablock symmetric Gauss-Seidel type preconditioner was also introduced in [224], whereas amultilevel 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 meansof a posteriori estimates of the error in the primal, mixed and hybrid variables, as well asin quantities of interest.Starting from the seminal works [104, 105] establishing reliability and efficiency of er-ror estimates for the HDG approximations of second-order elliptic equations, a posteriori estimates were developed for steady and unsteady scalar convection-diffusion problemsin [57] and [182], respectively, and for the vectorial case of incompressible Oseen [23] andBrinkman [22] flows. In addition, constant-free computable a posteriori error estimatesare 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 approximationsbased on non-uniform polynomial degrees [60,61] and the superconvergence property of thepostprocessed solution [77, 89] prompted the development of degree adaptive proceduresbased on superparametric HDG methods [152,153] and on isoparametric HDG-NEFEM ap-proaches [149, 239, 247]. Degree adaptivity is also applied to the simulation of cardiac elec-trophysiology in [159]. On the other hand, mesh adaptivity procedures to capture localisedabrupt changes in the solution were devised in [234] and [194] for the Grad-Shafranovequation and the phase-field model for brittle fracture, respectively. Octree-based mesh10efinement is performed in [230] for anisotropic inhomogeneous diffusion problems. Meshadaptivity driven by local error indicators is also employed in the context of second-orderFCFV approximations [150,262]. Concerning the error in quantities of interest, an adjoint-based method allowing to achieve superconvergent approximations of linear functionals isdescribed in [103] and goal-oriented mesh adaptation strategies are proposed in [136, 268].
The accuracy of high-order HDG approximations has been recently exploited to developefficient algorithms coupling different numerical methodologies in different regions of thecomputational domain.In [172], a strategy coupling HDG and a vertex-centred finite volume method is pro-posed to simulate transient inviscid flows using coarse meshes designed for steady-stateproblems. In addition, different couplings of HDG and continuous Galerkin (CG) dis-cretisations were explored in the literature. A strategy inspired by a non-overlapping DDmethod is presented in [208] in the context of incompressible Navier-Stokes flows cou-pled with conjugate heat transfer phenomena. An alternative minimally-intrusive couplingbased on a Nitsche’s formulation of the CG method was first introduced in [175] for lin-ear elastic problems involving nearly incompressible materials and was extended to FSIproblems with weakly compressible flows in [176].
In recent years, the accuracy of the HDG method and its flexibility to devise high-orderadaptive discretisation have been also employed to devise high-fidelity reduced and surro-gate models. In [260,261], a reduced order model to accelerate the Monte-Carlo simulationof stochastic elliptic PDEs is constructed coupling a high-order HDG method with a re-duced basis and empirical interpolation approach. The combination of an HDG solverfor time-harmonic Maxwell’s equations and a proper orthogonal decomposition (POD)strategy to design parametrised plasmonic nanogap structures is proposed in [259]. AnHGD-POD reduced order model (ROM) is also discussed in [249] for the fast simulationof the unsteady heat equation. More recently, an a priori
ROM based on HDG and theproper generalised decomposition was proposed to simulate Stokes flows in geometricallyparametrised domains [147, 240].
The success of hybrid discretisation methods led to the development of targeted open-source libraries and to their implementation in existing finite element libraries available11pen-source. To the best of the authors’ knowledge, the hybridised DG method based onprimal formulations is available in the following libraries: • MFEM [14, 20] • Netgen/NGSolve [11, 235]whereas the libraries • deal.II [12, 25] • Feel++ [4, 213] • Firedrake [6, 219] • Nektar++ [10, 43]provide implementations of the HDG method based on mixed formulations. Finally, theHHO method is available in • code aster [1, 211] • Code Saturne [2, 137] • DiSk++ [3, 71] • GetFEM [7, 221] • HArDCore [8]All above mentioned libraries rely on either Fortran or C/C++ implementations, whereasopen-source libraries implementing HDG in MATLAB include: • HDG3D [9, 143] • FESTUNG [5, 166]
In this section, the formulation of the HDG method for the Poisson equation is brieflyrecalled. Special attention is devoted to the identification of the building blocks of thenumerical scheme whose implementation will be detailed in section 8. Interested readersare referred to [89] for a complete theoretical introduction to the HDG method for Poissonequation and to [246] for a tutorial on its derivation.12et Ω ⊂ R n sd be an open bounded domain in n sd spatial dimensions such that itsboundary is ∂ Ω=Γ D ∪ Γ N and Γ D ∩ Γ N = ∅ . The strong form of the Poisson equation is − ∇ · ( κ ∇ u ) = s in Ω, u = u D on Γ D , n · κ ∇ u = g on Γ N , (1)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 theboundary. Consider a partition of Ω in n el disjoint subdomains such thatΩ = n el (cid:91) e =1 Ω e , Ω i ∩ Ω j = ∅ for i (cid:54) = j and define the mesh skeleton as Γ := (cid:34) n el (cid:91) e =1 ∂ Ω e (cid:35) \ ∂ Ω . Following the HDG rationale [89, 92, 200, 201, 204, 246], a mixed variable q = −√ κ ∇ u isintroduced and problem (1) is rewritten as a system of first-order equations element-by-element, that is, q + √ κ ∇ u = in Ω e , e =1 , . . . , n el , ∇ · ( √ κ q ) = s in Ω e , e =1 , . . . , n el , u = u D on Γ D , n ·√ κ q = − g on Γ N , (cid:74) u n (cid:75) = on Γ, (cid:74) n ·√ κ q (cid:75) = 0 on Γ, (2)where the jump operator (cid:74) · (cid:75) is defined as (cid:74) (cid:12) (cid:75) := (cid:12) i + (cid:12) j , being (cid:12) i and (cid:12) j the evaluations of the quantity (cid:12) in two neighbouring elements Ω i andΩ j sharing a given interface [191]. The last two conditions in (2), known as transmission onditions , enforce the continuity of the solution and of its normal flux across the internalmesh skeleton Γ.The HDG algorithm solves equation (2) in two stages. First, an independent hybridvariable ˆ u is introduced to represent the trace of the solution on ∂ Ω e \ Γ D and the primaland mixed variables ( u e , q e ) in each element Ω e , e =1 , . . . , n el are expressed as functions ofthe unknown ˆ u , namely q e + √ κ ∇ u e = in Ω e , e =1 , . . . , n el , ∇ · ( √ κ q e ) = s in Ω e , e =1 , . . . , n el , u e = u D on ∂ Ω e ∩ Γ D , u e = ˆ u on ∂ Ω e \ Γ D . (3) Remark 1.
Equation (3) represents the n el HDG local problems. This stage correspondsto the hybridisation of the mixed problem, see [138], and is equivalent to the static con-densation procedure in classical continuous Galerkin methods [157].Second, the hybrid variable is computed by solving the HDG global problem, which ac-counts for the transmission conditions on the mesh skeleton Γ and the Neumann boundarycondition on Γ N , that is, (cid:74) u n (cid:75) = on Γ, (cid:74) n ·√ κ q (cid:75) = 0 on Γ, n ·√ κ q = − g on Γ N . (4) Remark 2.
The first condition is automatically fulfilled owing to the Dirichlet boundarycondition u e =ˆ u on ∂ Ω e \ Γ D imposed in the local problem and to the uniqueness of thehybrid 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 bysolving n el independent problems, see equation (3), element-by-element. Following the rationale introduced in [246], the discrete functional spaces V h (Ω):= { v ∈ L (Ω) : v | Ω e ∈ P p (Ω e ) ∀ Ω e , e =1 , . . . , n el } , (5a) V (cid:86) h ( S ):= { ˆ v ∈ L ( S ) : ˆ v | Γ i ∈ P p (Γ i ) ∀ Γ i ⊂ S ⊆ Γ ∪ ∂ Ω } , (5b)are defined for the approximation of the element-based and face-based variables, respec-tively. In (5), P p (Ω e ) and P p (Γ i ) stand for the spaces of polynomial functions of completedegree at most p in Ω e and on Γ i , respectively.14or e =1 , . . . , n el , the weak form of the HDG local problem is: given u D on Γ D and ˆ u h on Γ ∪ Γ N , find ( u he , q he ) ∈ V h (Ω e ) × (cid:2) V h (Ω e ) (cid:3) n sd that satisfy − ( w , q he ) Ω e +( ∇ · ( √ κ w ) , u he ) Ω e = (cid:104) n ·√ κ w , u D (cid:105) ∂ Ω e ∩ Γ D + (cid:104) n ·√ κ w , ˆ u h (cid:105) ∂ Ω e \ Γ D , (6a)( v, ∇ · ( √ κ q he )) Ω e + (cid:104) v, τ u he (cid:105) ∂ Ω e = ( v, s ) Ω e + (cid:104) v, τ u D (cid:105) ∂ Ω e ∩ Γ D + (cid:104) v, τ ˆ u h (cid:105) ∂ Ω e \ Γ D , (6b)for all ( v, w ) ∈ V h (Ω e ) × (cid:2) V h (Ω e ) (cid:3) n sd , where ( · , · ) D and (cid:104)· , ·(cid:105) S denote the L inner productson a generic subdomain D ⊂ Ω and S ⊂ Γ ∪ ∂ Ω, respectively.
Remark 3.
In equation (6), τ represents a stabilisation parameter influencing accuracy,stability and convergence of the HDG method [89, 92, 200, 201, 204].Similarly, the weak form of the HDG global problem is: find ˆ u h ∈ V (cid:86) h (Γ ∪ Γ N ) thatsatisfies n el (cid:88) e =1 (cid:110) (cid:104) ˆ v, n ·√ κ q he (cid:105) ∂ Ω e \ Γ D + (cid:104) ˆ v, τ u he (cid:105) ∂ Ω e \ Γ D −(cid:104) ˆ v, τ ˆ u h (cid:105) ∂ Ω e \ Γ D (cid:111) = − n el (cid:88) e =1 (cid:104) ˆ v, g (cid:105) ∂ Ω e ∩ Γ N , (7)for all ˆ v ∈ V (cid:86) h (Γ ∪ Γ N ). An isoparametric formulation is considered for the primal, mixed and hybrid variables inthe discrete spaces (5), that is, u h = n en (cid:88) i =1 N i u i , q h = n en (cid:88) i =1 N i q i , ˆ u h = n fn (cid:88) i =1 ˆ N i ˆu i , (8)where u i , q i and ˆu i are the nodal values of the unknowns, N i and ˆ N i are the polynomialshape functions of degree p defined in a reference element and on a reference face, respec-tively 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 linearsystem of equations (cid:20) A uu A uq A Tuq A qq (cid:21) e (cid:26) uq (cid:27) e = (cid:26) f u f q (cid:27) e + (cid:20) A u ˆ u A q ˆ u (cid:21) e ˆu e , (9)from which the following solution is computed (cid:26) uq (cid:27) e = (cid:26) z fu z fq (cid:27) e + (cid:20) Z u ˆ u Z q ˆ u (cid:21) e ˆu e , (10a)with the matrices (cid:20) Z u ˆ u Z q ˆ u (cid:21) e := (cid:20) A uu A uq A Tuq A qq (cid:21) − e (cid:20) A u ˆ u A q ˆ u (cid:21) e (10b)15nd the vectors (cid:26) z fu z fq (cid:27) e := (cid:20) A uu A uq A Tuq A qq (cid:21) − e (cid:26) f u f q (cid:27) e . (10c)The corresponding discretisation of the global problem (7) leads to n el (cid:88) e =1 (cid:110) (cid:2) A Tu ˆ u A Tq ˆ u (cid:3) e (cid:26) uq (cid:27) e + [ A ˆ u ˆ u ] e ˆu e (cid:111) = n el (cid:88) e =1 [ f ˆ u ] e . (11)By plugging the local elemental solution (10a) into (11), the HDG problem Kˆu = ˆf involving only the globally-coupled degrees of freedom is obtained, where the matrix andthe right-hand side vector are obtained by assembling the elemental contributions (cid:98) K e := (cid:2) A Tu ˆ u A Tq ˆ u (cid:3) e (cid:20) Z u ˆ u Z q ˆ u (cid:21) e +[ A ˆ u ˆ u ] e , (12a)ˆ f e := [ f ˆ u ] e − (cid:2) A Tu ˆ u A Tq ˆ u (cid:3) e (cid:26) z fu z fq (cid:27) e . (12b)The expressions of the matrices and vectors introduced in this section are detailed inappendix A. Introduce the discrete functional space V h(cid:63) (Ω):= { v ∈ L (Ω) : v | Ω e ∈ P p +1 (Ω e ) ∀ Ω e , e =1 , . . . , n el } , (13)where P p +1 (Ω e ) denotes the space of polynomial functions of complete degree at most p +1in each element Ω e .The HDG postprocess procedure allows to compute a superconvergent approximation u (cid:63) of the primal variable by solving an independent local problem in each element, namely (cid:40) − ∇ · ( κ ∇ u (cid:63) ) = ∇ · ( √ κ q e ) in Ω e , n · κ ∇ u (cid:63) = − n ·√ κ q e on ∂ Ω e , (14)with the constraint ( u (cid:63) , Ω e = ( u e , Ω e (15)on the mean value of the solution in the element.16or each element Ω e , e = 1 , . . . , n el , the weak form of the postprocess procedure is:find u h(cid:63) ∈ V h(cid:63) (Ω e ) that satisfies( ∇ v (cid:63) , κ ∇ u h(cid:63) ) Ω e = − ( ∇ v (cid:63) , √ κ q he ) Ω e , (16a)( u h(cid:63) , Ω e = ( u he , Ω e , (16b)for all v (cid:63) ∈ V h(cid:63) (Ω e ).Using an isoparametric approximation for the functions in the space V h(cid:63) (Ω), the HDGlocal postprocess gives rise to the linear system (cid:20) A (cid:63)(cid:63) a (cid:63)λ a T(cid:63)λ (cid:21) e (cid:26) u (cid:63) λ (cid:27) e = (cid:20) (cid:63)q a T(cid:63)λ (cid:21) e (cid:26) I (cid:63) u I (cid:63) n sd q (cid:27) e , (17)where the saddle-point structure of the problem follows from the imposition of the con-straint (16b) via the Lagrange multiplier λ and I (cid:63) : V h → V h(cid:63) and I (cid:63) n sd : [ V h ] n sd → [ V h(cid:63) ] n sd denote the interpolation operators from the spaces of polynomial functions of degree p tothe ones of degree p +1 for scalar and vector-valued functions.The expressions of the matrices and vectors introduced in this section are detailed inappendix A. This section presents the formulation of the HDG method for the Stokes equations, ex-tending the framework previously introduced for the Poisson equation. For the sake ofsimplicity, the present work focuses on the velocity pressure formulation of the Stokesequations. For the Cauchy stress tensor formulation, a tutorial to devise an HDG methodbased on equal order approximation for all the variables and pointwise symmetric mixedvariable is presented in [151].The open bounded domain Ω ⊂ R n sd is characterised now by a boundary partitioned inthree portions disjoint by pairs such that ∂ Ω=Γ D ∪ Γ N ∪ Γ S , where Dirichlet, Neumann andslip conditions are imposed. The strong form of the Stokes equations is − ∇ · ( ν ∇ u − p I n sd ) = s in Ω, ∇ · u = 0 in Ω, u = u D on Γ D , n · ( ν ∇ u − p I n sd ) = g on Γ N , u · D + n · (cid:0) ν ∇ u − p I n sd (cid:1) E = on Γ S , (18)where the pair ( u , p ) denotes the unknown velocity and pressure fields, ν > n is the outward unit normal to the boundary, s is the vector of the bodyforces and u D and g represent the imposed velocity and pseudo-traction on the Dirichlet17nd Neumann boundaries, respectively. On the slip boundary, matrices D and E aredefined as D :=[ n , β t , . . . , β t n sd − ] and E :=[ α n , t , . . . , t n sd − ], the tangential vectors t j , j =1 , . . . , n sd − { n , t , . . . , t n sd − } 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 compat-ibility condition on the velocity field (cid:104) u D · n , (cid:105) Γ D + (cid:104) u · n , (cid:105) ∂ Ω \ Γ D = 0 . (19) Remark 5.
In case of a purely Dirichlet boundary value problem, that is ∂ Ω=Γ D , anadditional constraint needs to be introduced to retrieve uniqueness of the pressure field. Acommon approach relies on imposing the mean value of the pressure in the domain [120,218]1 | Ω | ( p, Ω = 0 , (20)or on the boundary of the domain [88, 98, 148, 199], namely1 | ∂ Ω | (cid:104) p, (cid:105) ∂ Ω = 0 . (21) Following the rationale presented in section 3, equation (18) is rewritten element-by-element as a system of first-order equations by introducing a mixed variable L = − √ ν ∇ u ,namely L + √ ν ∇ u = in Ω e , ∇ · (cid:0) √ ν L + p I n sd (cid:1) = s in Ω e , ∇ · u = 0 in Ω e , u = u D on Γ D , n · (cid:0) √ ν L + p I n sd (cid:1) = − g on Γ N , u · D − n · (cid:0) √ ν L + p I n sd (cid:1) E = on Γ S , (cid:74) u ⊗ n (cid:75) = on Γ, (cid:74) n · (cid:0) √ ν L + p I n sd (cid:1) (cid:75) = on Γ, (22)for e =1 , . . . , n el .First, the HDG algorithm performes the hybridisation step by expressing ( u e , p e , L e ) ineach element Ω e as functions of the unknown trace of the velocity (cid:98) u on the element faces18ia the HDG local problem L e + √ ν ∇ u e = in Ω e , ∇ · (cid:0) √ ν L e + p e I n sd (cid:1) = s in Ω e , ∇ · u e = 0 in Ω e , u e = u D on ∂ Ω e ∩ Γ D , u e = (cid:98) u on ∂ Ω e \ Γ D , (23a)for e =1 , . . . , n el . Note that equation (23a) is a purely Dirichlet boundary value problem.Hence, following remark 5, the constraint1 | ∂ Ω e | (cid:104) p e , (cid:105) ∂ Ω = ρ e , for e =1 , . . . , n el , (23b)is introduced, where ρ e is an independent variable representing the mean value of thepressure on the boundary ∂ Ω e . It is worth noting that the variable ρ was not present inthe HDG approximation of the Poisson equation and its treatment in the HDGlab code willbe detailed in section 9.The HDG global problem thus accounts for the transmission and non-Dirichlet bound-ary conditions, that is (cid:74) u ⊗ n (cid:75) = on Γ, (cid:74) n · (cid:0) √ ν L + p I n sd (cid:1) (cid:75) = on Γ, n · (cid:0) √ ν L + p I n sd (cid:1) = − g on Γ N , u · D − n · (cid:0) √ ν L + p I n sd (cid:1) E = on Γ S , (24a)where, following remark 2, the first condition is automatically fulfilled. In addition, theconstraint in remark 4 is rewritten element-by-element in terms of the hybrid variable (cid:98) u leading to (cid:104) u D · n , (cid:105) ∂ Ω e ∩ Γ D + (cid:104) (cid:98) u · n , (cid:105) ∂ Ω e \ Γ D = 0 , (24b)for e =1 , . . . , n el . The corresponding weak form of the HDG local problem (23) is: given u D on Γ D and (cid:98) u h on Γ ∪ Γ N ∪ Γ S , find ( L he , u he , p he ) ∈ [ V h (Ω e )] n sd × n sd × [ V h (Ω e )] n sd ×V h (Ω e ) that satisfy − ( W , L he ) Ω e +( ∇ · ( √ ν W ) , u he ) Ω e = (cid:104) n e ·√ ν W , u D (cid:105) ∂ Ω e ∩ Γ D + (cid:104) n e ·√ ν W , (cid:98) u h (cid:105) ∂ Ω e \ Γ D , ( w , ∇ · ( √ ν L he )) Ω e +( w , ∇ p he ) Ω e + (cid:104) w , τ u he (cid:105) ∂ Ω e =( w , s ) Ω e + (cid:104) w , τ u D (cid:105) ∂ Ω e ∩ Γ D + (cid:104) w , τ (cid:98) u h (cid:105) ∂ Ω e \ Γ D , ( ∇ q, u e ) Ω e = (cid:104) q, u D · n e (cid:105) ∂ Ω e ∩ Γ D + (cid:104) q, (cid:98) u h · n e (cid:105) ∂ Ω e \ Γ D , (cid:104) p he , (cid:105) ∂ Ω e = | ∂ Ω e | ρ he , (25)19or all ( W , v , q ) ∈ [ V h (Ω e )] n sd × n sd × [ V h (Ω e )] n sd ×V h (Ω e ). It is worth noting that, differentlyfrom the Poisson case, equation (25) provides ( L he , u he , p he ) in terms of two global unknowns, (cid:98) u h and ρ h := { ρ h , . . . , ρ h n el } T .Similarly, the weak form of the HDG global problem (24) is: find (cid:98) u h ∈ [ V (cid:86) h ] n sd and ρ h ∈ R n el such that n el (cid:88) e =1 (cid:110) (cid:104) (cid:98) v , n e ·√ ν L he (cid:105) ∂ Ω e \ (Γ D ∪ Γ S ) −(cid:104) (cid:98) v , n e ·√ ν L he E (cid:105) ∂ Ω e ∩ Γ S + (cid:104) (cid:98) v , p he n e (cid:105) ∂ Ω e \ (Γ D ∪ Γ S ) −(cid:104) (cid:98) v , p he n e · E (cid:105) ∂ Ω e ∩ Γ S + (cid:104) (cid:98) v , τ u he (cid:105) ∂ Ω e \ (Γ D ∪ Γ S ) −(cid:104) (cid:98) v , τ u he · E (cid:105) ∂ Ω e ∩ Γ S −(cid:104) (cid:98) v , τ (cid:98) u h (cid:105) ∂ Ω e \ (Γ D ∪ Γ S ) + (cid:104) (cid:98) v , (cid:98) u h · ( D + τ E ) (cid:105) ∂ Ω e ∩ Γ S (cid:111) = − n el (cid:88) e =1 (cid:104) (cid:98) v , g (cid:105) ∂ Ω e ∩ Γ N , (cid:104) (cid:98) u h · n e , (cid:105) ∂ Ω e \ Γ D = −(cid:104) u D · n e , (cid:105) ∂ Ω e ∩ Γ D , for e = 1 , . . . , n el (26)for all (cid:98) v ∈ [ V (cid:86) h ] n sd . The discretisation of the local problem (25) leads to the following linear system of equations A LL A Lu TLu A uu A Tpu
00 A pu ρp Tρp e Lup ζ e = f L f u f p e + A L ˆ u A u ˆ u A p ˆ u e ˆu e + e ρ e , (27)where the constraint (23b) is imposed using the Lagrange multiplier ζ . It is worth notingthat the Lagrange multiplier is required to guarantee that equation (27) is well-posed andthe computed pressure field is unique but is not utilised in the solution of the global HDGproblem. The resulting local elemental solution is given by Lup ζ e = z fL z fu z fp z fζ e + Z L ˆ u Z u ˆ u Z p ˆ u Z ζ ˆ u e ˆu e + z ρL z ρu z ρp z ρζ e ρ e , (28a)20ith the matrix and vectors defined as Z L ˆ u Z u ˆ u Z p ˆ u Z ζ ˆ u e := A LL A Lu TLu A uu A Tpu
00 A pu ρp Tρp − e A L ˆ u A u ˆ u A p ˆ u e , (28b) z fL z fu z fp z fζ e := A LL A Lu TLu A uu A Tpu
00 A pu ρp Tρp − e f L f u f p e , (28c) z ρL z ρu z ρp z ρζ e := A LL A Lu TLu A uu A Tpu
00 A pu ρp Tρp − e e . (28d)The discrete form of the global problem (26) reads as n el (cid:88) e =1 (cid:110) (cid:2) A ˆ uL A ˆ uu A ˆ up (cid:3) e Lup e + [ A ˆ u ˆ u ] e ˆu e (cid:111) = n el (cid:88) e =1 [ f ˆ u ] e , T [ A p ˆ u ] e ˆu e = − T [ f p ] e . (29)By inserting the solution (28a) into (29), the following system involving the globally-coupled unknowns ˆu and ρ is obtained (cid:20) (cid:98) K GG T (cid:21) (cid:26) ˆu ρ (cid:27) = (cid:26) ˆf ˆ u ˆf ρ (cid:27) , (30)where the matrices and vectors have the form (cid:98) K := A n el e =1 (cid:2) A ˆ uL A ˆ uu A ˆ up (cid:3) e Z L ˆ u Z u ˆ u Z p ˆ u Z ζ ˆ u e +[ A ˆ u ˆ u ] e , (31a) G T := T [ A p ˆ u ] T [ A p ˆ u ] · · · T [ A p ˆ u ] n el , (31b) ˆf ˆ u := A n el e =1 [ f ˆ u ] e − (cid:2) A ˆ uL A ˆ uu A ˆ up (cid:3) e z fL z fu z fp z fζ e , (31c) ˆf ρ := − T [ f p ] T [ f p ] · · · T [ f p ] n el . (31d)21he expressions of the matrices and vectors introduced above are detailed in ap-pendix B. Remark 6.
Differently from the Poisson case, the global problem (30) features a saddle-point structure, as classical in the context of incompressible flows [120]. The proof ofthe symmetry of the HDG matrix in equation (30) can be devised following the rationaledescribed in [151].
The HDG postprocess procedure presented in section 3.4 for the Poisson equation can beextended straightforwardly to the case of the vectorial variable u . Remark 7.
The postprocessing procedure utilised here is inspired by the work of Sten-berg [253] and was exploited in the framework of HDG to obtain an improved approxima-tion of the velocity field via the solution of an additional inexpensive element-by-elementproblem [199, 247]. Nonetheless, in the context of incompressible flows, it is often of in-terest retrieving an H (div)-conforming and exactly divergence-free approximation of thevelocity field. For this purpose, alternative postprocessing strategies inspired by the Brezzi-Douglas-Marini (BDM) projection operator [36] were proposed [90, 92]. It is worth notingthat the above mentioned procedures are suitable only for the velocity-pressure formulationof the incompressible flow equations. In case a formulation based on the Cauchy stresstensor is considered, an additional constraint is required to handle the rigid rotationalmodes as discussed in [148, 151, 245].A superconvergent velocity field u (cid:63) is thus obtained by solving an independent localproblem in each element, namely (cid:40) − ∇ · ( ν ∇ u (cid:63) ) = ∇ · ( √ ν L e ) in Ω e , n · ν ∇ u (cid:63) = − n ·√ ν L e on ∂ Ω e , (32)with the constraint ( u (cid:63) , Ω e = ( u e , Ω e (33)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 procedureis: find u h(cid:63) ∈ [ V h(cid:63) (Ω e )] n sd such that( ∇ v (cid:63) , ν ∇ u h(cid:63) ) Ω e = − ( ∇ v (cid:63) , √ ν L he ) Ω e , (34a)( u h(cid:63) , Ω e = ( u he , Ω e , (34b)for all v (cid:63) ∈ [ V h(cid:63) (Ω e )] n sd . 22sing an isoparametric approximation for the functions in the space [ V h(cid:63) (Ω)] n sd , theHDG local postprocess for the Stokes equations leads to (cid:20) A (cid:63)(cid:63) A (cid:63)λ A T(cid:63)λ (cid:21) e (cid:26) u (cid:63) λ (cid:27) e = (cid:20) (cid:63)L A T(cid:63)λ (cid:21) e (cid:26) I (cid:63) n sd u I (cid:63) n sd × n sd L (cid:27) e , (35)where λ is the vector of Lagrange multipliers imposing the constraint (34b) on the elementalmean value and I (cid:63) n sd × n sd :[ V h ] n sd × n sd → [ V h(cid:63) ] n sd × n sd denotes the interpolation operator for tensor-valued 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 ap-pendix B. HDGlab repository
The implementation of the HDG solver for the Poisson and Stokes equations has beenmade available as an open-source software, released under the terms of the GNU GeneralPublic License version 3.0 or any later version ( ) and isfreely 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 dimensionalmeshes in the directories meshes2D and meshes3D respectively, the pre-computed referenceelements in two and three dimensions in the directory refElem and the data structurerequired to postprocess high-order HDG solutions in the directory postprocess .The directory
Poisson contains the HDG solver for the Poisson problem as describedin section 3. The directory hdgPoisson contains the core HDG library for the Poissonequation. The directory testsPoisson is used to organise the functions that describe thesetup 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 theresults are saved after an execution of the Poisson solver.The directory
Stokes contains the HDG solver for the Stokes problem as described insection 4. The structure of this directory follows the same rationale as the one correspond-ing to the Poisson solver.The directory common contains a set of functions that are common to both the Poissonand the Stokes solvers.Finally, the directory importMesh contains a library that is provided to import a meshgenerated with the open source software
Gmsh [146] in
HDGlab . The core routines toconvert 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 ofthe imported mesh is stored in the directory meshFiles . This library is described in detailin Appendix C. 23 at meshes2Dpostprocessmeshes3DrefElemimportMesh GMSHmeshFilesexamplescommonPoisson hdgPoissontestsPoissonresPoissonStokes hdgStokestestsStokesresStokes
Figure 1: Structure of the
HDGlab repository.
Three data structures are used to manage the mesh, the reference element and the faceinformation 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 detaileddescription is provided in this section. To make the developed software more accessiblevariables of type struct are used in this work rather than class types.24 .1 Mesh
The variable mesh contains the following information: • nsd : Number of spatial dimensions. • optionNodes : Type of high-order nodal distribution, being an equally-spaced (0) ora Fekete (1) nodal set. • nOfNodes : Number of nodes. • X : Array of dimension nOfNodes × nsd containing the nodal coordinates of the mesh. • nOfElements : Number of elements. • indexT : Array of dimension nOfElements × • pElem : Array of dimension 1 × nOfElements containing the degree of approximationof each element. • matElem : Array of dimension 1 × nOfElements containing the material flag for eachelement. • nOfIntFaces : Number of interior faces (i.e. faces not on the boundary). • intFaces : Array of dimension nOfIntFaces ×
5. The first two columns contain thefirst element sharing this face and its local face number. The next two columnscontain the second element sharing this face and its local face number. The lastcolumn contains the local node number of the second face that matches with the firstlocal 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 theelement sharing this face and its local face number. The third column contains theboundary condition flag and the last column contains the flag of the boundary curveor 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 standardCG formulation, where only the field extFaces is needed to impose the boundary condi-tions. The last column of intFaces is needed to account for the different orientation ofan interior face as seen from the element on the left and on the right of a face. It is worthnoting that in two dimensions the information could be omitted because the local node25 Figure 2: Interface between two triangular elements, showing the orientation of each ele-ment and the local node number of each node on the face as seen from the element on theleft and on the right of the interface. Figure 3: Interface between two tetrahedral elements, showing the orientation of each faceand the local node number of each node on the face as seen from the element on the leftand on the right of the interface.number of the second face that matches with the first local node of the first face is alwaysequal to two, as illustrated in figure 2. However, in three dimensions this information isrequired as there are three possible rotations of the local face nodes that do not alter theorientation of the element/face, as depicted in figure 3.To illustrate the mesh data structure, a coarse mesh of the domain Ω = [0 , , with fourtriangular elements is considered. Figure 4 represents the mesh with the global numberingof 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 4is 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 functionshown in figure 7. This means that the mesh nodes in the array X are duplicated. Thisfunctionality is meant to help the handling of meshes with a non-uniform degree of approx-imation and to provide a seamless way to partition the mesh in case future users would liketo parallelise the code. If desired, a different array can be introduced for the connectivityinformation and the user only needs to redefine the function getElemConnectivity .26
23 456 78 91011 12 1 234 12 34 124 3 567 8 59 107 1314 6121181516
Figure 4: Mesh of Ω = [0 , with four triangular elements with the global numbering ofthe mesh elements (in green), the element nodes (in black), the mesh faces (in blue) andthe face nodes (in red). getElemConnectivity function Te = getElemConnectivity(mesh, iElem)Te = mesh.indexT(iElem,1):mesh.indexT(iElem,2); Mesh example >> disp(mesh)nsd: 2optionNodes: 0nOfNodes: 12X: [12x2 double]nOfElements: 4indexT: [4x2 double]pElem: [1 1 1 1]matElem: [1 1 1 1]nOfIntFaces: 4intFaces: [4x5 double]nOfExtFaces: 4extFaces: [4x4 double] Figure 5: Overview of the data contained in the mesh data structure for the example offigure 4.
As usual in an isoparametric finite element context, the information related to the approx-imation and numerical integration in an element is stored by means of a reference element,with local coordinates, ξ = ( ξ , . . . , ξ n sd ). To easily handle meshes with a non-uniform de-gree 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 eachcomponent of the refElem array, the following information is stored: • nsd : Number of spatial dimensions. • optionNodes : Type of high-order nodal distribution, being an equally-spaced (0) ora Fekete (1) nodal set. • p : Degree of approximation. • nOfVertices : Number of element vertices.27 > disp(mesh.X)0.0 0.01.0 0.00.5 0.51.0 0.01.0 1.00.5 0.51.0 1.00.0 1.00.5 0.50.0 1.00.0 0.00.5 0.5 >> disp(mesh.indexT)1 34 67 910 12 >> disp(mesh.intFaces)1 2 2 3 21 3 4 2 22 2 3 3 23 2 4 3 2 >> disp(mesh.extFaces)1 1 1 12 1 1 23 1 1 34 1 1 4 Figure 6: Detail of the fields X , indexT , intFaces and extFaces , corresponding to thedata structure mesh of figure 5. getElemConnectivity function Te = getElemConnectivity(mesh, iElem)Te = mesh.indexT(iElem,1):mesh.indexT(iElem,2); Mesh example >> disp(mesh)nsd: 2optionNodes: 0nOfNodes: 12X: [12x2 double]nOfElements: 4indexT: [4x2 double]pElem: [1 1 1 1]matElem: [1 1 1 1]nOfIntFaces: 4intFaces: [4x5 double]nOfExtFaces: 4extFaces: [4x4 double] Figure 7: Function to retrieve the connectivity of an element. • nOfNodes : Number of element nodes. • coordinates : Array of dimension nOfNodes × nsd containing the local nodal coor-dinates. • nOfFaces : Number of element faces. • face : Array of dimension 1 × nOfFaces . Each position is a structure containing thefollowing information about an element face: – nodes : Array containing the element local number of the nodes in the currentface. – nodesPerm : Array containing nsd + 1 permutations of the face nodes in the field nodes , using the element local number. Each row provides the permutation asrequired by the field intFaces in the mesh data structure. – nodesPermHDG : Array containing nsd + 1 permutations of the face nodes in thefield nodes , using the face local number. Each row provides the permutation asrequired by the field intFaces in the mesh data structure. • nOfGauss : Number of integration points.28
65 43 2
Figure 8: Reference triangular element for p = 2. • gaussWeights : Array of dimension nOfGauss ×
1, containing the integration weights. • shapeFunctions : Array of dimension nOfGauss × nOfNodes × ( nsd + 1). When thethird index is equal to 1, the array contains the value of the shape functions, for allthe nodes at all integration points. When the third index is equal to r >
1, the arraycontains the value of the derivative of the shape functions in the ξ k direction, for allthe nodes at all integration points. • shapeFunctionsNodesPPp1 : Array of dimension nOfNodes (cid:63) × nOfNodes , where nOfNodes (cid:63) denotes the number of nodes of an element with degree of approxima-tion p +1. It contains the value of the shape functions of the current element at thenodes of an element with one extra degree of approximation. This array is only usedto perform the HDG postprocess described in sections 3.4 and 4.4. • shapeFunctionsGaussPPp1 : Array of dimension nOfGauss (cid:63) × nOfNodes , where nOfGauss (cid:63) denotes the number of integration points of an element with degree ofapproximation p +1. It contains the value of the shape functions of the current ele-ment and its derivatives at the integration points of an element with one extra degreeof approximation. This array is only used to perform the HDG postprocess describedin sections 3.4 and 4.4.To illustrate the refElem data structure, figure 8 depicts the reference quadratic tri-angular element. An overview of the data contained in the refElem data structure for theexample of figure 8 is shown in figure 9.Similarly, figure 10 shows the reference quadratic tetrahedral element and figure 11depicts 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 poly-nomial approximation. However, it is straightforward for future users to change the basisused for the approximation by simply changing the reference element information. Withminimum effort it is also possible to incorporate other element types.29 > disp(mesh.X)0.0 0.01.0 0.00.5 0.51.0 0.01.0 1.00.5 0.51.0 1.00.0 1.00.5 0.50.0 1.00.0 0.00.5 0.5 >> disp(mesh.indexT)1 34 67 910 12 >> disp(mesh.intFaces)1 2 2 3 21 3 4 2 22 2 3 3 23 2 4 3 2 >> disp(mesh.extFaces)1 1 1 12 1 1 23 1 1 34 1 1 4 Reference element >> disp(refElem(2))nsd: 2optionNodes: 0p: 2nOfVertices: 3nOfNodes: 6coordinates: [6x2 double]nOfFaces: 3face: [1x3 struct]nOfGauss: 7gaussWeights: [7x1 double]shapeFunctions: [7x6x3 double]shapeFunctionsNodesPPp1: [10x6 double]shapeFunctionsGaussPPp1: [15x6x3 double] Figure 9: Overview of the data contained in the refElem data structure for the referencetriangular quadratic element of figure 8.Figure 10: Reference tetrahedral element for p = 2. Reference element2D >> disp(refElem(2))nsd: 2optionNodes: 0p: 2nOfVertices: 3nOfNodes: 6coordinates: [6x2 double]nOfFaces: 3face: [1x3 struct]nOfGauss: 7gaussWeights: [7x1 double]shapeFunctions: [7x6x3 double]shapeFunctionsNodesPPp1: [10x6 double]shapeFunctionsGaussPPp1: [15x6x3 double] Reference element 3D >> disp(refElem(2))nsd: 3optionNodes: 0p: 2nOfVertices: 4nOfNodes: 10coordinates: [10x3 double]nOfFaces: 4face: [1x4 struct]nOfGauss: 14gaussWeights: [14x1 double]shapeFunctions: [14x10x4 double]shapeFunctionsNodesPPp1: [20x10 double]shapeFunctionsGaussPPp1: [35x10x4 double] Figure 11: Overview of the data contained in the refElem data structure for the referencetetrahedral quadratic element of figure 10.
The information related to the approximation and numerical integration on a face is storedby means of a reference face, with local coordinates, η = ( η , . . . , η n sd − ). To easily handlemeshes with a non-uniform degree of approximation, the variable refFace is considered anarray of dimension p max × p max , where p max is the maximum degree of approximation used inall the elements. For each diagonal component of the refFace array, the information storedis a subset of the information stored in the refElem data structure. As an example, figure 12offers an overview of the data contained in the diagonal term of refFace corresponding to30 eference element 2D >> disp(refElem(2))nsd: 2optionNodes: 0p: 2nOfVertices: 3nOfNodes: 6coordinates: [6x2 double]nOfFaces: 3face: [1x3 struct]nOfGauss: 7gaussWeights: [7x1 double]shapeFunctions: [7x6x3 double]shapeFunctionsNodesPPp1: [10x6 double]shapeFunctionsGaussPPp1: [15x6x3 double] Reference element 3D >> disp(refElem(2))nsd: 3optionNodes: 0p: 2nOfVertices: 4nOfNodes: 10coordinates: [10x3 double]nOfFaces: 4face: [1x4 struct]nOfGauss: 14gaussWeights: [14x1 double]shapeFunctions: [14x10x4 double]shapeFunctionsNodesPPp1: [20x10 double]shapeFunctionsGaussPPp1: [35x10x4 double] Reference face 2D >> disp(refFace(2,2))nsd: 1optionNodes: 0p: 2nOfVertices: 2nOfNodes: 3coordinates: [-1 0 1]nOfGauss: 3gaussPoints: [3x1 double]gaussWeights: [3x1 double]shapeFunctions: [3x3x3 double] Figure 12: Overview of the data contained in the refFace data structure for the referenceface 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 informationassociated with the field shapeFunctions . This is only required when a mesh with a non-uniform 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 withdegree of approximation s . This is required when computing the integrals in an interiorface where the degree of approximation of the elements sharing this face is different. Itis worth noting that only the entries in the upper triangle of the array refFace containrelevant information because the degree of approximation used for the hybrid variable ina given face is defined as the maximum between the degree of approximation of the twoelements sharing the face. This section describes the preprocessing stage of the HDG solver. The implementation canbe found in the function hdgPreprocess , which produces as an output an updated versionof the data structures mesh and hdg .The data structure mesh is taken as an input, containing the fields described in sec-tion 6.1, and a new field, called indexTf is added. This field contains an array of dimension nOfFaces × nOfFaces isthe total number of mesh faces (i.e. interior faces plus exterior faces). The first columncontains the first degree of freedom of a face and the second column contains the lastdegree of freedom of the face. Figure 13 shows the data contained in indexTf , after thepreprocessing is performed, for the mesh of figure 4 and for a Poisson problem (i.e. scalarunknown). 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 hybridvariable contains an approximation of the trace of the solution, which is a scalar quantity. Incontrast, for the Stokes problem the global vector of unknowns contains an approximationof the trace of the velocity and an approximation of the mean pressure. Therefore, in two31 eference element 2D >> disp(refElem(2))nsd: 2optionNodes: 0p: 2nOfVertices: 3nOfNodes: 6coordinates: [6x2 double]nOfFaces: 3face: [1x3 struct]nOfGauss: 7gaussWeights: [7x1 double]shapeFunctions: [7x6x3 double]shapeFunctionsNodesPPp1: [10x6 double]shapeFunctionsGaussPPp1: [15x6x3 double] Reference element 3D >> disp(refElem(2))nsd: 3optionNodes: 0p: 2nOfVertices: 4nOfNodes: 10coordinates: [10x3 double]nOfFaces: 4face: [1x4 struct]nOfGauss: 14gaussWeights: [14x1 double]shapeFunctions: [14x10x4 double]shapeFunctionsNodesPPp1: [20x10 double]shapeFunctionsGaussPPp1: [35x10x4 double] Reference face 2D >> disp(refFace(2,2))nsd: 1optionNodes: 0p: 2nOfVertices: 2nOfNodes: 3coordinates: [-1 0 1]nOfGauss: 3gaussPoints: [3x1 double]gaussWeights: [3x1 double]shapeFunctions: [3x3x3 double] >> disp(mesh.indexTf)1 23 45 67 89 1011 1213 1415 16 >> disp(hdg)tau: 1problem: 'Poisson'faceInfo: [1x4 struct]nDOFglobal: 16vDOFtoSolve: [1 2 3 4 5 6 7 8] >> disp(hdg.faceInfo(1))local2global: [5 1 2]localNumFlux: [1 0 0]permutations: [0 0 0]pHat: [1 1 1] >> disp(hdg.faceInfo(2))local2global: [6 3 1]localNumFlux: [1 0 0]permutations: [0 0 2]pHat: [1 1 1] Figure 13: Data contained in mesh.indexTf data structure for the solution of the Poissonequation in the mesh of figure 4.
Reference element 2D >> disp(refElem(2))nsd: 2optionNodes: 0p: 2nOfVertices: 3nOfNodes: 6coordinates: [6x2 double]nOfFaces: 3face: [1x3 struct]nOfGauss: 7gaussWeights: [7x1 double]shapeFunctions: [7x6x3 double]shapeFunctionsNodesPPp1: [10x6 double]shapeFunctionsGaussPPp1: [15x6x3 double] Reference element 3D >> disp(refElem(2))nsd: 3optionNodes: 0p: 2nOfVertices: 4nOfNodes: 10coordinates: [10x3 double]nOfFaces: 4face: [1x4 struct]nOfGauss: 14gaussWeights: [14x1 double]shapeFunctions: [14x10x4 double]shapeFunctionsNodesPPp1: [20x10 double]shapeFunctionsGaussPPp1: [35x10x4 double] Reference face 2D >> disp(refFace(2,2))nsd: 1optionNodes: 0p: 2nOfVertices: 2nOfNodes: 3coordinates: [-1 0 1]nOfGauss: 3gaussPoints: [3x1 double]gaussWeights: [3x1 double]shapeFunctions: [3x3x3 double] indexTf Poisson >> disp(mesh.indexTf)1 23 45 67 89 1011 1213 1415 16 indexTf Stokes >> disp(mesh.indexTf)1 45 89 1213 1617 2021 2425 2829 32 >> disp(hdg)tau: 1problem: 'Poisson'faceInfo: [1x4 struct]nDOFglobal: 16vDOFtoSolve: [1 2 3 4 5 6 7 8] >> disp(hdg.faceInfo(1))local2global: [5 1 2]localNumFlux: [1 0 0]permutations: [0 0 0]pHat: [1 1 1] >> disp(hdg.faceInfo(2))local2global: [6 3 1]localNumFlux: [1 0 0]permutations: [0 0 2]pHat: [1 1 1] Figure 14: Data contained in mesh.indexTf data structure for the solution of the Stokesequation in the mesh of figure 4.dimensions each face contains 2( p + 1) degrees of freedom for the velocity and one extradegree of freedom per element is introduced for the mean pressure.The hdg data structure is also an input of the function hdgPreprocess , containing twouser defined parameters, namely: • tau : Value of the constant stabilisation parameter. • problem : String containing the name of the problem to be solved. The code providedcontains 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 informationof the mesh: – local2global : Array of dimension 1 × nOfElementFaces , containing the globalface number of the faces of the current element. – localNumFlux : Array of dimension 1 × nOfElementFaces , containing a flag forthe 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 numbercorresponds to the boundary condition to be imposed and it is linked to thethird column of the array extFaces of the mesh data structure described insection 6.1. 32 eference element 2D >> disp(refElem(2))nsd: 2optionNodes: 0p: 2nOfVertices: 3nOfNodes: 6coordinates: [6x2 double]nOfFaces: 3face: [1x3 struct]nOfGauss: 7gaussWeights: [7x1 double]shapeFunctions: [7x6x3 double]shapeFunctionsNodesPPp1: [10x6 double]shapeFunctionsGaussPPp1: [15x6x3 double] Reference element 3D >> disp(refElem(2))nsd: 3optionNodes: 0p: 2nOfVertices: 4nOfNodes: 10coordinates: [10x3 double]nOfFaces: 4face: [1x4 struct]nOfGauss: 14gaussWeights: [14x1 double]shapeFunctions: [14x10x4 double]shapeFunctionsNodesPPp1: [20x10 double]shapeFunctionsGaussPPp1: [35x10x4 double] Reference face 2D >> disp(refFace(2,2))nsd: 1optionNodes: 0p: 2nOfVertices: 2nOfNodes: 3coordinates: [-1 0 1]nOfGauss: 3gaussPoints: [3x1 double]gaussWeights: [3x1 double]shapeFunctions: [3x3x3 double] >> disp(mesh.indexTf)1 23 45 67 89 1011 1213 1415 16 >> disp(hdg)tau: 1problem: 'Poisson'faceInfo: [1x4 struct]nDOFglobal: 16vDOFtoSolve: [1 2 3 4 5 6 7 8] >> disp(hdg.faceInfo(1))local2global: [5 1 2]localNumFlux: [1 0 0]permutations: [0 0 0]pHat: [1 1 1] >> disp(hdg.faceInfo(2))local2global: [6 3 1]localNumFlux: [1 0 0]permutations: [0 0 2]pHat: [1 1 1] Figure 15: Data contained in the hdg data structure for the mesh of figure 4. – permutations : Array of dimension 1 × nOfElementFaces , containing a flag forthe permutation required to ensure that the ordering of the nodes in the globalface matches the ordering of the face nodes in the current element. – pHat : Array of dimension 1 × nOfElementFaces , containing the degree of ap-proximation used for the hybrid variable in the corresponding faces of the currentelement. • nDOFglobal : Number of global degrees of freedom. • vDOFtoSolve : Number of global degrees of freedom associated with nodes not on aDirichlet 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 Dirich-let boundary conditions. • columnsGlobalSystem : Number of columns in the global system, corresponding tothe number of unknowns given by the hybrid velocity and the mean pressure. • rowsGlobalSystem : Number of rows in the global system, including the constraintfor 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 infigure 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 infor-mation. The current implementation uses this data structure to carry the following infor-mation: • nOfMat : Number of materials in the domain. • charLength : A characteristic length, used to define the value of the stabilisationparameter of the HDG formulation. 33 eference element 2D >> disp(refElem(2))nsd: 2optionNodes: 0p: 2nOfVertices: 3nOfNodes: 6coordinates: [6x2 double]nOfFaces: 3face: [1x3 struct]nOfGauss: 7gaussWeights: [7x1 double]shapeFunctions: [7x6x3 double]shapeFunctionsNodesPPp1: [10x6 double]shapeFunctionsGaussPPp1: [15x6x3 double] Reference element 3D >> disp(refElem(2))nsd: 3optionNodes: 0p: 2nOfVertices: 4nOfNodes: 10coordinates: [10x3 double]nOfFaces: 4face: [1x4 struct]nOfGauss: 14gaussWeights: [14x1 double]shapeFunctions: [14x10x4 double]shapeFunctionsNodesPPp1: [20x10 double]shapeFunctionsGaussPPp1: [35x10x4 double] Reference face 2D >> disp(refFace(2,2))nsd: 1optionNodes: 0p: 2nOfVertices: 2nOfNodes: 3coordinates: [-1 0 1]nOfGauss: 3gaussPoints: [3x1 double]gaussWeights: [3x1 double]shapeFunctions: [3x3x3 double] >> disp(mesh.indexTf)1 23 45 67 89 1011 1213 1415 16 >> disp(hdg)tau: 1problem: 'Poisson'faceInfo: [1x4 struct]nDOFglobal: 16vDOFtoSolve: [1 2 3 4 5 6 7 8] >> disp(hdg.faceInfo(1))local2global: [5 1 2]localNumFlux: [1 0 0]permutations: [0 0 0]pHat: [1 1 1] >> disp(hdg.faceInfo(2))local2global: [6 3 1]localNumFlux: [1 0 0]permutations: [0 0 2]pHat: [1 1 1] Figure 16: Data contained in the hdg.faceInfo for the first two elements of the mesh offigure 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 thePoisson problem: • conductivity : Array of dimension 1 × nOfMat that contains the conductivity of eachmaterial 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 boundarycondition. • betaSlip : A scalar value describing the friction coefficient for the slip boundarycondition.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 facewhere a Dirichlet boundary condition is imposed. • iBC Neumann : A flag for the numerical flux function to be used on an exterior facewhere a Neumann boundary condition is imposed. • iBC Slip : A flag for the numerical flux function to be used on an exterior face wherea 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 bespecified 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.34 olve local problem Basic definitions
MeshReference element and face Loop on element facesCompute face integrals: A A A A
Depending on the boundary (vectorised loop on integration points)
Compute element integrals: A , A , f (vectorised loop on integration points) M − Compute and store , , , Element loop finished?YES NONO uu qq uq f u f q Assemble and solve global problem
Preprocess
Face connectivityLocal to global face mapping u Loop on elements u ˆ u ˆ u ˆ u Face loop finished?Loop on elements NOYESElement loop finished?Postprocess solutionLoop on elements NOYESElement loop finished?Compute ErrorVisualisation z fu z fq Z ˆ uu Z ˆ uq qq f u ˆ Figure 17: Code workflow diagram.
HDGlab
Poisson solver
A code workflow diagram of the
HDGlab
Poisson code is shown in figure 17. This sectionfocuses on the core part of the code that involves the assembly and solution of the globalsystem of equations, the element-by-element solution of the local problems and the localpostprocess to obtain a superconvergent solution.
The HDG global system of equations is assembled and solved in the hdg Poisson GlobalSystem function. For each element, hdg Poisson ElementalMatrices contains two parts corre-35
Element computation -----------------------------------------------------[N, dNx, wXY, Xg] = gaussElemCartesianInfo(refElem(pElem), Xe);Nw = bsxfun(@times, N, wXY');minusNN = -Nw*N';for iNsd = 1:nsdvElem = iNsd:nsd:ndofQ;Auq(:,vElem) = sqrt(kappa)*Nw*dNx(:,:,iNsd)';Aqq(vElem,vElem) = minusNN;endsourceTerm = poisson source(Xg,problemParams,matElem,nsd,problemParams.example);fu = Nw*sourceTerm; Figure 18: Extract of the function hdg Poisson ElementalMatrices that computes theelement integrals of the HDG formulation for the Poisson problem.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 theloop 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 leadto the matrices A uu , A u ˆ u , A q ˆ u and A ˆ u ˆ u . This computation distinguishes between interiorand 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 updatedand the vector f q is computed. For a Neumann face, the vector f ˆ u is computed.One of the distinctive parts of the HDG formulation is found in the loop on faces, wherea vector called indexFlip is used to ensure that the ordering of the degrees of freedomcorresponding to the hybrid variable, as seen from the current element, matches the globalordering 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 contributionsto the global system are prepared to be assembled, namely (cid:98) K e and ˆf e , as shown in theextract of figure 19. It is worth noting that at this stage the matrices Z u ˆ u and Z q ˆ u andthe vectors z fu and z fq , defined in equation (10), are stored in the data structure local inorder to be used during the solution of the element-by-element local problems. For largeproblems, the user might choose to save these matrices to disk before solving the globalsystem of equations. After the global system of equations is solved, the function hdg Poisson LocalProblem ,shown in figure 20, is called to solve the local problems. This function is just the imple-mentation of equation (10a). The only aspect that requires attention is the indexing of36
Elemental mapping -------------------------Z = [Auu Auq; Auq' Aqq] \ [Aul fu; Aql fq];vU = 1:ndofU;vQ = ndofU+1:ndofU+ndofQ;Zul = Z(vU,1:ndofUHat);zuf = Z(vU,ndofUHat+1);Zql = Z(vQ,1:ndofUHat);zqf = Z(vQ,ndofUHat+1);% Flipping due to the different numbering of% internal face nodes when seen from left or% right elementZql = Zql(:,indexFlip);Zul = Zul(:,indexFlip);Alq = Aql(:,indexFlip)';Alu = Aul(:,indexFlip)';All = All(indexFlip,indexFlip);fl = fl(indexFlip);% Elemental matricesKe = Alq*Zql + Alu*Zul + All;fe = fl - (Alq*zqf + Alu*zuf); Figure 19: Extract of the hdg Poisson ElementalMatrices function that computes theelemental matrix (cid:98) K e and vector ˆf e .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 com-ponents. It is also worth noting that this function accounts for the possibility to have anon-uniform degree of approximation. As discussed in section 3.4, once the primal and mixed variables are computed, it is possi-ble to perform a local, element-by-element, postprocess procedure to obtain a more accu-rate, 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 theelemental 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 infor-mation is known at this stage. Therefore, to compute the high-order nodal distribution ofdegree p +1, the nodal distribution in the reference element is mapped to the physical spaceby using the isoparametric mapping of degree p . This implies that, for curved elements,a subparametric formulation is considered. This formulation can lead to a suboptimalrate of convergence for the postprocessed solution as demonstrated in [239, 247], where aNURBS-enhanced implementation was proposed.37 unction [u,q] = ...hdg Poisson LocalProblem(mesh,refElem, ...refFace,hdg,uHat,local,nOfComponents)% Count the number of unknownsnOfDOFU = 0;for iElem = 1:mesh.nOfElementspElem = mesh.pElem(iElem);nOfElementNodes = refElem(pElem).nOfNodes;nOfDOFU = nOfDOFU + ...nOfElementNodes*nOfComponents;end% Initialisationu = zeros(nOfDOFU,1);q = zeros(nOfDOFU*mesh.nsd,1);indexUIni = 1;indexQIni = 1;for iElem = 1:mesh.nOfElementspElem = mesh.pElem(iElem);iGlobalFace = ...hdgElemToFaceIndex(mesh.indexTf, ...refFace,hdg.faceInfo(iElem), ...refElem(pElem).nOfFaces, ...nOfComponents);nOfElementNodes = refElem(pElem).nOfNodes;nDOFsElemU = nOfElementNodes*nOfComponents;nDOFsElemQ = nDOFsElemU*mesh.nsd;indexUEnd = indexUIni + nDOFsElemU - 1;indexQEnd = indexQIni + nDOFsElemQ - 1;u(indexUIni:indexUEnd) = ...local(iElem).Zul*uHat(iGlobalFace) ...+ local(iElem).zuf;q(indexQIni:indexQEnd) = ...local(iElem).Zql*uHat(iGlobalFace) ...+ local(iElem).zqf;indexUIni = indexUEnd + 1;indexQIni = indexQEnd + 1;end Figure 20: Function hdg Poisson LocalProblem to solve the element-by-element localproblems.
HDGlab
Stokes solver
In this section, the
HDGlab solver for the Stokes equations is presented. It is worth notingthat the code features the same structure introduced in the previous section for the Poissoncase. Hence, only the differences with respect to the Poisson solver will be detailed.
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.38 unction uStar = hdg Poisson LocalPostprocess(mesh,refElem,u,q,problemParams,nOfComponents)% InitialisationnOfDOFUStar = 0;for iElem = 1:mesh.nOfElementspElem = mesh.pElem(iElem);pElemStar = pElem + 1;nOfElementNodesStar = refElem(pElemStar).nOfNodes;nOfDOFUStar = nOfDOFUStar + nOfElementNodesStar*nOfComponents;enduStar = zeros(nOfDOFUStar,1);% ComputationindexIni = 1;for iElem = 1:mesh.nOfElementsTe = getElemConnectivity(mesh, iElem);Xe = mesh.X(Te,:);pElem = mesh.pElem(iElem);pElemStar = pElem + 1;NXeStar = refElem(pElem).shapeFunctionsNodesPPp1;XeStar = NXeStar*Xe;iMat = mesh.matElem(iElem);kappa = problemParams.conductivity(iMat);[ug,qg] = hdg Poisson LocalPostprocessInterpolation(NXeStar,u,q,refElem(pElemStar). ...nOfNodes,Te,mesh.nsd);[Ke,Bqe,intUStar,intU] = ...hdg Poisson LocalPostprocessElemMat(refElem(pElemStar),XeStar,ug,qg,kappa);% Constraint with Lagrange multipliersK = [Ke intUStar; intUStar' 0];f = [Bqe; intU];nOfDOFUStar = refElem(pElemStar).nOfNodes*nOfComponents;indexEnd = indexIni + nOfDOFUStar - 1;% Elemental solutionsol = K \ f;uStar(indexIni:indexEnd) = sol(1:end-1);% IndexingindexIni = indexEnd + 1;end Figure 21: Function hdg Poisson LocalPostprocess to perform a local, element-by-element, postprocess of the solution.The first difference with respect to the Poisson code is represented by the vectorialnature of the primal and hybrid variables and by the tensor-valued mixed variable. For thesake of computational efficiency, the representation of the second-order tensor L is writtenas a vector by rows, namely L = (cid:40) [L L L L ] T , in 2D , [L L L L L L L L L ] T , in 3D . Figure 22 reports the computation of the elemental matrices A LL , A Lu and A pu and theelemental vector f u . It is worth noting that the assembly of these matrices and this vectoraccount for the appropriate numbering of the vector-valued and tensor-valued unknowns.Similarly to the Poisson case, the loop on integration points is vectorised via the command bsxfun . 39 Initialisationmat.i = zeros(1, hdg.rowsGlobalSystem);mat.j = zeros(1, hdg.columnsGlobalSystem);mat.Kij = zeros(1, hdg.nDOFglobal);f = zeros(hdg.rowsGlobalSystem, 1);% Element computation -------------------------[N, dNx, wXY, Xg] = ...gaussElemCartesianInfo(refElem(pElem), Xe);Nw = bsxfun(@times, N, wXY');NwN = Nw*N';for iNsd = 1:nsd2vElem = iNsd:nsd2:ndofL;ALL(vElem,vElem) = -NwN;endsourceTerm = ...stokes source(Xg,problemParams,matElem, ...nsd,problemParams.example);for iNsd = 1:nsdwElem = iNsd:nsd:ndofU;Apu(:,wElem) = dNx(:,:,iNsd)*Nw';fu(wElem) = Nw*sourceTerm(:,iNsd);for jNsd = 1:nsdzElem = (iNsd-1)*nsd+jNsd:nsd2:ndofL;ALu(zElem,wElem) = ...sqrt(nu)*dNx(:,:,jNsd)*Nw';endend Figure 22: Extract of the hdg Stokes ElementalMatrices function that computes theelement integrals of the HDG formulation for the Stokes problem.
In the second part of hdg Stokes ElementalMatrices , the face integrals are computedwithin a loop on the element faces. More precisely, the matrices A L ˆ u , A uu , A u ˆ u and A p ˆ u are computed for the local problem, whereas the matrix A ˆ u ˆ u is computed for the globalproblem. In addition, on the Neumann faces the vector f ˆ u is computed, whereas on theDirichlet 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 matricesinvolved in the global problem are such that A ˆ uL = A TL ˆ u , (36a) A ˆ uu = A Tu ˆ u , (36b) A ˆ up = A Tp ˆ u . (36c)In case slip boundary conditions are also considered, the properties in equation (36)no longer hold. In this context, the matrices A ˆ uL , A ˆ uu and A ˆ up are computed in thefunction hdg Stokes ElementalMatrices . Figure 23 displays the initialisation of theabove elemental matrices which are then computed in the loop on faces when the flagin hdg.faceInfo.localNumFlux matches ctt.iBC Slip .A specific treatment of the case in which slip boundary conditions are considered is alsorequired for the definition of the elemental matrices of the global problem with appropriateordering of the degrees of freedom of the hybrid variable using the indexFlip vector, seefigure 24. 40 Initialisationmat.i = zeros(1, hdg.rowsGlobalSystem);mat.j = zeros(1, hdg.columnsGlobalSystem);mat.Kij = zeros(1, hdg.nDOFglobal);f = zeros(hdg.rowsGlobalSystem, 1);% Element computation -------------------------[N, dNx, wXY, Xg] = ...gaussElemCartesianInfo(refElem(pElem), Xe);Nw = bsxfun(@times, N, wXY');NwN = Nw*N';for iNsd = 1:nsd2vElem = iNsd:nsd2:ndofL;ALL(vElem,vElem) = -NwN;endsourceTerm = ...stokes source(Xg,problemParams,matElem, ...nsd,problemParams.example);for iNsd = 1:nsdwElem = iNsd:nsd:ndofU;Apu(:,wElem) = dNx(:,:,iNsd)*Nw';fu(wElem) = Nw*sourceTerm(:,iNsd);for jNsd = 1:nsdzElem = (iNsd-1)*nsd+jNsd:nsd2:ndofL;ALu(zElem,wElem) = ...sqrt(nu)*dNx(:,:,jNsd)*Nw';endend% Additional contributions if any of the ...faces features a slip BCif isAnySlipFaceAlL = zeros(ndofUhat,ndofL);Alu = zeros(ndofUhat,ndofU);Alp = zeros(ndofUhat,ndofP);end Figure 23: Extract of the hdg Stokes ElementalMatrices function that initialises theelemental matrices associated with the slip boundaries in the global problem for the Stokesproblem. % Initialisationmat.i = zeros(1, hdg.rowsGlobalSystem);mat.j = zeros(1, hdg.columnsGlobalSystem);mat.Kij = zeros(1, hdg.nDOFglobal);f = zeros(hdg.rowsGlobalSystem, 1);% Element computation -------------------------[N, dNx, wXY, Xg] = ...gaussElemCartesianInfo(refElem(pElem), Xe);Nw = bsxfun(@times, N, wXY');NwN = Nw*N';for iNsd = 1:nsd2vElem = iNsd:nsd2:ndofL;ALL(vElem,vElem) = -NwN;endsourceTerm = ...stokes source(Xg,problemParams,matElem, ...nsd,problemParams.example);for iNsd = 1:nsdwElem = iNsd:nsd:ndofU;Apu(:,wElem) = dNx(:,:,iNsd)*Nw';fu(wElem) = Nw*sourceTerm(:,iNsd);for jNsd = 1:nsdzElem = (iNsd-1)*nsd+jNsd:nsd2:ndofL;ALu(zElem,wElem) = ...sqrt(nu)*dNx(:,:,jNsd)*Nw';endend% Additional contributions if any of the ...faces features a slip BCif isAnySlipFaceAlL = zeros(ndofUhat,ndofL);Alu = zeros(ndofUhat,ndofU);Alp = zeros(ndofUhat,ndofP);end % Flipping due to the different numbering ...of internal face nodes when seen from ...left or right elementZLl = ZLl(:,indexFlip);Zul = Zul(:,indexFlip);Zpl = Zpl(:,indexFlip);if isAnySlipFaceAlL = AlL(indexFlip,:);Alu = Alu(indexFlip,:);Alp = Alp(indexFlip,:);elseAlL = ALl(:,indexFlip)';Alu = Aul(:,indexFlip)';Alp = Apl(:,indexFlip)';endAll = All(indexFlip,indexFlip);fl = fl(indexFlip);Arl = Arl(:,indexFlip); Figure 24: Extract of the hdg Stokes ElementalMatrices function that defines the el-emental matrices of the global problem for the Stokes case, depending on the boundaryconditions imposed.
A major difference between the Poisson and Stokes case is the structure of the local prob-lems (9) and (27). Despite both matrices are symmetric, the one arising from the discreti-sation of the Poisson equation is positive definite, whereas it is indefinite in the Stokescase. More precisely, the matrix in (27) features a saddle-point structure, as classical inthe context of finite element approximations of incompressible flow problems [120]. Inaddition, since the HDG local problem is a purely Dirichlet boundary value problem, theconstraint (23b) is introduced using a Lagrange multiplier ζ .The structure of the symmetric indefinite matrix involved in the local problem is dis-played on the left-hand side of figure 25. On the right-hand side, the blocks of the firstand last columns are associated with the contribution of ˆu and ρ , respectively, whereasthe second column is related to the independent term of the equation. The figure reportsthe hybridisation stage in which the elemental matrices and vectors defined in (28) arecomputed 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.41 Elemental mapping ---------------------------------------------------------------------------Z = [ALL ALu zeros(ndofL,ndofP) zeros(ndofL,1);ALu' Auu Apu' zeros(ndofU,1);zeros(ndofP,ndofL) Apu zeros(ndofP) Arp;zeros(1,ndofL) zeros(1,ndofU) Arp' 0] \ ...[ALl fL zeros(ndofL,1);Aul fu zeros(ndofU,1);Apl fp zeros(ndofP,1);zeros(1,ndofUhat) 0 1]; Figure 25: Extract of the hdg Stokes ElementalMatrices function that computes thematrices and vectors in equation (27). % Elemental mapping ---------------------------------------------------------------------------Z = [ALL ALu zeros(ndofL,ndofP) zeros(ndofL,1);ALu' Auu Apu' zeros(ndofU,1);zeros(ndofP,ndofL) Apu zeros(ndofP) Arp;zeros(1,ndofL) zeros(1,ndofU) Arp' 0] \ ...[ALl fL zeros(ndofL,1);Aul fu zeros(ndofU,1);Apl fp zeros(ndofP,1);zeros(1,ndofUhat) 0 1];% Elemental matricesif isPureDirichletKe = [All + (AlL*ZLl + Alu*Zul + Alp*Zpl), AlL*zLr + Alu*zur + Alp*zpr;Arl, 0;ArlExtra'*Zpl ArlExtra'*zpr];fe = [fl - (AlL*zLf + Alu*zuf + Alp*zpf);fn;ArlExtra'*zpf];elseKe = [All + (AlL*ZLl + Alu*Zul + Alp*Zpl), AlL*zLr + Alu*zur + Alp*zpr;Arl, 0];fe = [fl - (AlL*zLf + Alu*zuf + Alp*zpf);fn];end Figure 26: Extract of the hdg Stokes ElementalMatrices function that computes theblock matrix and vector in equation (30). % Initialisationmat.i = zeros(1, hdg.rowsGlobalSystem);mat.j = zeros(1, hdg.columnsGlobalSystem);mat.Kij = zeros(1, hdg.nDOFglobal);f = zeros(hdg.rowsGlobalSystem, 1);% Element computation -------------------------[N, dNx, wXY, Xg] = ...gaussElemCartesianInfo(refElem(pElem), Xe);Nw = bsxfun(@times, N, wXY');NwN = Nw*N';for iNsd = 1:nsd2vElem = iNsd:nsd2:ndofL;ALL(vElem,vElem) = -NwN;endsourceTerm = ...stokes source(Xg,problemParams,matElem, ...nsd,problemParams.example);for iNsd = 1:nsdwElem = iNsd:nsd:ndofU;Apu(:,wElem) = dNx(:,:,iNsd)*Nw';fu(wElem) = Nw*sourceTerm(:,iNsd);for jNsd = 1:nsdzElem = (iNsd-1)*nsd+jNsd:nsd2:ndofL;ALu(zElem,wElem) = ...sqrt(nu)*dNx(:,:,jNsd)*Nw';endend% Additional contributions if any of the ...faces features a slip BCif isAnySlipFaceAlL = zeros(ndofUhat,ndofL);Alu = zeros(ndofUhat,ndofU);Alp = zeros(ndofUhat,ndofP);end % Flipping due to the different numbering ...of internal face nodes when seen from ...left or right elementZLl = ZLl(:,indexFlip);Zul = Zul(:,indexFlip);Zpl = Zpl(:,indexFlip);if isAnySlipFaceAlL = AlL(indexFlip,:);Alu = Alu(indexFlip,:);Alp = Alp(indexFlip,:);elseAlL = ALl(:,indexFlip)';Alu = Aul(:,indexFlip)';Alp = Apl(:,indexFlip)';endAll = All(indexFlip,indexFlip);fl = fl(indexFlip);Arl = Arl(:,indexFlip);% Additional restriction for purely ...Dirichlet BCArlExtra = zeros(ndofP,1);analyticalPressure = ...stokes analyticalPressure(Xg,...problemParams,matElem,nsd,...problemParams.example);% If purely Dirichlet BCif isPureDirichletArlExtra = sum(Nw,2);intPressure = intPressure + ...sum(Nw*analyticalPressure,1);end Figure 27: Extract of the hdg Stokes ElementalMatrices function that computes thevector a ρρ for the pressure constraint.Finally, hdg Stokes ElementalMatrices computes the elemental contribution to thematrix in the global problem (30) as reported in figure 26. It is worth noting that thevariable hdg.pureDirichlet is utilised here to discriminate the construction of the globalmatrix of a purely Dirichlet boundary value problem. More precisely, besides the blocks (cid:98) K , G and G T , also the vector a ρρ (i.e. ArlExtra ) arising from the imposition of the globalconstraint (20) is taken into account. An extract of the hdg Stokes ElementalMatrices function computing such a vector is displayed in figure 27.42
Initialisationmat.i = zeros(1, hdg.rowsGlobalSystem);mat.j = zeros(1, hdg.columnsGlobalSystem);mat.Kij = zeros(1, hdg.nDOFglobal);f = zeros(hdg.rowsGlobalSystem, 1); Figure 28: Extract of the hdg Stokes Globalystem function that initialises the structuresrequired for the assembly of the global system.
Remark 9.
The assembly of the block matrix reported in figure 26 does not exploit thesymmetry 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 localproblem into the global problem (29), leading to G := A n el e =1 (cid:2) A ˆ uL A ˆ uu A ˆ up (cid:3) e z ρL z ρu z ρp z ρζ e . (37)The interested reader can thus employ the provided code to numerically verify that thematrix G defined in equation (37) is the transpose of the one introduced in (31). Theformal proof can be devised following the rationale described in [151]. As described in section 4, the global problem for the Stokes equations features a saddle-point structure. Hence, the assembly of the global system presents major differences withrespect to the Poisson case previously discussed. First, figure 28 reports the initialisationof 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 Dirich-let 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 isdisplayed in figure 29. According to the variable hdg.pureDirichlet , the above mentionedconstraint 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 theconstraint is imposed via a Lagrange multiplier (Fig. 30).
The structure of the code for the local problem in the Stokes case replicates the onepresented in figure 20 for the Poisson equation and is thus omitted. The only difference lies43
Assembly[indexGlobalFace, ⇠ ] = ...hdgElemToFaceIndex(mesh.indexTf, ...refFace, hdg.faceInfo(iElem), ...refElem(pElem).nOfFaces, ...ctt.nOfComponents);if hdg.pureDirichletiBlock = [indexGlobalFace, ...hdg.nDOFglobal+iElem, ...hdg.nDOFglobal+mesh.nOfElements+1];jBlock = [indexGlobalFace, ...hdg.nDOFglobal+iElem];elseiBlock = [indexGlobalFace, ...hdg.nDOFglobal+iElem];jBlock = iBlock;endnI =numel(iBlock);nJ =numel(jBlock);currentIndex = indexIni + (1:nJ);currentIndex2 = indexIni2 + (1:nI*nJ);for i=1:nImat.i(currentIndex) = iBlock(i);mat.j(currentIndex) = jBlock;currentIndex = currentIndex + nJ;endmat.Kij(currentIndex2) = ...reshape(Ke', 1, nI*nJ);if hdg.pureDirichletf(hdg.rowsGlobalSystem) = ...-f(hdg.rowsGlobalSystem)+intPressure;K = [K, [K(end,:)'; 0]];end Figure 29: Extract of the hdg Stokes Globalystem function that constructs the structuresrequired for the assembly of the global system. % Initialisationmat.i = zeros(1, hdg.rowsGlobalSystem);mat.j = zeros(1, hdg.columnsGlobalSystem);mat.Kij = zeros(1, hdg.nDOFglobal);f = zeros(hdg.rowsGlobalSystem, 1);% Element computation -------------------------[N, dNx, wXY, Xg] = ...gaussElemCartesianInfo(refElem(pElem), Xe);Nw = bsxfun(@times, N, wXY');NwN = Nw*N';for iNsd = 1:nsd2vElem = iNsd:nsd2:ndofL;ALL(vElem,vElem) = -NwN;endsourceTerm = ...stokes source(Xg,problemParams,matElem, ...nsd,problemParams.example);for iNsd = 1:nsdwElem = iNsd:nsd:ndofU;Apu(:,wElem) = dNx(:,:,iNsd)*Nw';fu(wElem) = Nw*sourceTerm(:,iNsd);for jNsd = 1:nsdzElem = (iNsd-1)*nsd+jNsd:nsd2:ndofL;ALu(zElem,wElem) = ...sqrt(nu)*dNx(:,:,jNsd)*Nw';endend% Additional contributions if any of the ...faces features a slip BCif isAnySlipFaceAlL = zeros(ndofUhat,ndofL);Alu = zeros(ndofUhat,ndofU);Alp = zeros(ndofUhat,ndofP);end % Flipping due to the different numbering ...of internal face nodes when seen from ...left or right elementZLl = ZLl(:,indexFlip);Zul = Zul(:,indexFlip);Zpl = Zpl(:,indexFlip);if isAnySlipFaceAlL = AlL(indexFlip,:);Alu = Alu(indexFlip,:);Alp = Alp(indexFlip,:);elseAlL = ALl(:,indexFlip)';Alu = Aul(:,indexFlip)';Alp = Apl(:,indexFlip)';endAll = All(indexFlip,indexFlip);fl = fl(indexFlip);Arl = Arl(:,indexFlip);% Additional restriction for purely ...Dirichlet BCArlExtra = zeros(ndofP,1);analyticalPressure = ...stokes analyticalPressure(Xg,...problemParams,matElem,nsd,...problemParams.example);% If purely Dirichlet BCif isPureDirichletArlExtra = sum(Nw,2);intPressure = intPressure + ...sum(Nw*analyticalPressure,1);endif hdg.pureDirichletf(hdg.rowsGlobalSystem) = ...-f(hdg.rowsGlobalSystem)+intPressure;K = [K, [K(end,:)'; 0]];end Figure 30: Extract of the hdg Stokes Globalystem function to impose the constraint inthe global system for purely Dirichlet boundary value problems. p(indexPIni:indexPEnd) = local(iElem).Zpl*uHat(indexGlobalFace) + local(iElem).zpr*rho(iElem) ...+ local(iElem).zpf;u(indexUIni:indexUEnd) = local(iElem).Zul*uHat(indexGlobalFace) + local(iElem).zur*rho(iElem) ...+ local(iElem).zuf;L(indexLIni:indexLEnd) = local(iElem).ZLl*uHat(indexGlobalFace) + local(iElem).zLr*rho(iElem) ...+ local(iElem).zLf; Figure 31: Extract of the hdg Stokes LocalProblem function where the elemental valuesof the primal and mixed unknowns are computed.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 isdisplayed in figure 31, focusing on the operations in equation (28).44 The use of a high-order functional approximation means that non-standard techniques arerequired to produce a reliable representation of the solution in each element. With the in-creased popularity of high-order methods in recent years, different strategies to efficientlydisplay high-order solutions have been proposed [188, 198, 220]. The
HDGlab provides thecapability to accurately display high-order solutions, including curved isoparametric ele-ments.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 thefollowing information: • resolution : Takes value 1 for a faster but less accurate representation of high-ordersolutions and value 2 for a slower but more accurate representation of high-ordersolutions. • 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 solutionto be represented. • alphaFace : Sets the transparency between 0 and 1.In addition, for the Stokes equation, two problem-specific options are available in the datastructure plotOpts : • componentsU : A boolean variable allowing the predefined visualisation of all thecomponents of the velocity vector. This functionality relies on the definition of thevector 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 withequally-spaced and Fekete nodal sets in the directory dat/postprocess . In two dimen-sions, the data structure postproc contains the following information: • nOfNodesPlot : Number of nodes used to display the high-order solution in eachelement. • nodesPlot : Array of dimension nOfNodesPlot ×
2. Coordinates of the nodes, in thereference element, used to display the high-order solution.45igure 32: Submesh of the reference element used to provide an accurate representation ofhigh-order solutions in each element. The left picture corresponds to resolution =1 andthe right picture to resolution =2. • nSubElemsOnePlot : Number of subelements used to display the high-order solutionin each element. • connecNodesPlot : Array of dimension nSubElemsOnePlot × • nOfEdges : Number of edges of the element. • edgeNodesSplit : Array of dimension × nOfEdges , where r is the resolution se-lected by the user in the data structure plotOpts . The i -th column contains thelocal 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 approximationused in all the elements. The component p of elem contains a field called N , ofdimension nOfNodesPlot × p ( p + 1) /
2, that stores the value of the shape functions oforder p at the positions given by nodesPlot . This information is used to interpolatethe solution at the nodes of the submesh, providing a more accurate representationof 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 functionsof order p at the positions of an edge. This information is used to interpolate thesolution at the edges of the submesh.The submeshes used for a triangular element with resolution =1 and resolution =2are 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 ofthe 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: 46igure 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 ineach element face. – connecPlot : Array of dimension nElemPlot × – nNodesPlot : Number of nodes used to display the high-order solution on eachelement face. It is given by (5 r + 1)(5 r + 2) /
2, where r is the resolution selectedby the user in the data structure plotOpts . – edgeNodesPlot : Array of dimension 1 × (15 r + 1), where r is the resolutionselected by the user in the data structure plotOpts . It contains the local numberof the face nodes that belong to the edges of the face. • Elem : Array of dimension 1 × p max , where p max is the maximum degree of approximationused in all the elements. The component p of Elem contains: – face : Array of dimension 1 × i -th component contains the localnumber of the vertices on the i -th face and the value of the element shapefunctions of order p at the nodal distribution used for plotting the solution onthe i -th face. – coord : Array of dimension p ( p +1)( p +2) / p . • faceVertices : Array of dimension 4 ×
3. The i -th row contains the local number ofthe vertices on the i -th face.The third data structure utilised during the postprocess stage is called visual andit is built by the function buildSubmeshPostprocess2D in two dimensions and by theroutine buildSubmeshPostprocess3D in three dimensions. This data structure containsthe following information: • X : Physical coordinates of all the nodes of the submesh used to display the high-ordersolution. 47 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 usedif the user sets fieldsWithNodes =1 in the data structure plotOpts . • edges : Structure containing the list of nodes of the submesh that form the high-orderrepresentation of the physical edges of the mesh.In a separate function, the data structure is updated by adding the field U that containsthe interpolated values of the solution on the submesh used to display the high-order solu-tion. 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 indepen-dent on the field to be represented and, therefore, it is only called once. Instead, the secondfunction that updates the data structure visual with the field U depends upon the fieldto 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 functionthat 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 representingthe solution only in a region of the computational domain. The user can set the valueof a string, called conditionPlot , that specifies a region in the physical space. Before visual is computed, the function selectFacesToPlot3D computes the list of faces inthe computational mesh that satisfy the condition given by conditionPlot . Then, thesubmesh and interpolation of the solution is only performed over the faces that satisfy thecondition specified by the user.Finally, the visualisation function also accounts for the need to represent the super-convergent 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 thedegree 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 solutioncan be used to postprocess the higher order superconvergent solution.
11 Numerical examples
In this section, several numerical examples showing the capabilities of the
HDGlab solversfor the Poisson and Stokes equations are presented. As mentioned in sections 3 and 4, thechoice of the stabilisation parameter τ is critical for the accuracy of the HDG approxima-tion. For the examples involving the Poisson equation, the definition τ = c P κ/(cid:96) is considered,48igure 34: Two-dimensional Wang flow. Convergence of the L (Ω) error of pressure, p ,mixed variable, L (left), primal, u , and postprocessed, u (cid:63) , velocities (right) as a functionof the characteristic mesh size h for polynomial degree of approxiomation p =1 , . . . , (cid:96) is a characteristic length of the domain and c P a scaling factor selected equal to1 [175]. Following [151], the stabilisation for the Stokes cases is defined as τ = c S ν/(cid:96) , (cid:96) , thescaling factor being c S =3. The optimal convergence properties of the proposed HDG implementation are presentedfor 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 , isconsidered. The analytical velocity field is u ( x ) = (cid:40) ax − bλ cos( λx ) exp {− λx } bλ sin( λx ) exp {− λx } (cid:41) , (38)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. Thesource term s and the boundary conditions are computed starting from the analyticalsolution above. More precisely, a pseudo-traction g is applied on the bottom surfaceΓ N := { ( x , x ) ∈ Ω | x =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 (Ω)norm, of the primal, mixed and postprocessed variables as a function of the characteristicmesh 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 thepostprocessed velocity u (cid:63) . 49igure 35: Three-dimensional manufactured Stokes flow. Convergence of the L (Ω) errorof pressure, p , mixed variable, L (left), primal, u , and postprocessed, u (cid:63) , velocities (right)as a function of the characteristic mesh size h for polynomial degree of approxiomation p =1 , . . . , , ,with the following manufactured solution u ( x ) = n + ( x − x ) sin( x − n ) m − x (cid:20)(cid:0) x − x (cid:1) cos( x − n ) + (cid:0) x − x (cid:1) cos( x − n ) (cid:21) n + ( x − x ) sin( x − n ) , (39) p ( x ) = x (1 − x ) + x (1 − x ) + x (1 − x ) , (40)where m =1 and n = . The kinematic viscosity is set to ν =1, a Neumann datum is imposedon the boundary Γ N := { ( x , x , x ) ∈ Ω | x =0 } , whereas Dirichlet conditions are prescribedon Γ D = ∂ Ω \ Γ N .The optimal convergence of order p +1 of the relative L (Ω) error for velocity, pres-sure and gradient of velocity and the superconvergence of the postprocessed velocity areconfirmed in the 3D case by the results in figure 35. The coaxial Couette flow [64] is considered to show the optimal convergence properties ofthe
HDGlab solver using a high-order isoparametric approximation in a domain featuringcurved boundaries.This test consists of an incompressible viscous flow within two coaxial circular cylindersof infinite length and radius R int =1 and R ext =5, respectively. The computational domain isdefined as a section of the 3D cylinders, that is, Ω= { ( x , x ) ∈ R | R int ≤ r ≤ R ext } , where50igure 36: First level of refinement of the third-order mesh used for the convergence studyof the two-dimensional Couette flow (left) and module of the computed velocity (right). r := (cid:112) x + x is the distance to the axis of the cylinders. Dirichlet boundary conditionsenforcing the value of the angular velocities ω int =0 and ω ext =1 /R ext are imposed on theinternal and external boundary, respectively. The analytical expression of the azimuthalcomponent of the velocity is u φ = R ω ext − R ω int R − R r + ( ω int − ω ext ) R R R − R r . (41)Of course, being a purely Dirichlet boundary value problem, the constraint on the meanvalue 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 con-structed using the strategy described in [212]. Figure 36 displays the first level of meshrefinement featuring 128 triangular elements of polynomial degree 3 and the module of thecomputed velocity.The convergence of the relative error, measured in the L (Ω) norm, of the primal, mixedand postprocessed variables is reported in figure 37 as a function of the characteristic meshsize. Optimal convergence of the primal and mixed variables and superconvergence of thepostprocessed variable is achieved also in presence of high-order curved meshes. In this section, the flexibility of
HDGlab to devise a non-uniform polynomial degree ap-proximation in the domain is presented. This case, inspired by the study on micromixersin [131], consists of the flow in a microchannel with five obstacles. The problem setup fea-tures a parabolic inlet velocity profile and homogeneous Dirichlet and Neumann conditionson the top/bottom walls and on the outlet, respectively.51igure 37: Two-dimensional Couette flow. Convergence of the L (Ω) error of pressure, p ,mixed variable, L (left), primal, u , and postprocessed, u (cid:63) , velocities (right) as a functionof the characteristic mesh size h for polynomial degree of approxiomation p =1 , . . . , , . × [ − . , .
5] and the obstacles, attached to the topand bottom walls have thickness 0 . .
5. A mesh with local element size rangingbetween 0 .
08 and 0 .
19 is generated without any specific a priori refinement. It is worthnoting 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-uniformpolynomial degree distribution is generated following the adaptivity strategy describedin [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 wherelocalised flow features appear, whereas low-order approximations are utilised in regionfurther away.The module of the velocity on the described mesh is also reported in figure 38. It isworth noting that the mesh structure featuring the non-uniform polynomial approximationis provided as a datum for this test case. The corresponding simulation is performedseamlessly in
HDGlab and no specific intervention is required to the user.
Finally, the Stokes flow past a sphere is considered. The domain Ω=[ − H, L ] × [ − H, H ] × [ − H, H ] \ B , is defined, with H =5, L =10 and B , being the sphere of radius 1 centredin the point (0 , , a) Degree distribution(b) Module of the velocity Figure 38: (a) Mesh and degree of approximation and (b) module of the velocity for theflow in a microchannel with obstacles using non-uniform polynomial degree p between 1and 7.A high-order mesh featuring 1,036 tetrahedral elements is generated via the solid me-chanics analogy described in [212, 269]. Figure 39 displays the module of the velocity andthe pressure field computed using an isoparametric approximation of degree 6.It is worth noting that the computations in sections 11.1, 11.2, 11.3 and 11.4 are per-formed using a unique HDGlab solver for the Stokes equations, independently on the numberof spatial dimensions of the problem under analysis. Indeed,
HDGlab provides a seamlessimplementation of the HDG method, in which all relevant information is extracted fromthe mesh structure and the user is required to specify only the physical parameters andthe boundary conditions to setup a test case.
The next example, taken from [187], shows the solution of an electrostatic problem gov-erned by the Poisson equation in three dimensions. The domain of interest corresponds tothe exterior of 11 conducting spheres and it is discretised with a mesh of 35,895 quadratictetrahedral elements, as shown in figure 40. This figure shows one of the implementedcapabilities of the postprocessing library to display the faces corresponding to the exteriorand interior faces of the mesh separately, with different colour and transparencies in eachcase. The first plot in figure 40 is produced by selecting the faces corresponding to the farfield boundary and using a transparency. In a second phase, the exterior faces correspond-ing to the conducting spheres are displayed with no transparency. It is worth noting thatthe far field boundary and the boundary corresponding to the conducting spheres couldbe distinguished using either the boundary condition flag or simply imposing a conditionon the faces to be displayed. The second plot of figure 40 is also obtained in two stages.53 a) Module of the velocity (b) Pressure
Figure 39: (a) Module of the velocity and (b) pressure field of the external flow past aquarter of a sphere using polynomial degree p =6.Figure 40: Two views of the quadratic tetrahedral mesh used to solve the electrostaticproblem. 54 a) Conducting spheres (b) Electric field Figure 41: (a) The 11 conducting spheres coloured according to the boundary conditionimposed. Red colour is used for the spheres where a positive potential is imposed, whereasblue is used for the spheres where a negative potential is imposed. (b) Magnitude of theelectric field on the surface of the 11 conducting spheres.First, the interior faces satisfying a condition corresponding to a positive x coordinate aredisplayed with a transparency. Second, the exterior faces corresponding to the conductingspheres is displayed with no transparency. When representing the interior and exteriorfaces, 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 computa-tional domain. A positive electrostatic potential of magnitude 5 is imposed on the centralsphere and five of the surrounding spheres, whereas a negative potential of magnitude − HDGlab , notonly the primal variable corresponding to the electrostatic potential is obtained but alsoits gradient, which corresponds to the electric field in this example. Figure 41(b) showsthe 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 di-mensional mechanical component. The domain is discretised with a mesh of 6,465 elementswith p = 4, as illustrated in figure 42. The first plot in figure 42 includes the high-ordernodal distribution and the second plot in figure 42 shows the possibilities offered by thepostprocessing library provided within HDGlab . In fact, it shows an extra functionality thatcan be added to display the intersection curves of the underlying CAD geometry when theuser have access to this data.Dirichlet and Neumann boundary conditions are imposed on the blue and red portionsof the boundary, respectively, as shown in figure 43(a). In the Dirichlet part of the bound-ary a temperature equal to 10 is imposed, whereas the Neumann boundary condition ishomogeneous, imposing that part of the boundary is perfectly insulated. The temperature55igure 42: Two views of the fourth order tetrahedral mesh used to solve the heat transferproblem. (a) Mechanical component (b) Temperature distribution
Figure 43: (a) The mechanical component coloured according to the boundary conditionimposed. Red colour denotes a homogeneous Neumann boundary condition and the bluecolour a Dirichlet boundary condition. (b) Temperature distribution on the surface of themechanical component. 56igure 44: Surface mesh of a generic UAV.Figure 45: Magnitude of the velocity field and pressure field on the surface of a genericUAV.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 un-manned 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 unitvelocity, is imposed on the inflow part of the boundary and a Dirichlet boundary conditioncorresponding to zero potential on the outflow. On the rest of the boundary a homoge-neous Neumann boundary condition is enforced to represent a physical wall. After solvingthe problem with HDGlab , the flow potential and the velocity field are obtained. Figure 45shows 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 involv-ing complex geometries, where the resulting global system has more than one million ofequations. Despite the code was not developed with computational efficiency in mind, itstill enables the interested users to solve relatively large problems and, with some addi-tional improvements in terms of performance, it can be used for larger three dimensionalproblems.
12 Concluding remarks
An open-source Matlab implementation of the HDG method for elliptic problems has beenpresented, the so-called
HDGlab . The code is capable of solving problems governed by57he Poisson and Stokes equations using high-order simplicial elements, including curvedelements, by means of an isoparametric formulation.
HDGlab provides a suitable environment to those interested in the programming ofhybrid methods in general and the HDG method in particular. The paper describes theHDG formulation employed for the Poisson and Stokes solvers and provides a very detaileddescription of the code, with particular emphasis in the data structures utilised. Thepresentation of the code involves the preprocess, computation and postprocess stages, andincludes detailed explanations of the most relevant functions. A set of examples is presentedto illustrate the use and the potential of
HDGlab . The examples go from simple test caseswith a known analytical solution, used to demonstrate the numerical properties of theHDG method, to more complex such as the computation of the potential flow around aUAV using high-order curved meshes.
HDGlab is available as an open-source software, released under the terms of the GNUGeneral Public License version 3.0 or any later version ( )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 formof the HDG local and global problems for the Poisson equation are presented.An isoparametric formulation is considered and the coordinates ξ = { ξ , . . . , ξ n sd } T ofa reference element (cid:101) Ω( ξ ) are mapped to the coordinates x = { x , . . . , x n sd } T of the localelement Ω( x ) by the transformation x ( ξ ) = n en (cid:88) i =1 N i ( ξ ) x i , where { x i } i =1 ,..., n en denotes the vector of the nodal coordinates of the element. Hence, theisoparametric approximations introduced in (8) are defined in a reference element (cid:101) Ω( ξ ) for u and q and on a reference face (cid:101) Γ( η ), η = { η , . . . , η n sd − } T for ˆ u , the corresponding shapefunctions being N ( ξ ) and ˆ N ( η ), respectively.For the sake of readability, introduce the compact forms of the shape functions and58heir derivatives, namely N := (cid:2) N N . . . N n en (cid:3) T , (42a) (cid:99) N := (cid:2) ˆ N ˆ N . . . ˆ N n fn (cid:3) T , (42b) N n := (cid:2) N n N n . . . N n en n (cid:3) T , (42c) N n sd := (cid:2) N I n sd N I n sd . . . N n en I n sd (cid:3) T , (42d) Q := (cid:2) ( J − ∇ N ) T ( J − ∇ N ) T . . . ( J − ∇ N n en ) T (cid:3) T , (42e)where n denotes the outward unit normal vector to a face and J is the Jacobian of theisoparametric mapping.The solution of the HDG local problem (9) involves the matrices and vectors[ A uu ] e := n efa (cid:88) f =1 τ f n fip (cid:88) g =1 N ( ξ fg ) N T ( ξ fg ) | J ( ξ fg ) | w fg , [ A uq ] e := √ κ n eip (cid:88) g =1 N ( ξ eg ) Q T ( ξ eg ) | J ( ξ eg ) | w eg , [ A qq ] e := − n eip (cid:88) g =1 N n sd ( ξ eg ) N T n sd ( ξ eg ) | J ( ξ eg ) | w eg , [ A u ˆ u ] e := n efa (cid:88) f =1 τ f (cid:18) n fip (cid:88) g =1 N ( ξ fg ) (cid:99) N T ( ξ fg ) | J ( ξ fg ) | w fg (cid:19)(cid:0) − χ f D (cid:1) , [ A q ˆ u ] e := √ κ n efa (cid:88) f =1 (cid:18) n fip (cid:88) g =1 N n ( ξ fg ) (cid:99) N T ( ξ fg ) | J ( ξ fg ) | w fg (cid:19)(cid:0) − χ f D (cid:1) , [ f u ] e := n eip (cid:88) g =1 N ( ξ eg ) s (cid:0) x ( ξ eg ) (cid:1) | J ( ξ eg ) | w eg + n efa (cid:88) f =1 τ f (cid:18) n fip (cid:88) g =1 N ( ξ fg ) u D (cid:0) x ( ξ fg ) (cid:1) | J ( ξ fg ) | w fg (cid:19) χ f D , [ f q ] e := √ κ n efa (cid:88) f =1 (cid:18) n fip (cid:88) g =1 N n ( ξ fg ) u D (cid:0) x ( ξ fg ) (cid:1) | J ( ξ fg ) | w fg (cid:19) χ f D , where n efa denotes the number of faces of the element Ω e , ξ eg are the n eip integration pointsdefined on the reference element and ξ fg are the n fip integration points of the referenceface. The corresponding weights for the integration points are denoted by w eg and w fg ,respectively. In addition, the indicator function χ f (cid:3) is defined as χ f (cid:3) := (cid:40) f belongs to Γ (cid:3) , . (43)59imilarly, for the HDG global problem (11) the following matrix and vector[ A ˆ u ˆ u ] e := − n efa (cid:88) f =1 τ f (cid:18) n fip (cid:88) g =1 (cid:99) N ( ξ fg ) (cid:99) N T ( ξ fg ) | J ( ξ fg ) | w fg (cid:19)(cid:0) − χ f D (cid:1) , [ f ˆ u ] e := − n efa (cid:88) f =1 (cid:18) n fip (cid:88) g =1 (cid:99) N ( ξ fg ) g (cid:0) x ( ξ fg ) (cid:1) | J ( ξ fg ) | w fg (cid:19) χ f N , are defined.To describe the terms appearing in the HDG postprocess (17), the compact forms N (cid:63) := (cid:2) N (cid:63) N (cid:63) . . . N (cid:63) n (cid:63) en (cid:3) T , (44a) N (cid:63) n sd := (cid:2) N (cid:63) I n sd N (cid:63) I n sd . . . N (cid:63) n (cid:63) en I n sd (cid:3) T , (44b) Q (cid:63) := (cid:2) ( J − (cid:63) ∇ N (cid:63) ) T ( J − (cid:63) ∇ N (cid:63) ) T . . . ( J − (cid:63) ∇ N (cid:63) n (cid:63) en ) T (cid:3) T , (44c)are introduced, where { N (cid:63)i } i =1 ,..., n (cid:63) en denote the shape functions for the discretisation of u (cid:63) , J (cid:63) is the Jacobian of the corresponding isoparametric transformation and n (cid:63) en is the numberof elemental nodes of the approximation.The following matrices and vector are thus defined as[ A (cid:63)(cid:63) ] e := κ n (cid:63) ip (cid:88) g =1 Q (cid:63) ( ξ (cid:63) g ) Q (cid:63) T ( ξ (cid:63) g ) | J (cid:63) ( ξ (cid:63) g ) | w (cid:63) g , [ a (cid:63)λ ] e := n (cid:63) ip (cid:88) g =1 N (cid:63) ( ξ (cid:63) g ) | J (cid:63) ( ξ (cid:63) g ) | w (cid:63) g , [ A (cid:63)q ] e := − √ κ n (cid:63) ip (cid:88) g =1 Q (cid:63) ( ξ (cid:63) g ) N (cid:63) T n sd ( ξ (cid:63) g ) | J (cid:63) ( ξ (cid:63) g ) | w (cid:63) g , where ξ (cid:63) g and w (cid:63) g are the n (cid:63) ip integration points and the weights associated with the higherorder approximation in the space V h(cid:63) . B Matrices and vectors for the Stokes solver
In this appendix, the expressions of the matrices and vectors appearing in the discrete formof the HDG local and global problems for the Stokes equations are presented. It is worthrecalling that in this problem, the primal, u , and hybrid, ˆu , variables are n sd -dimensionalvector-valued unknowns and the mixed variable L is a tensor-valued unknown of dimension n sd × n sd . 60n addition to the compact forms of the shape functions presented in equations (42)and (44), the following definitions are introduced for the vectorial problem N m sd := (cid:2) N I m sd N I m sd . . . N n en I m sd (cid:3) T , (45a) N (cid:63) m sd := (cid:2) N (cid:63) I m sd N (cid:63) I m sd . . . N (cid:63) n (cid:63) en I m sd (cid:3) T , (45b) (cid:99) N n sd := (cid:2) ˆ N I n sd ˆ N I n sd . . . ˆ N n fn I n sd (cid:3) T , (45c) (cid:99) N E := (cid:2) ˆ N E ˆ N E . . . ˆ N n fn E (cid:3) T , (45d) (cid:99) N D := (cid:2) ˆ N D ˆ N D . . . ˆ N n fn D (cid:3) T , (45e)where m sd := n sd . Moreover, for each component j =1 , . . . , n sd , the compact forms N n j := (cid:2) N n j N n j . . . N n en n j (cid:3) T , (45f) Q j := (cid:2) [ J − ∇ N ] j [ J − ∇ N ] j . . . [ J − ∇ N n en ] j (cid:3) T , (45g) Q (cid:63) j := (cid:2) [ J − (cid:63) ∇ N (cid:63) ] j [ J − (cid:63) ∇ N (cid:63) ] j . . . [ J − (cid:63) ∇ N (cid:63) n (cid:63) en ] j (cid:3) T , (45h)are defined, n j being the j -th component of the outward unit normal vector n to the face.The solution of the HDG local problem (27) involves the matrices and vectors[ A LL ] e := − n eip (cid:88) g =1 N m sd ( ξ eg ) N T m sd ( ξ eg ) | J ( ξ eg ) | w eg , [ A Lu ] e := √ ν n sd (cid:88) j =1 n eip (cid:88) g =1 Q j ( ξ eg ) N T n sd ( ξ eg ) | J ( ξ eg ) | w eg , [ A uu ] e := n efa (cid:88) f =1 τ f n fip (cid:88) g =1 N n sd ( ξ fg ) N T n sd ( ξ fg ) | J ( ξ fg ) | w fg , [ A pu ] e := n eip (cid:88) g =1 Q ( ξ eg ) N T n sd ( ξ eg ) | J ( ξ eg ) | w eg , [ a ρp ] e := n efa (cid:88) f =1 n fip (cid:88) g =1 N ( ξ fg ) | J ( ξ fg ) | w fg , [ A L ˆ u ] e := √ ν n efa (cid:88) f =1 (cid:18) n sd (cid:88) j =1 n fip (cid:88) g =1 N n j ( ξ fg ) (cid:99) N T n sd ( ξ fg ) | J ( ξ fg ) | w fg (cid:19)(cid:0) − χ f D (cid:1) , [ A u ˆ u ] e := n efa (cid:88) f =1 τ f (cid:18) n fip (cid:88) g =1 N n sd ( ξ fg ) (cid:99) N T n sd ( ξ fg ) | J ( ξ fg ) | w fg (cid:19)(cid:0) − χ f D (cid:1) , [ A p ˆ u ] e := n efa (cid:88) f =1 (cid:18) n fip (cid:88) g =1 N n ( ξ fg ) (cid:99) N T n sd ( ξ fg ) | J ( ξ fg ) | w fg (cid:19)(cid:0) − χ f D (cid:1) , f L ] e := √ ν n efa (cid:88) f =1 (cid:18) n sd (cid:88) j =1 n fip (cid:88) g =1 N n j ( ξ fg ) u D (cid:0) x ( ξ fg ) (cid:1) | J ( ξ fg ) | w fg (cid:19) χ f D , [ f u ] e := n eip (cid:88) g =1 N n sd ( ξ eg ) s (cid:0) x ( ξ eg ) (cid:1) | J ( ξ eg ) | w eg + n efa (cid:88) f =1 τ f (cid:18) n fip (cid:88) g =1 N n sd ( ξ fg ) u D (cid:0) x ( ξ fg ) (cid:1) | J ( ξ fg ) | w fg (cid:19) χ f D , [ f p ] e := n efa (cid:88) f =1 (cid:18) n fip (cid:88) g =1 N n ( ξ fg ) u D (cid:0) x ( ξ fg ) (cid:1) | J ( ξ fg ) | w fg (cid:19) χ f D . Similarly, the matrices and vectors for the global problem (29) are[ A ˆ uL ] e := √ ν n efa (cid:88) f =1 (cid:40) − (cid:18) n sd (cid:88) j =1 n fip (cid:88) g =1 (cid:99) N E ( ξ fg ) N Tn j ( ξ fg ) | J ( ξ fg ) | w fg (cid:19) χ f S + (cid:18) n sd (cid:88) j =1 n fip (cid:88) g =1 (cid:99) N n sd ( ξ fg ) N Tn j ( ξ fg ) | J ( ξ fg ) | w fg (cid:19)(cid:0) − χ f D (cid:1)(cid:0) − χ f S (cid:1)(cid:41) , [ A ˆ uu ] e := n efa (cid:88) f =1 τ f (cid:40) − (cid:18) n fip (cid:88) g =1 (cid:99) N E ( ξ fg ) N T n sd ( ξ fg ) | J ( ξ fg ) | w fg (cid:19) χ f S + (cid:18) n fip (cid:88) g =1 (cid:99) N n sd ( ξ fg ) N T n sd ( ξ fg ) | J ( ξ fg ) | w fg (cid:19)(cid:0) − χ f D (cid:1)(cid:0) − χ f S (cid:1)(cid:41) , [ A ˆ up ] e := n efa (cid:88) f =1 (cid:40) − (cid:18) n fip (cid:88) g =1 (cid:99) N E ( ξ fg ) N Tn ( ξ fg ) | J ( ξ fg ) | w fg (cid:19) χ f S + (cid:18) n fip (cid:88) g =1 (cid:99) N n sd ( ξ fg ) N Tn ( ξ fg ) | J ( ξ fg ) | w fg (cid:19)(cid:0) − χ f D (cid:1)(cid:0) − χ f S (cid:1)(cid:41) , [ A ˆ u ˆ u ] e := n efa (cid:88) f =1 (cid:40)(cid:18) n fip (cid:88) g =1 (cid:99) N D ( ξ fg ) (cid:99) N T n sd ( ξ fg ) | J ( ξ fg ) | w fg (cid:19) χ f S + τ f (cid:18) n fip (cid:88) g =1 (cid:99) N E ( ξ fg ) (cid:99) N T n sd ( ξ fg ) | J ( ξ fg ) | w fg (cid:19) χ f S − τ f (cid:18) n fip (cid:88) g =1 (cid:99) N n sd ( ξ fg ) (cid:99) N T n sd ( ξ fg ) | J ( ξ fg ) | w fg (cid:19)(cid:0) − χ f D (cid:1)(cid:0) − χ f S (cid:1)(cid:41) , [ f ˆ u ] e := − n efa (cid:88) f =1 (cid:18) n fip (cid:88) g =1 (cid:99) N n sd ( ξ fg ) g (cid:0) x ( ξ fg ) (cid:1) | J ( ξ fg ) | w fg (cid:19) χ f N . In addition, if a problem with purely Dirichlet boundary conditions is solved, the con-62traint (20) gives rise to the vector[ a ρρ ] e := n eip (cid:88) g =1 N ( ξ eg ) | J ( ξ eg ) | w eg , (46)which is used to enforce the constraint using an appropriately defined Lagrange multiplier.Finally, for the postprocess (32), the following matrices are required[ A (cid:63)(cid:63) ] e := ν n sd (cid:88) j =1 n (cid:63) ip (cid:88) g =1 Q (cid:63) j ( ξ (cid:63) g ) Q (cid:63) Tj ( ξ (cid:63) g ) | J (cid:63) ( ξ (cid:63) g ) | w (cid:63) g , [ A (cid:63)λ ] e := n (cid:63) ip (cid:88) g =1 N (cid:63) n sd ( ξ (cid:63) g ) | J (cid:63) ( ξ (cid:63) g ) | w (cid:63) g , [ A (cid:63)L ] e := − √ ν n sd (cid:88) j =1 n (cid:63) ip (cid:88) g =1 Q (cid:63) j ( ξ (cid:63) g ) N (cid:63) T m sd ( ξ (cid:63) g ) | J (cid:63) ( ξ (cid:63) g ) | w (cid:63) g . 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 insection 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 equally-spaced or Fekete points for the approximation.
Gmsh provides mesh discretisation basedon 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 concerningthe 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. • isPlotMesh : Boolean variable to activate the option for visualising the importedmesh with the nodal distribution. 63 (indexPIni:indexPEnd) = local(iElem).Zpl*uHat(indexGlobalFace) + local(iElem).zpr*rho(iElem) ...+ local(iElem).zpf;u(indexUIni:indexUEnd) = local(iElem).Zul*uHat(indexGlobalFace) + local(iElem).zur*rho(iElem) ...+ local(iElem).zuf;L(indexLIni:indexLEnd) = local(iElem).ZLl*uHat(indexGlobalFace) + local(iElem).zLr*rho(iElem) ...+ local(iElem).zLf;optionNodes = 0; % GMSH utilises uniform nodal distributions%% Option setupfileName = 'ringH1unifP2';nsd = 2; % Number of spatial dimensionspDegree = 2; % Polynomial degreeisPlotMesh = 1; % Boolean variable to plot the imported meshoutputPath = 'meshFiles';%% Construct the reference element[refElem, refFace] = getRefData(pDegree, nsd, optionNodes);refElem = refElem(1, pDegree);refFace = refFace(pDegree, pDegree);%% Store the coordinates of the nodes with duplication[X, T, faces, matElem] = scanMeshFileMSH(fileName, pDegree, nsd);%% Extract internal and external facesintFaces = getInternalFaces(T, refElem, refFace);extFaces = getBoundaryFaces(T, faces, refElem, refFace);%% Build the mat structure of the meshmesh = buildMeshStruct(X, T, matElem, intFaces, extFaces, refElem, optionNodes);%% Save mesh structure in .mat filemeshFile = sprintf('%s/%s', outputPath, fileName);save(meshFile, 'mesh');%% Visualise the imported meshif isPlotMeshvisualiseMeshend Figure 46: Function convertMSHtoMAT to convert the mesh file from .msh to .mat format. • outputPath : String that specifies the location where the imported mesh will bestored. 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 isextracted: • X : Array containing the coordinates of the mesh nodes. • T : Connectivity matrix featuring the identifiers of the nodes associated with eachelement. • faces : Array containing the information on the boundary faces, namely the identifierof the element the face belongs to, the list of vertices, the flags of the boundary andof the type of boundary condition imposed. • matElem : Array containing the flag of the region each element belongs to.64 Assembly[indexGlobalFace, ⇠ ] = ...hdgElemToFaceIndex(mesh.indexTf, ...refFace, hdg.faceInfo(iElem), ...refElem(pElem).nOfFaces, ...ctt.nOfComponents);if hdg.pureDirichletiBlock = [indexGlobalFace, ...hdg.nDOFglobal+iElem, ...hdg.nDOFglobal+mesh.nOfElements+1];jBlock = [indexGlobalFace, ...hdg.nDOFglobal+iElem];elseiBlock = [indexGlobalFace, ...hdg.nDOFglobal+iElem];jBlock = iBlock;endnI =numel(iBlock);nJ =numel(jBlock);currentIndex = indexIni + (1:nJ);currentIndex2 = indexIni2 + (1:nI*nJ);for i=1:nImat.i(currentIndex) = iBlock(i);mat.j(currentIndex) = jBlock;currentIndex = currentIndex + nJ;endmat.Kij(currentIndex2) = ...reshape(Ke', 1, nI*nJ);if hdg.pureDirichletf(hdg.rowsGlobalSystem) = ...-f(hdg.rowsGlobalSystem)+intPressure;K = [K, [K(end,:)'; 0]];end%% Node permutation from GMSH to Matlab ...orderingpermNodes = nodeRenumberingMSH(pDegree, nsd);T(:, 1:nOfElemNodes) = T(:, permNodes); Figure 47: Extract of the scanMeshFileMSH function that constructs the permutationvector for the appropriate renumbering of the mesh nodes.A detailed description of the structure of the .msh file as well as the notation utilised bythe mesh generator to identify element and face types is available in the official
Gmsh doc-umentation [13]. It is worth noting that the numbering of the nodes in
HDGlab differs fromthe one provided by
Gmsh . Hence, an appropriate permutation is performed as reported infigure 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 thedata structures intFaces and extFaces , respectively.To construct the structure intFaces , the getInternalFaces function explores themesh and, for each face of each element, it identifies the neighbouring element by comparingthe list of face nodes, see figure 48. In order to optimise this operation, a list of potentialneighbours 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 con-structing 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 infor-mation obtained by the mesh generator, namely X , T and matElem , and the structures ofthe interior and exterior faces, intFaces and extFaces , previously generated. More pre-cisely, the connectivity matrix is constructed via a loop on the mesh elements as reportedin figure 50.It is worth noting that the framework provided in HDGlab to import meshes gener-ated using
Gmsh can be easily extended to any mesh generator by introducing appropriatefunctions to construct the intFaces and extFaces data structures, as described above.65 unction [jElem, jFace, jPerm] = findNeighbourElement(T, iElem, faceNodes, ...listPotentialNeigh, elemFaces, nsd, nOfElemFaces, nOfFaceVertices)sortFaceNodes = sort(faceNodes);% InitialisationjElem = 0;jFace = 0;jPerm = 0;% Loop over potential neighboursnOfPotentialNeigh = length(listPotentialNeigh(1, :));for iPotentialNeigh = 1:nOfPotentialNeighjIndexElem = listPotentialNeigh(1, iPotentialNeigh);if ⇠ (jIndexElem == iElem)% Loop over faces of potential neighbourfor jFace = 1:nOfElemFacesjFaceNodes = T(jIndexElem, elemFaces(jFace).nodes);% Compare set of nodes after sortingsortJFaceNode = sort(jFaceNodes);if sum(abs(sortFaceNodes-sortJFaceNode)) == 0jElem = jIndexElem;jPerm = getFacePermutation(jFaceNodes, faceNodes(1), nsd, nOfFaceVertices);return;endendendend Figure 48: Function findNeighbourElement that identifies the element sharing a face ofa given element. % Assembly[indexGlobalFace, ⇠ ] = ...hdgElemToFaceIndex(mesh.indexTf, ...refFace, hdg.faceInfo(iElem), ...refElem(pElem).nOfFaces, ...ctt.nOfComponents);if hdg.pureDirichletiBlock = [indexGlobalFace, ...hdg.nDOFglobal+iElem, ...hdg.nDOFglobal+mesh.nOfElements+1];jBlock = [indexGlobalFace, ...hdg.nDOFglobal+iElem];elseiBlock = [indexGlobalFace, ...hdg.nDOFglobal+iElem];jBlock = iBlock;endnI =numel(iBlock);nJ =numel(jBlock);currentIndex = indexIni + (1:nJ);currentIndex2 = indexIni2 + (1:nI*nJ);for i=1:nImat.i(currentIndex) = iBlock(i);mat.j(currentIndex) = jBlock;currentIndex = currentIndex + nJ;endmat.Kij(currentIndex2) = ...reshape(Ke', 1, nI*nJ);if hdg.pureDirichletf(hdg.rowsGlobalSystem) = ...-f(hdg.rowsGlobalSystem)+intPressure;K = [K, [K(end,:)'; 0]];end%% Node permutation from GMSH to Matlab ...orderingpermNodes = nodeRenumberingMSH(pDegree, nsd);T(:, 1:nOfElemNodes) = T(:, permNodes); %% Store the list of elements containing ...each nodenodeToElem = zeros(nOfNodes, ...nOfMaxElemsOneNode);indexNode = zeros(nOfNodes, 1);for kElem = 1:nOfElementsTe = T(kElem, :);for iLocalNode = 1:nOfElemNodesiGlobalNode = Te(iLocalNode);indexNode(iGlobalNode) = ...indexNode(iGlobalNode) + 1;nodeToElem(iGlobalNode, ...indexNode(iGlobalNode)) = kElem;endend Figure 49: Extract of the getInternalFaces function that constructs the list of elementscontaining each node.
C.4 Some examples of high-order meshes using
Gmsh
In this section, several meshes of the ring domain introduced in section 11.2 are gen-erated using
Gmsh and tested to verify the optimal convergence rate of the code whenhigh-order meshes with equally-spaced nodal sets are considered. More precisely, fig-ure 51 displays the comparison of the third order meshes provided by
Gmsh using the
Mesh.HighOrderOptimize option either with an elastic approach or an optimisation algo-rithm [13]. It is worth noting that the elastic approach mainly induces the curvature ofthe boundary of the domain, whereas the optimisation strategy is responsible for curved66
Assembly[indexGlobalFace, ⇠ ] = ...hdgElemToFaceIndex(mesh.indexTf, ...refFace, hdg.faceInfo(iElem), ...refElem(pElem).nOfFaces, ...ctt.nOfComponents);if hdg.pureDirichletiBlock = [indexGlobalFace, ...hdg.nDOFglobal+iElem, ...hdg.nDOFglobal+mesh.nOfElements+1];jBlock = [indexGlobalFace, ...hdg.nDOFglobal+iElem];elseiBlock = [indexGlobalFace, ...hdg.nDOFglobal+iElem];jBlock = iBlock;endnI =numel(iBlock);nJ =numel(jBlock);currentIndex = indexIni + (1:nJ);currentIndex2 = indexIni2 + (1:nI*nJ);for i=1:nImat.i(currentIndex) = iBlock(i);mat.j(currentIndex) = jBlock;currentIndex = currentIndex + nJ;endmat.Kij(currentIndex2) = ...reshape(Ke', 1, nI*nJ);if hdg.pureDirichletf(hdg.rowsGlobalSystem) = ...-f(hdg.rowsGlobalSystem)+intPressure;K = [K, [K(end,:)'; 0]];end%% Node permutation from GMSH to Matlab ...orderingpermNodes = nodeRenumberingMSH(pDegree, nsd);T(:, 1:nOfElemNodes) = T(:, permNodes); %% Store the list of elements containing ...each nodenodeToElem = zeros(nOfNodes, ...nOfMaxElemsOneNode);indexNode = zeros(nOfNodes, 1);for kElem = 1:nOfElementsTe = T(kElem, :);for iLocalNode = 1:nOfElemNodesiGlobalNode = Te(iLocalNode);indexNode(iGlobalNode) = ...indexNode(iGlobalNode) + 1;nodeToElem(iGlobalNode, ...indexNode(iGlobalNode)) = kElem;endend%% Construct the connectivity matrixindexIni = 1;for iElem = 1:mesh.nOfElementsindexEnd = indexIni + refElem.nOfNodes ...- 1;mesh.indexT(iElem, 1) = indexIni;mesh.indexT(iElem, 2) = indexEnd;mesh.X(indexIni:indexEnd, :) = ...X(T(iElem,:), :);indexIni = indexEnd + 1;endmesh.nOfNodes = size(mesh.X,1); Figure 50: Extract of the buildMeshStruct function that constructs the connectivitymatrix.edges to appear in the interior of the domain as well (Fig. 51a-b). In both cases, the nodesare equally-spaced.In addition, the L (Ω) 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 thecharacteristic 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 post-processed variable u (cid:63) , showing the capability of HDGlab to work using both equally-spacedand Fekete nodal distributions. It is worth noting that the results obtained using Feketenodal sets (Fig. 37) provides up to one order of magnitude of extra accuracy for a givenmesh size compared to the equally-spaced nodes in figure 51.
Acknowledgements
This work was partially supported by the Spanish Ministry of Economy and Competitive-ness (Grant number: DPI2017-85139-C2-2-R). M.G. and A.H. are also grateful for the sup-port provided by the Spanish Ministry of Economy and Competitiveness through the SeveroOchoa programme for centres of excellence in RTD (Grant number: CEX2018-000797-S)and the Generalitat de Catalunya (Grant number: 2017-SGR-1278). R.S. also acknowl-edges the support of the Engineering and Physical Sciences Research Council (Grant num-ber: EP/T009071/1).
Conflict of interest
The authors declare that they have no conflict of interest.67 a) Elasticity: nodal distribution (b) Optimisation: nodal distribution(c) Elasticity: h -convergence of p and L (d) Optimisation: h -convergence of p and L (e) Elasticity: h -convergence of u and u (cid:63) (f) Optimisation: h -convergence of u and u (cid:63) Figure 51: Comparison of high-order meshes generated by
Gmsh using the elastic approach(left) and the optimisation algorithm (right) for the Stokes equation in a ring domain.(a-b) Nodal distributions for meshes of degree 3. Convergence of the L (Ω) error of (c-d)pressure, p , and mixed variable, L , and (e-f) primal, u , and postprocessed, u (cid:63) , velocitiesas a function of the characteristic mesh size h for polynomial degree of approxiomation p =1 , . . . ,
5. 68 eferences [1] code aster : Structures and Thermomechanics Analysis for Studies and Research. https://code-aster.org/ [2]
Code Saturne : EDF open-source software to solve computational fluid dynamicsapplications. [3]
DiSk++ : A C++ template library for Discontinuous Skeletal methods. https://github.com/wareHHOuse/diskpp [4]
Feel++ : A poweful, scalable and expressive finite element embedded library in C++. [5]
FESTUNG : A MATLAB/GNU Octave toolbox for the discontinuous Galerkin method. https://github.com/FESTUNG/FESTUNG [6]
Firedrake : An automated system for the solution of partial differential equationsusing the finite element method. [7]
GetFEM : An open-source finite element library. http://getfem.org [8]
HArDCore : Hybrid Arbitrary Degree::Core. https://github.com/jdroniou/HArDCore [9]
HDG3D : Matlab implementation of the Hybridizable Discontinuous Galerkin methodon general tetrahedrizations of polyhedra in three dimensional space. https://github.com/team-pancho/HDG3D [10]
Nektar++ : Spectral/hp element framework. [11]
Netgen/NGSolve : A high performance multiphysics finite element software. https://ngsolve.org [12] deal.II : An open source finite element library. [13]
Gmsh : A three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. https://gmsh.info [14]
MFEM : Modular finite element methods library. https://mfem.org [15] Abbas, M., Ern, A., Pignet, N.: Hybrid high-order methods for finite deformationsof hyperelastic materials. Computational Mechanics (4), 909–928 (2018)[16] Abbas, M., Ern, A., Pignet, N.: A hybrid high-order method for finite elastoplas-tic deformations within a logarithmic strain framework. International Journal forNumerical Methods in Engineering (3), 303–327 (2019)6917] Abbas, M., Ern, A., Pignet, N.: A hybrid high-order method for incremental asso-ciative plasticity with small deformations. Computer Methods in Applied Mechanicsand Engineering , 891–912 (2019)[18] Agullo, E., Giraud, L., Gob´e, A., Kuhn, M., Lanteri, S., Moya, L.: High order HDGmethod and domain decomposition solvers for frequency-domain electromagnetics.International Journal of Numerical Modelling: Electronic Networks, Devices andFields (2) (2020)[19] Ainsworth, M., Fu, G.: Fully computable a posteriori error bounds for hybridizablediscontinuous Galerkin finite element approximations. Journal of Scientific Comput-ing (1), 443–466 (2018)[20] Anderson, R., Andrej, J., Barker, A., Bramwell, J., Camier, J.S., Cerveny, J., Dobrev,V., Dudouit, Y., Fisher, A., Kolev, T., Pazner, W., Stowell, M., Tomov, V., Dahm,J., Medina, D., Zampini, S.: MFEM: a modular finite element methods library. Tech.rep., arXiv (2019). URL https://arxiv.org/abs/1911.09220 [21] Anderson, T.H., Civiletti, B.J., Monk, P.B., Lakhtakia, A.: Coupled optoelectronicsimulation and optimization of thin-film photovoltaic solar cells. Journal of Compu-tational Physics , 109,242 (2020)[22] Araya, R., Solano, M., Vega, P.: Analysis of an adaptive HDG method for theBrinkman problem. IMA Journal of Numerical Analysis (3), 1502–1528 (2019)[23] Araya, R., Solano, M., Vega, P.: A posteriori error analysis of an HDG method forthe Oseen problem. Applied Numerical Mathematics , 291–308 (2019)[24] Arbogast, T., Pencheva, G., Wheeler, M.F., Yotov, I.: A Multiscale Mortar MixedFinite Element Method. Multiscale Modeling & Simulation (1), 319–346 (2007)[25] Bangerth, W., Hartmann, R., Kanschat, G.: deal.II – a general purpose objectoriented finite element library. ACM Transactions on Mathematical Software (4),24/1–24/27 (2007)[26] Barrenechea, G.R., Bosy, M., Dolean, V., Nataf, F., Tournier, P.H.: Hybrid discon-tinuous Galerkin discretisation and domain decomposition preconditioners for theStokes problem. Computational Methods in Applied Mathematics (2018)[27] Bonaldi, F., Di Pietro, D.A., Geymonat, G., Krasucki, F.: A hybrid high-ordermethod for Kirchhoff-Love plate bending problems. ESAIM: Mathematical Modellingand Numerical Analysis (2), 393–421 (2018)[28] Bonnasse-Gahot, M., Calandra, H., Diaz, J., Lanteri, S.: Hybridizable discontinuousGalerkin method for the 2-D frequency-domain elastic wave equations. GeophysicalJournal International (1), 637–659 (2017)7029] Botti, L., Di Pietro, D.A.: Assessment of hybrid high-order methods on curvedmeshes and comparison with discontinuous Galerkin methods. Journal of Computa-tional Physics , 58–84 (2018)[30] Botti, L., Di Pietro, D.A., Droniou, J.: A hybrid high-order discretisation of theBrinkman problem robust in the Darcy and Stokes limits. Computer Methods inApplied Mechanics and Engineering , 278–310 (2018)[31] Botti, L., Di Pietro, D.A., Droniou, J.: A Hybrid High-Order method for the incom-pressible Navier-Stokes equations based on Temam’s device. Journal of Computa-tional Physics , 786–816 (2019)[32] Botti, M., Di Pietro, D.A., Guglielmana, A.: A low-order nonconforming method forlinear elasticity on general meshes. Computer Methods in Applied Mechanics andEngineering , 96–118 (2019)[33] Botti, M., Di Pietro, D.A., Le Maˆıtre, O., Sochala, P.: Numerical approximation ofporoelasticity with random coefficients using polynomial chaos and hybrid high-ordermethods. Computer Methods in Applied Mechanics and Engineering , 112,736(2020)[34] Botti, M., Di Pietro, D.A., Sochala, P.: A hybrid high-order method for nonlinearelasticity. SIAM Journal on Numerical Analysis (6), 2687–2717 (2017)[35] Botti, M., Di Pietro, D.A., Sochala, P.: A hybrid high-order discretization methodfor nonlinear poroelasticity. Computational Methods in Applied Mathematics (2),227–249 (2020)[36] Brezzi, F., Fortin, M.: Mixed and hybrid finite elements methods. Springer series incomputational mathematics. Springer-Verlag (1991)[37] Bui-Thanh, T.: From Godunov to a unified hybridized discontinuous Galerkin frame-work for partial differential equations. Journal of Computational Physics , 114–146 (2015)[38] Bui-Thanh, T.: Construction and analysis of HDG methods for linearized shallowwater equations. SIAM Journal on Scientific Computing (6), A3696–A3719 (2016)[39] Burman, E., Delay, G., Ern, A.: An unfitted hybrid high-order method for the Stokesinterface problem. IMA Journal of Numerical Analysis (2020)[40] Burman, E., Ern, A.: An unfitted hybrid high-order method for elliptic interfaceproblems. SIAM Journal on Numerical Analysis (3), 1525–1546 (2018)[41] Camargo, L., L´opez-Rodr´ıguez, B., Osorio, M., Solano, M.: An HDG method forMaxwell’s equations in heterogeneous media. Computer Methods in Applied Me-chanics and Engineering , 113,178 (2020)7142] Cangiani, A., Dong, Z., Georgoulis, E., Houston, P.: hp -Version DiscontinuousGalerkin Methods on Polygonal and Polyhedral Meshes. Springer International Pub-lishing (2017)[43] Cantwell, C.D., Moxey, D., Comerford, A., Bolis, A., Rocco, G., Mengaldo, G.,De Grazia, D., Yakovlev, S., Lombard, J.E., Ekelschot, D., Jordi, B., Xu, H., Mo-hamied, Y., Eskilsson, C., Nelson, B., Vos, P., Biotto, C., Kirby, R.M., Sherwin,S.J.: Nektar++: An open-source spectral/hp element framework. Computer PhysicsCommunications , 205 – 219 (2015)[44] Cascavita, K.L., Bleyer, J., Chateau, X., Ern, A.: Hybrid discretization methodswith adaptive yield surface detection for Bingham pipe flows. Journal of ScientificComputing (3), 1424–1443 (2018)[45] Castanon Quiroz, D., Di Pietro, D.A.: A Hybrid High-Order method for the incom-pressible Navier-Stokes problem robust for large irrotational body forces. Computersand Mathematics with Applications (9), 2655–2677 (2020)[46] Castillo, P., G´omez, S.: Conservative super-convergent and hybrid discontinuousGalerkin methods applied to nonlinear Schr¨odinger equations. Applied Mathematicsand Computation , 124,950 (2020)[47] Celiker, F., Cockburn, B., Shi, K.: Hybridizable discontinuous Galerkin methods forTimoshenko beams. Journal of Scientific Computing (1), 1–37 (2010)[48] Celiker, F., Cockburn, B., Shi, K.: A projection-based error analysis of HDG methodsfor Timoshenko beams. Mathematics of Computation (277), 131–151 (2011)[49] Cesmelioglu, A., Cockburn, B., Nguyen, N.C., Peraire, J.: Analysis of HDG methodsfor Oseen equations. Journal of Scientific Computing (2), 392–431 (2013)[50] Cesmelioglu, A., Cockburn, B., Qiu, W.: Analysis of a hybridizable discontinuousGalerkin method for the steady-state incompressible Navier-Stokes equations. Math-ematics of Computation (306), 1643–1670 (2017)[51] Cesmelioglu, A., Rhebergen, S., Wells, G.N.: An embedded-hybridized discontinuousGalerkin method for the coupled Stokes-Darcy system. Journal of Computationaland Applied Mathematics , 112,476 (2020)[52] Chave, F., Di Pietro, D.A., Formaggia, L.: A hybrid high-order method for passivetransport in fractured porous media. GEM - International Journal on Geomathe-matics (1) (2019)[53] Chave, F., Di Pietro, D.A., Marche, F., Pigeonneaux, F.: A Hybrid High-ordermethod for the Cahn-Hilliard problem in mixed form. SIAM Journal on NumericalAnalysis (3), 1873–1898 (2016) 7254] Chen, G., Cockburn, B., Singler, J., Zhang, Y.: Superconvergent Interpolatory HDGMethods for Reaction Diffusion Equations I: An HDG k Method. Journal of ScientificComputing (3), 2188–2212 (2019)[55] Chen, G., Cui, J., Xu, L.: Analysis of a hybridizable discontinuous Galerkin methodfor the Maxwell operator. ESAIM: Mathematical Modelling and Numerical Analysis (1), 301–324 (2019)[56] Chen, G., Monk, P., Zhang, Y.: An HDG method for the time-dependent drift-diffusion model of semiconductor devices. Journal of Scientific Computing (2019)[57] Chen, H., Li, J., Qiu, W.: Robust a posteriori error estimates for HDG method forconvection-diffusion equations. IMA Journal of Numerical Analysis (1), 437–462(2014)[58] Chen, H., Qiu, W., Shi, K.: A priori and computable a posteriori error estimates foran HDG method for the coercive Maxwell equations. Computer Methods in AppliedMechanics and Engineering , 287–310 (2018)[59] Chen, H., Qiu, W., Shi, K., Solano, M.: A Superconvergent HDG Method for theMaxwell Equations. Journal of Scientific Computing (3), 1010–1029 (2017)[60] Chen, Y., Cockburn, B.: Analysis of variable-degree HDG methods for convection-diffusion equations. Part I: General nonconforming meshes. IMA Journal of Numer-ical Analysis (4), 1267–1293 (2012)[61] Chen, Y., Cockburn, B.: Analysis of variable-degree HDG methods for convection-diffusion equations. Part II: Semimatching nonconforming meshes. Mathematics ofComputation (285), 87–111 (2014)[62] Chen, Y., Cockburn, B., Dong, B.: Superconvergent HDG methods for linear, sta-tionary, third-order equations in one-space dimension. Mathematics of Computation (302), 2715–2742 (2016)[63] Chen, Y., Dong, B., Jiang, J.: Optimally convergent hybridizable discontinuousGalerkin method for fifth-order Korteweg-de Vries type equations. ESAIM: Mathe-matical Modelling and Numerical Analysis (6), 2283–2306 (2018)[64] Childs, P.R.: Rotating flow. Elsevier (2010)[65] Cho, K., Moon, M.: Multiscale hybridizable discontinuous Galerkin method for el-liptic problems in perforated domains. Journal of Computational and Applied Math-ematics , 112,346 (2020)[66] Chouly, F., Ern, A., Pignet, N.: A Hybrid High-Order discretization combined withNitsche’s method for contact and Tresca friction in small strain elasticity. SIAMJournal on Scientific Computing (4), A2300–A2324 (2020)7367] Christophe, A., Descombes, S., Lanteri, S.: An implicit hybridized discontinuousGalerkin method for the 3D time-domain Maxwell equations. Applied Mathematicsand Computation , 395–408 (2018)[68] Chung, E., Cockburn, B., Fu, G.: The staggered DG method is the limit of a hy-bridizable DG method. SIAM Journal on Numerical Analysis (2), 915–932 (2014)[69] Chung, E., Cockburn, B., Fu, G.: The Staggered DG Method is the Limit of aHybridizable DG Method. Part II: The Stokes Flow. Journal of Scientific Computing (2), 870–887 (2016)[70] Chung, E., Efendiev, Y., Leung, W.T.: Generalized multiscale finite element methodswith energy minimizing oversampling. International Journal for Numerical Methodsin Engineering (3), 316–343 (2019)[71] Cicuttin, M., Di Pietro, D., Ern, A.: Implementation of discontinuous skeletal meth-ods on arbitrary-dimensional, polytopal meshes using generic programming. Journalof Computational and Applied Mathematics , 852 – 874 (2018)[72] Cicuttin, M., Ern, A., Gudi, T.: Hybrid high-order methods for the elliptic obstacleproblem. Journal of Scientific Computing (1), 8 (2020)[73] Cicuttin, M., Ern, A., Lemaire, S.: A hybrid high-order method for highly oscillatoryelliptic problems. Computational Methods in Applied Mathematics (2018)[74] Ciucˇa, C., Fernandez, P., Christophe, A., Nguyen, N.C., Peraire, J.: Implicit hy-bridized discontinuous Galerkin methods for compressible magnetohydrodynamics.Journal of Computational Physics: X , 100,042 (2020)[75] Cockburn, B.: Static condensation, hybridization, and the devising of the HDGmethods. In: G.R. Barrenechea, F. Brezzi, A. Cangiani, E. Georgoulis (eds.) BuildingBridges: Connections and Challenges in Modern Approaches to Numerical PartialDifferential Equations, pp. 129–177. Springer International Publishing, Cham (2016)[76] Cockburn, B.., Mustapha, K.: A hybridizable discontinuous Galerkin method forfractional diffusion problems. Numerische Mathematik (2), 293–314 (2015)[77] Cockburn, B., Dong, B., Guzm´an, J.: A superconvergent LDG-hybridizable Galerkinmethod for second-order elliptic problems. Mathematics of Computation (264),1887–1916 (2008)[78] Cockburn, B., Dong, B., Guzm´an, J.: A hybridizable and superconvergent discon-tinuous Galerkin method for biharmonic problems. Journal of Scientific Computing (1-3), 141–187 (2009) 7479] Cockburn, B., Dong, B., Guzm´an, J., Restelli, M., Sacco, R.: A hybridizable dis-continuous Galerkin method for steady-state convection-diffusion-reaction problems.SIAM Journal on Scientific Computing (5), 3827–3846 (2009)[80] Cockburn, B., Dubois, O., Gopalakrishnan, J., Tan, S.: Multigrid for an HDGmethod. IMA Journal of Numerical Analysis (4), 1386–1425 (2014)[81] Cockburn, B., Fu, G.: Devising superconvergent HDG methods with symmetricapproximate stresses for linear elasticity by M -decompositions. IMA Journal ofNumerical Analysis (2), 566–604 (2017)[82] Cockburn, B., Fu, G.: Superconvergence by M -decompositions. Part II: Constructionof two-dimensional finite elements. ESAIM: Mathematical Modelling and NumericalAnalysis (1), 165–186 (2017)[83] Cockburn, B., Fu, G.: Superconvergence by M -decompositions. Part III: Construc-tion of three-dimensional finite elements. ESAIM: Mathematical Modelling and Nu-merical Analysis (1), 365–398 (2017)[84] Cockburn, B., Fu, G., Qiu, W.: A note on the devising of superconvergent HDGmethods for Stokes flow by M -decompositions. IMA Journal of Numerical Analysis (2), 730–749 (2017)[85] Cockburn, B., Fu, G., Sayas, F.J.: Superconvergence by M -decompositions. PartI: General theory for HDG methods for diffusion. Mathematics of Computation (306), 1609–1641 (2017)[86] Cockburn, B., Fu, Z., Hungria, A., Ji, L., S´anchez, M.A., Sayas, F.J.: Stormer-Numerov HDG Methods for Acoustic Waves. Journal of Scientific Computing (2),597–624 (2018)[87] Cockburn, B., Gopalakrishnan, J.: A characterization of hybridized mixed methodsfor second order elliptic problems. SIAM Journal on Numerical Analysis (1), 283–301 (2004)[88] Cockburn, B., Gopalakrishnan, J.: The derivation of hybridizable discontinuousGalerkin methods for Stokes flow. SIAM Journal on Numerical Analysis (2),1092–1125 (2009)[89] Cockburn, B., Gopalakrishnan, J., Lazarov, R.: Unified hybridization of discontin-uous Galerkin, mixed, and continuous Galerkin methods for second order ellipticproblems. SIAM Journal on Numerical Analysis (2), 1319–1365 (2009)[90] Cockburn, B., Gopalakrishnan, J., Nguyen, N.C., Peraire, J., Sayas, F.J.: Analysisof HDG methods for Stokes flow. Mathematics of Computation (274), 723–760(2011) 7591] Cockburn, B., Karniadakis, G.E., Shu, C.W.: The development of discontinuousGalerkin methods. In: Discontinuous Galerkin methods (Newport, RI, 1999), Lect.Notes Comput. Sci. Eng. , vol. 11, pp. 3–50. Springer, Berlin (2000)[92] Cockburn, B., Nguyen, N.C., Peraire, J.: A comparison of HDG methods for Stokesflow. Journal of Scientific Computing (1-3), 215–237 (2010)[93] Cockburn, B., Quenneville-B´elair, V.: Uniform-in-time superconvergence of the HDGmethods for the acoustic wave equation. Mathematics of Computation (285), 65–85 (2014)[94] Cockburn, B., Sayas, F.J.: Divergence-conforming HDG methods for Stokes flows.Mathematics of Computation (288), 1571–1598 (2014)[95] Cockburn, B., Shen, J.: A hybridizable discontinuous Galerkin method for the p -Laplacian. SIAM Journal on Scientific Computing (1), A545–A566 (2016)[96] Cockburn, B., Shen, J.: An algorithm for stabilizing hybridizable discontinuousGalerkin methods for nonlinear elasticity. Results in Applied Mathematics (2019)[97] Cockburn, B., Shi, K.: Superconvergent HDG methods for linear elasticity withweakly symmetric stresses. IMA Journal of Numerical Analysis (3), 747–770 (2013)[98] Cockburn, B., Shi, K.: Devising HDG methods for Stokes flow: an overview. Com-puters & Fluids , 221–229 (2014)[99] Cockburn, B., Shu, C.W.: The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM Journal on Numerical Analysis (6),2440–2463 (1998)[100] Cockburn, B., Singler, J.R., Zhang, Y.: Interpolatory HDG Method for ParabolicSemilinear PDEs. Journal of Scientific Computing (3), 1777–1800 (2019)[101] Cockburn, B., Solano, M.: Solving Dirichlet boundary-value problems on curveddomains by extensions from subdomains. SIAM Journal on Scientific Computing (1), A497–A519 (2012)[102] Cockburn, B., Solano, M.: Solving convection-diffusion problems on curved domainsby extensions from subdomains. Journal of Scientific Computing (2), 512–543(2014)[103] Cockburn, B., Wang, Z.: Adjoint-based, superconvergent Galerkin approximationsof linear functionals. Journal of Scientific Computing (2-3), 644–666 (2017)[104] Cockburn, B., Zhang, W.: A posteriori error estimates for HDG methods. Journalof Scientific Computing (3), 582–607 (2012)76105] Cockburn, B., Zhang, W.: A posteriori error analysis for hybridizable discontinuousGalerkin methods for second order elliptic problems. SIAM Journal on NumericalAnalysis (1), 676–693 (2013)[106] Cockburn, B., Di Pietro, D. A., Ern, A.: Bridging the hybrid high-order and hy-bridizable discontinuous Galerkin methods. ESAIM: Mathematical Modelling andNumerical Analysis (3), 635–650 (2016)[107] Costa-Sol´e, A., Ruiz-Giron´es, E., Sarrate, J.: An HDG formulation for incompressibleand immiscible two-phase porous media flow problems. International Journal ofComputational Fluid Dynamics (4), 137–148 (2019)[108] Di Pietro, D., Ern, A.: Mathematical aspects of discontinuous Galerkin methods,vol. 69. Springer, Heidelberg (2012)[109] Di Pietro, D., Ern, A.: A hybrid high-order locking-free method for linear elasticityon general meshes. Computer Methods in Applied Mechanics and Engineering ,1–21 (2015)[110] Di Pietro, D., Ern, A., Lemaire, S.: An arbitrary-order and compact-stencil dis-cretization of diffusion on general meshes based on local reconstruction operators.Computational Methods in Applied Mathematics (4), 461–472 (2014)[111] Di Pietro, D.A., Droniou, J.: A Hybrid High-Order method for Leray-Lions ellip-tic equations on general meshes. Mathematics of Computation (307), 2159–2191(2017)[112] Di Pietro, D.A., Droniou, J.: The Hybrid High-Order Method for Polytopal Meshes.Modeling, Simulation and Applications series. Springer International Publishing(2020)[113] Di Pietro, D.A., Droniou, J., Ern, A.: A discontinuous-skeletal method for advection-diffusion-reaction on general meshes. SIAM Journal on Numerical Analysis (5),2135–2157 (2015)[114] Di Pietro, D.A., Droniou, J., Manzini, G.: Discontinuous Skeletal Gradient Dis-cretisation methods on polytopal meshes. Journal of Computational Physics ,397–425 (2018)[115] Di Pietro, D.A., Ern, A.: Hybrid high-order methods for variable-diffusion problemson general meshes. Comptes Rendus Mathematique (1), 31–34 (2015)[116] Di Pietro, D.A., Ern, A., Linke, A., Schieweck, F.: A discontinuous skeletal methodfor the viscosity-dependent Stokes problem. Computer Methods in Applied Mechan-ics and Engineering , 175–195 (2016)77117] Di Pietro, D.A., Krell, S.: A hybrid high-order method for the steady incompressibleNavier-Stokes problem. Journal of Scientific Computing (3), 1677–1705 (2018)[118] Diskin, B., Thomas, J.L.: Comparison of node-centered and cell-centered unstruc-tured finite-volume discretizations: inviscid fluxes. AIAA journal (4), 836–854(2011)[119] Diskin, B., Thomas, J.L., Nielsen, E.J., Nishikawa, H., White, J.A.: Comparison ofnode-centered and cell-centered unstructured finite-volume discretizations: viscousfluxes. AIAA journal (7), 1326 (2010)[120] Donea, J., Huerta, A.: Finite Element Methods for Flow Problems. Finite ElementMethods for Flow Problems. John Wiley & Sons (2003)[121] Dong, H., Wang, B., Xie, Z., Wang, L.L.: An unfitted hybridizable discontinuousGalerkin method for the Poisson interface problem and its error analysis. IMAJournal of Numerical Analysis (1), 444–476 (2016)[122] Du, S., Sayas, F.J.: A unified error analysis of hybridizable discontinuous Galerkinmethods for the static Maxwell equations. SIAM Journal on Numerical Analysis (2), 1367–1391 (2020)[123] Efendiev, Y., Lazarov, R., Moon, M., Shi, K.: A spectral multiscale hybridizable dis-continuous Galerkin method for second order elliptic problems. Computer Methodsin Applied Mechanics and Engineering , 243 – 256 (2015). Special Issue on Ad-vances in Simulations of Subsurface Flow and Transport (Honoring Professor MaryF. Wheeler)[124] Egger, H., Sch¨oberl, J.: A hybrid mixed discontinuous Galerkin finite-elementmethod for convection-diffusion problems. IMA Journal of Numerical Analysis (4),1206–1234 (2009)[125] Egger, H., Waluga, C.: hp -analysis of a hybrid DG method for Stokes flow. IMAJournal of Numerical Analysis (2), 687–721 (2012)[126] Egger, H., Waluga, C.: A hybrid mortar method for incompressible flow. Interna-tional Journal of Numerical Analysis and Modeling (4), 793–812 (2012)[127] Fabien, M.S.: A GPU-accelerated hybridizable discontinuous Galerkin method forlinear elasticity. Communications in Computational Physics (2), 513–545 (2020)[128] Fabien, M.S.: A high-order implicit HDG method for the Benjamin-Bona-Mahonyequation. International Journal for Numerical Methods in Fluids (2020)[129] Fabien, M.S., Knepley, M.G., Mills, R.T., Riviere, B.M.: Manycore parallel com-puting for a hybridizable discontinuous Galerkin nested multigrid method. SIAMJournal on Scientific Computing (2), C73–C96 (2019)78130] Fabien, M.S., Knepley, M.G., Rivi`ere, B.M.: A hybridizable discontinuous Galerkinmethod for two-phase flow in heterogeneous porous media. International Journal forNumerical Methods in Engineering (3), 161–177 (2018)[131] Farahinia, A., Zhang, W.J.: Numerical investigation into the mixing performance ofmicro T-mixers with different patterns of obstacles. Journal of the Brazilian Societyof Mechanical Sciences and Engineering (11), 491 (2019)[132] Fernandez, P., Christophe, A., Terrana, S., Nguyen, N.C., Peraire, J.: Hybridizeddiscontinuous Galerkin methods for wave propagation. Journal of Scientific Comput-ing (3), 1566–1604 (2018)[133] Fernandez, P., Nguyen, N.C., Peraire, J.: The hybridized discontinuous Galerkinmethod for implicit large-eddy simulation of transitional turbulent flows. Journal ofComputational Physics , 308 – 329 (2017)[134] Fidkowski, K.J.: A hybridized discontinuous Galerkin method on mapped deformingdomains. Computers & Fluids , 80 – 91 (2016)[135] Fidkowski, K.J.: Comparison of hybrid and standard discontinuous Galerkin meth-ods in a mesh-optimisation setting. International Journal of Computational FluidDynamics (1-2), 34–42 (2019)[136] Fidkowski, K.J., Chen, G.: Output-based mesh optimization for hybridized and em-bedded discontinuous Galerkin methods. International Journal for Numerical Meth-ods in Engineering (5), 867–887 (2020)[137] Fournier, Y., Bonelle, J., Moulinec, C., Shang, Z., Sunderland, A.G., Uribe, J.C.:Optimizing Code Saturne computations on Petascale systems. Computers & Fluids (1), 103 – 108 (2011)[138] Fraeijs de Veubeke, B.: Displacement and equilibrium models in the finite elementmethod. In: O.C. Zienkiewicz and G.S. Holister (ed.) Stress Analysis, pp. 145–197.John Wiley & Sons (1965)[139] Franciolini, M., Fidkowski, K.J., Crivellini, A.: Efficient discontinuous Galerkin im-plementations and preconditioners for implicit unsteady compressible flow simula-tions. Computers & Fluids , 104,542 (2020)[140] Fu, G.: Arbitrary Lagrangian-Eulerian hybridizable discontinuous Galerkin methodsfor incompressible flow with moving boundaries and interfaces. Computer Methodsin Applied Mechanics and Engineering , 113,158 (2020)[141] Fu, G., Cockburn, B., Stolarski, H.: Analysis of an HDG method for linear elastic-ity. International Journal for Numerical Methods in Engineering (3-4), 551–575(2015) 79142] Fu, G., Jin, Y., Qiu, W.: Parameter-free superconvergent H(div)-conforming HDGmethods for the Brinkman equations. IMA Journal of Numerical Analysis (2),957–982 (2018)[143] Fu, Z., Gatica, L.F., Sayas, F.J.: Algorithm 949: MATLAB Tools for HDG in ThreeDimensions. ACM Transactions on Mathematical Software (3) (2015)[144] Gander, M.J., Hajian, S.: Analysis of Schwarz methods for a hybridizable discontinu-ous Galerkin discretization: the many-subdomain case. Mathematics of Computation (312), 1635–1657 (2018)[145] Gatica, G.N., Sequeira, F.A.: Analysis of an augmented HDG method for a classof quasi-Newtonian Stokes flows. Journal of Scientific Computing (3), 1270–1308(2015)[146] Geuzaine, C., Remacle, J.F.: Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methodsin Engineering (11), 1309–1331 (2009)[147] Giacomini, M., Borchini, L., Sevilla, R., Huerta, A.: Separated response surfacesfor flows in parametrised domains: comparison of a priori and a posteriori PGDalgorithms. Tech. rep., arXiv (2020). URL https://arxiv.org/abs/2009.02176 .Submitted.[148] Giacomini, M., Karkoulias, A., Sevilla, R., Huerta, A.: A superconvergent HDGmethod for Stokes flow with strongly enforced symmetry of the stress tensor. Journalof Scientific Computing (3), 1679–1702 (2018)[149] Giacomini, M., Sevilla, R.: Discontinuous Galerkin approximations in computationalmechanics: hybridization, exact geometry and degree adaptivity. SN Applied Sci-ences , 1047 (2019)[150] Giacomini, M., Sevilla, R.: A second-order face-centred finite volume method ongeneral meshes with automatic mesh adaptation. International Journal for NumericalMethods in Engineering (2020)[151] Giacomini, M., Sevilla, R., Huerta, A.: Tutorial on Hybridizable DiscontinuousGalerkin (HDG) Formulation for Incompressible Flow Problems. In: L.D. Loren-zis, A. D¨uster (eds.) Modeling in Engineering Using Innovative Numerical Methodsfor Solids and Fluids, CISM International Centre for Mechanical Sciences , vol. 599,pp. 163–201. Springer International Publishing (2020)[152] Giorgiani, G., Fern´andez-M´endez, S., Huerta, A.: Hybridizable discontinuousGalerkin p -adaptivity for wave propagation problems. International Journal for Nu-merical Methods in Fluids (12), 1244–1262 (2013)80153] Giorgiani, G., Fern´andez-M´endez, S., Huerta, A.: Hybridizable discontinuousGalerkin with degree adaptivity for the incompressible Navier–Stokes equations.Computers & Fluids , 196–208 (2014)[154] G¨urkan, C., Kronbichler, M., Fern´andez-M´endez, S.: eXtended Hybridizable Discon-tinuous Galerkin with Heaviside enrichment for heat bimaterial problems. Journalof Scientific Computing (2), 542–567 (2017)[155] G¨urkan, C., Kronbichler, M., Fern´andez-M´endez, S.: eXtended hybridizable discon-tinuous Galerkin for incompressible flow problems with unfitted meshes and inter-faces. International Journal for Numerical Methods in Engineering (7), 756–777(2019)[156] G¨urkan, C., Sala-Lardies, E., Kronbichler, M., Fern´andez-M´endez, S.: eXtendedHybridizable Discontinous Galerkin (X-HDG) for void problems. Journal of ScientificComputing (3), 1313–1333 (2016)[157] Guyan, R.: Reduction of stiffness and mass matrices. AIAA Journal (2), 380–380(1965)[158] Hesthaven, J., Warburton, T.: Nodal High-Order Methods on Unstructured Grids:I. Time-Domain Solution of Maxwell’s Equations. Journal of Computational Physics (1), 186–221 (2002)[159] Hoermann, J.M., Bertoglio, C., Kronbichler, M., Pfaller, M.R., Chabiniok, R., Wall,W.A.: An adaptive hybridizable discontinuous Galerkin approach for cardiac electro-physiology. International Journal for Numerical Methods in Biomedical Engineering (5) (2018)[160] Horv´ath, T.L., Rhebergen, S.: A locally conservative and energy-stable finite-elementmethod for the Navier-Stokes problem on time-dependent domains. InternationalJournal for Numerical Methods in Fluids (12), 519–532 (2019)[161] Horv´ath, T.L., Rhebergen, S.: An exactly mass conserving space-time embedded-hybridized discontinuous Galerkin method for the Navier-Stokes equations on movingdomains. Journal of Computational Physics , 109,577 (2020)[162] Huang, J., Huang, X.: A hybridizable discontinuous Galerkin method for Kirchhoffplates. Journal of Scientific Computing (1), 290–320 (2019)[163] Huerta, A., Angeloski, A., Roca, X., Peraire, J.: Efficiency of high-order elements forcontinuous and discontinuous Galerkin methods. International Journal for NumericalMethods in Engineering (9), 529–560 (2013)[164] Hungria, A., Prada, D., Sayas, F.J.: HDG methods for elastodynamics. Computers& Mathematics with Applications (11), 2671 – 2690 (2017)81165] Huynh, L.N.T., Nguyen, N.C., Peraire, J., Khoo, B.C.: A high-order hybridizablediscontinuous Galerkin method for elliptic interface problems. International Journalfor Numerical Methods in Engineering (2), 183–200 (2013)[166] Jaust, A., Reuter, B., Aizinger, V., Sch¨utz, J., Knabner, P.: FESTUNG: A MAT-LAB/GNU Octave toolbox for the discontinuous Galerkin method. Part III: Hy-bridized discontinuous Galerkin (HDG) formulation. Computers & Mathematicswith Applications (12), 4505 – 4533 (2018)[167] Kabaria, H., Lew, A.J., Cockburn, B.: A hybridizable discontinuous Galerkin for-mulation for non-linear elasticity. Computer Methods in Applied Mechanics andEngineering , 303–329 (2015)[168] Kang, S., Bui-Thanh, T., Arbogast, T.: A hybridized discontinuous Galerkin methodfor a linear degenerate elliptic equation arising from two-phase mixtures. ComputerMethods in Applied Mechanics and Engineering , 315–336 (2019)[169] Kang, S., Giraldo, F.X., Bui-Thanh, T.: IMEX HDG-DG: A coupled implicit hy-bridized discontinuous Galerkin and explicit discontinuous Galerkin approach forshallow water systems. Journal of Computational Physics , 109,010 (2020)[170] Kirby, R., Sherwin, S.J., Cockburn, B.: To CG or to HDG: A comparative study.Journal of Scientific Computing (1), 183–212 (2011)[171] Kirk, K.L.A., Rhebergen, S.: Analysis of a pressure-robust hybridized discontinuousGalerkin method for the stationary Navier-Stokes equations. Journal of ScientificComputing (2), 881–897 (2019)[172] Komala-Sheshachala, S., Sevilla, R., Hassan, O.: A coupled HDG-FV scheme for thesimulation of transient inviscid compressible flows. Computers & Fluids , 104,495(2020)[173] Kronbichler, M., Schoeder, S., M¨uller, C., Wall, W.A.: Comparison of implicit andexplicit hybridizable discontinuous Galerkin methods for the acoustic wave equation.International Journal for Numerical Methods in Engineering (9), 712–739 (2016)[174] Kronbichler, M., Wall, W.A.: A performance comparison of continuous and discon-tinuous Galerkin methods with fast multigrid solvers. SIAM Journal on ScientificComputing (5), A3423–A3448 (2018)[175] La Spina, A., Giacomini, M., Huerta, A.: Hybrid coupling of CG and HDG dis-cretizations based on Nitsche’s method. Computational Mechanics (2), 311–330(2020) 82176] La Spina, A., Kronbichler, M., Giacomini, M., Wall, W., Huerta, A.: A weaklycompressible hybridizable discontinuous Galerkin formulation for fluid-structure in-teraction problems. Computer Methods in Applied Mechanics and Engineering ,113,392 (2020)[177] Lederer, P.L., Lehrenfeld, C., Sch¨oberl, J.: Hybrid discontinuous Galerkin methodswith relaxed H (div)-conformity for incompressible flows. Part I. SIAM Journal onNumerical Analysis (4), 2070–2094 (2018)[178] Lederer, P.L., Lehrenfeld, C., Sch¨oberl, J.: Hybrid discontinuous Galerkin methodswith relaxed H (div)-conformity for incompressible flows. Part II. ESAIM: Mathe-matical Modelling and Numerical Analysis (2), 503–522 (2019)[179] Lederer, P.L., Lehrenfeld, C., Sch¨oberl, J.: Divergence-free tangential finite elementmethods for incompressible flows on surfaces. International Journal for NumericalMethods in Engineering (11), 2503–2533 (2020)[180] Lee, J.J., Shannon, S.J., Bui-Thanh, T., Shadid, J.N.: Analysis of an HDG methodfor linearized incompressible resistive MHD equations. SIAM Journal on NumericalAnalysis (4), 1697–1722 (2019)[181] Lehrenfeld, C., Sch¨oberl, J.: High order exactly divergence-free hybrid discontinuousGalerkin methods for unsteady incompressible flows. Computer Methods in AppliedMechanics and Engineering , 339–361 (2016)[182] Leng, H., Chen, Y.: Adaptive hybridizable discontinuous Galerkin methods for non-stationary convection-diffusion problems. Advances in Computational Mathematics (4) (2020)[183] Li, G., Shi, K.: Upscaled HDG methods for Brinkman equations with high-contrastheterogeneous coefficient. Journal of Scientific Computing (3), 1780–1800 (2018)[184] Li, L., Lanteri, S., Mortensen, N.A., Wubs, M.: A hybridizable discontinuousGalerkin method for solving nonlocal optical response models. Computer PhysicsCommunications , 99–107 (2017)[185] Li, L., Lanteri, S., Perrussel, R.: A hybridizable discontinuous Galerkin methodcombined to a Schwarz algorithm for the solution of 3D time-harmonic Maxwell’sequation. Journal of Computational Physics , 563–581 (2014)[186] Li, L., Lanteri, S., Perrussel, R.: A class of locally well-posed hybridizable discon-tinuous Galerkin methods for the solution of time-harmonic Maxwell’s equations.Computer Physics Communications , 23 – 31 (2015)[187] Liu, Y.: Fast multipole boundary element method: theory and applications in engi-neering. Cambridge university press (2009)83188] Loseille, A., Feuillet, R.: Vizir: High-order mesh and solution visualization usingOpenGL 4.0 graphic pipeline. In: 2018 AIAA Aerospace Sciences Meeting, p. 1174(2018)[189] Lu, P.., Chen, H., Qiu, W.: An absolutely stable hp-HDG method for the time-harmonic Maxwell equations with high wave number. Mathematics of Computation (306), 1553–1577 (2017)[190] McLachlan, R.I., Stern, A.: Multisymplecticity of hybridizable discontinuousGalerkin methods. Foundations of Computational Mathematics (1), 35–69 (2020)[191] Montlaur, A., Fern´andez-M´endez, S., Huerta, A.: Discontinuous Galerkin methodsfor the Stokes equations using divergence-free approximations. International Journalfor Numerical Methods in Fluids (9), 1071–1092 (2008)[192] Moon, M., Lazarov, R., Jun, H.K.: Multiscale HDG model reduction method forflows in heterogeneous porous media. Applied Numerical Mathematics , 115–133(2019)[193] Moro, D., Nguyen, N.C., Peraire, J.: Navier-Stokes solution using hybridizable dis-continuous Galerkin methods. In: 20th AIAA Computational Fluid Dynamics Con-ference. AIAA (2011)[194] Muix´ı, A., Rodr´ıguez-Ferran, A., Fern´andez-M´endez, S.: A hybridizable discontinu-ous Galerkin phase-field model for brittle fracture with adaptive refinement. Inter-national Journal for Numerical Methods in Engineering (6), 1147–1169 (2020)[195] Muralikrishnan, S., Bui-Thanh, T., Shadid, J.N.: A multilevel approach for tracesystem in HDG discretizations. Journal of Computational Physics , 109,240(2020)[196] Muralikrishnan, S., Tran, M., Bui-Thanh, T.: An improved iterative HDG approachfor partial differential equations. Journal of Computational Physics , 295 – 321(2018)[197] Mustapha, K., Nour, M., Cockburn, B.: Convergence and superconvergence analysesof HDG methods for time fractional diffusion problems. Advances in ComputationalMathematics (2), 377–393 (2016)[198] Nelson, B., Liu, E., Kirby, R.M., Haimes, R.: Elvis: A system for the accurate andinteractive visualization of high-order finite element solutions. IEEE transactions onvisualization and computer graphics (12), 2325–2334 (2012)[199] Nguyen, N., Peraire, J., Cockburn, B.: A hybridizable discontinuous Galerkinmethod for Stokes flow. Computer Methods in Applied Mechanics and Engineer-ing (9-12), 582–597 (2010) 84200] Nguyen, N.C., Peraire, J., Cockburn, B.: An implicit high-order hybridizable dis-continuous Galerkin method for linear convection-diffusion equations. Journal ofComputational Physics (9), 3232–3254 (2009)[201] Nguyen, N.C., Peraire, J., Cockburn, B.: An implicit high-order hybridizable dis-continuous Galerkin method for nonlinear convection-diffusion equations. Journal ofComputational Physics (23), 8841–8855 (2009)[202] Nguyen, N.C., Peraire, J., Cockburn, B.: High-order implicit hybridizable discontin-uous Galerkin methods for acoustics and elastodynamics. Journal of ComputationalPhysics (10), 3695–3718 (2011)[203] Nguyen, N.C., Peraire, J., Cockburn, B.: Hybridizable discontinuous Galerkin meth-ods for the time-harmonic Maxwell’s equations. Journal of Computational Physics (19), 7151–7175 (2011)[204] Nguyen, N.C., Peraire, J., Cockburn, B.: An implicit high-order hybridizable discon-tinuous Galerkin method for the incompressible Navier-Stokes equations. Journal ofComputational Physics (4), 1147–1170 (2011)[205] Nguyen, N.C., Peraire, J., Cockburn, B.: A class of embedded discontinuous Galerkinmethods for computational fluid dynamics. Journal of Computational Physics ,674–692 (2015)[206] Oikawa, I.: A hybridized discontinuous Galerkin method with reduced stabilization.Journal of Scientific Computing (1), 327–340 (2015)[207] Oikawa, I.: Analysis of a reduced-order HDG method for the Stokes equations. Jour-nal of Scientific Computing (2), 475–492 (2016)[208] Paipuri, M., Tiago, C., Fern´andez-M´endez, S.: Coupling of continuous and hybridiz-able discontinuous Galerkin methods: Application to conjugate heat transfer prob-lem. Journal of Scientific Computing (1), 321–350 (2019)[209] Peraire, J., Nguyen, N.C., Cockburn, B.: A hybridizable discontinuous Galerkinmethod for the compressible Euler and Navier-Stokes equations. AIAA paper ,2010 (2010)[210] Peters, E., Evans, J.: A divergence-conforming hybridized discontinuous Galerkinmethod for the incompressible Reynolds-averaged Navier-Stokes equations. Interna-tional Journal for Numerical Methods in Fluids , 112–133 (2019)[211] Pignet, N.: Hybrid High-Order methods for nonlinear solid mechanics. PhD thesis,Universit´e Paris-Est Marne la Vall´ee (2019). TEL 0231815785212] Poya, R., Sevilla, R., Gil, A.J.: A unified approach for a posteriori high-order curvedmesh generation using solid mechanics. Computational Mechanics (3), 457–490(2016)[213] Prud’homme, C.: A Domain Specific Embedded Language in C++ for AutomaticDifferentiation, Projection, Integration and Variational Formulations. Scientific Pro-gramming , 150,736 (2006)[214] Qiu, W., Shen, J., Shi, K.: An HDG method for linear elasticity with strong sym-metric stresses. Mathematics of Computation (309), 69–93 (2018)[215] Qiu, W., Shi, K.: A superconvergent HDG method for the incompressible Navier-Stokes equations on general polyhedral meshes. IMA Journal of Numerical Analysis (4), 1943–1967 (2016)[216] Qiu, W., Shi, K.: Analysis on an HDG method for the p -Laplacian equations. Journalof Scientific Computing (2), 1019–1032 (2019)[217] Qiu, W., Solano, M., Vega, P.: A high order HDG method for curved-interface prob-lems via approximations from straight triangulations. Journal of Scientific Comput-ing (3), 1384–1407 (2016)[218] Quarteroni, A.: Numerical models for differential problems, MS&A. Modeling, Sim-ulation and Applications , vol. 16. Springer, Cham (2017)[219] Rathgeber, F., Ham, D.A., Mitchell, L., Lange, M., Luporini, F., Mcrae, A.T.T.,Bercea, G.T., Markall, G.R., Kelly, P.H.J.: Firedrake: Automating the finite elementmethod by composing abstractions. ACM Transactions on Mathematical Software (3) (2016)[220] Remacle, J.F., Chevaugeon, N., Marchandise, E., Geuzaine, C.: Efficient visualiza-tion of high-order finite elements. International Journal for Numerical Methods inEngineering (4), 750–771 (2007)[221] Renard, Y., Poulios, K.: GetFEM: Automated FE modeling of multiphysics problemsbased on a generic weak form language. Tech. rep., HAL (2020). URL https://hal.archives-ouvertes.fr/hal-02532422 [222] Rhebergen, S., Cockburn, B.: A space-time hybridizable discontinuous Galerkinmethod for incompressible flows on deforming domains. Journal of ComputationalPhysics (11), 4185–4204 (2012)[223] Rhebergen, S., Wells, G.: A hybridizable discontinuous Galerkin method for theNavier-Stokes equations with pointwise divergence-free velocity field. Journal of Sci-entific Computing (3), 1484–1501 (2018)86224] Rhebergen, S., Wells, G.: Preconditioning of a hybridized discontinuous Galerkinfinite element method for the Stokes equations. Journal of Scientific Computing (3), 1936–1952 (2018)[225] Rhebergen, S., Wells, G.N.: An embedded-hybridized discontinuous Galerkin finiteelement method for the Stokes equations. Computer Methods in Applied Mechanicsand Engineering , 112,619 (2020)[226] Rivi`ere, B.: Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equa-tions. Society for Industrial and Applied Mathematics (2008)[227] Rocha, B.M., dos Santos, R. W., Igreja, I., Loula, A.F.D.: Stabilized hybrid dis-continuous Galerkin finite element method for the cardiac monodomain equation.International Journal for Numerical Methods in Biomedical Engineering (7) (2020)[228] Samii, A., Dawson, C.: An explicit hybridized discontinuous Galerkin method forSerre-Green-Naghdi wave model. Computer Methods in Applied Mechanics andEngineering , 447 – 470 (2018)[229] Samii, A., Kazhyken, K., Michoski, C., Dawson, C.: A comparison of the explicitand implicit hybridizable discontinuous Galerkin methods for nonlinear shallow waterequations. Journal of Scientific Computing (3), 1936–1956 (2019)[230] Samii, A., Michoski, C., Dawson, C.: A parallel and adaptive hybridized discontinu-ous Galerkin method for anisotropic nonhomogeneous diffusion. Computer Methodsin Applied Mechanics and Engineering , 118 – 139 (2016)[231] Samii, A., Panda, N., Michoski, C., Dawson, C.: A hybridized discontinuous Galerkinmethod for the nonlinear Korteweg-de Vries equation. Journal of Scientific Comput-ing (1), 191–212 (2016)[232] S´anchez, M.A., Ciuca, C., Nguyen, N.C., Peraire, J., Cockburn, B.: SymplecticHamiltonian HDG methods for wave propagation phenomena. Journal of Computa-tional Physics , 951–973 (2017)[233] S´anchez-Vizuet, T., Solano, M.E.: A hybridizable discontinuous Galerkin solver forthe Grad-Shafranov equation. Computer Physics Communications , 120–132(2019)[234] S´anchez-Vizuet, T., Solano, M.E., Cerfon, A.J.: Adaptive Hybridizable Discontinu-ous Galerkin discretization of the Grad-Shafranov equation by extension from polyg-onal subdomains. Computer Physics Communications , 107,239 (2020)[235] Sch¨oberl, J.: C++11 Implementation of Finite Elements in NGSolve. Tech.Rep. ASC-30/2014, Institute for Analysis and Scientific Computing - TU Wien(2014). URL (2), 969–1006 (2018)[237] Schoeder, S., Sticko, S., Kreiss, G., Kronbichler, M.: High-order cut discontinuousGalerkin methods with local time stepping for acoustics. International Journal forNumerical Methods in Engineering (13), 2979–3003 (2020)[238] Sch¨utz, J., Aizinger, V.: A hierarchical scale separation approach for the hybridizeddiscontinuous Galerkin method. Journal of Computational and Applied Mathematics , 500 – 509 (2017)[239] Sevilla, R.: HDG-NEFEM for two dimensional linear elasticity. Computers & Struc-tures , 69–80 (2019)[240] Sevilla, R., Borchini, L., Giacomini, M., Huerta, A.: Hybridisable discontinuousGalerkin solution of geometrically parametrised Stokes flows. Computer Methods inApplied Mechanics and Engineering , 113,397 (2020)[241] Sevilla, R., Fern´andez-M´endez, S., Huerta, A.: NURBS-enhanced finite elementmethod (NEFEM). International Journal for Numerical Methods in Engineering (1), 56–83 (2008)[242] Sevilla, R., Fern´andez-M´endez, S., Huerta, A.: 3D NURBS-enhanced finite elementmethod (NEFEM). International Journal for Numerical Methods in Engineering (2), 103–125 (2011)[243] Sevilla, R., Giacomini, M., Huerta, A.: A face-centred finite volume method forsecond-order elliptic problems. International Journal for Numerical Methods in En-gineering (8), 986–1014 (2018)[244] Sevilla, R., Giacomini, M., Huerta, A.: A locking-free face-centred finite volume(FCFV) method for linear elastostatics. Computers & Structures , 43–57 (2019)[245] Sevilla, R., Giacomini, M., Karkoulias, A., Huerta, A.: A superconvergent hybridis-able discontinuous Galerkin method for linear elasticity. International Journal forNumerical Methods in Engineering (2), 91–116 (2018)[246] Sevilla, R., Huerta, A.: Tutorial on Hybridizable Discontinuous Galerkin (HDG) forsecond-order elliptic problems. In: J. Schr¨oder, P. Wriggers (eds.) Advanced FiniteElement Technologies, CISM International Centre for Mechanical Sciences , vol. 566,pp. 105–129. Springer International Publishing (2016)[247] Sevilla, R., Huerta, A.: HDG-NEFEM with degree adaptivity for Stokes flows. Jour-nal of Scientific Computing (3), 1953–1980 (2018)88248] Sheldon, J.P., Miller, S.T., Pitt, J.S.: A hybridizable discontinuous Galerkin methodfor modeling fluid-structure interaction. Journal of Computational Physics , 91– 114 (2016)[249] Shen, J., Singler, J.R., Zhang, Y.: HDG-POD reduced order model of the heatequation. Journal of Computational and Applied Mathematics , 663–679 (2019)[250] Solano, M., Vargas, F.: A high order HDG method for Stokes flow in curved domains.Journal of Scientific Computing (3), 1505–1533 (2019)[251] Soon, S.C., Cockburn, B., Stolarski, H.K.: A hybridizable discontinuous Galerkinmethod for linear elasticity. International Journal for Numerical Methods in Engi-neering (8), 1058–1092 (2009)[252] Stanglmeier, M., Nguyen, N.C., Peraire, J., Cockburn, B.: An explicit hybridizablediscontinuous Galerkin method for the acoustic wave equation. Computer Methodsin Applied Mechanics and Engineering , 748–769 (2016)[253] Stenberg, R.: Some new families of finite elements for the Stokes equations. Nu-merische Mathematik (8), 827–838 (1990)[254] Su, W., Wang, P., Zhang, Y., Wu, L.: A high-order hybridizable discontinuousGalerkin method with fast convergence to steady-state solutions of the gas kineticequation. Journal of Computational Physics , 973–991 (2019)[255] Terrana, S., Nguyen, N.C., Bonet, J., Peraire, J.: A hybridizable discontinuousGalerkin method for both thin and 3D nonlinear elastic structures. Computer Meth-ods in Applied Mechanics and Engineering , 561 – 585 (2019)[256] Terrana, S., Vilotte, J., Guillot, L.: A spectral hybridizable discontinuous Galerkinmethod for elastic-acoustic wave propagation. Geophysical Journal International (1), 574–602 (2017)[257] Vidal-Codina, F., Mart´ın-Moreno, L., Cirac`ı, C., Yoo, D., Nguyen, N.C., Oh, S.H.,Peraire, J.: Terahertz and infrared nonlocality and field saturation in extreme-scalenanoslits. Optics Express (6), 8701–8715 (2020)[258] Vidal-Codina, F., Nguyen, N., Oh, S.H., Peraire, J.: A hybridizable discontinuousGalerkin method for computing nonlocal electromagnetic effects in three-dimensionalmetallic nanostructures. Journal of Computational Physics , 548 – 565 (2018)[259] Vidal-Codina, F., Nguyen, N., Peraire, J.: Computing parametrized solutions forplasmonic nanogap structures. Journal of Computational Physics , 89 – 106(2018) 89260] Vidal-Codina, F., Nguyen, N.C., Giles, M.B., Peraire, J.: A model and variancereduction method for computing statistical outputs of stochastic elliptic partial dif-ferential equations. Journal of Computational Physics , 700–720 (2015)[261] Vidal-Codina, F., Nguyen, N.C., Giles, M.B., Peraire, J.: An empirical interpo-lation and model-variance reduction method for computing statistical outputs ofparametrized stochastic partial differential equations. SIAM-ASA Journal on Uncer-tainty Quantification (1), 244–265 (2016)[262] Vieira, L.M., Giacomini, M., Sevilla, R., Huerta, A.: A second-order face-centredfinite volume method for elliptic problems. Computer Methods in Applied Mechanicsand Engineering , 112,655 (2020)[263] Vila-P´erez, J., Giacomini, M., Sevilla, R., Huerta, A.: Hybridisable discontinuousGalerkin formulation of compressible flows. Tech. rep., arXiv (2020). Submitted.[264] Wang, C.Y.: Exact solutions of the steady-state Navier-Stokes equations. AnnualReview of Fluid Mechanics (1), 159–177 (1991)[265] Wildey, T., Muralikrishnan, S., Bui-Thanh, T.: Unified geometric multigrid algo-rithm for hybridized high-order finite element methods. SIAM Journal on ScientificComputing (5), S172–S195 (2019)[266] Williams, D.M.: An entropy stable, hybridizable discontinuous Galerkin method forthe compressible Navier-Stokes equations. Mathematics of Computation (309),95–121 (2018)[267] Woopen, M., Balan, A., May, G., Sch¨utz, J.: A comparison of hybridized and stan-dard DG methods for target-based hp -adaptive simulation of compressible flow. Com-puters & Fluids , 3 – 16 (2014)[268] Woopen, M., May, G., Sch¨utz, J.: Adjoint-based error estimation and mesh adapta-tion for hybridized discontinuous Galerkin methods. International Journal for Nu-merical Methods in Fluids (11), 811–834 (2014)[269] Xie, Z.Q., Sevilla, R., Hassan, O., Morgan, K.: The generation of arbitrary ordercurved meshes for 3D finite element analysis. Computational Mechanics , 361–374(2013)[270] Yang, Y., Shi, K., Fu, S.: Multiscale hybridizable discontinuous Galerkin methodfor flow simulations in highly heterogeneous media. Journal of Scientific Computing (3), 1712–1731 (2019)[271] Yoo, D., Vidal-Codina, F., Cirac`ı, C., Nguyen, N.C., Smith, D.R., Peraire, J., Oh,S.H.: Modeling and observation of mid-infrared nonlocality in effective epsilon-near-zero ultranarrow coaxial apertures. Nature Communications10