Target-Matrix Optimization Paradigm (TMOP)

• TMOP is a variational framework for high-order mesh optimization that precisely controls element geometry using spatially varying target matrices.
• It employs modular quality metrics and adaptive control mechanisms—supporting r-, h-, and hr-adaptivity—to efficiently fine-tune mesh elements based on application-specific criteria.
• The paradigm leverages high-performance techniques including matrix-free assembly and GPU acceleration to significantly speed up complex simulations and finite element analyses.

The Target-Matrix Optimization Paradigm (TMOP) is a variational framework for high-order mesh optimization. TMOP enables precise and flexible control over mesh element geometry by encoding desired sizes, shapes, aspect ratios, skewness, and orientations through a spatially varying target matrix. The approach is algebraic, element-local, and compatible with arbitrary polynomial order and element topology, making it extensively applicable in mesh quality control, adaptive mesh refinement, interface fitting, and PDE-driven mesh movement for both 2D and 3D finite element applications (Barrera et al., 2022, Camier et al., 2022, Dobrev et al., 2020, Dobrev et al., 2020, Dobrev et al., 2018).

1. Mathematical Foundation and General Formulation

A high-order mesh in TMOP is defined by a set of curved elements, each constructed via an isoparametric map $\Phi_E:\bar E\to\mathbb{R}^d$:

$$x(\bar x) = \Phi_E(\bar x) = \sum_{i=1}^{N_p} x_{E,i}\, \bar w_i(\bar x)$$

where $\{x_{E,i}\}$ are the control-point coordinates and $\{\bar w_i\}$ are basis functions on the reference element $\bar E$. The local Jacobian at a point $\bar x$ is

$$A(\bar x) = \frac{\partial \Phi_E}{\partial \bar x} = \sum_{i=1}^{N_p} x_{E,i} [\nabla \bar w_i(\bar x)]^T.$$

At every quadrature node, the user supplies a target Jacobian matrix $W(\bar x)\in \mathbb{R}^{d \times d}$ encoding ideal geometry. The core deviation transformation is

$T(\bar x) = A(\bar x) W(\bar x)^{-1}$

A scalar distortion metric $\mu(T)\ge 0$ quantifies the deviation from the ideal; $\mu(T)=0$ iff $T=I$ (i.e., the mesh matches the target exactly).

The general mesh optimization functional is

$$F(x) = \sum_{E \in \mathcal{M}} \int_{E_t} \mu\big(T(x_t)\big) dx_t$$

where $E_t$ is the “target” realization of element $E$, and integrations are performed using sufficient-order quadrature.

2. Target Matrix Construction and Adaptive Control

TMOP’s adaptability arises from the construction of $W(x)$ at each quadrature location. Following Knupp’s decomposition, $W$ can synthesize scalar size $\zeta$, arbitrary rotation $R\in SO(d)$, symmetric positive-definite skew $Q$, and diagonal aspect ratio $D$:

$W = \zeta R Q D$

Adaptive control is realized by driving $W(x)$ from problem-dependent indicators—solution gradients, shock positions, material interfaces, etc.—often encoded as finite element fields. For mesh anisotropy near interfaces,

$$W(x) = \text{vol}^{1/d} \,\times\, R(\theta) \,\times\, \text{skew}(\phi) \,\times\, \text{aspect ratio}(r)$$

where the fields may be reconstructed from, e.g., $\nabla\eta(x)$ (interface indicator), user-defined orientation, and prescribed skewness.

This construction is updated at each nonlinear solve iteration, including necessary interpolation or advection when the computational mesh changes (Dobrev et al., 2020).

3. Choice of Quality Metrics and Objective Functional

TMOP supports modular selection of distortion metrics $\mu(T)$:

• Shape metric (rotation and volume invariant): $\mu_2(T) = \frac{\|T\|_F^2}{2\,\det(T)} - 1$
• Size metric (element volume): $\mu_{77}(T) = \frac12 (\det T - \frac1{\det T})^2$
• Shape+Size metric: $$\mu_{80}(T) = \gamma\,\mu_2(T) + (1-\gamma)\,\mu_{77}(T)$$

Multiple metrics can be mixed with spatially-varying weights:

$$J = \sum_{E \in \mathcal{M}} \sum_{q} \sum_{i=1}^m \alpha_i(\bar x_q) \det(W_E^{(i)}(\bar x_q)) \mu_i(T_E^{(i)}(\bar x_q))$$

Metric normalization is optionally employed for invariance under mesh refinement or scaling, ensuring that each metric begins at unit value (Dobrev et al., 2020).

4. Optimization Procedure: Variational Minimization and Solvers

TMOP minimizes the global objective by varying all nodal positions $x$ (and, in adaptive cases, mesh topology). Gradients and Hessians of $F(x)$ are derived element- and quadrature-locally via careful chain-rule application:

$$\frac{\partial F}{\partial x_{a,i}} = \sum_{E, q} w_q \det W \sum_{n=1}^d \frac{\partial \mu}{\partial T_{a n}} \frac{\partial w_i(\xi_q)}{\partial \xi_n}$$

This produces block-sparse algebraic systems suitable for Newton-type methods, typically solved via global line search or trust-region techniques. Krylov solvers (e.g., MINRES with Jacobi preconditioning) are employed for the linearized systems, with globalization strategies designed to ensure energy decrease, element orientation preservation ($\det A > 0$), and (in interface fitting) geometric accuracy (Barrera et al., 2022, Camier et al., 2022).

5. Mesh Adaptivity: $r$-, $h$-, and $hr$-Variants

TMOP naturally supports:

• $r$-adaptivity: Nodal movement to minimize $F(x)$ towards encoded targets.
• $h$-adaptivity: Nonconforming mesh refinement/derefinement driven by local energy reduction indicators:

$$\Delta F_E^\gamma = F_E^{\gamma=0} - \frac{1}{N_c} \sum_{i=1}^{N_c} F_{E_i^{(\gamma)}}$$

An element is refined if $\max_\gamma \Delta F_E^\gamma > 0$.

• $hr$-adaptivity: Alternating $r$- (node movement) and $h$- (topological adjustment) steps until no further $h$-changes occur, accelerating convergence to geometric targets and reducing degrees of freedom compared to pure $h$- or $r$-adaptivity (Dobrev et al., 2020).

Metrics guide anisotropic vs. isotropic refinements: pure size metrics permit only isotropic patterns, pure shape metrics require anisotropic patterns, and combined metrics select the refinement pattern that maximally decreases TMOP energy.

6. Interface and Boundary Fitting: Level-Set Augmentation

TMOP is extended for implicit boundary/interface morphing by augmenting the objective functional with a penalty enforcing alignment to zero-isocontours of a level-set function $\sigma(x)$:

$$F(x) = F_\mu(x) + F_\sigma(x),\quad F_\sigma(x) = w_\sigma \sum_{s \in S} \sigma(x_s)^2$$

$S$ is a set of nodes/faces marked for fitting (e.g., those adjacent to “fictitious” material labels). Adaptive penalty weights are tuned per Newton iteration to ensure geometric accuracy without over-constraining:

$$w_{k+1,\sigma} = \alpha_\sigma w_{k,\sigma} \quad \text{if relative error reduction falls below threshold}$$

Level-set functions may reside on finer or nonconforming source meshes; high-order interpolation transfers $\sigma(x)$ and its gradients to the computational mesh at each iteration. All derivatives for penalty terms are closed-form, enabling efficient Newton solves (Barrera et al., 2022).

7. Computational Aspects and High-Performance Implementations

TMOP’s discrete objective and gradient assembly is recast in standard finite element operator form, fully compatible with modern software infrastructures and hardware acceleration:

• Partial assembly and tensor-product factorization exploit locality and reduce computational complexity from $O(p^{2d})$ to $O(p^{d+1})$.
• Matrix-free global assembly enables efficient GPU mapping using libraries such as MFEM: local quadrature operations, residual assembly, and Hessian actions are fused in GPU kernels with optimized memory access and register usage.
• Significant speed-ups (30$\times$–40$\times$ for V100 GPUs at $p\geq 2$) have been reported for high-order mesh optimization (Camier et al., 2022).

8. Illustrative Applications and Impact

TMOP underpins a wide range of mesh optimization challenges:

• Robust morphing of meshes to fit explicit or implicit boundaries/interfaces (2D circle, 3D sphere, complex CSG domains).
• Simulation-driven mesh adaptation in multi-material Arbitrary Lagrangian-Eulerian (ALE) hydrodynamics, sharply resolving interfaces and supporting large Lagrangian steps between remeshes.
• Coupling with automatic triggers for remap in ALE simulations to balance interface accuracy and computational cost:
  • Size-adaptive with automatic triggers yielded the lowest interface error at moderate remap cost (Dobrev et al., 2020).
• Preconditioning for shape optimization, accurately morphing initial meshes to target configurations before PDE-constrained optimization (Barrera et al., 2022).
• $hr$-adaptivity achieves similar solution accuracy with significantly fewer degrees of freedom compared to pure $h$- or $r$-adaptivity (Dobrev et al., 2020).

No explicit topological splitting is required except optional local conforming splits (e.g., for quad meshes with multiple faces marked for alignment); all mesh operations reuse the existing high-order finite element infrastructure (Barrera et al., 2022).

9. Summary Table: TMOP Component Features

Component	Description	Paper Reference
Target matrix $W$	Size, shape, aspect, orientation encoding	(Camier et al., 2022, Dobrev et al., 2018)
Quality metrics	Shape/size invariants, combined penalties	(Barrera et al., 2022, Dobrev et al., 2020)
Level-set penalty	Alignment to implicit interfaces	(Barrera et al., 2022)
$hr$-adaptivity	Combined node movement and refinement	(Dobrev et al., 2020)
GPU acceleration	Matrix-free, tensor contractions, partial assembly	(Camier et al., 2022)

TMOP constitutes a highly general, modular, and robust mesh optimization paradigm, directly enabling high-fidelity geometric control in fast-moving simulation contexts and advanced multi-material finite element workflows.