A unified approach for a posteriori high-order curved mesh generation using solid mechanics

The paper presents a unified approach for the a posteriori generation of arbitrary high-order curvilinear meshes via a solid mechanics analogy. The approach encompasses a variety of methodologies, ranging from the popular incremental linear elastic approach to very sophisticated non-linear elasticity. In addition, an intermediate consistent incrementally linearised approach is also presented and applied for the first time in this context. Utilising a consistent derivation from energy principles, a theoretical comparison of the various approaches is presented which enables a detailed discussion regarding the material characterisation (calibration) employed for the different solid mechanics formulations. Five independent quality measures are proposed and their relations with existing quality indicators, used in the context of a posteriori mesh generation, are discussed. Finally, a comprehensive range of numerical examples, both in two and three dimensions, including challenging geometries of interest to the solids, fluids and electromagnetics communities, are shown in order to illustrate and thoroughly compare the performance of the different methodologies. This comparison considers the influence of material parameters and number of load increments on the quality of the generated high-order mesh, overall computational cost and, crucially, the approximation properties of the resulting mesh when considering an isoparametric finite element formulation.


Introduction
The performance of high-order discretisation methods for the simulation of various problems in science and engineering has been the object of intensive research during the last two decades [32,41,44,67]. These methods have the potential to offer an increased level of accuracy with a reduced number of degrees of freedom and, more importantly, a reduced computational cost [16,35,63].
The potential of high-order unstructured methods has been intensively studied by the computational fluid dynamics (CFD) community in the last decade due to their inherent ability to accurately predict the behaviour of complex high Reynolds number flows [37,45,48,74]. It is also well known that low-order methods are highly dissipative and extremely refined meshes are required to properly resolve the propagation of vortices over long distances. The advantages of high-order methods have also attracted the attention of researchers working in wave propagation problems (e.g. acoustics and electromagnetics) due to their low dispersion and dissipation compared to low-order methods [2,7,31,46,47,64]. In particular, the high-order discontinuous Galerkin method has become popular in this area due to its ability to propagate waves over long periods of time with a reduced computational cost compared to alternative low-order methods [13,15,38,40,42].
The use of curved elements is nowadays accepted to be crucial in order to fully exploit the advantages of high-order discretisation methods [5,19,43,49,61,62,70,78], but until relatively recently, the challenge of automatically generating high-order curvilinear meshes has been an obstacle for the widespread application of high-order methods [73]. Methods to produce high-order curvilinear meshes are traditionally classified into direct methods and a posteriori methods [20,21]. Direct methods build the curvilinear highorder mesh directly from the CAD boundary representation of the domain whereas a posteriori approaches rely on mature low-order mesh generation algorithms to produce an initial mesh that is subsequently curved using different techniques, such as local modification of geometric entities [20,21,50,65,66], solid mechanics analogies [58,77] or optimisation [26,71].
Within the category of a posteriori approaches, the solid mechanics analogy first proposed in [58] has become increasingly popular. The main idea is to consider the initial, loworder, mesh as the undeformed configuration of an elastic solid. High-order nodal distributions are then inserted into all of the elements and then the nodes over element edges/faces in contact with the curved parts of the boundary are projected onto the true CAD boundary. The displacement required to move the nodes onto the true boundary is interpreted as an essential boundary condition within the solid mechanics analogy. The solution of the elastic problem provides the desired curvilinear mesh as the deformed configuration. The initial approach proposed in [58] used a non-linear neo-Hookean constitutive model. Several attempts to reduce the computational cost of this approach have been proposed based on a linear elastic analogy, see [1,77]. It is clear that when large deformations are induced to produce the deformed curvilinear high-order mesh, a linear elastic model can result in non-valid elements due to the violation of the hypothesis of small deformations. In order to alleviate this problem, it is possible to split the desired (potentially large) displacement of boundary nodes into smaller load increments. Other approaches to increase the robustness of the linear elastic analogy have been recently introduced, see for instance [56], where pseudo thermal effects are introduced. It is worth noting that mesh moving strategies based on an elastic analogy have also been proposed and successfully used with a proven track record of robustness in the low-order context [39,68,69].
Although some approaches have been demonstrated to be capable of producing curvilinear high-order meshes of highly complex geometrical configurations, including anisotropic boundary layer meshes around a full aircraft configuration [77], a comparison of the proposed solid mechanics analogies has not been investigated. With this comparison in mind, in this paper a unified theoretical and computational solid mechanics approach is proposed. This formulation encompasses the linear and non-linear formulations proposed in [58] and [77], respectively. In addition, a new incremen-tally linearised elasticity formulation, not previously applied to generate curvilinear high-order meshes, is proposed within this unified approach. Different distortion measures are considered in order to evaluate the quality of the generated meshes for the different formulations, analysing the effect of material parameters, load increments, computational cost and, more importantly, the approximation properties of the resulting high-order mesh.
The paper is organised as follows. In Sect. 2, the fundamentals of non-linear continuum mechanics are briefly revisited by following some recent developments in [11,12], where the kinematics of the non-linear continua and the principle of virtual work for a displacement-based formulation, in material and spatial settings, are presented. The new consistent incrementally linearised approach is detailed in Sect. 3 and the material characterisation for all the different formulations is described in detail in Sect. 4. Using the derivation of all the formulations from an energy principle, a range of quality measures are proposed in Sect. 5 and their relations with existing quality indicators is briefly discussed. Finally, Sect. 6 presents a number of numerical examples both in two and three dimensions and an extensive comparison of performances of the different formulations is presented. The examples include geometries appearing in a range of areas of computational mechanics, e.g. computational solid mechanics, CFD and computational eletromagnetics. Meshes are produced for a variety of degrees of approximation and for interior and exterior domains, illustrating the potential of the proposed approach.

Non-linear continuum mechanics 2.1 Kinematics
Let us consider the motion of a continuum from its initial undeformed (or material) configuration Ω 0 ⊂ R d , with boundary ∂Ω 0 and outward unit normal n 0 , into its final deformed (or spatial) configuration Ω ⊂ R d , with boundary ∂Ω and outward unit normal n, where d represents the number of spatial dimensions. In the context of curved mesh generation, the initial (undeformed) configuration Ω 0 represents a linear mesh with planar faces (edges in two dimensions) and the final (deformed) configuration Ω represents the final curved high-order mesh, as illustrated in Fig. 1. The motion is described by a mapping φ which links a material particle from material configuration X to spatial configuration x according to x = φ(X).
The following well-known strain measures can be introduced. Firstly, the two-point deformation gradient tensor or fibre-map F, which relates a fibre of differential length from the material configuration d X to the spatial configuration d x, namely where ∇ 0 denotes the gradient with respect to material coordinates. Secondly, the two-point co-factor or adjoint tensor or area map H which relates an element of differential area from material configuration d A (co-linear with n 0 ) to spatial configuration d a (co-linear with n), namely where denotes the tensor cross product between two second order tensors, introduced in [18] and recently brought into the context of solid mechanics in [10][11][12]30]. 1 Finally, the Jacobian J or volume-map of the deformation which relates differential volume elements in the material configuration dΩ 0 and the spatial configuration dΩ, is introduced as The fundamental strain measures {F, H, J }, also illustrated in Fig. 1, encode the essential modes of deformation.

The principle of virtual work in material and spatial settings
While a myriad of methodologies can be applied to solve for the deformation of a continuum described by the motion map φ, such as optimisation and projection techniques [20,21,25,1 The cross product ( ) between two tensors is defined as where ξ is the third order permutation tensor, see [11,12] for a list of properties. 28,50,51,57,65,66], in the context of solid mechanics, the deformation of a continuum from its undeformed configuration to its deformed configuration can be posed as a problem of minimisation of the total potential energy Π , subjected to certain desired constraints [9,34,52]. In other words, the displacement of a deformable body can be obtained by finding the stationary condition of the total potential energy, also called the principle of virtual work (or variational principle), of an assumed strain energy density. A valid expression for strain energy density must fulfill certain requirements to guarantee the existence of minimisers [4,14]. In addition, to satisfy objectivity (i.e. invariance with respect to rotations or coordinate transformation), the strain energy density must be constructed in terms of the right Cauchy-Green strain tensor C = F T F [9]. 2 Hence, in our setting, we can write where is worth noting that, as usual in the context of a posteriori curved mesh generation, no external forces are considered and Dirichlet boundary conditions are imposed on the whole boundary, corresponding to the displacement needed to place the high-order nodes on the CAD boundary entities, at positionx.
The stationary condition of the total potential energy in (4) obtained after linearisation with respect to virtual displacements v, leads to the principle of virtual work in the material configuration where S is the symmetric second order Piola-Kirchhoff stress tensor and DC [v] denotes the first directional derivative of the right Cauchy-Green strain tensor along the direction of virtual displacements v. Equation (5) results in a system of non-linear equations that need to be solved through an iterative scheme. The tangent operator required to facilitate convergence of the non-linear iterative scheme (e.g. Newton-Rapshon) can be determined by computing the second directional derivative of the total potential energy where the material fourth order tangent elasticity tensor is defined as and D 2 C[v; w] denotes the second directional derivatives of the right Cauchy-Green strain tensor.
In order to establish a unified approach for various solid mechanics approaches discussed in this work, it is convenient to re-express (5) in the spatial configuration where σ is the Cauchy stress tensor, ε(a) is the small strain tensor given by ε(a) = 1 2 (∇ a + (∇ a) T ) and ∇ is the gradient operator in the spatial configuration such that ∇ 0 (a) = ∇(a)F. Analogously, the second directional derivative of the total potential energy (6) can be re-expressed in the spatial configuration as where c is the spatial fourth order tangent elasticity tensor obtained as In this setting, the first and second terms in the right hand side of (9) yield the constitutive and the geometric/initial stiffness components, respectively.

A consistent incrementally linearised solid mechanics approach
The standard non-linear solid mechanics methodology described in the previous section can be proven costly (due to the iterative nature of the solution finding process) when the ultimate goal is solely to deform a mesh in order to conform to the exact geometry. Alternative solid mechanics methodologies have been developed in the past, based on a variety of linearised elasticity approaches [1,56,77].
x n = φ n (X) It is worth emphasising that, to guarantee and/or maintain previously mentioned mathematical requirements for the linearised strain energy density, a linearised solid mechanics approach must emanate from an underlying non-linear variational principle, as the notion of objectivity and polyconvexity cannot be invoked in small strains. This is typically achieved by consistent linearisation of the total potential energy (4) through a Taylor series expansion. To illustrate this, let us consider the total potential energy in (4), cast in the form of an iterative (Newton-Raphson) scheme x n+1 on ∂Ω 0 and x n+1 = φ n+1 (X) is the position vector of the material points at increment n+1, which can be evaluated through an incremental displacement u superimposed on the deformed configuration at increment n, i.e. x n+1 = x n + u, as illustrated in Fig. 2. At increment n, the current position vector x n , the state of deformation gradient F n and, subsequently, the Cauchy-Green strain C n are fully known. In a non-linear regime, the motion of the continuum from n to n+1 is solved iteratively, as the amount of displacements, the state of deformation gradient F n+1 and the Cauchy-Green strain C n+1 cannot be determined explicitly. However, in the context of high-order curved mesh generation, it is convenient to approximate (11) through a Taylor series expansion of the form embedded in the definition of the new total potential energy (Π lin ) in (12) are the first and second directional derivatives of the non-linear total potential energy (4), where the virtual and incremental variations, namely v and w, are now replaced with u. Hence, unlike (5), (12) is fully and consistently linearised in u. Furthermore, notice that the first term in the integrand in (12) is a constant term describing the state of strain energy density at increment n, which vanishes at the moment of computing the stationary point of (12). In fact, the stationary condition of the linearised total potential energy (12) is identical to the stationary condition of the following potential energy Substituting for (8)(9) in (13), we obtain the spatial form of the linearised total potential energy as where the subscript n denotes the state of deformation, stresses, tangent elasticity and the volume at increment n, namely ε n , σ n , c n and Ω n . In addition, ∇ n represents the spatial gradient operator at increment n. The stationary condition of (14), obtained after the linearisation with respect to the virtual displacement v, leads to the principle of virtual work It is worth noting that, in the right hand side of (15), the first term R n corresponds to the residual stresses, the second term C n to the linearised constitutive stiffness term and the last term G n to the geometric stiffness term. The emergence of the geometric stiffness term is due to consistent linearisation of the non-linear total potential energy, which would not have appeared, had the starting point not been chosen to correspond to a non-linear total potential energy, as also described in [36], p. 104 and in [52]. As will be seen in the numerical examples, in the context of high-order curved mesh generation, the geometric stiffness term, stiffens the interior elements of the computational mesh against severe distortion, hence producing meshes with better quality. Note that unlike in the non-linear analysis, since (15) is linear in u, a further linearisation is not required. It is evident that, if a single increment is used to reach the final configuration, i.e. when n = 0, the equations of classical linear elasticity are recovered. It should be emphasised, that in the context of linearised approaches there are a multitude of heuristic formulations which do not necessarily come from a variational principle. In fact, apart from the case when all the terms {R n , C n , G n } are present in the principle of virtual work, all the other formulations based on the combinations of these terms lead to non-consistent formulations. For the sake of completeness, the four distinct linearised cases used in the literature of curved mesh generation are identified in the following; [see also Table 1 for a schematic comparison of these cases and their associated computational requirement]: 1. When the state of deformation at increment n + 1 is obtained by computing the residual stresses, the constitutive stiffness and the geometric stiffness at the previous deformed configuration i.e. based on {R n , C n , G n }. This consistent incrementally linearised methodology for high-order curved mesh generation, is presented in this paper for the first time. 2. When the residual stresses at increment n+1 are obtained from the previous deformed configuration based on R n , but the constitutive stiffness is evaluated at the initial undeformed (or stress-free) configuration i.e. C 0 and the geometric stiffness term is absent from the formulation. The technique developed by [56] falls into this category. 3. When both the residual stresses and constitutive stiffness at increment n + 1 are computed based on the initial undeformed configuration i.e. R 0 and C 0 and the geometric stiffness term is absent from the formulation, but the geometry itself is updated incrementally such that x n+1 = x n +u. This approach has been pursued in [1,77] and from here onwards we will refer to this approach as the incremental linear elastic approach. 4. When n = 0 or in other words, when the residual stresses, constitutive stiffness and geometric stiffness are evaluated once at the initial undeformed configuration {R 0 , C 0 , G 0 } i.e. particularisation to the case of classical linear elasticity.

Hyperelasticity and material characterisation for different solid mechanics formulations
In the previous section, the principle of virtual work was introduced in terms of the stresses in the body. These stresses result from the deformation of the body expressed in terms of the fundamental kinematic measures {F, H, J }. To close the system of equations, introducing a relationship between stresses and the kinematic measures is necessary. Known as the constitutive relations, these relations depend on the type of material under consideration (i.e. isotropic, anisotropic, compressible, nearly incompressible and so on). Typically, these constitutive relations are obtained from the linearisation of the strain energy density of the material. Apart from satisfying the requirements of existence of minimisers, objectivity (and favourably polyconvexity), the choice of material generally imposes further physical requirements on the strain energy density. For an isotropic material, the strain energy density can be written as a function of three independent invariants as Ψ (C) = Ψ iso (I 1 , I 2 , I 3 ), I 1 = I : C, where the three isotropic invariants I 1 , I 2 and I 3 are indeed further related to the fundamental kinematic measures as 3 Furthermore, for a transversely isotropic material, the strain energy density can be expressed as Ψ (C) = Ψ aniso (I 1 , I 2 , I 3 , I 4 , I 5 ), I 4 = N · C N, 3 Note that for plane strain problems, the first two isotropic invariants are identical i.e. I 1 = I 2 .
where the two transversely isotropic invariants I 4 and I 5 are indeed related to the fundamental kinematic measure F as (19) where N is the unit Lagrangian vector characterising the direction of transverse isotropy. It is worth emphasising that (17) and (19) represent a set of independent invariants that can be used to construct isotropic and transversely isotropic strain energy density expressions, respectively. 4 Any other invariant used to construct the strain energy density, would be a combination of the aforementioned invariants. As will be discussed later, these invariants play a key role in quantifying the quality of generated curvilinear meshes.
To establish a unified approach for the different solid mechanics formulations presented in the previous section and to further facilitate a comparison among them, it is essential that material parameters are calibrated such that the strain energies (and consequently the stresses and the constitutive tensors) for the different formulations are identical in the initial undeformed configuration. To this end, in this section, we discuss characterisation of material constants through an example of a hyperelastic neo-Hookean model. More sophisticated strain energies accounting for near incompressibility and transverse isotropy are also considered, which for the purpose of brevity are only listed in Table 2.

The nonlinear hyperelastic case
Let us consider a compressible neo-Hookean model which is also considered, for instance in [58] in the context of high-order curvilinear mesh deformation. The strain energy density of the material is given by where μ and λ are two material constants. Note that this model can be obtained by setting β = 0 and α = μ 2 in the more sophisticated compressible Mooney-Rivlin model presented in Table 2. Following the procedure outlined in the previous section, the Cauchy stress tensor and the tangent elasticity tensor can be obtained as Nearly incompressible Mooney-Rivlin (NI-MR) Transversely Isotropic Hyperelastic (TI) a For all the analyses we fix α = β b For all the analyses we fix α = μlin 2 c C 11 , C 13 , C 33 , G A , E and ν are the components of transversely isotropic linear elastic material, where the subscript A denotes the direction of anisotropy. For all the analyses we fix α = β ↔ implies relationship between non-linear and linear material constants after calibration. Also note that where b = F F T is the left Cauchy-Green strain tensor and I is the symmetric fourth order identity tensor [I] i jkl = δ ik δ jl + δ il δ jk , where δ mn denotes the Kronecker delta.

The classical linear elastic case
For the classical compressible linear elastic constitutive model, the strain energy is defined as where λ lin and μ lin represent the well-known Lamé parameters. The Cauchy stress tensor and the tangent elasticity tensor are then obtained as σ lin = λ lin trε I + 2μ lin ε (24) c lin = λ lin I ⊗ I + μ lin I. (25) It is worth noting that the approach pursued in [56,77] corresponds to a linear elastic approach with only the geometry being updated at each increment x = x n .
Comparison of the tangent elasticity tensor of the neo-Hookean model (22) evaluated at the origin (i.e. J = 1) to its linear elastic counterpart (25), a relationship between material constants can be defined as In the case of the neo-Hookean material model the relationship between material constants is one-to-one, whereas for more complex material models dependent upon more than two constants, one can arrive at correlations between those and the Lamé constants, preferably having to fix some of those constants. The calibration of material constants of a nonlinear or linearised energy functional against the linear elastic model is key to the comparison of the approaches. In practice, normally the Young's modulus (E) and Poisson's ratio (ν) of the material are provided, which for the threedimensional and plane strain isotropic cases are related to the Lamé constants as To facilitate a comparison between material models for producing the higher-order curvilinear meshes, all isotropic materials described in Table 2 can be expressed in terms of the Poisson's ratio ν using (26). For transversely isotropic material, the relationship between material constants and ν is given in Table 2, cf. [8].

Mesh quality measures
Quality or distortion measures are traditionally used both in a low and high-order finite element context in order to quantify the approximation properties induced by a computational mesh. In a standard high-order finite element formulation, measures involving the Jacobian of the isoparametric mapping have been extensively used [29,77], in particular the so-called scaled Jacobian. This measure only quantifies volumetric deformations and alternative measures that exploit different modes of deformation and account for shape, skewness and degeneracy of elements have only been recently considered [27,28]. However, it is worth noting that not all of these quality measures can be regarded as independent quantities.
In the unified solid mechanics approach proposed here, due to the derivation of all the approaches from an energy principle, we propose five quality measures that are defined in terms of the invariants of (17) and (19), used to construct the strain energy. The quality measures for a generic element e are where R denotes the reference element employed in the isoparametric formulation, with local coordinates ξ . If necessary, further quality measures can be obtained through a linear combination of the invariants I j which will be independent of the geometrical parametrisation.
In practice, the invariants are evaluated at a discrete set of points within the reference element, usually the quadrature points that will be employed during a computational simulation. For the numerical examples presented here, a quadrature rule is used that integrates polynomials of degree up to 2 p, where p is the order of approximation.
In order to obtain a representative quality measure for a given computational mesh, a variety of statistical data can be reported, such as the mean quality or the standard deviation. However, in the numerical examples presented here, the mesh quality is defined by computing the minimum over all the elements, namely Q j = min e {Q e j }. Despite this being the least favourable choice, it is well known that a few low quality elements can substantially deteriorate the overall quality of a finite element simulation, specially if these elements are near a curved boundary. Several numerical examples in two and three dimensions are used in the next section to evaluate the performance of different approaches for a posteriori mesh generation. The objective is to produce meshes where the minimum quality is as high as possible as this will provide the better approximation properties of a high-order finite element solver.
From the mesh distortion point of view, the first quality measure Q 1 , quantifies fibre deformation (for instance, distortion of the edges of an element), the second quality measure Q 2 , quantifies surface deformation (for instance, distortion of the faces of an element) and the third quality measure Q 3 , quantifies volumetric deformation (distortion of the element itself). In fact, it is worth noting that the scaled Jacobian corresponds to the quality measure Q 3 . For simplicial elements (i.e. triangles and tetrahedra), this measure is identical to the Jacobian of the deformation gradient tensor J because the isoparametric mapping for a simplicial elements with planar faces (or edges in two dimensions) is affine. This result is valid because, in the context of a posteriori high-order mesh generation, the undeformed configuration typically corresponds to a mesh formed by elements with planar faces (or edges).
The quality measures Q 4 and Q 5 , based on the two anisotropic invariants, quantify the distortion in the direction of anisotropy. These measures can only be utilised when the internal energy of the material is anisotropic and, since in the context of curved mesh generation this is not often the case, their usage remains limited. Moreover, anisotropic quality measures are typically dependent on the geometrical parametrisation.
Finally, it should be emphasised that, in contrast to the non-linear approach, the solution of the incrementally linearised problem in (13), does not correspond to the minimisation of the total potential energy (11) with respect to the fundamental strain measures {F, H, J } per se, but rather with respect to the incrementally linearised versions of these quantities. Furthermore, it is easy to identify that for plane strain problems, the first and the second invariants are indeed identical i.e. I 1 = I 2 which in turn translates into the corresponding quality measures being identical Q 1 = Q 2 . This is only true for two-dimensional plane strain problems.

Examples
This section presents a detailed comparison of the various solid mechanics formulations considered in this work (refer to Table 1) for the a posteriori generation of high-order curvilinear meshes. The comparison focuses on the advantages and disadvantages of the various formulations, the influence of the material parameters, the degree of approximation obtained by using two and three dimensional examples and the monitoring of different quality measures. In this work, the only material parameter that is varied is the Poisson's ratio (ν). Notice that as detailed in [77], the Young's modulus has no real effect on the resulting high-order meshes because only Dirichlet boundary conditions are considered. Therefore, in all the examples we consider E = 10 5 , E A = 5E 2 and G A = E 2 and the Poisson's ratio is selected within the interval [0.001,0.495].
To simplify the presentation, the following acronyms are utilised: for incremental linear elastic and consistent incrementally linearised formulations, respectively. When these acronyms are not used in conjunction with a material model, the formulation should be assumed to correspond to a fully non-linear analysis. For the sake of brevity the names of the following two material models are also shortened to -Nearly Incompressible Mooney-Rivlin (NI-MR) -Transversely Isotropic (TI) In all the examples, the linear system of equations resulting from the finite element discretisation, is solved using UMFPACK [17] and the Multi-frontal Massively Parallel Solver (MUMPS) [3] for the systems that result in twodimensional and three-dimensional examples, respectively. The non-linear systems are solved using a standard Newton-Raphson algorithm where the tolerance is set to 10 −5 . Finally, it is worth noting that a standard isoparametric finite element formulation is considered throughout this work, using Lagrange polynomials with optimal distribution of nodes for interpolation [77] and also the optimal quadrature points for triangles and tetrahedra reported in [75] are considered for numerical integration.
The developed code, called PostMesh, has been released as an open-source software under MIT license and is available through the repository https://github.com/romeric/ PostMesh.

Mesh of a mechanical component
The first example considers an isotropic mesh of a mechanical component. The geometry is given by 28 NURBS curves describing the boundary of the domain as depicted in Fig. 3.
The initial linear triangular mesh is shown in Fig. 4a, having 192 elements, 129 nodes and 68 boundary edges. The produced mesh using the ILE approach for a degree of approximation p = 5 is shown in Fig. 4b, having 2569 nodes.
A detailed view of four high-order meshes produced using the same ILE approach is shown in Fig. 5, showing the large distortion that is induced by the projection of the boundary nodes over the true boundary. In addition, the better approximation of the true boundary shown in Fig. 3 as the polynomial order is increased, can be clearly observed. Figure 6 shows a comparison of the quality of the generated curvilinear meshes using linear, incrementally linearised and non-linear approaches. In all cases, the deformation of the boundary has been imposed using five increments and the minimum scaled Jacobian is used as a quality measure.  It can be observed that the quality of the meshes produced with the ILE isotropic and CIL neo-Hookean approaches is almost identical, although the CIL neo-Hookean approach seems to provide better quality for high-order (e.g. p = 5) approximations and for values of the Poisson's ratio near the incompressible limit. Despite this difference, both approaches are able to produce high quality meshes for any degree of approximation tested. In contrast, the (nonlinear) neo-Hookean approach fails to produce a high-order mesh for high-order approximations, except for a few cases where a low-quality mesh for p = 4 is produced. The quality of the produced meshes for lower order approximations (i.e., p = 3,4) is similar to the quality produced by the ILE isotropic and CIL neo-Hookean approaches, but it is worth noting that the (non-linear) neo-Hookean approach also fails in the nearly incompressible region, whereas the ILE isotropic and CIL neo-Hookean approaches produce the best quality meshes in this scenario.
It should also be noted that, unlike the linearised approaches wherein the internal nodes of the mesh move proportionally to the boundary nodes, in the non-linear approach the internal nodes can move arbitrarily within the element, and this can in turn affect the quality and approximation property of the produced meshes. In a purely displacement-based formulation, it is not feasible to restrain the movement of internal nodes to a desired proportion. In this context, higher order gradient theories [6,23,55] and more elaborate mixed formulations [11,60], offer a potential future research direction.
Next, we compare the effect of material models presented in Table 2 on the quality of generated meshes, for all the three approaches. Figure 7 shows the quality (min- For the transversely isotropic model, the negative x-axis is chosen as the direction of anisotropy for the interior elements. For the elements in the boundary, the direction of anisotropy is computed to be perpendicular to the unit normal to boundary edge. This technique is customary in the field of fibre-reinforced composites.
The results show that the quality displayed with neo-Hookean, Mooney-Rivlin and nearly incompressible materials is almost identical for any value of the Poisson's ratio, whereas a different behaviour is obtained for the transversely isotropic model. The best quality is obtained with the ILE TI model and with a Poisson's ratio near 0.5. However, it is worth emphasising that a small variability of the quality is obtained in all cases as all simulations provide a high-order mesh with quality belonging to [0.67,0.77].
Next, the same comparison is performed for higher order approximations, but the results with a Mooney-Rivlin and nearly incompressible models are omitted because, in all cases, the results are almost identical to those obtained with a neo-Hookean model. Figure 8 shows the quality (minimum scaled Jacobian) as a function of the Poisson's ratio when a polynomial approximation of degree p = 3 is considered.
A different trend is observed, when comparing the results with p = 3 to the results with p = 2 displayed in Fig. 7. With p = 3 the quality of the mesh improves as the Poisson's ratio is increased, providing the best results always when the incompressible limit is approached. This behaviour is expected in general because when the Poisson's ratio is taken near 0.5, the imposed displacement on the boundary induces a larger displacement of the interior nodes. In contrast, when a value of the Poisson's ratio near 0 is considered, the imposed displacement on the boundary induces small displacement on the interior nodes, resulting in more distorted elements (i.e. reduced quality elements). The reason why this expected behaviour was not obtained with p = 2 is attributed to the lack of resolution of the displacement field when the coarse mesh considered here, see Fig. 4, is employed with (a) (b) Fig. 8 Minimum scaled Jacobian of the generated meshes as a function of the Poisson's ratio for p = 3 using the ILE isotropic, CIL and non-linear approaches with different material models. a Neo-Hookean. b TI a quadratic approximation. In fact, further simulations not reported here for brevity confirm that with a finer mesh the expected trend is obtained even with a degree of approximation p = 2.
In addition, the results show that the quality of the meshes produced with ILE isotropic, CIL and non-linear approaches is almost identical if a neo-Hookean model is considered, whereas the use of a transversely isotropic model reveals some differences between the three approaches. The results demonstrate the significance of chosing a well-defined material model like neo-Hookean (with a quality reported near 0.85), in contrast with a transversely isotropic model (quality reported below 0.6 for any value of the Poisson's ratio), whose limitations would be discussed shortly. It is worth emphasising that the quality obtained with the Mooney-Rivlin and the NI-MR models is almost identical to that produced by the neo-Hookean model, so that any of the three models is equally suitable to produce high quality meshes in this example, since all these material models are mathematically well-defined.
Finally, Fig. 9 shows the quality (minimum scaled Jacobian) as a function of the Poisson's ratio when a polynomial approximation of degree p = 6 is considered. The conclusions that are implied by the results are similar to those obtained from the simulation with p = 3. First, the quality obtained with the neo-Hookean model is similar for the ILE and CIL approaches whereas some differences are observed for the transversely isotropic model. However, in this case the non-linear approach is not able to converge, as already mentioned and shown in Fig. 6. The quality of the produced meshes increases as the Poisson's ratio approaches 0.47 and the best results are obtained when a neo-Hookean (equivalently compressible or nearly incompressible Mooney-Rivlin) model is considered. It is worth mentioning that this example shows a slight drop in the quality of the mesh as the Poisson's ratio increases from 0.47 to 0.49. It should be noted that imposing the material to be incompressible in this example is not physically possible because the initial and deformed configuration have a pre-defined (and non-equal) volume as shown in Fig. 6. Therefore, the results suggest that the Poisson's ratio should be carefully selected near the incompressible limit, but preferably of a value to ensure that some level of compressibility is allowed, for instance 0.45. This behaviour is only observed with p = 6 because for lower order approximations there is a lack of resolution to capture the displacement field.
The analysis for the different approaches and material models is summarised in Fig. 10. This figure shows the mean and standard deviation of the scaled Jacobian for the ILE isotropic, CIL and non-linear approaches with different materials and degrees of approximation.
It can be concluded that the choice of material model does not play a major role, as long as the model is well-defined. As hinted before, unlike the other material models, the transversely isotropic (TI) material defined in Table 2, does not correspond to a polyconvex energy functional, or more specifically, the invariant N · C 2 N = (F T F N) · (F T F N), is not convex with respect to F, and hence under highly large deformations, the model experiences instabilities in the form of loss of ellipticity which can manifest through shear-bands, fibre kinking if under compression or fibre de-bonding if under stretch; cf. [53] and [54], for an intensive study on the loss of ellipticity for this invariant. The latter two phenomena (fibre kinking and de-bonding) also hold true for the transversely isotropic linear materials. As a consequence, it can be observed that the mean quality of the high-order meshes generated with a transversely isotropic material deteriorates Overall, the ILE isotropic approach is found to be the most robust, providing the best or near the best mean quality for all orders of approximation. Also, it is worth noting that for all material models the standard deviation grows as the order of approximation is increased, implying that a good choice of the Poisson's ratio is more important as the order of approximation is increased.
Next, the effect of the Poisson's ratio, the different approaches and material models on the condition number of the system matrix is illustrated in Fig. 11. The condition number is computed using the lower bound one-norm estimate of Higham [33].
Again, the results show that the condition number with neo-Hookean, Mooney-Rivlin and NI-MR is almost identical for any value of the Poisson's ratio, whereas a different behaviour is obtained for the transversely isotropic model. In all cases the condition number increases as the Poisson's ratio approaches the incompressible limit but it is worth noting that a slightly lower condition number is obtained when the transversely isotropic model is considered, irrespective of the use of ILE isotropic, CIL and non-linear approaches. This is inherently due to anisotropic nature of the model, as the deformation is not homogenous in every direction and hence the effect of Poisson's ratio is not equally pronounced for this model. The results with other degrees of approximation are omitted, as exactly the same behaviour is observed. The last analysis is aimed to compare the computational cost of each formulation with different material models and different orders of approximation p. As it is not feasible to a priori know the number of iterations required by the non-linear approach to converge, a comparison of the actual computing time is considered here.
The computational code is written in Python using Python scientific stack such as NumPy and SciPy, with the core routines implemented in C, C++ and Cython. With the compiled code, we utilise X86-64 explicit vectorisation by dispatching most of the computation at quadrature points to SIMD intrinsic calls. A parallel map model is used for computation of local elemental matrices and multi-core assembly of the system of equations, as the interprocess communication, at this stage is minimal. A similar strategy for assembling the elemental matrices is adopted and discussed in detail in [72]. The code is linked against the optimised Open-BLAS library. For the purpose of comparison, we exclude the pre-processing and post-processing and only measure the run-time of the assembly of tangent stiffness, application of Dirichlet boundary conditions and the solution of the system of equations.
Whilst theoretically, the non-linear approach should cost number of iterations × number of increments times more than the linear model, in practice, due to differences in the sparsity pattern and condition number of the system as well as CPU warm-up and pipelining, this is not often the case. In fact, comparison of non-linear against linear approaches is analogous to cold versus hot benchmarking, in that with a higher number of iterations, the processor becomes progressively more accurate with branch prediction and guessing jmp operations, which helps improve processor pipelining. On the other hand, for highly non-linear problems, with every iteration of the non-linear analysis the condition number increases, hence impacting the run-time. With this in mind, we report the geometrical mean of 100 runtimes, excluding the timing for the first 10 runs. An in-house tool similar to Google Benchmark is used for time measurements. For all time measurements, parallelisation has been turned off. Material data and p are deliberately chosen such that the non-linear analysis would converge. The analysis corresponds to ν = 0.4 with other parameters remaining constant as before. Figure 12 shows the CPU time using the three formulations and different material models when the boundary displacement is imposed using five load increments using a polynomial approximation of degree p = 2 and p = 3. The CPU time has been normalised with respect to that of classical linear elasticity (i.e. one increment). It is important to note that, due to the small size of the problem, all systems are solved using UMFPACK, the cost of solver is negligible and the condition number of the system does not adversely affect quadratic convergence of Newton-Raphson. Furthermore, for a problem of this size, a portion of the computational time corresponds to the overhead of function calls i.e. argument parsing at interpreted level and push/pop operations at compiled level, even after some aggressive inlining.
The computational cost associated with the different material models is clearly related to their tangent operator. In the case of linear elasticity, the tangent operator can be computed at the pre-processing stage. For a neo-Hookean model the two fourth order identity tensors I ⊗ I (δ i j δ kl ) and I (δ ik δ jl + δ il δ jk ), appearing in the tangent operator, are compile time constants. For Mooney-Rivlin, nearly incompressible Mooney-Rivlin and transversely isotropic hyperelastic models, the dyadic products in the tangent operators [see Table 2] are run-time variables and their computation is not always cache-friendly due to unavoidable strided indexing. In fact for complex material models identifying the optimal networks of tensor contraction is not trivial [24,59]. Using Voigt notation and further permutations, these dyadic products can be transformed to further gemm calls, which eventually may or may not be beneficial. Moreover, the nearly incompressible and transversely isotropic hyperelastic models require computation of co-factors H and H H T , at every quadrature point which are all O(n 3 ) in computational complexity.
In this example, ILE isotropic approach is found to be the most competitive. This allows to conclude that, for this example, the ILE isotropic approach provides both the best quality and the lowest computational cost compared to other approaches and material models. Furthermore, one should note that the mesh qualities reported here are not indicative of the maximum quality that can be obtained, as the number of load increments is rather kept fixed to facilitate an impartial comparison between different approaches. Finally, it is worth mentioning that although CPU time measurements are always dependent on the implementation, the results reported here provide a qualitative indication of the higher cost associated to a non-linear approach. The CPU timing, together with the already discussed convergence difficulties of non-linear approaches for high-order approximations, clearly provides an indication of the limited scope of such an approach for a posteriori high-order mesh generation.

Mesh around the SD7003 aerofoil
The second example considers anisotropic boundary layer meshes around the SD7003 aerofoil with different levels of stretching in the boundary layer. The detailed view near the The produced mesh using the ILE approach for a degree of approximation p = 5 is shown in Fig. 13 (b), having 27,410 nodes.
It should be mentioned that, similar to the previous example, it was found that the choice of material model does not have an effect on the quality of the curved meshes and that the transversely isotropic material shows a similar pattern of loss of ellipticity. In the light of these findings, we abandon the comparison of material models for the present example and unless otherwise stated, we only utilise the neo-Hookean model with its linearised version. In contrast, due to high level of stretching of the meshes considered here, the effect of the number of load increments on the quality of generated meshes will be investigated. Figure 14 shows the quality of the high-order meshes, measured as the minimum scaled Jacobian, as a function of the Poisson's ratio and the number of load increments for the ILE isotropic, CIL neo-Hookean and (non-linear) neo-Hookean approaches.
Once more, the non-linear approach is not able to provide a solution in all cases (i.e. for all values of the Poisson's ratio and number of load increments). In fact, when it converges, the quality of the non-linear approach is generally lower than the quality of the ILE isotropic and the CIL neo-Hookean approaches. It can also be observed that the quality of the meshes produced with the ILE isotropic and the CIL neo-Hookean approaches is almost identical, for any value of the Poisson's ratio and for any number of load increments. Finally, the results show that the best quality is obtained for value of the Poisson's ratio near the incompressible limit and ten load increments approximately. A further increase of the number of load increments does not improve substantially the quality of the meshes but it enables to obtain high quality meshes for slightly lower values of the Poisson's ratio. Figure 15 shows the same analysis for meshes with significantly higher level of stretching, namely 100 and 800, for the same degree of approximation, p = 2.
For these meshes, the non-linear approach is only able to provide a result in a few cases. In fact, further numerical experiments show that the higher the stretching, the more cases would display no convergence of the non-linear approach. In addition, the quality of the meshes produced with the ILE isotropic and the CIL neo-Hookean approaches is, again, almost identical, for any value of the Poisson's ratio and for any number of load increments, showing that the conclusions presented do not strongly depend on the level of stretching within the boundary layer.
For the incremental linear elastic and the consistent incrementally linearised approaches, the number of load increments generally improves the mesh quality but the same is not true for the non-linear approach. If phenomena such as buckling, snap-back and snap through are not expected, the non-linear approach provides the same mesh quality irrespective of the number of increments. However, in the presence of buckling, it is possible to jump through snapback/snap-through region with fewer load increments, but as the number of load increments is increased the buckling (i.e. snap-back/snap-through regions) cannot be avoided, which in the absence of an arc-length technique leads to nonconvergence of the Newton-Raphson method. Furthermore, it is possible for the Newton-Rapshon scheme to converge just prior to the onset of buckling, at the cost of losing quadratic rate of convergence due to ill-conditioning of the system which essentially emanates from nearly zero Jacobian(s).
Next, the same analysis is performed for higher orders of approximation. Figures 16 and 17 show the quality of the high-order meshes as a function of the Poisson's ratio and the number of load increments for the ILE isotropic, CIL neo-Hookean and (non-linear) neo-Hookean approaches, for a degree of approximation p = 4 and p = 6 respectively. For p = 4 the non-linear approach is unable to converge in the majority of cases. Only for a relatively low stretching, such as 50, and using one load increment, this approach provides a solution for any value of the Poisson's. ratio. When the stretching is increased to 400, this approach fails to converge even when one increment is used if the Poisson's ratio is selected near the incompressible limit. For the ILE isotropic and CIL neo-Hookean approaches the quality of the produced meshes is, once more, almost identical for any value of the Poisson's ratio, number of increments and stretching. It worth noting that for this order of approximation, the increase in stretching translates into a significant decrease in the maximum quality that can be obtained with the ILE isotropic and CIL neo-Hookean approaches.
If the order of approximation is further increased to p = 6, the problem becomes substantially more challenging and the quality of the produced meshes with either the ILE isotropic and CIL neo-Hookean approaches is significantly lower, as observed in Fig. 17. For a stretching factor of 100, a higher number of load increments is required (approximately 40) compared to previous examples and a value of the Poisson's ratio near the incompressible limit is mandatory to obtain the best quality meshes. It is also worth noting that this example shows, for the first time, a subtle difference between the ILE isotropic and CIL neo-Hookean approaches. For a value of the Poisson's ratio near the incompressible limit, the CIL neo-Hookean approach requires more load increments than the ILE isotropic approach to obtain similar quality. The conclusions for a stretching factor of 800 are similar but, as it can be observed in Fig. 17, both the ILE isotropic and CIL neo-Hookean approaches can only provide a maximum quality near 0.3. In this example, the non-linear approach shows once more the inability to converge in the majority of simulations.
To summarise, Fig. 18 shows the ratio of the scaled Jacobian with 50 load increments over the classical linear elasticity (i.e. single increment), in a logarithmic scale. Note that, due to the logarithmic nature of this measure, a factor of zero implies no improvement and, furthermore, a slight improvement in terms of this factor can imply a significant change in terms of percentage value. For instance, for p = 6, and stretching of 1600, the scaled Jacobian improves from 0.0011 for a single increment to 0.3382 for 50 increments. It can be observed, that at high p, it is crucial to increase the number of load increments to obtain good quality meshes, specially if the stretching is also high. In contrast, for low-order approximations the gain obtained by increasing the number of load increments is marginal.
To further illustrate the improvement induced by an increase on the number of load increments in the quality of the generated meshes, Fig. 19 shows a histogram of the quality for two different values of the Poisson's ratio, namely A marginal difference is observed between the ILE isotropic and CIL neo-Hookean approaches both at lower values of Poisson's ratio as well as for values near the incompressible limit. As discussed earlier, this figure shows that the non-linear model is only able to converge when a few load increments are considered. However, note that due to the presence of the geometric stiffness term in the CIL approach, the interior elements are stiffened against heavy distortion and hence the CIL approach typically produces meshes with a slightly better distribution of the quality over the computational mesh, irrespective of the minimum value for quality measures. The results also show the improvement induced by an increase of the Poisson's ratio. For instance, Fig. 19 (a) shows that the mesh contains a significant number of elements of quality 0.45 when the Poisson's ratio is 0.11 whereas the minimum quality of the mesh associated to Fig. 19 (c), with a Poisson's ratio of 0.44, the minimum quality is near 0.75.
Next, the quality of the generated meshes in terms of different measures is studied, namely the measures defined in Equation (27) that are defined in terms of the invariants in Equation (17). The two anisotropic mesh quality measures Q 4 and Q 5 are dropped from the comparison because they are only valid for the transversely isotropic material model, which has been shown to produce low quality meshes in the examples considered. It is worth emphasising that the qualities Q 1 and Q 2 are the same for two-dimensional plane strain problems. Figure 20 shows the quality Q 1 as a function of the Poisson's ratio and the number of load increments for the mesh with a stretching factor 100 and p = 6.
The results show that the quality is substantially improved when the number of load increments is increased, as previously observed with the minimum scaled Jacobian as quality Fig. 20 Quality Q 1 of the generated meshes with p = 6 as a function of the Poisson's ratio and the number of load increments with a stretching of 100. a ILE isotropic. b CIL neo-Hookean. c (Non-linear) neo-Hookean measure. However, in this example, increasing the Poisson's ratio near the incompressible limit induces a lower quality except if a large number of load increments is considered. As shown in previous examples, the ILE isotropic approach performs slightly better than the CIL approach when the Poisson's ratio is selected near the incompressible limit and the non-linear approach fails to converge in the majority of cases. Finally, by comparing Figs. 20 and 17, it is clearly observed that a value closer to one is obtained when using Q 1 instead of the minimum scaled Jacobian. This behaviour is expected because, for this problem, the deformation is primarily volumetric and the deviatoric contribution is minimal. Figure 21 shows the three quality measures Q 1 , Q 2 and Q 3 as a function of the Poisson's ratio for the mesh with p = 2, a stretching factor of 25 and using five load increments. The results confirm, numerically, that the quality measures Q 1 and Q 2 are the same for two-dimensional plane strain problems. It can also be observed that the ILE and CIL approaches produce meshes of the same quality, irrespectively of the measure used. In addition, the results illustrate that the quality measure Q 1 (and Q 2 ) is less influenced by changes on the Poisson's ratio, compared to Q 3 . Finally, the results confirm, once more, the lower quality obtained with the non-linear approach compared to the ILE and CIL approaches, irrespective of the measure used.
In Fig. 22, the effect of the stretching factor on the different quality measures is illustrated for the ILE isotropic approach using the mesh with p = 4 and by introducing five load increments. Almost identical results are obtained with the CIL neo-Hookean approach whereas the non-linear approach fails to converge in the majority of the cases.
The results show that the quality measure Q 1 (and Q 2 ) are less influenced by an increase in the stretching factor, compared to the minimum scaled Jacobian Q 3 . In all cases, and for all values of the Poisson's ratio, the value of Q 1 (and Q 2 ) is approximately 0.9, whereas the quality Q 3 can vary from 0.4 to 0.9 depending on the value of the Poisson's ratio and the level of stretching. When the quality Q 3 is considered, the optimal value of the Poisson's ratio is clearly dependent on the level of stretching. For low to moderate stretching factors, a Poisson's ratio near the incompressible limit provides the highest quality whereas for very high stretching factors it is better to consider values in between 0.3 and 0.4.
Finally, Fig. 23 shows the three quality measures Q 1 , Q 2 and Q 3 as a function of the Poisson's ratio for the mesh with p = 6, a stretching factor of 200 and using five load increments. The results correspond to the ILE isotropic and CIL neo-Hookean approaches because the non-linear approach fails to converge in all cases due to the high stretching and high-order considered in this example. This example, shows a different behaviour of the ILE isotropic and CIL neo-Hookean approaches. The CIL approach shows a significant deterioration of the quality measure Q 1 (and Q 2 ) near the incompressible limit, whereas the ILE isotropic approach maintains a high quality for all values of the Poisson's ratio.
Next, under the same setting as in the previous problem, the computational time is analysed, for different formulations and using different material models. Figure 24 shows the CPU time using the three formulations and different material models when the boundary displacement is imposed using five load increments. The Poisson's ratio is ν = 0.4 and the order of approximation is p = 2. These values are deliberately chosen such that the non-linear analysis converges for most material models.
Compared to the previous problem, there is a significant increase in the degrees of freedom and hence the overhead of function calls is insignificant compared to the actual cost of computation. For highly stretched meshes, the Newton-Raphson scheme loses quadratic convergence. The increased number of iterations required and the higher cost of each iteration, due to ill-conditioning, makes the cost of the non-linear approach significantly higher. The ILE (isotropic and TI) approaches are found to be the most competitive. This allows to conclude that, as in the previous example, the ILE approaches provide both the best quality and the lowest computational cost compared to other approacher and material models.
The last study for this example, involves a p-convergence analysis in order to illustrate the optimal approximation properties of the produced meshes. Given a smooth function defined in Cartesian coordinates, the strategy consists on computing the exact value of the solution at the mesh nodes. Then, the error between the approximated solution, interpolated from the nodal values, and the exact solution is computed at each integration point to compute the error in the L 2 (Ω) norm. Figure 25 shows the approximation error in the L 2 (Ω) norm as a function of the square root of the number of degrees of freedom for two different levels of stretching and for a degree of approximation ranging from p = 2 up to p = 9.
The results show the expected exponential convergence in the approximation of a smooth function. In addition, it is interesting to observe that the error is almost identical for the ILE and CIL approaches. This conclusion is in line with the previous analysis where it was shown that the quality of the meshes produced with the ILE and CIL approaches is almost identical, except in some extreme cases considering In contrast, the CIL TI approach, which was shown to produce lower quality for high-order approximation shows a deterioration in the convergence rate, which illustrates the importance of producing high quality meshes for finite element analysis. Finally, the results also show the ability to preserve the approximation properties independently on the level of stretching.

Mesh around the NASA almond
The next example considers a tetrahedral mesh around the NASA almond, a popular geometry for benchmarking 3D radar cross section computations in computational electromagnetics [22,76]. Figure 26 shows the linear surface mesh of the almond, the high-order surface mesh corresponding to a degree of approximation p = 4 and a cut of the highorder volume mesh. The linear mesh contains 6247 elements, 1288 nodes and 688 faces on the almond. The corresponding high-order mesh with p = 6 contains 233,205 nodes and 16,420 nodes to be projected over the true almond geometry to obtain the Dirichlet boundary conditions of the solid mechanics problem.
Similar to the previous examples, the effect of the Poisson's ratio on the quality of the generated meshes is investigated first, for different degrees of approximation. Figure 27 shows the quality measure Q 1 for the linear, incrementally linear and non-linear approaches. In all cases the imposed displacement on the boundary has been introduced using 10 load increments.
Compared to the two-dimensional results of the isotropic meshes in Sect. 6.1, similar conclusions are derived here. First, the quality of both the meshes produced with the ILE isotropic and CIL neo-Hookean approaches is similar, although the ILE isotropic provides better quality near the incompressible limit and, for some particular choices of the approximation degree, for the whole range of values of the Poisson's ratio (e.g., for p = 5). As shown in previous examples, the non-linear approach produces good quality meshes for low-order approximations (i.e., p = 2,3). For p = 4 a valid mesh is only obtained for values of the Poisson's ratio between 0.1 and 0.4, and no convergence is obtained if the order of approximation is further increased.
Similar conclusions are obtained if other quality measures are utilised. For instance, Figs. 28 and 29 show the same analysis in terms of the quality measures Q 2 and Q 3 respectively.
Although the actual value of the quality is different, depending on the selected measure, the qualitative behaviour is the same compared to the quality Q 1 . As reported earlier with the two dimensional examples, the quality measure that produces a lower absolute value is the scaled Jacobian, Q 3 , traditionally used by the high-order mesh generation community. This is attributed to the motion resulting from an imposed boundary displacement that results from projecting the high-order nodes to the true CAD surface. In this scenario, the volumetric deformation related to Q 3 , is much more important than the deformations related to Q 1 and Q 2 . Figure 30 shows the three quality measures Q 1 , Q 2 and Q 3 as a function of the Poisson's ratio for the mesh with p = 3 and using 10 load increments.
The results confirm that, contrary to two dimensional plane strain problems, the quality measures Q 1 and Q 2 are different. It can be observed that the ILE and CIL approaches produce meshes of similar quality, irrespective of the measure considered. In addition, the results illustrate that the quality measure Q 1 is less influenced by the Poisson's ratio whereas the quality Q 3 shows a major dependence on this material parameter. Finally, the results shows that for loworder approximations the non-linear approach can produce meshes of slightly better scaled Jacobian compared to the ILE and CIL approaches although when the quality measures Q 1 and Q 2 are used, the non-linear approach produce the lowest quality meshes compared to the ILE and CIL approaches. This is again due to the non-proportional movement of the nodes in the non-linear approach, which results in distortion of edges and faces of the element, despite a reasonable volumetric deformation being maintained. If a higher order of approximation is considered, say p = 5, the non-linear approach fails to converge for any value of the Poisson's ratio, as illustrated in Fig. 27. A comparison of the different quality measures for the ILE and CIL approaches is shown in Fig. 31.
The results reveal important differences between the ILE and CIL approaches and illustrate the robustness of the ILE approach as the quality is significantly less dependent on the value of the Poisson's ratio selected, compared to the CIL approach. In fact, the results show that high quality meshes can be obtained for the ILE approach with any value of the Poisson's ratio, even with a value near 0, whereas a substan-tial decrease in the quality is observed if a Poisson's ratio near 0 is selected for the CIL approach.
Next, the computational time is analysed. Figure 32 shows the CPU time using the three formulations and different material models when the boundary displacement is imposed using five load increments.
As done in previous examples, the CPU time has been normalised with respect to that of classical linear elasticity and the geometrical mean of 100 run-times, excluding the timing for the first 10 runs, is reported. Compared to previous twodimensional examples, the number of degrees of freedom is now significantly larger for a single core and, therefore, the cost of actual computation dominates over secondary effects such as inlining and branch prediction. The systems of linear equations are now solved using the Multi-frontal Massively Parallel Solver (MUMPS). It is interesting to observe that, despite these differences compared to the two-dimensional examples, similar conclusions are obtained from the CPU time analysis. Once more, both the ILE approaches are found to be the most competitive and the non-linear approaches the most computationally expensive. These results, together with the quality study presented in this section, enables to conclude that the ILE and CIL approaches are recommended for producing high-order curvilinear meshes from an initial linear mesh.
To conclude, p-convergence analysis of the interpolation error is performed to illustrate the optimal approximation properties of the produced meshes. Following the strategy presented in Sect. 6.2, Fig. 33 shows the approximation error in the L 2 (Ω) norm as a function of the cubic root of the number of degrees of freedom for a degree of approximation ranging from p = 2 up to p = 6.
The results show the expected exponential convergence in the approximation of a smooth function. In addition, it is interesting to observe that the error is almost identical for the ILE and CIL approaches. Once more, the CIL TI approach shows a slight deterioration in the rate of convergence for high-order approximations due to the lower quality of the meshes produced with this approach. This result in fact pinpoints the importance of choosing a well-defined polyconvex material model, in the context of a posteriori mesh generation.

Meshes around full aircraft configurations
The next examples consider meshes around two full aircraft configurations, showing the capability of the proposed unified approach for generating meshes around realistic First, a tetrahedral mesh around a generic Falcon aircraft is considered. The linear mesh has 185,191 elements, 35,875 vertices and 16,922 triangular faces on the aircraft to be projected on the true CAD geometry to obtain the Dirichlet boundary condition for the solid mechanics problem. The corresponding CAD geometry has 54 surfaces with 240 intersection curves. For an interpolation degree of p = 3, there  Figure 34 shows the linear surface mesh of the aircraft, the high-order surface mesh corresponding to a degree of approximation p = 3 and a cut of the high-order volume mesh. The problem is solved using the CIL Mooney-Rivlin approach with ν = 0.45 and 20 load increments. The minimum Scaled Jacobian for this mesh is Q 3 = 0.337 and there are 181,251 elements (i.e. 97.87 percent of the total number of elements) for which Q 3 > 0.9. The minimum values of the other two quality measures, accounting for fibre and surface deformations, are Q 1 = 0.605 and Q 2 = 0.467.
Next, a tetrahedral mesh around the DLR-F6 transport configuration is considered. The linear mesh has 68,571 elements, 31,080 vertices and 31,836 tetrahedral faces on the aircraft to be projected on the true CAD geometry to obtain the Dirichlet boundary conditions for the solid mechanics problem. The corresponding CAD geometry has 128 surfaces with 634 intersection curves. For an interpolation degree of p = 4, there are 1,601,015 nodes on the domain and 255,584 nodes on the aircraft that require projection. Figure 35 shows the linear surface mesh of the aircraft, the high-order surface mesh corresponding to a degree of approximation p = 4 and a cut of the high-order volume mesh.
The problem is solved using the ILE isotropic approach with ν = 0.45 and 100 load increments. The minimum values of the three quality measures for this mesh are Q 1 = 0.482, Q 2 = 0.377 and Q 3 = 0.329. Moreover there are only 11 elements with a quality Q 3 < 0.9.
Finally, a boundary layer tetrahedral mesh around the DLR-F6 transport configuration with a stretching of 317 is considered. The boundary layer has been constructed such that the final mesh is suitable for a compressible Navier-Stokes simulation up to a Reynolds number of approximately Re = 4 × 10 7 . The linear mesh has 4,482,662 elements, 787,712 vertices and 110,458 triangular faces on the aircraft. Two curved boundary layer meshes are generated for this geometry with p = 3 and p = 4, respectively. The resulting high-order mesh with p = 3 has 20,434,689 nodes with 498,590 nodes on the aircraft and the p = 4 mesh has 48,279,087 nodes with 885,712 nodes on the aircraft. Both meshes are produced using the ILE isotropic approach. The number of increments are chosen such that a balance is kept between computational cost and final quality of the computational mesh. This corresponds to 60 and 30 increments with a minimum scaled Jacobian of 0.06 and 0.02 for p = 3 and p = 4, respectively. However, for both meshes, 99.5 % of the elements have scaled Jacobian above 0.8. Figure 36 shows the surface mesh of the aircraft for p = 4 and cuts of the high-order volume mesh for p = 3 and p = 4.

Meshes of complex mechanical components
Two complex three-dimensional mechanical components are considered in this section. The Poisson's ratio for all the examples considered in this section is chosen as ν = 0.45.
The first example considers a mechanical valve where the CAD geometry has 45 surfaces and 260 intersection curves. The linear mesh has 16,509 elements 4176 nodes and 5364 triangular faces on the boundary. The resulting high order mesh with p = 5 has 377,994 nodes and 67,047 nodes on the CAD surfaces. The problem is solved using the ILE isotropic approach with 5 load increments and the resulting minimum quality measures are Q 1 = 0.917, Q 2 = 0.841 and Q 3 = 0.768. Moreover, there are only 4 elements for which Q 3 < 0.9. Figure 37 shows two views of the gen-erated high-order curved surface mesh corresponding to a degree of approximation p = 5.
The last example considers a more complex mechanical component and it has been selected to illustrate the robustness and potential of the proposed approach when dealing with complex geometries formed by a large number of surfaces. The corresponding CAD geometry has 638 surfaces with 3459 intersection curves. The linear mesh has 64,599 elements 17,025 vertices and 23,506 triangular faces on the CAD boundary. The resulting high order mesh with p = 4 has 784,670 nodes and 187,903 nodes on the boundary. The problem is solved using the ILE isotropic approach with 200 load increments. The minimum values of the three quality measures for this mesh are Q 1 = 0.719, Q 2 = 0.605 and Q 3 = 0.451, respectively. Moreover, there are only 6 elements for which Q 3 < 0.9. Figure 38 shows different profiles of the high-order curved surface mesh corresponding to a degree of approximation p = 4. Observe that in this mesh, there are 221 planar surfaces and the nodes lying on these surfaces require in-plane translations, which the proposed unified approach is capable of resolving.

Conclusions
A unified approach for the generation of high-order curvilinear meshes derived via a solid mechanics analogy has been presented. This proposed theoretical and computational approach encompasses the incremental linear elastic approach (wherein only the geometry is updated incrementally) and the fully non-linear approach, both previously applied in the context of a posteriori high-order mesh generation. In addition, the new incrementally linearised elasticity formulation (wherein the geometry, the tangent operator and the stresses are updated incrementally), not previously applied to generate curvilinear high-order meshes, is included within this unified approach. The material parameters are calibrated such that the tangent operators of all the aforementioned approaches with various material models are identical in the reference configuration, i.e. for the (undeformed) mesh with planar faces or edges. The derivation of all the approaches, based on energy principles, is used to propose mesh quality measures based on independent invariants of the strain energy density. The relation of the proposed qual- Several numerical examples are presented in both two and three dimensions, including realistic geometries of interest to the solids, fluids and electromagnetics communities. A detailed comparison of all the methodologies is made, including the quality of the generated high-order meshes, the influence of material parameters and load increments on the resulting meshes, the computational cost and the approximation properties of the meshes when applied to an isoparametric finite element formulation.
In terms of the material parameters, the use of a Poisson's ratio near the incompressible limit is generally advised in order to maximise the quality of the resulting mesh. For isotropic meshes, a low number of increments (e.g. five increments) is typically sufficient to obtain the maximum possible quality, whereas for highly stretched meshes and for highorders of approximation (i.e. p > 4) a higher number (e.g. 40 increments) is needed to obtain good quality meshes. Both factors are in fact related as the results show that a higher number of increments is needed when the Poisson's ratio approaches the incompressible limit.
When the material parameters are kept the same, all the linearised approaches, in particular, the incremental linear elastic and the consistent incrementally linearised approach produce meshes of very similar quality and only small differences are observed for highly stretched meshes when high-orders of approximation are used and the Poisson's ratio approaches the incompressible limit. In contrast, the non-linear approach has been found to produce poor quality elements when a high-order approximation is utilised. The non-proportional displacement of interior nodes with respect to the imposed displacement of boundary nodes has a significant negative impact on the convergence of the non-linear solver. Only for low-order approximations has the non-linear approach shown robustness and the ability to produce good quality meshes. The importance of having a well-defined internal energy for the non-linear material model has been illustrated using the transversely isotropic hyperelastic material. For highly stretched meshes, buckling can be expected in the non-linear analysis and the Dirichlet-driven nature of the problems demands a sophisticated and expensive arc-length technique to guarantee convergence, hindering its practical use in an a posteriori mesh generation framework.
The three quality measures proposed for isotropic materials, namely, Q 1 related to fibre maps, Q 2 related to surface maps and Q 3 related to volume maps, show a similar trend with respect to the material parameters. In fact, the first two quality measures are identical for two dimensional plane strain problems. For all the examples considered, Q 3 is the most impactful indicator, which corresponds to the so-called scaled Jacobian traditionally used by the high-order mesh generation community.
In terms of the computational cost, the non-linear approach is much more expensive than the linearised approaches. For highly stretched meshes, where the Newton-Raphson scheme may lose its quadratic convergence due to illconditioning of the system, a higher number of iterations is required and the solver time is drastically increased. The linearised approaches are not only much more economical but, in addition, more robust and produce better quality meshes.
Finally, the approximation properties of the resulting meshes have been assessed and the results show that a similar quality of mesh (as indicated by Q 1 , Q 2 and Q 3 ) translates in similar interpolation errors (i.e. the quality indicators have been shown to be well chosen).