tucanos and pytucanos

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

The remeshing algorithm is a sequence of steps applied to the edges or vertices of the mesh

• collapse
• split
• (edge) swap
• smooth

Collapse

The collapse operation aims at removing "small" edges. It is applied to edges whose length (in metric space) is smaller than $l_0<1/\sqrt{2}$.

It can however introduce

• "long" edges
• element of low quality (including invalid elements)
• a poor representation of the geometry that will be difficult to recovered later

These 3 criteria have to be assessed to determine if a collapse operation is accepted or not

Split

The collapse operations aims at splitting "long" edges. It is applied to edges whose length (in metric space) is larger than $l_0>\sqrt{2}$.

It can however introduce

• "short" edges
• element of low quality (including invalid elements)

These 2 criteria have to be assessed to determine if a split operation is accepted or not

When introducing new vertices on boundaries, a projection step is required to ensure the consistency with the CAD model.

Edge swap

The swap operation aims at improving the quality of the elements.

It can however introduce

• "long" or "short" edges
• a poor representation of the geometry that will be difficult to recovered later
• inconsistent tagging

These 3 criteria have to be assessed to determine if a swap operation is accepted or not

Smoothing

The smoothing operation aims at improving the quality of the elements by moving vertices to some average of the locations of its neighbors.

In order to have a consistent smoothing on the boundaries of the computational domain, only the neighbors tagged on the same topological entity or one of its children are considered for smoothing. A projection step is still required for boundary vertices to ensure the consistency with the CAD model.

TODO: smoothing algorithms

Remeshing algorithm

The remeshing algorithm is iterative, and consists is multiple passes over the mesh edges / vertices to apply collapse, split, swap and smoothing.

The detail of the algorithm (order of the steps, number of iterations, criteria for acceptable operation, ...) controls the quality of remeshing as well as computational cost.

While most remeshing libraries share the same basic steps, they use different empirical rules to build the overall algorithm leading to significantly different results.

Remeshing quality

Enforcing mesh validity at each remeshing step ensures that the output of the remeshing algorithm is a valid mesh.

Given the remeshing algorithm (split and collapse operations), the target is not to reach a final mesh where all edges are of unit length in metric space, but only that these lengths lie in the interval $(1/\sqrt{2},\sqrt{2})$

The quality of remeshing is usually considered in terms of the histograms of edge lengths and element qualities.