Robust untangling of curvilinear meshes
Thomas Toulorge1 , Christophe Geuzaine2 , Jean-François Remacle1 andJonathan Lambrechts1,3
1

Université catholique de Louvain, Institute of Mechanics, Materials and Civil
Engineering (iMMC), Bâtiment Euler, Avenue Georges Lemaître 4, 1348
Louvain-la-Neuve, Belgium
2
Université de Liège, Department of Electrical Engineering and Computer Science,Montefiore Institute B28, Grande Traverse 10, 4000 Liège, Belgium
3
Fonds National de la Recherche Scientifique, rue d’Egmond, Bruxelles.

Abstract
This paper presents a technique that allows to untangle high order/curvilinear
meshes. The technique makes use of unconstrained optimization where element Jacobians are constrained to lie in a prescribed range through moving
log-barriers. The untangling procedure starts from a possibly invalid curvilinear mesh and moves mesh vertices with the objective of producing elementsthat all have bounded Jacobians. Bounds on Jacobians are computed using
results of papers [1, 2]. The technique is applicable to any kind of elements,
both for surface, volume, hybrid or boundary layer meshes. A series of examples demonstrate both the robustness and the efficiency of the technique.
The final example, involving a time explicit computation, shows that it ispossible to control the stable time step of the computation for curvilinear
meshes through an alternative element deformation measure.
Keywords: High order methods, finite elements, mesh generation

1. Introduction
There is a growing consensus in the computational mechanics communitythat state of the art solver technology requires, and will continue to require
too extensive computational resources to provide the necessary resolution fora broad range of demanding applications, even at the rate that computational
Email address: thomas.toulorge@uclouvain.be (Thomas Toulorge1 , Christophe
Geuzaine2 , Jean-François Remacle1 and Jonathan Lambrechts1,3 )Preprint submitted to Journal Of Computational Physics

February 26, 2013

power increases. The requirement for high resolution naturally leads us toconsider methods which have a higher order of grid convergence than theclassical (formal) 2nd order provided by most industrial grade codes. This
indicates that higher-order discretization methods will replace at some point
the current finite volume and finite element solvers, at least for part of theirapplications.
The development of high-order numerical technologies for engineering
analysis has been underway for many years now. For example, Discontinuous
Galerkin methods (DGM) have been largely studied in the literature, initiallyin a theoretical context [3], and now from the application point of view[4]. In many contributions, it is shown that the accuracy of the methodstrongly depends on the accuracy of the geometrical discretization [5–7].
Consequently, it is necessary to address the problem of generating the highorder meshes that are needed to fully benefit from high-order methods.
Modern mesh generation procedures take as input CAD1 models composed of model entities. Each model entity Gdi has a geometry (or shape) [8,9], that depends on the solid modeler for its underlying representation. Solidmodelers usually provide a parametrization of the shapes, that is, a mapping ξ ∈ Rd 7→ x ∈ R3 . The geometry of a model vertex G0i is simply
its 3D location, xi = (xi , yi , zi ). The geometry of a model edge G1i is its
underlying curve with its parametrization x(t), t ∈ [t1 , t2 ]. The geometry ofa model face G2i is its underlying surface with its parametrization x(u, v),(u, v) ∈ S ⊂ R2 . The geometry associated with a model region G3i is R3 .
There are also four kind of mesh entities: mesh vertices Mi0 , mesh edges Mi1 ,mesh faces Mi2 and mesh regions Mi3 that are said to be classified on model
entities2 . Each mesh entity is classified on the model entity of the smallest
dimension that contains it.
The way of building a high order mesh is to first generate a straight sided
mesh. Then, mesh entities that are classified on the curved boundaries ofthe domain are curved accordingly (see Figure 1). For mesh edges that areclassified on model edges (for example M21 @ G11 on Figure 1), additional
high order mesh points are added on the geometry of the model edge. Then,high order points are added on mesh edges that are classified on model faces(for example M11 @ G21 on Figure 1). Finally, high order points may be
added on mesh faces that are classified on model faces. The position of thehigh order points can be chosen in such a way that the geometrical error,
1
2Computer Aided Design.
We use the symbol @ for indicating that a mesh entity is classified on a model entity.

2

Mesh Curving

Mesh edge M11 @ G21

Model face G21
Model edge G11

Mesh edge M21 @ G11

Figure 1: Straight sided mesh (left) and curvilinear (cubic) mesh (right).

i.e. the distance between the CAD model and the mesh, is minimized. In
this paper however, the high order points are simply defined by orthogonal
projection from the straight mesh entity onto the curved geometry definedby the CAD model.
The naive curving procedure described above does not ensure that all theelements of the final curved mesh are valid. Figure 2 gives an illustrationof this important issue. Some of the curved triangles are tangled: they selfintersect after having been curved. It is important to note that this problem
is not related to the accuracy of the geometrical discretization: in Figure 2,
the mesh would not be valid even if the curved edge was assigned the exact geometry (blue curve). Invalid elements may be detected by exploiting

First-order mesh

High-order vertices

Untangling

Figure 2: Straight sided mesh (left) basic curvilinear (quadratic) mesh with tangled elements (center) and untangled mesh (right).

specific properties of the Jacobian of the transformations that map a fixedreference element onto each element in the physical domain, as shown in
Figure 3. In particular, a change of sign in the Jacobian over an element
indicates that the corresponding map is not injective, so the numerical meth3

x3
X2

x5

x(X)

X5
X3

Y

x(ξ)

X(ξ)

x4

ξ3

X6
X

x2

x6

η

X4

y

ξ5

ξ6

x1
x

X1

ξ
ξ1

ξ4

ξ2

Figure 3: Reference unit triangle in local coordinates ξ = (ξ, η) and the mappings x(ξ),X(ξ) and X(x).

ods that rely on integration over elements cannot be used. This is the casefor self-intersecting elements. In two recent papers [1, 2], a general formulation has been developed for computing robust estimates of the geometricalvalidity of a curvilinear element. Provable bounds on element Jacobians
can be computed for high order triangles, quads, tetrahedra, hexahedra andprisms.
Figure 2 indicates that, without refining the mesh, the only way of generating a valid high order mesh is to curve not only mesh entities classified
on curved model entities, but also those that are initially straight-sided inside the computational domain. It is thus necessary to somehow propagate
the curvature inside the domain. Several smoothing schemes have been proposed in the literature to this effect: linear smoothing techniques such as
Laplacian smoothing [10], Winslow smoothing [11] or linear elasticity withvarying stiffness [10, 12]. Even though such simple techniques may often leadto interesting results, there is no guarantee whatsoever that applying such a
linear smoother will result in an untangled mesh. Nonlinear smoothing techniques have also been investigated, in particular using a nonlinear elasticity
analogy [13], which results in a valid curvilinear mesh. Yet, computing a nonlinear mechanics problem including large deformations on a high-order andhighly stretched mesh is numerically very complex (more than solving theNavier-Stokes equations on the same grid, for instance). The computational
cost may then be excessive, given that the meshing time is expected to be4

a fraction of the overall CPU time. Besides linear and nonlinear smoothing
techniques, other authors [14–16] make use of mesh adaptation techniquesby eliminating invalid elements by a combination of local mesh refinements,
edge and face swaps, and node relocations.
In this paper, we propose a robust smoothing scheme that makes it possible to build a curvilinear mesh for which the validity of every element
is guaranteed. This new untangling procedure relies on an optimization
algorithm rather than a mechanical analogy or local mesh refinement: it
specifically targets element Jacobians and modifies node locations in such a
way that Jacobian values lie in a predefined range. This approach is loosely
related to the optimization schemes proposed by several authors for untangling straight-sided meshes [17–20], but with a focus on arbitrarily high-orderelements.
The paper is organized as follows. In Section 2, we briefly recall the
results of paper [1, 2] on Jacobian bounds. Section 3 is dedicated to thepractical computation of both Jacobian bounds and their derivatives withrespect to the motion of mesh vertices. Section 4 is the kernel of the paper.
An objective function that specifically targets invalid Jacobians is describedin Section 4.1. Constraints on Jacobian positivity are imposed through logbarriers, allowing the use of unconstrained optimization procedures. The
optimization strategy is described in Section 4.2. The optimization starts
with an invalid mesh and the asymptote in the log barrier is progressively
moved into the valid region. Some examples of mesh untangling are presented in Section 5 where both the robustness and the efficiency of the new
methodology are demonstrated. Some 3D examples are presented in Section 5.1 with timings and statistics. In Section 5.2, the example of a highorder boundary layer mesh is presented. A new term is added to the objective function with the aim of controlling the maximum Jacobian as well.
Finally, Section 6 presents a simple idea for controlling the impact of meshcurvature on the explicit time step of computations.
2. Validity estimates of curvilinear meshes
Let us introduce the following notations. We call ne and nv the numberof elements and vertices of the mesh. Each element e of the mesh contains
Np vertices. Note that in this context, the word “vertex” denotes any nodeof the mesh, whether it is a geometrical vertex of a straight-sided element or
a high-order node of a curved element.
We note X ei the position of the ith vertex of the element e in the straightsided configuration and xei the location of the same vertex, yet in some5

deformed configuration.
The shape of an element e is defined geometrically through its vertices xei ,i = 1 . . . Np that are nodes associated to a set of Lagrange shape functions(p)
Li (ξ), i = 1 . . . Np at order p. This allows us to map a reference element
to the real one:NpX
(p)x(ξ) =Li (ξ) xei .
(1)i=1Consider now the transformation x(X) that maps straight sided elements
onto curvilinear elements (see Figure 3). The mapping x(X) should be bijective, i.e. it should admit an inverse. This implies that the determinantof the Jacobian det x,X has to be non-zero, for every value of ξ in the reference element. In practice, this condition reduces to strict positivity for 3D
elements, but an ambiguity remains for 2D elements, in which the Jacobian
is defined with respect to a local normal. The orientation of planar 2D elements can be chosen in such a way that the direction of the constant normal
ensures the positivity of the Jacobian. For curved surfaces in 3D however,the definition of the Jacobian is more complex, so that the condition of strict
positivity that we retain in this paper may be overly restrictive.
It is possible to write this determinant in terms of the ξ coordinates as:
det x,X =det x,ξ
J(ξ)= e ,det X ,ξJ0

where J0e is the strictly positive and constant3 Jacobian of the straight sided
element. The function x(X) is called the distortion mapping. Its determinant det x,X , that we call the scaled Jacobian, should be as close to 1 aspossible in order not to degrade the quality of the straight sided element e.
In [1, 2] it is shown that it is possible to reliably detect invalid elements.
In other words, it is possible to find reliable bounds to Jmin = minξ J andto Jmax = maxξ J over the whole element. The Jacobian J is a polynomialin ξ. It can then be interpolated exactly as a linear combination of Bézier
(q)
polynomials Bi at a certain order q ≥ p over the element. Provable bounds
for Jmin and Jmax are then computed using an interesting property of the
Bézier polynomials, namely that their value is contained within the range of3

Straight sided element Jacobians are constant for simplical elements only, i.e. triangles
in 2D and tetrahedra in 3D.

6their coefficients. Assuming that J is a polynomial of order q in ξ, we writeJ(ξ) =

Nq
X
i=1

(q)

Bi (ξ)Bi(2)

and bounds can be computed as
min J(ξ) ≥ min Bi and max J(ξ) ≤ max Bi .
ξ

i

ξ

i

Although they are reliable, the bounds obtained in this manner may not
be accurate: mini Bi may be much lower than minξ J(ξ) and maxi Bi may
be much greater than maxξ J(ξ). The adaptive procedure described in [1, 2]
overcomes this problem, but it implies a higher computational cost, as wellas additional complexity in the formulation. Therefore, we favor the direct
evaluation provided by Eq. (2) in the present work, even though it yields anoverly restrictive estimation of the Jacobian bounds.
The following section is dedicated to the practical computation of the coefficients Bi as well as their derivatives with respect to the vertex coordinatesxei .
3. Computation of Bézier coefficients and their derivatives
The aim of our method is to untangle both surfacic and volume meshes.
For that, we assume that a point x has always 3 coordinates x = {x, y, z}.
Local coordinates ξ = {ξ, η, ζ} are also supposed to be three dimensional.
Yet, for surface meshes, we assume that vector∂x ∂y ∂zn=, ,∂ζ ∂ζ ∂ζis the constant unit normal vector to the straight sided element. With that
hypothesis, it is possible to compute the determinant of the Jacobian at
every Lagrange node ξk = (ξk , ηk , ζk ) at order q:
Jk = J(ξk ) =∂x ∂y ∂z ∂z ∂x ∂y ∂y ∂z ∂x++−∂ξ ∂η ∂ζ∂ξ ∂η ∂ζ∂ξ ∂η ∂ζ∂z ∂y ∂x ∂x ∂z ∂y ∂y ∂x ∂z−−.
∂ξ ∂η ∂ζ∂ξ ∂η ∂ζ∂ξ ∂η ∂ζConsidering that
x=

Np
X
i=1

(p)

xei Li (ξk ),7(3)

it is possible to compute the sensitivity of the Jacobian at point k withrespect to the x coordinate of the vertex i:
∂Jk
∂xei

(p)

(p)

=

(p)

∂Li ∂y ∂z ∂z ∂Li ∂y ∂y ∂z ∂Li++∂ξ ∂η ∂ζ∂ξ ∂η ∂ζ∂ξ ∂η ∂ζ
(p)

(p)

(p)

∂z ∂y ∂Li∂ξ ∂η ∂ζ−−∂Li ∂z ∂y ∂y ∂Li ∂z−.
∂ξ ∂η ∂ζ∂ξ ∂η ∂ζ

(4)

∂Jk
k
The same computation can be done for ∂J∂yie and ∂zie . In practice, the following matrix J of size Nq × (3Np + 1) is computed for every element e:
 ∂J

∂J1
∂J1
∂J1
∂J1
∂J1
1
. . . ∂x. . . ∂y. . . ∂z
J1
e
e
e
e
e
∂xe1
∂y
∂z
1
1
Np
Np
Np


..
..
..
..
..
.. 
 ..
J= .
.
.
.
.
.
. 
 ∂JN

∂JNq
∂JNq
∂JNq
∂JNq
∂JNq
q
.
.
.
.
.
.
.
.
.
J
e
e
e
e
e
e
Nq
∂x
∂x
∂y
∂y
∂z
∂z
1

Np

1

Np

1

Np

(q)

q
= Bl (ξk ) is the transformation matrix that enables toAssuming that Tlk
compute Bézier coefficients Bl using Lagrange coefficients Jk , the matrix
 ∂B

∂B1
∂B1
∂B1
∂B1
∂B1
1
.
.
.
.
.
.
.
.
.
B
e
e
e
e
e
e
1
∂x1
∂xNp
∂y1
∂yNp
∂z1
∂zNp


..
..
..
..
..
.. 
 ..
B= .
.
.
.
.
.
. .
 ∂BN

∂B
∂B
∂B
∂B
∂B
N
N
N
N
N
q
q
q
q
q
q
.
.
.
.
.
.
.
.
.
B
e
e
e
e
e
e
N
q
∂x
∂x
∂y
∂y
∂z
∂z
1

Np

1

Np

1

Np

(5)
that contains both the Bézier coefficients Bl as well as their gradients withrespect to the position of nodes of element e is calculated through a singleqmatrix-matrix product: Blj = Tlk
Jkj .
It is then possible to use the Bi ’s and their gradients in a gradient-based
optimization procedure. In what follows, an objective function that contains
the coefficients Bi is constructed in order to control the quality of elements.
4. Curvilinear mesh untangling
4.1. An objective function
This section describes the objective function f (xei ) that is used to untangle invalid curved elements through an unconstrained optimization procedure. We design a function
f =E +F
8

that is composed of two parts E and F.
Our assumption is that the method is provided with a straight-sided meshof high quality. This mesh has potentially been defined to satisfy multiple
criteria, such as a predetermined size field, or anisotropic adaptation. Whencurving such meshes, it is desirable to preserve as much as possible all thesefeatures, which implies keeping the nodes as close as possible to their initial
positions in the straight sided mesh.
To this end, we want to introduce some kind of energy E associated withthe displacement of the nodes xei − X ei , i.e. a positive quadratic form thatis a measure of the distance between the straight sided nodes X ei and theirposition xei in the curved mesh:
Np Np

1 XXX eE(xi , K) =(xi − X ei )Keij (xej − X ej ) ≥ 0
2 e

(6)

i=1 j=1

where K is a symmetric positive matrix of size 3nv × 3nv and Keij is of size
3 × 3. In this paper, we choose K as the diagonal matrix with entries wie /L2 ,where wie are user-defined weights used to set the balance between the twoparts E and F of the objective function f , and L is a length scale representative of the problem. We choose L as the maximum distance, among allvertices of the initial tangled mesh, between a node and its counterpart in thestraight-sided mesh. As for the non-dimensional weights wie , we usually set
wie = 105 if the node i of element e lies on the boundary, and wie = 100 otherwise. However, we observe in practice that the value of wie has little influenceon the convergence of the optimization method. The presence of E prevents
the problem from being under-determined, and it orients the optimization
procedure towards a solution that tends to preserve the straight-sided mesh,but the term F dominates for the most problematic (tangled) elements thatdrive the mesh deformation.
The second part F of the functional deals with Jacobian positivity. We
use a log barrier [21] in order to avoid Jacobians that are too small, and aquadratic function to penalize Jacobians that are too high:F(xi , ) =with
Fle (xei , ) =




log

Nq
ne X
X

Fle (xei , )

e=1 l=1

Ble (xei ) − J0e
J0e − J0e
9

2


+

Ble (xei )
−1
J0e

2
,

(7)

Be

F ( J el )015F (x) = (x − 1)2 + log



x−
1−

2

10

5

0

0 

1

2

3

Ble
J0e

Figure 4: Plot of the barrier function F (Jle ) for  = 0.2.

that is defined in such a way that F blows up when Ble = J0e , but still
vanishes whenever Ble = J0e . Barrier methods are among the most powerfulclass of algorithms available for attacking general nonlinear optimization
problems. These techniques converge to at least a local minimum in mostcases, regardless of the convexity characteristics of the objective function andconstraints [22]. Figure 4 shows a plot of the barrier function in Equation (7)for  = 0.2.
The objective function f being smooth, it is possible to compute its
gradient ∇f with respect to the positions of the element vertices xei . This
gradient is used in a gradient-based optimization procedure.
Mesh vertices can be classified on mesh entities of various dimensions.
Some of the mesh vertices Mi0 @ G1j are classified on model edges. Sucha vertex can only be moved along G1j , i.e. its position only depends on one
curve parameter t. We have thereforedf∂f dxei=·dt∂xei dtdxewith dti the tangent vector to the curve at point t, that is obtained fromthe CAD model.
Other vertices, that are classified on model faces Mi0 @ G2j , can only be
moved along the surface. Two parameters u and v are associated to thosevertices. We have therefore∂f∂f ∂xei∂f∂f ∂xei
=
·
and
=
·
∂u
∂xei ∂u∂v∂xei ∂v
10

∂xe

∂xe

with ∂ui and ∂vi the two tangent vectors to the surface at point (u, v).
Those can be computed using the CAD model.
Vertices that are classified on model regions have a complete freedom tomove in every direction of the 3D space, in which case the physical coordinates {x, y, z} are used. Finally, mesh vertices that are classified on model
vertices have no freedom to move, and are excluded from the optimization
problem.
4.2. Optimization Strategy
The problem of untangling curvilinear meshes is defined as
min f (xi , K, ), i = 1, . . . , nv .
xi

There is a variety of methods that can be used to solve unconstrained
minimization problems. We have tested several alternatives: interior point
methods implemented in the software package IPOPT [23], as well as LBFGS [24] and conjugate gradients [25] algorithms provided by ALGLIB [25].
In the end, the use of conjugate gradients seemed to be the best choice interms of computational efficiency.
The most important part of the optimization strategy is indeed to define
the right sequence of optimization problems. Some preliminary remarks canbe made at this point:• The optimization should not be applied to the whole mesh but locally.
Blobs of elements that surround an invalid element are constructed.
Mesh vertices that lie on the boundary of the blob are fixed.
• It is necessary to apply preconditioning to the optimization problem,
because the scale of parametric or physical coordinates of different
mesh vertices can differ by orders of magnitude, depending on themodel entity on which they are classified. We found that a simple
dxe
diagonal preconditioner, based on the norm of the tangent vectors dti∂xe∂xe
for vertices classified on model edges (respectively ∂ui and ∂vi forvertices classified on model faces), and unity for vertices classified onmodel regions, yields fast and robust convergence.
• In order for the variables to remain in the domain of definition of thebarrier function F, the value of the barrier  must be lower than theworst scaled Jacobian of all elements in one given blob. In particular, has to be negative for an initially tangled mesh. Therefore, we compute a sequence of optimization problems with “moving barriers”:  is
11

F0

Iterations

F1

F2
• J min

F

•

60

20
15

•

30

10
5

0 

0

•

−2

−1

1

0 2 0.5

12

0Figure 5: Optimization process: three successive series of (maximum) 30 conjugate gradient iterations, with their respective log barriers.

increased between each optimization problem, until the desired barrier
value ¯ is reached. This procedure is illustrated in Figure 5.
The optimization strategy is described in Algorithm 1.
Algorithm 1: Optimization strategy
1 Compute non overlapping element blobs Bk , k = 1 . . . NB ;2 for k = 1 to NB do
3
4

repeatBe
compute κ = mine minl J el , e ∈ Bi , l ∈ [1, Nq ];
0

5
6
7
8

if κ < ¯ thenset  = κ − 0.1 |κ|;solve minxi f (xi , K, ) for all elements of blob Bk ;Be

recompute κ = mine minl J el , e ∈ Bi , l ∈ [1, Nq ];0910until κ ≥ ¯;As an example, consider a coarse 3D tetrahedral mesh of a sphere, aspresented in Figure 6. The surface of the sphere is described in the CAD
system as one single patch that covers the whole range of spherical coordinates. In order to challenge our optimization strategy, high order nodes
classified on the surface are added along lines in the parameter plane. The
resulting mesh that is presented in the middle image of Figure 6 is clearly12Figure 6: Examples of mesh untangling with mesh vertices motions on manifolds. Upper
figures describe the untangling of a 3D coarse mesh of a sphere. The straight sided mesh(left) is made quadratic (center) and is subsequently untangled (right).

wrong. Our untangling strategy is then successfully applied to the invalid
mesh: the final valid mesh that is presented in Figure 6 has all elements withscaled Jacobians in the range [0.9, 1.1]. Less than one second is required forconverging. The other example in Figure 6 is a quartic surface mesh that is
untangled using the same procedure.
5. Examples of mesh untangling
5.1. 3D Mechanical parts
Figure 7 shows images of the mesh of a mechanical part before and after
the untangling process. The mesh is composed of 5,605 quadratic tetrahedra
e /J e < 0.3. The untangling process was performedof which 215 have a Jmin0
considering 13 blobs of elements. The process for constructing one blob can
be described as follows. One invalid element is initially inserted into theblob. N layers of elements around the invalid element are inserted in the
blob. Vertices that lie on the boundary of the blob are fixed in order not totake any risk of invalidating elements outside the blob in the optimization
process. In this process, blobs may of course overlap. A value of N = 2 isusually sufficient for ensuring the convergence of the optimization process.
Note that this way of building the blobs is not adequate for untangling highly
stretched meshes. This issue will be addressed in the next subsection.
13

Initial quadratic mesh

Invalid tetrahedra

Some blobs of elements

Final untangled mesh

Figure 7: Mesh of a mechanical part that is composed of 5,605 quadratic tetrahedra.
Figures show a surface view of the initial quadratic mesh, a volume view of the 215 invalid
tetrahedra, a view of 10 among 13 blobs of elements used in the optimization process andthe final untangled mesh.

14

Figure 8: Mesh of a mechanical part that is composed of 2,988 quadratic tetrahedra.
Figures show a surface view of the initial quadratic mesh, a volume view of the 131 invalid
tetrahedra and the final untangled mesh.

In the example of Figure 7, we take the value N = 2 and the total CPU
time for untangling that mesh is 18 seconds on a standard laptop. A global
optimization procedure involving one single blob of 5,605 elements takes 60
seconds to converge.
Another example is presented on Figure 8. The same parameters as in
the first example is used. The total time needed to untangle the 3 blobs ofelements is now 23 seconds.
5.2. High-lift airfoil
We consider the high-lift airfoil shown in Figure 12. Two straight-sided
meshes M1 and M2 , featuring boundary layer-type stretching around every
surface of the airfoil, are generated. M1 is made out of 9,814 triangles, while
M2 contains 3,178 triangles and 3,318 quads. Both meshes are converted toorder 2 to 5 and optimized.
The objective function described in Section 4.1 is designed to untangle
invalid elements that are created by the conversion of straight-sided meshes
to higher order, which happens mainly near convex parts of the geometry.
However, we note that the conversion can also generate elements of lower
quality close to concave parts, as seen in Figure 9. These elements are notstrictly invalid, because their Jacobian is clearly positive everywhere, butthey do not maintain the boundary layer-type size progression imposed in thestraight-sided mesh. In order to remedy this inconvenience, we supplement
the procedure described in Section 4 with a second optimization pass that
controls the maximum Jacobian in the mesh. In the first pass, the objectivefunction f introduced in Section 4.1 is used in the “moving barrier” procedure
15

Figure 9: Detail of the mesh on the suction side of the slat in the high-lift airfoil case:primary second-order mesh (left), optimized mesh (right).
Be

F ( J el )015F (x) = (x − 1)2 + log



x−
1−

2

10

5

0

0

1

2

 3

Ble
J0e

Figure 10: Plot of the barrier function F (Jle ) for  = 2.75.

mentioned in Section 4.2 to reach the desired minimum Jacobian ¯min . In
the second pass, we switch to a different objective function f1 :
f1 (xi , K, ¯min , ¯max ) = E(xi , K) + F(xi , ¯min ) + F(xi , ¯max )where the term F(xi , ¯min ) maintains the minimum Jacobian constraint atthe value ¯min , while the term F(xi , ¯max ) (plotted in Figure 10) is handled
as a moving barrier that is decreased to reach the desired maximum Jacobian ¯max . As a result, the elements located near the concave parts of thegeometry also curve to adapt the high-order boundary, as shown in Figure 9.
In boundary layer-type meshes, where elements have high aspect ratio,
the boundary deformation responsible for the mesh tangling is typically large
16Figure 11: Detail of the mesh near the leading edge of the slat in the high-lift airfoil case:tangled second-order mesh (left), blob definition with a circle representing the geometricalcriterion (center), and untangled mesh (right).

compared to the normal size of the tangled element, but small compared to
its tangential size (see Figure 11). In order to untangle it, several layers mustbe curved in the direction normal to the boundary, whereas modifying theelements that are adjacent in the tangential direction is useless. Therefore,
the blobs in which we apply the optimization procedure can be constructedby adding a geometrical criterion to the procedure described in Section 5.1,
as illustrated in Figure 11: among the N layers of elements surrounding theinvalid element, only those located within a certain distance are retained.
This distance is defined by multiplying the distance between the straightsided and high-order boundaries by a user-defined factor. Given the element
aspect ratio and grid stretching in the boundary layer zone for this test case,
we obtain good results with blobs of N = 6 layers of quads (or N = 12 layersof triangles) and a distance factor of 12.
We optimize meshes M1 and M2 with ¯min = 0.4 and ¯max = 1.6 atorder p = 2 . . . 5. Information about the results is provided in Table 1. The
number of operations needed to compute f and its derivatives depends on
the size of matrix B of Equation (5), i.e. it is a function of Np × Nq . In
table 1, we report Np and Nq for the different mesh orders p. It must be
noted that, for mesh M2 , all elements in the boundary layer are quads, sothat Np and Nq correspond to quads. The number of Bézier points Nq for
triangles corresponds to order q = 2(p − 1) while the Jacobians of quads areof order q = 2p − 1. It is therefore not surprising that the computation time
for untangling quartic meshes is about 10 times larger than to untangle aquadratic mesh. Quintic meshes (p = 5) however are disappointingly expensive to compute for quads. This problem may be essentially related to thebad conditioning of the optimization problem. In our implementation, we
use Lagrange equidistant shape functions, which are known to be highly os-

17

Mesh
M1
M1
M1
M1
M2
M2
M2
M2

p
2
3
4
5
2
3
4
5

q
2
4
6
8
3
5
7
9

Np
6
10
15
21
9
16
25
36

Nq
6
15
28
45
16
36
64
100

Ninvalid
51
56
64
74
32
34
31
31

nv
20, 751
45,452
79,967124,296
20,75145,452
79,967124,296t(sec)
3.0
8.2
20.0
65.2
1.4
7.0
36.3
330.7

e /J e
Jmin
0
0.40
0.40
0.40
0.41
0.55
0.50
0.46
0.44

Table 1: Statistics for the high-lift airfoil test case.

cillatory for high orders. The use of polynomial basis with a better Lebesgue
constant such as those of Ref. [26] may improve the results. Nevertheless,our procedure remains quite efficient for high order meshes of orders thatare required for engineering analysis, i.e. p = 2 and p = 3 [27, 28].
Table 1 may look surprising: the target for the minimum scaled Jacobiane /J e is sometimes
has been taken as ¯min = 0.4, but the final value of Jmin0quite higher than the threshold. Indeed, in the optimization procedure,
Jacobians are scaled by a constant value J0e computed on the initial mesh.
Yet, corner vertices of elements may be changed in the optimization process,leading to possible changes in values of straight sided Jacobians J0e . Thus, thee /J e could be lower than the threshold, but the positivity
final value of Jmin0of the Jacobian is not threatened as long as the optimization procedure is
successful.
Finally, Figure 12 shows a comparison between the linear elastic analogy
and the present procedure for the mesh M1 at p = 2. It illustrates thefailure of the simple elastic analogy approach in presence of highly stretchedboundary layer elements [13], whereas the optimization procedure manages
to untangle the mesh by propagating the high-order boundary deformation
through the neighboring layers of elements.
6. Mesh curvature impact on the explicit time step: a simple idea
6.1. Another measure of deformation
It is shown in Ref. [29] that the conditioning of the operators derivingfrom high-order spatial discretization methods can be affected by the curvature of the elements. When using implicit time integration in simulations,18

Figure 12: Quadratic triangular mesh for the high-lift airfoil case. Zoom at the leadingedge of the slat: comparison between the mesh obtained by elastic analogy (upper frame)and by the present optimization procedure (lower frame).

19

(a)

(b)

(c)

Figure 13: Typical optimization of a two-layer patch around an invalid element of theGreat Barrier Reef 3rd order coarse mesh: (a) initial patch and original linear mesh (indashed lines), (b) first optimization based on the Jacobian determinant and (c) second
optimization based on the minimum eigenvalue of the metric matrix M. The blue line isthe coastline given by the CAD model.

this effect has a negative impact on the convergence of iterative linear algebra solvers. With explicit time stepping schemes, it reduces the maximum
time step allowed by the stability restrictions, as explained more thoroughly
in Section 6.2.
The optimization procedure described in Section 4 is very efficient at
producing meshes with scaled Jacobians that are all positive and close to 1.
In mathematical terms, the procedure tends to generate mappings x(X) thatare one-to-one and close to be equiareal (in 2D). As mentioned in Section 4.1,
this is due to the term F that dominates the objective function f for invalid
elements. In comparison, the part E is small, so that the vertices can be
moved enough to fix the tangled elements. In other words, the optimization
procedure tries to untangle high-order elements by recovering the area in 2D(or the volume in 3D) that they have in the initial straight-sided mesh.
In some cases however, the approximate conservation of the area or volume may not prevent an untangled element to end up with a wildly different
shape from the original straight-sided one, as illustrated in Figure 13. This
may affect the numerical behaviour of the simulations using such meshes.
For instance, consider the three quadratic triangles shown in Figure 14, thathave a scaled Jacobian strictly constant and equal to one. It is clear thatthose triangles have very different interpolation properties in a finite element
computation. Prescribing the Jacobian of the elements is thus required toensure the mesh validity, but it is clearly not sufficient for controlling themesh quality. The mesh optimization procedure should also be able to control element lengths, in order to avoid problematic cases like the the one
shown in Figure 13.
20

Figure 14: Three second-order triangles with the same area and the same constant Jacobian determinant but obviously leading to different stable time step.

Consider the two triangles of Figure 15. Let A = x,X be the constant
Jacobian matrix of the transformation between the two triangles. The untransformed triangle and the transformed triangle have the same surface(det A = 1). The first fundamental form or metric tensor of the mapping
M = AA∗ .

(8)

allows to compute variations of lengths. Assume U to be a vector that defines
a direction in the undeformed space (see Figure 15). Vector U transforms
in AU as it is depicted on the figure. The ratio between the lengths of thetwo vectors is written asrkAU k
U ∗ MU
l=
=
.
kU k
U ∗U

A=

2 1
0 1/2




,M=

5 1/2
1/2 1/4



1
U1

AU1

1/2
AU2

U2
1
Untransformed

2
Transformed√

1k
Figure 15: Two triangles that have the same area. We have here kAU= 2√372 ' 2.15
kU1 k
√
√
p
21+ 3772k√
which is close to λmax (M) =' 2.247 and kAU= √15 ' 0.447 which is close
kU2 k
2 2
√
√
p
21− 377√
to λmin (M) =' 0.445.
2 2

21

Then, we havekAU k=
kU k>0 kU klmin = minqλmin (M)where λmin (M) is the smallest eigenvalue of M.
Thus, considering the eigenvalues of the metric tensor in addition to theJacobian gives a better characterization of the mesh quality. This problem
is well posed in the sense that the only mapping that result in a metrictensor that has three equal eigenvalues (e.g. the isometry) does not deformthe element. In this perspective, we extend our approach not only to obtain
positive Jacobians but also to control the minimum eigenvalue of M. We
thus look locally at how element lengths are transformed and try to control
the minimum value λmin on all the element.
In two dimensions, M is expressed as:kx,X k2 x,X · y,XM=.
x,X · y,X ky,X k2
And its smallest eigenvalue, λmin can be obtained easily:qλmin (M) = kx,X k2 + ky,X k2 − (kx,X k2 − ky,X k2 )2 + 4(x,X · y,X )2
For each element blob, we introduce a second stage to the optimization
process described in Section 4. After the optimization of the Jacobian in
the first stage, the same procedure is repeated with a modified objectivefunction:f2 (xi , K, ¯, ζ) = E(xi , K) + F(xi , ¯) + G(xi , ζ)
The term F(xi , ¯) maintains the minimum Jacobian constraint ¯, while theadditional term G controls the minimum eigenvalue of the metric tensor:
G(xi , ζ) =

Nq
ne X
X

Gel (xei , ζ)

e=1 l=1

  min e22
λ (xi ) − ζe eGl (xi , ζ) = log+ λmin (xei ) − 1 .
1−ζUnfortunately, due to the square root, λmin is not polynomial. Therefore, itcannot be expressed as a sum of Bézier polynomials to guarantee that it isabove the objective value everywhere on the element. In practice, the valueof λmin is sampled at the same points as the Jacobian and these additional22

constraints, combined with those on the Jacobian, are enough to preserve
the quality of the elements.
Assuming the usual expansion of Eq. (1) for x, the derivative of λminwith respect to nodal positions can be computed as!
2 − ky k2 )x
(kx
k
+
2(x
·
y
)y
,X
,X
,X
,X
,X ,X
(p)
λmin (M),xe = 2 x,X − p· Li,Xi(kx,X k2 − ky,X k2 )2 + 4(x,X · y,X )2!
(ky,X k2 − kx,X k2 )y,X + 2(x,X · y,X )x,X(p)min
· Li,X .
λ (M),ye = 2 y,X − p
2
2
2
2
i
(kx,X k − ky,X k ) + 4(x,X · y,X )Figure 13 shows how optimization on λmin affects the resulting meshes.
Left Figure shows the initial invalid mesh. Algorithm 1 is then applied,
resulting in the mesh of the central figure. This mesh, while perfectly valid,
has elements that are dramatically shrinked in some directions. Finally,
optimization on λmin is applied, resulting on the mesh of the right Figure.
This mesh has element shapes which are very close to those of the straightsided mesh.
6.2. Impact of the deformation with explicit time stepping
It is mentioned in Section 6.1 that the deformation of high-order elements
may affect the numerical methods used in simulations through the mappingx(ξ). Hereafter, we explain how these effects can be characterized by the
deformation measure introduced in Section 6.1, in the particular case ofnumerical methods using explicit time stepping schemes to resolve unsteadyphenomena.
With explicit time stepping schemes, the curvature of the elements reduces the maximum time step allowed by the stability restrictions throughthe Jacobian matrix x,ξ of the mapping x(ξ), as shown in Ref. [29]. Thisis illustrated in Figure 16, where the evolution of the maximum time step isplotted while progressively deforming a grid, for a simple advection problem
solved with a DG spatial discretization and an explicit Runge-Kutta time
integrator. Thus, setting a scaled Jacobian barrier ¯ closer to 1 tends tomitigate the adverse effects of high-order mesh curvature.
However, in the simple example of Figure 16, the elements are deformed
in a constant direction. The general results of Ref. [29], as well as our experience with high-order computations, suggest that there is no simple relation
between the conditioning of the semi-discrete operator and the Jacobian
det x,ξ .

23

1
0.9
0.8

∆t/∆t0

0.7
0.6
k=1
k=2
k=3
k=4
k=5
k=6
k=7
k=8
k=9
k=10

0.5
0.4
0.3
0.2
0.1
0

0.1

0.2

0.3

0.4

0.5
0.6
Jmin/J0

0.7

0.8

0.9

1

Figure 16: Impact of the progressive deformation of a structured grid for a scalar advection
problem, solved with the classic explicit 4-stage Runge-Kutta time integrator and a DG
spatial discretization of order k ranging from 1 to 10: deformation of the pattern of
triangular elements of which the structured grid is made up (left), evolution of the time
step ∆t (scaled by the time step ∆t0 for straight-sided elements) as a function of thescaled Jacobian Jmin/J0 (right). More details can be found in Ref. [29].

Consider a simple scalar advection equation, as a model for more general
hyperbolic systems: find u(x, t) solution of∂u+ V · ∇x u = 0.
∂t

(9)The explicit integration in time of the ordinary differential equation resulting
from the discretization of the spatial operator V · ∇x is subject to conditional stability. Linear stability conditions are usually more restrictive than
non-linear stability conditions, so the stability results obtained for the conservation law (9) can be applied to non-linear partial differential equations(PDE’s). These conditions can further be generalized to systems of hyperbolic PDE’s by considering that u corresponds to a characteristic variableof the system and V corresponds to the associated characteristic velocity.
The stability of the system is then driven by the characteristic leading to themost restrictive stability condition. However, the systems of multidimensional PDE’s of practical interest often have an infinite set of characteristics
forming a Monge cones instead of degenerating into lines. Thus, we are interested in a stability condition that takes all possible directions of V into
consideration, i.e. a condition on the norm kV k of V .
This so-called Courant-Friedrichs-Levy (CFL) condition can be definedas:kV k ∆t ≤ C ∆x
(10)24

where ∆t is the stable time step, ∆x is a length that represents the element
size, and C a constant that depends on the numerical method [30]. This condition should be fulfilled everywhere on the element. In the case of straightsided triangular elements, ∆x is usually chosen as the inner radius of thetriangle.
How does this size definition extend to curvilinear meshes? In order tofind a stability condition of the advection problem on the curved element,
we transform it into a problem on a straight-sided element for which we canuse the stability condition (10). The conservation law (9) is written in asystem of coordinates, x(X), which can be seen as the transformation ofanother one, X corresponding to a straight-sided element. The Jacobianof this transformation is A = x,X . In the system of coordinates X, theconservation law (9) reads:∂u+ V · A−1 ·∇X u = 0.
∂t | {z0 }

(11)

V

The problems (9) and (11) being equivalent, the stability condition (10) canbe written in the system of coordinates X as:
kV 0 k ∆t = kV · A−1 k ∆t < C ∆X
where ∆X is the size of the straight-sided element. Let lmin be the smallest
1
singular value of A so that lminis the largest singular value of A−1 . We have:kV A−1 k <

kV k
.
lmin

The condition:

∆X
.
(12)
lmin
is then sufficient to satisfy the stability condition (12). Condition (10) can beseen as the usual CFL condition (10) with an element size ∆x = ∆X/lmin
depending on the smallest singular value lmin of the mapping A, i.e. thesmallest eigenvalue of the metric tensor M = AAT , as defined in Section 6.1.
kV k ∆t < C

6.3. Application on a large scale example
In order to show the interest of the deformation measure introduced in
Section 6.1 and Section 6.2 for practical applications, we apply the mesh
optimization procedure to a large-scale oceanic simulation.
The patch of Figure 13 is a very small part of a larger mesh of the Great
Barrier Reef (GBR), in Australia. The GBR is composed of more than 2500
25

reefs and islands. This complex topography generates circulation patterns
on a wide range of scales and makes the generation of a curved mesh of thisdomain particularly challenging. Second-, third- and fourth-order versions of
two linear meshes of this region are considered. The first mesh is composedof about 330,000 triangles and the second one, even coarser, is composed ofabout 88,000 triangles. Those meshes are presented on Figure 17. In order totest our method on a difficult case, we use in both meshes a coarse resolution
compared to the small features of the underlying geometrical model, so thathighly curved elements are generated near the boundaries.
On each mesh, we determine the maximum stable time step to solve
the shallow water equations using a discontinuous Galerkin method of thesame polynomial order as the mesh with a fourth-order explicit Runge-Kutta
temporal scheme. The details of the parameterization and the boundary conditions can be found in [31]. For each mesh, an estimation of the maximum
time step is progressively refined through a simple bisection method that
determines the stability by running the simulation for a few iterations andchecking the global norm of the solution.
Table 2 shows the different time steps obtained on linear meshes, on
curvilinear meshes optimized with respect to the scaled Jacobians only andon meshes further optimized with respect to λmin . By chance, the Jacobianbased optimization is enough to recover the same time step on the seconde /J e = 0.5.
order mesh as on the linear mesh with a requirement of Jmin0
However, this is not the case for the other meshes, even when a high scaled
Jacobian is imposed. On the contrary, on all meshes optimized first withrespect to the scaled Jacobians and then with respect to λmin (using an
objective value of 0.5 in both optimizations), the maximum stable time stepis very close to the maximum stable time step of the same finite element
method on the linear mesh. The CPU time taken by the optimization process
ranges from 5 seconds for the second-order coarse mesh to about 10 minutes
for the fourth-order fine mesh.
7. Conclusions
The ability of generating high order/curvilinear mesh is a prerequisite to
the generalization of high order numerical methods in engineering analysis.
This paper offers a framework that allows to generate such meshes in a robustand efficient way.
The methodology consists in solving a sequence of optimization problems in which the deformation of high-order mesh elements, compared totheir linear counterpart, is minimized. The deformation measure is based26

2nd order coarse mesh

4th order coarse mesh

4th order fine mesh

Figure 17: Top: fourth-order coarse curvilinear mesh (88k triangles) of the Great Barrier
Reef, Australia and detail on the Whitsunday islands. Below: closer zoom on Hamilton
Island for different meshes.

27

p

Coarse

Fine

2
3
4
2
3
4

e /J eJmin0
0.1 0.5
1.7 3.0
0.8 0.9
0.5 0.5
0.8 1.0
0.4 0.5
0.4 0.5

λmin
0.1 0.5
2.5 3.1
2.0 2.0
1.2 1.2
1.4 3.1
1.4 2.0
0.9 1.3

Linear
3.0
1.9
1.2
3.0
2.0
1.4

Table 2: Maximum stable time step to solve the shallow water with discontinuous Galerkinmethod of polynomial order 2 to 4 and an explicit fourth order Runge-Kutta, on linear
meshes, curvilinear meshes optimized for the Jacobians only and curvilinear meshes further
optimized for the metric. The optimization process stops when all the Jacobians or the
eigenvalues of the metric are greater than 0.1 or 0.5.

on the Jacobian of the transformation that maps a high-order element onto
the straight-sided one, which enables to strictly controls the mesh validity.
Optionally, the minimum eigenvalue of the corresponding metric can be used
in a second step, with positive impact on simulations involving explicit time
stepping. Constraints on the maximum deformation are imposed throughmoving barriers, in order to guarantee the quality of the high-order meshgenerated. The restriction of the optimization process to blobs of elements
surrounding a tangled element, instead of the whole mesh, make the methods
computationally efficient.
We have tested the new untangling procedure on various other examples
than those presented in this paper, and our conclusions can be summarized
as follows:• The use of moving log-barriers for imposing Jacobian constraints isclearly a better solution than using constrained optimization, both interms of computational cost and of robustness.
• The energy part (6) of the objective function is not a critical parameterin the untangling process, its only purpose being most of the time toensure a unique solution to the optimization process.
• It is of paramount importance to allow mesh vertices to move on themodel entities they are classified on. This is especially true in 3D where
both surface and volume meshes should be untangled at once.
• The untangling of quadratic meshes using the new technique is very28

fast, robust and not sensitive to any parameter that we have introduced. We consider that issue to be well addressed.
All the material that is presented in this paper is readily available in Gmsh,
the open source mesh generator[9]. We are confident that users will help us
to improve further Gmsh’s high-order mesh generation capabilities. Nevertheless, some direct extensions of this work are already planned.
First, the objective function should somehow take into account the geometrical accuracy of the high order mesh representations. High order vertices
should allow to reduce some norm of the distance between the mesh and the
CAD model.
Then, the influence of element shapes on the numerical solutions should
be investigated. In this paper, we have shown that it is possible to computea length quantity representing the size of a curvilinear element and use it tocompute the CFL condition. We should now investigate the influence of meshcurving on the quality of numerical solutions, especially in computational
fluid dynamics.
Acknowledgement
This research project was funded in part by the Walloon Region under
WIST 3 grant 1017074 DOMHEX (Dominant Hexahedral Mesh Generation)and in part by the European Union under the 7th PCRD grant IDIHOM(Industrialisation of High-Order Methods - A Top-Down Approach).
[1] A. Johnen, J.-F. Remacle, C. Geuzaine, Geometrical validity of curvilinear finite elements, in: W. R. Quadros (Ed.), Proceedings of the 20th
International Meshing Roundtable, Springer Berlin Heidelberg, 2012,
pp. 255–271. doi:10.1007/978-3-642-24734-7_14.
[2] A. Johnen, J.-F. Remacle, C. Geuzaine, Geometric validity of curvilinear
finite elements, Journal of Computational Physics.
[3] B. Cockburn, G. Karniadakis, C.-W. Shu (Eds.), Discontinuous
Galerkin Methods, Vol. 11 of Lecture Notes in Computational Scienceand Engineering, Springer, Berlin, 2000.
[4] N. Kroll, H. Bieler, H. Deconinck, V. Couaillier, H. Van Der Ven,
K. Sorensen, ADIGMA – A European Initiative on the Developmentof Adaptive Higher-Order Variational Methods for Aerospace Applications: Results of a Collaborative Research Project Funded by the29European Union, 2006-2009, Notes on Numerical Fluid Mechanics andMultidisciplinary Design, Springer, 2010.
[5] F. Bassi, S. Rebay, High-order accurate discontinuous finite element
solution of the 2D Euler equations, J. Comput. Phys. 138 (2) (1997)
251–285. doi:10.1006/jcph.1997.5454.
[6] P.-E. Bernard, J.-F. Remacle, V. Legat, Boundary discretization for
high-order discontinuous Galerkin computations of tidal flows aroundshallow water islands, International Journal for Numerical Methods in
Fluids 59 (5) (2009) 535–557. doi:10.1002/fld.1831.
[7] T. Toulorge, W. Desmet, Curved boundary treatments for the discontinuous Galerkin method applied to aeroacoustic propagation, AIAA J.
48 (2) (2010) 479–489. doi:10.2514/1.45353.
[8] J.-F. Remacle, M. S. Shephard, An algorithm oriented mesh database,
International Journal for Numerical Methods in Engineering 58 (2)(2003) 349–374. doi:10.1002/nme.774.
[9] C. Geuzaine, J.-F. Remacle, Gmsh: a three-dimensional finite element
mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (11) (2009)
1309–1331.
[10] P. L. George, H. Borouchaki, Construction de maillages de degré 2 partie
3 : Tétraèdre p2, Tech. rep., INRIA (France) (2011).
[11] P. M. Knupp, Winslow smoothing on two-dimensional unstructured meshes, Engineering with Computers 15 (1999) 263–268.
doi:10.1007/s003660050021.
[12] Z. Q. Xie, R. Sevilla, O. Hassan, K. Morgan, The generation of arbitrary order curved meshes for 3d finite element analysis, Computational
Mechanics.
[13] P.-O. Persson, J. Peraire, Curved mesh generation and mesh refinement
using lagrangian solid mechanics, in: Proceedings of the 47th AIAA
Aerospace Sciences Meeting and Exhibit, Orlando (FL), USA, 5-9 January 2009, 2009.
[14] X. Li, M. Shephard, M. Beall, Accounting for curved domains in mesh
adaptation, International Journal for Numerical Methods in Engineering
58 (2) (2003) 247–276.
30

[15] X. Luo, M. Shephard, R. O’Bara, R. Nastasia, M. Beall, Automaticp-version mesh generation for curved domains, Engineering with Computers 20 (3) (2004) 273–285.
[16] M. Shephard, J. Flaherty, K. Jansen, X. Li, X. Luo, N. Chevaugeon,
J. Remacle, M. Beall, R. O’Bara, Adaptive mesh generation for curved
domains, Applied Numerical Mathematics 52 (2) (2005) 251–271.
[17] L. A. Freitag, P. Plassmann, Local optimization-based simplicial mesh
untangling and improvement, International Journal for Numerical
Methods in Engineering 49 (1) (2000) 109–125.
[18] E. J. López, N. M. Nigro, M. A. Storti, J. A. Toth, A minimal element distortion strategy for computational mesh dynamics, International Journal for Numerical Methods in Engineering 69 (2007) 1898–
1929.
[19] E. J. López, N. M. Nigro, M. A. Storti, Simultaneous untangling andsmoothing of moving grids, International Journal for Numerical Methods in Engineering 76 (7) (2008) 994–1019.
[20] L. Diachin, P. Knupp, T. Munson, S. Shontz, A comparison of two
optimization methods for mesh quality improvement, Engineering with
Computers 22 (2) (2006) 61–74.
[21] L. Freitag, P. Knupp, T. Munson, S. Shontz, A comparison of optimization software for mesh shape-quality improvement problems., Tech. rep.,
Argonne National Lab., IL (US) (2002).
[22] A. Fiacco, G. McCormick, Nonlinear programming: sequential unconstrained minimization techniques, Vol. 4, Society for Industrial Mathematics, 1990.
[23] A. Waechter, C. Laird, F. Margot, Y. Kawajir, Introduction to ipopt:
A tutorial for downloading, installing, and using ipopt (2009).
[24] D. Liu, J. Nocedal, On the limited memory bfgs method for large scale
optimization, Mathematical programming 45 (1) (1989) 503–528.
[25] R. Fletcher, C. Reeves, Function minimization by conjugate gradients,
The computer journal 7 (2) (1964) 149–154.
[26] J. Hesthaven, T. Warburton, Nodal discontinuous Galerkin methods:
algorithms, analysis, and applications, Vol. 54, Springer-Verlag New
York Inc, 2008.
31

[27] P. Vos, S. Sherwin, R. Kirby, From h to p efficiently: Implementing finite
and spectral hp element methods to achieve optimal performance forlow-and high-order discretisations, Journal of Computational Physics
229 (13) (2010) 5161–5181.
[28] P. Bernard, J. Remacle, R. Comblen, V. Legat, K. Hillewaert, Highorder discontinuous galerkin schemes on general 2d manifolds applied tothe shallow water equations, Journal of Computational Physics 228 (17)(2009) 6514–6535.
[29] T. Toulorge, W. Desmet, Spectral properties of discontinuous Galerkinspace operators on curved meshes, in: J. S. Hesthaven, E. M. Rønquist, T. J. Barth, M. Griebel, D. E. Keyes, R. M. Nieminen, D. Roose,
T. Schlick (Eds.), Spectral and High Order Methods for Partial Differential Equations, Vol. 76 of Lecture Notes in Computational Science and Engineering, Springer Berlin Heidelberg, 2011, pp. 495–502.
doi:10.1007/978-3-642-15337-2_48.
[30] B. Seny, J. Lambrechts, R. Comblen, V. Legat, J.-F. Remacle, Multirate
time stepping for accelerating explicit discontinuous galerkin computations with application to geophysical flows, International Journal for
Numerical Methods in Fluids in press. doi:10.1002/fld.3646.
[31] J. Lambrechts, E. Hanert, E. Deleersnijder, P.-E. Bernard, V. Legat,
J.-F. Remacle, E. Wolanski, A multiscale model of the whole Great
Barrier Reef hydrodynamics, Estuarine, Coastal and Shelf Science 79(2008) 143–151. doi:10.1016/j.ecss.2008.03.016.

32