Tucanos

Tucanos is a 2D and 3D metric-based adaptation library based on Four-Dimensional Anisotropic Mesh Adaptation for Spacetime Numerical Simulations by Philip Claude Caplan.

The public repository containing the sources can be found here.

This document briefly presents the main theoretical aspects and functionnalities that are offered by this remesher.

Mesh

Metric-based mesh adaptation applies only to simplex meshes, i.e. meshes only composes of simplices. In practice, such elements are triangles in 2D and tetrahedras in 3D.

Simplex meshes

Mesh adaptation can only be applied to conformal meshes made only of simplex elements, i.e. triangles in 2D and tetrahedra in 3D. For other types of meshes, the elements must initially be split into simplices: algorithms for standard elements are given in How to Subdivide Pyramids, Prisms and Hexahedra into Tetrahedra, Julien Dompierre Paul Labbé Marie-Gabrielle Vallet Ricardo Camarero and implemented in tucanos.

Mesh

A mesh consists of

• Vertices ${\mathbf{x}}_i\in {\mathbb{R}}^d$
• Elements $K_i$ and element tags $T_{K,i}$
• Tagged faces $F_i$ and face tags $T_{F,i}$

It is characterized by a spatial dimension $d$ as well as a cell dimension $d_{cell}\le d$ (3 for tetrahedra, 2 for triangles, ...). The set of mesh entities of dimension $d$ is denoted $\mathcal{E}(d)$.

A mesh is valid if

• All the elements have a volume $|K_i|>0$,
• All the faces must be either either
  • tagged if
    • it belongs to only one element. In this case, they must be oriented outwards
    • it belongs to two elements with different tags
    • it belongs to more than two elements
  • untagged

The mesh also has to be consistent with the geometry (i.e. a CAD-like model) if all the vertices of the tagged faces lie on the corresponding geometry surfaces.

Topology

This module provides data structures and algorithms to represent, build, and query the topological relationships of a mesh. It extracts hierarchical adjacency information from mesh elements (e.g., tetrahedrons, triangles) and faces, mapping out how entities of different dimensions (volumes, faces, edges, vertices) connect to one another.

Core Concepts

• Dim: Represents the dimension of an entity (e.g., 0 for vertices, 1 for edges, 2 for faces, 3 for volumes).
• Tag: A unique integer identifier for a specific entity within a given dimension.
• TopoTag: A tuple of (Dim, Tag) uniquely identifying any topological entity across the entire mesh.
• Hierarchy: In this system, "parents" are strictly higher-dimensional entities that bound a "child" entity. For example, a 3D Volume (Dim 3) is a parent to its bounding 2D Faces (Dim 2). A 2D Face is a parent to its 1D Edges.

Core Data Structures

TopoNode

Represents a single topological entity in the mesh and its immediate relationships.

pub struct TopoNode {
    pub tag: TopoTag,            // The unique identifier (Dim, Tag)
    pub children: HashSet<Tag>,  // Tags of lower-dimensional boundary entities (Dim - 1)
    pub parents: HashSet<Tag>,   // Tags of higher-dimensional bounding entities (Dim + 1)
}

Topology

The central registry containing all TopoNodes, categorized by dimension, and caching their hierarchical relationships.

pub struct Topology {
    dim: Dim,
    entities: Vec<Vec<TopoNode>>, 
    parents: FxHashMap<(TopoTag, TopoTag), TopoTag>,
}

• entities: A 2D vector where entities[d] holds all the TopoNodes of dimension d.
• parents: A cached mapping used to quickly find the closest common parent of a given pair of TopoTags.

MeshTopology

A high-level wrapper that binds a Topology to a specific Mesh instance. It manages the Topology and stores the resolved topological tags for all vertices (vtags).

pub struct MeshTopology {
    topo: Topology,
    vtags: Vec<TopoTag>,
}

API Reference: Topology

Initialization and Serialization

• new(dim: Dim) -> Self: Initializes an empty topology supporting up to the specified dimension.
• from_json(fname: &str) -> Result<Self>: Deserializes a topology from a JSON file and automatically recomputes the parent caches.
• to_json(&self, fname: &str) -> Result<()>: Serializes the topology to a pretty-printed JSON file.

Querying Tags and Entities

• ntags(&self, dim: Dim) -> usize: Returns the number of tags registered in a specific dimension.
• tags(&self, dim: Dim) -> Vec<Tag>: Returns a list of all tags present in a specific dimension.
• first_available_tag(&self, dim: Dim) -> Tag: Calculates the next available unique tag for a given dimension (the maximum absolute tag value + 1).
• get(&self, tag: TopoTag) -> Option<&TopoNode>: Retrieves a reference to a TopoNode by its TopoTag.

Hierarchy Construction

• insert(&mut self, tag: TopoTag, parents: &[Tag]): Manually inserts a new entity. It automatically updates the children set of the specified parent entities (dimension tag.0 + 1).
• get_from_parents(&self, dim: Dim, parents: &[Tag]) -> Option<&TopoNode>: Finds a node in the specified dimension that has the exact set of given parent tags.
• parent(&self, topo0: TopoTag, topo1: TopoTag) -> Option<TopoTag>: Retrieves the closest common parent of two topology entities, leveraging the internal FxHashMap cache.

Automatic Topology Inference

The module contains robust algorithms to automatically deduce lower-dimensional topology from a mesh's elements and boundary faces:

• update_from_elems_and_faces: The primary engine used by MeshTopology. It iterates over the mesh elements and faces, infers the tags for intermediate dimensions (like edges in a 3D mesh), resolves tagging conflicts at corners/junctions (check_and_fix), and returns a fully resolved vector of vertex tags (vtags).
• clear<F>(&mut self, filter: F): Filters out topology nodes based on a closure, cleaning up dangling child/parent relationships along the way.

API Reference: MeshTopology

• new<const D: usize, M: Mesh<D>>(msh: &M) -> Self: Constructs a full topological mapping from a generic mesh M. It analyzes element tags (etags) and face tags (ftags), generating unique tags for unmarked internal boundaries, edges, and vertices.
• topo(&self) -> &Topology: Grants immutable access to the underlying Topology.
• vtags(&self) -> &[TopoTag]: Returns a slice containing the computed TopoTag for every vertex in the underlying mesh, indexed by the vertex ID.

Usage Example

While usually constructed automatically via MeshTopology::new(&mesh), the topology can be constructed manually (which is often useful for testing and custom mesh generation):

let mut t = Topology::new(3);

// Insert a 3D volume (Dim 3, Tag 1) with no parents
t.insert((3, 1), &[]);

// Insert 2D faces (Dim 2) whose parent is the volume (Tag 1)
t.insert((2, 1), &[1]);
t.insert((2, 2), &[1]);

// Insert a 1D edge (Dim 1) bounded by both faces
t.insert((1, 1), &[1, 2]); // Edge shared by Face 1 and Face 2

// Automatically resolve all parent-child hierarchy caches
t.compute_parents();

// Fetch an edge bounded by faces 1 and 2
let edge = t.get_from_parents(1, &[1, 2]).unwrap();
assert_eq!(edge.tag.1, 1);

Geometry

This module provides the structural foundation for associating computational meshes with their underlying geometric domains. It facilitates projecting mesh vertices onto true geometric surfaces (e.g., during mesh adaptation or order elevation), computing geometric deviations (like distance and normal angles), and handling surface curvature.

Core Trait: Geometry<const D: usize>

The Geometry trait defines the essential interface for any geometric representation in a D-dimensional space. It is designed to be thread-safe (Send + Sync).

Key Methods

• check(&self, topo: &Topology) -> Result<()> Validates that the geometric model is consistent with the provided topological hierarchy.
• project(&self, pt: &mut Vertex<D>, tag: &TopoTag) -> f64 Projects a given vertex pt onto the geometric entity identified by tag. It modifies the vertex coordinates in-place to the projected coordinates and returns the projection distance.
• angle(&self, pt: &Vertex<D>, n: &Vertex<D>, tag: &TopoTag) -> f64 Computes the angle (in degrees) between a given normal vector n and the true geometric normal at the projection of pt on the specified geometric surface.
• project_vertices<M>(&self, mesh: &mut M, topo: &MeshTopology) -> f64 Iterates through all vertices in a mesh and projects boundary vertices (entities with a dimension < D) onto the geometry, returning the maximum projection distance.
• max_distance<M>(&self, mesh: &M) -> f64 Computes the maximum geometric deviation (distance) between the mesh's face/element centers and the true geometry.
• max_normal_angle<M>(&self, mesh: &M) -> f64 Computes the maximum angle deviation between the mesh's discrete normals and the true continuous geometry normals.
• to_quadratic_triangle_mesh & to_quadratic_edge_mesh Elevates a linear mesh (straight edges/flat faces) to a quadratic mesh. It achieves this by inserting mid-edge nodes and projecting those new nodes onto the true geometry using project_vertices, thereby curving the mesh to match the physical boundaries.

Implementations

1. NoGeometry

A dummy implementation representing the absence of an underlying geometric model.

• project() leaves the vertex unchanged and returns 0.0.
• angle() always returns 0.0.

2. MeshedGeometry<const D: usize, M: Mesh<D>>

A discrete, STL-like representation of a geometry constructed from a high-resolution boundary mesh. It relies on spatial indices (ObjectIndex, a bounding volume hierarchy) to efficiently find the closest projection points on the discrete surface.

Internal Patch Structures

To manage different topological dimensions, the geometry is split into discrete patches:

• MeshedPatchGeometry: Represents a specific topological edge/curve (dimension D-1). It holds the elements belonging to a specific topological tag and uses a spatial index to compute projections.
• MeshedPatchGeometryWithCurvature: Represents a specific topological surface patch (dimension D). In addition to projecting points, it pre-computes and stores principal curvature directions (u and optionally v for 3D flows) for anisotropic refinement.

Attributes of MeshedGeometry

• patches: A hash map of surface patches, keyed by their topological tag.
• edges: A hash map of bounding edges/curves, keyed by their topological tag.
• edge2faces & edge_map: Topological mapping structures. edge_map maps the tags of a working simulation mesh back to the original baseline tags of this geometric representation (crucial if the simulation mesh undergoes splitting or re-tagging).

Workflow & Usage

1. Initialization: Instantiated via MeshedGeometry::new(&boundary_mesh). The provided boundary mesh must have correctly tagged faces and edges (e.g., using mesh.fix()).
2. Topology Mapping: Use set_topo_map(&mesh_topo) to reconcile edge tags if the simulation mesh's topology differs slightly from the baseline geometry.
3. Projection: When project() is invoked, the geometry inspects the dimension of the TopoTag to route the point:
   • Dim == D: Projects the point onto the corresponding 2D surface patch.
   • Dim == D - 1: Projects the point onto the corresponding 1D boundary curve/edge.
   • Dim == 0: Represents a hard corner/vertex; the point remains fixed in space.

Utility Functions

• curvature(&self, pt: &Vertex<D>, tag: Tag) Extracts the principal curvature directions at a specific vertex on a surface patch.
• write_curvature(&self, fname: &str) Exports the computed curvature vector fields as .vtu files (appending the topological tag to the filename). This is extremely useful for verifying the curvature field visually in tools like ParaView.

Metric

Metric and Euclidean metric space

An Euclidean metric space $({\mathbb{R}}^d,\mathcal{M})$ is a finite-dimensional vector space where the dot product is defined by means of a symmetric definite positive matrix $\mathcal{M}$:

$$⟨\mathbf{u},\mathbf{v}{⟩}_{\mathcal{M}}=⟨\mathbf{u},\mathcal{M}\mathbf{v}⟩={\mathbf{u}}^{\top }\mathcal{M}\mathbf{v},\text{ for }(\mathbf{u},\mathbf{v})\in {\mathbb{R}}^d\times {\mathbb{R}}^d$$

The properties of $\mathcal{M}$ ensure that it defines a dot product. We will now refer to $\mathcal{M}$ as a metric tensor or metric.

Canonical Euclidean metric space $({\mathbb{R}}^3,{\mathcal{I}}_3)$

The most familiar example of Euclidean space is the one defined by an identity matrix. The 2D canonical one is defined by the 2x2 identity matrix ${\mathcal{I}}_2$ and the 3D by the 3x3 identity matrix ${\mathcal{I}}_3$.

Geometric definitions in Euclidean spaces

The dot product defined by $\mathcal{M}$ allows to define the distance notion in ${\mathbb{R}}^d$. In the following $d$ denotes the dimension and is equal either to 2 or 3. We can then define norm and distance definitions:

• $\forall \mathbf{u}\in {\mathbb{R}}^d,||\mathbf{u}|{|}_{\mathcal{M}}=\sqrt{⟨\mathbf{u},\mathcal{M}\mathbf{u}⟩}$
• $\forall (\mathbf{a},\mathbf{b})\in {\mathbb{R}}^d\times {\mathbb{R}}^d,d_{\mathcal{M}}(\mathbf{a},\mathbf{b})=||\mathbf{b}-\mathbf{a}|{|}_{\mathcal{M}}=||\mathbf{ab}|{|}_{\mathcal{M}}$

The length $l_{\mathcal{M}}$ of a segment $\mathbf{ab}=[\mathbf{a},\mathbf{b}]$ is then given by the distance between its extremities:

$$l_{\mathcal{M}}(\mathbf{ab})=d_{\mathcal{M}}(\mathbf{a},\mathbf{b})=||\mathbf{ab}|{|}_{\mathcal{M}}$$

Moreover, given a bounded subset $K$ of ${\mathbb{R}}^d$, the volume $|K{|}_{\mathcal{M}}$ is given by:

$$|K{|}_{\mathcal{M}}={\int }_K\sqrt{\det \mathcal{M}}d\mathbf{x}=\sqrt{\det \mathcal{M}}|K{|}_{{\mathcal{I}}_d}$$

where $|K{|}_{{\mathcal{I}}_d}$ is the Euclidean volume of $K$.

Finally, as metric tensor $\mathcal{M}$ is symmetric, it is diagonalizable. It thus admits the following spectral decomposition:

$$\mathcal{M}=\mathcal{R}\Lambda {\mathcal{R}}^{\top }$$

where $\mathcal{R}$ is an orthonormal matrix composed of the eigenvectors $({\mathbf{v}}_{\mathbf{i}}{)}_{i=1,d}$ of $\mathcal{M}$:

$$\mathcal{R}=({\mathbf{v}}_1\dots {\mathbf{v}}_{\mathbf{d}})$$

verifying $\mathcal{R}{\mathcal{R}}^{\top }={\mathcal{R}}^{\top }\mathcal{R}={\mathcal{I}}_d$. $\Lambda =\text{diag}({\lambda }_i)$ is a diagonal matrix composed of the eigenvalues of $\mathcal{M}$, denoted $({\lambda }_i{)}_{i=1,d}$, which are strictly positive.

Natural metric mapping

The metric $\mathcal{M}$ can be seen as a mapping from the unit ball ${\mathcal{B}}_{{\mathcal{I}}_d}$ of identity matrix ${\mathcal{I}}_d$ onto the unit ball ${\mathcal{B}}_{\mathcal{M}}$ of $\mathcal{M}$, where these unit balls are defined as:

$$\begin{array}{rl}{\mathcal{B}}_{{\mathcal{I}}_d}& =\left\{\mathbf{x}\in \mathcal{V}(\mathbf{a})\mid (\mathbf{x}-\mathbf{a}{)}^{\top }{\mathcal{I}}_d(\mathbf{x}-\mathbf{a})=1\right\}\\ {\mathcal{B}}_{\mathcal{M}}& =\{\mathbf{x}\in \mathcal{V}(\mathbf{a})\mid (\mathbf{x}-\mathbf{a}{)}^{\top }\mathcal{M}(\mathbf{x}-\mathbf{a})=1\}\end{array}$$

Where $\mathcal{V}(\mathbf{a})$ denotes the vicinity of point $\mathbf{a}$. A visual representation of these unit balls is given in Figure 1. Notice that ${\mathcal{B}}_{{\mathcal{I}}_d}$ defines the unit circle in 2D and the unit sphere in 3D.

Figure 1: Geometric interpretation of the unit ball ${\mathcal{B}}_{\mathcal{M}}$ of $\mathcal{M}$. ${\mathbf{v}}_i$ are the eigenvectors of $\mathcal{M}$ and ${\lambda }_i=h_i^{-2}$ are the eigenvalues of $\mathcal{M}$. (Loseille 2011)

${\mathcal{M}}^{-\frac{1}{2}}$ is defined by the spectral decomposition:

$${\mathcal{M}}^{-\frac{1}{2}}=\mathcal{R}{\Lambda }^{-\frac{1}{2}}{\mathcal{R}}^{\top }\text{where}{\Lambda }^{-\frac{1}{2}}=\text{diag}({\lambda }_i^{-\frac{1}{2}})$$

Where $\text{diag}$ denotes a 2D or 3D diagonal matrix. ${\Lambda }^{-\frac{1}{2}}$ defines the diagonal matrix of the sizes $h_i$ along directions as $h_i={\lambda }_i^{-\frac{1}{2}}$. This gives the relation between the metric $\mathcal{M}$ and the size it enforces.

The mapping defined above is written in the form of the following application:

$$\begin{array}{rl}{\mathcal{M}}^{-\frac{1}{2}}:{\mathbb{R}}^d& ⟼{\mathbb{R}}^d\\ \mathbf{x}& ⟶{\mathcal{M}}^{-\frac{1}{2}}\mathbf{x}\end{array}$$

Figure 2: Natural mapping associated with metric $\mathcal{M}$. (Loseille 2008)

for more details, see A. Loseille and F. Alauzet, Continuous mesh framework. Part I: well-posed continuous interpolation error, SIAM J. Numer. Anal., Vol. 49, Issue 1, pp. 38-60, 2011. PDF

Riemannian metrics

Riemannian metric space

A Riemannian metric space $({\mathcal{R}}^d,\mathcal{M}(\mathbf{x}))$ is defined by a continuous metric field $\mathcal{M}(\cdot )$ in ${\mathcal{R}}^d$. The length of an edge $\mathbf{ab}$ is computed using the path $\gamma (t)=\mathbf{a}+t\mathbf{ab}$ where $t\in [0,1]$:

$$l_{\mathcal{M}}(\mathbf{ab})={\int }_0^1||{\gamma }^{{}^{\prime }}(t)||dt={\int }_0^1\sqrt{{\mathbf{ab}}^{\top }\mathcal{M}(\mathbf{a}+t\mathbf{ab})\mathbf{ab}}dt$$

The figure below shows length computation in different metric spaces. We see that the iso-values of distance in the Riemannian space are "distorted" with the loss of symmetry when compared to the Euclidean spaces.

Figure 1: Iso-values of the distance w.r.t canonical Euclidean space (left), euclidean space defined by $\mathcal{M}$ and Riemannian space defined by $\mathcal{M}(\cdot )$ (Loseille 2008)

Given a bounded subset $K$ of ${\mathbb{R}}^d$, the volume $|K{|}_{\mathcal{M}(\cdot )}$ is given by:

$$|K{|}_{\mathcal{M}(\cdot )}={\int }_K\sqrt{\text{det}\mathcal{M}(\mathbf{x})}d\mathbf{x}$$

Unit elements and meshes

Let $K$ be a simplex element in $\Omega \subset {\mathbb{R}}^n$ where $n$ denotes its dimension. $K$ is unit w.r.t the euclidean metric field $\mathcal{M}$ when all its edges are equal to 1 in this metric. Mathematically speaking, if we define the list of edges of $K$ by $({\mathbf{e}}_i{)}_{i=1,n(n+1)/2}$, then

$$\forall i=1,...,\frac{n(n+1)}{2},l_{\mathcal{M}}({\mathbf{e}}_i)=1$$

And, the volume of $K$ is given by:

$$|K{|}_{\mathcal{M}}=\frac{\sqrt{n+1}}{2^{n/2}n!}\text{ and }|K{|}_{{\mathcal{I}}_n}=\frac{\sqrt{n+1}}{2^{n/2}n!}\text{det}({\mathcal{M}}^{-\frac{1}{2}})$$

Figure 2: Geometric representation of unit tetrahedrons in metric ${\mathcal{I}}_3$ (left) and an anisotropic metric $\mathcal{M}$ (right) (Loseille 2008)

The metric $\mathcal{M}$ for which $K$ is unit is unique and referred to as the implied metric of $K$. Conversely, when given a metric $\mathcal{M}$, the set of unit elements is not empty. Similarly, the unit notion can be extended to a mesh composed of simplexes with a continuous metric field $\mathcal{M}(\cdot )$. The mesh could then be said unit, if all its simplex elements are unit in $\mathcal{M}(\cdot )$. However, it is well known that the equilateral tetrahedron does not fill a subset $\Omega$ of ${\mathbb{R}}^3$. Following the proposition in the literature (Loseille 2011), the unit notion for simplex elements is relaxed to the quasi-unit element notion. $K$ is said to be quasi-unit if:

$$\forall i=1,...,\frac{n(n+1)}{2},l_{\mathcal{M}}({\mathbf{e}}_i)\in \left[\frac{1}{\sqrt{2}},\sqrt{2}\right]$$

Yet, some examples show that elements with a null volume (i.e. degenerated elements) can satisfy the above property. The authors propose a strategy to control the ratio between the volume $|K{|}_{\mathcal{M}}$ and the sum of the square length of the edges $l_{\mathcal{M}}^2({\mathbf{e}}_i)$. This ratio defines the quality function ${\mathcal{Q}}_{\mathcal{M}}$ that measures the proximity of $K$ to the perfectly unit element in $\mathcal{M}$:

$${\mathcal{Q}}_{\mathcal{M}}(K)=\frac{n(n!{)}^{2/n}}{(n+1{)}^{\frac{1-n}{n}}}\frac{|K{|}_{\mathcal{M}}^{\frac{2}{n}}}{{\sum }_{i=1}^{\frac{n(n+1)}{2}}{l}_{\mathcal{M}}^2({\mathbf{e}}_i)}\in [0,1]$$

The constant $\frac{n(n!{)}^{2/n}}{(n+1{)}^{\frac{1-n}{n}}}$ is chosen so that ${\mathcal{Q}}_{\mathcal{M}}$ is equal to 1 when $K$ is perfectly unit in $\mathcal{M}$. For a null volume ${\mathcal{Q}}_{\mathcal{M}}$ is 0. A constraint on the quality is added to the definition of quasi-unit element to avoid describing degenerated elements.

Then, $K$ is said to be quasi-unit if:

$$\forall i=1,...,\frac{n(n+1)}{2},l_{\mathcal{M}}({\mathbf{e}}_i)\in \left[\frac{1}{\sqrt{2}},\sqrt{2}\right]\text{ and }{\mathcal{Q}}_{\mathcal{M}}(K)\in [\alpha ,1]\text{ with }\alpha >0$$

Finally, the simplex mesh ${\mathcal{T}}_h$ is said unit with respect to the metric field $(\mathcal{M}(\mathbf{x}){)}_{\mathbf{x}\in \Omega }$ when all its elements are quasi-unit in $(\mathcal{M}(\mathbf{x}){)}_{\mathbf{x}\in \Omega }$. This defines the duality between the discrete and continuous mesh, in the sense that, given a metric field there exists a unit mesh with respect to that metric field $(\mathcal{M}(\mathbf{x}){)}_{\mathbf{x}\in \Omega }$.

Complexity

We define the complexity $\mathcal{C}$ of metric field $(\mathcal{M}(\mathbf{x}){)}_{\mathbf{x}\in \Omega }$:

$$\mathcal{C}(\mathcal{M}(\cdot ))={\int }_{\Omega }\sqrt{\text{det}(\mathcal{M}(\mathbf{x}))}d\mathbf{x}$$

This real value parameter is useful to quantify the global level of accuracy of $(\mathcal{M}(\mathbf{x}){)}_{\mathbf{x}\in \Omega }$. It can also be interpreted as the continuous counterpart of the number of vertices of a discrete mesh. There is a direct relationship between the prescribed metric complexity and the number of elements of the generated mesh. Let ${\mathcal{T}}_h$ be a simplex mesh of $\Omega \subset {\mathbb{R}}^d$ composed of simplex elements $K$. We can describe this mesh with the metric $\mathcal{M}(\mathbf{x}{)}_{\mathbf{x}\in \Omega }$ in which ${\mathcal{T}}_h$ is unit. Then, the continuous mesh complexity and the mesh number of element $N_{elts}$ are linked using the unit element properties:

$$\mathcal{C}\approx \sum _K\sqrt{\text{det}\mathcal{M}}|K|\approx \sum _K\frac{\sqrt{d+1}}{2^{d/2}d!}=\frac{\sqrt{d+1}}{2^{d/2}d!}{N}_{elts}$$

After remeshing, this relation is not exact because most of the elements are not perfectly unit.

Embedded Riemannian spaces

Two embedded Riemannian spaces have the same anisotropic ratios and orientations. They differ only from their complexity. Mathematically speaking, two Riemannian spaces $(\mathcal{M}(\mathbf{x}){)}_{\mathbf{x}\in \Omega }$ and $(\mathcal{N}(\mathbf{x}){)}_{\mathbf{x}\in \Omega }$ are embedded if a constant $c$ exists such that:

$$\forall \mathbf{x}\in \Omega ,\mathcal{N}(\mathbf{x})=c\mathcal{M}(\mathbf{x})$$

Conversely, from $\mathbf{M}=(\mathcal{M}(\mathbf{x}){)}_{\mathbf{x}\in \Omega }$ we can deduce $\mathbf{N}=(\mathcal{N}(\mathbf{x}){)}_{\mathbf{x}\in \Omega }$ of complexity $N$ having the same anisotropic properties (anisotropic orientations and ratios) by considering:

$$\mathcal{N}(\mathbf{x})={\left(\frac{N}{\mathcal{C}(\mathbf{M})}\right)}^{\frac{2}{3}}\mathcal{M}(\mathbf{x})$$

In the context of error estimation, this notion enables convergence order study with respect to an increasing complexity $N$. Consequently, the complexity $\mathcal{C}(\mathbf{M})$ is also the continuous counterpart of the classical parameter $h$ used for uniform meshes while studying convergence.

Figure 3: Different unit elements where only the density increases from left to right. (Loseille 2008)

Local decomposition of the metric field in ${\mathbb{R}}^n$

Let $\mathcal{M}(\mathbf{x})$ be a metric defined locally on point $\mathbf{x}$ of domain $\Omega \subset {\mathbb{R}}^n$. Then, by definition, this metric can be written in diagonal form yielding:

$$\begin{array}{rl}\mathcal{M}(\mathbf{x})& ={\mathcal{R}}^{\top }(\mathbf{x})\left(\begin{array}{ccc}{\lambda }_1(\mathbf{x})& & \\ & \ddots & \\ & & {\lambda }_d(\mathbf{x})\end{array}\right)\mathcal{R}(\mathbf{x})\\ & ={\mathcal{R}}^{\top }(\mathbf{x})\left(\begin{array}{ccc}{h}_1^{-2}(\mathbf{x})& & \\ & \ddots & \\ & & h_n^{-2}(\mathbf{x})\end{array}\right)\mathcal{R}(\mathbf{x})\end{array}$$

In this form, the optimal orientations are given by the rotation matrix $\mathcal{R}(\mathbf{x})$ and the sizes are given by $h_i(\mathbf{x})$, given by the eigenvalues ${\lambda }_i^{-\frac{1}{2}}(\mathbf{x})$. However, we want to write a form that locally decouples the anisotropic information from the global size. Following (Loseille 2011b), the following local parameters can be defined:

• Anisotropic quotients $r_i$ defined as:

$$r_i=h_i{\left(\prod _{k=1}^n{h}_k\right)}^{-\frac{1}{n}}$$

• A density $d$ defined as:

$$d={\left(\prod _{k=1}^n{\lambda }_k\right)}^{\frac{1}{2}}={\left(\prod _{k=1}^n{h}_k\right)}^{-1}$$

Locally, the metric writes:

$$\mathcal{M}(\mathbf{x})=d^{\frac{2}{n}}(\mathbf{x}){\mathcal{R}}^{\top }(\mathbf{x})\left(\begin{array}{ccc}{r}_1^{-2}(\mathbf{x})& & \\ & \ddots & \\ & & r_n^{-2}(\mathbf{x})\end{array}\right)\mathcal{R}(\mathbf{x})$$

This decomposition is used in the Hessian scaling technique described in the section on Hessian Scaling.

Geometric invariants for unit elements

Unit element's geometric properties can be connected to the linear algebra properties of metric tensors thanks to geometric invariants:

• Invariant related to the Euclidean volume $|K|$: $$|K|=\frac{\sqrt{3}}{4}\text{det}({\mathcal{M}}^{-\frac{1}{2}})\text{ in 2D and }|K|=\frac{\sqrt{2}}{12}\text{det}({\mathcal{M}}^{-\frac{1}{2}})\text{ in 3D.}$$

• Invariant related to the square length of the edges for all symmetric matrix $H$:

  $$\begin{array}{rl}\sum _{i=1}^3{\mathbf{e}}_i^{\top }H{\mathbf{e}}_i& =\frac{3}{2}\text{tr}({\mathcal{M}}^{-\frac{1}{2}}H{\mathcal{M}}^{-\frac{1}{2}})\text{ in 2D,}\\ \sum _{i=1}^6{\mathbf{e}}_i^{\top }H{\mathbf{e}}_i& =2\text{tr}({\mathcal{M}}^{-\frac{1}{2}}H{\mathcal{M}}^{-\frac{1}{2}})\text{ in 3D,}\end{array}$$

Element-implied metric computation

Practically we compute the implied metric ${\mathcal{M}}_K$ by first computing ${\mathcal{M}}_K^{-\frac{1}{2}}$ from a two-step transformation. Let ${\phi }_0$ be the jacobian transformation matrix from $K_{\perp }$ to $K_{\Delta }$ and $\phi$ be the jacobian transformation matrix from $K_{\perp }$ to $K$.

We use the triangle mapping definition to compute ${\phi }_0$ and $\phi$. The element-implied metric becomes:

$${\mathcal{M}}_K=(\mathbf{\Phi }{\mathbf{\Phi }}^{\top }{)}^{-1}=(\phi {\phi }_0^{-1}{\phi }_0^{-\top }{\phi }^{\top }{)}^{-1}$$

Figure 4: Relationship between the physical $K$ and reference elements ($K_{\Delta }$ and $K_{\perp }$ ) for the calculation of the element-implied metric. (CAPLAN2020102915)

Numerical computation of lengths

Let $\mathbf{a}$ and $\mathbf{b}$ be two points of the domain. The length of segment $\mathbf{ab}$ can be numerically approximated using a $k$-points Gaussian quadrature with weights $({\omega }_i{)}_{i=\mathrm{1...}k}$ and barycentric coefficients $({\alpha }_i{)}_{i=\mathrm{1...}k}$:

$$l_{\mathcal{M}}(\mathbf{ab})={\int }_0^1\sqrt{{\mathbf{ab}}^{\top }\mathcal{M}(\mathbf{a}+t\mathbf{ab})\mathbf{ab}}dt\approx \sum _{i=1}^k{\omega }_i\sqrt{{\mathbf{ab}}^{\top }\mathcal{M}(\mathbf{a}+{\alpha }_i\mathbf{ab})\mathbf{ab}}$$

where $\mathcal{M}(\mathbf{a}+{\alpha }_i\mathbf{ab})$ is the metric at the $i^{th}$ Gauss point. In the context of meshing, we are interested in the evaluation of edge length in the metric. This is a recurrent operation which could be time consuming for the mesh generator. Fortunately, this operation is at the edge level where the metric is known at the endpoints. The problem can thus be simplified by considering a variation law of the metric along the edge and computing analytically the edge length in the metric.

Formally speaking, let $\mathbf{e}={\mathbf{p}}_1{\mathbf{p}}_2$ be an edge of the mesh of Euclidean length $||\mathbf{e}|{|}_2$, and $\mathcal{M}({\mathbf{p}}_1)$ and $\mathcal{M}({\mathbf{p}}_2)$ be the metrics at the edge extremities ${\mathbf{p}}_1$ and ${\mathbf{p}}_2$. We denote by $l_i(\mathbf{e})=\sqrt{{\mathbf{e}}^{\top }\mathcal{M}({\mathbf{p}}_{\mathbf{i}})\mathbf{e}}$ the length of the edge in metric $\mathcal{M}({\mathbf{p}}_i)$. We assume $l_1(\mathbf{e})>l_2(\mathbf{e})$ and we set $a=\frac{l_1(\mathbf{e})}{l_2(\mathbf{e})}$.

Now, we consider the metric variation law associated with the interpolation operator on metrics.

Logarithmic interpolation is used to interpolate between metrics: $$\mathcal{M}(t)=\exp \{(1-t)\log ({\mathcal{M}}_0)+t\log ({\mathcal{M}}_1)\}$$ as it verifies a maximum principle, i.e. if $\det ({\mathcal{M}}_0)\le \det ({\mathcal{M}}_1)$ then $\det ({\mathcal{M}}_0)\le \det (\mathcal{M}(t))\le \det ({\mathcal{M}}_1)$ for $0\le t\le 1$

When interpolating a discrete metric field at $\mathbf{x}=\sum {\alpha }_i{\mathbf{x}}_i$ is approximated as $$\mathcal{M}(\mathbf{x})=\exp \left(\sum {\alpha }_i\log ({\mathcal{M}}_i)\right)$$ which is consistent with the assumption of geometric progression of sizes.

Why not linear interpolation?

Linear interpolation between ${\mathcal{M}}_0$ and ${\mathcal{M}}_1$ is given by $$\mathcal{M}(t)=(1-t){\mathcal{M}}_0+t{\mathcal{M}}_1$$ and does not respect the maximum principle.

As such, logarithmic interpolation is considered along the edge. It can be shown that it is equivalent to a geometric interpolation of the physical length prescribed by the metric F. Alauzet, leading to the following Riemannian distance computation:

$$l_{\mathcal{M}}(\mathbf{e})=l_1(\mathbf{e})\frac{a-1}{a\text{ln}(a)}$$

Numerical computation of volumes

The evaluation of a tetrahedron volume in a Riemannian metric space consists in computing the volume integral numerically. If a first order approximation in the log-Euclidean framework is considered, we get:

$$|K{|}_{\mathcal{M}}={\int }_K\sqrt{\text{det}\mathcal{M}}d\mathbf{x}\approx \sqrt{\text{det}\left(\text{exp}\left(\frac{1}{4}\sum _{i=1}^4\text{log}({\mathcal{M}}_i)\right)\right)}|K{|}_{{\mathcal{I}}_3}$$

Higher order approximation can be obtained by using Gaussian quadrature. For instance, if one considers a $k$-point Gaussian quadrature with weights $({\omega }_j{)}_{j=\mathrm{1...}k}$ and barycentric coefficients $({\beta }_j^1,{\beta }_j^2,{\beta }_j^3,{\beta }_j^4{)}_{j=\mathrm{1...}k}$, it yields:

$$|K{|}_{\mathcal{M}}\approx |K{|}_{{\mathcal{I}}_3}\sum _{j=1}^k{\omega }_j\sqrt{\text{det}\left(\text{exp}\left(\sum _{i=1}^4{\beta }_j^i\text{log}({\mathcal{M}}_i)\right)\right)}$$

Metric operations

Building a metric that is suitable for remeshing may require combining information from different sources or imposing constraints (e.g. minimum or maximum edge lengths, smooth variation with respect to the input mesh, ... ).

Metric Interpolation

In practice, depending on the CFD solver, the metric is either stored on the vertex or the cell center. The definition of an interpolation procedure on metrics is therefore mandatory to be able to compute the metric at any point of the domain. When more than two discrete metrics needs to be interpolated, it is important to use a consistent interpolation framework. For instance, the interpolation needs to be commutative (i.e. the resulting metric does not depend on the order of the interpolation operations between metrics). To that end, the log-Euclidean framework introduced by Arsigny et al. (2006) is used.

Log-Euclidean framework

The authors first define the notion of metric logarithm and matrix exponential. The metric logarithm is defined on the set of metric tensors. For metric tensor $\mathcal{M}=\mathcal{R}\Lambda {\mathcal{R}}^{\top }$, it is given by:

$$\ln (\mathcal{M}):=\mathcal{R}\ln (\Lambda ){\mathcal{R}}^{\top }$$

where $\ln (\Lambda )=\text{diag}(\ln ({\lambda }_i))$. The matrix exponential is defined on the set of symmetric matrices. For any symmetric matrix $S=\mathcal{Q}\Sigma {\mathcal{Q}}^{\top }$, it is given by:

$$\exp (S):=\mathcal{Q}\exp (\Sigma ){\mathcal{Q}}^{\top }$$

where $\exp (\Sigma )=\text{diag}(\exp ({\xi }_i))$. We can now define the logarithmic addition $\oplus$ and the logarithmic scalar multiplication $\odot$:

$${\mathcal{M}}_1\oplus {\mathcal{M}}_2:=\exp (\ln ({\mathcal{M}}_1)+\ln ({\mathcal{M}}_2))$$

$$\alpha \odot \mathcal{M}:=\exp (\alpha \cdot \ln (\mathcal{M}))={\mathcal{M}}^{\alpha }$$

The logarithmic addition is commutative and coincides with matrix multiplication whenever the two tensors ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ commute in the matrix sense.

Metric interpolation in the log-Euclidean framework

Let $({\mathbf{x}}_i{)}_{i=1,...k}$ be a set of cell centers and $(\mathcal{M}({\mathbf{x}}_i){)}_{i=1,...k}$ their associated metrics. Then, for a point $\mathbf{x}$ of the domain such that:

$$\mathbf{x}=\sum _{i=1}^k{\alpha }_i\cdot {\mathbf{x}}_i\text{ with }\sum _{i=1}^k{\alpha }_i=1$$

the interpolated metric is defined by:

$$\mathcal{M}(\mathbf{x})=\underset{i=1}{\overset{k}{⨁}}{\alpha }_i\odot \mathcal{M}({\mathbf{x}}_i)=\exp \left(\sum _{i=1}^k{\alpha }_i\ln (\mathcal{M}({\mathbf{x}}_i))\right)$$

This interpolation is commutative, but its bottleneck is to perform $k$ diagonalizations and to request the use of the logarithm and the exponential functions which are CPU consuming. However, this procedure is essential to define continuously the metric map on the entire domain. Also, since CODA is a cell centered solver and Tucanos is vertex based, this interpolation definition is crucial to get the metric field at the vertexes of the mesh. Moreover, it has been demonstrated (Arsigny et al., 2006) that this interpolation preserves the maximum principle, i.e., for an edge $\mathbf{ab}$ with endpoint metrics $\mathcal{M}(\mathbf{a})$ and $\mathcal{M}(\mathbf{b})$ such that $\text{det}(\mathcal{M}(\mathbf{a}))<\text{det}(\mathcal{M}(\mathbf{b}))$ then we have:

$$\text{det}(\mathcal{M}(\mathbf{a}))<\text{det}(\mathcal{M}(\mathbf{a}+t\mathbf{ab}))<\text{det}(\mathcal{M}(\mathbf{b})),\forall t\in [0,1]$$

Metric Intersection

In metric-based adaptation, the metric field computed from the error estimation does not take into account geometric or gradation constraints, that would be defined using different metric fields. Metric intersection is used to get a unique metric field. It consists in keeping the most restrictive size constraint in all directions imposed by this set of metrics.

Figure 1: Left, view illustrating the metric intersection procedure with the simultaneous reduction in three dimensions. In red, the resulting metric of the intersection of the blue and green metrics. Right, metric interpolation along a segment where the endpoints metrics are the blue and violet ones. (Loseille 2008)

Formally speaking, let ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ be two metric tensors given at a point. The metric tensor ${\mathcal{M}}_{1\cap 2}$ corresponding to the intersection of ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ is the one prescribing the largest possible size under the constraint that the size in each direction is always smaller than the sizes prescribed by ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$. Let us give a geometric interpretation of this operator. Metric tensors are geometrically represented by an ellipse in 2D and an ellipsoid in 3D. But the intersection between two metrics is not directly the intersection between two ellipsoids as their geometric intersection is not an ellipsoid. Therefore, we seek for the largest ellipsoid representing ${\mathcal{M}}_{1\cap 2}$ included in the geometric intersection of the ellipsoids associated with ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$, cf. the Figure above (left). The ellipsoid (metric) verifying this property is obtained by using the simultaneous reduction of two metrics.

Simultaneous reduction

The simultaneous reduction enables to find a common basis $({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3)$ such that ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ are congruent to a diagonal matrix in this basis, and then to deduce the intersected metric. To do so, the matrix $\mathcal{N}={\mathcal{M}}_1^{-1}{\mathcal{M}}_2$ is introduced. $\mathcal{N}$ is diagonalizable with real-eigenvalues. The normalized eigenvectors of $\mathcal{N}$ denoted by ${\mathbf{e}}_1$, ${\mathbf{e}}_2$ and ${\mathbf{e}}_3$ constitute a common diagonalization basis for ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$. The entries of the diagonal matrices, that are associated with the metrics ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ in this basis, are obtained with the Rayleigh formula:

$${\lambda }_i={\mathbf{e}}_i^{\top }{\mathcal{M}}_1{\mathbf{e}}_i\text{ and }{\mu }_i={\mathbf{e}}_i^{\top }{\mathcal{M}}_2{\mathbf{e}}_i$$

Remark: ${\lambda }_i$ and ${\mu }_i$ are not the eigenvalues of ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$. They are the spectral values associated with basis $({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3)$.

Let $\mathcal{P}=({\mathbf{e}}_1{\mathbf{e}}_2{\mathbf{e}}_3)$ be the matrix with the columns that are the eigenvectors $\{{\mathbf{e}}_i{\}}_{i=\mathrm{1...3}}$ of $\mathcal{N}$. $\mathcal{P}$ is invertible as $({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3)$ is a basis of ${\mathbb{R}}^3$. We have:

$${\mathcal{M}}_1={\mathcal{P}}^{-\top }\left(\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_2& 0\\ 0& 0& {\lambda }_3\end{array}\right){\mathcal{P}}^{-1}\text{ and }{\mathcal{M}}_2={\mathcal{P}}^{-\top }\left(\begin{array}{ccc}{\mu }_1& 0& 0\\ 0& {\mu }_2& 0\\ 0& 0& {\mu }_3\end{array}\right){\mathcal{P}}^{-1}$$

Computing the metric intersection

The resulting intersected metric ${\mathcal{M}}_{1\cap 2}$ is then analytically given by:

$${\mathcal{M}}_{1\cap 2}={\mathcal{M}}_1\cap {\mathcal{M}}_2={\mathcal{P}}^{-\top }\left(\begin{array}{ccc}\max ({\lambda }_1,{\mu }_1)& 0& 0\\ 0& \max ({\lambda }_2,{\mu }_2)& 0\\ 0& 0& \max ({\lambda }_3,{\mu }_3)\end{array}\right){\mathcal{P}}^{-1}$$

The ellipsoid associated with ${\mathcal{M}}_{1\cap 2}$ is the largest ellipsoid included in the geometric intersection region of the ellipsoids associated with ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$.

Numerically, to compute ${\mathcal{M}}_{1\cap 2}$, the real-eigenvalues of $\mathcal{N}$ are first evaluated with a Newton algorithm. Then, the eigenvectors of $\mathcal{N}$, which define $\mathcal{P}$, are computed using the algebra notions of image and kernel spaces.

See proof

Metric gradation

Metric fields may have huge variations or may be quite irregular when evaluated from numerical solutions that present discontinuities or steep gradients. This makes the generation of a unit mesh difficult or impossible, thus leading to poor quality anisotropic meshes. Generating high-quality anisotropic meshes requires to smooth the metric field by bounding its variations in all directions. It also helps with the convergence of the implicit solver. Mesh gradation (Alauzet, 2010) consists in reducing in all directions the size prescribed at any points if the variation of the metric field is larger than a fixed threshold.

Spanning a metric field

Let $\mathbf{p}$ be a point of a domain $\Omega$ supplied with a metric ${\mathcal{M}}_{\mathbf{p}}$ and $\beta$ the specified gradation. Two laws governing the metric growth in the domain are proposed. In the first one, the metric growth is homogeneous in the Euclidean metric field defined by ${\mathcal{M}}_{\mathbf{p}}$. The second metric growth is homogeneous in the physical space, i.e., the classical Euclidean space.

The first law associates for any point $\mathbf{x}$ of the domain a unique scale factor with the metric given by:

$${\eta }^2(\mathbf{px})=(1+l_p(\mathbf{px})\cdot \ln (\beta ){)}^{-2}=(1+\sqrt{{\mathbf{px}}^{\top }{\mathcal{M}}_{\mathbf{p}}\mathbf{px}}\cdot \ln (\beta ){)}^{-2}$$

With this formulation, each pointwise metric ${\mathcal{M}}_{\mathbf{p}}$ spans a global continuous smooth metric field all over domain $\Omega$ parametrized by the given gradation value $\beta$:

$$({\mathcal{M}}_{\mathbf{p}}(\mathbf{x}){)}_{\mathbf{x}\in \Omega }\text{with }{\mathcal{M}}_{\mathbf{p}}(\mathbf{x})={\eta }^2(\mathbf{px}){\mathcal{M}}_{\mathbf{p}}$$

In this case, the resulting metric field grows homogeneously in the Euclidean metric space defined by ${\mathcal{M}}_{\mathbf{p}}$ as the scale factor depends on the length of segment $\mathbf{px}$ with respect to ${\mathcal{M}}_{\mathbf{p}}$. As a result, the shape of the metric is kept unchanged while growing. This law conserves the same anisotropic ratio.

For the second law, we associate independently a growth factor with each eigenvalue of ${\mathcal{M}}_{\mathbf{p}}=\mathcal{R}\Lambda {\mathcal{R}}^{\top }$, with $\Lambda =\text{diag}({\lambda }_i{)}_{i=1,3}$:

$${\eta }_i^2(\mathbf{px})=(1+\sqrt{{\lambda }_i}||\mathbf{px}|{|}_2\cdot \ln (\beta ){)}^{-2}$$

the grown metric at $\mathbf{x}$ is given by:

$${\mathcal{M}}_{\mathbf{p}}(\mathbf{x})=\mathcal{RN}(\mathbf{px})\Lambda {\mathcal{R}}^{\top }\text{where }\mathcal{N}(\mathbf{px})=\left(\begin{array}{ccc}{\eta }_1^2(\mathbf{px})& 0& 0\\ 0& {\eta }_2^2(\mathbf{px})& 0\\ 0& 0& {\eta }_3^2(\mathbf{px})\end{array}\right)$$

The resulting metric field grows homogeneously in the physical space. Indeed, each eigenvalue grows similarly in all directions, as the factor ${\eta }_i$ depends only on the distance (in the physical space) from the original point. Consequently, the shape, i.e., the anisotropic ratio of the metric, is no more preserved as the eigenvalues are growing separately and differently. This law gradually makes the metric more and more isotropic as it gradually propagates in the domain.

Alauzet et al. (2021) suggest to mix these two laws to achieve an efficient metric gradation algorithm. For this new law, a growth factor is associated independently with each eigenvalue of ${\mathcal{M}}_{\mathbf{p}}$:

$${\eta }_i^2(\mathbf{px})=((1+\sqrt{{\lambda }_i}||\mathbf{px}|{|}_2\cdot \ln (\beta ){)}^t(1+l_p(\mathbf{px})\cdot \ln (\beta ){)}^{1-t}{)}^{-2}$$

The authors consider $t=\frac{1}{8}$ within their numerical examples. For the applications in the work, we will use $t=1.$

Metric reduction

The reduced metric at a point $\mathbf{x}$ of the domain $\Omega$ is given by the strongest size constraint imposed by the metric at $\mathbf{x}$ and by the spanned metrics (parametrized by the given size gradation) of all the other points of the domain at $\mathbf{x}$:

$$\widetilde{\mathcal{M}(\mathbf{x})}=\left(\bigcap _{\mathbf{p}\in \Omega }{\mathcal{M}}_{\mathbf{p}}(\mathbf{x})\right)\cap \mathcal{M}(\mathbf{x})$$

Practical implementation are given in Alauzet (2010).

Thresholds on the local sizes

In order to ensure that the metric field will generate meshes with edge sizes between $h_{min}$ and $h_{max}$, the metric $\mathcal{M}=R^{T}diag(s_1,\cdots s_d)R$ with $R=({\mathbf{v}}_1|\cdots |{\mathbf{v}}_d)$ is modified as $$\mathcal{T}(\mathcal{M},h_{min},h_{max})=R^{T}diag(s_1,\cdots s_d)R$$ with ${\tilde{s}}_i=\min \left(\max \left(s_i,\frac{1}{h_{max}^2}\right),\frac{1}{h_{min}^2}\right)$.

Element-implied metric

Given an element $K$, defined by its vertices $({\mathbf{x}}_0,\cdots ,{\mathbf{x}}_d)$, the element-implied metric is the metric for which the element is equilateral

If $J_{K,\Delta }$ is the jacobian of the transformation from the reference equilateral element to the physical element, the element-implied metric is $${\mathcal{M}}_K=(J_{K,\Delta }{J}_{K,\Delta }^T{)}^{-1}$$

In practice $J_{K,\Delta }$ is computed as $J_{K,\Delta }=J_{K,\perp }{J}_{\perp ,\Delta }$ where - $J_{K,\perp }=({\mathbf{x}}_1-{\mathbf{x}}_0|\cdots |{\mathbf{x}}_d-{\mathbf{x}}_0)$ is the jacobian of transformation from the reference orthogonal element to the physical element - $J_{\perp ,\Delta }$ is the jacobian of transformation from the reference equilateral element to the reference orthogonal element (independent of $K$ )

Controlling the step between two metrics

In the remeshing process, it may be interesting to control the step between the element-implied and target metrics. Such control is available within the remesher Tucanos and specified by the parameter $f$.

Given two metric fields ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$, the objective is to find $\mathcal{M}=\mathcal{L}({\mathcal{M}}_1,{\mathcal{M}}_2,f)$ as close as possible from ${\mathcal{M}}_2$ such that, for all edges $e$:

$$\frac{1}{f}\le \sqrt{\frac{e^T\mathcal{M}e}{e^T{\mathcal{M}}_{1}e}}\le f$$

i.e. to have:

$$\frac{1}{f}\le {\lambda }_{\min }({\mathcal{M}}_1^{-\frac{1}{2}}\mathcal{M}{\mathcal{M}}_1^{-\frac{1}{2}})\le {\lambda }_{\max }({\mathcal{M}}_1^{-\frac{1}{2}}\mathcal{M}{\mathcal{M}}_1^{-\frac{1}{2}})\le f$$

Practically, we define as close as possible as minimizing the Froebenius norm $\Vert {\mathcal{M}}_1^{-\frac{1}{2}}(\mathcal{M}-{\mathcal{M}}_2){\mathcal{M}}_1^{-\frac{1}{2}}{\Vert }_F$. The optimal $\mathcal{M}$ is then computed as follows:

• compute $\mathcal{N}:={\mathcal{M}}_1^{-\frac{1}{2}}{\mathcal{M}}_2{\mathcal{M}}_1^{-\frac{1}{2}}$,
• compute the eigenvalue decomposition $QD{Q}^T=\mathcal{N}$, with $D=\text{diag}({\lambda }_i)$
• compute ${\mathcal{N}}^{\ast }:=Q\text{diag}({\hat{\lambda }}_i)Q^T$ where ${\hat{\lambda }}_i:=\min (\max ({\lambda }_i,\frac{1}{f}),f)$,
• compute ${\mathcal{M}}^{\ast }:={\mathcal{M}}_1^{\frac{1}{2}}{\mathcal{N}}^{\ast }{\mathcal{M}}_1^{\frac{1}{2}}$.

Figure 2: Representation of the controlled metric $\mathcal{L}({\mathcal{M}}_1,{\mathcal{M}}_2,f)$ by stepping with stepping parameter $f=2$. The thin blue lines represent $2{\mathcal{M}}_1$ and $\frac{1}{2}{\mathcal{M}}_1$, X. Garnaud

See proof

Metric scaling

• The "ideal" metric $\mathcal{M}$ is modified using a series of constraints
  • Minimum and maximum edge sizes $h_{min}$ and $h_{max}$,
  • Maximum gradation $\beta$,
  • Maximum step $f$ with respect to the element implied metric ${\mathcal{M}}_i$,
  • A minimum curvature radius to edge lenth ratio $r$ (and optionally $y_{target}^{+}$), giving a fixed metric ${\mathcal{M}}_f$
• A scaling $\alpha$ is applied, such that the actual metric used for remeshing is $$\tilde{\mathcal{M}}=\mathcal{L}(\mathcal{T}(\alpha \mathcal{M},h_{min},h_{max})\cap {\mathcal{M}}_f,{\mathcal{M}}_i,f)$$
• $\alpha$ is chosen to have a given complexity $$\mathcal{C}(\tilde{\mathcal{M}})=\mathcal{N}$$
• A simple bisection is used to compute $\alpha$ if all the constraints can be met.

Maths

Notations

• One denotes the Frobenius norm of a matrix $A\in {\mathbb{R}}^{n\times n}$ by $\Vert A{\Vert }_F:=\mathrm{tr}(A^{T}A{)}^{1/2}=({\sum }_{i,j}{A}_{i,j}^2{)}^{1/2}$. This norm is associated with the scalar product $(A,B{)}_F:=\mathrm{tr}(A^{T}B)$.
• One denotes by $\le$ the partial order in the symmetric matrices space ${\mathcal{S}}_n(\mathbb{R})$: $S_1\le S_2$ if $S_2-S_1\in {\mathcal{S}}_n^{+}(\mathbb{R})$, i.e. for all $x\in {\mathbb{R}}^n$, $x^T{S}_{1}x\le x^T{S}_{2}x$.

Projection on a closed convex subset of the symmetric matrices

The following theorem has been proven in Computing a nearest symmetric positive semidefinite matrix, Higham Nicholas J, 1988:

Theorem 1 Let ${\mathcal{S}}_{\delta }:=\{S\in {\mathcal{S}}_n({\mathbb{R}}^n):\delta I_n\le S\}$. Then, ${\mathcal{S}}_{\delta }$ is a closed convex subset of the space ${\mathcal{S}}_n({\mathbb{R}}^n)$ and then there exists a projection application $P_{\delta }:{\mathcal{S}}_n({\mathbb{R}}^n)\to {\mathcal{S}}_{\delta }$ such that $$\forall A\in {\mathcal{S}}_n(\mathbb{R}),\Vert A-P_{\delta }(A){\Vert }_F=\inf \{\Vert A-S{\Vert }_F:S\in {\mathcal{S}}_{\delta }\}.>$$ Let $A\in {\mathcal{S}}_n(\mathbb{R})$ and $A=Q^{T}DQ$ its diagonalisation in an orthonormal basis: $Q\in {\mathcal{O}}_n(\mathbb{R})$, $D=\mathrm{diag}({\lambda }_i)$. Then, $$P_{\delta }(A)=Q\hat{D}{Q}^T,\text{where}\hat{D}:=\mathrm{diag}(\max ({\lambda }_i,\delta )).$$

Proof

Let $A\in {\mathcal{S}}_n(\mathbb{R})$ and its diagonalisation in a orthonomal basis $P^{T}AP=D=\mathrm{diag}({\lambda }_i)$. As $\Vert .{\Vert }_F$ is invariant by linear isometry ($\Vert PAQ{\Vert }_F^2=\mathrm{tr}((PAQ{)}^{T}PAQ)=\mathrm{tr}(Q^T{A}^T{P}^{T}PAQ)=\mathrm{tr}(Q{Q}^T{A}^{T}A)=\Vert A{\Vert }_F^2$, one gets $\Vert A-S{\Vert }_F^2=\Vert D-Q^{T}SQ{\Vert }_F^2$ and then the change of variable $S\in {\mathcal{S}}_n(\mathbb{R})\mapsto X=Q^{T}SQ\in {\mathcal{S}}_n(\mathbb{R})$ together with $\delta I\le S\iff \delta I\le Q^{T}SQ$, leads to the following problem $$\inf \{\Vert D-X{\Vert }_F^2:X\in {\mathcal{S}}_n(\mathbb{R}),\delta I\le X\}.$$ This problem can be computed explicitly. Since $X_{ii}=e_i^{T}X{e}_i\ge \delta e_i^T{e}_i=\delta$, ($e_i=({\delta }_{ij}{)}_j$ the canonical basis vectors) $$\Vert D-X{\Vert }_F^2=\sum _{i=1}^n({\lambda }_i-X_{ii}{)}^2+2\sum _{i<j}{X}_{ij}^2\ge \sum _{i=1}^n({\lambda }_i-X_{ii}{)}^2\ge \sum _{i:{\lambda }_i<\delta }({\lambda }_i-\delta {)}^2.$$ This lower bound is reached for the following diagonal matrix $X:=\hat{D}:=\mathrm{diag}({\hat{\lambda }}_i)$, where ${\hat{\lambda }}_i:=\max ({\lambda }_i,\delta )$, and moreover, $\hat{D}\ge \delta I$ and is in ${\mathcal{S}}_n(\mathbb{R})$, which ends the proof.

Generalization

One can generalize the Theorem above as follow:

Theorem 2 Given two symmetric real numbers $0<\alpha <\beta$, let ${\mathcal{S}}_{\alpha ,\beta }:=\{S\in {\mathcal{S}}_n(\mathbb{R}):\alpha I_n\le S\le \beta I_n\}$. ${\mathcal{S}}_{\alpha ,\beta }$ is a compact convex subset of ${\mathcal{S}}_n(\mathbb{R})$ and then let us denote $P_{\alpha ,\beta }:{\mathcal{S}}_n(\mathbb{R})\to {\mathcal{S}}_{\alpha ,\beta }$ the associated projection application: $$\forall A\in {\mathcal{S}}_n(\mathbb{R}),\Vert A-P_{\alpha ,\beta }(A){\Vert }_F=\inf \{\Vert A-S{\Vert }_F:S\in {\mathcal{S}}_{\alpha ,\beta }\}.>$$ Given $A\in {\mathcal{S}}_n(\mathbb{R})$ and $A=QD{Q}^T$ its diagonalisation in an orthonormal basis: $Q\in {\mathcal{O}}_n(\mathbb{R})$, $D=\mathrm{diag}({\lambda }_i)$. Then, $$P_{\alpha ,\beta }(A)=Q\hat{D}{Q}^T,\text{where}\hat{D}:=\mathrm{diag}({\hat{\lambda }}_i),{\hat{\lambda }}_i:=\min \left(\max \right({\lambda }_i,\alpha ),\beta ).$$

Proof

The proof follows the one above very closely. The lower bound is got using the inequality: $\alpha e_i^T{I}_n{e}_i\le e_i^{T}X{e}_i\le \beta e_i^T{I}_n{e}_i\Rightarrow \alpha \le X_{ii}\le \beta$: $$\Vert D-X{\Vert }_F^2\ge \sum _{i=1}^n({\lambda }_i-X_{ii}{)}^2+2\sum _{i<j}{X}_{ij}^2$$ $$\ge \sum _{i=1}^n({\lambda }_i-X_{ii}{)}^2$$ $$\ge \sum _{i:{\lambda }_i<\alpha }({\lambda }_i-\alpha {)}^2+\sum _{i:{\lambda }_i>\beta }({\lambda }_i-\beta {)}^2,$$ which is reached for $\hat{D}:=\mathrm{diag}({\hat{\lambda }}_i)$ where ${\hat{\lambda }}_i=\alpha$ if ${\lambda }_i<\alpha$, ${\hat{\lambda }}_i=\beta$ if ${\lambda }_i>\beta$ and otherwise ${\hat{\lambda }}_i={\lambda }_i$. Again, the matrix $\hat{D}$ is symmetric and indeed in ${\mathcal{S}}_{\alpha ,\beta }$ (because since $\hat{D}$ is diagonal, $\alpha \le {\hat{\lambda }}_i\le \beta \iff \alpha I_n\le \hat{D}\le \beta I_n$) and thus is the unique solution of the problem (existence and uniqueness being given by the projection on closed convex sets in a Hilbert space).

Metric optimization problem

Given two metrics ${\mathcal{M}}_1,{\mathcal{M}}_2\in {\mathcal{S}}_n^{++}(\mathbb{R})$. The deformation ratio between a metric $\mathcal{M}\in S_n^{++}(\mathbb{R})$ and ${\mathcal{M}}_0$ is defined as a function of x∈Rn</span>{0} r(M,M0​;x):=(xTM0​xxTMx​)1/2.

The problem.

Given two metrics ${\mathcal{M}}_1,{\mathcal{M}}_2\in {\mathcal{S}}_n^{++}(\mathbb{R})$, and two real numbers $0<\alpha <\beta$, let ${\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)$ be the set of metrics which are not to stretched with regard to ${\mathcal{M}}_0$: Rα,β​(M0​):={M∈Sn​(R):∀x∈Rn</span>{0},α≤r(M,M0​;x)2≤β}. Again, this set is compact and convex. Given a matrix norm $$|||.|||$$, we want to find $${\mathcal{M}}^{\ast }\in {\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)$$ which minimize the distance to $${\mathcal{M}}_1$$, i.e. find $${\mathcal{M}}^{\ast }\in {\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)$$ such that: $$|||{\mathcal{M}}_1-{\mathcal{M}}^{\ast }|||=\inf \{|||{\mathcal{M}}_1-\mathcal{M}|||:\mathcal{M}\in {\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)\}$$

Lemma : Let $\mathcal{M},{\mathcal{M}}_0\in {\mathcal{S}}_n^{++}(\mathbb{R})$ two metrics > and $\alpha ,\beta \in \mathbb{R}$. The function x∈Rn</span>{0}↦r(M,M0​;x)2 is bounded and reaches its bounds which are the minimum and maximum of the eigenvalues of M0>​−1/2MM0−1/2​. Thus, >∀x∈Rn</span>{0},>α≤r(M,M0​;x)≤β> ⇔>αIn​≤M0−1/2​MM0−1/2​≤β>In​> ⇔>αM0​≤M≤βM0​>

Thanks to this lemma, the admissible set ${\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)$ can be re-written as $${\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)=\{\mathcal{M}\in {\mathcal{S}}_n(\mathbb{R}):\alpha {\mathcal{M}}_0\le \mathcal{M}\le \beta {\mathcal{M}}_0\}$$ $$=\{\mathcal{M}\in {\mathcal{S}}_n(\mathbb{R}):\alpha I_n\le {\mathcal{M}}_0^{-1/2}\mathcal{M}{\mathcal{M}}_0^{-1/2}\le \beta I_n\}.$$