A comparative numerical study of meshing functionals for variational mesh adaptation^††thanks: Supported in part by NSF (U.S.A.) through grant DMS-1115118 and NSERC (Canada) through grant A8781.

Variational mesh adaptation is an important type of mesh adaptation method and has received considerable attention from scientists and engineers; e.g., see the books [15, 19, 23, 24] and references therein. It also serves as the base of a number of commonly used adaptive moving mesh methods (e.g., see [5, 12, 14, 22]). In the variational approach, an adaptive mesh is generated as the image of a reference mesh under a coordinate transformation and such a coordinate transformation is determined as a minimizer of a certain meshing functional. A number of meshing functionals have been developed in the past (cf. the above mentioned books). For example, Winslow [25] proposed an equipotential method based on variable diffusion. Brackbill and Saltzman [3] developed a method by combining mesh concentration, smoothness, and orthogonality. Dvinsky [6] used the energy of harmonic mappings as his meshing functional, while Knupp [20] and Knupp and Robidoux [21] developed functionals based on the idea of conditioning the Jacobian matrix of the coordinate transformation. More recently, Huang [7] and Huang and Russell [15] proposed functionals based on the so-called equidistribution and alignment conditions.

With variational mesh adaptation, the mesh concentration is typically controlled through a scalar or a matrix-valued function, often referred to as the metric tensor or monitor function and defined based on some error estimates and/or physical considerations. While most of the meshing functionals have been developed with physical or geometric intuitions and have various levels of success in the adaptive numerical solution of partial differential equations (PDEs) and other applications, there is only a limited understanding on how the metric tensor affects the behavior of the mesh. An attempt to alleviate this lack of understanding was made by Cao et al. [4] for a generalization of Winslow’s variable diffusion functional. They showed that a significant change in an eigenvalue of the metric tensor along the corresponding eigendirection (first increasing and then decreasing, or vice versa) will result in adaptation of coordinate lines along this direction, although this adaptation competes with far more complicated effects, including those from changes in eigenvectors and other eigenvalues and the effects of the shapes of the physical and computational domains and the mesh point distribution on the boundaries. In [7, 15], two functionals have been developed based directly on the equidistribution and alignment conditions. These two conditions provide a complete characterization of the mesh elements through the metric tensor [7]. Minimizing the functionals leads to meshes which tend to satisfy the conditions in an integral sense. Nevertheless, this characterization is only qualitative, and how closely the mesh satisfies the conditions depends on the boundary correspondence between the computational and physical domains and the mesh point distribution on the boundaries. Thus, numerical studies, especially comparative ones, are useful, and often necessary, in understanding how the mesh adaptation for those meshing functionals is controlled precisely by the metric tensor. There do exist a few comparative numerical studies for meshing functionals. For example, a gallery of (adaptive and non-adaptive) meshes is given in [19] for a number of meshing functionals. Some comparative meshes are given in [15] for the harmonic mapping functional [6] and the subsequent functional based on equidistribution and alignment [7].

The main objective of this work is to present a comparative study for three of the most appealing meshing functionals, a generalization of Winslow’s variable diffusion functional (cf. Eq.˜6) and two functionals based on equidistribution and alignment (cf. Eqs.˜11 and 13). They are selected because Eqs.˜6 and 11 have been known to work well for many problems (e.g., see [1, 2, 7, 13, 14, 22]) while Eq.˜13 is similar to Eq.˜11 and has some very nice theoretical properties (cf. Section˜3.2). Another motivation is to present some three dimensional numerical results for those functionals. Critical for our study is to perform the substantial computations using the improved implementation of the variational methods introduced in [11]. In a sharp contrast to the situation in two dimensions, very little work has been done with variational mesh adaptation and adaptive moving mesh methods in three dimensions (e.g., see [15, 22]). It is particularly interesting to see how the functionals perform in this case.

An outline of the paper is given as follows. We describe the basic ideas of the variational mesh adaptation and its direct discretization (that is, first to discretize and then optimize) in Section˜2. In Section˜3 we introduce the three functionals to be studied for the numerical comparison, a generalization of Winslow’s variable diffusion functional and two functionals based on equidistribution and alignment. Numerical results and example adaptive meshes are given in Section˜4, followed by conclusions in Section˜5.

In variational mesh adaptation, an adaptive mesh is generated as the image of a reference mesh under a coordinate transformation. Denote the physical domain by $\Omega\subset\mathbb{R}^{d}$ ($d\geq 1$), and assume that we are given a computational domain $\Omega_{c}\subset\mathbb{R}^{d}$ and a quasi-uniform mesh $\hat{\mathcal{T}_{h}}_{c}$ thereon (in this work we consider only simplicial meshes). In many situations we can choose $\Omega_{c}$ to be the unit square/cube or simply $\Omega$. Denote the coordinate transformation by $\boldsymbol{x}=\boldsymbol{x}(\boldsymbol{\xi})\colon\Omega_{c}\to\Omega$ and its inverse by $\boldsymbol{\xi}=\boldsymbol{\xi}(\boldsymbol{x})\colon\Omega\to\Omega_{c}$. Such a coordinate transformation (more precisely, its inverse) is determined as a minimizer of a meshing functional. Most of the existing meshing functionals can be cast in a general form as

where $G$ is a smooth function, $\mathbb{J}$ is the Jacobian matrix of $\boldsymbol{\xi}=\boldsymbol{\xi}(\boldsymbol{x})$, $\det(\mathbb{J})$ denotes the determinant of $\mathbb{J}$, and $\mathbb{M}=\mathbb{M}(\boldsymbol{x})$ is the metric tensor supplied by the user to control mesh concentration. We assume that $\mathbb{M}=\mathbb{M}(\boldsymbol{x})$ is a symmetric and uniformly positive definite $d$-by-$d$ matrix-valued function on $\Omega$. Notice that Eq.˜1 is formulated in terms of the inverse coordinate transformation. One reason for this is that this form is less likely to produce singular meshes [6]. Another reason is that $\mathbb{M}$ is a function of $\boldsymbol{x}$ and thus finding the functional derivative of $I[\boldsymbol{\xi}]$ will not directly involve the derivatives of $\mathbb{M}$. This is convenient since in practice $\mathbb{M}$ is known only at the vertices of a mesh and its derivatives are not cheap to find. The main disadvantage of the formulation in this form is that $\boldsymbol{\xi}=\boldsymbol{\xi}(\boldsymbol{x})$ (or its numerical approximation) does not give the physical mesh directly. This is remedied either by interchanging the roles of the independent and dependent variables in the Euler-Lagrange equation of $I[\boldsymbol{\xi}]$ (e.g., see [15]) or, in a recently developed implementation (see below), computing the new physical mesh from a computational one using linear interpolation.

A minimizer of Eq.˜1 can be found numerically in the MMPDE (moving mesh PDE) framework. A conventional implementation [15] is to find the functional derivative of Eq.˜1 and then define the MMPDE as the gradient flow equation of the functional. Having been transformed by interchanging the roles of the dependent and independent variables, the MMPDE can be discretized on $\hat{\mathcal{T}_{h}}_{c}$ and a system of equations for the nodal velocities is obtained. Finally, the new mesh is obtained by integrating the mesh equation over a time step.

A much simpler implementation was proposed recently by Huang and Kamenski [11]. Instead of utilizing the MMPDE directly, the new implementation first discretizes the functional on the current mesh $\mathcal{T}_{h}$ and then, following the idea of the MMPDE, defines the mesh equation as the gradient equation of the discretized functional (with respect to the computational coordinates of the vertices). To be specific, denote by $\boldsymbol{x}_{i}$, $\hat{\boldsymbol{\xi}}_{i}$, and $\boldsymbol{\xi}_{i}$, $i=1,\dotsc,N_{v}$ the coordinates of the vertices of the current physical mesh ($\mathcal{T}_{h}$), the reference mesh ($\hat{\mathcal{T}_{h}}_{c}$), and the computational mesh (${\mathcal{T}_{h}}_{c}$), respectively. We assume that these meshes have the same numbers of the elements and vertices and the same connectivity. For any element $K\in\mathcal{T}_{h}$ (with vertices $\boldsymbol{x}_{i}^{K}$, $i=0,\dotsc,d$), the corresponding element in ${\mathcal{T}_{h}}_{c}$ is denoted by $K_{c}$ (with vertices $\boldsymbol{\xi}_{i}^{K}$, $i=0,\dotsc,d$). The edge matrices for $K$ and $K_{c}$ are defined as

Let $\omega_{i}$ be the element patch associated with vertex $\boldsymbol{x}_{i}$ (i.e., the collection of the elements containing $\boldsymbol{x}_{i}$ as a vertex). Then, the equation for the nodal velocities reads as

where $|K|$ is the volume of $K$, $\boldsymbol{v}_{i_{K}}^{K}$ is the local mesh velocity associated with vertex $\boldsymbol{x}_{i}$ in $K$, $i_{K}$ denotes the local index of $\boldsymbol{x}_{i}$ in $K$, $\tau>0$ is a constant parameter used to adjust the time scale of mesh movement, and $P=(P_{1},\dotsc,P_{N_{v}})$ is a positive function used to make the mesh equation to have desired invariance properties. Although the parameter $\tau$ can be absorbed in $P$, using two parameters has the advantage that the role of each parameter is clear: $\tau$ for the time scale of mesh movement while $P$ for the invariance properties. Ideally, $\tau$ should be chosen such that the mesh movement has the same scale as the physical equation. Unfortunately, there is no theoretical analysis for this yet and trial and error is still the most practical way to choose $\tau$. Numerical experience shows that a value in the range $[0.01,0.1]$ seems to work well for most problems [15]. In our computation, we choose $P$ such that the equation is invariant under the scaling transformation $\mathbb{M}\to c\,\mathbb{M}$ for all non-zero constants $c$ (cf. Section˜3): in variational mesh adaptation it is the relative distribution of $\mathbb{M}$ over the physical domain (instead of its absolute distribution) that determines the variation of the mesh density and therefore it is essential for the moving mesh equation to be invariant under the scaling transformation of $\mathbb{M}$.

The local velocities are given by

where the derivatives of $G$ with respect to $\mathbb{J}$ and $\det(\mathbb{J})$ (see Eqs.˜7, 12 and 14 below) are evaluated as

The second equation in Eq.˜3 is an inherent property resulting directly from the differentiation of the discretized meshing functional; it states that the centroid of $K$ stays fixed if only the contribution from $K$ is taken into account.

The above mesh equation should be modified properly for boundary vertices. For example, if $\boldsymbol{\xi}_{i}$ is a fixed boundary vertex, we replace the corresponding equation by

When $\boldsymbol{\xi}_{i}$ is allowed to move on a boundary curve/surface represented by

then the mesh velocity $\frac{\partial{}\boldsymbol{\xi}_{i}}{\partial{}t}$ needs to be modified such that its normal component along the curve or surface is zero, i.e.,

Mesh equation Eq.˜2 defines the movement of the nodes of the computational mesh ${\mathcal{T}_{h}}_{c}$ starting from the reference mesh $\hat{\mathcal{T}_{h}}_{c}$ at $t_{n}$. The equation can be integrated in time to obtain the computational mesh at $t_{n+1}$. For notational simplicity, we denote the computational mesh at $t_{n+1}$ by ${\mathcal{T}_{h}}_{c}$ as well. Notice that the physical mesh $\mathcal{T}_{h}$ is fixed during the time integration from $t_{n}$ to $t_{n+1}$. Meshes ${\mathcal{T}_{h}}_{c}$ and $\mathcal{T}_{h}$ define a correspondence

The new physical mesh is computed by means of linear interpolation as

where $\hat{\mathcal{T}_{h}}_{c}$ is the reference mesh on $\Omega_{c}$. The computational mesh plays the role of an intermediate variable.

Recall that the mesh concentration in variational mesh adaptation is controlled through the metric tensor $\mathbb{M}=\mathbb{M}(\boldsymbol{x})$. Such a metric tensor can be defined based on physical or geometric considerations or some error estimates. For example, for the $L^{2}$ norm of the error of piecewise linear interpolation on simplicial meshes, the optimal metric tensor [9, 16] (also see [15, (5.192)]) is

where $H(u)$ is the Hessian of function $u$, $|H(u)|$ is the eigen-decomposition of $H(u)$ with the eigenvalues being replaced by their absolute values, and the regularization parameter $\alpha>0$ is chosen such that

In practical computation, $u$ is typically unknown, and only the approximations to its values at the vertices are available. For this reason (and even in situations where an analytical expression for $u$ is available), the Hessian in Eq.˜4 is replaced by an approximation obtained by a Hessian recovery technique from the nodal values of $u$ or the approximations of the nodal values of $u$. A number of such techniques are known to produce nonconvergent recovered Hessians from a linear finite element approximation (e.g., see Kamenski [17]). Nevertheless, it is shown by Kamenski and Huang [18] that a linear finite element solution of an elliptic BVP converges at a second order rate as the mesh is refined if the recovered Hessian used to generate the adaptive mesh satisfies a closeness assumption. Numerical experiment shows that this closeness assumption is satisfied by the approximate Hessian obtained with commonly used Hessian recovery methods. We use a Hessian recovery method based on a least squares fit: a quadratic polynomial is constructed locally for each vertex via least squares fitting to neighboring nodal function values and an approximate Hessian at the vertex is then obtained by differentiating the polynomial.

Here we introduce the three meshing functionals used in the numerical study. A generalization of Winslow’s variable diffusion functional and the two functionals based on equidistribution and alignment are selected because they are reasonably simple, have nice theoretical properties, and are known to work well for many problems.

In the following we consider a number of examples in two and three dimensions. For a given function we consider $\mathbb{M}$ defined in Eq.˜4 which is optimal for minimizing the $L^{2}$ norm of the linear interpolation error of this function and compare meshes obtained from using the functionals of Winslow Eq.˜6 (W), Huang Eq.˜11 (H), and Huang and Russell Eq.˜13 (HR).

To assess the quality of the generated meshes, we compare the $L^{2}$ norm of the linear interpolation error and the equidistribution and alignment mesh quality measures, which describe how far the mesh is from being uniform in the metric defined by $\mathbb{M}$. The element-wise quality measures are based on Eq.˜8 and Eq.˜9 and defined as

while for the overall mesh quality measures we take their root-mean-squared values,

If the mesh is uniform with respect to $\mathbb{M}$, then $Q_{eq}=Q_{ali}=1$; if the mesh is far from being uniform with respect to $\mathbb{M}$, then $Q_{eq}$ and $Q_{ali}$ will become large. In other words, these quality measures describe how well the volume (measured by $Q_{eq}$) and the shape and orientation (measured by $Q_{ali}$) of mesh elements correspond to the desired size and shape prescribed by $\mathbb{M}$ (see [8] for more details on the mesh quality measures).