A surface moving mesh method based on equidistribution and alignment

Journal Pre-proof A surface moving mesh method based on equidistribution and alignment

Avary Kolasinski, Weizhang Huang

PII:

S0021-9991(19)30802-2

DOI:

https://doi.org/10.1016/j.jcp.2019.109097

Reference:

YJCPH 109097

To appear in:

Journal of Computational Physics

Received date:

27 January 2019

Revised date:

22 October 2019

Accepted date:

3 November 2019

Please cite this article as: A. Kolasinski, W. Huang, A surface moving mesh method based on equidistribution and alignment, J. Comput. Phys. (2019), 109097, doi: https://doi.org/10.1016/j.jcp.2019.109097.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

Highlights • Presents a surface moving mesh method for general surfaces with or without explicit pa- rameterization. • Develops the equidistribution and alignment conditions for surface meshes and establishes a meshing energy function based on the conditions. • Formulates a new surface moving mesh algorithm based on the gradient system of the en- ergy function, with the nodal mesh velocities being projected onto the underlying surface. • Proves that any mesh trajectory generated by the method remains nonsingular if it is so initially.

A surface moving mesh method based on equidistribution and alignment Avary Kolasinski∗†

Weizhang Huang‡

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 based on the current mesh. 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.

∗

Prepared by LLNL under Contract DE-AC52-07NA27344. LLNL-JRNL-782518. Department of Mathematics, the University of Kansas, Lawrence, KS 66045 ([email protected]). ‡ Department of Mathematics, the University of Kansas, Lawrence, KS 66045 ([email protected]). †

1

1. Introduction We are interested in methods that can directly move simplicial meshes on general surfaces with or without analytical expressions. Consider the four-petal rose (a curve or referred to as a ‘surface’ in two dimensions) shown in Fig. 1 as an example. An initial mesh is given as in Fig. 1(a). We want to develop a method that can use the information given by the initial mesh (mainly the location and connectivity of its vertices) but does not require any prior knowledge about the underlying curve including its analytical expression to move the mesh around for better uniformity or concentration. Meshes in Fig. 1(b,c) are obtained by such a method presented in this work for being more uniform/equidistant or more concentrated in places with larger curvature. The mesh movement is curve blind in the sense that it uses no information about the curve including its analytical expression although it is available for this example. 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, 27]. 0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

-0.2

-0.2

-0.2

-0.4

-0.4

-0.4

-0.6

-0.6

-0.6

-0.8

-0.8

-1

-0.5

0

0.5

(a) Initial Mesh

1

-0.8

-1

-0.5

0

0.5

1

(b) Final Mesh, equidistant

-1

-0.5

0

0.5

1

(c) Final Mesh, curvature-based

Figure 1: Meshes for four-petal rose.

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. [26] solve a Monge-Amp´ere 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. [24] 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. MacKenzie et al. [25] perform a moving mesh method to solve a forced curve shortening goemetric evolution equation. In this, the mesh is controlled via a tangential mesh velocity which is derived from the equidistribution principle with respect to an adaptivity criterion or a monitor function. Dassi et al. [6] generalize the higher embedding approach proposed in [1]. They modify the embedding map between the underlying

2

surface and R6 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 R6 . 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, 23] 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 (denoted by Ih , see (16) in Section 2.4), with the nodal mesh velocities being projected onto the underlying surface, i.e.,       ∂Ih T dxi ∂Ih T Pi =− − · ni ni , i = 1, ..., Nv (1) dt τ ∂xi ∂xi where Nv is the number of vertices, ni is the unit normal to the surface at the node xi , Pi is the value of a given positive scalar function evaluated at xi , and τ > 0 is a constant parameter. The choice of P and τ is discussed in Section 3. Note that (1) forms a system of ordinary differential equations (ODEs), with the right-hand side representing the nodal velocities for mesh movement. Moreover, it does not require any spatial discretization because it is already in spatially discrete form. 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. With these expressions, (1) can be rewritten as dxi Pi  K = v iK , i = 1, ..., Nv (2) dt τ K∈ωi

where K is an arbitrary element of the element patch ωi associated with xi and v K iK is the local mesh velocity given analytically in matrix form in (28), (29), and (36). ODE system (2) can be integrated in time using any scheme such as MATLAB’s ODE solvers ode45 or ode15s for meshes at discrete time instants. 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

3

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 based on the current mesh, the new method can apply to general surfaces with or without explicit parameterization. A selection of twoand 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 S be a bounded surface in Rd (d ≥ 2). Assume that we have a mesh Th = {K} on S and let N and Nv be the number of its elements and vertices, respectively. The elements are surface simplexes in Rd , i.e., they are (d − 1)-dimensional simplexes in a d-dimensional space. Notice that their area in d ˆ dimensions is equivalent to their volume in (d − 1) dimensions. Assume that the reference element K has been chosen to be a (d − 1)-dimensional equilateral and unitary simplex in a (d − 1)-dimensional ˆ and any element K ∈ Th let FK : K ˆ ⊂ Rd−1 → K ⊂ Rd be the affine mapping between space. For K  be its Jacobian matrix. Denote the vertices of K by xK ∈ Rd , j = 1, . . . , d and the them and FK j d−1 ˆ vertices of K by ξ j ∈ R , j = 1, . . . , d. Then xK j = FK (ξ j ), From this, we have or

j = 1, . . . , d.



K  xK j − x1 = FK ξ j − ξ 1 , 

j = 2, . . . , d

K K K  xK 2 − x1 , . . . , xd − x1 = FK [ξ 2 − ξ 1 , . . . , ξ d − ξ 1 ] ,

 =E E ˆ are the edge matrices for K and K, ˆ i.e., ˆ −1 , where EK and E which gives FK K  K K K ˆ = [ξ 2 − ξ 1 , . . . , ξ d − ξ 1 ] . E EK = xK 2 − x1 , . . . , xd − x1 ,

ˆ is a (d − 1) × (d − 1) square matrix and its inverse exists since K ˆ is not degenerate. Notice that E  ∈ Rd×(d−1) are not square. This makes the However, unlike the bulk mesh case, matrices EK , FK 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.

4

In the following we can see that the area of the physical element K ∈ Th can be determined using or EK .

 FK

Lemma 2.1. For any surface simplex K, there holds 1/2



|K|  T  FK , = det FK ˆ |K|

(3)

ˆ denote the area of the simplexes K and K, ˆ respectively, and det(·) denotes the where |K| and |K| determinant of a matrix.  =E E ˆ −1 , we have Proof. From FK K

1/2 1/2



 T  T ˆ −T EK ˆ −1 FK = det E EK E det FK T

1/2 ˆ −1 det EK = det(E) EK

1/2 T 1 = EK , det EK ˆ (d − 1)! |K|

ˆ = where we have used |K|

1 (d−1)!

ˆ Let the QR-decomposition of EK ∈ Rd×(d−1) be given by det(E).   RK E K = QK , 0

where QK ∈ Rd×d is a unitary matrix, RK ∈ R(d−1)×(d−1) is an upper triangular matrix, and 0 ∈ R1×(d−1) is a row vector of zeros. This decomposition indicates that K is formed by rotating the  convex hull with edges formed by the column vectors of 

 |K| = area(EK ) = area QK

RK 0

RK . We have 0



 = area

RK 0

 ,

where we have used the fact  that  rotation, QK , does not change the area. Since the convex hull formed by the column vectors of

RK lies on the x(1) – · · · – x(d−1) – plane, its area is equal to the (d − 1)0

dimensional volume of the convex hull formed by the column vectors of RK in (d − 1)-dimensions. Then, 1 1 T det(RK ) = det(RK RK )1/2 . |K| = volume(d−1) (RK ) = (d − 1)! (d − 1)! On the other hand, we have ⎛   T ⎞    R  T



T R R K K K T ⎠ = det RT 0 det EK EK = det ⎝ QK QK RK . = det RK K 0 0 0 Therefore,



1/2  T  det FK FK =

1 ˆ (d − 1)! |K|

1/2 T EK = det EK

5

1/2 |K| T RK = det RK . ˆ ˆ (d − 1)! |K| |K| 1

2.2. Area of surface elements in a Riemannian metric  or E . The formula We now formulate the area of a surface element in a Riemannian metric using FK K is needed later in the development of algorithms for mesh adaptation. First we consider a symmetric, uniformly positive definite metric tensor M(x) which satisfies

mI ≤ M(x) ≤ mI,

∀x ∈ S

(4)

where m and m are positive constants, I 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 M over K as  1 MK = M(x)dx. |K| K Recall that the distance in the Riemannian metric, MK , is defined as   T     1/2  1/2 1/2 T MK x MK x = MK x , x MK = x MK x =

(5)

where · denotes the standard Euclidean norm. This implies that the geometric properties of K in 1/2 the metric MK can be obtained from those of MK K in the Euclidean metric. Lemma 2.2. For any surface simplex K, there holds 1/2



|K|MK  T  M K FK , = det FK ˆ |K|

(6)

where |K|MK denotes the area of K in the metric MK . ˆ to M1/2 K is given by Proof. The Jacobian matrix of the affine mapping from K K

 1/2   ˆ −1 = M1/2 FK FM = MK EK E . K K ,K From the discussion following (5), we know that the area of K in the metric MK is equal to the area 1/2 of MK K in the Euclidean metric. Thus, from Lemma 2.1 we have 1/2 1/2



T  |MK K| |K|MK  = = det FMK ,K FMK ,K ˆ ˆ |K| |K| 1/2 



 1/2 1/2  T 1/2   T  = det MK FK MK FK = det FK M K FK .

The following lemma gives a lower bound for the area of K with respect to the metric MK in terms of the minimum altitude of K with respect to MK . Lemma 2.3. Let aK,MK denote the minimum altitude of K with respect to MK . Then, |K|MK ≥

1 (d − 1)

d−1 2

6

(d − 1)!

ad−1 K,MK .

(7)

 =E E ˆ −1 , we have Proof. From Lemma 2.2 and FK K 1/2



 T  ˆ det FK |K|MK = |K| M K FK

ˆ

1/2 T |K| M K EK det EK ˆ det(E)  T 1/2 1 1/2 1/2 det MK EK = . M K EK (d − 1)! =

1/2

Let the QR-decomposition of MK EK be denoted as 1/2

M K E K = QK



 RK , 0

where QK ∈ Rd×d is a unitary matrix, RK ∈ Rd−1×d−1 is an upper triangular matrix, and 0 is a (d − 1)-dimensional row vector of zeros. This gives  T 1/2 1 1/2 1/2 det MK EK M K EK |K|MK = (d − 1)!   1/2 1 R K T 0T ]QTK QK det [RK = (d − 1)! 0 =

1 T det(RK RK )1/2 (d − 1)!

=

d−1  1 si , (d − 1)! i=1

where si , i = 1, . . . , d − 1 denote the singular values of RK . Additionally, by [2, Lemma 5.12] we have that aR si ≥ √ K , d−1 where aRK denotes the minimum altitude of the simplex formed by the columns of RK . Combining these, we get 1 ad−1 |K|MK ≥ d−1 RK . 2 (d − 1) (d − 1)! Since QK is a rotational matrix, the minimum altitude of K with respect to the metric MK is the same as the minimum altitude of the convex hull formed by the columns of RK i.e., aK,MK = aRK . Thus, we have obtained (7). 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 (d = 2), K is a line segement and both |K|MK and aK,MK represent the length of K in the metric MK and are equal. In this case, the inequality (7) reduces to |K|MK ≥ aK,MK ,

7

which is very sharp. For d = 3, (7) becomes |K|MK ≥

1 2 a , 4 K,MK

which is not as sharp as in two dimensions. Indeed, when K is equilateral with respect to MK , we have [11] 1 |K|MK = √ a2K,MK . 3

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 |K|MK =

σh , N

∀K ∈ Th

(8)

where, as before, | · |MK denotes the area of the surface with respect to the metric MK and  ˆ = 1, we have σh = K∈Th |K|MK . Using Lemma 2.2 and recalling |K| |K|MK = det



 FK

T

 M K FK

1/2

,

σh =



det



 FK

T

 M K FK

1/2

.

K∈Th

Thus, the equidistribution condition (8) becomes det



 FK

T

 M K FK

1/2

=

σh , N

∀K ∈ Th .

(9)

The alignment condition, on the other hand, requires that all of the elements K ∈ Th be similar ˆ Notice that any element K is similar to K ˆ if and only if FK : K ˆ →K to the reference element K.  is composed by dilation, rotation, and translation, or equivalently, FK is composed by dilation and  can therefore be expressed as rotation. Mathematically, FK   I  = αU V T, (10) FK 0 where α is a constant representing dilation and U ∈ Rd×d and V ∈ R(d−1)×(d−1) are orthogonal matrices representing rotation. Thus,   I 1/2 1/2  V T. MK FK = α MK U 0 It can be verified that the above equation is equivalent to 



 1



1 d−1  T   T  = det FK tr FK M K FK M K FK , d−1

8

∀K ∈ Th

(11)

which is referred to as the alignment condition. Here, tr(·) 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 (11) and note that an equivalent condition is 1   −1  −1  d−1  T  T 1   tr FK MK FK = det F K M K FK . d−1

Notice that the left- and right-hand sides are the arithmetic mean and the geometric mean of the  T

 −1 , respectively. The inequality of arithmetic and geometric eigenvalues of the matrix (FK ) M K FK means gives 1   −1  −1  d−1  T  T 1   tr FK MK FK ≥ det FK MK FK , (12) d−1 with equality if and only if all of the eigenvalues are equal. From (12), for any general mesh which does not necessarily satisfy (11), we have  tr

−1  T  F K M K FK

d−1

≥ (d − 1)d−1 det



 FK

T

 M K FK

−1

,

and therefore  tr

−1  T  F K M K FK

 p(d−1) 2

− (d − 1)

p(d−1) 2

det



 FK

T

 M K FK

− p 2

≥ 0,

where p > 0 is a dimensionless parameter which has been added to agree with the equidistribution energy function below. Then, we define the alignment energy function as Iali =



ˆ det |K|



 T FK

K∈Th

− (d − 1)

p(d−1) 2

 M K FK





1 2



tr

ˆ det |K|



 T FK

 FK

T

 M K FK

 M K FK

−1  p(d−1) 2

 1−p 2

,

(13)

K∈Th

whose minimization will result in a mesh that closely satisfies the alignment condition (11). One may 1

2 T   ˆ = |K|MK has been added as a weight. notice that |K| det (FK ) MK FK Similarly, we consider the equidistribution condition (9). From H¨older’s inequality, for any p > 1 we have ⎞1 ⎛ p  |K|M  |K|M 1 1 ⎟ ⎜ K K · · (14)





1/2 ≤ ⎝ p/2 ⎠ , σh σh T    T M F K∈Th K∈Th det FK MK FK det FK K K with equality if and only if det



−1/2  T  FK M K FK = constant,

9

∀ K ∈ Th .

That is, minimizing the difference between the left-hand side and the right-hand side of (14) tends

−1/2  )T M F  to make det (FK constant for all K ∈ Th . Noticing that the left-hand side of (14) is K K simply N/σh , we can rewrite this inequality as  p



 1−p  N 2  T  ˆ det FK · σh ≤ |K| M K FK . σh K∈Th  det(M(x))1/2 dx and hence it only weakly depends on the We can consider σh constant since σh ≈ S

mesh. Therefore, we define the equidistribution energy function as



 1−p p(d−1)  2  T  ˆ det FK |K| M K FK , Ieq = (d − 1) 2

(15)

K∈Th

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 (13) and (15), that is, define Ih = θIali + (1 − θ)Ieq for θ ∈ [0, 1]. This leads to  1 

−1  p(d−1)



 2 2  T   T  ˆ Ih = θ |K| det FK MK FK tr F K M K FK K∈Th

+ (1 − 2θ)(d − 1)

p(d−1) 2



ˆ det |K|



 FK

T

 M K FK

 1−p 2

,

(16)

K∈Th

where p > 1 and θ ∈ [0, 1] are dimensionless parameters, with the latter balancing the equidistribution and alignment conditions for which full alignment is achieved when θ = 1 and full equidistribution is achieved when θ = 0. We can write (16) as



1  2  T  ˆ det FK ˜K , Ih = |K| M K FK (17) G K∈Th

where ˜ K = θ tr G



 T FK



 M K FK

−1  p(d−1) 2

+ (1 − 2θ)(d − 1)

p(d−1) 2

det



 FK

T

 M K FK

− p 2

.

(18)

˜ K is coercive; see the definition of coercivity in Section We remark that for 0 < θ ≤ 12 and p > 1, G 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 (16) is very similar to a Riemann sum of the meshing function developed in [13] for bulkmeshes based on equidistribution and alignment. One of the main  )T M F  differences is that (FK K K cannot be simplified in (16) since it is not a square matrix as it is in the bulk mesh case. Additionally, the constant terms and exponents that contain d are (d − 1) in (16) instead of d in the bulk mesh case. The functional of [13] has been proven to work well for a variety of problems [20].

10

3. Moving mesh equations for surface meshes In this section we describe the surface MMPDE method used to find the minimizer of (16). In principle, we can directly minimize it. However, this minimization problem can be extremely difficult to solve since (16) 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 (16), we consider meshing functions in a general form (17), i.e.,



1   2  T  ˆ ˜K ≡ |K| det FK MK FK G G(JK , rK ), Ih = K∈Th

˜ K is a given smooth function of where G



−1  T  M K FK , JK = FK ˜ K , rK ), and ˜ K = G(J that is, G

(19)

K∈Th

rK = det



 FK

T

 M K FK

−1

,

1

ˆ −2 G ˜K . G(JK , rK ) = |K|r K

(20)

˜ K can be chosen differently. Moreover, both JK and rK depend Indeed, a special example is (18) but G on the coordinates of the vertices of the physical element K and hence G is a function of them, i.e., G (JK , rK ) can be expressed as

K (21) G (JK , rK ) = GK xK 1 , . . . , xd , d where xK i ∈ R for i = 1, . . . , d are the coordinations of the vertices of K. As a consequence, the sum in (19) is a function of the coordinates of all vertices of the physical mesh Th , i.e., 

K (22) GK x K Ih (x1 , . . . , xNv ) = 1 , . . . , xd , K∈Th

where xi ∈ Rd for i = 1, . . . , Nv 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 Ih with respect to the physical K coordinates x1 , . . . , xNv which requires elementwise derivatives of GK with respect to xK 1 , . . . , xd . That is,  ∂GK  ∂GK ∂Ih = = , i = 1, . . . , Nv (23) ∂xi ∂xi ∂xK iK K∈T K∈ω h

i

where iK denotes the local index of vertex xi in K and ωi is the element patch associated with xi . In order to calculate the necessary derivatives, we recall some definitions and properties of scalar-bymatrix differentiation (cf. [16] for details). Let f = f (A) be a scalar function of a matrix A ∈ Rm×n . Then the scalar-by-matrix derivative of f with respect to A is defined as ⎡ ∂f ⎤ · · · ∂A∂fm1   ∂A11 ∂f ∂f ∂f ⎢ .. ⎥ . . =⎣ . or = . (24) ⎦ . ∂A ∂A i,j ∂Aj,i ∂f ∂f ∂A1n · · · ∂Amn n×m

11

The chain rule of differentiation with respect to t is    ∂f ∂Aj,i   ∂f  ∂Aj,i ∂f ∂A ∂f = = = tr . ∂t ∂Aj,i ∂t ∂A i,j ∂t ∂A ∂t ij

(25)

ij

With this and when A is a square matrix, the following properties have been proven in [16],   ∂ det(A) ∂A−1 ∂tr (A) −1 ∂A −1 −1 ∂A = I, = −A A , = det(A) tr A . ∂A ∂t ∂t ∂t ∂t

(26)

∂G Using the above, we can find the expressions for ∂G ∂J and ∂r which are needed to compute (23). For the function (16), the first derivatives of G are given by ⎧ p(d−1)−2 ∂G θp(d − 1) ˆ − 1 ⎪ ⎪ = |K|r 2 tr(J) 2 I, ⎪ ⎪ ⎨ ∂J 2 (27) ⎪ ⎪ ⎪ p(d−1) p(d−1) ⎪ ˆ − 32 tr(J) 2 + p − 1 (1 − 2θ)(d − 1) 2 |K|r ˆ p−3 ⎩ ∂G = − θ |K|r 2 . ∂r 2 2

3.2. Derivatives of the meshing function with respect to the physical coordinates From (23), we can see that we will need ∂GK /∂xK iK to compute ∂Ih /∂xi . The former can be obtained once we know the derivatives of GK with respect to the coordinates of all vertices of K, i.e., ⎡ ⎤ ∂GK ⎢ ∂xK ⎥ ⎢ 1 ⎥ ∂GK ⎢ .. ⎥ = ⎢ . ⎥. K K K ⎥ ∂[x1 , x2 , . . . , xd ] ⎢ ⎣ ∂GK ⎦ ∂xK d The derivation of these derivatives is given in Appendix A. They are given as ⎤ ⎡ ∂GK ⎢ ∂xK ⎥ ⎢ 2 ⎥ T

−1 T ∂GK ⎢ .. ⎥ T T ˆ K ˆ E(E M K EK MK EK )−1 EK MK E ⎢ . ⎥ = − 2 EK ⎥ ⎢ ∂J ⎣ ∂GK ⎦ ∂xK d ˆ 2

−1 T ∂GK T det(E) T

E K M K EK −2 EK M K det EK MK EK ∂r ⎤ ⎡ ∂φj,K  ∂x  d ⎢ . ⎥ ∂GK 1 . ⎥ + tr Mj,K ⎢ ⎣ . ⎦, d ∂M j=1

K

∂φj,K ∂x

  d d  ∂φj,K ∂GK ∂GK ∂GK  , =− + tr Mj,K K K ∂MK ∂x ∂x1 ∂x j j=2 j=1

12

(28)

(29)

where ∂GK /∂J and ∂GK /∂r are given in (27), φj,K is the linear basis function associated with xK j , K Mj,K = M(xj ), and ∂GK T T T ˆ K ˆ T ∂GK E(E = −EK (EK MK EK )−1 E MK EK )−1 EK ∂MK ∂J ˆ 2 T

−1 T ∂GK det(E)

T E K EK − M K EK EK , det EK MK EK ∂r ⎤ ⎡ ∂φ2,K d  ⎢ ∂.x ⎥ ∂φ1,K ∂φj,K T −1 T ⎢ .. ⎥ = (EK = − . E ) E , K K ⎦ ⎣ ∂x ∂x ∂φd,K ∂x

(30)

(31)

j=2

Having computed ∂GK /∂xK j (j = 1, ..., d) for all elements using (28) and (29), we can obtain ∂Ih /∂xi from (23).

3.3. Surface moving mesh equations As mentioned above, we employ a surface MMPDE method to minimize the meshing function (16) or a more general form (19). 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 ∂Ih , i = 1, · · · , Nv . ∂xi Let Φ(x) = 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 Φ(xi ) = 0, i = 1, ..., Nv or at least dxi · ∇Φ(xi ) = 0, dt

i = 1, ..., Nv

(32)

xi is the nodal mesh velocity. Following the MMPDE approach, we would define the mesh where ddt equation as the gradient system of Ih , i.e.,   dxi Pi ∂Ih T =− , i = 1, ..., Nv (33) dt τ ∂xi where Pi is a positive scalar function used to make the equation have desired invariance properties and τ > 0 is a constant parameter used for adjusting the time scale of mesh movement. Obviously, this does not satisfy (32). Here we propose to project the velocities in (33) onto the surface and define the surface moving mesh equation as       dxi ∂Ih T ∂Ih T Pi =− − · ni ni , i = 1, ..., Nv (34) dt τ ∂xi ∂xi where ni = ∇Φ(xi )/ ∇Φ(xi ) is the unit normal to the surface at xi and the difference inside the square bracket is the projection of the vector ∂Ih /∂xi onto the tangential plane of the surface at xi .

13

Notice that this surface MMPDE inherently ensures that (32) be satisfied or, in words, the nodes stay on the surface during the mesh movement. Moreover, it is important to note that (34) only utilizes the unit normal vectors of the surface whose computation does not require explicit parameterization or analytical expression of the surface. In Section 5, we will present numerical results obtained with the normal vectors computed numerically based on the current mesh and their comparison with those obtained with the normal vectors computed using analytical formulas. Using (23) we can rewrite the above equation in a compact form as dxi Pi  K = v iK , i = 1, ..., Nv (35) dt τ K∈ωi

d K where v K iK ∈ R is the local mesh velocities contributed by K to xiK and has the expressions ⎞ T ⎛ T  ∂GK ∂GK +⎝ · n iK ⎠ niK , vK iK = − K ∂xiK ∂xK iK

(36)

and ∂GK /∂xK iK is given in (28) and (29). The surface MMPDE (35) must be modified properly for boundary vertices when S has a boundary. For fixed boundary vertices, the corresponding equation is replaced by dxi = 0. dt The velocities for other boundary vertices should be modified such that they slide on the boundary. It should be noted that (35) forms a system of ODEs. Its right-hand side represents the nodal velocities for mesh movement and is given analytically in matrix form in (28), (29), and (36), with the former making the implementation of the method relative easy and robust. Moreover, the system is already in spatially discrete form so it does not need any spatial discretization. With proper modification of the boundary vertices, the system (35) can be integrated in time. To do so, one first ˆ for the reference element. One starts by calculating the edge matrices EK for all elements and E can then readily calculate (27) which is needed for (28) and (29). Then one can integrate (35) 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 t → ∞ for the MMPDE (35).

4.1. Equivalent measure of minimum height

−1  )T M F  We begin the theoretical analysis by establishing the relation between (FK and the K K minimum altitude of K with respect to MK .

14

Lemma 4.1. There holds a ˆ2 a2K,MK

 −1   (d − 1)2 a   T ˆ2  ≤ ≤ F M F , K K K   2 aK,MK

(37)

ˆ and aK,M is the minimum altitude of K with respect to the metric MK . where a ˆ is the altitude of K K Proof. First of all, we have    −1    T  −1   −T T    T −1 T  ˆ EK E K E ˆ EK ˆ −1 ˆ   F K FK = E = E E E  . K     Now, consider the QR decomposition of EK



E K = QK

 RK , 0

where QK ∈ Rd×d is a unitary matrix, RK ∈ R(d−1)×(d−1) is an upper triangular matrix, and 0 is a (d − 1)-dimensional row vector of zeros. With this we have     −1        R −1 T  ˆ T K T T T T ˆ [RK 0 ]QK QK ˆ  = E ˆ  E E E EK EK   0      ˆ −1 −T ˆ T  = ER K RK E   −1 −T    −1 −1 ˆ ˆ . = RK E  RK E  By [16, Lemma 4.1] we have  −1 −T   (d − 1)2 a  ˆ2 a ˆ2 −1 −1 ˆ ˆ ≤  R ≤ R , E E K K   a2RK a2RK where aRK is the minimum altitude of the simplex formed by the columns of RK . Since QK is a rotation matrix, aRK is the same as aK , the minimum altitude of K with respect to the Euclidean metric. Combining the above results, we get     T  −1  (d − 1)2 a a ˆ2 ˆ2   ≤ ≤ F F . K K   a2K a2K The inequality (37) follows from this and the observation that the geometric properties of K with 1/2 respect to the metric MK are the same as those of MK K with respect to the Euclidean metric. ˆ is chosen to satisfy |K| ˆ = O(1) then Lemma 4.1 indicates that if K  −1     T   ∼ a−2 .  FK MK FK K,MK  

15

(38)

4.2. Mesh nonsingularity We now consider the MMPDE (35). 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 ˆ is taken to satisfy |K| ˆ = 1 instead of being unitary that for theoretical purposes, we assume that K N as we have been considering thus far. This change does not affect the actual computation. However, ˆ = 1 will likely lead to F  = O(1) and since typically we expect |K| = O(1/N ), the assumption |K| K N ˆ =1 thus Ih (Th (0)) (the value of Ih on the initial mesh Th (0)) stays O(1). On the other hand, if |K|  (unitary), we have FK = O(1/N ) and Ih (Th (0)) will depend strongly on N . In the following analysis, the mesh at time t is denoted by Th (t) = (x1 (t), . . . , xNv (t)). Theorem 4.1. Assume that the meshing function in the form (19) satisfies the coercivity condition   −1 q  T  ˜ (J, det (J) , x) ≥ α tr G FK M K FK − β, ∀x ∈ S (39) ˆ is equilateral and where q > (d − 1)/2, α > 0, and β ≥ 0 are constants. We also assume that K 1 ˆ = . Then if the elements of the mesh trajectory of the MMPDE (35) have positive areas |K| N initially, they will have positive areas for all time. Moreover, their minimum altitudes in the metric MK and their areas in the Euclidean metric are bounded below by aK,MK ≥ C1 |K| ≥ C2 where

⎛ C1 = ⎝

αd

q(d−2) d−1



¯ d/2 |S| Ih (Th (0)) + β m



¯ d/2 |S| Ih (Th (0)) + β m

(d − 1)!

(d − 1)

2q−d+1 d−1

d−1+2q 2

⎞

−

1 2q−d+1

−

d−1 2q−d+1

N

2q − (d−1)(2q−d+1)

N

2q − 2q−d+1

1 2q−d+1

⎠

,

C2 =

m

,

− d2

(40)

,

C1d−1 (d − 1)

d−1 2

(d − 1)!

Proof. From (34) we have       Pi ∂Ih  ∂Ih T dIh  ∂Ih dxi ∂Ih T = =− − · ni n i dt ∂xi dt τ ∂xi ∂xi ∂xi i i ⎡ 2 ⎤ 2    T  Pi  ∂Ih  ∂Ih  − ⎣ =− · ni ⎦ τ  ∂xi  ∂xi i

≤ 0.

16

(41)

.

(42)

This implies Ih (Th (t)) ≤ Ih (Th (0)) for all t. From coercivity (39) and Lemma 4.1, we get −1 q



1/2  

   T  T  ˆ Ih (Th (t)) ≥ α tr F K M K FK − βm ¯ d/2 |S| |K| det FK MK FK K∈Th

≥α



ˆ det |K|

K∈Th

≥α



ˆ det |K|





 T FK  FK

T

 M K FK  M K FK

1/2  −1    T q   F K M K FK  − βm ¯ d/2 |S|   1/2

a ˆ2q

K∈Th

1/2 ˆ det (F  )T MK F  By Lemma 12, |K| = |K|MK ≥ K K Ih (Th (t)) + β m ¯ d/2 |S| ≥ and therefore

(d − 1)

d−1 2

1

(d−1)

d−1 2 (d−1)!

(d − 1)

d−1 2

ad−1 K,MK , thus



αˆ a2q (d − 1)!

1

2q−d+1 K∈Th aK,MK

αˆ a2q

a2q−d+1 K,MK ≥

− βm ¯ d/2 |S|.

a2q K,MK

(d − 1)!

Ih (Th (0)) + β m ¯ d/2 |S|

ˆ is equilateral and |K| ˆ = Moreover, from the assumption that K √

1 N

−1

,

(43)

.

(44)

it follows that

1

1 d (d − 1)! d−1 a ˆ= √ N − d−1 . 1 d − 1 d 2(d−1)

(45)

Combining (44) and (45) we get ⎛ aK,MK ≥ ⎝

αd

q(d−2) d−1

(d − 1)!

(d − 1)

2q−d+1 d−1

d−1+2q 2

⎞

1 2q−d+1

⎠



Ih (Th (0)) + β m ¯ d/2 |S|

−

1 2q−d+1

N

2q − (d−1)(2q−d+1)

,

(46)

which gives (40). Furthermore, we have ad−1 K,MK (d − 1)

d−1 2

(d − 1)!

ˆ det ≤ |K|MK = |K| ≤m

d/2

ˆ det |K|





 FK

 FK

T

T

 FK

 M K FK

1/2

1/2

=m

d/2

|K|.

Then (41) follows from the above inequality and (40). Finally, from (28) and (29) it is not difficult to see that the magnitude of the mesh velocities is bounded from above when |K| is bounded from below. As a consequence, the mesh vertices will move continuously with time and |K| cannot jump over the bound to become negative. Hence, |K| 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 (35) by design while the

17

latter is an assumption for the meshing function. We emphasize that the result holds for any function satisfying the coercivity condition (39). On the other hand, the condition (39) is satisfied by the meshing function (16) for 0 < θ ≤ 12 and p > 1 (with q = (d − 1)p/2 and β = 0 in Theorem 4.1). It is interesting to point out that the role of the parameter p can be explained from (40). Indeed, for this case the inequality (40) becomes 1 −  p 1 (d−1)(p−1) − aK,MK ≥ C1 Ih (Th (0)) + β m ¯ d/2 |S| N (d−1)(p−1) → C1 N − d−1 , p → ∞. (47) 1

Since N − d−1 represents the average diameter of the elements, the above inequality implies that the mesh becomes more uniform as p gets larger. We take p = 1.5, 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 {Th (t), t > 0} of MMPDE (35) has the following properties. 1. Ih (Th (t)) has a limit as t → ∞, i.e., lim Ih (Th (t)) = L.

(48)

t→∞

2. The mesh trajectory has limiting meshes, all of which are nonsingular and satisfy (40) and (41). 3. The limiting meshes are critical points of Ih , i.e., they satisfy ∂Ih = 0, ∂xi

i = 1, . . . , Nv .

(49)

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 Ih (Th (t)) and the compactness of S. 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 (35) 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 (35). 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 (Qeq ) and alignment (Qali ) mesh quality measures which are defined as 

  T

1  )T M F  −1  tr (F 2 K K K det (FK ) MK FK Qeq = max and Qali = max

− 1 . (50)  K∈Th K∈Th σh /N d−1 (d − 1) det (F )T M F  K

18

K

K

These measures are indications of how closely the mesh satisfies the equidistribution condition (9) and the alignment condition (11), respectively. The closer these quality measures are to 1, the closer they are to a uniform mesh with respect to the metric MK . It should be noted that the alignment condition does not apply to the two-dimensional case where a “surface” is actually a curve. Mathematically,  )T M F  is a number and hence (11) is always satisfied. when d = 2, (FK K K For all computations we use p = 3/2 and θ = 1/3 in the meshing function (16). 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 τ = 0.01, dt = 0.01, and Pi = det (M(xi ))

p(d−1)−d 2

.

The latter is to ensure that the MMPDE (35) be invariant under scaling transformations of M. For all of the results, we run to a final time of 1.0. We choose two main forms of MK . The first is MK = I, 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 MK = (kK + ) I, where kK is the mean curvature and is machine epsilon. The mean curvature is defined (e.g., see [12]) for a curve Φ(x, y) = 0 in R2 as ) ) ) ) 2 2 ) Φxx Φy − 2Φxy Φx Φy + Φx Φyy ) ) ) k=)

3 ) 2 2 2 ) ) Φ x + Φy and for a surface Φ(x, y, z) = 0 in R3 as

) ) )D + D + D − D ) ) 1 2 3 4) k=) ), ) 2 Φ2 + Φ2 + Φ2 3/2 ) x

y

z

where D1 = Φx (Φx Φxx + Φy Φxy + Φz Φxz ) , D2 = Φy (Φx Φxy + Φy Φyy + Φz Φyz ) , D3 = Φz (Φx Φxz + Φy Φyz + Φz Φzz ) ,

D4 = Φ2x + Φ2y + Φ2z (Φxx + Φyy + Φzz ) . To further illustrate the surface moving mesh method and show how the mesh adapts accordingly to the metric tensor provided, we present a two-dimensional and a three-dimensional example with a gradient-based metric tensor (see Examples 5.3 and 5.8). We would like to explore more metric tensors in future work but will focus on these for this paper. In our computation we have used the analytical formula and numerical approximation for normal vectors. In the later case, normal vectors are computed based on the current mesh. They are obtained first for elements, and then for vertices by simply averaging those on neighboring elements. This approximation method is simple and straightforward however, our limited experience shows that it works well with all of the examples we have tested. In the future, this could be made more accurate using spline functions. We should point out, though, that the surface moving mesh method based on numerically approximated normal vectors can fail when the initial mesh is too coarse (to give a

19

reasonable representation of the underlying surface). To explain this, we consider the cardioid defined by

2 Φ(x, y) = x2 + y 2 + 4x(x2 + y 2 ) − 4y 2 . (51) We take MK = I and fix the node x1 = (0, 0). Fig. 2 compares the final meshes using analytical normal vectors to those using numerically approximated ones. We see that the mesh adapts from a very nonuniform initial mesh (Fig. 2(a) and (d)) to a uniform final mesh (Fig. 2(b), (c), (e), and (f)), which is consistent with the fact that the metric tensor corresponds to the Euclidean metric. However, for N = 10, the nodes move off the curve when numerically approximated normal vectors are used (Fig. 2(c)) whereas when analytical normal vectors are used, the nodes remain on the curve (Fig. 2(b)). For N = 40, the nodes remain on the curve for both cases, as shown in Fig. 2(e) and (f). This example shows that the initial mesh must be fine enough in order for the nodes to remain on the surface during movement when numerically approximated normal vectors are used. When the initial mesh is sufficiently fine, however, this is not a problem and the final meshes when using approximate or analytical normal vectors are almost identical. The same has been observed for all of the examples we have tested. If, however, in application the initial mesh is too coarse, one can use a surface fitting (e.g. using spline functions) to obtain a finer mesh and proceed with the surface MMPDE method proposed. An alternative approach is to employ a more stabilized version of the moving mesh PDE. More specifically, the MMPDE (34) is formulated directly using Lagrange multipliers imposing the condition ∇Φ ·

dxi = 0. dt

(52)

However, we can impose a stronger condition to obtain a stabilized MMPDE. That is, we impose Φ(xi (t + τ δ)) = 0,

(53)

where δ ≥ 0 is a constant and τ is the same parameter used in (33). In this, we are not requiring that the nodes be directly projected onto the surface at time t, i.e., Φ(xi (t)) = 0, but instead be projected onto the surface at some time later t + τ δ. This is called a delayed projection which is a stronger imposition than (52) and results in the following moving mesh PDE       ∂Ih T ∂Ih T dxi Pi Φ =− ni . − · ni n i − (54) dt τ ∂xi ∂xi τ δ ∇Φ The role of the last term can be seen through the different values of δ. That is, from (53) we can see that if δ = 0 then Φ(xi (t)) = 0 and hence we project the nodes directly onto the surface. If δ = ∞ then (54) gives (34), i.e., we project the nodes onto the tangential space. For all values 0 < δ < ∞, (54) uses a delayed projection of the nodes onto the surface. It should be noted that we use τ δ as the delay parameter so that the dimension of τ in all of the terms in (54) agrees. Note that (54) requires we know Φ or at least its approximation. For more details including examples displaying the effectiveness of the delayed projection method, see [22]. It is remarked that we use analytical normal vectors for some examples and numerical normal vectors for the others to demonstrate the ability of the surface moving mesh method to work with both. Nevertheless, we observe that both types of normal vectors lead to almost the same results for all of

20

2.5

2.5

2

2

2

1.5

1.5

1.5

2.5

1

1

1

0.5

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-1

-1

-1

-1.5

-1.5

-1.5

-2

-2

-2

-2.5

-2.5

-2.5

-5

-4

-3

-2

-1

0

-5

1

(a) Initial Mesh.

-4

-3

-2

-1

0

-5

1

(b) Final Mesh, analytical n.

2

2

2

1.5

1.5

1.5

1

1

1

0.5

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-1

-1

-1

-1.5

-1.5

-1.5

-2

-2

-2

-2.5

-2.5

-2.5

-3

-2

-1

0

(d) Initial Mesh.

1

-2

-1

0

1

2.5

2.5

-4

-3

(c) Final Mesh, approximate n.

2.5

-5

-4

-5

-4

-3

-2

-1

0

1

(e) Final Mesh, analytical n.

-5

-4

-3

-2

-1

0

1

(f) Final Mesh, approximate n.

Figure 2: Meshes of N = 10 and N = 40 for the cardioid using analytical and approximate normal vector with the Euclidean metric tensor.

21

the examples, as shown in Fig. 2.

Example 5.1. For the first example, we generate adaptive meshes for the unit circle in two dimensions, Φ(x, y) = x2 + y 2 − 1. We take N = 80 and fix the node x1 = (1, 0). Analytical normal vectors are used in this example. Fig. 3 shows the meshes for this example. Studying the figures we see that the initial mesh Fig. 3(a) is very nonuniform but the final meshes Fig. 3(b) and (c) have adapted to be equidistant along the curve. Moreover, the final meshes for both MK = I (Fig. 3 (b)) and MK = (kK + )I (Fig. 3(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 Qeq improves from 7.509 to 1.000 for both cases of MK . The fact that Qeq ≈ 1 indicates that the mesh is close to satisfying the equidistribution condition (9) and hence the mesh is almost uniform with respect to the metric tensor MK . 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. 0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

-0.2

-0.2

-0.2

-0.4

-0.4

-0.4

-0.6

-0.6

-0.6

-0.8

-0.8

-1

-0.5

0

0.5

(a) Initial Mesh

1

-0.8

-1

-0.5

0

0.5

1

(b) Final Mesh, MK = I

-1

-0.5

0

0.5

1

(c) Final Mesh, MK = (kK + )I

Figure 3: Example 5.1. Meshes of N = 80 are obtained for Φ(x, y) = x2 + y 2 − 1.

Example 5.2. The second two-dimensional example is the ellipse defined by x2 + y 2 − 1. 64 In this example we take N = 60 and fix the node x1 = (8, 0). Numerically approximated normal vectors are used in this example. The initial nodes (Fig. 4(a)) are randomly distributed through the curve. However, for MK = I, the final mesh (Fig. 4(b)) is equidistant along the ellipse providing a much more uniform mesh. This can also be seen in Qeq which improves from 5.497 initially to 1.026 in the final mesh. Φ(x, y) =

22

Now considering the curvature-based metric tensor (Fig. 4(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 x = −8, 8 and almost 0 for x ∈ (−2, 2). From Fig. 4(c) we can see that the adaptation with MK = (kK + )I provides a mesh that represents the shape of the curve much better than other two meshes. The improvement of Qeq from 5.126 to 1.015 indicates that the final mesh is almost uniform with respect to the curvature-based metric tensor. 6

6

6

4

4

4

2

2

2

0

0

0

-2

-2

-2

-4

-4

-4

-6

-6

-6 -6

-4

-2

0

2

4

(a) Initial Mesh

6

8

-6

-4

-2

0

2

4

6

(b) Final Mesh, MK = I

8

-6

-4

-2

0

2

4

6

8

(c) Final Mesh, MK = (kK + )I

Figure 4: Example 5.2. Meshes of N = 60 are obtained for Φ(x, y) =

x2 + y 2 − 1. 64

Example 5.3. For the next two-dimensional example, we generate adaptive meshes for the sine curve defined by Φ(x, y) = 4 sin(x) − y. In this example we take N = 60 and fix the end nodes x1 = (0, 0) and x61 = (2π, 0). We use analytical normal vectors in this example. Fig. 5 shows the meshes for this example. From Fig. 5(a) and (b) we see that for MK = I, the mesh becomes much more uniform. This is consistent with the fact that for M = I, 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 Qeq measure improves from 4.183 to 1.002 indicating that the final mesh satisfies the equidistribution condition (9) closely. To show mesh compression (in Euclidean metric) more clearly, we also consider the ratio of the area of the largest element to the smallest one, i.e., ν=

maxK∈Th |K| . minK∈Th |K|

(55)

Obviously, the larger ν is, the more compressive and nonuniform the mesh elements are. On the other hand, ν = 1 implies a perfectly uniform mesh. In this example, ν changes from 25.1 to 1.05 indicating that the initial mesh is nonuniform whereas the final mesh is close to being uniform.

23

Now studying Fig. 5(c) where MK = (kK + )I is used, we see that there is a high concentration of mesh elements in regions with large curvature, i.e., the hill at y = 4 and cup at y = −4, which is consistent with the use of the curvature-based metric tensor. Moreover, the equidistribution measure Qeq improves from 6.254 to 1.007. This indicates that although the mesh may seem nonuniform in the Euclidean metric, it is almost uniform in the metric MK . Considering the area ratio as defined in (55), we see an change from 25.1 for the initial mesh to 243.2 for the final mesh. This implies that the final mesh is extremely compressive in the Euclidean metric. 3

3

3

2

2

2

1

1

1

0

0

0

-1

-1

-1

-2

-2

-2

-3

-3

-3

-1

0

1

2

3

4

5

(a) Initial Mesh

6

7

8

-1

0

1

2

3

4

5

6

7

(b) Final Mesh, MK = I

8

-1

0

1

2

3

4

5

6

7

8

(c) Final Mesh, MK = (kK + )I

Figure 5: Example 5.3. Meshes of N = 60 are obtained for Φ(x, y) = 4 sin(x) − y. As discussed in Section 4, theoretically we know that the value of Ih is decreasing and |K| is bounded below. To see these numerically, we plot Ih and |K|min as functions of t in Fig. 6, where |K|min denotes the minimum area of K over all elements in Th . The numerical results are shown to be consistent with the theoretical predictions. Specifically, for MK = I, Fig. 6(a) shows that Ih is decreasing and bounded below by 9.535. Additionally, Fig. 6(b) suggests that |K|min is bounded below by 0.235 which is the value of |K|min of the initial mesh. As we see, |K|min first increases and then converges to about 0.285 ≈ |S| N . The reason is because in the final mesh, the elements are close to being uniform with respect to the Euclidean metric and thus |K| ≈ |S| N for all K. Since the initial mesh is nonuniform, we expect an increase in |K|min as the mesh is becoming more uniform. Moreover, as the mesh reaches the limiting mesh trajectory around t = 0.05, we see that |K|min converges as shown in Fig. 6(b). For the case with MK = (kK + )I, the numerical results are again consistent with the theoretical predictions. Fig. 6(c) shows that Ih is decreasing for all time and bounded below by 15.5. This figure also shows that at around t = 0.15, Ih begins to converge. In Fig. 6 (d), |K|min has similar properties to Fig. 6(b). That is, we see an initial increase in |K|min after which, the value converges to 0.11 starting at around t = 0.15. Furthermore, Fig. 6(d) suggests that |K|min is bounded below by the initial value of 0.045. We show another example of the surface moving mesh method for the sine curve where now we use a gradient-based metric tensor defined as MK = (gK + )I, where Φ2x gK =  y 2 + 100

24

is based on the gradient of Φ and is machine epsilon. We take N = 60 and fix the end nodes x1 = (0, 0) and x61 = (2π, 0). Studying the final mesh, Fig. 7(b), we see a concentration of mesh elements in the regions of the curve with larger gradient and close to y = 0. This adaptation is consistent with the gradient-based metric tensor which has the large determinant in those regions. The accurate adaptation can also be seen in the improvement of the equidistribution quality measure from Qeq = 6.892 to Qeq = 1.052. While the final mesh is close to being uniform in the metric MK , it is very nonuniform in the Euclidean metric as indicated by ν = 399. Finally, we explore a case of the sine curve for which the limiting mesh is not unique. More specifically, Theorem 4.2 implies that although Ih is monotonically decreasing, the mesh trajectory may not converge to a single final mesh. To demonstrate this we take N = 11 and fix the end nodes as before. For this example we use the curvature-based metric tensor. Fig. 8(d) shows the norm of the difference between consecutive meshes, Thn − Thn−1 ∞ ≡ max |xni − xn−1 |, i i

as a function of t, where Thn denotes the mesh at iteration n for n = 1, 2, . . . . Studying Fig. 8(d) we see that the difference converges to about 0.3 indicating that the mesh does not converge to a unique final mesh but instead alternates between two, Fig. 8(b) and Fig. 8(c). The only difference between

Fig. 8(b) and Fig. 8(c) is the coordinate positions of the two nodes closest to π2 , 4 however, this does not seem to affect the mesh quality. The equidistribution quality measure improves from 2.644 to 1.177 for the final mesh in Fig. 8(b) and 1.177 for the one in Fig. 8(c). This implies that although the mesh for this example is not convergent, it does not affect the mesh quality improvement. It is interesting to point out that the convergence improves as the mesh becomes finer. Fig. 9(b) shows the final mesh of N = 22 for this example. The mesh trajectory converges to a single limiting mesh as confirmed in Fig. 9(c) where the difference decreases to the round-off error level. We have also observed similar convergence behavior for sufficiently fine meshes for all of other examples presented in this paper although this is not clear from theory (particularly Theorem 4.2). Example 5.4. This two-dimensional example is the rose defined by     θ θ x(θ) = cos cos(θ), y(θ) = cos sin(θ), θ ∈ [0, 8π]. 4 4 We take N = 220 and use numerically approximated normal vectors. Fig. 10 shows the meshes for this example. We can see that the initial mesh Fig. 10(a) is very nonuniform but the final mesh (Fig. 10(b)) is almost equidistant along the curve which is consistent with the use of the Euclidean metric. This quality improvement is also seen in Qeq from 6.310 to 1.030. Consider now the curvature-based metric tensor (Fig. 10(c)). We can see a high concentration of mesh elements in the regions of the rose with large curvature, i.e., the inner most smaller circles, which is consistent with the use of the curvature-based metric tensor. The improvement of Qeq from 5.942 to 1.001 also indicates that the final mesh is almost uniform with respect to the curvature-based metric tensor.

25

Moreover, this example illustrates that crossings in a curve do not affect the mesh movement. As Fig. 10(c) indicates, the nodes remain on the correct path of the curve and do not become tangled even at the areas where the curve crosses itself. This is consistent with the theoretical property given in Section 4 that states the mesh does not become tangled if it is not so initially. Example 5.5. As the final two-dimensional example, we generate adaptive meshes for the lemniscate defined by Φ(x, y) = (x2 + y 2 )2 − 4(x2 − y 2 ). In this example adapt the mesh on the curve for both N = 60 and N = 120. In both situations, we fix the node x1 = (2, 0). Analytical normal vectors are used in this example. From Fig. 11 we see that for N = 60 the mesh adapts from a very nonuniform initial mesh (Fig. 11(a)) to a uniform final mesh (Fig. 11(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.083 for the initial mesh to 1.002 for the final mesh. We see a similar result when the curvature-based metric tensor is used (Fig. 11(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.364 to 1.001 indicating that the final mesh is much more uniform with respect to the curvature-based MK than the initial mesh. For N = 120, Fig. 12 shows similar findings. When considering the Euclidean metric, we see the mesh, Fig. 12(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.552 to 1.002 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. 12(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. 12(c) with N = 120 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.023 for the initial mesh to 1.002 for the final mesh. Example 5.6. Let us now consider surfaces in R3 . In this first example, we consider adaptive meshes for the torus defined by

2  Φ(x, y, z) = 2 − x2 + y 2 + z 2 − 1, where x, y ∈ [−3, 3], and z ∈ [−1, 1]. We take N = 3200. Numerically approximated normal vectors are used for this example.

26

Fig. 13 shows the meshes for this example in two different views. Studying Fig. 13(a), the initial mesh, and Fig. 13(b), the final mesh with MK = I, 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 M = I 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.50 and for the final mesh 1.332. Similarly, the initial alignment quality measure is 30.63 compared to that of the final mesh which is 1.920. 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. 13(c), looks identical to the final mesh for the Euclidean metric, Fig. 13(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.7. The second three-dimensional example is the cylinder defined by Φ(x, y, z) = x2 + y 2 − 1, where z ∈ [−2, 2]. For this example we take N = 3200. In this example, normal vectors are calculated using an analytical formula. Two boundary nodes were fixed, x1 = (0, 1, −2) and x1 = (0, 1, 2), but the remaining boundary nodes were allowed to slide along the boundary. Although the cylinder has constant curvature like Example 5.6, this example shows the adaptation on a surface with a boundary. Fig. 14 shows the adaptive meshes for the cylinder in two different views. For both MK = I and MK = (kK + )I, 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.07 to 1.054 and the alignment quality measure from 23.35 to 1.192. The fact that the final quality measures for both conditions are close to 1 indicates that the final meshes are close to satisfying conditions (9) and (11). Example 5.8. Our next example is the sine surface in three dimensions defined by Φ(x, y, z) = sin(x + y) − z. For this example we take N = 3200 and fix the boundary nodes. Numerically approximated normal vectors are used in this example. Fig. 15 shows the adaptive meshes for this examples in two different views. It is clear in Fig. 15, when MK = I, the mesh becomes much more uniform with respect to the Euclidean metric from the initial mesh Fig. 15(a) to the final mesh Fig. 15(b). The top view of the surface, Fig. 15(d) and (e), further confirms this observation. It is also supported by the improvement of the quality measures from Qeq = 4.234 and Qali = 6.643 to Qeq = 1.669 and Qali = 1.702. The area ratio ν improves from 13.9 to 1.11, indicating that the final mesh is close to being uniform with respect to the Euclidean metric. When MK is curvature-based, we see a similar result to Example 5.3. That is, Fig. 15(c) and (f) show that the elements are more concentrated in those regions of the surface with larger curvature,

27

i.e., the dip when z = −1 and the hill when z = 1. The quality measures with respect to the metric tensor improve from Qeq = 21.69 to Qeq = 1.634 and Qali = 6.527 to Qali = 2.586. The final quality measure for the equidistribution condition close to 1 hence indicating that the final mesh is close to satisfying (9). The final quality measure for the alignment condition is not as close to 1 as the equidistribution condition. Recall that θ in the meshing function (16) balances equidistribution and alignment and the choice θ = 1/3 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 (16) for the curvature-based metric tensor for this example. We would like to take a look at the changes of Ih and |K|min along the mesh trajectory. As we recall from Section 4, |K| is bounded from below and Ih is decreasing. These can be seen numerically for MK = I in Fig. 16(a) and Fig. 16(b). Similar to what we saw in Example 5.3, Fig. 16(a) shows that Ih is always decreasing and at around t = 0.10 begins to converge. In Fig 16(b) we see an initial increase in the |K|min value and then it begins to converge to 4.64×10−3 ≈ |S| N at t = 0.10. 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 |K|min can be very small whereas when the mesh is adapted, the mesh becomes more uniform and hence the values of |K| ≈ |S| N become almost identical. This implies that the value of |K|min is likely to increase as the mesh adapts. For the case with MK = (kK + )I, Fig. 16(c) and (d) show similar findings. In 16(c) we see that Ih is decreasing for all time and converging beginning at around t = 0.15. Fig. 16(d) shows |K|min initially increases then begins to converge to about 2.0×10−4 . Furthermore, |K|min is bounded below by the initial |K|min value of 0.90 × 10−4 . These numerical results for the curvature based metric tensor further support the theoretical predictions. Finally, to further validate the surface moving mesh method, we show another example for the sine surface now with a gradient-based metric tensor defined as MK = (gK + )I, where gK = ∇Φ 2 is the norm of the gradient and is machine epsilon. We take N = 3200 and fix the boundary nodes. Fig. 17 shows the adaptive mesh for this example in two different views. Studying Fig. 17(b) and (d), the final mesh, we see similar results to Example 5.3 with the gradient-based metric tensor. More specifically, Fig. 17(b) and (d) show that the mesh elements are more concentrated in those regions of the surface with steeper slope and fewer elements in the more planar regions, i.e., the dip when z = −1 and the hill when z = 1. This concentration is further supported by the quality measures which improved from Qeq = 11.39 and Qali = 7.209 to Qeq = 1.700 and Qali = 1.897. The final quality measures are closer to 1 which indicates that the mesh is closer to being uniform with respect to the gradient-based metric tensor. On the other hand, a change of ν defined in (55) from 13.9 to 53.5 implies that the mesh becomes more nonuniform in the Euclidean metric. Example 5.9. In this example we consider the wave surface in three dimensions defined by

  sin x2 + y 2 + 16  Φ(x, y, z) = . x2 + y 2 + 16 For this example we take N = 3872 and fix the boundary nodes. Numerically approximated normal vectors are used in this example.

28

Fig. 18 shows the adaptive meshes for the wave in two different views. For MK = I, the mesh becomes much more uniform throughout the surface as can be seen in Fig. 18(b) and (c). The final mesh quality is further supported in the measures which improve from Qeq = 7.012 and Qali = 9.119 to Qeq = 1.201 and Qali = 1.892. When MK is curvature-based, we see the elements concentrating in the regions with larger curvature, i.e., the peaks and dips of the wave. The quality measures with respect to the curvature-based metric improve from Qeq = 8.671 to Qeq = 1.769 and Qali = 8.996 to Qali = 1.501. These final mesh quality measures are closer to 1 indicating that the final mesh is closer to being uniform in the curvature-based metric tensor. Example 5.10. Our final example explores the sphere and ellipsoid defined by an icosahedral initial mesh (see [28] for more details). We begin with the sphere Φ(x, y, z) = x2 + y 2 + z 2 − 1. For this example we take N = 1280. Analytical normal vectors are used in this example. As we see from Fig. 19(a), the initial mesh is close to being uniform however, there is a very slight difference in the final mesh Fig. 19(b) when MK = I. Indeed, this slight adaptation can be seen in the quality measures which change from Qeq = 1.068 and Qali = 1.025 for the initial mesh to Qeq = 1.289 and Qali = 1.025 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 explore a shock structure on the unit sphere. More specifically, we define the metric tensor using a sine wave near the equator by M = (0.01 + ∇f 2 ) I, where

(56)

 

 y  . f (x, y, z) = tanh 30 z + sin 0.2 sin 6 arctan x Fig. 20 shows the meshes for this example in two different views. One can see in Fig. 20(d) that the mesh elements adapt in accordance with the gradient-based metric tensor as defined above, and particularly, a large concentration of elements occurs along the wave near the equator on the unit sphere. The equidistribution measure changes from Qeq = 1.405 to Qeq = 1.026 whereas the alignment quality measure from Qali = 3.648 to Qali = 13.81. One may notice that the mesh alignment quality decreases. This is because the metric tensor (56) is scalar (which favors isotropic elements) and the icosahedral initial mesh has already a very good alignment quality in the metric tensor. (The L2 norm of the alignment quality measure is 1.129 for the initial mesh and 2.195 for the final mesh.) As a result, any deformation from this initial mesh will almost guarantee a decrease in the alignment quality. Moreover, the surface moving mesh method, balancing between equidistribution and alignment (c.f. (16)), is likely to move the mesh from the initial mesh and improve the equidistribution measure. As we can see in Fig. 20(d), this happened and the elements concentrate around the shock structure. The element compression can be seen to change from ν = 1.299 to ν = 102.5.

29

Now consider adaptive meshes for the ellipsoid defined by Φ(x, y, z) = x2 + y 2 +

z2 − 1. 4

We move the mesh on the surface for both N = 1280 and N = 5120. First considering N = 1280, Fig. 21 shows the meshes for this example in two different views. Studying Fig. 21(a) and Fig. 21(d), the initial mesh, and Fig. 21(b) and Fig. 21 (e), the final mesh with MK = I, 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 Qeq = 1.724 and Qali = 1.453 for the initial mesh to Qeq = 1.571 and Qali = 1.102 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 MK = (kK + )I 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 curvaturebased metric tensor to concentrate more mesh elements at the tips of the ellipsoid by redefining MK as   1 1 ˜ K = MK +  M + I. (57) (zK − 2)2 +

(zK + 2)2 +

Fig. 21(c) and Fig. 21(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 1.374 initially to 1.967 whereas the alignment quality measure from 1.453 to 1.262. Similar results are seen with a finer mesh in Fig. 22 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 (6) 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 (9) and (11), 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 (39) for θ ∈ (0, 1/2] and p > 1. We defined the surface MMPDE (34) 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 (28) and (29) give explicit, compact formulas for the mesh velocities making the time integration of the surface MMPDE (35) 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

30

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 (34) only depends on surface normal vectors which can be computed based on the current mesh. 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 main cases of the metric tensor MK = I and MK = (kK + )I as well as the gradient-based metric tensor. Furthermore, the MMPDE method showed to accurately adapt the mesh to a shock structure on the unit sphere. The successful adaptation to the final mesh in all of the examples was supported by the mesh quality measures. Additionally, the compression measure defined in (55) as the ratio of the area of the largest element to the smallest further indicated the degree to which the method adapted the mesh. The theoretical properties were numerically verified in Example 5.3 and Example 5.8 as we showed that Ih is decreasing and |K| is bounded below. Finally, we explored the case for which the limiting mesh is not unique. These results showed that although the mesh did not converge, it did not affect the mesh quality improvement. The future goal is to investigate the use of spline functions to more accurately approximate the normal vectors required for this method. 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.10.

Appendices A. Derivation of derivatives of the meshing function with respect to the physical coordinates Recall from Section 3.2 that GK = G (JK , rK ) and our objective is to compute the derivatives ∂GK . K K ∂[x1 , xK 2 , . . . , xd ] K K Let t be an entry of [xK 1 , x2 , . . . , xd ]. Using the chain rule we have     ∂GK ∂EK ∂GK ∂MK ∂GK = tr + tr . ∂t ∂EK ∂t ∂MK ∂t

Denote

∂GK (I) = tr ∂t



31

∂GK ∂EK ∂EK ∂t

 ,

(58)

∂GK (II) = tr ∂t



∂GK ∂MK ∂MK ∂t

 .

(59)

K K K K To begin, consider (58). When t is an entry of [xK 2 , . . . , xd ], recalling that EK = [x2 − x1 , . . . , xd − xK 1 ], we have  K  ∂GK ∂[xK ∂GK 2 , . . . , xd ] (I) = tr , ∂t ∂EK ∂t

which implies ∂GK (I) K ∂[x2 , . . . , xK d ]

(1) Moreover, for t = xK (the first component 1 ⎛ ⎡ −1 ⎜ ⎢ ⎜ ∂GK ⎢ 0 ∂GK ⎢ K (1) (I) = tr ⎜ ⎜ ∂EK ⎢ .. ⎝ ⎣ . ∂ x 1

0

=

∂GK . ∂EK

of xK 1 ), we have ⎤⎞ −1 · · · −1 ⎥⎟   ∂GK  0 · · · 0 ⎥⎟ ⎥ ⎟ . .. ⎥⎟ = − .. ∂EK i,1 . ⎦⎠ . i 0 ··· 0

(j) for j = 2, . . . , d. This gives We can have similar expressions for xK 1 ∂GK ∂GK (I) = −eT K ∂EK ∂x1 where eT = [1, . . . , 1] ∈ R1×(d−1) . For (59), we assume that M = M(x) is a piecewise linear function d  K Mj,K φK defined on the current mesh, i.e., M = j , where φj is the linear basis function associated j=1

(i)

(i) and x , with the vertex xK j for all j = 1, . . . , d. Denote the ith components of x and xK by x K K , . . . , xK ], we have respectively. Then, for any entry t of [xK , x 1 2 d

∂GK (II) = tr ∂t

 ⎛

d ∂GK  ∂MK ∂MK ∂x(i) i=1



(i)

∂xK ∂t

⎞ d  d (i)  ∂φ ∂G j,K ⎠ ∂xK K = tr ⎝ Mj,K ∂MK ∂t ∂x(i) =

d  j=1

 tr

i=1 j=1

∂GK Mj,K ∂MK



∂φj,K ∂xK , ∂x ∂t

∂ xK where we notice that ∂ j,K x and ∂t are a row and a column vector, respectively, and thus K is a dot product. From this and the identity xK = (xK 1 + · · · + xd )/d, we get ⎤ ⎡ ∂φj,K   d ⎢ ∂.x ⎥ ∂GK ∂GK 1 ⎢ .. ⎥ (II) = tr M j,K ⎦ ⎣ d ∂M ∂[xK , . . . , xK ] ∂φ

2

d

j=1

32

K

∂φj,K ∂x

∂φj,K ∂ xK ∂x ∂t

and

∂GK 1 (II) = tr d ∂xK 1 j=1 d



∂GK Mj,K ∂MK



∂φj,K . ∂x

Summarizing the above results, we have ∂GK K ∂[x2 , ..., xK d ]

=

d 

⎡

⎤

∂φj,K ⎢ ∂x ⎥

∂GK .. ⎥ , Mj,K ⎢ . ⎦ ⎣ ∂MK ∂φj,K j=1 ∂x   d ∂φj,K ∂GK 1 . + tr Mj,K d ∂MK ∂x

∂GK 1 + ∂EK d

∂GK ∂GK = −eT K ∂EK ∂x1





tr

(60)

(61)

j=1

Notice that (61) can be rewritten as   d d  ∂φj,K ∂GK ∂GK  ∂GK , =− + tr Mj,K K K ∂MK ∂x ∂x1 ∂xj j=2 j=1

(62)

which gives (29). Next, we establish the relations between ∂GK , ∂EK

∂GK ∂MK

and

∂GK , ∂J

∂GK . ∂r

 =E E ˆ −1 , thus First recall that FK K

J=

−1



T

−1 T  T  ˆ EK ˆ . M K FK =E M K EK FK E

(63)

Let EK = EK (t). Then we have ⎛

−1 ⎞ −1 T    )T M F  ∂ (F ∂ det (F ) M F K K K K K K ∂GK ⎜ ∂GK ⎟ ∂GK = tr ⎝ ⎠+ ∂t ∂J ∂t ∂r ∂t  = tr



−1

−1 T T ∂ det E M K EK M E ∂G ∂GK ˆ ∂ EK K K K K T 2 ˆ ˆ E E . + det(E) ∂J ∂t ∂r ∂t

Consider the first term of (64). Using the properties of matrix derivatives (26) we get   T

−1 M K EK ∂GK ˆ ∂ EK ˆT E E tr ∂J ∂t  

T

−1 ∂ EK

−1 T M K EK T ∂GK ˆ T ˆ = −tr E E K M K EK E K M K EK E ∂J ∂t   T  

−1 ∂EK

−1 T ∂GK ˆ T ∂EK T T ˆ E E K M K EK MK E K + EK M K EK MK EK . = −tr E ∂J ∂t ∂t

33

(64)

Since

∂GK ∂J ,

T

−1 MK , and EK M K EK are all symmetric, it follows that   T

−1 M K EK ∂GK ˆ ∂ EK ˆT tr E E ∂J ∂t   T

−1 T ∂GK T

∂EK T ˆ EK MK EK −1 EK ˆ = −2tr EK E . M K EK MK E ∂J ∂t

The second term of (64) is then T

−1 ∂ det E M E ∂G K K K K 2 ˆ det(E) ∂r ∂t 

−1  T 2 ˆ T

∂ EK M K EK det(E) ∂GK T

= tr EK MK EK ∂r ∂t det EK M K EK  

TM E ˆ 2

−1 ∂ EK ∂GK det(E) K K T T

=− tr E K M K EK ∂r ∂t det EK M K EK  T   ˆ 2

−1 ∂EK ∂GK ∂EK T det(E) T

T =− tr MK E K + EK M K EK MK EK ∂r ∂t ∂t det EK M K EK   ˆ 2 T

−1 T ∂GK ∂EK det(E)

= −2 tr EK MK EK . EK M K TM E ∂r ∂t det EK K K Therefore T

−1 T ∂GK ∂GK T T ˆ K ˆ E(E = −2 EK M K EK MK EK )−1 EK MK E ∂EK ∂J ˆ 2

−1 T ∂GK T det(E) T

E K M K EK −2 EK M K . det EK MK EK ∂r Combining this with (61) we obtain (28). The identity (30) can be obtained similarly. Finally, we derive the relations between ∂φj,K , ∂x

j = 1, ..., d

and

EK .

First, note that the basis functions satisfy d 

φi,K = 1

and

i=1

d 

xK i φi,K = x.

i=1

Eliminating the xK 1 yields x − xK 1 =

d 

K (xK i − x1 )φi,K .

i=2

Then differentiating with respect to ek =

x(k)

∂(x − xK ∂ 1 ) = (k) ∂x ∂x(k)

gives  d  d   K K ∂φi,K = (xK − x )φ (xK , i,K i 1 i − x1 ) ∂x(k) i=2 i=2

34

(65)

where ek is the k th unit vector in Rd . Hence we have ⎡

⎤

∂φ2,K ⎢ ∂x ⎥

I = EK ⎢ ⎣ which gives

.. .

∂φd,K ∂x

⎥, ⎦

⎤

⎡

∂φ2,K ⎢ ∂x ⎥

T EK EK ⎢ ⎣

.. .

∂φd,K ∂x

T ⎥ = EK ⎦

and thus (31).

References [1] B. L´evy and N. Bonneel. Variational anisotropic surface meshing with Voronoi parallel linear enumeration. The 21st International Meshing Roundtable, Springer-Verlag, 344-366, (2012). [2] J.-D. Boissonnat, F. Chazal, and M. Yvinec. Geometric and Topological Inference. Cambridge University Press, (2018). [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´ere type equation. J. Comput. Phys., 308:102-123, (2016). [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] B. Crestel, R. D. Russell, and S. Ruuth. Moving mesh methods on parametric surfaces. Procedia Engineering, 124:148-160, (2015). [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] 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] G. Dziuk. Finite elements for the Beltrami operator on arbitrary surfaces. Springer, New York. Lecture Notes in Math., Vol. 1357, (2006). [9] G. Dziuk and C. Elliott. Finite element methods for surface PDEs. Acta Numerica, 22:289-396, (2013). [10] G. Dziuk and C. Elliott. Surface finite elements for parabolic equations. J. Comput. Math., 25:385-407, (2007).

35

[11] J. Emert and R. Nelson. Volume and surface area for polyhedra and polytopes. Math. Mag., 70:365-371, (1997). [12] F. Goldman. Curvature formulas for implicit curves and surfaces. Comp. Aided Geometric Design., 22:632-658, (2005). [13] W. Huang. Variational mesh adaptation: isotropy and equidistribution. J. Comput. Phys., 174:903-924, (2001). [14] W. Huang. Metric tensors for anisotropic mesh generation. J. Comput. Phys., 204:633-665, (2005). [15] W. Huang and L. Kamenski. On the mesh nonsingularity of the moving mesh PDE method. Math. Comp., 87:1887-1911, (2018). [16] W. Huang and L. Kamenski. A geometric discretization and a simple implementation for variational mesh generation and adaptation. J. Comput. Phys., 301:322–337, (2015). [17] W. Huang, L. Kamenski, and R. D. Russell. A comparative numerical study of meshing functions for variational mesh adaptation. J. Math. Study, 48:168-186, (2015). [18] W. Huang, L. Kamenski, and H. Si. Mesh smoothing: an MMPDE approach. Research note at the 24th International Meshing Roundtable (2015). [19] W. Huang and R. D. Russell. Moving mesh strategy based upon a gradient flow equation for two dimensional problems. SIAM J. Sci. Comput., 20:998-1015, (1999). [20] W. Huang and R. D. Russell. Adaptive Moving Mesh Methods. Springer, New York. Applied Mathematical Sciences Series, Vol. 174, (2011). [21] W. Huang, Y. Ren, and R. D. Russell. Moving mesh partial differential equations (MMPDEs) based upon the equidistribution principle. SIAM J. Numer. Anal., 31:707-730, (1994). [22] A. Kolasinski. Surface and bulk moving mesh methods based on equidistribution and alignment. Department of Mathematics, University of Kansas, PhD dissertation, (2019). [23] A. Kolasinski and W. Huang. A new function for variational mesh generation and adaptation based on equidistribution and alignment conditions. Comput. Math. Appl., 75:2044-2058, (2018). [24] G. MacDonald, J. A. Mackenzie, M. Nolan, and R. H. Insall. A computational method for the coupled solution of reaction-diffusion equations on evolving domains and manifolds: Application to a model of cell migration and chemotaxis. J. Comput. Phys., 309:207-226, (2016). [25] J. A. Mckenzie, M. Nolan, C. F. Rowlatt, and R. H. Insall. An adaptive moving mesh method for forced curve shortening flow. SIAM J. Sci. Comput., 41:A1170-A1200 (2019). [26] A. T. T. McRae, C. J. Cotter, and C. J. Budd. Optimal-transport-based mesh adaptivity on the plane and sphere using finite elements. SIAM J. Sci. Comput., 40:A1121-A1148, (2018).

36

[27] N. Tuncer and A. Madzvamuse. Projected finite elements for systems of reaction-diffusion equations on closed evolving spheroidal surfaces. Comm. Comp. Phys., 21:718-747, (2017). [28] N. Wang and J. Lee. Geometric properties of the icosahedral-hexagonal grid on the two-sphere. SIAM J. Sci. Comput., 33:2536-2559, (2011).

37

9.61 0.28

9.6

0.275

9.59 0.27

9.58 0.265

9.57

0.26

9.56

0.255

9.55

0.25 0.245

9.54

0.24

9.53

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

(a) Ih , MK = I

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.8

0.9

1

(b) |K|min , MK = I

18.5 0.1

18

0.09

0.08

17.5

0.07

17

0.06

16.5

0.05

16

15.5

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

(d) |K|min , MK = (kK + )I

(c) Ih , MK = (kK + )I

Figure 6: Example 5.3. Ih and Kmin plotted as functions of t for Φ(x, y) = 4 sin(x) − y.

38

3

3

2

2

1

1

0

0

-1

-1

-2

-2

-3

-3

-1

0

1

2

3

4

5

6

7

-1

8

0

1

2

3

4

5

6

7

8

(b) Final Mesh, MK = (gK + )I

(a) Initial Mesh

Figure 7: Example 5.3. Meshes of N = 60 for the curve Φ(x, y) = 4 sin(x) − y with a gradient-based metric tensor.

3

3

3

2

2

2

1

1

1

0

0

0

-1

-1

-1

-2

-2

-2

-3

-3

-3

-1

0

1

2

3

4

5

6

(a) Initial Mesh

7

8

-1

0

1

2

3

4

5

6

7

10

-1

8

(b) Final Mesh 1

0

1

2

3

4

5

6

7

(c) Final Mesh 2

8

10

0

-1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

(d) Thn − Thn−1 ∞

Figure 8: Example 5.3. Meshes of N = 11 are obtained for Φ(x, y) = 4 sin(x) − y with MK = (kK + )I and Thn − Thn−1 ∞ plotted as functions of t.

39

10 0

3

3

2

2

1

1

0

0

-1

-1

-2

10

-10

-2

-3

-3

10

-12

10

-14

10 -2

10 -4

10 -6

-1

0

1

2

3

4

5

6

7

10

8

-1

0

(a) Initial Mesh

1

2

3

4

5

6

7

8

-8

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

(c) Thn − Thn−1 ∞

(b) Final Mesh

Figure 9: Example 5.3. Meshes of N = 22 are obtained for Φ(x, y) = 4 sin(x) − y with MK = (kK + )I and Thn − Thn−1 ∞ plotted as functions of t.

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

-0.2

-0.2

-0.2

-0.4

-0.4

-0.4

-0.6

-0.6

-0.6

-0.8

-0.8

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

(a) Initial Mesh

0.8

1

-0.8

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

(b) Final Mesh, MK = I

-1

1

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

(c) Final Mesh, MK = (kK + )I

Figure 10: Example 5.4 Meshes for the flower defined by x(θ) = cos for θ ∈ [0, 8π] with N = 220.

40

-0.8

θ

4

cos(θ), y(θ) = cos

θ

4

sin(θ)

1.5

1.5

1.5

1

1

1

0.5

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-1

-1

-1

-1.5

-1.5 -1.5

-1

-0.5

0

0.5

1

1.5

-1.5

-1.5

2

-1

-0.5

0

0.5

1

1.5

(b) Final Mesh, MK = I

(a) Initial Mesh

-1.5

2

-1

-0.5

0

0.5

1

1.5

2

(c) Final Mesh, MK = (kK + )I

Figure 11: Example 5.5. Meshes of N = 60 are obtained for the lemniscate Φ(x, y) = (x2 + y 2 )2 − 4(x2 − y 2 ).

1.5

1.5

1.5

1

1

1

0.5

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-1

-1

-1

-1.5

-1.5

-1.5

-1.5

-1

-0.5

0

0.5

1

(a) Initial Mesh

1.5

2

-1.5

-1

-0.5

0

0.5

1

1.5

(b) Final Mesh, MK = I

2

-1.5

-1

-0.5

0

0.5

1

1.5

2

(c) Final Mesh, MK = (kK + )I

Figure 12: Example 5.5. Meshes of N = 120 are obtained for the lemniscate Φ(x, y) = (x2 + y 2 )2 − 4(x2 − y 2 ).

41

4

3

3 2 1 0

0.5 0

-2

-1

2 1

0

0.5

-1

0

0

0.5

-1

0

-2

-2 -0.5

-0.5

2 1

-3

-0.5

-3

-4 -4

-3

-2

-1

0

1

2

3

-3

4

-2

-1

0

1

2

-3

3

(b) Final Mesh, MK = I

(a) Initial Mesh

-2

-1

0

1

2

3

(c) Final Mesh, MK = (kK + )I

1 0.5

0.5

0

0

0

-0.5

-0.5

-0.5

0.5

-1 -4

-3

-2

-1

0

1

2

(d) side view of (a)

3

4

-3

-2

-1

0

1

2

(e) side view of (b)

3

-3

-2

-1

0

1

2

3

(f) side view of (c)

Figure 13: Example 5.6. Meshes of N = 3200 are obtained for the surface Φ(x, y, z) = (2 −  2 2 x + y )2 + z 2 − 1.

42

2

2

2

1.5

1.5

1.5

1

1

1

0.5

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-1

-1

-1

-1.5

-1.5

-1.5

-2

-2

-2

0.5

0.5

0.5

0

0

0 -0.5 -0.5

0

-0.5

0.5

(a) Initial Mesh

-0.5

0

-0.5

0.5

(b) Final Mesh MK = I

-0.5

2

2

1.5

1.5

1.5

1

1

1

0.5

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-1

-1

-1

-1.5

-1.5

-1.5

-2

-2

-0.5

0

0.5

(d) side view of (a)

0.5

(c) Final Mesh MK = (kK + )I

2

-2

0

-0.5

0

0.5

(e) side view of (b)

-0.5

0

0.5

(f) side view of (c)

Figure 14: Example 5.7. Meshes of N = 3200 are plotted for Φ(x, y, z) = x2 + y 2 − 1.

43

1

1

1

0.5

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-1

-1 0

-1 2

2.5

0

-1

1 3

3.5

4

2

2.5

2

4.5

3.5

4

4.5

2

2

(b) Final Mesh, MK = I

(a) Initial Mesh

-1 0

-1

1 3

4.5

4.5

4

4

4

3.5

3.5

3.5

3

3

3

2.5

2.5

2.5

2

2

2

-1.5

-1

-0.5

0

0.5

1

1.5

(d) top view of (a)

2

-2

-1.5

-1

-0.5

0

0.5

1

1.5

(e) top view of (b)

1 3

3.5

4

4.5

2

(c) Final Mesh, MK = (kK + )I

4.5

-2

2.5

2

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

(f) top view of (c)

Figure 15: Example 5.8. Meshes of N = 3200 for the surface Φ(x, y, z) = sin(x + y) − z.

44

10

-3

1981.6

4.62 4.6

1981.4

4.58 1981.2

4.56 4.54

1981

4.52 4.5

1980.8

4.48 1980.6

4.46 4.44

1980.4

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

(a) Ih , MK = I

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.8

0.9

1

(b) |K|min , MK = I 10

6400

-4

1.8 6200

1.7 1.6

6000

1.5

5800

1.4 5600

1.3 5400

1.2 5200

1.1 5000

1

4800 4600

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

(d) |K|min , MK = (kK + )I

(c) Ih , MK = (kK + )I

Figure 16: Example 5.8. Ih and Kmin are plotted as functions of t for Φ(x, y, z) = sin(x + y) − z.

45

1

1

0.5

0.5

0

0

-0.5

-0.5

-2

-1

-1

0

-1 2

2.5

-1

1 3

3.5

4

0 2

2

4.5

4.5

4.5

4

4

3.5

3.5

3

3

2.5

2.5

2

2

-1.5

-1

-0.5

0

0.5

1 3

3.5

4

4.5

2

(b) Final Mesh, MK = (gK + )I

(a) Initial Mesh

-2

2.5

1

1.5

-2

2

(c) top view of (a)

-1.5

-1

-0.5

0

0.5

1

1.5

2

(d) top view of (b)

Figure 17: Example 5.8. Meshes of N = 3200 for the surface Φ(x, y, z) = sin(x + y) − z with a gradient-based metric tensor.

46

0.15

0.15

0.15

0.1

0.1

0.1

0.05

0.05

0.05

0

0

0

-0.05

-0.05

-0.05

-0.1

-0.1

-0.1

-0.15

-0.15

-0.15

-0.2

-0.2

-0.2

-0.25 20

-0.25 20

-0.25 20

10

10 0

10 0

-10 -20

-5

-10

-15

0

5

15

10

0 -10 -20

-5

-10

-15

0

5

15

10

-20

(b) Final Mesh, MK = I

(a) Initial Mesh

-10

15

15

10

10

10

5

5

5

0

0

0

-5

-5

-5

-10

-10

-10

-10

-5

0

5

10

(d) aerial view of (a)

15

-15 -15

-10

-5

0

5

10

15

0

5

-15 -15

-10

-5

0

(e) aerial view of (b)

5

10

(f) aerial view of (c) √  sin x2 +y 2 +16 Figure 18: Meshes of N = 3872 are obtained for Φ(x, y, z) = √ 2 2 . x +y +16

47

15

10

(c) Final Mesh, MK = (kK + )I

15

-15 -15

-5

-10

-15

15

1

0.8 0.6

0.5

0.4 0.2 0

0 -0.2 -0.4

-0.5

-0.6 -0.8

-1 1 0.5

1 0

0.5

0.5

0

0

-0.5

-0.5 -1

0.5 0

-0.5

-1

-0.5

(b) Final Mesh MK = I

(a) Initial Mesh 1 0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

-0.2

-0.2

-0.4

-0.4

-0.6

-0.6

-0.8

-0.8

-1 -1

-0.5

0

0.5

1

-0.8

(c) top view of (a)

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

(d) top view of (b)

Figure 19: Example 5.10. Meshes of N = 1280 are plotted for Φ(x, y, z) = x2 + y 2 + z 2 − 1.

48

1

1

0.5

0.5

0

0

-0.5

-0.5

-1 1

-1 1 0.5

0.5

1 0.5

0 -0.5

0

-0.5

-0.5 -1

0.5

0

0

-0.5 -1

-1

(b) Final Mesh M = (0.01 + ∇f 2 ) I

(a) Initial Mesh 1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

-0.2

-0.2

-0.4

-0.4

-0.6

-0.6

-0.8

-0.8 -1

-1 -1

-0.5

0

0.5

-0.8

1

(c) side view of (a)

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

(d) side view of (b)

Figure 20: Example 5.10. Meshes of N = 5120 are plotted for Φ(x, y, z) = x2 + y 2 + z 2 − 1.

49

2

2

1.5

1.5

1

1

0.5

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-1

-1

-1

-1.5

-1.5

-1.5

-2 1

-2

1.5 1

1

0

0.5

0.5

0

0 -1

-0.5

-1

0

-0.5

(b) Final Mesh MK = I

(a) Initial Mesh

0

0.5

-0.5

-0.5

0

0.5

˜ K in (57) (c) Final Mesh M

1 0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

-0.2

-0.2

-0.2

-0.4

-0.4

-0.4

-0.6

-0.6

-0.6

-0.8

-0.8

-1 -1

-0.5

0

0.5

1

(d) top view of (a)

-0.8

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

(e) top view of (b)

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

(f) top view of (c)

Figure 21: Example 5.10. Meshes of N = 1280 are plotted for Φ(x, y, z) = x2 + y 2 +

50

z2 − 1. 4

2 2

2

1.5

1.5

1

1

0.5

0.5

0

0

-0.5

-0.5

-1

-1

-1

-1.5

-1.5

-1.5

-2 1

-2

-2 1

1.5

1

0

1 0.5 0 -0.5

0.5

-1

1

0

0

0

0

-0.5

-1

-1

-1

(b) Final Mesh MK = I

(a) Initial Mesh

-1

0.5

˜ K in (57) (c) Final Mesh M 1

1 0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

-0.2

-0.2

-0.2

-0.4

-0.4

-0.4

-0.6

-0.6

-0.6

-0.8

-0.8

-1 -1

0

-0.5

-0.5

0

0.5

1

(d) top view of (a)

-1

-0.8

-0.5

0

0.5

1

(e) top view of (b)

-1 -1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

(f) top view of (c)

Figure 22: Example 5.10. Meshes of N = 5120 are plotted for Φ(x, y, z) = x2 + y 2 +

51

z2 − 1. 4