Meshing, Geometric Modeling and Numerical Simulation, Volume 2: Metrics, Meshes and Mesh Adaptation

By Paul Louis George, Houman Borouchaki, Frederic Alauzet and 3 more

Triangulations, and more precisely meshes, are at the heart of many problems relating to a wide variety of scientific disciplines, and in particular numerical simulations of all kinds of physical phenomena. In numerical simulations, the functional spaces of approximation used to search for solutions are defined from meshes, and in this sense these meshes play a fundamental role. This strong link between meshes and functional spaces leads us to consider advanced simulation methods in which the meshes are adapted to the behaviors of the underlying physical phenomena. This book presents the basic elements of this vision of meshing.

These mesh adaptations are generally governed by a posteriori error estimators representing an increase of the error with respect to a size or metric. Independently of this metric of calculation, compliance with a geometry can also be calculated using a so-called geometric metric. The notion of mesh thus finds its meaning in the metric of its elements.

• Language: English
• Publisher: Wiley
• Release date: Jan 25, 2019
• ISBN: 9781119384366

Book preview

Foreword

The Geometric Modeling and Applications series is made up of five volumes of research on geometric modeling. The act of geometric modeling is an age-old one. To look at a few key points in its history: geometrists in antiquity devoted themselves to it, though, of course, it had a different form from that which we know today; Descartes, in the mid-17th Century, worked on partitioning of planes; and Voronoi established concepts relating to this toward the end of the 19th Century.

But geometric modeling, as we understand it now, can be considered to be linked to the use of computers, especially since this tool has become ubiquitous. We are tempted to fix the scientific recognition of this activity to the date of the First Conference on Computer Aided Geometric Design, which took place in March 1974 at the University of Utah. This was organized by R. Barhnill and R. Riesenfeld, whose presentations may be found in the book Computer Aided Geometric Design, Academic Press, 1974, which clearly associates geometric modeling with the computer.

Since then, many detailed and high-quality books have been published in this field. One of the pioneering works was, undoubtedly, that authored by D. Faux and M. Pratt, Computational Geometry for Design and Manufacture, Ellis Horwood Publishers, 1979. Some books went through several editions, taking note of the advances made in research in the field. In this context, bringing out a new general book on the subject would mark no significant scientific progress.

Moreover, for any object or entity that is to be modeled, there is no single geometric model. Instead, there are several geometric models, each being adapted to the processes to be carried out. For example, the model to finely analyze the shape of an object is necessarily different from models that are used for mechanical computations on this same object. We have, therefore, chosen to focus on specific points of research in fields that we find particularly important. Thus, the five references given here, written by specialists, bring you up to date on advances in research in five domains:

• "Constraints and geometric modeling", by Dominique Michelucci, Pascal Schreck and Pascal Mathis, looks at geometric problems that can be approached by describing the constraints, geometric or more complex, that a solution must verify, rather than directly giving the geometric definition of such a solution.

• "Geometric modeling of fractal forms CAD", by Christian Gentil, Dmitry Sokolov and Gilles Gouaty, discusses the definition of fractal forms: this formalism makes it possible to reuse smooth forms of geometric models that are classic today to define forms that have not yet been explored till present.

• "Meshing, geometric modeling and numerical simulation", in two volumes, by Houman Borouchaki and Paul Louis George, details methodologies advanced for the construction of meshes and geometric modeling for numerial simulation with applications, especially in mechanics.

• "Geometric and Topological Mesh Feature Extraction for 3D Shape Analysis", by Jean-Luc Mari, Franck Hétroy-Wheeler and Gérard Subsol, studies the analysis of three-dimensional-meshed shapes via the extraction of characteristics, both geometric and topological. The fields of application are also detailed to illustrate the transdisciplinary nature of this work.

• "Relevant Triangulations and Samplings for 3D Shapes", by Raphaëlle Chaine and Julie Digne, discusses sampling, meshing and compression of free forms from a large quantity of real data from a virtual, interactive deformation process.

These are books that present recent research. The reader can find, in the general books on geometric modeling that were mentioned before, elements to fill any gaps they may have and to enable them to read the books in this series.

Happy reading.

Marc DANIEL

Introduction

Triangulations and, more precisely, meshes, with the subtle difference that separates these two entities, lie at the heart of any number of problems that arise from varied scientific disciplines. A triangulation or a mesh is a discrete representation, using simple geometric elements (triangle, quadrilateral, tetrahedron, etc., any arbitrary polygon or polyhedron), of a domain that may be an object or a region of which we want a discrete, spatial description. There are, thus, many applications, including numerical simulations of any kind of physical problem, though not restricted to these. In particular, a discrete representation of a (volume) object or a surface may simply be seen as a geometric modeling problem as is. This book adopts a double point of view, as indicated by its title, and we will look both at the use of meshes in numerical simulations (the finite element method, especially) with, of course, the underlying constraints, as well as the use of these meshes for the (discrete) modeling of geometry.

The literature on triangulations and meshes may be classified in two chief categories: one more purely mathematical and geometric, the other more oriented toward industrial applications (numerical simulations) with, of course, though not always, relations between these categories.

The first point of view is covered by the computation geometry community, which studies (among others) Delaunay triangulations in all cuts, definitions, properties, construction algorithms and the complexity of these algorithms. They also study some applications of these triangulations. Nonetheless, relatively recently, mesh generation problems are also being studied, but under a more theoretical angle, generally relating to situations that allow for the use of Delaunay triangulations for which robust construction methods as well as interesting geometric properties have been known for a long time. This (somewhat monotheistic) philosophy necessarily imposes limits on the nature of the problems worked (workable). The first reference book on these subjects is that of Preparata and Shamos [Preparata, Shamos-1985] published in 1985. This was followed by several others, among which we cite two by Edelsbrunner [Edelsbrunner-1987] and [Edelsbrunner-2001], that of Yvinnec and Boissonnat [Boissonnat, Yvinec-1997], by Dey [Dey-2007] and winding up the list with the book by Cheng et al. [Cheng et al. 2012], which was published in 2012. With a few exceptions, the orientation chosen by these references is not always guided by the preoccupations of mathematicians, engineers and the world of numerical simulation in general.

It is this very need for numerical simulations, in particular (and historically) in solid or fluid mechanics, that has led to the emergence of a meshing community among mathematicians and engineers. Without fear of contradiction, we can state that the very first book on meshes that was seen from this point of view is that of Thompson et al. [Thompson et al. 1985], which was also published in 1985. This book essentially discusses structured meshes or grids1 while the first truly generalist book, discussing all types of meshes, structured or not, is that by George [George-1991], which dates back to 1991. Over the years and with the evolution of computers and methods as well as the increasingly complex nature of the meshes to consider, other books have been published. Let us mention books by George and Borouchaki [George, Borouchaki-1998], published in 1997, which revisits meshing methods based on Delaunay algorithms; by Carey [Carey-1997], also in 1997, which gives a very pedagogic presentation of methods; by Topping and co-authors [Topping et al. 2004], in 2004, cited in order to be exhaustive; by Frey and George [Frey, George-2008], in 2000, with a second edition in 2008, which returns to the general view and incorporates newly appeared techniques; by Löhner [Löhner-2008], in 2002, with a second edition that came out in 2008, which abounds in innovative details; up to the recent book by Lo [Lo-2015], published in 2015, which, among others, sheds new light on some problems. In this context, we may wonder what motivated the writing of this book. We will attempt to answer this question (and also try, by doing so, to arouse the reader’s curiosity).

The first observation is that the few books on this subject are already either a little old or rather classic in the way they discuss the approaches and the methods concerning the many questions linked to meshes. The second remark is that, evidently, new questions have appeared relatively recently, e.g. everything related to (finite) elements or patches of higher degrees, metrics and their links to interpolation errors, geometric modeling by meshes, numerical simulation via advanced methods (adaptation of meshes) and so on. This remark immediately calls forth another: we will resolutely adopt a double point of view by considering finite elements (or the constitutive elements of a mesh) as geometric patches and vice versa, thus establishing a link between the Lagrange world of finite elements and the Bézier world of CAD. This choice, as we will see, noticeably changes the manner in which we perceive the problems and, consequently, the manner in which we try to find solutions. In this spirit, this book in two volumes, beyond a set of subjects that we may qualify as classic (and, thus, essential), approaches many subjects that are very recent or even completely novel. Moreover, and this may be more surprising (and certainly rare), some methods (which we find in the literature) have been deliberately left out as they were deemed to not be of great interest and, a contrario, some methods are described in a manner that is, at the very least, nuanced and even critical. Indeed, it seemed pertinent to share our perception, even though this is subjective, up to a point, of these methods, perhaps to prevent the reader from going down paths that we think are dead ends. Having specified this, this book provides the necessary basis for a course on meshes and related problems at the masters level and it may serve as a technical reference for engineering students in science and, more generally, engineers using numerical simulations. We will give a brief description of the content of the chapters in both volumes.

⋄ Volume 1:

We first give a chapter-wise overview of the first volume of this book. In the first three chapters, this volume introduces the basic concepts related to finite elements seen through their shape functions either as finite elements or as patches. In addition to the classic expression via Lagrange polynomials (complete, reduced or rational) (Chapter 1), we give the equivalent formulation based on Bézier forms (Chapter 2), which makes it possible to easily find the conditions for the geometric validity of these finite elements (or patches) whether they are straight or curved and whatever their degree (Chapter 3). Triangulation problems are the subject of the next three chapters, where we specify vocabulary and the basic concepts that make it possible to describe the different construction methods for any triangulation (Chapter 4), Delaunay triangulation (Chapter 5) as well as the triangulation (of a domain) with constraints (Chapter 6). In the next two chapters, we use the concepts introduced earlier to discuss the vast problem of geometric modeling of a domain in its various aspects. Geometric modeling, from our point of views (meshes and numerical simulations), consists of construction of a discrete representation from a continuous representation and vice versa. Different methods are described in Chapters 7, while Chapter 8 gives several significant examples for the application of these methods. Chapter 9 brings together some basic algorithms and formulas related to finite elements (patches) and triangulations, to persuade the reader to go beyond a theoretical viewpoint, to a practical viewpoint, which we believe is essential.

⋄ Volume 2:

We share the contents of Volume 2 of this book, for which four new authors have joined us. This volume (finally) approaches mesh problems and begins with a detailed description of the concept of metric, a concept which will be seen to be fundamental in all that follows. The first two chapters, Chapters 1 and 2, thus focus on metrics. We introduce the concept of metric, we demonstrate their properties and their relations with the geometry of elements of a mesh and how to control interpolation errors during the resolution of a problem (by finite elements). The construction methods of meshes and their optimization are the focus of the following five chapters: Chapter 3 (mesh for a curve), Chapter 4 (simplicial meshes), Chapter 5 (non-simplicial meshes), Chapter 6 (meshes of higher degree) and Chapter 7 (optimizing meshes). We will then discuss the large subject of the adaptation of meshes, controlling the solutions via error estimates and corresponding interpolation metrics in Chapter 8. The use of a dose or a certain dose of parallelism (based on its multiple forms) is seen through Chapter 9. Chapter 10 illustrates, using concrete examples, a series of applications of the methods and methodologies introduced and described through the different chapters of both volumes of this book. As in Volume 1, the final chapter, Chapter 11, gives practical observations on some of the algorithms from Volume 2, especially how to use and manipulate metrics.

⋄ Finally, at the end of this second volume, we sketch out future prospects.

Acknowledgments.

We are especially grateful to two then-doctoral students who willingly helped in developing this book. Nicolas Barral helped us through his expertise with Maple to resolve systems that we encountered during the construction of Serendipity elements. Loïc Frazza, Rémi Feuillet and Julien Vanharen have helped us regarding high-order elements and finite volume methods. Victorien Menier is the author of most of the figures that illustrate this book. We also wish to thank Peha for his magnificent free-hand images, illustrating the impossibility of triangulating certain polyhedra without adding internal vertices. And, of course, we would also like to thank everyone in the common Inria2 and UTT3, Gamma3 team: researchers, engineers and doctoral students, whose work contributed to the development and consolidation of various subjects discussed in this book.

1. The grid-type meshes are generated by solving partial derivative equations, for example elliptical equations. It must be noted that other, more recent books on this theme, followed this pioneering work, but have not been mentioned here.

2. Institut National de Recherche en Informatique et en Automatique.

3. Université de Technologie de Troyes.

Chapter 1

Metrics, Definitions and Properties

The concept of metric has been outlined in Chapter 5 of Volume 1 concerning anisotropic triangulations in which directions are specified as well as lengths according to these directions. It is now essential to clarify what a metric is (a metric field) because metrics are going to be ever-present in the following areas: curve and surface meshing – construction is governed by means of control on curvature(s), meshing in general – construction is governed through the control (at best) of the shape and size of elements and finally, adapted meshing – construction is governed through the control, for example, of interpolation errors resulting in directives concerning shape, directionality and element size. In effect, one will explicitly rely on adequate metrics to orientate the creation or the modification of the meshing involved. It should be noted that we are now referring to meshes, whereas so far we were talking about discretizations in the generic sense of the term or simply about triangulations.

A mesh created ex nihilo or obtained by modification of an existing mesh (remeshing) is characterized by the nature of its elements. The nature of an element, taken individually, is obviously linked to its size and orientation that themselves are completely dependent on the length and orientation of its edges. The nature of a mesh, thereby of a set of elements, depends on its elements; their directions give the directional properties of the mesh seen as a whole and their sizes specifically reflect the density of its points. The most important thing is thus to observe that if edges are controlled, elements are controlled, one by one, and thereby the whole mesh is controlled.

The length of a segment is calculated using the usual dot product, and if u designates the vector corresponding to this segment, we have:

where 〈.,.〉 denotes the dot product. Let O be a given point, the locus of points P such that ||OP||² = 〈OP, OP〉 = 1 is the circle (the sphere) of radius r = 1 centered in O, namely the unit ball (classic). In other words, the distance between O and any one of these points P is the same. The idea of metrics is to change the way lengths are calculated by modifying the definition of the dot product. If is a positive definite symmetric matrix, it can be defined that:

then, O being fixed, the locus of points P such that is the ellipse (ellipsoid) centered at O whose axes have the same directions as those of the eigenvectors of the matrix and whose radii, denoted by ri, are such that , where λi is the ith eigenvalue (normalized via eigenvector normalization) of . Here again, we have the unit ball relative to the underlying metric. The result is that O being fixed, the Euclidean distance from a point P such that at O is different from that of another point satisfying the same relation. This simple observation indicates that the notion of length in a metric contains both length information (Euclidean), but also directional information. This property is the basis itself of all the techniques for the construction of anisotropic meshes (or remeshing) necessary to properly address the problems where direction has a strong impact (curvature direction, directionality of the physical phenomenon under study, etc.). Intuitively, it can be seen that from the matrix (written, see below, as ), associated with the underlying metric¹ , is extracted, which is the transformation that allows shifting from a circle (a sphere) to an ellipse (an ellipsoid). This transformation is interpreted as a change of scale coupled with a rotation (see Figure 1.1). This naive introduction being made, we will introduce and discuss the concept of metric in a formal manner.

First, the definition of a metric is introduced and its properties are indicated. The next step is then to look into ways to combine two or more metrics, either by interpolation or by intersection. The origin of the metrics being diverse, we look into what metrics of geometric origin are (useful for (re)meshing curves and surfaces) and then we show how to find the intrinsic metric of a given element, and, finally, we indicate what calculus metrics are, in the process of solving systems of equations with partial derivatives, obtained via estimators of errors, a topic that will be the subject of Chapter 2. In the end, it is quickly shown how to calculate the length of a segment in the presence of a discrete metric field (that is known only at the vertices of a mesh), and the way to calculate the length of an edge and, finally, we also describe how angles, areas and volumes are to be computed in the presence of a metric field, and all this will be revisited in detail in the last chapter.

1.1. Definitions and properties

A metric represents the unit measure in every direction. It is, therefore, linked to the notion of distance or length. Distances and lengths are defined based on the result of a dot product. We will recall what such a product is before giving a general definition allowing the introduction of the concept of metric.

Let u and v be two vectors (in ℝd with, in our case, d = 2 or 3), the dot product of u and v designated by 〈u, v〉 is a bilinear form, verifying the following three properties:

– [-Pr1] symmetry: 〈u, v〉 = 〈v, u〉;

– [-Pr2] positivity: 〈u, u〉 ≥ 0, for all u =0;

– [-Pr3] definite: 〈u, u〉 = 0 if and only if u = 0.

The second property allows the definition of the norm or vector length u (or still the distance between its two ends) as being . We can then define the unit ball (representing the unit measure in every direction) by ||u|| = 1 with , O fixed and P varying, describing this ball centered in O.

The usual dot product in canonical coordinates is written as:

with ui (respectively, vi) the coordinates of u (respectively, v) within this frame of reference. It is easy to verify that this form satisfies the three properties above (the origin of this dot product dates back to Pythagoras!). It should be noted that the above sum only involves products of coordinates having the same index, and bearing this in mind, the cross-products of the coordinates are not taken into account and, as a result, coordinates of different indices are decorrelated. In this case, a vector has the same norm as any vector corresponding to one of its rotations. This is, of course, natural for the usual canonical space within which distances and lengths are isotropic and the unit ball is a circle (respectively, a sphere) in ℝ² (respectively, ℝ³).

In order to generalize the usual dot product (by taking now the cross-products of the coordinates into account), the general bilinear form is introduced and defined by:

where mij are weighting reals of the coordinate products. By considering the matrix , this generalized dot product is written as:

The matrix must verify the three properties of a dot product, which corresponds to the fact that this matrix has to be symmetric positive definite. The length of a vector u (or its norm) is expressed as . The unit ball with , O fixed and P varying, is an ellipse (respectively, an ellipsoid) in ℝ² (respectively, in ℝ³), centered in O. As a matter of fact, since the matrix is real symmetric, it is diagonalizable and its eigenvectors are orthogonal two by two. Let be the orthogonal matrix whose columns are composed of these vectors; once normalized, it follows that:

with Λ the diagonal matrix composed of the eigenvalues of . Therefore, matrix is written as:

because . We have:

with . As a result, the unit ball can be written in the new basis defined by the columns of as ||v||Λ = 1. Since the coefficients λi of the matrix Λ are strictly positive (because is positive definite), this equation is that of an ellipse (of an ellipsoid) whose main axes are of length . Therefrom, the unit of measure in the direction of the ith eigenvector of (also called ith principal direction of ) is equal to a measure of in the common Euclidean space. The unit of measure in other directions is defined via the elliptical equation ||v||Λ = 1. In the case where λi are identical, the metric is said to be isotropic. In this case, the unit of measure in all directions is equal to the size in the usual metric and the matrix can be written as ( denoting the identity). This relation yields the link with a size in the usual metric and an isotropic metric defined by this size. If λi are different, the metric is said to be anisotropic and, in general, the link between the coefficients of and underlying said to be anisotropic and sizes is not directly visible.

A space equipped with a constant metric (therefore a generalized dot product) is called Euclidean space. The usual Euclidean space is a special case of Euclidian spaces in which the metric is , the identity. The quadratic form associated with , thus , makes it possible to define the unit ball within this space. We have:

with . Thereby, the quadratic form associated with , after a change in variables, can be reduced to the usual quadratic form (that of the usual metric, the identity; Figure 1.1). It can be derived thereof that a Euclidean space endowed with a metric is equivalent to a Euclidean space endowed with the usual identity metric deformed by the linear application:

[1.1]

verifying . This mapping makes it possible, among others, to undo the metric to obtain an identity metric. It is composed of the change in basis followed by the homothety .

This process of deformation of spaces in order to undo a metric is a means to extend the concepts defined in the usual Euclidean space to the Euclidean space equipped with any metric. In particular, in mesh problems, the interest is more specifically in the angle formed by two vectors, in the determinant of a matrix (representing an area or a volume) and in the length of a curved segment.

Let u and v be two vectors of the Euclidean space ℝ² (or ℝ³) equipped with the metric , the angle denoted by (for ) formed by these two vectors is defined as the angle formed by the vectors and , which for their part are defined in the usual space. This angle is thus defined by:

Figure 1.1. The ellipse of the usual Euclidean space (on the left) and the unit ball of the space associated with the metric , on the right. The transformation allowing shifting from one space to another is or according to the direction

or still:

[1.2]

Let two vectors u and v be of the Euclidean space ℝ² (or three vectors u, v and w of ℝ³) equipped with the metric , one considers the second-order matrix (or third-order) whose columns are these vectors. The determinant of this matrix in this space is defined as being the determinant of the matrix formed by the two column vectors and of the space ℝ² (or by the three columns vectors , and of the space ℝ³) equipped with the usual metric. Let:

However²,:

since , thereby, in ℝ²:

[1.3]

Similarly, in ℝ³:

[1.4]

In ℝ² (respectively, ℝ³) equipped with the metric , the quantity (respectively, ) is the area (respectively, the volume) of the parallelogram (respectively, of the parallelepiped) based on u and v (respectively, on u, v and w).

Finally, let Γ be a curved segment defined by the continuous mapping γ of an interval I = [a, b] ⊂ R onto the space ℝ² or ℝ³ equipped with the metric :

The length of Γ is calculated by considering the deformed space associated with the usual metric:

or still:

[1.5]

In the particular case where Γ is a straight segment, defined in the interval I = [0, 1], we have:

which represents the length of the vector in the metric .

The metric was supposed to be constant and has allowed us to define a Euclidean space. A Riemannian space is a space equipped with a variable metric denoted³ by (.), thus dependent on the position. In this case, the scalar product defined by the bilinear form associated with (.) changes from one point to another, implying that every vector must be considered with the metric at its origin. As a result, the dot product of two vectors only makes sense if both vectors have the same origin and is expressed based on the metric at this origin. Under this assumption, the calculations of angles and determinants seen above can be applied. As stated in the Euclidean case, the length of a curved segment Γ depends on the norm in the metric of the vector γ′(t). Since the metric varies from one point to another, this norm also varies from one point to another. It therefore yields:

[1.6]

In the particular case where Γ is a straight segment defined in the interval I = [0, 1] by , one gets:

which does obviously not represent the length of the vector in the metric (A).

The length of a curved segment, even in the case where this segment is a straight segment, involves an elliptic integral that generally has no analytical resolution. The computation of such an integral is therefore based on integral calculus approximation methods as it will be seen at the end of this chapter.

1.2. Metric interpolation and intersection

In this section, the focus is on interpolation problems of a discrete metric field through the meshing of a domain, this being seen as a simplicial covering of the domain. More specifically, the metrics are known at the nodes (which amount to the vertices in the case of a first-degree meshing) of the mesh and the problem consists of defining the metric at any point of the domain. As it will be seen later, this problem occurs both in mesh construction methods as in mesh adaptation methods. In addition, we take a look at metric intersection problems where several discrete metric fields are given, the goal being to define a unique field that best reflects initial fields.

1.2.1. Metric interpolation

A discrete metric field is given over a domain, which is namely known at the nodes of the elements of a mesh of this domain, and different interpolation methods are discussed to find the value of this metric field at any point of the domain. The interpolation on the domain is carried out through its meshing, element by element.

Initially, one considers the case of a mesh with simplicial elements of the first degree. Let K be such an element (a segment, a triangle or a tetrahedron) and Pi its vertices, then any point P of this element is written as , where the αi are the barycentric coordinates⁴ of P. The interpolation problem involves defining the metric (P) at P from the metrics (Pi) at vertices Pi of K and from αi.

The interpolation that we seek, beyond its character of pure interpolation, must verify two paramount properties:

– it must be monotonic from the point of view of the variation in size;

– it must verify a maximum principle from the point of view of the unit balls associated with the metrics.

The monotonicity indicates that the variation in size in each direction is an increasing or decreasing monotonic function. The maximum principle states that, at any P of K, the volume of the unit ball associated with the metric of P is comprised between the minimum and maximum volumes of the balls associated with the vertices of K. It should be noted that the volume of the unit ball associated with the metric (P) is proportional to . Furthermore, from the point of view of each element, if the metric is constant on K, the volume of this element is , where |K| is the usual Euclidean volume of K and the term can be seen as a correction factor for the calculation of a volume in a metric. Therefore, the second property gives control over the change in volume of the elements.

In the isotropic case, the metric at a point P is represented by the matrix , where is the matrix identity and h(P) is the size at P. It is, thus, easy to find the metric at any P of K according to the metrics at its vertices, since it suffices to build a proper interpolation function of sizes. The simplest is a size linear interpolation, that is . Another choice is based on geometric interpolation of size, that is , and the metric in P can be expressed as with either one of these choices for h(P). It should be noted that geometric interpolation gives preference to small sizes, and thereby, meshes that have to satisfy such a metric will be finer than those defined by linear interpolation where, at any point, sizes are larger. Obviously, geometric and linear interpolations verify both monotonicity and maximum properties. Finally, let us mention the existence of other interpolation functions, for example harmonic interpolation defined by .

The anisotropic case is not so straightforward and many approaches of more or less complexity (simplicity) have been proposed and their properties verified. Since we want to see the two aforementioned properties satisfied, we will examine different solutions in light of this double criterion.

⋄ Coefficient interpolation

The most simplistic idea is to interpolate term by term the coefficients of the matrices representing the metrics. In general, we designate by mkl the coefficients of metric . For linear interpolation, we have or even . For geometric interpolation, one would have but since the coefficients mkl can be negative, this interpolation is not valid. Coefficient linear interpolation is generally not satisfactory because the coefficients mkl do not directly reflect the size of the element. However, the first property of monotonicity is verified in every direction and despite the second property not being satisfied, this approach remains a simple possible solution. Moreover, any metric is written as , with a matrix of change in basis and Λ a diagonal matrix whose coefficients are the reciprocals of the squared size, that is to say in the ith direction defined by the change in basis matrix (whose vectors are the columns⁵ of the matrix); termwise interpolation of two metrics and only makes sense if ; therefore, one finds the fact that this approach is legitimate in the isotropic case only.

⋄ Simultaneous reduction interpolation

The idea is to find a common basis (eventually non-orthogonal) in which two given matrices can be written in a diagonal form and, thus, to be able to combine them in this basis using any particular interpolation formula. Before moving on, it should be noted that this approach makes it possible to interpolate in a unique manner the metric on a segment but not in an element as it will be seen further on.

Let two matrix metrics be and , since is invertible, the matrix is well-defined. It is that is symmetric and therefore diagonalizable. Let be the matrix whose columns are the eigenvectors ei (forming a basis of ℝd) of and λi the eigenvalues associated, one thus has . It should be noted that, if there are no multiple eigenvalues, the basis is unique up to a homothety. Otherwise, there are an infinite number of non-homothetic bases for each eigensubspace associated with a multiple eigenvalue. Let X be a vector of ℝd; the coordinates X′ this vector in the basis verify . We have:

Therefore, the matrix of the quadratic form in the basis is . The coefficient ij of this matrix is, therefore, . Similarly, the coefficient ij of the matrix of the quadratic form in the basis is . We are going to show that for i ≠ j, a basis can be chosen (even in the case of multiple eigenvalues) in which coefficients ij are zero for i ≠ j. This means that and are diagonal in the database . Let us consider the expression for i ≠ j; it follows that:

because , we also have:

because is symmetric and is a scalar equal to its transposed. Similarly, since , it follows that:

If λi ≠ λj, one gets , and subsequently, . Otherwise, λi or λj is a multiple eigenvalue and one can choose, in its own subspace (by means of the Schmidt or-thogonalization process), a basis -orthogonal, which thus guarantees that , and subsequently, also . Note that the basis in which and are diagonal must be composed of -orthogonal eigenvectors for the eigensubspaces associated with the multiple eigenvalues. By designating μ1,i (respectively, μ2,i) the quantities (respectively, ), the quadratic form (respectively, ) is written in the diagonal form as Λ1 (respectively, Λ2), a diagonal matrix having μ1,i (respectively, μ2,i) as coefficients. Coefficients μ1,i (respectively, μ2,i) can also be written as (respectively, ), where h1,i (respectively, h2,i) represents the size according to the direction ei.

Therefore, one considers a segment [P1, P2], the metrics at these points, and , and one seeks to define the metric at a point P(t) = (1 – t)P1 + tP2, where the parameter t varies in the range [0, 1], with P(0) = P1 and P(1) = P2. Metrics and , in the basis , are defined by the diagonal matrices Λ1 and Λ2 of coefficients and . The interpolated , in the basis , is defined by a diagonal matrix Λt of coefficients , where ht,i is an interpolation function of h1,i and h2,i. In other words, the interpolated interpolates the sizes following the directions ei.

As already seen, there are several possible choices for the interpolation function. The simplest ones are the following:

[1.7]

[1.8]

Given that the metric is defined in the basis , by the diagonal matrix Λt, its expression in the usual basis is written as:

where –tƤ is an abbreviated notation for , the inverse of the transposed.

In both interpolation cases, since sizes are interpolated only according to d preferred directions (the ei), the monotonicity property cannot be globally verified. The determinant of the matrix is that of the diagonal matrix found above and is therefore equal to . Using a simple counterexample, it is easy to see that the maximum principle is violated in the linear interpolation case (set h1,1 = 1, h2,1 = 3, h1,2 = 3, h2,2 = 1 and evaluate the determinant at ). On the other hand, in the geometric case, one successively gets:

therefore the determinant itself follows a geometric interpolation and, consequently, the maximum principle is verified.

In conclusion, for this method based on simultaneous reduction (with generally non-proportional metrics), linear interpolation, being least expensive, does not verify the maximum principle, whereas geometric interpolation, which is more expensive, has the right properties. In the presence of more than two metrics, it is impossible to build a common basis in which all of these metrics are diagonal. However, since associativity is not guaranteed, by imposing an order of consideration of these metrics, one can gradually build a process of interpolation by taking into account the barycentric coordinates of the point in the element containing it.

⋄ Interpolation using matrix functions

One should bear in mind that interpolation of metric coefficients is not satisfactory but that interpolation on (somehow representing sizes), which is also a metric, leads to considering the mapping of a scalar function f to a matrix, here . Let be a metric represented by its matrix. It can be written as where Λ is the diagonal matrix of the eigenvalues and the basis of the eigenvectors of , each eigenvalue being inversely proportional to the square of a size. The matrix f( ) is defined as .

Let Ƥ be a point of the domain. There is an element K of the mesh of vertices Pi in which it is contained such that with, for all i, 0 ≤ αi ≤ 1. If one considers a linear interpolation on , it then follows that:

On the other hand, considering a geometric interpolation yields the expression:

however, the product of metrics is usually not a metric. It should be preferable to rewrite (Pi) as:

As a result:

or still:

By making use of the exponential and the logarithm [Arsigny et al. 2006], the product of matrices has been transformed into a sum of matrices in which each one is a metric, thereby properly yielding a metric.

As in the case of interpolation based on simultaneous reduction, the monotonicity property is not verified for any of these interpolations; and obviously, linear interpolation does not verify the maximum principle. On the other hand, the latter is properly verified via geometric interpolation. In effect:

where tr(.) refers to the trace of a matrix. Therefore, is a geometric interpolation of .

⋄ Interpolation using natural bases

The previous methods require matrix diagonalizations, and the last one requires costly evaluations of transcendent functions (log and exp). The main idea is to define a common basic form for every metric in which it is diagonal or even identity. In the isotropic case, a metric can be written as , h being the size prescribed by the metric. In the case of an anisotropic metric, this decomposition represents the Cholesky decomposition, namely , where is a lower triangular matrix with strictly positive diagonal coefficients [Laug, Borouchaki-2013], [Laug, Borouchaki-2017]. The fact that matrix acts as suggests that interpolation with respect to should be considered. It therefore yields:

In this case, only the linear interpolation of the coefficients of is possible because they may be negative.

Another possible interpolation can be built based on a factorization of and written as , where is also a lower triangular matrix but with 1 as diagonal coefficients and the coefficients of matrix are inversely proportional to the sizes. Specifically, every metric is written as:

Then, for a linear interpolation, we have:

For geometric interpolation, only the matrix is interpolated assuming a geometric interpolation of the coefficients of (this makes sense because is diagonal and its coefficients are positive). In this case, given that and that is diagonal, this interpolation is a geometric interpolation of the determinants and the maximum principle is verified.

It should be noted that the Cholesky decomposition is by construction dependent in the system of coordinates, whereas the previous method (interpolation of ) is not. Clearly, this means that this method is not really suitable in cases where there is strong anisotropy or fast variation of the metric. On the other hand, this method is less costly than the others.

To go back to the initial point, the construction of a continuous metric field from the metric data at the vertices of a mesh is achieved piece