tucanos and pytucanos

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

Pytucanos is the python wrapper and it can be found here.

This document briefly presents the theory behind the remesher, as some aspects are new to other remesher.

Mesh

Mesh adaptation applies only to simplex meshes, i.e. triangles in 2D and tetrahedra 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
    • belong to only one element. In this case, they must be oriented outwards
    • belong to two elements with different tags
    • (belong to more than two elements)
  • untagged

and 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

Metric

A metric field is used in order to prescribe the element sizes and orientations to be used during remeshing. The continuous mesh framework establishes a duality between this (continuous) metric field and the resulting discrete mesh : 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 & unit mesh

A Riemanian metric is a symmetric positive definite (SPD) tensor field $\mathcal{M}(\mathbf{x})$, defined everywhere in the computational domain, such that the length of an edge $\mathbf{y}-\mathbf{x}$ is given by $${\ell }_{\mathcal{M}}(\mathbf{y}-\mathbf{x})=\Vert \mathbf{y}-\mathbf{x}{\Vert }_{\mathcal{M}}={\int }_0^1\sqrt{(\mathbf{y}-\mathbf{x}{)}^T\mathcal{M}((1-t)\mathbf{x}+t\mathbf{y})(\mathbf{y}-\mathbf{x})}dt$$

This field can then be used to define the size and orientations of mesh elements through the generation of a (quasi) unit mesh, i.e. a mesh for which all edges $\mathbf{e}$ have a unit length: $${\ell }_{\mathcal{M}}(\mathbf{e})\approx 1$$

Local sizes and directions

The metric is a SPD tensor: it can be characterized locally by its eigenvalues $(s_1,\cdots ,s_d)$ and an orthonormal basis of eigenvectors $({\mathbf{v}}_1,\cdots ,{\mathbf{v}}_d)$. The edges $${\mathbf{e}}_i=h_i{\mathbf{v}}_i$$ with a local size $h_i=1/\sqrt{s_i}$ are unit edges in metric space (if $\mathcal{M}$ is constant over ${\mathbf{e}}_i$ ).

In order to generate an isotropic mesh with local size $h(x)$, the metric should be $\mathcal{M}(x)=diag(1/h^2(x),\cdots ,1/h^2(x))$

Metric complexity

The physical volume of and equilateral element with unit edges in the (uniform) metric space (i.e. ideal elements) in $d$ dimensions is $$\frac{v_{\Delta }^{(d)}}{\sqrt{\det (\mathcal{M})}}\equiv v_{\Delta }^{(d)}\mathcal{V}(\mathcal{M})$$ where $\mathcal{V}(\mathcal{M})$ is the volume associated with metric $\mathcal{M}$ and $v_{\Delta }^{(d)}=\frac{\sqrt{d+1}}{d!\sqrt{2^d}}$.

The number of ideal elements per unit physical volume is thefore $1/\left(v_{\Delta }^{(d)}\mathcal{V}(\mathcal{M})\right)$ so the ideal elements to fill a domain $\Omega$ is $$\frac{1}{v_{\Delta }^{(d)}}{\int }_{\Omega }\sqrt{\det (\mathcal{M})}dx$$

The complexity is $$\mathcal{C}(\mathcal{M})={\int }_{\Omega }\sqrt{\det (\mathcal{M})}dx$$

Discrete metric field

In the context of remeshing, only a discrete metric field is known. A metric ${\mathcal{M}}_i$ is defined at each vertex ${\mathbf{x}}_i$ of the input mesh.

Metric interpolation

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:

Edge length in metric space

In order to compute the length of an edge ${\mathbf{e}}_{i,j}={\mathbf{x}}_j-{\mathbf{x}}_i$, 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: $$l_{\mathcal{M}}(e_{i,j})=l_i\frac{a-1}{a\ln (a)}$$ with $l_k=\sqrt{e_{i,j}^T{\mathcal{M}}_k{e}_{i,j}}$ and $a=l_j/l_i$.

Element quality in metric space

The quality of an element in metric space is given by $$q(K)={\beta }_d\frac{|K{|}_{\mathcal{M}}^{2/d}}{{\sum }_e||e|{|}_{\mathcal{M}}^2}$$ where the normalization factor ${\beta }_d$ is a such that $q=1$ of equilateral elements:

• ${\beta }_2=1/(6\sqrt{2})$
• ${\beta }_3=\sqrt{3}/4$

The same metric is used to compute all the edge lengths: among the metrics defined at all the vertices of element $K$, the one with minimum volume is used in the formula above

Invalid elements have a quality $q<0$

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 intersection

Given two metrics ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$, the intersection ${\mathcal{M}}_{1\cap 2}$ corresponds to a metric that imposes the largest sizes that are smaller than those imposed by ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$: $${\mathcal{M}}_{1\cap 2}=\arg \min \{\det (\mathcal{M})|\mathcal{M}\in {\mathcal{S}}^{+}s.t.\mathcal{M}\ge {\mathcal{M}}_1,\mathcal{M}\ge {\mathcal{M}}_2\}$$

It is computed as follows: let $\mathcal{P}=(e_0|...|e_d)$ be the generalized eigenvectors of $({\mathcal{M}}_0,{\mathcal{M}}_1)$: $${\mathcal{M}}_0\mathcal{P}=\Lambda {\mathcal{M}}_1\mathcal{P}$$ then $${\mathcal{M}}_i={\mathcal{P}}^{-1,T}{\Lambda }^{(i)}{\mathcal{P}}^{(-1)}$$ with ${\Lambda }_{jk}^{(i)}=e_j^T{\mathcal{M}}_i{e}_k{\delta }_{jk}$. The intersection is then $${\mathcal{M}}_0\cap {\mathcal{M}}_1={\mathcal{P}}^{-1,T}{\Lambda }^{(i,j)}{\mathcal{P}}^{-1}$$ with ${\Lambda }_{jk}^{(i,j)}=\max ({\Lambda }_{jk}^{(i)},{\Lambda }_{jk}^{(j)})$.

See proof

Metric gradation

It is sometimes desirable to control how fast edge lengths grow from from element to the next. This idea can be translated into constraints on the variation of the metric field on the edges of the mesh (see Size gradation control of anisotropic meshes, F. Alauzet, 2010 pdf) summarized as follows:

• The gradation on an edge ${\mathbf{e}}_{i,j}={\mathbf{x}}_j-{\mathbf{x}}_i$ is (with $a=l_{{\mathcal{M}}_i}(e_{i,j})/l_{{\mathcal{M}}_j}(e_{i,j})$ ) $$\max {\left(a,\frac{1}{a}\right)}^{\frac{1}{l_{\mathcal{M}}(e_{i,j})}}$$
• Given a metric $\mathcal{M}(\mathbf{x})=R^{T}diag(s_1^2,\cdots s_d^2)R$, is is possible to span a field at any location $\mathbf{y}$ as $${\mathcal{M}}_s(\mathbf{x},\mathbf{y})=R^{T}diag({\eta }_1^2{s}_1^2,\cdots {\eta }_d^2{s}_d^2)R$$
  with ${\eta }_i=1+s_i\Vert \mathbf{y}-\mathbf{x}{\Vert }_2\log (\beta )$ so that the gradation along $\mathbf{y}-\mathbf{x}$ is $\beta$
• In practice a maximum gradation is enforces along edge ${\mathbf{e}}_{i,j}={\mathbf{x}}_j-{\mathbf{x}}_i$, by modifying ${\mathcal{M}}_k$ to $$\widetilde{{\mathcal{M}}_j}={\mathcal{M}}_j\cap {\mathcal{M}}_s({\mathbf{x}}_i,{\mathbf{x}}_j)$$
  and $$\widetilde{{\mathcal{M}}_i}={\mathcal{M}}_i\cap {\mathcal{M}}_s({\mathbf{x}}_j,{\mathbf{x}}_i)$$
• Achieving a maximum gradation over all the edges in the mesh would require ${\#}_{edges}^2$ operations; in practices, the above operations is only applied on all the edges of the mesh a small number of times.

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$ )

Controling the step between two metrics

In the remeshing process, it may be interesting to limit the step between the element-implied and target metrics

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$ $$1/f\le \sqrt{\frac{e^T\mathcal{M}e}{e^T{\mathcal{M}}_{1}e}}\le f$$ i.e. to have $$1/f\le {\lambda }_{min}({\mathcal{M}}_1^{-1/2}\mathcal{M}{\mathcal{M}}_1^{-1/2})\le {\lambda }_{max}({\mathcal{M}}_1^{-1/2}\mathcal{M}{\mathcal{M}}_1^{-1/2})\le f$$

• If as close as possible is defined as minimizing the Froebenius norm $\Vert {\mathcal{M}}_0^{-1/2}(\mathcal{M}-{\mathcal{M}}_2){\mathcal{M}}_0^{-1/2}{\Vert }_F$, the optimal $\mathcal{M}$ is computed as follows
  • compute ${\mathcal{N}}_1:={\mathcal{M}}_0^{-1/2}{\mathcal{M}}_1{\mathcal{M}}_0^{-1/2}$,
  • compute the eigenvalue decomposition $QD{Q}^T={\mathcal{N}}_1$,
  • compute ${\mathcal{N}}^{\ast }:=Q\mathrm{diag}({\hat{\lambda }}_i)Q^T$ where ${\hat{\lambda }}_i:=\min (\max ({\lambda }_i,1/f),f)$,
  • compute ${\mathcal{M}}^{\ast }:={\mathcal{M}}_0^{1/2}{\mathcal{N}}^{\ast }{\mathcal{M}}_0^{1/2}$.

(the thin blue lines represent $2{\mathcal{M}}_1$ and $1/2{\mathcal{M}}_1$)

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 beeing 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-writen 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\}.$$

In the ${\mathcal{M}}_0$-twisted Frobenius norm

First, one studies the case of choosing $$|||\mathcal{M}|||:=\Vert {\mathcal{M}}_1-\mathcal{M}{\Vert }_{{\mathcal{M}}_0,F}:=\Vert {\mathcal{M}}_0^{-1/2}\mathcal{M}{\mathcal{M}}_0^{-1/2}{\Vert }_F.$$ The above problem reads as follows. Find ${\mathcal{M}}^{\ast }\in {\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)$ such that: $$|||{\mathcal{M}}_1-{\mathcal{M}}^{\ast }|||=\inf \{||{\mathcal{M}}_0^{-1/2}({\mathcal{M}}_1-\mathcal{M}){\mathcal{M}}_0^{-1/2}|{|}_F:\mathcal{M}\in {\mathcal{S}}_n(\mathbb{R}),{\mathcal{M}}_0^{-1/2}\mathcal{M}{\mathcal{M}}_0^{-1/2}\in {\mathcal{S}}_{\alpha ,\beta }\}.$$ It is then clear that the change of variable $\mathcal{M}\in {\mathcal{S}}_n(\mathbb{R})\to \mathcal{N}:={\mathcal{M}}_0^{-1/2}\mathcal{M}{\mathcal{M}}_0^{-1/2}\in {\mathcal{S}}_n(\mathbb{R})$ leads to a problem of the form of Theorem 2. Denoting ${\mathcal{N}}_1:={\mathcal{M}}_0^{-1/2}{\mathcal{M}}_1{\mathcal{M}}_0^{-1/2}$ the problem reads find ${\mathcal{M}}^{\ast }={\mathcal{M}}_0^{1/2}{\mathcal{N}}^{\ast }{\mathcal{M}}_0^{1/2}\in {\mathcal{R}}_{\alpha ,\beta }({\mathcal{M}}_0)$ such that $$\Vert {\mathcal{N}}_1-{\mathcal{N}}^{\ast }{\Vert }_F=\inf \{\Vert {\mathcal{N}}_1-\mathcal{N}{\Vert }_F:\mathcal{N}\in {\mathcal{S}}_{\alpha ,\beta }\}.$$ Then, from Theorem 2, the solution is $${\mathcal{M}}^{\ast }={\mathcal{M}}_0^{1/2}{P}_{\alpha ,\beta }({\mathcal{N}}_1){\mathcal{M}}_0^{1/2}$$ and can be computed as follows:

• Compute ${\mathcal{N}}_1:={\mathcal{M}}_0^{-1/2}{\mathcal{M}}_1{\mathcal{M}}_0^{-1/2}$,
• Find $Q\in {\mathcal{O}}_n(\mathbb{R})\text{ and }D=\mathrm{diag}({\lambda }_i)$ such that $QD{Q}^T={\mathcal{N}}_1$.
• Compute ${\mathcal{N}}^{\ast }:=P_{\alpha ,\beta }({\mathcal{N}}_1)=Q\mathrm{diag}({\hat{\lambda }}_i)Q^T,\text{where }{\hat{\lambda }}_i:=\min \left(\max \right({\lambda }_i,\alpha ),\beta )$.
• Compute ${\mathcal{M}}^{\ast }:={\mathcal{M}}_0^{1/2}{\mathcal{N}}^{\ast }{\mathcal{M}}_0^{1/2}$

Remark Maybe a more suitable formulation would minimizing the Frobenius norm $|||{\mathcal{M}}_1-\mathcal{M}|||:=\Vert {\mathcal{M}}_1-\mathcal{M}{\Vert }_F$. This leads to a different problem. For instance, if $${\mathcal{M}}_0=\left(\begin{array}{cc}1& 0\\ 0& 1/\epsilon \end{array}\right)$$ where $0<\epsilon \ll 1$. One gets, denoting the error $\mathcal{E}:=(e_{ij}):={\mathcal{M}}_1-{\mathcal{M}}^{\ast }$ $$\Vert \mathcal{E}{\Vert }_{{\mathcal{M}}_0,F}^2=e_{11}^2+2\epsilon e_{12}^2+{\epsilon }^2{e}_{22}^2,$$ instead of $$\Vert \mathcal{E}{\Vert }_F^2=e_{11}^2+2e_{12}^2+e_{22}^2.$$ The two distance measures give very different results... does it matter in practice?

Intersection of metrics

Problem.

Given two metrics ${\mathcal{M}}_1,{\mathcal{M}}_2\in {\mathcal{S}}_n^{++}(\mathbb{R})$, we want to find a metric ${\mathcal{M}}_{1\cap 2}\in {\mathcal{S}}_n^{++}(\mathbb{R})$ which respects the most restrictive length of ${\mathcal{M}}_1$ and ${\mathcal{M}}_1$. Moreover, we want that this metric beeing maximal. More formally, the lenght restriction can be expressed by ∀x∈Rn</span>{0},xTM1∩2​x≤xTM1​xandxTM1∩2​x≤xTM2​x, indeed, $$\sqrt{x^T\mathcal{M}x}$$ is the length of $$x$$ associated with a metric $$\mathcal{M}$$. The local volume associated with a metric $$\mathcal{M}$$ is its determinant: $$\det \mathcal{M}$$, so that the optimization problem reads $$\max \det {\mathcal{M}}_{1\cap 2}$$ such that $${\mathcal{M}}_{1\cap 2}\in {\mathcal{S}}_n(\mathbb{R}),$$ $$0\le {\mathcal{M}}_{1\cap 2}\le {\mathcal{M}}_1,$$ $$0\le {\mathcal{M}}_{1\cap 2}\le {\mathcal{M}}_2.$$

Solution of the problem.

First, one multiplies by ${\mathcal{M}}_1^{-1/2}\in {\mathcal{S}}_n^{++}(\mathbb{R})$ on each side of the constraints, which leads to $$0\le {\mathcal{M}}_1^{-1/2}{\mathcal{M}}_{1\cap 2}{\mathcal{M}}_1^{-1/2}\le I_n$$ $$0\le {\mathcal{M}}_1^{-1/2}{\mathcal{M}}_{1\cap 2}{\mathcal{M}}_1^{-1/2}\le {\mathcal{M}}_1^{-1/2}{\mathcal{M}}_2{\mathcal{M}}_1^{-1/2}.$$ Then, denoting ${\overline{\mathcal{M}}}_2:={\mathcal{M}}_1^{-1/2}{\mathcal{M}}_2{\mathcal{M}}_1^{-1/2}$ and applying the change of variable $\mathcal{M}\mapsto \mathcal{N}:={\mathcal{M}}_1^{-1/2}\mathcal{M}{\mathcal{M}}_1^{-1/2}$, one gets the simplified problem: $$\max \det \left({\mathcal{M}}_1^{1/2}\mathcal{N}{\mathcal{M}}_1^{1/2}\right)$$ such that $$\mathcal{N}\in {\mathcal{S}}_n(\mathbb{R}),$$ $$0\le \mathcal{N}\le I_n,$$ $$0\le \mathcal{N}\le {\overline{\mathcal{M}}}_2.$$

This problem can be even more simplified using $\det \left({\mathcal{M}}_1^{1/2}\mathcal{N}{\mathcal{M}}_1^{1/2}\right)=\det {\mathcal{M}}_1\det \mathcal{N}$ (does not change the $\max$ objective as $\det {\mathcal{M}}_1>0$) and by diagonalization of ${\overline{\mathcal{M}}}_2$ in an orthonomal basis: $${\overline{\mathcal{M}}}_2=P^{T}DP,\text{with }{P}^{T}P=I_n.$$ Indeed, $0\le \mathcal{N}\le P^{T}DP\Rightarrow 0\le P\mathcal{N}{P}^T\le D$ and $0\le \mathcal{N}\le I_n\Rightarrow 0\le P\mathcal{N}{P}^T\le I_n$, and let us change again the variable $\mathcal{X}:=P\mathcal{N}{P}^T$, with $\det \mathcal{N}=\det \mathcal{X}$, the problem is finally: $$\max \det \mathcal{X}$$ such that $$\mathcal{X}\in {\mathcal{S}}_n(\mathbb{R}),$$ $$0\le \mathcal{X}\le I_n,$$ $$0\le \mathcal{X}\le D.$$

This problem can now to be solved explicitly by finding an upper bound. Let $\mathcal{Y}\in {\mathcal{S}}_n(\mathbb{R})$ such that $\mathcal{X}={\mathcal{Y}}^2$, then Hadamard inequality provides $\det \mathcal{X}=(\det \mathcal{Y}{)}^2\le {\prod }_j|{\mathcal{Y}}_{\bullet i}{|}_2^2$. Moreover, ${\mathcal{X}}_{jj}={\sum }_k{\mathcal{Y}}_{jk}{\mathcal{Y}}_{kj}=|{\mathcal{Y}}_{\bullet j}{|}_2^2$, from which $\det \mathcal{X}\le {\prod }_i{\mathcal{X}}_{ii}$ and then, using the two constraints: $0\le \mathcal{X}\le I_n\text{ and }0\le \mathcal{X}\le D\Rightarrow 0\le {\mathcal{X}}_{ii}\le \min (D_{ii},1)$, one gets the upper bound $$\det \mathcal{X}\le \prod _{i=1}^n{\mathcal{X}}_{ii}\le \prod _{i=1}^n\min (D_{ii},1),$$ reached for the diagonal matrix $D^{\ast }:=\mathrm{diag}(\min (D_{ii},1))$. Moreover, this matrix is symmetric and verifies $0\le D^{\ast }\le I_n$ and $0\le D^{\ast }\le D$.

Finally, applying the two changes of variables, the unique solution of the problem is $${\mathcal{M}}_{1\cap 2}={\mathcal{M}}_1^{1/2}{P}^T{D}^{\ast }P{\mathcal{M}}_1^{1/2}$$

Hadamard inequality

Given a matrix $M=[v_1|\dots |v_n]\in {\mathbb{C}}^{n\times n}$, one has $$\left|\det M\right|\le \prod _{i=1}^n|v_i{|}_2,$$ with equality iff the $(v_1,\dots ,v_n)$ are orthogonal.

The result if clear if $M$ is singular, if not let $\hat{M}:=[{\hat{v}}_1|\dots |{\hat{v}}_n]$, ${\hat{v}}_i=v_i/|v_i{|}_2$ so that one has to show that $|\det \hat{M}|\le 1$. ${\hat{M}}^{\ast }\hat{M}\in {\mathcal{H}}_n(\mathbb{C})$ is hermitian hence diagonalizable in an orthonomal basis, denoting $({\lambda }_1,\dots ,{\lambda }_n)$ its eigenvalues, one has $$|\det \hat{M}{|}^2=\det {\hat{M}}^{\ast }\hat{M}=\prod _{i=1}^n{\lambda }_i$$ $$\le {\left(\frac{1}{n}\sum _{i=1}^n{\lambda }_i\right)}^n={\left(\frac{1}{n}\mathrm{tr}\left({\hat{M}}^{\ast }M\right)\right)}^n={\left(\frac{1}{n}\sum _{i=1}^n\sum _{k=1}^n\overline{{\hat{M}}_{ki}}{\hat{M}}_{ki}\right)}^n$$ $$={\left(\frac{1}{n}\sum _{i=1}^n\sum _{k=1}^n|{\hat{M}}_{ki}{|}^2\right)}^n={\left(\frac{1}{n}\sum _{i=1}^n|{\hat{v}}_i{|}_2^2\right)}^n$$

Curvature metric

A metric field can be defined on the boundaries of the computational domain to ensure an accurate representation of the curvature

• the curvature tensor is computed on an accurate geometrical model, giving principal directions ${\mathbf{t}}_0$ and ${\mathbf{t}}_1$ and the associated curvature radii $r_0$ and $r_1$.
• along direction ${\mathbf{t}}_{\mathbf{i}}$, a size $h_i=r_i/c$ is imposed (with a given value of $c$, typically $c\sim 5$)
• along the normal direction $\mathbf{n}={\mathbf{t}}_0\times {\mathbf{t}}_1$ the edge length $h_n$ can be chosen arbitrarily. It the representation of the curvature is the only objective, one can chose $h_n=\min (h_0,h_1)$.

This metric field, defined on the boundary of the computational domain, can then be extended into the domain using a fixed gradation

boundary-normal size for turbulent boundary layers

For the simulation of turbulent boundary layers, it may be interesting to ensure that the mesh will be fine enough to capture the linear (or log if a wall model is available) layer. In this case, the boundary normal size $h_n$ can be chosen to have a target value of $y^{+}$.

In order to estimate the wall shear stress, required to estimate the target edge length in the normal direction, a wall model should be used (as the initial mesh is in practive much too coarse to resolve the turbulent boundary layer)

Remesher

Input metric	Adapted mesh

The idea is to perform a series of local modification on an input mesh to have a

• edge lengths that are of approximately unit length in the metric space
• element qualities as high as possible

In order to ensure the validity of the final mesh, validity is enforced at each step.

Cavity operators

Remesher implementation: the cavity operator

Given a mesh entity $e$ (a vertex or an edge in the present case), the cavity $C(e)$ is the set of mesh elements that contain the entity $e$.

A cavity $C(e)$ can be filled from a vertex $v$ as the set of elements created from $\mathbf{v}$ and the faces of the boundary of the cavity $\partial C(e)$ (oriented outwards). $$\mathcal{F}(\mathbf{v},C(e))=\{K=(\mathbf{v},{\mathbf{g}}_1,\cdots ,{\mathbf{g}}_d)|g=({\mathbf{g}}_1,\cdots ,{\mathbf{g}}_d)\in \partial C(e),\mathbf{v}\notin g\}$$

Edge collapse using cavity operators

The collapse of edge $e_{i,j}=({\mathbf{x}}_i,{\mathbf{x}}_j)$ can be performed as follow

• check the tags
  • if $dim(e_{i,j})>dim({\mathbf{x}}_i)$ and $dim(e_{i,j})>dim({\mathbf{x}}_j)$, do not collapse
  • if $dim(e_{i,j})>dim({\mathbf{x}}_i)$ and $dim(e_{i,j})=dim({\mathbf{x}}_j)$, swap indices $i$ and $j$
• build the vertex cavity $C({\mathbf{x}}_i)$
• fill the vertex cavity from ${\mathbf{x}}_j$: $F=\mathcal{F}({\mathbf{x}}_j,C({\mathbf{x}}_i))$
• if $F$ passes all the criteria for swaps
  • remove all the elements from $C({\mathbf{x}}_i)$ from the mesh
  • insert all the elements from $F$ to the mesh

What to do with the metric ${\mathcal{M}}_i$?

Edge split using cavity operators

The collapse of edge $e_{i,j}$ can be performed as follow

• build the edge cavity $C(e_{i,j})$
• create the midpoint $\mathbf{x}=({\mathbf{x}}_i+{\mathbf{x}}_j)/2$, project it onto the geometry if needed, and compute the metric $\mathcal{M}=\exp ((\log ({\mathcal{M}}_i)+\log ({\mathcal{M}}_j))/2)$
• fill the vertex cavity from $\mathbf{x}$: $F=\mathcal{F}(\mathbf{x},C({\mathbf{e}}_{i,j}))$
• if $F$ passes all the criteria for swaps
  • remove all the elements from $C({\mathbf{e}}_{i,j})$ from the mesh
  • insert verted $\mathbf{x}$ (with metric $\mathcal{M}$) to the mesh
  • insert all the elements from $F$ to the mesh

Edge swap using cavity operators

The swap of edge $e_{i,j}$ can be performed as follow

• build the edge cavity $C(e_{i,j})$
• compute $q_{min}={\min }_{K\in C(e_{i,j})}q(K)$
• if $q_{min}$ is too low, loop over all the vertices $\mathbf{x}$ on the boundary of C (that are neither ${\mathbf{x}}_i$ nor ${\mathbf{x}}_j$)
  • fill the vertex cavity from $\mathbf{x}$: $F=\mathcal{F}(\mathbf{x},C({\mathbf{e}}_{i,j}))$
  • if $F$ passes all the criteria, it is a valid candidate
• if valid candidates are available, let $F_0$ be the best one (the one with for which the minimum element quality is maximum)
  • remove all the elements from $C({\mathbf{e}}_{i,j})$ from the mesh
  • insert all the elements from $F_0$ to the mesh

Smoothing using cavity operators

The smoothnig of vertex ${\mathbf{x}}_i$ can be performed as follow

• build the vertex cavity $C({\mathbf{x}}_i)$
• compute $q_{min}={\min }_{K\in C({\mathbf{x}}_i)}q(K)$
• compute the subset of the vertices in $C({\mathbf{x}}_i)$ that are tagged in the same entity as ${\mathbf{x}}_i$ or one of its children, and average the positions of these vertices to obtain the new position $\mathbf{x}$
• if ${\mathbf{x}}_i$ is on a boundary, project $\mathbf{x}$ on the geometry
• fill the vertex cavity from $\mathbf{x}$: $F=\mathcal{F}(\mathbf{x},C({\mathbf{x}}_i))$
• if ${\min }_{K\in F}q(K)<q_{min}$:
  • locate the element of $C({\mathbf{x}}_i)$ that contains $\mathbf{x}$ (or is the closest), and interpolate the metric in $\mathbf{x}$ using the barycentric coordinates
  • update the location and metric of vertex $i$

Remeshing quality

Remeshing loop

Benchmarks

Isotropic remeshing in a 2D square domain

Anisotropic remeshing in a 2D square domain

Isotropic remeshing in a 3D cubic domain

Anisotropic remeshing in a 3D cubic domain

polar-1 benchmark from the UGAWG