A surface moving mesh method based on equidistribution and alignment
Avary Kolasinski, Weizhang Huang
TL;DR
This paper introduces a novel surface moving mesh method based on equidistribution and alignment conditions, capable of handling general surfaces without relying on surface parameterization, and demonstrates its robustness through various examples.
Contribution
It extends the moving mesh PDE method to general surfaces without explicit parameterization, providing a robust, easy-to-implement approach using surface normals.
Findings
Method works well for 2D and 3D surface examples.
Mesh trajectories remain nonsingular if initially so.
Applicable to surfaces with or without explicit parameterization.
Abstract
A surface moving mesh method is presented for general surfaces with or without explicit parameterization. The method can be viewed as a nontrivial extension of the moving mesh partial differential equation method that has been developed for bulk meshes and demonstrated to work well for various applications. The main challenges in the development of surface mesh movement come from the fact that the Jacobian matrix of the affine mapping between the reference element and any simplicial surface element is not square. The development starts with revealing the relation between the area of a surface element in the Euclidean or Riemannian metric and the Jacobian matrix of the corresponding affine mapping, formulating the equidistribution and alignment conditions for surface meshes, and establishing a meshing energy function based on the conditions. The moving mesh equation is then defined as…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 18
Figure 19
Figure 20
Figure 21
Figure 22
Figure 23
Figure 24
Figure 25
Figure 26
Figure 27
Figure 28
Figure 29
Figure 30
Figure 31
Figure 32
Figure 33
Figure 34
Figure 35
Figure 36
Figure 37
Figure 38
Figure 39
Figure 40Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
A surface moving mesh method based on equidistribution and alignment
Avary Kolasinski Department of Mathematics, the University of Kansas, Lawrence, KS 66045 ([email protected]).
Weizhang Huang Department of Mathematics, the University of Kansas, Lawrence, KS 66045 ([email protected]).
Abstract
A surface moving mesh method is presented for general surfaces with or without explicit parameterization. The method can be viewed as a nontrivial extension of the moving mesh partial differential equation method that has been developed for bulk meshes and demonstrated to work well for various applications. The main challenges in the development of surface mesh movement come from the fact that the Jacobian matrix of the affine mapping between the reference element and any simplicial surface element is not square. The development starts with revealing the relation between the area of a surface element in the Euclidean or Riemannian metric and the Jacobian matrix of the corresponding affine mapping, formulating the equidistribution and alignment conditions for surface meshes, and establishing a meshing energy function based on the conditions. The moving mesh equation is then defined as the gradient system of the energy function, with the nodal mesh velocities being projected onto the underlying surface. The analytical expression for the mesh velocities is obtained in a compact, matrix form, which makes the implementation of the new method on a computer relatively easy and robust. Moreover, it is analytically shown that any mesh trajectory generated by the method remains nonsingular if it is so initially. It is emphasized that the method is developed directly on surface meshes, making no use of any information on surface parameterization. It utilizes surface normal vectors to ensure that the mesh vertices remain on the surface while moving, and also assumes that the initial surface mesh is given. The new method can apply to general surfaces with or without explicit parameterization since the surface normal vectors can be computed even when the surface only has a numerical representation. A selection of two- and three-dimensional examples are presented.
AMS 2010 Mathematics Subject Classification. 65M50, 65N50
Key Words. surface mesh movement, surface mesh adaptation, moving mesh PDE, mesh nonsingularity, surface parameterization
Abbreviated title. A surface moving mesh method.
1 Introduction
We are interested in methods that can directly move simplicial meshes on general surfaces with or without analytical expressions. Such surface moving mesh methods can be used for adaptation and/or quality improvements of surface meshes and thus are useful for computational geometry and numerical solutions of partial differential equations (PDEs) defined on surfaces; e.g., see [7, 10, 25].
There has been some work done on mesh movement and adaptation for surfaces. For example, Crestel et al. [5] present a moving mesh method for parametric surfaces by generalizing Winslow’s meshing functional to Riemannian manifolds and taking into consideration the Riemannian metric associated with the manifolds. The method is simplified and implemented on a two-dimensional domain for surfaces that accept certain parameterizations. Weller et al. [3] and McRae et al. [24] solve a Monge-Ampére type equation on the surface of the sphere to generate optimally transported meshes that become equidistributed with respect to a suitable monitor function. MacDonald et al. [23] devise a moving mesh method for the numerical simulation of coupled bulk-surface reaction-diffusion equations on an evolving two-dimensional domain. They use a one-dimensional moving mesh equation in arclength to concentrate mesh points along the evolving domain boundary. Dassi et al. [6] generalize the higher embedding approach proposed in [1]. They modify the embedding map between the underlying surface and to include more information associated with the physical solution and its gradient. The idea behind this mapping is that it essentially approximates the geodesic length on the surface via a Euclidean length in . The mesh adapts in the Euclidean space and then is mapped back to the physical domain.
The objective of this paper is to present a surface moving mesh method for general surfaces with or without explicit parameterization. The method can be viewed as a nontrivial extension of the moving mesh PDE (MMPDE) method that has been developed for bulk meshes and demonstrated to work well for various applications; e.g. see [19, 20, 21]. The main challenges in the development of surface mesh movement come from the fact that the Jacobian matrix of the affine mapping between the reference element and any simplicial surface element is not square. To overcome these challenges, we start by connecting the area of the surface element in the Euclidean metric or a Riemannian metric with the Jacobian matrix. This connection allows us to formulate the equidistribution and alignment conditions and ultimately, form a meshing energy function for surface meshes. This meshing function is similar to a discrete version of Huang’s functional [13, 17, 22] for bulk meshes which has been proven to work well in a variety of problems. Following the MMPDE approach, we define the surface moving mesh equation as the gradient system of the meshing function, with the nodal mesh velocities being projected onto the underlying surface. The analytical expression for the mesh velocities is obtained in a compact, matrix form, which makes the implementation of the new method on a computer relatively easy and robust. Several theoretical properties are obtained for the surface moving mesh method. In particular, it is proven that a surface mesh generated by the method remains nonsingular for all time if it is so initially. Moreover, the element altitudes and areas of the physical mesh are bounded below by positive constants depending only on the initial mesh, the number of elements, and the metric tensor that is used to provide information on the size, shape, and orientation of the elements throughout the surface. Furthermore, limiting meshes exist and the meshing function is decreasing along each mesh trajectory. These properties are verified in numerical examples.
It is emphasized that the new method is developed directly on surface meshes, making no use of any information on surface parameterization. It utilizes surface normal vectors to ensure that the mesh vertices remain on the surface while moving, and also assumes that the initial surface mesh is given. Since the surface normal vectors can be computed even when the surface only has a numerical representation, the new method can apply to general surfaces with or without explicit parameterization. A selection of two- and three-dimensional examples are presented.
This paper is organized as follows. In Section 2, the area formula for a surface element and the equidistribution and alignment conditions for surface meshes are established. The surface moving mesh equation is described in Section 3 and its theoretical analysis is given in Section 4. Numerical examples are then provided in Section 5 followed by conclusions and further remarks in Section 6. Appendix A contains the derivation of derivatives of the meshing function with respect to the physical coordinates.
2 Equidistribution and alignment for surface meshes
In this section we formulate the equidistribution and alignment conditions for a surface mesh. These conditions are used to characterize the size, shape, and orientation of the elements and develop a meshing function for surface mesh generation and adaptation. The function is similar to the one [13, 16] used in bulk mesh generation and adaptation and also based on mesh equidistribution and alignment.
2.1 Area and affine mappings for surface elements
Let be a bounded surface in (). Assume that we have a mesh on and let and be the number of its elements and vertices, respectively. The elements are surface simplexes in , i.e., they are -dimensional simplexes in a -dimensional space. Notice that their area in dimensions is equivalent to their volume in dimensions. Assume that the reference element has been chosen to be a -dimensional equilateral and unitary simplex in a -dimensional space. For and any element let be the affine mapping between them and be its Jacobian matrix. Denote the vertices of by \mbox{\boldmathx}_{j}^{K}\in\mathbb{R}^{d}, and the vertices of by \mbox{\boldmath\xi}_{j}\in\mathbb{R}^{d-1}, . Then
[TABLE]
From this, we have
[TABLE]
or
[TABLE]
which gives where and are the edge matrices for and , i.e.,
[TABLE]
Notice that is a square matrix and its inverse exists since is not degenerate. However, unlike the bulk mesh case, matrices are not square. This makes the formulation of adaptive mesh methods more difficult for surface than bulk meshes. Nevertheless, the approach is similar for both situations, as will be seen below.
In the following we can see that the area of the physical element can be determined using or .
** Lemma 2.1****.**
For any surface simplex , there holds
[TABLE]
where and denote the area of the simplexes and , respectively, and denotes the determinant of a matrix.
Proof.
From , we have
[TABLE]
where we have used . Let the QR-decomposition of be given by
[TABLE]
where is a unitary matrix, is an upper triangular matrix, and \mbox{\boldmath0}\in\mathbb{R}^{1\times(d-1)} is a row vector of zeros. This decomposition indicates that is formed by rotating the convex hull with edges formed by the column vectors of \left[\begin{matrix}R_{K}\\ \mbox{\boldmath0}\end{matrix}\right]. We have
[TABLE]
where we have used the fact that rotation, , does not change the area. Since the convex hull formed by the column vectors of \left[\begin{matrix}R_{K}\\ \mbox{\boldmath0}\end{matrix}\right] lies on the \mbox{\boldmathx}^{(1)} – – \mbox{\boldmathx}^{(d-1)} – plane, its area is equal to the -dimensional volume of the convex hull formed by the column vectors of in -dimensions. Then,
[TABLE]
On the other hand, we have
[TABLE]
Therefore,
[TABLE]
∎
2.2 Area of surface elements in a Riemannian metric
We now formulate the area of a surface element in a Riemannian metric using or . The formula is needed later in the development of algorithms for mesh adaptation. First we consider a symmetric, uniformly positive definite metric tensor \mathbb{M}(\mbox{\boldmathx}) which satisfies
[TABLE]
where and are positive constants, is the identity matrix, and the less-than-or-equal-to sign is in the sense of negative semi-definiteness. We define the average of over as
[TABLE]
Recall that the distance in the Riemannian metric, , is defined as
[TABLE]
where denotes the standard Euclidean norm. This implies that the geometric properties of in the metric can be obtained from those of in the Euclidean metric.
** Lemma 2.2****.**
For any surface simplex , there holds
[TABLE]
where denotes the area of in the metric .
Proof.
The Jacobian matrix of the affine mapping from to is given by
[TABLE]
From the discussion following (3), we know that the area of in the metric is equal to the area of in the Euclidean metric. Thus, from Lemma 2.1 we have
[TABLE]
∎
The following lemma gives a lower bound for the area of with respect to the metric in terms of the minimum altitude of with respect to .
** Lemma 2.3****.**
Let denote the minimum altitude of with respect to . Then,
[TABLE]
Proof.
From Lemma 2.2 and , we have
[TABLE]
Let the -decomposition of be denoted as
[TABLE]
where is a unitary matrix, is an upper triangular matrix, and [math] is a -dimensional row vector of zeros. This gives
[TABLE]
where , denote the singular values of . Additionally, by [2, Lemma 5.12] we have that
[TABLE]
where denotes the minimum altitude of the simplex formed by the columns of . Combining these, we get
[TABLE]
Since is a rotational matrix, the minimum altitude of with respect to the metric is the same as the minimum altitude of the convex hull formed by the columns of i.e., . Thus, we have obtained (5). ∎
The relationship given in the above lemma between the area and minimum height will be used in the proof of the nonsingularity for surface meshes in §4. It is instructional to note that in two dimensions (), is a line segement and both and represent the length of in the metric and are equal. In this case, the inequality (5) reduces to
[TABLE]
which is very sharp. For , (5) becomes
[TABLE]
which is not as sharp as in two dimensions. Indeed, when is equilateral with respect to , we have [11]
[TABLE]
2.3 Equidistribution and alignment conditions
We can now define the equidistribution and alignment conditions characterizing a general nonuniform, simplicial surface mesh. Notice that any nonuniform mesh can be viewed as a uniform one in some metric tensor. Specifically, a mesh is uniform in some metric if all of the elements in the mesh have the same size and are similar to a reference element with respect to that metric. In this point of view, the equidistribution condition requires that all of the elements in the mesh have the same size. Mathematically, this can be expressed as
[TABLE]
where, as before, denotes the area of the surface with respect to the metric and
. Using Lemma 2.2 and recalling , we have
[TABLE]
Thus, the equidistribution condition (6) becomes
[TABLE]
The alignment condition, on the other hand, requires that all of the elements be similar to the reference element . Notice that any element is similar to if and only if is composed by dilation, rotation, and translation, or equivalently, is composed by dilation and rotation. Mathematically, can therefore be expressed as
[TABLE]
where is a constant representing dilation and and are orthogonal matrices representing rotation. Thus,
[TABLE]
It can be verified that the above equation is equivalent to
[TABLE]
which is referred to as the alignment condition. Here, denotes the trace of a matrix.
With these two conditions, we can now formulate the meshing function. To do so, we first consider the alignment condition (9) and note that an equivalent condition is
[TABLE]
Notice that the left- and right-hand sides are the arithmetic mean and the geometric mean of the eigenvalues of the matrix , respectively. The inequality of arithmetic and geometric means gives
[TABLE]
with equality if and only if all of the eigenvalues are equal. From (10), for any general mesh which does not necessarily satisfy (9), we have
[TABLE]
and therefore
[TABLE]
where is a dimensionless parameter which has been added to agree with the equidistribution energy function below. Then, we define the alignment energy function as
[TABLE]
whose minimization will result in a mesh that closely satisfies the alignment condition (9). One may notice that has been added as a weight.
Similarly, we consider the equidistribution condition (7). From Hölder’s inequality, for any we have
[TABLE]
with equality if and only if
[TABLE]
That is, minimizing the difference between the left-hand side and the right-hand side of (12) tends to make constant for all . Noticing that the left-hand side of (12) is simply , we can rewrite this inequality as
[TABLE]
We can consider constant since \sigma_{h}\approx\displaystyle\int_{S}\det(\mathbb{M}(\mbox{\boldmathx}))^{1/2}d\mbox{\boldmathx} and hence it only weakly depends on the mesh. Therefore, we define the equidistribution energy function as
[TABLE]
whose minimization will result in a mesh that closely satisfies the equidistribution condition.
2.4 Energy function for combined equidistribution and alignment
We now have two functions, one for each of equidistribution and alignment. Our goal is to formulate a single meshing function for which minimizing will result in a mesh that closely satisfies both conditions. One way to ensure this is to average (2.3) and (13), that is, define for . This leads to
[TABLE]
where and are dimensionless parameters, with the latter balancing the equidistribution and alignment conditions for which full alignment is achieved when and full equidistribution is achieved when . We can write (14) as
[TABLE]
where
[TABLE]
We remark that for and , is coercive; see the definition of coercivity in Section 4. As can be seen therein, coercivity is an important property when proving mesh nonsingularity. It is also instructional to point out that the function (14) is very similar to a Riemann sum of the meshing function developed in [13] for bulk meshes based on equidistribution and alignment. One of the main differences is that cannot be simplified in (14) since it is not a square matrix as it is in the bulk mesh case. Additionally, the constant terms and exponents that contain are in (14) instead of in the bulk mesh case. The functional of [13] has been proven to work well for a variety of problems [20].
3 Moving mesh equations for surface meshes
In this section we describe the surface MMPDE method used to find the minimizer of (14). In principle, we can directly minimize it. However, this minimization problem can be extremely difficult to solve since (14) is highly nonlinear in general. We employ here the MMPDE approach (a time transient approach) to find the minimizer and define the moving mesh equation as the gradient system of the meshing function.
3.1 Gradient of meshing energy
Motivated by the function (14), we consider meshing functions in a general form (15), i.e.,
[TABLE]
where is a given smooth function of
[TABLE]
that is, , and
[TABLE]
Indeed, a special example is (16) but can be chosen differently. Moreover, both and depend on the coordinates of the vertices of the physical element and hence is a function of them, i.e., can be expressed as
[TABLE]
where \mbox{\boldmathx}_{i}^{K}\in\mathbb{R}^{d} for are the coordinations of the vertices of . As a consequence, the sum in (17) is a function of the coordinates of all vertices of the physical mesh , i.e.,
[TABLE]
where \mbox{\boldmathx}_{i}\in\mathbb{R}^{d} for are the coordinates of the vertices of the mesh with global indices. One of the underlying keys to our approach is to find the derivatives of with respect to the physical coordinates \mbox{\boldmathx}_{1},\dots,\mbox{\boldmathx}_{N_{v}} which requires elementwise derivatives of with respect to \mbox{\boldmathx}_{1}^{K},\dots,\mbox{\boldmathx}_{d}^{K}. That is,
[TABLE]
where denotes the local index of vertex \mbox{\boldmathx}_{i} in and is the element patch associated with \mbox{\boldmathx}_{i}.
In order to calculate the necessary derivatives, we recall some definitions and properties of scalar-by-matrix differentiation (cf. [16] for details). Let be a scalar function of a matrix . Then the scalar-by-matrix derivative of with respect to is defined as
[TABLE]
The chain rule of differentiation with respect to is
[TABLE]
With this and when is a square matrix, the following properties have been proven in [16],
[TABLE]
Using the above, we can find the expressions for and which are needed to compute (21). For the function (14), the first derivatives of are given by
[TABLE]
3.2 Derivatives of the meshing function with respect to the physical coordinates
From (21), we can see that we will need \partial G_{K}/\partial\mbox{\boldmathx}_{i_{K}}^{K} to compute \partial I_{h}/\partial\mbox{\boldmathx}_{i}. The former can be obtained once we know the derivatives of with respect to the coordinates of all vertices of , i.e.,
[TABLE]
The derivation of these derivatives is given in Appendix A. They are given as
[TABLE]
where and are given in (25), is the linear basis function associated with \mbox{\boldmathx}_{j}^{K}, \mathbb{M}_{j,K}=\mathbb{M}(\mbox{\boldmathx}_{j}^{K}), and
[TABLE]
Having computed \partial G_{K}/\partial\mbox{\boldmathx}_{j}^{K} () for all elements using (26) and (27), we can obtain \partial I_{h}/\partial\mbox{\boldmathx}_{i} from (21).
3.3 Surface moving mesh equations
As mentioned above, we employ a surface MMPDE method to minimize the meshing function (14) or a more general form (17). An MMPDE is a mesh equation that involves mesh speed. There are various formulations of MMPDEs; we focus here on the approach where the surface MMPDE is defined as the modified gradient system of the meshing function. A distinct feature for surface meshes, other than bulk meshes, is that the nodes need to stay on the surface. By Section 3.2 we may assume that we have the matrix
[TABLE]
Let \Phi(\mbox{\boldmathx})=0 denote the surface, where can be defined through an analytical expression or a numerical representation such as by spline functions. Then for the vertices to stay on the surface we should have \Phi(\mbox{\boldmathx}_{i})=0,\;i=1,...,N_{v} or at least
[TABLE]
where \frac{d\mbox{\boldmathx}_{i}}{dt} is the nodal mesh velocity. Following the MMPDE approach, we would define the mesh equation as the gradient system of , i.e.,
[TABLE]
where is a positive scalar function used to make the equation have desired invariance properties and is a constant parameter used for adjusting the time scale of mesh movement. Obviously, this does not satisfy (30). Here we propose to project the velocities in (31) onto the surface and define the surface moving mesh equation as
[TABLE]
where \mbox{\boldmathn}_{i}={\nabla\Phi(\mbox{\boldmathx}_{i})}/{\|\nabla\Phi(\mbox{\boldmathx}_{i})\|} is the unit normal to the surface at \mbox{\boldmathx}_{i} and the difference inside the square bracket is the projection of the vector {\partial I_{h}}/{\partial\mbox{\boldmathx}_{i}} onto the tangential plane of the surface at \mbox{\boldmathx}_{i}. Notice that this surface MMPDE inherently ensures that (30) be satisfied or, in words, the nodes stay on the surface during the mesh movement. Moreover, it is important to note that (32) only utilizes the unit normal vectors of the surface whose computation does not require explicit parameterization or analytical expression of the surface. As mentioned, for surfaces without an analytical expression, spline functions may be used to approximate the gradient for (32). Although the numerical examples presented in this work have explicit parameterizations, it is a goal for our future work to study the use of spline functions to approximate the unit normal vectors for surfaces represented by simplicial background meshes.
Using (21) we can rewrite the above equation in a compact form as
[TABLE]
where \mbox{\boldmathv}_{i_{K}}^{K}\in\mathbb{R}^{d} is the local mesh velocities contributed by to \mbox{\boldmathx}_{i_{K}}^{K} and has the expressions
[TABLE]
and {\partial G_{K}}/{\partial\mbox{\boldmathx}_{i_{K}}^{K}} is given in (26) and (27).
The surface MMPDE (33) must be modified properly for boundary vertices when has a boundary. For fixed boundary vertices, the corresponding equation is replaced by
[TABLE]
The velocities for other boundary vertices should be modified such that they slide on the boundary.
With proper modification of the boundary vertices, the system (33) can be integrated in time. To do so, one first starts by calculating the edge matrices for all elements and for the reference element. One can then readily calculate (25) which is needed for (26) and (27). Then one can integrate (33) in time. For this work we use Matlab’s ODE solvers ode45 and ode15s. The explicit scheme, ode45, implements a 4(5)-order Runge-Kutta method with a variable time step. The implicit scheme, ode15s, is a variable time step and variable-order solver based on the numerical differentiation formulas of orders 1 to 5. All of the numerical examples in this paper use ode45 although both ode45 and ode15s have been tested and proven to work very well in computation.
4 Nonsingularity of surface moving meshes
In this section we study the nonsingularity of the mesh trajectory and the existence of limiting meshes as for the MMPDE (33).
4.1 Equivalent measure of minimum height
We begin the theoretical analysis by establishing the relation between and the minimum altitude of with respect to .
** Lemma 4.1****.**
There holds
[TABLE]
where is the altitude of and is the minimum altitude of with respect to the metric .
Proof.
First of all, we have
[TABLE]
Now, consider the QR decomposition of
[TABLE]
where is a unitary matrix, is an upper triangular matrix, and [math] is a -dimensional row vector of zeros. With this we have
[TABLE]
By [16, Lemma 4.1] we have
[TABLE]
where is the minimum altitude of the simplex formed by the columns of . Since is a rotation matrix, is the same as , the minimum altitude of with respect to the Euclidean metric. Combining the above results, we get
[TABLE]
The inequality (35) follows from this and the observation that the geometric properties of with respect to the metric are the same as those of with respect to the Euclidean metric. ∎
Lemma 4.1 indicates that if is chosen to satisfy then
[TABLE]
4.2 Mesh nonsingularity
We now consider the MMPDE (33). Recall that the velocities for the boundary vertices need to be modified in order for them to stay on the boundary. However, the analysis is similar with or without modifications. Hence, for simplicity we do not consider modifications in the analysis. We also note that for theoretical purposes, we assume that is taken to satisfy instead of being unitary as we have been considering thus far. This change does not affect the actual computation. However, since typically we expect , the assumption will likely lead to and thus (the value of on the initial mesh ) stays . On the other hand, if (unitary), we have and will depend strongly on .
In the following analysis, the mesh at time is denoted by \mathcal{T}_{h}(t)=\left(\mbox{\boldmathx}_{1}(t),\dots,\mbox{\boldmathx}_{N_{v}}(t)\right).
** Theorem 4.1****.**
Assume that the meshing function in the form (17) satisfies the coercivity condition
[TABLE]
where , , and are constants. We also assume that is equilateral and . Then if the elements of the mesh trajectory of the MMPDE (33) have positive areas initially, they will have positive areas for all time. Moreover, their minimum altitudes in the metric and their areas in the Euclidean metric are bounded below by
[TABLE]
where
[TABLE]
Proof.
From (32) we have
[TABLE]
This implies for all . From coercivity (37) and Lemma 4.1, we get
[TABLE]
By Lemma 6, , thus
[TABLE]
and therefore
[TABLE]
Moreover, from the assumption that is equilateral and it follows that
[TABLE]
Combining (42) and (43) we get
[TABLE]
which gives (38).
Furthermore, we have
[TABLE]
Then (39) follows from the above inequality and (38).
Finally, from (26) and (27) it is not difficult to see that the magnitude of the mesh velocities is bounded from above when is bounded from below. As a consequence, the mesh vertices will move continuously with time and cannot jump over the bound to become negative. Hence, will stay positive if so initially. ∎
From the proof we have seen that the key points are the energy decreasing property and the coercivity of the meshing function. The former is satisfied by the MMPDE (33) by design while the latter is an assumption for the meshing function. We emphasize that the result holds for any function satisfying the coercivity condition (37).
On the other hand, the condition (37) is satisfied by the meshing function (14) for and (with and in Theorem 4.1). It is interesting to point out that the role of the parameter can be explained from (38). Indeed, for this case the inequality (38) becomes
[TABLE]
Since represents the average diameter of the elements, the above inequality implies that the mesh becomes more uniform as gets larger. We take , which has been found to work well for all examples we have tested.
** Theorem 4.2****.**
Under the assumptions of Theorem 4.1, for any nonsingular initial mesh, the mesh trajectory of MMPDE (33) has the following properties.
* has a limit as , i.e.,*
[TABLE] 2. 2.
The mesh trajectory has limiting meshes, all of which are nonsingular and satisfy (38) and (39). 3. 3.
The limiting meshes are critical points of , i.e., they satisfy
[TABLE]
Proof.
The proof is very much the same as that for [15, Theorem 4.3] for the bulk mesh case. The key ideas to the proof are the monotonicity and boundedness of and the compactness of . With these holding for the surface mesh case, one can readily prove the three properties. ∎
It is remarked that the above two theorems have been obtained for the MMPDE (33) which is semi-discrete in the sense that it is discrete in space and continuous in time. A fully discrete scheme can be obtained by applying a time-marching scheme to (33). Then similar results can be obtained for the fully discrete scheme under similar assumptions and for a sufficiently small but not diminishing time step. Since the analysis is similar to that for the bulk mesh case, the interested reader is referred to [15] for the detailed discussion.
5 Numerical examples
Here we present numerical results for a selection of two- and three-dimensional examples to demonstrate the performance of the surface moving mesh method described in the previous sections. The main focus will be on showing how our method can be used for mesh smoothing and concentration. To assess the quality of the generated meshes, we compare the equidistribution () and alignment () mesh quality measures which are defined as
[TABLE]
These measures are indications of how closely the mesh satisfies the equidistribution condition (7) and the alignment condition (9), respectively. The closer these quality measures are to 1, the closer they are to a uniform mesh with respect to the metric . It should be noted that the alignment condition does not apply to the two-dimensional case where a “surface” is actually a curve. Mathematically, when , is a number and hence (9) is always satisfied.
For all computations we use and in the meshing function (14). This choice has been known to work well in bulk mesh applications. Interestingly, we have found that it also works well for all surface mesh examples we have tested. We take , , and
[TABLE]
The latter is to ensure that the MMPDE (33) be invariant under scaling transformations of . For all of the results, we run to a final time of 1.0.
We choose two forms of . The first is , which will ensure the mesh move to become as uniform as possible with respect to the Euclidean norm. The second is a curvature-based metric tensor defined as a scalar matrix , where is the mean curvature and is machine epsilon. The mean curvature is defined (e.g., see [12]) for a curve in as
[TABLE]
and for a surface in as
[TABLE]
where
[TABLE]
We would like to explore more metric tensors in future work but will focus on these two for this paper.
* Example 5.1**.*
For the first example, we generate adaptive meshes for the unit circle in two dimensions,
[TABLE]
We take and fix the node \mbox{\boldmathx}_{1}=(1,0).
Fig. 1 shows the meshes for this example. Studying the figures we see that the initial mesh Fig. 1(a) is very nonuniform but the final meshes Fig. 1(b) and (c) have adapted to be equidistant along the curve. Moreover, the final meshes for both (Fig. 1 (b)) and (Fig. 1(c)) adapt the mesh in the same manner. This is consistent with the fact that the curvature of a circle is constant thus the nodes do not concentrate in one particular region of the curve. The final meshes in both cases provide good size adaptation and are more uniformly distributed along the curve when compared with the initial mesh. This can be further supported assessing the mesh quality measure for which improves from 7.509604 to 1.000004 for both cases of . The fact that indicates that the mesh is close to satisfying the equidistribution condition (7) and hence the mesh is almost uniform with respect to the metric tensor . It can also be seen that the nodes remain on the curve , which, as mentioned earlier, is an inherent feature of the new surface moving mesh method and indeed an important one when adapting a mesh on a curve.
* Example 5.2**.*
The second two-dimensional example is the ellipse defined by
[TABLE]
In this example we take and fix the node \mbox{\boldmathx}_{1}=(8,0).
The initial nodes (Fig. 2(a)) are randomly distributed through the curve. However, for , the final mesh (Fig. 2(b)) is equidistant along the ellipse providing a much more uniform mesh. This can also be seen in which improves from 5.497002 initially to 1.026912 in the final mesh.
Now considering the curvature-based metric tensor (Fig. 2(c)), we can see a high concentration of elements near the regions of the ellipse with large curvature. This is consistent with the equidistribution principle which requires higher concentration in the regions with larger determinant of the metric tensor (larger mean curvature in the current situation). The mean curvature is large in the regions of the ellipse close to and almost 0 for . From Fig. 2(c) we can see that the adaptation with provides a mesh that represents the shape of the curve much better than other two meshes. The improvement of from 5.126216 to 1.015848 indicates that the final mesh is almost uniform with respect to the curvature-based metric tensor.
* Example 5.3**.*
For the next two-dimensional example, we generate adaptive meshes for the sine curve defined by
[TABLE]
In this example we take and fix the end nodes \mbox{\boldmathx}_{1}=(0,0) and \mbox{\boldmathx}_{61}=(2\pi,0).
Fig. 3 shows the meshes for this example. From Fig. 3(a) and (b) we see that for , the mesh becomes much more uniform. This is consistent with the fact that for , the minimization of the meshing function will make the mesh more uniform with respect to the Euclidean norm. The observation can be further supported by assessing the mesh quality measures for which measure improves from 4.183312 to 1.002906 indicating that the final mesh satisfies the equidistribution condition (7) closely.
Now studying Fig. 3(c) where is used, we see that there is a high concentration of mesh elements in regions with large curvature, i.e., the hill at and cup at , which is consistent with the use of the curvature-based metric tensor. Moreover, the equidistribution measure improves from 6.254755 to 1.007493. This indicates that although the mesh may seem nonuniform in the Euclidean metric, it is almost uniform in the metric .
As discussed in Section 4, theoretically we know that the value of is decreasing and is bounded below. To see these numerically, we plot and as functions of in Fig. 4, where denotes the minimum area of over all elements in . The numerical results are shown to be consistent with the theoretical predictions. Specifically, for , Fig. 4(a) shows that is decreasing and bounded below by 9.535. Additionally, Fig. 4(b) suggests that is bounded below by 0.235 which is the value of of the initial mesh. As we see, first increases and then converges to about 0. The reason is because in the final mesh, the elements are close to being uniform with respect to the Euclidean metric and thus for all . Since the initial mesh is nonuniform, we expect an increase in as the mesh is becoming more uniform. Moreover, as the mesh reaches the limiting mesh trajectory around , we see that converges as shown in Fig. 4(b).
For the case with , the numerical results are again consistent with the theoretical predictions. Fig. 4(c) shows that is decreasing for all time and bounded below by 15.5. This figure also shows that at around , begins to converge. In Fig. 4 (d), has similar properties to Fig. 4(b). That is, we see an initial increase in after which, the value converges to 0.11 starting at around . Furthermore, Fig. 4(d) suggests that is bounded below by the initial value of 0.045.
* Example 5.4**.*
As the final two-dimensional example, we generate adaptive meshes for the lemniscate defined by
[TABLE]
In this example adapt the mesh on the curve for both and . In both situations, we fix the node \mbox{\boldmathx}_{1}=(2,0).
From Fig. 5 we see that for the mesh adapts from a very nonuniform initial mesh (Fig. 5(a)) to a uniform final mesh (Fig. 5(b)) when considering the metric tensor corresponding to the Euclidean metric. The nodes are equidistant apart while remaining on the curve. This improvement in uniformity can be further supported by the equidistribution quality measure which improves from 2.083287 for the initial mesh to 1.002549 for the final mesh.
We see a similar result when the curvature-based metric tensor is used (Fig. 5(c)). A higher concentration of nodes occurs in the circular regions with larger curvature compared to the cross section which has smaller curvature (i.e., the linear regions). It is not a significant difference in concentration but this is consistent with the fact that the curvature of the lemniscate is close to but not exactly constant. The equidistribution quality measure improves from 3.364232 to 1.001011 indicating that the final mesh is much more uniform with respect to the curvature-based than the initial mesh.
For , Fig. 6 shows similar findings. When considering the Euclidean metric, we see the mesh, Fig. 6(a), is very nonuniform initially and adapts to a equidistant spacing of the nodes along the curve. This is further supported in the quality measures for which the equidistribution measure improves from 2.552134 to 1.002167 indicating that the final mesh is close to satisfying the equidistribution condition.
The curvature-based metric tensor results in a similar adaptation as before. Fig. 6(c) shows the final mesh has adapted in such a way where there is a higher concentration of nodes in those regions of the curve with larger curvature, i.e., circular regions. Comparatively, there are fewer nodes in the cross section which has smaller curvature. The difference in concentration can be clearly seen in Fig. 6(c) with nodes. The adaptation is consistent with the curvature of the lemniscate, which is close to but not exactly constant. This improvement in uniformity can be further supported by the equidistribution quality measure which improves from 8.023253 for the initial mesh to 1.001855 for the final mesh.
* Example 5.5**.*
Let us now consider surfaces in . In this first example, we consider adaptive meshes for the torus defined by
[TABLE]
where , and . We take .
Fig. 7 shows the meshes for this example in two different views. Studying Fig. 7(a), the initial mesh, and Fig. 7(b), the final mesh with , we can see that the final mesh provides a more uniform distribution of the nodes. That is consistent with the use of the metric tensor whose goal is to make the mesh as uniform as possible in the Euclidean norm. This can also be confirmed from the equidistribution and alignment quality measures. The equidistribution measure for the initial mesh is 15.50150 and for the final mesh 1.332488. Similarly, the initial alignment quality measure is 30.63276 compared to that of the final mesh which is 1.920701.
For the curvature-based metric tensor, we see similar results to that of the Euclidean metric. That is, the final mesh for the curvature-based metric tensor, Fig. 7(c), looks identical to the final mesh for the Euclidean metric, Fig. 7(b). This is because the absolute value of the mean curvature of a torus is close to constant and hence, the elements do not concentrate in any particular region of the surface.
* Example 5.6**.*
The second three-dimensional example is the cylinder defined by
[TABLE]
where . For this example we take . Two boundary nodes were fixed, \mbox{\boldmathx}_{1}=(0,1,-2) and \mbox{\boldmathx}_{1}=(0,1,2), but the remaining boundary nodes were allowed to slide along the boundary. Although the cylinder has constant curvature like Example 5.5, this example shows the adaptation on a surface with a boundary.
Fig. 8 shows the adaptive meshes for the cylinder in two different views. For both and , the mesh becomes much more uniform and identical. This is consistent with the constant curvature of the cylinder hence the nodes do not concentrate in any specific region of the surface. The equidistribution quality measure improves from 19.07656 to 1.054857 and the alignment quality measure from 23.35403 to 1.192268. The fact that the final quality measures for both conditions are close to 1 indicates that the final meshes are close to satisfying conditions (7) and (9).
* Example 5.7**.*
Our next example is the sine surface in three dimensions defined by
[TABLE]
For this example we take and fix the boundary nodes.
Fig. 9 shows the adaptive meshes for this examples in two different views. It is clear in Fig. 9, when , the mesh becomes much more uniform with respect to the Euclidean metric from the initial mesh Fig. 9(a) to the final mesh Fig. 9(b). The top view of the surface, Fig. 9(d) and (e), further confirms this observation. It is also supported by the improvement of the quality measures from and to and .
When is curvature-based, we see a similar result to Example 5.3. That is, Fig. 9(c) and (f) show that the elements are more concentrated in those regions of the surface with larger curvature, i.e., the dip when and the hill when . The quality measures with respect to the metric tensor improve from to and to . The final quality measure for the equidistribution condition close to 1 hence indicating that the final mesh is close to satisfying (7). The final quality measure for the alignment condition is not as close to 1 as the equidistribution condition. Recall that in the meshing function (14) balances equidistribution and alignment and the choice has been used in the computation. Further computations show that increasing will improve the alignment quality but worsen the equidistribution quality, and vice versa. This suggests that a perfectly uniform mesh cannot be obtained by minimizing (14) for the curvature-based metric tensor for this example.
Finally, we would like to take a look at the changes of and along the mesh trajectory. As we recall from Section 4, is bounded from below and is decreasing. These can be seen numerically for in Fig. 10(a) and Fig. 10(b). Similar to what we saw in Example 5.3, Fig. 10(a) shows that is always decreasing and at around begins to converge. In Fig 10(b) we see an initial increase in the value and then it begins to converge to 4.64 at . This initial increase, as discussed above, is due to the nonuniformity of the initial mesh. That is, the initial mesh is very nonuniform and therefore can be very small whereas when the mesh is adapted, the mesh becomes more uniform and hence the values of become almost identical. This implies that the value of is likely to increase as the mesh adapts.
For the case with , Fig. 10(c) and (d) show similar findings. In 10(c) we see that is decreasing for all time and converging beginning at around . Fig. 10(d) shows initially increases then begins to converge to about 2.0. Furthermore, is bounded below by the initial value of . These numerical results for the curvature based metric tensor further support the theoretical predictions.
* Example 5.8**.*
Our final example explores the sphere and ellipsoid defined by an icosahedral initial mesh (see [26] for more details). We begin with the sphere
[TABLE]
For this example we take .
As we see from Fig. 11(a), the initial mesh is close to being uniform however, there is a very slight difference in the final mesh Fig. 11(b) when . Indeed, this slight adaptation can be seen in the quality measures which change from and for the initial mesh to and for the final mesh. The difference in the quality measures indicates that the initial icosahedral mesh is almost uniform and so the moving mesh method does not affect the mesh significantly.
We further this example to consider adaptive meshes for the ellipsoid defined by
[TABLE]
We move the mesh on the surface for both and .
First considering , Fig. 12 shows the meshes for this example in two different views. Studying Fig. 12(a) and Fig. 12(d), the initial mesh, and Fig. 12(b) and Fig. 12 (e), the final mesh with , we can see that the final mesh adapts to provide a higher concentration of elements in the middle region of the ellipsoid and fewer elements near the tips. The quality measures improve from and for the initial mesh to and for the final mesh. Although the initial mesh is close to uniform, the final mesh adapts in such a way to satisfy the equidistribution and alignment condition on the surface. However, this is not an accurate representation of the shape thus we consider a curvature-based metric tensor.
In our numerical experiments, when is used, we saw the mesh adapt in a similar way as with the Euclidean metric. This is because the curvature of the ellipsoid does not change significantly at the tips thus not many nodes move there. With this in mind, we altered the curvature-based metric tensor to concentrate more mesh elements at the tips of the ellipsoid by redefining as
[TABLE]
Fig. 12(c) and Fig. 12(f) show the final mesh using this altered metric tensor. As we can see, the mesh elements have concentrated at the tips of the ellipsoid better representing the shape of the surface. The equidistribution quality measure changes from initially to whereas the alignment quality measure from to . Similar results are seen with a finer mesh in Fig. 13 for both the Euclidean metric and altered curvature-based metric.
6 Conclusions and further comments
We have proposed a direct approach for surface mesh movement and adaptation that can be applied to a general surface with or without analytical expressions. We did so by first proving the relation (4) between the area of a surface element in a Riemannian metric and the Jacobian matrix of the affine mapping between the reference element and any simplicial surface element. From this we formulated the equidistribution and alignment conditions as given in (7) and (9), respectively. These two conditions enabled us to formulate a surface meshing function that is similar to a discrete version of Huang’s functional for bulk meshes [13]. The surface function satisifies the coercivity condition (37) for and .
We defined the surface MMPDE (32) as the gradient system of the meshing function, which utilizes surface normal vectors to inherently ensure that the mesh vertices remain on the surface during movement. Equations (26) and (27) give explicit, compact formulas for the mesh velocities making the time integration of the surface MMPDE (33) relatively easy to implement. Moreover, we showed that this surface MMPDE satisfies the energy decreasing property, which is one of the keys to proving Theorem 4.1. This theorem is an important theoretical result as it states that the surface mesh remains nonsingular for all time if it is so initially. Finally, we proved Theorem 4.2 that states the mesh has limiting meshes, all of which are nonsingular.
A point of emphasis is that the new method is developed directly on surface meshes thus, making no use of any information on surface parameterization. As mentioned, the MMPDE (32) only depends on surface normal vectors which can be computed even when the surface has a numerical representation. This allows the new method to be applied to general surfaces with or without explicit parameterization.
The numerical results presented in this work demonstrated that this new approach to surface mesh movement is successful. In all of the examples, the final mesh was seen to be much more uniform with respect to both cases of the metric tensor and which was supported by the mesh quality measures. Moreover, the theoretical properties were numerically verified in Ex. 5.3 and Ex 5.7 as we showed that is decreasing and is bounded below.
The future goal is to develop this algorithm for any surface with or without analytical expression. Even though we only presented examples which have analytical expressions, we should emphasize that the MMPDE (33) uses only the normal direction of the surface. Since these derivatives can be obtained numerically from the initial mesh or a background mesh, the method developed in this paper should work in principle for surfaces without explicit expressions. A practical difficulty is that the initial mesh or a background mesh typically does not represent the underlying surface accurately and approximate partial derivatives of obtained directly from the mesh may not lead to acceptable results. Our next step in the research is to investigate the use of spline approximations of surfaces for this purpose. Moreover, the monitor functions we used in the examples are limited to simple scalar metric tensors. It will be interesting to see how an anisotropic metric tensor such as one based on the shape map affects mesh movement and quality.
Acknowledgement. We would like to thank Dr. Lei Wang at the University of Wisconsin-Madison for providing us her code to generate initial icosahedral meshes for Example 5.8.
Appendix A Derivation of derivatives of the meshing function with respect to the physical coordinates
Recall from Section 3.2 that and our objective is to compute the derivatives
[TABLE]
Let be an entry of [\mbox{\boldmathx}_{1}^{K},\mbox{\boldmathx}_{2}^{K},\dots,\mbox{\boldmathx}_{d}^{K}]. Using the chain rule we have
[TABLE]
Denote
[TABLE]
[TABLE]
To begin, consider (50). When is an entry of [\mbox{\boldmathx}_{2}^{K},\dots,\mbox{\boldmathx}_{d}^{K}], recalling that E_{K}=[\mbox{\boldmathx}_{2}^{K}-\mbox{\boldmathx}_{1}^{K},\dots,\mbox{\boldmathx}_{d}^{K}-\mbox{\boldmathx}_{1}^{K}], we have
[TABLE]
which implies
[TABLE]
Moreover, for t=\left(\mbox{\boldmathx}_{1}^{K}\right)^{(1)} (the first component of \mbox{\boldmathx}_{1}^{K}), we have
[TABLE]
We can have similar expressions for \left(\mbox{\boldmathx}_{1}^{K}\right)^{(j)} for . This gives
[TABLE]
where \mbox{\boldmathe}^{T}=[1,\dots,1]\in\mathbb{R}^{1\times(d-1)}. For (51), we assume that \mathbb{M}=\mathbb{M}(\mbox{\boldmathx}) is a piecewise linear function defined on the current mesh, i.e., , where is the linear basis function associated with the vertex \mbox{\boldmathx}_{j}^{K} for all . Denote the th components of and \mbox{\boldmathx}_{K} by \mbox{\boldmathx}^{(i)} and \mbox{\boldmathx}_{K}^{(i)}, respectively. Then, for any entry of [\mbox{\boldmathx}_{1}^{K},\mbox{\boldmathx}_{2}^{K},\dots,\mbox{\boldmathx}_{d}^{K}], we have
[TABLE]
where we notice that \frac{\partial\phi_{j,K}}{\partial\mbox{\boldmathx}} and \frac{\partial\mbox{\boldmathx}_{K}}{\partial t} are a row and a column vector, respectively, and thus \frac{\partial\phi_{j,K}}{\partial\mbox{\boldmathx}}\frac{\partial\mbox{\boldmathx}_{K}}{\partial t} is a dot product. From this and the identity \mbox{\boldmathx}_{K}=(\mbox{\boldmathx}_{1}^{K}+\dots+\mbox{\boldmathx}_{d}^{K})/d, we get
[TABLE]
and
[TABLE]
Summarizing the above results, we have
[TABLE]
Notice that (53) can be rewritten as
[TABLE]
which gives (27).
Next, we establish the relations between
[TABLE]
First recall that , thus
[TABLE]
Let . Then we have
[TABLE]
Consider the first term of (56). Using the properties of matrix derivatives (24) we get
[TABLE]
Since , , and are all symmetric, it follows that
[TABLE]
The second term of (56) is then
[TABLE]
Therefore
[TABLE]
Combining this with (53) we obtain (26).
The identity (28) can be obtained similarly.
Finally, we derive the relations between
[TABLE]
First, note that the basis functions satisfy
[TABLE]
Eliminating the \mbox{\boldmathx}_{1}^{K} yields
[TABLE]
Then differentiating with respect to \mbox{\boldmathx}^{(k)} gives
[TABLE]
where \mbox{\boldmathe}_{k} is the unit vector in . Hence we have
[TABLE]
which gives
[TABLE]
and thus (29).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] B. Lévy and N. Bonneel. Variational anisotropic surface meshing with Voronoi parallel linear enumeration. The 21st International Meshing Roundtable, Springer-Verlag, 344-366, (2012).
- 2[2] J.-D. Boissonnat, F. Chazal, and M. Yvinec. Geometric and Topological Inference. Cambridge University Press, (2018).
- 3[3] P. Browne, C. J. Budd, M. Cullen, and H. Weller. Mesh adaptation on the sphere using optimal transport and the numerical solution of a Monge-Ampére type equation. J. Comput. Phys. , 308:102-123, (2016).
- 4[4] J. Cavalcante-Neto, M. Freitas, D. Siqueira, and C. Vidal. An adaptive parametric surface mesh generation method guided by curvature. Proceedings of the 22nd International Meshing Roundtable, Springer International Publishing. 425-443, (2014).
- 5[5] B. Crestel, R. D. Russell, and S. Ruuth. Moving mesh methods on parametric surfaces. Procedia Engineering , 124:148-160, (2015).
- 6[6] F. Dassi, S. Perotto, H. Si, and T. Streckenbach. A priori anisotropic mesh adaptation driven by a higher dimensional embedding. Computer-Aided Design , 85:111-122, (2017).
- 7[7] A. Demlow and G. Dziuk. An adaptive finite element method for the Laplace-Beltrami Operator on implicitly defined surfaces. SIAM J. Numer. Anal. , 45:421-442, (2007).
- 8[8] G. Dziuk. Finite elements for the Beltrami operator on arbitrary surfaces. Springer, New York. Lecture Notes in Math., Vol. 1357, (2006).