HDGNEFEM with degree adaptivity for Stokes flows
Abstract
The NURBSenhanced 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 HDGNEFEM with degree adaptation is illustrated.
Keywords: Hybridisable Discontinuous Galerkin, NURBSenhanced finite element method, degree adaptivity, Stokes
1 Introduction
Early work on mesh and degree adaptivity schemes for the finite element method [65, 44, 30] 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 loworder 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 [63, 14, 25], the interest in degree adaptive (or the combination of mesh and degree adaptivity) processes has increased [4, 23, 31, 24].
One of 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, 43, 58, 26, 21, 15]. 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 loworder 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, 33, 27]. The ability to substantially reduce the number of degrees of freedom combined with the possibility to obtain a postprocessed 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 [38, 12, 10, 56]. Moreover, this is achieved while preserving the wellknown advantages of DG for stabilizing convection and circumventing the socalled LadyzhenskayaBabuškaBrezzi (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 [1]. Error indicators and error estimators are typically employed to asses the error of a simulation with an adaptive framework [29]. Error indicators are computationally inexpensive but they are problem dependent whereas error estimators are considerably more expensive but more general [46, 19, 47]. A cheap, general and reliable error estimator was proposed in [23, 24] by exploiting the ability of the HDG method to construct a postprocessed 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, 31, 24]. In 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 socalled NURBSenhanced 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 regeneration of nodal distributions for curved elements.
The second approach proposed here considers, for the first time, the combination of the socalled NURBSenhanced 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. Finally, Section 8 summarises the main conclusions of the work that has been presented.
2 Hybridisable discontinuous Galerkin for Stokes flow
2.1 Problem statement
Let us consider an open bounded domain with boundary , where 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 NavierStokes equations [20]. The socalled velocitypressure formulation is obtained by invoking the Stoke’s law and results in
(1) 
where is the velocity vector, is the kinematic viscosity, denotes the dynamic pressure, is a body force, is the imposed velocity on the Dirichlet boundary , is the outward unit normal vector to and is the pseudotraction vector imposed on the Neumann boundary . The disjoint boundaries and satisfy .
In what follow, denotes the scalar product in a generic subdomain , that is
for scalars, vectors and second order tensors respectively. Analogously, denotes the scalar product in any domain .
The free divergence condition in Equation (1) induces the compatibility condition
(2) 
2.2 HDG weak formulation
The domain is assumed partitioned in disjoint subdomains with boundaries , which define an internal interface
(4) 
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
(5) 
for , where is the identity tensor of dimension , 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 .
The free divergence condition in Equation (5) induces the compatibility condition
(6) 
The last two equations in (5) impose the continuity of velocity and continuity of the normal component of the pseudostress 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, and , namely
The HDG method solves problem (5) in two stages, see for instance [11, 8, 40, 41, 42]. First, an elementbyelement problem is defined with as unknowns. This is the socalled local problem and is given by
(7) 
for , where the last equation in (7) has been introduced to remove the indeterminacy of the pressure and denotes the mean pressure on the boundary of element . The local problem is used to obtain the solution in each element, , and , for , in terms of and along the interface .
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
(8) 
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 on , as imposed in the local problem.
The weak formulation for each element equivalent to (7) is as follows: for , given on and on , find that satisfies
(9a)  
(9b)  
(9c)  
(9d) 
for all , where the numerical trace of the normal flux is defined as
(10) 
with being a stabilization parameter, whose selection has an important effect on the stability, accuracy and convergence properties of the resulting HDG method. The influence of the stabilization parameter has been studied extensively by Cockburn and coworkers, see for instance [11, 8, 40, 41, 39, 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 , find such that
(11a)  
(11b)  
(11c)  
(11d) 
for all .
For the global problem, the weak formulation equivalent to (8) is: find and that satisfies
(12a)  
(12b) 
for all .
Introducing the definition of the numerical trace of the normal flux in Equation (12) leads to the weak form of the global problem: find and such that, for all ,
(13a)  
(13b) 
3 Spatial discretisation
This section presents the discretisation of the HDG weak forms derived in the previous section. Both the standard isoparametric and the socalled 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.
3.1 Isoparametric elements
Standard isoparametric formulations map each element and face in the physical domain into a reference element, , and a reference face, , where polynomial functional approximations characterize the discrete finite dimensional spaces. Namely, and are the spaces of polynomial functions of degree at most and in the reference element and the reference face respectively. Finally, the approximations for each variable are defined as
(14)  
(15)  
(16)  
(17) 
where , , and are nodal values, are polynomial shape functions of order in the reference element, is the number of nodes per element, are polynomial shape functions of order in the reference face and 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
(18)  
where are the nodal coordinates of the computational element .
It is worth noting that in general, when the physical element is curved, the isoparametric mapping is nonlinear and the approximation defined in the reference element do not induce a polynomial interpolation in the physical space. In addition, the computational element is just an approximation of , see [53] for a detailed discussion.
Similarly, an isoparametric mapping is used to link the reference face and the computational face
(19)  
where 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
(20) 
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 , and from Equation (20) in the global problem of Equation (13), a global system of equations is obtained
(21) 
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.
3.2 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 , where is the NURBS boundary parametrisation and and are the parametric coordinates (in the parametric space of the NURBS) of the end points of .
The discrete approximations are defined now as:
(22)  
(23)  
(24)  
(25) 
where , , and are nodal values, are polynomial shape functions of order in the physical element, is the number of nodes per element, are polynomial shape functions of order in and 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].
In 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]. For instance, in two dimensions, the following mapping is introduced between a reference rectangle and the physical element
(26)  
where and is the internal vertex of .
Using the mapping in Equation (26) and the NURBS boundary representation given by , 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).
4 Error estimation and adaptivity
In HDG, the possibility to obtain a postprocessed solution [11] that converges at a higher rate (i.e. ) 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 [24, 23]. In 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]
(27) 
where the normalisation becomes critical when meshes with different element sizes are considered [18].
For 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]
(28) 
for , where is a constant, the characteristic element size of and the degree of approximation used in . 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
(29) 
for , 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 is lower than the desired error .
4.1 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 highorder approximations of both the solution and the geometry are required to exploit the full potential of a highorder method [2, 17, 34, 55]. This aspect is usually ignored as degree adaptive procedures are applied to problems involving polygonal boundaries, see for instance [16, 45, 61, 22]. Here 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, 35, 24]. This 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 regenerate 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 () and the polynomial degree of approximation of the solution increases from to . The second row shows a similar situation where the boundary of the computational domain is described using quadratic polynomials () and the degree of approximation of the solution is increased from to . The boundary of the computational domain is denoted by 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 to .
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 regenerate the 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 [49, 6] for more details. For the second strategy, proposed here, convergence to the physical solution is guaranteed because no geometrical error is introduced [54]. Finally, 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 are bounded by , for [49, 6], where 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 [6].
5 Validation 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 with Dirichlet boundary conditions. The viscosity is considered as and source and boundary conditions are taken such that the analytical solution is given by
(30) 
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 to .
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 , whereas for the computations both NEFEM and isoparametric meshes are employed.
The results of the convergence study are presented in Figures 5 and 6. Figure 5 shows the error of the solution and the postprocessed solution in the norm for isoparametric and NEFEM elements for .
The optimal rate of convergence is obtained in all cases for both the solution (rate ) and the postprocessed solution (rate ).
In Figure 6 a similar analysis is conducted, but now the error is measured for the dual variable and the pressure , also in the norm and for isoparametric and NEFEM elements.
Again, the optimal rate of convergence is obtained in all cases for both and (rate ).
6 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 [50]. More precisely, the curved part of the boundary is given by the curve .
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.
6.1 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 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 norm of the error is represented as a constant value in each element, showing a good agreement between the estimated and the exact error.
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 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 (). 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 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 is represented in Figure 11.
The results show that even with a more accurate geometric approximation, the exact error in the elements 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 preadapted 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.
6.2 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 communication the CAD model is required as the exact boundary representation is already used by the NEFEM solver.
The process starts with a degree of approximation 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 () 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 regeneration 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 , whereas Figure 10 (b) shows the same study but with a desired error of . 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 regenerated the highorder 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.
7 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 highorder isoparametric and NEFEM meshes are generated using the techniques described in [64, 48] and [57] respectively.
7.1 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 (noslip) 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 . 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 . 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 , 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 , 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 preadaptation 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 preadaptation 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 . 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, , 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 is enough to obtain the required error, illustrating why cubic isoparametric elements outperformed the use of quadratic elements in the previous computations.
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 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 . 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 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 preadapted 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 .
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 . 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 preadaptation 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 regenerated 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.
7.2 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 nanofluidics 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]. The 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 HDGNEFEM 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.
7.3 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.
7.4 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 ultracoarse 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 approximation is employed in elements with almost identical geometric complexity due to the different complexity of the solution.
8 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 postprocessed 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 regenerate 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.
Acknowledgements
The authors gratefully acknowledge the financial support of the Ministerio de Economía y Competitividad (grant number DPI201451844C22R). The first author also gratefully acknowledges the financial support provided by the Sêr Cymru National Research Network for Advanced Engineering and Materials (grant number NRN045). The second author is also grateful for the financial support provided by the Generalitat de Catalunya (grant number 2014SGR1471).
References
 [1] W. Bangerth and R. Rannacher. Adaptive finite element methods for differential equations. Birkhäuser, 2013.
 [2] F. Bassi and S. Rebay. Highorder accurate discontinuous finite element solution of the 2D Euler equations. Journal of Computational Physics, 138(2):251–285, 1997.
 [3] C. Bernardi. Optimal finiteelement interpolation on curved domains. SIAM Journal on Numerical Analysis, 26(5):1212–1240, 1989.
 [4] A. Burbeau and P. Sagaut. A dynamic adaptive discontinuous Galerkin method for viscous flow with shocks. Computers & fluids, 34(4):401–417, 2005.
 [5] Q. Chen and I. Babuška. Approximate optimal points for polynomial interpolation of real functions in an interval and in a triangle. Computer Methods in Applied Mechanics and Engineering, 128(3–4):405–417, 1995.
 [6] P. G. Ciarlet. The finite element method for elliptic problems. SIAM, 2002.
 [7] F. Cirak, M. Ortiz, and P. Schroder. Subdivision surfaces: a new paradigm for thinshell finiteelement analysis. International Journal for Numerical Methods in Engineering, 47(12):2039–2072, 2000.
 [8] B. Cockburn, B. Dong, and J. Guzmán. A superconvergent LDGhybridizable Galerkin method for secondorder elliptic problems. Mathematics of Computation, 77(264):1887–1916, 2008.
 [9] B. Cockburn, B. Dong, J. Guzmán, M. Restelli, and R. Sacco. A hybridizable discontinuous Galerkin method for steadystate convectiondiffusionreaction problems. SIAM Journal of Scientific Computing, 31(5):3827–3846, 2009.
 [10] B. Cockburn and J. Gopalakrishnan. The derivation of hybridizable discontinuous Galerkin methods for Stokes flow. SIAM Journal on Numerical Analysis, 47(2):1092–1125, 2009.
 [11] B. Cockburn, J. Gopalakrishnan, and R. Lazarov. Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems. SIAM Journal on Numerical Analysis, 47(2):1319–1365, 2009.
 [12] B. Cockburn and K. Shi. Devising HDG methods for Stokes flow: an overview. Computers & Fluids, 98:221–229, 2014.
 [13] B. Cockburn and C.W. Shu. Runge–Kutta discontinuous Galerkin methods for convectiondominated problems. Journal of scientific computing, 16(3):173–261, 2001.
 [14] G. C. Cohen and Q. H. Liu. Higherorder numerical methods for transient wave equations. The Journal of the Acoustical Society of America, 114(1):21–21, 2003.
 [15] J. D. De Basabe, M. K. Sen, and M. F. Wheeler. The interior penalty discontinuous Galerkin method for elastic wave propagation: grid dispersion. Geophysical Journal International, 175(1):83–93, 2008.
 [16] L. Demkowicz. Fully automatic –adaptivity for Maxwells equations. Computer methods in applied mechanics and engineering, 194(2):605–624, 2005.
 [17] S. Dey, M. S. Shephard, and J. E. Flaherty. Geometry representation issues associated with version finite element computations. Computer Methods in Applied Mechanics and Engineering, 150(14):39–55, Dec. 1997.
 [18] P. Díez and A. Huerta. A unified approach to remeshing strategies for finite element adaptivity. Computer Methods in Applied Mechanics and Engineering, 176(14):215–229, 1999.
 [19] P. Díez, J. J. Ródenas, and O. C. Zienkiewicz. Equilibrated patch recovery error estimates: simple and accurate upper bounds of the error. International Journal for Numerical Methods in Engineering, 69(10):2075–2098, 2007.
 [20] J. Donea and A. Huerta. Finite element methods for flow problems. John Wiley & Sons, Chichester, 2003.
 [21] M. Dumbser and M. Käser. An arbitrary highorder discontinuous Galerkin method for elastic waves on unstructured meshes– II. The threedimensional isotropic case. Geophysical Journal International, 167(1):319–336, 2006.
 [22] M. Dumbser, M. Käser, and E. F. Toro. An arbitrary highorder discontinuous Galerkin method for elastic waves on unstructured meshes – V. Local time stepping and adaptivity. Geophysical Journal International, 171(2):695–717, 2007.
 [23] G. Giorgiani, S. FernándezMéndez, and A. Huerta. Hybridizable discontinuous Galerkin adaptivity for wave propagation problems. International Journal for Numerical Methods in Fluids, 72(12):1244–1262, 2013.
 [24] G. Giorgiani, S. FernándezMéndez, and A. Huerta. Hybridizable discontinuous Galerkin with degree adaptivity for the incompressible Navier–Stokes equations. Computers & Fluids, 98:196–208, 2014.
 [25] J. S. Hesthaven. Highorder accurate methods in timedomain computational electromagnetics: A review. Advances in Imaging and Electron Physics, 127:59–123, 2003.
 [26] J. S. Hesthaven and T. Warburton. Nodal highorder methods on unstructured grids: I. Timedomain solution of Maxwell’s equations. Journal of Computational Physics, 181(1):186–221, 2002.
 [27] A. Huerta, A. Angeloski, X. Roca, and J. Peraire. Efficiency of highorder elements for continuous and discontinuous Galerkin methods. International Journal for Numerical Methods in Engineering, 96(9):529–560, 2013.
 [28] A. Huerta and P. Díez. Error estimation including pollution assessment for nonlinear finite element analysis. Computer Methods in Applied Mechanics and Engineering, 181(1):21–41, 2000.
 [29] A. Huerta, A. RodríguezFerran, P. Díez, and J. Sarrate. Adaptive finite element strategies based on error assessment. International Journal for Numerical Methods in Engineering, 46(10):1803–1818, 1999.
 [30] C. Johnson and P. Hansbo. Adaptive finite element methods in computational mechanics. Computer methods in applied mechanics and engineering, 101(13):143–181, 1992.
 [31] P. Karban, F. Mach, and I. Doležel. Advanced adaptive algorithms in 2D finite element method of higher order of accuracy. COMPELThe international journal for computation and mathematics in electrical and electronic engineering, 32(3):834–849, 2013.
 [32] S. Khuri. Stokes flow in curved channels. Journal of computational and applied mathematics, 187(2):171–191, 2006.
 [33] R. Kirby, S. J. Sherwin, and B. Cockburn. To CG or to HDG: A comparative study. Journal of Scientific Computing, 51(1):183–212, 2011.
 [34] L. Krivodonova and M. Berger. Highorder accurate implementation of solid wall boundary conditions in curved geometries. Journal of Computational Physics, 211(2):492–512, Jan. 2006.
 [35] L. Li, S. Lanteri, and R. Perrussel. Numerical investigation of a high order hybridizable discontinuous Galerkin method for 2D timeharmonic Maxwell’s equations. COMPELThe international journal for computation and mathematics in electrical and electronic engineering, 32(3):1112–1138, 2013.
 [36] Y. Liu. A new fast multipole boundary element method for solving 2D Stokes flow problems based on a dual BIE formulation. Engineering Analysis with Boundary Elements, 32(2):139–151, 2008.
 [37] A. Montlaur, S. FernándezMéndez, and A. Huerta. Discontinuous Galerkin methods for the Stokes equations using divergencefree approximations. International Journal for Numerical Methods in Fluids, 57(9):1071–1092, 2008.
 [38] N. Nguyen, J. Peraire, and B. Cockburn. A hybridizable discontinuous Galerkin method for Stokes flow. Computer Methods in Applied Mechanics and Engineering, 199(9):582–597, 2010.
 [39] N. Nguyen, J. Peraire, and B. Cockburn. A hybridizable discontinuous Galerkin method for Stokes flow. Comput. Methods Appl. Mech. Eng., 199(912):582–597, 2010.
 [40] N. C. Nguyen, J. Peraire, and B. Cockburn. An implicit highorder hybridizable discontinuous Galerkin method for linear convectiondiffusion equations. Journal of Computational Physics, 228(9):3232–3254, 2009.
 [41] N. C. Nguyen, J. Peraire, and B. Cockburn. An implicit highorder hybridizable discontinuous Galerkin method for nonlinear convectiondiffusion equations. Journal of Computational Physics, 228(23):8841–8855, 2009.
 [42] N. C. Nguyen, J. Peraire, and B. Cockburn. An implicit highorder hybridizable discontinuous Galerkin method for the incompressible NavierStokes equations. Journal of Computational Physics, 230(4):1147–1170, 2011.
 [43] N. C. Nguyen, P.O. Persson, and J. Peraire. RANS solutions using high order discontinuous Galerkin methods. AIAA Paper, 914:2007, 2007.
 [44] J. T. Oden, L. Demkowicz, W. Rachowicz, and T. Westermann. Toward a universal adaptive finite element strategy, Part 2. A posteriori error estimation. Computer Methods in Applied Mechanics and Engineering, 77(12):113–180, 1989.
 [45] D. Pardo. Multigoaloriented adaptivity for –finite element methods. Procedia Computer Science, 1(1):1953–1961, 2010.
 [46] N. Parés, P. Díez, and A. Huerta. Subdomainbased fluxfree a posteriori error estimators. Computer Methods in Applied Mechanics and Engineering, 195(4):297–323, 2006.
 [47] N. Parés, P. Díez, and A. Huerta. Exact bounds for linear outputs of the advectiondiffusionreaction equation using fluxfree error estimates. SIAM journal on scientific computing, 31(4):3064–3089, 2009.
 [48] R. Poya, R. Sevilla, and A. J. Gil. A unified approach for a posteriori highorder curved mesh generation using solid mechanics. Computational Mechanics, 58(3):457–490, 2016.
 [49] P.A. Raviart, J.M. Thomas, P. G. Ciarlet, and J. L. Lions. Introduction à l’analyse numérique des équations aux dérivées partielles, volume 2. Dunod Paris, 1998.
 [50] M. Scholle. Creeping Couette flow over an undulated plate. Archive of Applied Mechanics, 73(1112):823–840, 2004.
 [51] R. Sevilla and S. FernándezMéndez. Numerical integration over 2D NURBSshaped domains with applications to NURBS–enhanced FEM. Finite Elements in Analysis and Design, 47(10):1209–1220, 2011.
 [52] R. Sevilla, S. FernándezMéndez, and A. Huerta. 3D NURBSenhanced finite element method (NEFEM). International Journal for Numerical Methods in Engineering, 88(2):103–125, 2011.
 [53] R. Sevilla, S. FernándezMéndez, and A. Huerta. Comparison of highorder curved finite elements. International Journal for Numerical Methods in Engineering, 87(8):719–734, 2011.
 [54] R. Sevilla, S. FernándezMéndez, and A. Huerta. NURBSenhanced finite element method (NEFEM). Archives of Computational Methods in Engineering, 18(4):441–484, 2011.
 [55] R. Sevilla, O. Hassan, and K. Morgan. An analysis of the performance of a highorder stabilised finite element method for simulating compressible flows. Computer Methods in Applied Mechanics and Engineering, 253:15–27, 2013.
 [56] R. Sevilla and A. Huerta. Tutorial on hybridizable discontinuous galerkin (HDG) for secondorder elliptic problems. In Advanced Finite Element Technologies, pages 105–129. Springer, 2016.
 [57] R. Sevilla, L. Rees, and O. Hassan. The generation of triangular meshes for NURBS–enhanced FEM. International Journal for Numerical Methods in Engineering, 108(8):941–968, 2016.
 [58] C.W. Shu. Discontinuous Galerkin methods for timedependent convection dominated problems: Basics, recent developments and comparison with other methods. In Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations, pages 369–397. Springer, 2016.
 [59] S. Sirivithayapakorn and A. Keller. Transport of colloids in saturated porous media: A porescale observation of the size exclusion effect and colloid acceleration. Water Resources Research, 39(4):WR001583, 2003.
 [60] S. Soghrati and R. A. Merel. NURBS enhanced HIFEM: A fully meshindependent method with zero geometric discretization error. Finite Elements in Analysis and Design, 120:68–79, 2016.
 [61] P. Šolín, J. Červenỳ, and I. Doležel. Arbitrarylevel hanging nodes and automatic adaptivity in the –FEM. Mathematics and Computers in Simulation, 77(1):117–132, 2008.
 [62] C. Wang. On Stokes slip flow through a transversely wavy channel. Mechanics Research Communications, 38(3):249–254, 2011.
 [63] Z. J. Wang, K. Fidkowski, R. Abgrall, F. Bassi, D. Caraeni, A. Cary, H. Deconinck, R. Hartmann, K. Hillewaert, H. T. Huynh, et al. Highorder CFD methods: current status and perspective. International Journal for Numerical Methods in Fluids, 72(8):811–845, 2013.
 [64] Z. Q. Xie, R. Sevilla, O. Hassan, and K. Morgan. The generation of arbitrary order curved meshes for 3D finite element analysis. Computational Mechanics, pages 1–14, 2013.
 [65] O. Zienkiewicz, J. Zhu, and N. Gong. Effective and practical version adaptive analysis procedures for the finite element method. International Journal for Numerical Methods in Engineering, 28(4):879–891, 1989.