Geometrical Validity of Curvilinear Finite Elements
A. Johnen1 , J.-F. Remacle2 and C. Geuzaine1
1

Université de Liège, Department of Electrical Engineering and Computer Science,Montefiore Institute B28, Grande Traverse 10, 4000 Liège, Belgium
2
Université catholique de Louvain, Institute of Mechanics, Materials and Civil
Engineering (iMMC), Place du Levant 1, 1348 Louvain-la-Neuve, Belgium

Abstract
In this paper, we describe a way to compute accurate bounds on Jacobian determinants of curvilinear polynomial finite elements. Our condition enables
to guarantee that an element is geometrically valid, i.e., that its Jacobian
determinant is strictly positive everywhere in its reference domain. It also
provides an efficient way to measure the distortion of curvilinear elements.
The key feature of the method is to expand the Jacobian determinant usinga polynomial basis, built using Bézier functions, that has both propertiesof boundedness and positivity. Numerical results show the sharpness of ourestimates.
Keywords: Finite element method; high-order methods; mesh generation;Bézier functions
1. Introduction
There is a growing consensus in the Finite Element community thathigher-order discretization methods will replace at some point the solvers
of today, at least for part of their applications. These high-order methods
require a good accuracy of the geometrical discretization to be accurate—
in other words, such methods will critically depend on the availability ofhigh-quality curvilinear meshes.
The usual way of building such curvilinear meshes is to first generate astraight sided mesh. Then, mesh entities that are classified on the curved
boundaries of the domain are curved accordingly [1–3]. Some internal mesh
Email address:a.johnen@ulg.ac.be,jean-francois.remacle@uclouvain.be,cgeuzaine@ulg.ac.be(A. Johnen1 , J.-F. Remacle2 and C. Geuzaine1 )Preprint submitted to Journal Of Computational Physics

July 4, 2012

entities may be curved as well. If we assume that the straight sided meshis composed of well shaped elements, curving elements introduces a “shape
distortion” that should be controlled so that the final curvilinear mesh is also
composed of well shaped elements. The optimization of the shape distortion
is a computationally expensive operation, especially when applied globallyover the full mesh. It is thus crucial to be able to get fast and accuratebounds on the distortion in order to 1) evaluate the quality of the elements
during the optimization process; and 2) reduce the sets of elements to be
optimized, so that the optimization can be applied locally, i.e., only where
it is necessary.
In this paper we present a method to analyze curvilinear meshes in termsof their elementary Jacobian determinants. The method does not deal withthe actual generation/optimization of the high order mesh. Instead, it provides an efficient way to guarantee that each curvilinear element is geometrically valid, i.e., that its Jacobian determinant is strictly positive everywherein its reference domain. It also provides a way to measure the distortionof the curvilinear element. The key feature of the method is to adaptively
expand the elementary Jacobian determinants in a polynomial basis thathas both properties of boundedness and positivity. Bézier functions areused to generate these bases in a recursive manner. The proposed method
can be either used to check the validity and the distortion of an existing
curvilinear mesh, or embedded in the curvilinear mesh generation procedure to assess the validity and the quality of the elements on the fly. The
algorithm described in this paper has been implemented in the open sourcemesh generator Gmsh [4], where it is used in both ways.
2. Curvilinear Meshes, Distortion and Bounds on Jacobian Determinants
Let us consider a mesh that consists of a set of straight-sided elementsof order p. Each element is defined geometrically through its nodes xi ,(p)
i = 1, . . . , Np and a set of Lagrange shape functions Li (ξ), i = 1, . . . , Np .
The Lagrange shape functions (of order p) are based on the nodes xi andallow to map a reference unit element onto the real one:x(ξ) =Np
X(p)

Li (ξ) xi .

(1)

i=1

The mapping x(ξ) should be bijective, which means that it should admit aninverse. This implies that the Jacobian determinant det x,ξ has to be strictly
2

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 1: Reference unit triangle in local coordinates ξ = (ξ, η) and the mappings x(ξ),X(ξ) and X(x).

positive. In all what follows we will always assume that the straight-sided
mesh is composed of well-shaped elements, so that the positivity of det x,ξ
is guaranteed. This standard setting is presented on Figure 1 (left) for thequadratic triangle.
Let us now consider a curved element obtained after application of thecurvilinear meshing procedure, i.e., after moving some or all of the nodes ofthe straight-sided element. The nodes of the deformed element are called
Xi , i = 1 . . . Np , and we haveX(ξ) =Np
X(p)

Li (ξ) Xi .

(2)

i=1

Again, the deformed element is assumed to be valid if and only if the Jacobian determinant J(ξ) := det X,ξ is strictly positive everywhere over theξ reference domain. The Jacobian determinant J, however, is not constantover the reference domain, and computing Jmin := minξ J(ξ) is necessary toensure positivity.
The approach that is commonly used is to sample the Jacobian determinant on a very large number of points. Such a technique is however both
expensive and not fully robust since we only get a necessary condition. In
this paper we follow a different approach: because the Jacobian determinant
J is a polynomial in ξ, J can be interpolated exactly as a linear combinationof specific polynomial basis functions over the element. We would then like3

to obtain provable bounds on Jmin by using the properties of these basis
functions.
In addition to guaranteeing the geometrical validity of the curvilinear
element, we are also interested in quantifying its distortion, i.e., the deformation induced by the curving. To this end, let us consider the transformation X(x) that maps straight sided elements onto curvilinear elements (seeFigure 1). It is possible to write the determinant of this mapping in termsof the ξ coordinates as:
det X,x =det X,ξ
J(ξ)=.
det x,ξ
det x,ξ

(3)

We call X(x) the distortion mapping and its determinant δ(ξ) := det X,xthe distortion. The distortion δ should be as close to δ = 1 as possible inorder not to degrade the quality of the straight sided element. Elementsthat have negative distortions are of course invalid but elements that havedistortions δ  1 or δ  1 lead to some alteration of the conditioning of thefinite element problem. In order to guarantee a reasonable distortion it isthus necessary to find a reliable bound on Jmin and Jmax := maxξ J(ξ) over
the whole element.
Note that many different quality measures can be defined based on the
Jacobian determinant J. For example, one could look at the Jacobian determinant divided by its average over the element instead of looking at the
distortion. Obtaining bounds on Jmin and Jmax is thus still the main underlying challenge. Other interesting choices are presented and analyzed in[5].
3. Bounds for Second Order Planar Triangles
We start our analysis with the particular case of second order planar
triangles for which a direct computation of Jmin is relatively easy. The
determinant J(ξ) = J(ξ, η) for a planar triangle at order p is a polynomialin ξ and η of order at most 2 (p − 1). For quadratic planar triangles, J(ξ, η)
is therefore quadratic at most in ξ and η.
The geometry of the six-node quadratic triangle is shown in Figure 1.
Inspection reveals two types of nodes: corners (1, 2 and 3) and midside nodes(4, 5 and 6). If Ji is defined as J(ξ, η) evaluated at node i, it is possible towrite the Jacobian determinant exactly as a finite element expansion whose4coefficients are the Jacobian determinants at the nodes:
J(ξ, η) = J1 (1 − ξ − η)(1 − 2ξ − 2η) + J2 ξ(2ξ − 1) + J3 η(2η − 1) +| {z }| {z }
|
{z
}
(2)

(2)

L1 (ξ,η)L2 (ξ,η)J4 4(1 − ξ − η)ξ + J5
|
{z
}
(2)

4ξη
|{z}

(2)

+ J6 4(1 − ξ − η)η .
|{z}

L5 (ξ,η)L4 (ξ,η)(2)

L3 (ξ,η)(4)(2)

L6 (ξ,η)(2)

In equation (4), the functions Li (ξ, η) are the equidistant quadratic Lagrange shape functions that are commonly used in the finite element community [6].
We first show how to compute the exact minimal Jacobian determinant
Jmin . Then we examine different bounds that can be provided on Jmin by
exploiting the properties of the basis used in the expansion.
3.1. Exact Computation of Jmin
From equation (4), the stationnary point of J can be computed by solving
∂J
∂J
=
= 0,
∂ξ
∂η

(5)

which leads to the following linear system of two equations and two unknowns ξsta and ηsta :



ξsta
4(J1 + J2 − 2J4 )4(J1 − J4 + J5 − J6 )(6)4(J1 − J4 + J5 − J6 ) 4(J1 + J3 − 2J6 )ηsta
−(−3J1 − J2 + 4J4 ).
=−(−3J1 − J3 + 4J6 )Algorithm 1 allows to compute the minimal Jacobian determinant over one
quadratic planar element exactly. If the minimum of the function is outsideof the element, it computes the minimum on its border assuming a function
MINQ(a, b, c) that computes
MINQ(a, b, c) = min a x2 + b x + c.
x∈[0,1](7)

Although Algorithm 1 is quite simple, applying similar techniques forhigher order elements would become extremely expensive computationally.
For example, for a third order triangle, the Jacobian determinant is of order
4 and the algorithm requires the solution of a system of cubic equations; atorder 4, it requires the solution of a system of quintic equations. Insteadof trying to evaluate Jmin directly, we should try to compute (the sharpest
possible) bounds in a computationally efficient manner.
5

Algorithm 1: Exact computation of Jmin over a quadratic triangle
1 compute nodal Jacobian determinants Ji , i = 1, . . . , 6;
2 compute ξsta , ηsta as in equation (6);3 if ηsta > 0 and ξsta > 0 and 1 − ξsta − ηsta > 0 then
4

Jmin = min(J(ξsta , ηsta ), J1 , J2 , J3 );5 else6789m1 = MINQ(2(J1 + J2 − 2J4 ), −3J1 − J2 + 4J4 , J1 );m2 = MINQ(2(J1 + J3 − 2J6 ), −3J1 − J3 + 4J6 , J1 );m3 = MINQ(2(J2 + J3 − 2J5 ), −3J2 − J3 + 4J5 , J2 );Jmin = min (m1 , m2 , m3 );10 return Jmin ;3.2. The Principle for Computing Bounds on Jmin
It is obvious that a necessary condition for having J(ξ, η) > 0 everywhere is that Ji > 0, i = 1, . . . , 6. Yet, this condition is not sufficient.
The expression (4) does not give more information because the quadratic(2)
Lagrange shape functions Li (ξ, η) change sign on the reference triangle.
What polynomial basis should we chose to obtain usable bounds?
The first idea is to expand (4) into monomials, which gives:
J(ξ, η) = J1 + (−3J1 − J2 + 4J4 )ξ + (−3J1 − J3 + 4J6 )η+4(J1 − J4 + J5 − J6 )ξη + 2(J1 + J2 − 2J4 )ξ 2 + 2(J1 + J3 − 2J6 )η 2 . (8)Every monomial being positive on the reference triangle, we have now a setof sufficient conditions that can be written as
4J4 ≥ 3J1 +J2 , 4J6 ≥ 3J1 +J3 , J1 +J5 ≥ J4 +J6 , J1 +J2 ≥ 2J4 , J1 +J3 ≥ 2J6 .
However these constraints do not provide a usable bound on Jmin and breakthe symmetry of the expression with respect to a rotation of corner nodes.
A second idea is to expand (4) in terms of the second order hierarchicalbasis functions ψi (ξ, η), i = 1, . . . , 6, which are also positive on the triangle[7]:
J(ξ, η) = J1 (1 − ξ − η) +J2 ξ +J3 η +(4J4 −2J1 −2J2 ) (1 − ξ − η)ξ +
|{z}
|{z}
| {z }|{z}ψ1 (ξ,η)ψ2 (ξ,η)ψ3 (ξ,η)ψ4 (ξ,η)(4J5 − 2J3 − 2J2 ) ξη +(4J6 − 2J1 − 2J3 ) (1 − ξ − η)η . (9)
|{z}
|
{z
}
ψ5 (ξ,η)ψ6 (ξ,η)6

This last expression has the right symmetry, and leads to the following
validity conditions:J1 ≥ 0, J2 ≥ 0, J3 ≥ 0, 4J4 ≥ 2J1 +2J2 , 4J5 ≥ 2J2 +2J3 , 4J6 ≥ 2J3 +2J1 .
(10)J(ξ, η) is a degree two polynomial, therefore it has an expression in the(ΨP6i ) basis. Let Ki denote the coefficients in this basis. Writing J(ξ, η) :=i=1 ψi (ξ, η)Ki , we have!
!
XX
ψi (ξ, η) min Ki = min Ki ,min J(ξ, η) = min
ψi (ξ, η)Ki ≥ min
ξ,η

ξ,η

ξ,η

i

i

i

i

P
because i ψi = 1 + ξ + η − ξ 2 − η 2 − ξη has its minimum on the corner
nodes (where its value is equal to 1). And since Ki , i = 1, . . . , 3 are values of
the Jacobian determinant (at the three corners), they form an upper boundon it. Thus, expansion (9) leads to the following estimate for the minimumof the Jacobian determinant over the triangle:
Jmin ≥ min{J1 , J2 , J3 , 4J4 − 2J1 − 2J2 , 4J5 − 2J2 − 2J3 , 4J6 − 2J3 − 2J1 }≤ min{J1 , J2 , J3 }.

(11)

It is easy to see that the estimate is however of very poor quality: foran element that has a constant and positive J, (11) simply tells us that
Jmin ≥ 0.
In order to find a sharper estimate, instead of the hierarchical quadratic
functions ψi (ξ, η), we can use the quadratic triangular Bézier functions(2)
B2 (ξ, η) [8]:J(ξ, η) = J1 (1 − ξ − η)2 +J2 ξ 2 +J3 η 2 +
|{z}
|{z}
|
{z
}
(2)

B1 (ξ,η)(2)

B2 (ξ,η)(2)

B3 (ξ,η)





1
1
2J4 − (J2 + J1 ) 2 ξ (1 − ξ − η) + 2J5 − (J3 + J2 )2 ξη +
|{z}
|
{z
}
2
2
(2)

(2)

B5 (ξ,η)B4 (ξ,η)
12J6 − (J1 + J3 ) 2 η (1 − ξ − η) . (12)
|
{z
}
2
(2)

B6 (ξ,η)

Since

(2)
(2)
i=1 Bi (ξ, η) = 1 and Bi (ξ, η) ≥ 0, we obtain the following esti-

P6

7

mate



J1 + J2
J2 + J3
J3 + J1
Jmin ≥ min J1 , J2 , J3 , 2J4 −, 2J5 −
, 2J6 −
2
2
2
≤ min {J1 , J2 , J3 } .
(13)
One can show that this estimate is always better than the one using the hierarchical basis. It provides two conditions on the geometrical validity of the triangle: a sufficient condition (if min{J1 , J2 , J3 , 2J4 −J1 +J2J2 +J3J3 +J12 , 2J5 −2 , 2J6 −2 } > 0, the element is valid) and a necessary condition (if min{J1 , J2 , J3 } < 0, the element is invalid). However, these two conditions are sometimes insufficient to determine the validity of the element, as the bound (13) is often not sharp enough (having
J2 +J3J3 +J12
min{2J4 − J1 +J
2 , 2J5 −2 , 2J6 −2 } < 0 does not imply that theelement is invalid).
A sharp necessary and sufficient condition on the geometrical validityof an element can be achieved in a general way by refining the Bézier estimate adaptively so as to achieve any prescribed tolerance—and thus provide
bounds as sharp as necessary for a given application.
4. Adaptive Bounds for Arbitrary Curvilinear Finite Elements
In order to explain the adaptive bound computation let us first focus on
the one-dimensional case, for “line” finite elements. Since Bézier functions
can be generated for all types of common elements (triangles, quadrangles,
tetrehedra, hexahedra and prisms), the generalization to 2D and 3D elements will be straightforward.
4.1. The One-Dimensional Case
In 1D the Bézier functions are the Bernstein polynomials:
 
n
(n)
Bk (ξ) =(1 − ξ)n−k ξ k(ξ ∈ [0, 1] ; k = 0, ..., n)
(14)
k

n!
where nk = n!(n−k)!
is the binomial coefficient. The Bézier interpolation
requires n + 1 control values bi . We haveJ(ξ) =Nn
X(n)

Bk (ξ) bk .

(15)

k=0

Bernstein-Bézier functions have the nice following properties : (i) they formP(n)a partition of unity which means that nk=0 Bk (ξ) = 1 for all ξ ∈ [0, 1] and8(n)(ii) they are positive which means that Bk (ξ) ≥ 0 for all ξ ∈ [0, 1]. This
leads to the well known property of Bézier interpolations:min J(ξ) ≥ bmin = min bi andmax J(ξ) ≤ bmax = max bi .

iξ∈[0,1]iξ∈[0,1](16)

Moreover, the control values related to the corner nodes of the element are
equal to the values of the interpolated function. In what follows we assume
that these “corner” control values are always ordered at the Kf first indices.
We then havemin J(ξ) ≤ min bi andξ∈[0,1]i<Kf

max J(ξ) ≥ max bi .
i<Kf

ξ∈[0,1](17)

Since Lagrange and Bézier functions span the same function space, computation of the Bézier values bi from the nodal values Ji (and convertly) is(n)done by a transformation matrix. The tranformation matrix T B→L , which
computes nodal values from control values, is created by evaluating Bézierfunctions at sampling points:
 (n)

(n)
B0 (ξ0 ) . . . Bn (ξ0 ) (n) B0 (ξ1 ) . . . Bn(n) (ξ1 ) (n).
T B→L = 
..
..
..

.
.
.


(n)

(n)

B0 (ξn ) . . . Bn (ξn )Those sampling points are taken uniformly, i.e. at the location of the nodes(n)(n) −1of the element of order n. The inverse transformation is T L→B = T B→L
and from the expression of the interpolation of the Jacobian determinant(15), we can write(n)J = T B→L B(n)

B = T L→B J ,(18)

where B and J are the vectors containing respectively the bi ’s and the Ji ’s.
4.2. Adaptive Subdivision
Let us assume that the domain [0, 1] is subdivided into Q parts. Theinterpolation J [q] (ξ [q] ) on the q th subdomain [a, b] (0 ≤ a < b ≤ 1) mustverify[q][q]

J (ξ ) =

Nn
X
k=0

(n)
[q]
Bk (ξ [q] ) bk =

Nn
X

(n)

Bk (ξ(ξ [q] )) bk

k=09(ξ [q] ∈ [0, 1]),(19)

with ξ(ξ [q] ) = a + (b − a) ξ [q] .
[q][q]Considering the nodes ξk such that ξk = ξk (k = 0, . . . , n) (i.e., suchthat they are ordered like the sampling points), the expression (19) reads
 (n)

(n)
B0 (a + (b − a) ξ0 ) . . . Bn (a + (b − a) ξ0 )
 (n)
 B0 (a + (b − a) ξ1 ) . . . Bn(n) (a + (b − a) ξ1 ) 
[q]
(n)
[q]

 B = T (n)T B→L B = 
B,
..
.
B→L
..

..
.
.


(n)
(n)
B0 (a + (b − a) ξn ) . . . Bn (a + (b − a) ξn )where B [q] is the vector containing the control values of the q th subdomain.
This implies that(n)(n) [q][q]
B = T L→B T B→L
B = M[q] B.
(20)
Each set of new control values bounds the Jacobian determinant on its own
subdomain and we have:[q]b0min = min bi
i,q

and

[q]

max bii<Kf ,q

[q]

(21)

[q]

(22)

≤ Jmin ≤ min bi
i<Kf ,q≤ Jmax ≤ b0max = max bi .
i,qIf an estimate is not sufficiently sharp, we can thus simply subdivide theappropriate parts of the element. This leads to a simple adaptive algorithm,
exemplified in Figure 2. In this particular case the original estimate (16)-(17)is not sharp enough (Jmin ∈ [−3, 1]). After one subdivision, the Jacobian
determinant is proved to be positive on the second subdomain. The first
subdomain is thus subdivided once more, which proves the validity. In
practice, as will be seen in Section 5, a few levels of refinement lead to thedesired accuracy. The subdivision has quadratic speed of convergence [9, 10].
Note that in a practical implementation (in finite precision arithmetic)
we must take care of a tricky situation. If the minimum of the Jacobian determinant is too close to zero but positive, then the upper bound is positive
while the lower bound might never get positive. In order to avoid this situation, we limit the number of consecutive subdivisions that can be applied.
The undetermined elements are then considered as invalid. Another way ofgetting rid of this issue is to relax the condition of rejection as explained in
section 4.4.

10

5

1

4
0.8

x

3
2

0.6

J

0.4

1
0
ï1

0.2

ï2
0

0

0.2

0.4

0.6

0.8

ï3

1

0

0.2

j

0.4

0.6

0.8

1

j

Jmin

2.

1.

[-3,1][-0.25, 0.88][0.88, 0.88]

3. [0.22, 0.65] [0.65, 0.65]

Figure 2: Top left: One-dimensional element mapping x(ξ). Top right: Exact Jacobian
determinant J(ξ) (black), control values on the original control points (green) and two
adaptive subdivisions (blue and red). Bottom: Estimates of Jmin at each step in theadaptive subdivision process.

11

4.3. Extension to Higher Dimensions
The extension of the method to higher dimensions is straightforward,
provided that Bézier functions can be generated and that a subdivision
scheme is available. Jacobian determinants J are polynomials of ξ, η in 2D
and of ξ, η, ζ in 3D.
For high order triangles, the Bézier triangular polynomials are defined
as
 

p p−i
(p)
Ti,j (ξ, η) =
ξ i η j (1 − ξ − η)p−i−j(i + j ≤ p).
ij
It is possible to interpolate any polynomial function of order at most pon the unit triangle ξ > 0, η > 0, ξ + η < 1 as an expansion into Bézier
triangular polynomials. Recalling that, for a triangle at order p, its Jacobian
determinant J(ξ, η) is a polynomial in ξ and η at order at most n = 2(p − 1),
we can writeX(n)J(ξ, η) =bij Ti,j (ξ, η).
i+j≤n

It is also possible to compute J in terms of Lagrange polynomialsX(n)J(ξ, η) =Ji Li (ξ, η)iwhere the Ji are the Jacobian determinants calculated at Lagrange points.
n such that
It is then easy to find a transformation matrix TLBn
B = TLB
J,

where B and J are the vectors containing respectively the control values of
the Jacobian determinant bij and the Ji ’s. As an example, for quadratic
triangles we obtain
1
0
0 0 0 001
0 0 0 0 00
1 0 0 0 
2

,
(23)
TLB = 0 2 0 0 
 −1/2 −1/2

0 −1/2 −1/2 0 2 0 −1/2
0 −1/2 0 0 2
which directly provides the estimate (13).
Other element shapes can be treated similarly. For quadrangles, tetrahedra, prisms and hexahedra, the Bézier are functions respectively:
(p)

(p)

(p)

Qi,j (ξ, η) = Bi (ξ) Bj (η)
12(i ≤ p, j ≤ p),

Line
Triangle
Quadrangle
Tetrahedron
Prism
Hexahedron

Order (n) of J
p−1
2(p − 1)
2p − 13(p − 1)
3p − 1
3p − 1

Number of coefficientsn+1(n + 1)(n + 2)/2
(n + 1)2(n + 1)(n + 2)(n + 3)/6
(n + 1)2 (n + 2)/2
(n + 1)3

Table 1: Order of the Jacobian determinant and number of coefficients in the expansion
for an element of order p.

(p)Ti,j,k (ξ, η, ζ) =

 


p p−i p−i−jξ i η j ζ k (1 − ξ − η − ζ)p−i−j−kijk(i + j + k ≤ p),

(p)

(p)

(p)

(i + j ≤ p, k ≤ p)(p)(i ≤ p, j ≤ p, k ≤ p).

Pi,j,k (ξ, η, ζ) = Ti,j (ξ, η) Bk (ζ)
and
(p)

(p)

(p)

Hi,j,k (ξ, η, ζ) = Bi (ξ) Bj (η) Bk (ζ)

Matrices of change of coordinates can then be computed inline for every
polynomial order, and bounds of Jacobian determinants computed accordingly. Table 1 summarizes the order of the Jacobian determinant and thenumber of coefficients in its expansion for all common element types. In
all cases the subdivision scheme works exactly in the same way as for lines.
Figure 3 shows the first level of subdivision for a third-order triangle.
4.4. Implementation
As mentioned in Section 2, the bounds on the Jacobian determinant canbe used to either make the distinction between valid and invalid elementswith respect to a condition on Jmin , or to measure the quality of the elements
by systematically computing Jmin and Jmax with a prescribed accuracy.
In both cases the same operations are executed on each element. First,
the Jacobian determinant is sampled on a determined number of points Ns ,equal to the dimension of the Jacobian determinant space, and so to thenumber of Bézier functions. Second, Bézier values are computed. Then
adaptive subdivision is executed if necessary. Algorithms 2 and 3 show inpseudo-code the algorithm used to determine whether the Jacobian determinant of the element is everywhere positive or not.

13

1
0.8
0.6
y
0.4
0.2
0
ï0.2
ï0.2

0

0.2

0.4

0.6

0.8

1

1.2

x

2

2

1.5

1.5

1
J

J

1

0.5
0.5

0
ï0.5
0

0
0

0
0.5

0.5
1

0
0.5

0.5

1

d

1
j

1

d

j

Figure 3: Top: third-order planar triangle. Bottom left: exact Jacobian determinant andcontrol values (dots) on the original control points; the validity of the element cannot be
asserted. Bottom right: exact Jacobian determinant and control values (dots) after one
subdivision; the element is provably correct.

Algorithm 2: Check if an element is valid or invalid
Input: a pointer to an element.
Output: true if the element is valid, false if the element is invalid
1 set sampling points Pi , i = 1, . . . , Ns ;2 compute Jacobian determinants Ji at points Pi ;3 for i = 1 to Ns do
4

if Ji <= 0 then return false;

5 compute Bézier coefficients bi , i = 1, . . . , Ns using (18);6 i = 1;
7 while i ≤ Ns and bi > 0 do
8i = i + 1;

9 if i > Ns then return true;
10 call algorithm 3 with bi as arguments and return output;

14

Algorithm 3: Compute the control values of the subdivisions
Input: Bézier coefficients bi , i = 1, . . . , Ns
Output: true if the Jacobian determinant on the domain iseverywhere positive, false if not[q]
1 compute new Bézier coefficients bi , q = 1, . . . , Q as in equation (20);
2 for q = 1 to Q do
3
4for i = 1 to Kf do[q]
if bi <= 0 then return false;

5 for q = 1 to Q do
6

i = 1;

7

while i ≤ Ns and bi > 0 do
i = i + 1;

8
9
10
11

[q]

if i ≤ Ns then[q]call algorithm 3 with bi as arguments and store output;
if output = false then return false;

12 return true;

Algorithm 3 could be further improved by optimzing the loop on line
5, by first selecting q for which we have the best chance to have a negative
Jacobian determinant (line 4, algo 3). In practice this improvement is not
significant since the only case for which we can save calculation is for invalid
elements—and the proportion of them which require subdivision in order tobe detected is usually small (about 3% for the mesh depicted in Figure 4).
Note that we may also want to find, for example, all the elements for which
the Jacobian determinant is somewhere smaller than 20% of its average. We
then just have to compute this average and replace the related lines (4 and
7 for algorithm 2).
Another possible improvement is to relax the condition of rejection. We
could accept elements for which all control values are positive but reject anelement as soon as we find a Jacobian determinant smaller than a defined
percent of the average Jacobian determinant. The computational gain canbe significant, since elements that were classified as good and which neededa lot of subdivisions (and have a Jacobian determinant close to zero) will be
instead rapidly be detected as invalid.
More interestingly, the computation of sampled Jacobian determinants
and the computation of Bézier control values in algorithm 2 can easily be
15executed for a whole groups of elements at the same time. This allows to
use efficient BLAS 3 (matrix-matrix product) functions, which significantly
speeds up the computations.
The algorithm used for all the tests in the next section is implemented in
the open source mesh generator Gmsh [4] as the AnalyseCurvedMesh plugin.
5. Numerical Results
We start by comparing the new adaptive computation of bounds on
Jacobian determinants with the brute-force sampling of the Jacobian determinant for the detection of invalid high-order triangles.
The points at which we sample the Jacobian determinant for the bruteforce method are taken as the nodes of an element of order k, for k =1, . . . , 65, leading to a number of sampling points comprised between 3 and2211. In order to make the comparison as fair as possible, we have implemented the brute-force computation as efficiently as possible, i.e., fork (> n) sufficiently large we sample the Jacobian determinant on the points
computed for an element at order n (the order of the Jacobian determinant)and then compute the desired Jacobian determinant values by a matrixvector product, just like in our own adaptive method.
We consider the two-dimensional microstructure with oval holes depicted
in Figure 4, meshed with 82,009 sixth-order triangles. In this mesh 38,978
triangles are curved, and 8,112 are invalid. The new algorithm successfully
detects all the 8,112 invalid elements in 2.146 s on a 2.4GHz Core 2 Duo
laptop computer1 . Some elements needed as much as 7 levels of subdivisions
in order to be classified: see Table 2. The brute-force approach required
1, 596 sample points per triangle (k = 55) in order to detect all the invalid
elements, and took more than 7 times longer. But far worse, increasing thenumber of sampling points beyond 1, 596 can actually lead to a decreased
accuracy of the prediction, as shown in Figure 5.
Let us now examine the use of the adaptive bounds in the curvilinear
meshing algorithm. We consider the boundary layer mesh of the B-Spline
CAD model of a tri-part wing: see Figure 6. The cubic triangular meshis generated as follows. We first generate a straight-sided mesh (Figure
6/(a)). Then, every mesh edge that is classified on a model edge is curved
by snapping their center vertices on the model edge. High order nodes are
1
Note that for completeness the algorithm also analyzes straight-sided elements, whichis unnecessary in practice.

16

δ
ï1.5

ï1

ï0.5

0

0.5

1

Figure 4: Two-dimensional mesh with sixth order triangles; 47.5% of the elements are
curved. The straight elements are in green and the curved ones are colored in function of
the minimum of the distortion. The valid elements (δmin > 0) are colored in blue. The
invalid ones are colored in red if δmin is near 0 and in black if δmin < −1.5.

17

Curved Element Classification

First stage
1 subdiv.
2 subdiv.
3 subdiv.
4 subdiv.
5 subdiv.
6 subdiv.
7 subdiv.
Subtotal
Total

Valid
Invalid
elements
elements
29,715
8,039
+787
+0
+285
+17
+56
+15
+16
+22
+5
+15
+1
+2
+1
+2
30,866
8,112
38,978Underterminedelements
1,224
437
135
64
26
6
3
0

# elements
analysed atgiven stage

CPU time for
given stage [s]38,978
1,224
437
135
64
26
6
3

1.865
1.16e-1
8.40e-2
4.02e-22.40e-2
1.10e-2
4.34e-3
1.47e-3
2.146

Undetected invalid elements

Table 2: On the left, number of curved elements detected as valid or invalid at each stageof the adaptive algorithm; At the first stage, 8,309 elements can be classified as invalid due
to a negative value of at least one of the 66 sampling points. Then Bézier coefficients arecomputed and 29,715 elements are classified as valid because those coefficients are positive.
The 1,224 (3.14 %) remaining elements need to be subdivided adaptively. On the right,
computation time; most of the time is spent on sampling the Jacobian determinant andcomputing the first Bézier coefficients.

2

10

1

10

0

10

2

3

10

10
Sampling points

Figure 5: Number of undetected invalid elements using brute-force sampling of the Jacobian determinant. The five red data points correspond to the correct result, i.e., when no
invalid triangle is left undetected.

18

then inserted on every edge and in the middle of every face (Figure 6/(b)).
This simple procedure does not guarantee that the final mesh is valid. In
our case, 66 elements are invalid. Then, an optimization is applied globally(Figure 6/(c)). The final curvilinear mesh contains about 31% of curved
elements. During the meshing process, the adaptive bound computation
allowed to detect all invalid elements (the worst minimum of the distortion
that was observed was δmin = −27.72). After optimization, the final meshis composed of elements that have δmin > 0.64.
Finally, a 3D mesh is considered. The CAD model of an A319 plane ismeshed with 168,884 p3 tetrahedra (see figure 7 and 8). Without executing
any optimization, 76 elements are invalid. The new algorithm detects themin 9.88 s on a 2.4GHz Core 2 Duo laptop computer. The worst elements interm of their Jacobian determinant are located around leading edges, where
the curvature is the most important. 999 elements have δmin ≤ .7 and theworst distortion is δmin = −7.74.
6. Conclusions and Perspectives
In this paper we presented a way to compute accurate bounds on Jacobian determinants of curvilinear finite elements, based on the efficient
expansion of these Jacobian determinants in terms of Bézier functions.
The overall idea can be summarized as follows:
1. The Jacobian determinant of a polynomial element is also a polynomial(of higher order).
2. The Bézier polynomial basis satisfies the convex hull property, whichmeans that any polynomial expressed in this basis is bounded by thevalues at the control nodes.
3. By a change of polynomial basis (e.g. from Lagrange to Bézier), one
naturally gets the bounds of the Jacobian determinant.
4. If the bounds are not accurate enough, one can subdivide the element,
and once again, with a change of polynomial basis, obtain a moreaccurate bound in each sub-element (and so on and so forth until
enough precision is reached).
The proposed algorithm can either be used to determine the validity
or invalidity of curved elements, or provide an efficient way to measure
their distortion. Triangles, quadrangles, tetrahedra, prisms and hexahedra
can be analyzed using the same algorithm, which is available in the open
source mesh generator Gmsh. Numerical tests show that the method is19

(a) Straight boundary mesh(a) Straight-sided mesh(b) Raw curvilinear mesh (c) After optimization

δ
ï1.5

ï1

ï0.5

0

0.5

(b) Raw curvilinear mesh
50

300

Elements

30

ï1.5

20

ï1

ï0.5

Elements

δ

40

10
50

30
20

200

0

0.5

1

100
300

ï20

ï10
Distorsion

0

0

10
Elements

Elements

0
40
ï30

1

(c) After optimization

200

0.7

0.8

0.9
1
1.1
Distorsion

1.2

1.3

(d) Distortion before and after optimization
100

Figure10 6: Distortion of the curvilinear mesh of a wing (p3 triangles) before and after
00
optimization.
ï30
ï20
ï10
0
10
0.7 0.8 0.91
1.1 1.2 1.3
Distorsion

Distorsion(e) Distortion before and after optimization
Figure 5: Distortion of the curvilinear mesh of a wing (p3 triangles) before and after
optimization.

20
19

δ
ï1.5

ï1

ï0.5

0

0.5

1

Figure 7: The geometry of a A319 plane is meshed with p3 tetrahedron without executing
optimization. The curved elements are all shown in the figure at the bottom. For the
figure at the center, only elements for which δ ≤ 0.7 are shown. On 168,884 elements,
24,691 are curved and 76 are invalid.

21

δ
ï1.5

ï1

ï0.5

0

0.5

1

Figure 8: The top figure shows a cut of the mesh of the A319 plane (figure 7). The twobottom figures show same invalid elements. On the left, elements are colored in functionof their distortion. Here, we can clearly see that some faces intersect each other. In somecases, the Jacobian determinant is negative only inside the element so as it is not possible
to see it visually.

22

robust, and a user-defined error tolerance permits to adjust the accuracy vs.
computational time ratio.
Perspectives for future work are numerous. We are currently investigating two related areas: first, we are working on a generalization of thebounds presented in this paper to the case of surface meshes embedded in
3D (curved surfaces). Second, we are investigating the use of various optimization strategies to generate meshes with provably good qualities.
Acknowledgement
This research project was funded in part by the Walloon Region under
WIST 3 grant 1017074 (DOMHEX).
Authors gratefully thank E. Bechet from the University of Liège and K.
Hillewaert from Cenaero for insightful discussions about Bézier functions
and curvilinear mesh generation. Authors also thank V. D. Nguyen forproviding the microstructure geometry used in Figure 4.
References[1] S. Dey, R. M. O’Bara, M. S. Shephard, Curvilinear mesh generation in
3D, Computer Aided Geom. Design 33 (2001) 199–209.
[2] M. S. Shephard, J. E. Flaherty, K. E. Jansen, X. Li, X. Luo, N. Chevaugeon, J.-F. Remacle, M. W. Beall, R. M. O’Bara, Adaptive mesh generation for curved domains, Applied Numerical Mathematics 52 (2005)
251–271.
[3] S. J. Sherwin, J. Peiró, Mesh generation in curvilinear domains using
high-order elements, International Journal for Numerical Methods inEngineering 53 (2002) 207–223.
[4] 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.
[5] X. Roca, A. Gargallo-Peiró, J. Sarrate, Defining Quality Measures forHigh-Order Planar Triangles and Curved Mesh Generation, Proceedings of the 20th International Meshing Roundtable (2012) 365–383.
[6] T. Hughes, The Finite Element Method, Dover, 2003.

23

[7] I. Babus̆ka, B. Szabó, R. L. Actis, Hierarchic models for laminated composites, International Journal for Numerical Methods in Engineering 33(1992) 503–535.
[8] G. E. Farin, Curves and surfaces for CAGD: a practicle guide, MorganKaufmann, 2002.
[9] J. M. Lane, R. F. Riesenfeld, A theoretical development for the computer generation and display of piecewise polynomial surfaces, IEEE
Transactions on Pattern Analysis and Machine Intelligence 2 (1) (1980)
35–46.
[10] E. Cohen, L. L. Schumacker, Rates of convergence of control polygons,
Computer Aided Geometric Design 2 (1985) 229–235.

24