HDG-NEFEM with Degree Adaptivity for Stokes Flows

The NURBS-enhanced finite element method (NEFEM) combined with a hybridisable discontinuous Galerkin (HDG) approach is presented for the first time. The proposed technique completely eliminates the uncertainty induced by a polynomial approximation of curved boundaries that is common within an isoparametric approach and, compared to other DG methods, provides a significant reduction in number of degrees of freedom. In addition, by exploiting the ability of HDG to compute a postprocessed solution and by using a local a priori error estimate valid for elliptic problems, an inexpensive, reliable and computable error estimator is devised. The proposed methodology is used to solve Stokes flow problems using automatic degree adaptation. Particular attention is paid to the importance of an accurate boundary representation when changing the degree of approximation in curved elements. Several strategies are compared and the superiority and reliability of HDG-NEFEM with degree adaptation is illustrated.


Introduction
Early work on mesh and degree adaptivity schemes for the finite element method 30,44,65 already showed the advantages of adaptive schemes to achieve a required accuracy in an economic manner.The use of mesh adaptive methods is substantially more extended due to the popularity of low-order methods in the computational mechanics community.This is largely due, as discussed later, to the fact that mesh adaptation is easier to implement, compared to degree adaptivity, in standard finite element codes.But, with recent needs on high fidelity simulations for fluids and wave propagation phenomena, 14,25,63 the interest in degree adaptive (or the combination of mesh and degree adaptivity) processes has increased. 4,23,24,31 Onf the main reasons for the increasing popularity of degree adaptive schemes in the last years is the rise of discontinuous Galerkin (DG) methods as a viable alternative for convection dominated flow and wave propagation problems.13,15,21,26,43,58 In a standard continuous Galerkin framework, the implementation of variable degree of approximation is cumbersome, whereas its application in a DG context is straightforward due to the weak imposition of the continuity of the solution by means of numerical fluxes.Despite traditional DG methods have not been able to consistently prove its superiority against low-order techniques traditionally employed in industry (e.g.finite volume methods), the recently proposed hybridisable DG (HDG) 11 has shown its superiority compared to traditional DG methods.9,27,33 The ability to substantially reduce the number of degrees of freedom combined with the possibility to obtain a post-processed solution that converges at a faster rate to the exact solution are the two main properties of HDG methods behind its superiority compared to other DG methods.10,12,38,56 Moreover, this is achieved while preserving the well-known advantages of DG for stabilizing convection and circumventing the so-called Ladyzhenskaya-Babuška-Brezzi (LBB) condition in the incompressible limit.
A key aspect in any adaptive scheme is the ability to devise cheap and reliable error measures for a given numerical solution in order to decide the regions where a more accurate solution is required. 1Error indicators and error estimators are typically employed to asses the error of a simulation with an adaptive framework. 29Error indicators are computationally inexpensive but they are problem dependent whereas error estimators are considerably more expensive but more general. 19,46,47 Acheap, general and reliable error estimator was proposed in 23,24 by exploiting the ability of the HDG method to construct a post-processed solution, more accurate than the HDG solution.
One of the aspects that is normally ignored when devising degree adaptive schemes is the geometric representation of domains with curved boundaries.Despite it is now well known that a poor representation of the geometry can have an important effect on the results of a finite element simulation, 2,7,54,60 the most extended practice consists on maintaining the shape of the elements during the degree adaptive process. 23,24,31 I the majority of cases, a polynomial representation of the boundary is selected whereas the polynomial degree of the functional approximation changes at each iteration of the degree adaptive scheme.This work analyses and discusses three approaches to perform a degree adaptive process in domains with curved boundaries.The first one corresponds to the approach typically employed in practice, consisting of fixing the shape of the curved elements and changing the degree of the functional approximation as dictated by the degree adaptivity procedure.The second approach proposed in this work is to employ the so-called NURBS-enhanced finite element method (NEFEM) that enables to exactly represent the geometry of the computational domain, given by a CAD model, irrespectively of the degree of the polynomials used to approximate the solution.The third approach, despite not considered useful from a practical point of view, consists of changing the geometry representation of the computational domain to represent with the same degree of polynomials both the geometry and the solution at each iteration of the degree adaptive process.This approach is not considered of interest from a practical point of view because it requires communication with the CAD model at each iteration and re-generation of nodal distributions for curved elements.
The second approach proposed here considers, for the first time, the combination of the so-called NURBS-enhanced finite element method (NEFEM) and the HDG rationale.The resulting method combines all the advantages of both methods, that is the efficiency of HDG and the ability of NEFEM to decouple the functional approximation from the geometric representation, usually tied in traditional isoparametric implementations.
A number of numerical examples is considered in order to compare the different degree adaptivity approaches.Furthermore, this work presents a simple idea to verify computational methods that are able to use different degrees of approximation for the solution in different elements.The idea is based on an existing local a priori error estimator developed in 18 for elliptic problems.
The remainder of the paper is organised as follows.Sections 2 briefly presents the model problem considered (i.e.Stokes flows) and the HDG formulation.The spatial discretisation of the HDG weak formulation is presented in Section 3 for both isoparametric and NEFEM, with particular emphasis on the differences between both formulations.The details about the proposed error estimator and degree adaptivity process proposed are presented in Section 4, including a discussion of the three approaches considered to perform a degree adaptive process.In Section 5 a simple technique to verify the implementation of a solver with variable degree of approximation is presented and used to test the implementation of the HDG code for Stokes flows with isoparametric and NEFEM.Section 6 presents a comparison of the different degree adaptive approaches and a number of numerical examples are used in Section to show the potential of the proposed approach.The so-called velocity-pressure formulation is obtained by invoking the Stoke's law and results in where u is the velocity vector, ν is the kinematic viscosity, p denotes the dynamic pressure, s is a body force, u D is the imposed velocity on the Dirichlet boundary Γ D , n is the outward unit normal vector to ∂Ω and t is the pseudo-traction vector imposed on the Neumann boundary Γ N .The disjoint boundaries Γ D and Γ N satisfy In what follow, •, • D denotes the L 2 scalar product in a generic subdomain D ⊂ Ω, that is for scalars, vectors and second order tensors respectively.Analogously, •, • S denotes the L 2 scalar product in any domain S ⊂ Γ ∪ ∂Ω.
The free divergence condition in Equation ( 1) induces the compatibility condition It is worth noting that, if only Dirichlet boundary conditions are considered (i.e.Γ N = ∅), an additional constraint on the pressure must be imposed to avoid its indeterminacy.It is common 10,12,38 to impose the mean pressure on the element boundary, namely p, 1 ∂Ω = 0. (3)

HDG weak formulation
The domain Ω is assumed partitioned in n el disjoint subdomains Ω e with boundaries ∂Ω e , which define an internal interface Γ The corresponding strong form of the Stokes system given in Equation ( 1) can be written in mixed form and in the broken computational domain as for e = 1, . . ., n el , where I is the identity tensor of dimension n sd , L = −∇u is a new variable (the second order velocity gradient tensor) which is introduced after splitting the second order momentum conservation equation in two first order equations and û is an independent variable representing the trace of the solution in ∂Ω e \ Γ D .
The free divergence condition in Equation ( 5) induces the compatibility condition The last two equations in (5) impose the continuity of velocity and continuity of the normal component of the pseudo-stress across the interior faces respectively, where the jump • operator has been introduced following the definition in, 37 such that, along each portion of the interface Γ it sums the values from the element on the left and right of say, Ω e and Ω l , namely = e + l .
[42] First, an element-by-element problem is defined with (L, u, p) as unknowns.This is the so-called local problem and is given by for e = 1, . . ., n el , where the last equation in (7) has been introduced to remove the indeterminacy of the pressure and ρ e denotes the mean pressure on the boundary of element Ω e .The local problem is used to obtain the solution in each element, L, u and p, for e = 1, . . ., n el , in terms of û and ρ along the interface Γ ∪ Γ N .
Second, a global problem is defined to determine the traces of the velocity and the mean pressure, denoted by û and ρ, on the element boundaries.This is given by where the first equation is automatically satisfied due to the unique definition of the hybrid variable û on each face of the mesh skeleton and the condition u = û on Γ, as imposed in the local problem.
The weak formulation for each element equivalent to (7) is as follows: for e = 1, . . ., n el , given u D on Γ D and û on Γ ∪ Γ N , find (L, u, p) , where the numerical trace of the normal flux is defined as elsewhere, (10)   with τ being a stabilization parameter, whose selection has an important effect on the stability, accuracy and convergence properties of the resulting HDG method.1][42] Introducing the definition of the numerical trace of the normal flux in Equation ( 9) leads to the weak form of the local problem: for e = 1, . . ., n el , find (L, u, p) For the global problem, the weak formulation equivalent to (8) Introducing the definition of the numerical trace of the normal flux in Equation ( 12) leads to the weak form of the global problem:

Spatial discretisation
This section presents the discretisation of the HDG weak forms derived in the previous section.Both the standard isoparametric and the so-called NEFEM formulations are presented.Special attention is paid to the differences between both formulations as this represents the first time NEFEM is considered in an HDG framework.

Isoparametric elements
Standard isoparametric formulations map each element Ω e and face Γ e in the physical domain into a reference element, Ω, and a reference face, Γ, where polynomial functional approximations characterize the discrete finite dimensional spaces.Namely, P k ( Ω) and P k( Γ) are the spaces of polynomial functions of degree at most k ≥ 1 and k ≥ 1 in the reference element and the reference face respectively.Finally, the approximations for each variable are defined as where u j , p j , L j and ûj are nodal values, N j are polynomial shape functions of order k in the reference element, n en is the number of nodes per element, Nj are polynomial shape functions of order k in the reference face and n fn is the corresponding number of nodes per face.Note that equal interpolation is used for all element variables (i.e.velocity, pressure and gradient of velocity).Recall that HDG allows for equal interpolation because of the numerical fluxes and the stabilization parameter τ .They ensure solvability and stability, see, 8 without the need of an enriched space for the gradient variable, or a reduced space for the trace variable.
An isoparametric mapping is used to link the reference element Ω and the computational element where x j are the nodal coordinates of the computational element Ω h e .It is worth noting that in general, when the physical element Ω e is curved, the isoparametric mapping is non-linear and the approximation defined in the reference element do not induce a polynomial interpolation in the physical space.In addition, the computational element Ω h e is just an approximation of Ω e , see 53 for a detailed discussion.Similarly, an isoparametric mapping is used to link the reference face Γ and the com- where x j denote the face nodal coordinates.
Using the mappings in Equations ( 18) and ( 19), the integrals appearing in the weak form of the local problems are transformed to the reference element and reference face/edge respectively.Then, the nodal interpolations given by Equations ( 14) to ( 17) are introduced, leading to a system of equations for each element with the following structure where ζ is the Lagrange multiplier corresponding to the constraint of Equation (11d).
Analogously, using the isoparametric mappings, the nodal interpolations of the corresponding variables and introducing the expression of L, u and p from Equation (20) in the global problem of Equation ( 13), a global system of equations is obtained where the vector of unknowns Û contains the nodal values of the trace of the velocity on the elementa faces and the mean pressure within each element.

NEFEM elements
In NEFEM, the boundary of the computational domain ∂Ω is exactly represented by NURBS.In what follows, in order to simplify the presentation and without loss of generality, the NURBS are restricted to two dimensional problems, see 52 for a detailed description of the three dimensional case.An edge is given by Γ e := C([λ e a , λ e b ]), where C is the NURBS boundary parametrisation and λ a and λ b are the parametric coordinates (in the parametric space of the NURBS) of the end points of Γ e .
The discrete approximations are defined now as: where u j , p j , L j and ûj are nodal values, N j are polynomial shape functions of order k in the physical element, n en is the number of nodes per element, Nj are polynomial shape functions of order k in [λ e a , λ e b ] and n fn is the corresponding number of nodes per face.The main differences of NEFEM with respect to the isoparametric formulation are: • The exact description of the computational domain is considered by means of its NURBS boundary representation.
• The approximation of the elemental variables directly in the physical space, with Cartesian coordinates.
• The approximation of the trace of the velocity is defined in the parametric space of the NURBS.It is worth noting that other options could be considered such as defining the approximation directly in the physical space.The main advantage of defining the approximation in the parametric space of the NURBS is that the number of unknowns remains the same as in the isoparametric formulation.In contrast, if the approximation of this variable is selected in the physical space it would require further degrees of freedom. 53,54 addition, from the computational point of view, NEFEM uses specifically designed numerical quadratures that provide a more efficient alternative to standard quadratures defined in a reference triangle. 51,52 or instance, in two dimensions, the following mapping is introduced between a reference rectangle and the physical element where R = [λ e a , λ e b ] × [0, 1] and x I is the internal vertex of Ω e .Using the mapping in Equation ( 26) and the NURBS boundary representation given by C, the integrals appearing in the weak form of the local problems are transformed to the reference rectangle and the parametric space of the NURBS respectively.Then, the nodal interpolations given by Equations ( 22), ( 23), ( 24) and ( 25) are introduced, leading to a system of equations similar to Equation (20).Analogously, the global problem with NEFEM leads to a global system of equations similar to Equation (21).

Error estimation and adaptivity
In HDG, the possibility to obtain a postprocessed solution 11 that converges at a higher rate (i.e.k + 2) than the HDG solution, not only provides a higher accurate solution to the problem at hand but it can also be used to build an inexpensive, reliable and computable error estimator. 23,24 n this section, particular attention is paid to the fact that, when the degree of approximation is changed in a curved element, a choice must be made regarding the geometric definition of the element.
An element by element measure of the error is defined by employing the HDG solution and the postprocessed solution as proposed in 24 where the normalisation becomes critical when meshes with different element sizes are considered. 18r elliptic problems, and by using that the influence of pollution errors becomes negligible if the mesh is sufficiently refined in the area where the pollution error is generated, 28 the following a priori error estimate was derived in 18 for e = 1, . . ., n el , where C is a constant, h e the characteristic element size of Ω e and k e the degree of approximation used in Ω e .It is worth noting that the error estimate of Equations ( 28) was initially derived for the standard finite element method but its extension to the HDG method is straightforward.
By applying a standard Richardson extrapolation, it is possible to predict the required change in the degree of approximation in order to ensure that the error in each element is lower than a desired accuracy , namely for e = 1, . . ., n el , where • denotes the ceiling function.
The adaptive procedure consist on solving the Stokes problem using the HDG formulation as described in Section 3 and estimating the required degree of approximation in each element according to Equation (29).The process is repeated until convergence is achieved, meaning that the error in each element ε e is lower than the desired error .

Geometry update
The technique described to drive a degree adaptive process only focuses on the degree of approximation used for the functional approximation, but in the presence of curved boundaries it is known that high-order approximations of both the solution and the geometry are required to exploit the full potential of a high-order method. 2,17,34,55 Ths aspect is usually ignored as degree adaptive procedures are applied to problems involving polygonal boundaries, see for instance. 16,22,45,61 Hee three options are discussed and assessed and compared later using numerical examples.
The first, and the one typically considered in a degree adaptive process, technique consists of defining a polynomial representation of the curved boundaries that is maintained during the adaptive process, irrespectively of the degree of approximation used for the solution. 23,24,35 Tis option is attractive because when the degree of approximation is changed in an element, there is no need to communicate with a CAD library to re-generate the nodal distributions in curved elements at each iteration of the degree adaptive process.The strategy is illustrated in Figure 1.The first row of plots show triangular elements where the geometric approximation is linear (q = 1) and the polynomial degree of approximation  of the solution increases from k = 1 to k = 3.The second row shows a similar situation where the boundary of the computational domain is described using quadratic polynomials (q = 2) and the degree of approximation of the solution is increased from k = 2 to k = 4.The boundary of the computational domain is denoted by ∂Ω h whereas the exact boundary is denoted by ∂Ω.

∂Ω ∂Ω ∂Ω
The second alternative, proposed in this work, consists of using NEFEM, where the exact boundary representation of the computational domain is considered irrespectively of the degree of approximation considered for the solution.As NEFEM encapsulates the necessary information to define the approximation and perform the numerical integration in curved elements in contact with a NURBS boundary, communication with a CAD library is avoided.The strategy is illustrated in Figure 2, showing a NEFEM element where the exact boundary representation is considered and a degree of approximation for the solution being increased from k = 1 to k = 3.
A third alternative, not considered in practice, consists of communicating with the CAD model after each iteration of the degree adaptive process in order to re-generate the Figure 3: Illustration of a degree adaptation in an element where the same degree of approximation is used for both the solution and the geometry.
nodal distribution of curved elements by placing the nodes over the true boundary.The strategy is illustrated in Figure 3, showing a triangular element where both the degree of the functional approximation and the degree of the polynomials used to approximate the solution are updated at each iteration.This strategy has not been considered in practical applications due to the cost associated to communicating with the CAD model at each iteration.
Remark 1.It is important to note that the first strategy, where the geometry remains unchanged, does not guarantee the convergence of the numerical solution to the physical solution in domains with curved boundaries because the distance between the computational domain and the physical domain does not converge to zero with as the degree of approximation is increased, see 6,49 for more details.For the second strategy, proposed here, convergence to the physical solution is guaranteed because no geometrical error is introduced. 54Finally, for the third approach, convergence is also guaranteed if the distance between the computational boundary and the physical boundary tends to zero as the order of the approximation is increased and the derivatives of the isoparametric mapping up to order k + 1 are bounded by h s , for s = 2, . . ., k + 1, 6,49 where h denotes the characteristic element size.It is worth noting that a specifically designed nodal distribution for curved elements is required in the third approach to guarantee that the second hypothesis is fulfilled. 6Validation of the HDG formulation with variable degree of approximation The first example provides a novel and simple technique to fully validate a solver that employs variable different degree of approximation in different elements for the solution of elliptic problems.The idea consists of utilising the local a priori error estimate of Equation ( 28) that states how the error, measured in an element, decreases when the mesh is refined.
To illustrate the proposed technique and validate the HDG isoparametric and NEFEM implementations with variable degree of approximation, the Stokes equations are solved in a circle of radius 0.5 centred at (0.5, 0.5) with Dirichlet boundary conditions.The viscosity is considered as ν = 1 and source and boundary conditions are taken such that the analytical solution is given by Six triangular meshes of the domain are generated using nested refinement.The first three meshes are shown in Figure 4, where the colour of each element represents the degree of approximation used, ranging from k = 1 to k = 6.In each mesh, there is one element per degree of approximation highlighted with a thicker line and darker colour, representing the region where the error is measured to test the local a priori error estimate.It is worth noting that the meshes shown in Figure 4 are NEFEM meshes, as the exact boundary representation is always employed, even for k = 1, whereas for the computations both NEFEM and isoparametric meshes are employed.
The results of the h-convergence study are presented in Figures 5 and 6. Figure 5 shows the error of the solution u and the postprocessed solution u in the L 2 (Ω e ) norm for isoparametric and NEFEM elements for k = 1, . . ., 5. The optimal rate of convergence is obtained in all cases for both the solution u (rate k + 2) and the postprocessed solution u (rate k + 3).
In Figure 6 a similar analysis is conducted, but now the error is measured for the dual variable L and the pressure p, also in the L 2 (Ω e ) norm and for isoparametric and NEFEM elements.Again, the optimal rate of convergence is obtained in all cases for both L and p (rate k + 2).(Ω) )

Comparison of degree adaptivity strategies
The same model problem employed in the previous example is utilised to compare the strategies described in Section 4.1 to update the geometry during a degree adaptive process.The computational domain selected, shown in Figure 7, features an oscillatory boundary and represents a common problem encountered in biological transport applications, see for instance. 50More precisely, the curved part of the boundary is given by the curve f (x) = (1 + cos(5πx))/10.
Dirichlet boundary conditions are imposed on the polygonal part of the boundary whereas a Neumann boundary condition, corresponding to the exact traction derived from the solution in Equation ( 30) is imposed on the oscillatory part of the boundary.

No geometric update
First, the degree adaptive process with no communication with the CAD model is studied, as illustrated in Figure 1.The process starts with a degree of approximation k = 1 in all elements.At each iteration the degree of the functional approximation is adapted according to the strategy presented in Section 4 whereas a linear approximation of the geometry is kept irrespectively of the degree of the functional approximation.Figure 8 shows the original mesh, the estimated error and the exact error, computed by using the known analytical solution.The L 2 norm of the error is represented as a constant value in each element, showing a good agreement between the estimated and the exact error.: Sixth iteration of the degree adaptivity procedure with no geometric update using HDG isoparametric elements and q = 1.
After six iterations of the adaptivity process, the degree of approximation is adapted in each element as shown in Figure 9 (a) but a linear geometric approximation of the curved boundary is still considered.The estimated error in each element, shown in Figure 9 (b), is below the desired error which is 0.5 × 10 −2 in this example but the computation of the exact error, shown in Figure 9 (c), reveals a significant disparity when compared to the estimated error.
To better analyse the results, Figure 10 shows the evolution of the maximum estimated error in each element and the maximum exact in each element for different geometric approximations of the curved boundary.Figure 10 (a) corresponds to the case illustrated in Figures 8 and 9, where a linear approximation of the geometry is considered (q = 1).It can be clearly observed that, as the degree adaptive process evolves, the difference between estimated and exact error becomes more and more sizeable.In the sixth iteration, the adaptive process shows convergence, with an estimated error less than the desired error 0.5 × 10 −2 but this is two orders of magnitude lower than the exact error, clearly indicating that the error estimator is not reliable because the estimator assumes that the geometry is exactly represented.
As a linear approximation of the geometry is well known to be not suitable when high order functional approximations are considered, 2, 7, 54, 60 the same experiment is repeated by using a more accurate boundary representation.The plots in Figure 10 (b), (c) and (d) show the evolution of the maximum estimated error in each element and the maximum exact in each element for quadratic, cubic and quartic approximation of the geometry.In all cases it is clearly observed that the error estimator is not reliable because the adaptive process converges but the exact error is more than one order of magnitude higher than the desired error.
The degree of approximation, estimated error and exact error obtained in the last iteration of the adaptive process for q = 4 is represented in Figure 11.The results show that even with a more accurate geometric approximation, the exact error in the elements Figure 10: Evolution of the estimated and exact errors during a degree adaptivity process for different degrees of the polynomials used to approximate the geometry (q).close to the curved boundary is much higher than the estimated error.There is clear evidence that, if no communication with a CAD model is undertaken during the degree adaptive process, the original mesh must be pre-adapted manually in order to ensure that the geometric error is small enough in order to ensure that the error estimator is reliable, clearly compromising the robustness of the whole adaptivity process.

NEFEM HDG
The strategy proposed in this work consists of utilising NEFEM, where the geometry is always given by its CAD boundary representation, irrespective of the degree of the functional approximation.In the context of a degree adaptive process, this means that no  The process starts with a degree of approximation k = 1 in all elements.At each iteration the degree of the functional approximation is adapted according to the strategy presented in Section 4 and a new nodal distribution is generated for each curved element.
Figure 12 shows the original mesh, the estimated error and the exact error, computed by using the known analytical solution.It is worth emphasising that, even when the degree of the functional approximation used is linear (k = 1) the exact boundary representation is considered, as shown in Figure 12 (a).The results show a very similar distribution for the estimated and exact errors.
In this case, after only three iterations of the degree adaptive process convergence is achieved.The degree of approximation used in each element, the estimated and the exact     errors in each element are represented in Figure 13.
Remark 2. As discussed in Section 4.1, an alternative, not employed in practice due to the high cost induced by the re-generation of the mesh at each iteration of the adaptive process, consists of changing both the degree of approximation for the solution and for the geometry during the adaptivity process, as illustrated in Figure 3.
To illustrate the superiority of NEFEM, Figure 14 shows the evolution of the maximum estimated error in each element and the maximum exact error in each element for isoparametric and NEFEM approaches and for two magnitudes of the desired error.Figure 10 (a) corresponds to the case previously illustrated, where the desired error is 0.5×10 −2 , whereas Figure 10 (b) shows the same study but with a desired error of 0.5 × 10 −3 .In both cases the superiority of NEFEM is clear as the desired error is achieved by substantially reducing the computational cost (i.e. the number of iterations of the degree adaptive process and the number of degrees of freedom required to achieve the required accuracy).
In addition, it is worth emphasising that the isoparametric approach requires communication with the CAD model in each iteration to re-generated the high-order nodal distribution.These nodal distributions in curved elements must be specifically designed to ensure optimal convergence of the isoparametric approach, 3, 5 while for NEFEM the Cartesian approximation of the solution ensures that the accuracy of the approximation is much less sensitive to the quality of the nodal distribution.

Numerical Examples
This section presents four numerical examples to illustrate the potential of NEFEM when combined with HDG to perform a degree adaptive process.The examples involve geometries with curved boundaries and where coarse meshes are considered to show the robustness of the proposed methodology.In all the examples the high-order isoparametric and NEFEM meshes are generated using the techniques described in 48, 64 and 57 respectively.

Flow in a channel with randomly distributed ellipses
The first example, similar to a test case presented in, 36 considers the flow around a set of randomly distributed set of 25 ellipses in a channel.Dirichlet boundary conditions are considered in the whole domain corresponding to a parabolic velocity profile on the left (inflow) and right (outflow) boundaries and zero (no-slip) Dirichlet boundary condition on the top and bottom walls and on the boundary of the ellipses, as illustrated in Figure 15.
A coarse mesh with 2,443 triangular elements is first considered.As no analytical solution is available, a reference solution is computed in a much finer mesh with 28,150 elements and by employing a degree of approximation k = 4.This reference solution is used to measure the accuracy of the adaptive computations performed in much coarse   meshes.
An adaptive process is performed using a quadratic approximation of the curved boundaries and standard isoparametric elements with a desired error of 0.5 × 10 −3 .Figure 16 (a) shows the computational mesh and degree of approximation after five iterations.In the vicinity of the ellipses the majority of elements have a cubic degree of approximation where the elements in contact with the ellipses need a higher order of approximation to capture all the flow features.The highest order of approximation is k = 6, used, as expected, with the regions with higher curvature of the boundary.Figure 16 (b) and (c) show the estimated and the reference errors after the adaptive process converged.The discrepancy between the estimated and the reference errors is clearly observed.Despite the adaptive process converges, meaning that all elements have an estimated error below the desired error, a total of 408 elements have a reference error above the desired tolerance.
When the same computation is performed by considering a cubic approximation of the geometry (not reported for brevity), the adaptive process converges again in five iterations. he highest order of approximation used in a few elements is now k = 7, indicating that a different geometric representation leads to a different degree of approximation required to achieve convergence.In addition, the error estimator is again not reliable as there are 15 elements where the reference error is above the desired error.
It is apparent that an adaptive computation with isoparametric elements requires an initial pre-adaptation of the mesh and the degree of approximation used to approximate the geometry in order to obtain a reliable error estimator.Next, the finer mesh with 4,048 triangular elements depicted in Figure 17 (a) is considered.When a quadratic approximation of the geometry is considered, the adaptive process converges in four iterations but the results are still not satisfactory as there are 155 elements where the estimated error, shown in Figure 17 (b), is above the desired error, represented in Figure 17 (c).
Finally, if a cubic approximation of the geometry is employed the adaptive process converges in only three iterations with one element still showing an error above the desired tolerance.The results clearly indicate that in the presence of curved boundaries the level of pre-adaptation required negates all the advantages of an an automatic adaptive process as the initial mesh has to be designed in such a way so that the error in the first computation is near the desired error.
To show the potential of NEFEM in this scenario, an adaptive process is performed employing the coarse mesh with 2,443 triangular elements and starting with a degree of approximation k = 1.The adaptive process converges in four iterations.The final degree of approximation used in each element is shown in Figure 18 (a), with two elements having the maximum degree of approximation, k = 6, required to achieve the desired error.The estimated and reference errors, depicted in Figures 18 (b) and (c), respectively, shows a consistent behaviour that illustrates the reliability of the proposed strategy to estimate the error due to the use of the exact boundary representation.It is worth noting that in the majority of the elements surrounding the ellipses a degree of approximation k = 3 is enough to obtain the required error, illustrating why cubic isoparametric elements outperformed the use of quadratic elements in the previous computations.
The velocity field computed with NEFEM on the mesh shown in Figure 18 (a) is depicted in Figure 19.
To better analyse the effect of an accurate boundary representation in a degree adaptive procedure, Figure 20 shows the evolution of the estimated and exact errors as a function of the number of iterations of the adaptive process by using the coarse mesh with 2,443 elements.For a desired error of 0.5 × 10 −3 a quadratic approximation of the geometry prevents convergence of the exact error whereas better results are obtained with a cubic  representation of the geometry.It is worth noting that in both cases the exact error stagnates, indicating that the geometric error dominates over the interpolation error.Even when a cubic representation of the geometry is considered, the exact error is not decreasing after the third iteration.
To further illustrate the limitations of an isoparametric formulation during a degree adaptivity procedure, the same analysis is repeated with a lower desired error, namely 0.5 × 10 −4 .Figure 21 shows the evolution of the estimated and exact errors as a function of the number of iterations of the adaptive process by using the coarse mesh with 2,443 elements.The results show that a cubic approximation of the geometry is not enough    Figure 20: Evolution of the estimated and exact errors during the degree adaptivity process for q = 2 and q = 3 with isoparametric HDG in the coarse mesh with 2,443 elements and a desired error of 0.5 × 10 −3 .Figure 21: Evolution of the estimated and exact errors during the degree adaptivity process for q = 2 and q = 3 with isoparametric HDG in the coarse mesh with 2,443 elements and a desired error of 0.5 × 10 −4 .
because the difference between the estimated and exact error in the final computation with cubic elements is almost an order of magnitude.This suggests, once more that the initial mesh has to be pre-adapted to achieve a reliable degree adaptive process.
Finally, Figure 22 shows the results obtained with NEFEM using the coarse mesh with 2,443 elements and starting with a linear approximation of the solution k = 1.The robustness of the proposed approach is clearly illustrated as convergence of both the estimated  and exact errors is achieved in the coarse mesh even for a desired error of 0.5 × 10 −4 .
It is worth emphasising that with NEFEM the adaptive process provides a reliable error estimator even when the desired error is several orders of magnitude lower than the error of the computation in the first mesh.The results clearly indicate that no pre-adaptation of the mesh is required with NEFEM as the geometry is exactly represented irrespectively of the spatial discretisation.Therefore, the adaptive process is purely driven by the functional approximation and not by the geometric error as it happens with an isoparametric formulation.
Further numerical experiments, not reported here, indicate that NEFEM is also superior to an isoparametric approach where the mesh is re-generated at each iteration of the adaptive procedure by using the same degree of the approximation for both the geometry and the solution.In all cases, not only the time required by NEFEM is lower (due to the extra time required to communicate with the CAD model and the mesh generator) but also due to the fact that more iterations of the adaptive process are required with an isoparametric formulation.

Flow in a channel with wavy boundaries
The next example considers the flow in a channel with oscillatory boundaries.This problem, of interest to the micro and nano-fluidics community, is often considered to study the flow structure induced by the different phase of the oscillations of the top and bottom boundaries. 32,62 he two extreme cases are considered here, where the oscillations are exactly on phase and completely out of phase.Figure 23 (a) shows the coarse computational mesh employed with HDG-NEFEM and the final degree of approximation obtained after ten iterations of the degree adaptive process for the case where the boundary oscillations are on phase.Similarly, Figure 23 (b) shows the mesh and degree of approximation obtained after eight iterations of the degree adaptive process for the case where the boundary oscillations are completely out of phase.The velocity fields obtained for both cases are represented in Figure 24 showing the ability of the proposed approach to capture the different flow structure induced by the oscillatory boundaries.

Flow in a porous media
The next example, taken from, 59 considers the flow in the interstices of a porous media.The geometry consists of the surroundings of a large number of particles in the porous media.Figure 25 shows the mesh and degree of approximation after eight iterations of the degree adaptivity procedure.It is worth noting that a linear degree of approximation is used in many elements in contact with curved boundaries.This shows that the proposed adaptivity strategy is completely driven by the complexity of the solution and not by the complexity of the geometry.

Flow in a channel with thin obstacles
The last example shows a further benefit of using NEFEM by demonstrating its unique ability to obtain accurate solutions with ultra-coarse meshes even when geometric features, smaller than the element size, are present in the boundary representation of the computational domain.
The flow in a channel with a number of thin obstacles is considered.The thickness of the obstacles is approximately 0.08 whereas the minimum element size of the mesh that has been generated, using the technique proposed in, 57 is 0.32.The degree adaptive process is started, as in previous examples, with a linear approximation of the solution and convergence is achieved in four iterations.The final degree of approximation in each element is represented in Figure 26 (a). Figure 26 shows a detailed view of a region in the channel showing the elements near the end of some of the obstacles.This plot shows not only that the element size is independent on the geometric complexity but it also demonstrates the robustness of the proposed degree adaptive technique.The adaptivity process is clearly driven only by the complexity of the solution as a different degree of  The velocity field, obtained on the mesh shown in Figure 26, is represented in Figure 27.

Concluding remarks
A new degree adaptive methodology that combines the advantages of the HDG formulation and NEFEM has been presented.The proposed method results in a cheap and reliable error estimator due to the cheap computation of a post-processed solution provided by HDG and the ability to exactly represent curved boundaries irrespectively of the polynomial degree used for the functional approximation that is characteristic of NEFEM.
The proposed approach is compared against two alternative options to perform degree adaptivity.The first approach, broadly used in practice, consists of keeping the shape of the curved elements during the degree adaptivity process.It is found that this approach leads to an unreliable error estimator.The numerical examples show that even when the estimated error is below the required tolerance, the exact error can be orders of magnitude higher.The second approach, not used in practice, consists of changing the shape of the curved elements during the adaptive process.The main drawback is its high cost due to the need to constantly communicate with the CAD model and re-generate the nodal distributions for curved elements.
The proposed approach considers, for the first time, the implementation of the NEFEM rationale in an HDG framework.A number of numerical examples have been presented to compare the performance of the proposed methodology and to shows its superiority on a number of problems involving domains with curved boundaries.

Figure 1 :
Figure 1: Illustration of a degree adaptation in an element with linear (top) and quadratic (bottom) approximation of the geometry.

Figure 2 :
Figure 2: Illustration of a degree adaptation in a NEFEM element.

Figure 4 :
Figure 4: First three NEFEM meshes where the colour indicates the degree of approximation used in each element and the highlighted element represents the region where the error is measured for each degree of approximation.

Figure 5 :
Figure 5: Error of the solution u and the postprocessed solution u in the L 2 (Ω e ) norm for different degrees of approximation in each element.

Figure 6 :
Figure 6: Error of the dual variable L and the pressure p in the L 2 (Ω e ) norm for different degrees of approximation in each element.

Figure 7 :
Figure 7: Computational domain for the test problem used to compare the different geometry update options in a degree adaptive process.

Figure 8 :
Figure 8: First iteration of the degree adaptivity procedure with HDG isoparametric elements.

Figure 9
Figure9: Sixth iteration of the degree adaptivity procedure with no geometric update using HDG isoparametric elements and q = 1.

Figure 11 :
Figure 11: Sixth iteration of the degree adaptivity procedure with no geometric update using HDG isoparametric elements and q = 4.

Figure 12 :
Figure 12: First iteration of the degree adaptivity procedure with geometric update using HDG NEFEM elements.

Figure 13 :
Figure 13: Third iteration of the degree adaptivity procedure with geometric update using HDG NEFEM elements.

Figure 14 :
Figure 14: Estimated and exact errors for isoparametric and NEFEM.

Figure 15 :
Figure 15: Computational domain and boundary conditions for the problem involving a flow in a channel with randomly distributed ellipses.

Figure 16 :
Figure 16: Final degree distribution, estimated and reference errors for an adaptive computation with isoparametric HDG and quadratic approximation of the curved boundaries in a coarse mesh with 2,443 elements.

Figure 17 :
Figure 17: Final degree distribution, estimated and reference errors for an adaptive computation with isoparametric HDG and quadratic approximation of the curved boundaries in a finer mesh with 4,048 elements.

Figure 18 :
Figure 18: Final degree distribution, estimated and reference errors for an adaptive computation with HDG-NEFEM in a coarse mesh with 2,443 elements.

Figure 19 :
Figure 19: Magnitude of velocity and isolines computed with HDG-NEFEM on the mesh shown in Figure 18 (a) after four iterations of the degree adaptive process.

Figure 22 :
Figure 22: Evolution of the estimated and exact errors during the degree adaptivity process for NEFEM HDG in the coarse mesh with 2,443 elements.
(a) On phase oscillations (b) Out of phase oscillations

Figure 23 :
Figure 23: Mesh and degree of approximation of the converged degree adaptive procedure with HDG-NEFEM for the computation of the flow in a channel with oscillations of the top and bottom boundaries.
(a) On phase oscillations (b) Out of phase oscillations

Figure 24 :
Figure 24: Magnitude of velocity computed with HDG-NEFEM after convergence of the degree procedure on the meshes shown in Figure 23.

Figure 25 :
Figure 25: Mesh and degree of approximation of the converged degree adaptive procedure with HDG-NEFEM and velocity field.

Figure 26 :
Figure 26: Mesh and degree of approximation of the converged degree adaptive procedure with HDG-NEFEM for the computation of the flow in a channel with high curvature.

Figure 27 :
Figure 27: Velocity computed with HDG-NEFEM after convergence of the degree adaptive procedure in the of Figure 26.
Let us consider an open bounded domain Ω ∈ R nsd with boundary ∂Ω, where n sd the number of spatial dimensions.The strong form of the stationary Stokes problem is obtained by neglecting the transient and convective effects in the full incompressible Navier-Stokes equations.