Geometrical Validity of High-Order Triangular
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{a.johnen, cgeuzaine}@ulg.ac.be2
Institute of Mechanics, Materials and Civil Engineering, Université catholique deLouvain, Avenue Georges-Lemaı̂tre 4, 1348 Louvain-la-Neuve, Belgium
jean-francois.remacle@uclouvain.be
Summary. This paper presents a method to compute accurate bounds on Jacobian
determinants of high-order (curvilinear) triangular finite elements. This method canbe used to guarantee that a curvilinear triangle 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 quality the triangles. The key feature of themethod is to expand the Jacobian determinant using a polynomial basis, built usingBézier functions, that has both properties of boundedness and positivity. Numericalresults show the sharpness of our estimates.

Key words: Finite element method; high-order methods; mesh generation;Bézier functions

1 Introduction
High-order finite element solvers critically depend on the availability of highquality curvilinear meshes. Such curvilinear meshes are usually built by curving an initial, straight sided mesh [1, 2, 3]. If we assume that the straightsided mesh is 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.
In this paper we present a method to analyze the quality of curvilinear
meshes in terms of their elementary Jacobian determinants. The method does
not deal with the actual generation of the high order mesh. Instead, it provides
an efficient way to guarantee that each curvilinear element is geometricallyvalid, i.e., that its Jacobian determinant is strictly positive everywhere in its
reference domain. It also provides a way to measure the distortion of thecurvilinear element. The key feature of the method is to adaptively expand

2

A. Johnen, J.-F. Remacle, and C. Geuzaine

the elementary Jacobian determinants in a polynomial basis that has bothproperties of boundedness and positivity. Bézier functions are used to generate these bases in a recursive manner. Contrary to existing approaches usingBézier functions [4, 5], which can only provably detect valid elements (thevalidity of some triangles remaining uncertain), the adaptive scheme allows toprovably detect both valid and invalid elements. The proposed method can beeither used to check the validity and the distortion of an existing curvilinear
mesh, or embedded in the curvilinear mesh generation procedure to assess thevalidity and the quality of the elements on the fly. The algorithm described inthis paper has been implemented in the open source mesh generator Gmsh [6],where it is used in both ways.

2 Curvilinear Meshes, Jacobian Determinants andDistortion
Let us consider a mesh that consists of a set of straight-sided triangles of orderp. Each triangle is defined geometrically through its nodes xi , i = 1, . . . , Np(p)and a set of Lagrange shape functions Li (ξ), i = 1, . . . , Np [7]. The Lagrange
shape functions (of order p) are based on the nodes xi and allow to map areference unit triangle onto the real one:x(ξ) =Np
X(p)

Li (ξ) xi .

(1)

i=1

The mapping x(ξ) should be bijective, which implies that the Jacobian determinant det x,ξ (which is constant over the reference domain since the triangle
is straight-sided) has to be strictly positive. In all what follows we will always
assume that the straight-sided mesh is composed of well-shaped elements, sothat the positivity of det x,ξ is guaranteed. This standard setting is presented
on Figure 1 for the quadratic triangle.
Let us now consider a curved triangle obtained after application of thecurvilinear meshing procedure, i.e., after moving some or all of the nodes ofthe straight-sided triangle. The nodes of the deformed triangle are called Xi ,i = 1 . . . Np , and we haveX(ξ) =Np
X(p)

Li (ξ) Xi .

(2)

i=1

Again, the deformed triangle is valid if the Jacobian determinant J(ξ) :=det X,ξ is strictly positive everywhere over the ξ reference domain. The Jacobian determinant J, however, is not constant anymore over the reference
domain, and computing Jmin := minξ J(ξ) is necessary to ensure 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

Geometrical Validity of High-Order Triangular Finite Elements

3

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

Fig. 1. Reference unit triangle in local coordinates ξ = (ξ, η) and the mappingsx(ξ), X(ξ) and X(x).

expensive and not fully robust since we only get a necessary condition. In thispaper we follow a different approach: because the Jacobian determinant J is apolynomial in ξ, J can be interpolated exactly as a linear combination of specific polynomial basis function over the element. We would then like to obtain
provable bounds on Jmin by using the properties of these basis functions.
In addition, we are also interested in measuring the distortion induced by
the transformation X(x) that maps straight sided elements onto curvilinear
elements. It is possible to write its determinant in terms of the ξ coordinates
as:
det X ,ξJ(ξ)
det X ,x ==.
(3)
det x,ξ
det x,ξ
We call the determinant of the mapping δ(ξ) := det X ,x the distortion and it
should be as close to δ = 1 as possible in order not to degrade the quality ofthe straight sided element. 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.

3 Jacobian Determinant Bounds using Bezier Functions
Computing exact bounds on the Jacobian determinant is only possible forsimple cases. For example, the minimum Jacobian determinant of the planar
2nd order triangle is located either at its unique, easily computable stationary
point, or on the boundary of the triangle.
For arbitrarily high-order curvilinear triangles, computing exact bounds iscomputationally impractical: the Jacobian determinant of a p-th order trian-4

A. Johnen, J.-F. Remacle, and C. Geuzaine

gle is a polynomial function of order at most 2(p − 1). In order to compute
approximate bounds, we propose to expand the Jacobian determinant using apolynomial basis that has both properties of boundedness and positivity—sothat one can verify the positivity of the Jacobian determinant by checking thepositivity of the coefficients in the expansion.
A first candidate for the expansion is the hierarchical basis [8], but quick
analysis shows that the bounds are of poor quality. A better candidate is theBézier basis [9]. Bézier functions are positive and also have the nice property
that they form a partition of unity. As a consequence, the minimum (resp.
maximum) of the coefficients in the expansion is always smaller (resp. greater)
than the Jacobian determinant. In order to compute these coefficients for thehigh-order Lagrange triangles used in standard FEM codes [7], we must firstcompute the change of basis from the Lagrange to the Bézier basis.
3.1 From a Lagrange to a Bezier Basis
The order-p Bézier functions for a triangle in terms of ξ and η are defined ina subspace ET of R2 , with ET = {(ξ, η) : ξ ≥ 0, η ≥ 0, ξ + η ≤ 1} [9]:
 

p p−i
(p)
Ti,j (ξ, η) :=ξ i η j (1 − ξ − η)p−i−j(i + j ≤ p), (4)ija!
is the binomial coefficient. The Jacobian determinant of
where ab = b!(a−b)!
triangles is a polynomial of order at most n = 2(p − 1). Given the expansionof the Jacobian determinant J in terms of Lagrange polynomialsX(n)J(ξ, η) =Ji Li (ξ, η),iwhere the Ji ’s are the values of the Jacobian determinant calculated at Lagrange points, it is possible to express it as an expansion into Bézier triangularpolynomials:X(n)J(ξ, η) =bij Ti,j (ξ, η),
(5)
i+j≤n

where the coefficients bij are the control values of the Bézier expansion.
Thanks to the properties of the Bézier functions, i.e. positivity and partition of unity, the control values bound the Jacobian,min(ξ,η)∈ET

J(ξ, η) ≥ bmin := min bij
i,j

and

max J(ξ, η) ≤ bmax := max bij .

(ξ,η)∈ETi,j(6)
Moreover, there is only one non-zero function at each corner, as for example
at (0, 1)((n)
T0,n (0, 1) = 1(n)Ti,j (0, 1) = 0, ∀ (i, j) 6= (0, n).

Geometrical Validity of High-Order Triangular Finite Elements

5This implies that the 3 control values b00 , b0n and bn0 , are values of theinterpolated function, and we havemin(ξ,η)∈ET

J(ξ, η) ≤ min(b00 , b0n , bn0 ) andmax J(ξ, η) ≥ max(b00 , b0n , bn0 ).

(ξ,η)∈ET(7)
Since Lagrangian and Bézier functions span the same function space, computation of the Bézier values bi from the nodal values Ji (and conversely) is(n)done by a transformation matrix. The transformation matrix T B→L , which
computes nodal values from control values, is created by evaluating Bézierfunctions at sampling points:

 (n)
(n)
T0,0 (0, 0) . . . Tn,0 (0, 0)

 (n)
(n)
 T0,0 (0, 1/n) . . . Tn,0 (0, 1/n) 
(n)
.

(8)
T B→L = 
..
..
..

.


.
.
(n)

(n)

T0,0 (1, 0) . . . Tn,0 (1, 0)
Those sampling points are taken uniformly, i.e., at the location of the nodes(n)(n) −1of the triangle of order n. The inverse transformation T L→B is T B→Lfrom the interpolation of the Jacobian determinant (5), we can write

and

(n)

J = T B→L B,(n)

B = T L→B J ,(9)where B and J are the vectors containing respectively the control values of theJacobian determinant bij and the Ji ’s. For example, for quadratic triangles
we obtain


1
0
0
0 0 0
 0
1
0
0 0 0
001
0 0 0
(2)
.
(10)
T L→B =  −1/2 −1/2 0
2 0 0
 0 −1/2 −1/2 0 2 0 −1/2 0 −1/2 0 0 2
3.2 Adaptive Subdivision
Subdivision consists in interpolating the Jacobian determinant on subdomainsof the triangle. By defining 4 triangular subdomains as shown in Figure 2 (topleft) all the initial domain is recovered, and the subdivision can be carried out
in a simple, recursive manner.
Consider the one-to-one linear map from the reference space to itself thatassociates a point (ξ [q] , η [q] ) ∈ ET [q] of the q th subdomain to every point(ξ, η) ∈ ET of the reference domain, where ET [q] is the subspace of R2 that
corresponds to the q th subdomain. The new interpolation must verify

6

A. Johnen, J.-F. Remacle, and C. Geuzaine

J(ξ [q] , η [q] ) :=X

(n)Ti,j (ξ [q] , η [q] ) bij =

i+j≤n

X

(n)

[q]

(ξ [q] , η [q] ) ∈ ET [q] ,Ti,j (ξ, η) bij

i+j≤n

(11)
[q]
where the coefficients bij are the control values corresponding to the q th subdomain.
Considering the uniform sampling points (ξk , ηk ) as in equation (8) and[q] [q]their image (ξk , ηk ), the expression (11) reads

(n)
(n)
T0,0 [ξ [q] (0, 0), η [q] (0, 0)] . . . Tn,0 [ξ [q] (0, 0), η [q] (0, 0)] (n) [q](n)
 T0,0 [ξ (0, 1/n), η [q] (0, 1/n)] . . . Tn,0 [ξ [q] (0, 1/n), η [q] (0, 1/n)] 
B

..
..
..


.


.
.
(n) [q](n) [q][q][q]
T0,0 [ξ (1, 0), η (1, 0)] . . . Tn,0 [ξ (1, 0), η (1, 0)](n) [q]:= T B→L(n)

B = T B→L B [q] ,(12)

where B [q] is the vector containing the control values of the related subdomain. This implies that we have to compute only one matrix per subdomain:(n)(n) [q]B [q] = T L→B T B→L
B = M[q] B.
(13)Each set of new control values bounds the Jacobian determinant on its own
subdomain and we have:
[q]

[q]

[q]

[q]

b0min := min bij ≤ Jmin ≤ min min(b00 , b0n , bn0 )
q

i,j,q

(14)

and
[q]

[q]

[q]

[q]

max max(b00 , b0n , bn0 ) ≤ Jmax ≤ b0max := max bij .
q

i,q

(15)

If 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 (6)-(7) is
not sharp enough. After one subdivision, the Jacobian determinant is proved
to be positive on all subdomains. In practice, as will be seen in Section 4, a
few levels of refinement lead to the desired accuracy. The convergence of thesubdivision can be proven to be quadratic [10, 11].
3.3 Implementation
As mentioned in Section 2, the Jacobian determinant bounds can be used
to either make the distinction between valid and invalid elements with respect to a condition on Jmin , or to measure the quality of the elements bysystematically 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 =

Geometrical Validity of High-Order Triangular Finite Elements

7

1
0.8

3

0.6
y
0.4

4

0.2

1

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.5

0.5
1
d

0
0

0
0.5

0.5

1

1
j

d

1
j

Fig. 2. Top left: 4 subdomains of the triangle. Top right: third-order planar triangle. Bottom left: exact Jacobian determinant and control values (dots) on theoriginal control points (before subdivision); 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.
(n+1)(n+2), equal to the dimension of the Jacobian determinant space, and so2

to the number of Bézier functions. Second, Bézier values are computed. Then
adaptive subdivision is executed if necessary. Algorithm 1 shows in pseudocode the algorithm used to determine whether the Jacobian determinant of
the element is everywhere positive or not. The algorithm assumes that thecorner points of the triangle are ordered at the 3 first indices.
Algorithm 2 could be further improved by optimizing 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 2). 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 that require subdivision in order to be
detected is usually small. 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 1).
Another possible improvement is to relax the condition of rejection. We
could accept elements for which all control values are positive but reject an

8

A. Johnen, J.-F. Remacle, and C. Geuzaine

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

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

9

if i > Ns then return true;

10

call algorithm 2 with bi as arguments and return its output;

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

for q = 1 to 4 do
i = 1;[q]7
while i ≤ Ns and bi ≥ 0 do
8i = i + 1;

5

6

9
10
11
12

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

element as soon as we find a Jacobian determinant value smaller than a defined
percent of the average Jacobian determinant. The computational gain can besignificant, since elements that were classified as good and that needed a 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 andthe computation of Bézier control values in algorithm 1 can easily be executedfor a whole groups of elements at the same time. This allows to use efficientGeometrical Validity of High-Order Triangular Finite Elements

9

BLAS 3 (matrix-matrix product) functions, which significantly speeds up thecomputations.
The algorithm used for all the tests in the next section is implemented in
the open source mesh generator Gmsh [6] as the AnalyseCurvedMesh plugin.

4 Numerical Results
We start by comparing the new adaptive computation of Jacobian determinantbounds with the brute-force sampling of the Jacobian determinant for thedetection 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. We started thetest for k = 1 and we incremented k until the brute-force approach detected
all the invalid elements. We still executed the algorithm 10 times (while incrementing k) so as to plot the change in the number of invalid element detected.
In order to make the comparison as fair as possible, we have implemented thebrute-force computation as efficiently as possible, i.e., for k (> 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 computethe desired Jacobian determinant values by a matrix-vector product, just likein our own adaptive method.

1Fig.
the elements
elements
Fig.1.3.Two-dimensional
Two-dimensionalmesh
mesh with
with sixth
sixth order
order triangles;
triangles; 23.6%
23.6% ofof theareThearecurved.
curved.The
Thestraightstraightelements
elementsare
are inin blueblue andand thethe invalid
invalid are in dark red. The
other
are
coloredinfunctionoftheminimumof
theirdistortion.
They
are
green
other are colored in function of the minimum of their distortion.
green
if ifthey/J
distortedtheyarearenearly
nearlystraightstraight(J
(Jmin
/Jmax
' 1)
1) andand rather
rather light
light red if really distorted
min
max �
(J(J
/J/J
0).
0).
min
max
min
max�'

We consider the two-dimensional microstructure with circular holes depicted in Figure 3, meshed with 331,050 sixth-order triangles. In this mesh

10

A. Johnen, J.-F. Remacle, and C. Geuzaine

78,180 triangles are curved, and 45,275 are invalid. The new algorithm successfully detects all the 45,275 invalid elements in 6.194s. Some elements needed
as much as 8 levels of subdivisions in order to be classified: see Table 1. The
brute-force approach required 666 sample points per triangle in order to detect
all the invalid elements, and took 4 times longer. But far worse, increasing
the number of sampling points beyond 666 can actually lead to a decreased
accuracy of the prediction, as shown in Figure 4.
Curved Element Classification
Valid curved elements Invalid curved elements
First stage
29303
44967
1 subdivision2436
2 subdivisions
1119
299
3 subdivisions
23
4 subdivisions
10
4
5 subdivisions
92
6 subdivisions
5
7 subdivisions2
8 subdivisions
1
Subtotal
32,90545,275
Total
78,180Undetected invalid elements

Table 1. Number of elements detected as valid or invalid at each stage of theadaptive algorithm; 5 % of the curved elements had to be subdivided adaptively.

4

10

2

10

0

10

1

10

2

10
Sampling points

3

10

Fig. 4. Number of undetected invalid elements using brute-force sampling of theJacobian determinant. The three data points not displayed correspond to the correct
result, i.e., when no invalid triangle is left undetected.

Geometrical Validity of High-Order Triangular Finite Elements

11

Let us now examine the use of the adaptive Jacobian determinant bounds
in the curvilinear meshing algorithms as implemented in Gmsh. We consider
the mesh of a rather coarse version of the world ocean. In our CAD model,
shorelines are described using cubic B-splines: for example, Europe and Asia
are discretized by only one B-spline with about 3,500 control points. The description of this kind of meshing procedure is described in [12]. The quadratic
triangular mesh is generated as follows. We first generate a straight sided
mesh (see Figure 5/(a)). Then, every mesh edge that is classified on a model
edge is curved by snapping its center vertex on the model edge. High ordernodes are then inserted in the middle of every edge that is classified on amodel face (see Figure 5/(b)). This simple procedure does not guarantee thatthe final mesh is valid. In our case, 175 elements are invalid. Then, a global
elasticity analogy is applied to the quadratic mesh that enables to reduce thenumber of invalid elements to 70 (see Figure 5/(c)). Then local optimizations
are performed to remove all invalid elements (see Figure 5/(d)). The final
curvilinear mesh contains about 30% of curved elements. During the meshingprocess, the adaptive Jacobian determinant bound computation allowed todetected all invalid elements. After optimization, the final mesh is composedof elements that have a distortion > 0.1.

(a) Straight sided mesh(b) Raw curvilinear mesh(c) After elastic analogy(d) Final mesh without invalid elements

Fig. 5. The four stages of the curvilinear mesh procedure for the world ocean,meshed with second order triangles.

12

A. Johnen, J.-F. Remacle, and C. Geuzaine

The validity check presented in this article applies directly to triangles
in the plane as shown in the previous examples. When meshing 3D (curved)
surfaces, an additional mapping between the parametric plane where the meshgeneration is performed [6] and the surface has to be taken into account.
This mapping depends on the underlying geometrical representation of thesurface: B-Spline or NURBS-based for CAD models, piece-wise polynomial
for reparametrized STL models [13, 14], ...
In Figure 6 we used the adaptive bounds in the parametric plane for generating a fourth-order mesh of the STEP model of a rotor. After generatinga first order mesh (Figure 6(a)) and snapping the high-order vertices on thegeometrical model (Figure 6(b)), the adaptive Jacobian determinant boundcomputation allowed to pinpoint invalid elements. The final, locally optimized
mesh is displayed in Figure 6(c).
This procedure however does not guarantee the general validity of the 3D
mesh. Indeed, a valid element in the parametric plane can be made invalid
after mapping to the curved surface since only the nodes of the triangle aremapped (and problems can arise between these nodes depending on the actual
mapping). The converse is also true: an invalid element in the parametric plane
can become valid after mapping. This is the subject of ongoing research, as is
the optimization of the quality of general curvilinear finite elements.

5 Conclusion
In this paper we presented a way to compute accurate bounds on Jacobian
determinants of high-order triangular finite elements, based on the efficient
expansion of these Jacobian determinants in terms of Bézier functions. The
proposed algorithm can either be used to determine the validity or invalidityof curved elements, or provide an efficient way to measure their distortion.
Numerical tests show that the method is robust, and a user-defined error
tolerance permits to adjust the accuracy vs. computational time ratio. The
extension to all classical finite element shapes (quadrangles, tetrahedra, prisms
and hexahedra) is currently underway in the open source finite element meshgenerator Gmsh.

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 andcurvilinear mesh generation. Authors also thank V. D. Nguyen for providing
the microstructure geometry used in Figure 3.

Geometrical Validity of High-Order Triangular Finite Elements(a) Straight sided mesh(b) Raw curvilinear mesh(c) After elastic analogy
Fig. 6. Curvilinear mesh of a rotor using fourth-order curved triangles.

13

14

A. Johnen, J.-F. Remacle, and C. Geuzaine

References
1. S. Dey, R. M. O’Bara, and M. S. Shephard. Curvilinear mesh generation in 3D.
Computer Aided Geom. Design, 33:199–209, 2001.
2. M. S. Shephard, J. E. Flaherty, K. E. Jansen, X. Li, X. Luo, N. Chevaugeon,
J.-F. Remacle, M. W. Beall, and R. M. O’Bara. Adaptive mesh generation for
curved domains. Applied Numerical Mathematics, 52:251–271, 2005.
3. S. J. Sherwin and J. Peiró. Mesh generation in curvilinear domains using highorder elements. International Journal for Numerical Methods in Engineering,
53:207–223, 2002.
4. X.-J. Luo, M. S. Shephard, R. M. O’Bara, R. Nastasia, and M. W. Beall. Automatic p-version mesh generation for curved domains. Engineering with Computers, 20:273–285, 2004.
5. O. Sahni, X.J. Luo, K.E. Jansen, and M.S. Shephard. Curved boundary layermeshing for adaptive viscous flow simulations. Finite Elements in Analysis andDesign, 46:132–139, 2010.
6. C. Geuzaine and J.-F. Remacle. Gmsh: a three-dimensional finite element meshgenerator with built-in pre- and post-processing facilities. International Journal
for Numerical Methods in Engineering, 79(11):1309–1331, 2009.
7. T. Hughes. The Finite Element Method. Dover, 2003.
8. I. Babuška, B. Szabò, and R. L. Actis. Hierarchic models for laminated composites. International Journal for Numerical Methods in Engineering, 33:503–535,
1992.
9. G. E. Farin. Curves and surfaces for CAGD: a practicle guide. MorganKaufmann, 2002.
10. J. M. Lane and 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):35–46, 1980.
11. E. Cohen and L. L. Schumacker. Rates of convergence of control polygons.
Computer Aided Geometric Design, 2:229–235, 1985.
12. J. Lambrechts, R. Comblen, V. Legat, C. Geuzaine, and J.-F. Remacle. Multiscale mesh generation on the sphere. Ocean Dynamics, 58:461–473, 2008.
13. J.-F. Remacle, C. Geuzaine, G. Compère, and E. Marchandise. High-quality
surface remeshing using harmonic maps. International Journal for Numerical
Methods in Engineering, 83(4):403–425, 2010.
14. E. Marchandise, C. Carton de Wiart, W. G. Vos, C. Geuzaine, and J.-F.
Remacle. High quality surface remeshing using harmonic maps. Part II: Surfaces
with high genus and of large aspect ratio. International Journal for Numerical
Methods in Engineering, 86(11):1303–1321, 2011.