A local rezoning and remapping method for unstructured mesh

Computer Physics Communications 182 (2011) 1361–1376

Contents lists available at ScienceDirect

Computer Physics Communications www.elsevier.com/locate/cpc

A local rezoning and remapping method for unstructured mesh ✩ Zhiwei Lin ∗ , Shaoen Jiang, Shunchao Wu, Longyu Kuang Laser Fusion Research Center, Chinese Academy of Engineering Physics, Mianyang, Sichuan 621900, China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 22 January 2010 Received in revised form 21 June 2010 Accepted 8 November 2010 Available online 3 March 2011 Keywords: ALE methods Rezoning Remapping Mesh closure Mesh union Unstructured mesh ICF simulation

This paper presents explicit description of the rezoning and remapping method for unstructured mesh. The rezoning and remapping constitute two of the three phases of Arbitrary Lagrangian–Eulerian (ALE) methods, while the other one is Lagrangian phase. Our method is local (vs. global) in order that the majority of mesh remains unchanged, and that the error caused by remapping can be confined within the changed areas. This is achieved by specifying the local worse areas of mesh based on the meshclosure and mesh-union concepts. After the computing of worse areas, they are rezoned and united with the remaining areas to compose the new mesh. Then within the new mesh only the worse areas are remapped using four methods with different orders of accuracy, including linear interpolation, particle remapping, first order integral remapping and high order ENO (essentially non-oscillatory) remapping. (However, our program only uses the first order integral remapping method for the moment.) Our method’s application in the inertial confinement fusion (ICF) simulation is introduced at the end. Program summary Program title: Maxis Catalogue identifier: AEIP_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEIP_v1_0.html Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 285 253 No. of bytes in distributed program, including test data, etc.: 13 384 369 Distribution format: tar.gz Programming language: Matlab Computer: PC Operating system: Windows/Linux/Unix RAM: 50 MB Word size: 32 bits Classification: 12, 19.7 Subprograms used: Cat Id Title Reference AECV_v1_0 MULTI2D CPC 180 (2009) 977 Nature of problem: Although the great deformation of mesh elements moved with the fluid, which causes the severe shortening of computing time step and the increasing of computing error, makes it indispensable that the rezoning of distorted mesh and the remapping of physical variables such as density or velocity are performed, it still bothers people that the rezoning and remapping also bring in the loss of Lagrangian motion information of old mesh. Therefore the essential problem is how to minimize the error caused by remapping while improving the mesh quality. Solution method: The best method would be based on the remapping error estimation. Unfortunately, in practice, this method would be quite intricate. Instead of the error estimation we adopt the local rezoning and remapping technique to confine the error within the worse or ’not very important’ areas, that is, only the worse areas are rezoned and remapped which are determined by the mesh-closure or mesh-union method. Restrictions: It is sometimes annoying that there is no uniform criterion to distinguish the worse mesh elements from the good ones.

✩ This paper and its associated computer program are available via the Computer Physics Communications homepage on ScienceDirect (http://www.sciencedirect.com/ science/journal/00104655). Corresponding author. E-mail address: [email protected] (Z. Lin).

*

0010-4655/$ – see front matter doi:10.1016/j.cpc.2010.11.034

©

2011 Elsevier B.V. All rights reserved.

1362

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

Additional comments: Since MULTI2D created by R. Ramis is used to fulfill the Lagrangian phase of the simulation of the ICF process, while our programs are called to perform the rezoning and remapping phases, it is required for non-MULTI2D users that our routines should be slightly reprogrammed. Running time: 1967 seconds for my example. © 2011 Elsevier B.V. All rights reserved.

1. Introduction

In numerical simulations of fluid movement, nowadays there are three popular addressing methods: a Lagrangian method in which the mesh moves with the fluid together, an Eulerian method in which the fluid moves through the fixed mesh, and an Arbitrary Lagrangian–Eulerian (ALE) method in which the best properties of the Lagrangian and Eulerian methods are combined [1]. The ALE method comprises three phases: a Lagrangian phase in which the mesh is repositioned corresponding to the Lagrangian motion of fluid and the physical variables are updated, a rezoning phase in which the nodes of the computational mesh are moved to more optimal positions to improve the mesh quality, and a remapping phase in which the physical variables are transferred to the new mesh. In this paper we will focus on the rezoning and remapping phases. In the literatures on rezoning methods there are two broad classes of them: velocity-based and coordinate (mesh)-based. Also the rezoning methods can be divided into occasional and continuous rezoning depending on the frequency with which the rezoning phase is performed. We refer to [4] for more detailed review on the rezoning methods which also demonstrates an excellent strategy for the continuous rezoning method. Although the occasional rezoning is less accurate than the continuous one, our rezoning method is occasional mainly based on the following considerations:

Fig. 1. The flowchart of our method.

2. The criteria for worse areas (1) The rezoning and remapping can bring in the loss of Lagrangian motion information of old mesh and cause computing error. Therefore they are not executed unless necessary. (2) Usually the continuous rezoning is much more expensive in the calculating time aspect than the occasional one because the former needs to be implemented every several time steps. (3) The ALE methods often apply to the situation in which multifluids are inter-pushed and the interfaces between them should be well maintained. It is easier for the occasional rezoning to keep the boundaries. Our rezoning method also employs the idea of remeshing, and the mesh connectivities can be different before and after the rezoning. The related literatures on remeshing can be found in [11–15]. The remapping algorithms abound in literatures [5–10,16,21– 25] and are maturer than the rezoning ones. Essentially they all interpolate the physical variables from the old mesh onto the new one, but with different orders of accuracy. Four remapping methods are investigated and compared in this paper, including linear interpolation, particle remapping, first order integral remapping and high order ENO remapping. Since any two-dimensional polygon can be partitioned into the union of some triangles, also for convenience we will only consider triangular meshes in the following sections, where we first build the criteria for determining the worse areas within a given mesh, then discuss issues on specifying, rezoning and remapping these worse areas, and finally apply our methods to the ICF simulation. Fig. 1 illustrates the flowchart of our method.

Whichever the continuous rezoning or occasional one we choose, the following problems always arise: when to rezone, what to rezone, how to rezone and how to remap. Since our aim of performing rezoning and remapping is to improve the accuracy and efficiency of the simulation when the mesh is distorted, the ultimate criterion must be tied to the solution error and the computing time step, that is, the error between the new solution after the rezoning and remapping and the solution to the original equations is below some acceptable level, and the new computing time step after the rezoning and remapping is satisfyingly long. However, the remapping error estimation can be so intricate that we won’t follow this approach. It is the ‘good’ geometrical properties of the new mesh and the fact that our remapping is only performed within the worse areas that we rely on to guarantee the small remapping error and long computing time step. Our treatment to the problem about when to rezone is simple: if the computing time step is less than some constant, then the mesh needs to be rezoned.1 The problem about what to rezone is equivalent to the one whether a triangle within a given mesh is good or bad. Our criteria for labeling the triangular elements within a mesh can be manual or automatic, introduced respectively as follows. The further specification of worse areas that need to be rezoned is provided in the next section, followed by another two sections which answer the remaining two problems about how to rezone and how to remap.

1

For example, dt < 10−14 s, where dt stands for the computing time step.

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

1363

2.1. The polygon-given criterion The basic idea about the polygon-given criterion is that people can easily tell the worse areas within a given heavily-distorted mesh and so can encircle them with a simple polygon (see Definition 3.2.1).2 Let Ω denote the simple polygon, then we can make the following rule:

if a triangle q1 q2 q3 ⊆ Ω,

then q1 q2 q3 is bad.

The advantage of this criterion is that the worse areas are often encircled exactly, while the disadvantage is that it is manual. Fig. 2. A simple triangular mesh.

2.2. The CFL-condition and shape-size metric criterion Besides the polygon-given criterion, we can propose other miscellaneous automatic criteria in practice. The CFL-condition and shape-size metric criterion is examined in our method. The CFL-condition (short for Courant–Friedrichs–Lewy) calculates the longest allowable computing time step for each triangle as follows:

dt C (q1 q2 q3 ) = λ(q1 q2 q3 )/C s (q1 q2 q3 ), (1)  2 2 12 where λ(q1 q2 q3 ) = I + − ( I − + I xy ) (see [31]), I + = ( I x + I y )/2, I − = ( I x − I y )/2, I x = x221 + x231 − xc2 /3, I y = y 221 + y 231 − y c2 /3, I xy = x21 y 21 + x31 y 31 − xc y c /3, x21 = x2 − x1 , x31 = x3 − x1 , xc = x21 + x31 , y 21 = y 2 − y 1 , y 31 = y 3 − y1 , y c = y 21 + y 31 , if ∂p ∂ρ

qi = (xi , y i ), i = 1, 2, 3, and C s (q1 q2 q3 ) =

+

∂p p ∂ e /ρ 2 if the

pressure p, the density ρ , the mass-average internal energy e, and the partial derivatives of q1 q2 q3 can be obtained from the EOS (equation of state) p = p (ρ , e ). The shape quality metric and relative size metric measure the shape and size qualities for each triangle respectively as follows [2,3]: let f shape (q1 q2 q3 , q01 q02 q03 ) represent the shape quality of q1 q2 q3 with the reference triangle q01 q02 q03 , then





 





f shape q1 q2 q3 , q01 q02 q03 = 2/  AW −1  ·  W A −1  , where



A=

x2 − x1 y2 − y1

 W =

x3 − x1 y3 − y1



x02 − x01

x03 − x01

y 02 − y 01

y 03 − y 01

if qi = (xi , y i ),

q0j

=

(x0j ,

(2)

,

y 0j ),

,

and

T  =



t i2j

12

f size q1 q2 q3 , q01 q02 q03

1i , j 2

i , j = 1, 2, 3;



    = min S (q1 q2 q3 )/ S q01 q02 q03 , 1 ,

(3)

where S (q1 q2 q3 ) = | det( A )|/2, S (q01 q02 q03 ) = | det( W )|/2. Let  be an equilateral triangle, then we can make the following rules: (1) if dt C (q1 q2 q3 ) < c 1 , then q1 q2 q3 is bad; 2

In Matlab the function ginput () can fulfill the object.

where c 1 , c 2 , c 3 are constants selected in [0, 1] and q01 q02 q03 can take the triangular element of the initial mesh. The advantage of this criterion is that it is automatic, while the disadvantage is that sometimes the worse areas are unsatisfactorily labeled. 3. The specifying of worse areas In this section we will offer two methods for the specifying of worse areas with the polygonal appearances3 corresponding to the two criteria in the previous section. The first one can be regarded as a particular case of the second while the second one as the natural generalization of the first. But just before stretching out our argument it is necessary to review the data structure for triangular meshes which are properly defined, that is, any two overlapping triangles of the mesh can only have either a common edge or a common vertex. We will only consider the properly defined triangular meshes in the following parts. 3.1. The data structure for triangular meshes For the sake of saving computer memory and convenience, a triangular mesh is usually stored in three vectors: x, in the format [x1 x2 x3 · · ·], which is a real number array consisting of x-coordinates of all nodes of the mesh, y, in the format [ y 1 y 2 y 3 · · ·], which consists of y-coordinates, and tn (short for triangle-node connectivity), in the format [qt11 qt21 qt31 qt12 qt22

qt32 · · ·], which is an integer array consisting of indices into x and y of the coordinates of three vertices of each mesh triangular element. For the triangular mesh showed in Fig. 2, we can have that x = [1 1.5 0 0], y = [0 1.5 1 0] and tn = [1 2 3 1 3 4].

,

let f size (q1 q2 q3 , q01 q02 q03 ) express the size quality of q1 q2 q3 with the reference size of q01 q02 q03 , then



(2) if f shape (q1 q2 q3 , ) < f shape (q01 q02 q03 , ) and f shape (q1 q2 q3 , q01 q02 q03 ) < c 2 , then q1 q2 q3 is bad; (3) if f size (q1 q2 q3 , q01 q02 q03 ) < c 3 , then q1 q2 q3 is bad;

3.2. The mesh-closure of a simple polygon within a triangular mesh This subsection provides the method for the specifying of worse areas using the polygon-given criterion. For convenience, the polygon (or line segment) in the following parts means the union of its interior and boundary, while the triangular mesh means the set composed by its all triangular elements (triangles). Definition 3.2.1. A polygon is simple if it is homeomorphous to a circle. Simply speaking, a simple polygon Ω = p 1 p 2 p 3 · · · pn (more rigorous notation is Ω = p 1 p 2 p 3 · · · pn p 1 ) meets the conditions 3

The vertices of the polygon are exactly the given mesh nodes.

1364

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

Fig. 4. A special case in the mesh-closure-contracting algorithm where q3 = q1 → prev and q3 = q2 → next and q3 ∈ Vertex(Ω1 ) and q1 q2 q3 ∩ int(Ω) = ∅.

Fig. 3. An example for the mesh-closure.

that p 1 , p 2 , p 3 , . . . , pn are different and that two inconsecutive edges do not intersect. Also in the following discussion Edge(Ω) = { the line segment p i p i +1 : p i ∈ Vertex(Ω)} and Edge( M ) = ∈ M Edge(), Vertex( M ) = ∈ M Vertex() where M is a triangular mesh. Definition 3.2.2. A polygon Ω1 is the mesh-closure of a given simple polygon Ω within a triangular mesh M if it satisfies: (1) (2) (3) (4)

Ω1 is simple; Ω1 ⊇ Ω ; Edge(Ω1 ) ⊆ Edge( M ); Ω1 is the minimum (in the sense of inclusion or size) among the polygons meeting the above three conditions.

The mesh-closure exists uniquely if Ω ⊆ ∂ M (the boundary ¯ M . Fig. 3 provides an example polygon of M), and is denoted by Ω| for mesh-closure. Note the difference between the mesh-closure and the closure of a simple polygon, whereas

the latter is the ¯ = ( l∈Edge(Ω) l) ∪ int(Ω), union of its interior and boundary (e.g., Ω where int(Ω) is the set of points which lie strictly inside Ω ). The algorithm for computing the mesh-closure of a given simple polygon within a triangular mesh is proposed in pseudo-code. Mesh-closure problem Input: A simple polygon Ω = p 1 p 2 p 3 · · · pn and a triangular mesh M such that Ω ⊆ ∂ M and ∂ M is also simple; ¯ M = Output: The mesh-closure of Ω within M which is Ω| q 1 q 2 q 3 · · · qm . MESH-CLOSURE-CONTRACTING( M , Ω) (1) Determine the relations between Ω and every edge of M as follows: A 1 A 2 ∈ Edge( M ), flag( A 1 A 2 ) = 0 if and only if A 1 A 2 ∩ int(Ω) = ∅ (int(Ω) is for the interior of Ω ); (2) let Ω1 = ∂ M (in counterclockwise order) and q1 be the Y – X ordering minimum of Vertex(Ω1 )4 ; (3) let ptr = q1 and n(q) = 0, ∀q ∈ Vertex( M ); (4) while (∃q ∈ Vertex(Ω1 ), n(q) = 0) (5) do q1 = ptr, q2 = q1 → next (in counterclockwise order within Vertex(Ω1 )), and q3 takes the unique point q such that q ∈ Vertex( M ), Edge(q1 q2 q) ⊆ Edge( M ) and q1 q2 q ⊆ Ω1 5 ; (there are four circumstances discussed below)

4 5

Actually the initial value of q1 can be any point within Vertex(Ω1 ). See Fig. 8 for the method of choosing q3 .

Fig. 5. A special case in the mesh-closure-contracting algorithm where q3 ∈ / Vertex(Ω1 ) and q1 q2 q3 ∩ int(Ω) = ∅.

(6) case 1: q3 = q1 → prev and q3 = q2 → next and q3 ∈ Vertex(Ω1 ),6 then: if q1 q2 q3 ∩ int(Ω) = ∅ (usually ⇔ flag(q1 q2 ) = flag(q2 q3 ) = flag(q1 q3 ) = 0), then Ω1 can be partitioned into two simple sub-polygons Ωi , Ωii besides q1 q2 q3 , where Ωi = [q1 , q3 , . . . , q1 → prev], Ωii = [q3 , q2 , q2 → next, . . .], and there is only one of them that covers Ω denoted by Ω  ,7 therefore let Ω1 = Ω  , and ptr = q1 , if Ω  = Ωi , or ptr = q2 , if Ω  = Ωii ; else (which means q1 q2 q3 ∩ int(Ω) = ∅) let n(q1 ) = 1, ptr = q2 ; (7) case 2: q3 ∈ / Vertex(Ω1 ), then: if q1 q2 q3 ∩ int(Ω) = ∅, modify Ω1 by inserting q3 into between q1 and q2 8 ; else let n(q1 ) = 1, ptr = q2 ; (8) case 3: q3 = q2 → next, then: if q1 q2 q3 ∩ int(Ω) = ∅, modify Ω1 by deleting q2 from between q1 and q3 9 ; else let n(q1 ) = n(q2 ) = 1, ptr = q3 ; (9) case 4: q3 = q1 → prev, then: if q1 q2 q3 ∩ int(Ω) = ∅, modify Ω1 by deleting q1 from between q3 and q2 ,10 with letting ptr = q3 ; else let n(q3 ) = n(q1 ) = 1, ptr = q2 ; (10) end; (the end of while loop) (11) return Ω1 . Fig. 8 demonstrates the initial several computing steps of meshclosure showed in Fig. 3 using the mesh-closure-contracting algorithm. Now we are going to prove its correctness based on the following lemma. Lemma 3.2.1. Let M be a triangular mesh and Ω, Ω1 be two simple polygons such that Ω1 ⊇ Ω , Edge(Ω1 ) ⊆ Edge( M ) and Ω ⊆ ∂ M (also 6 See Fig. 4 where the given simple polygon Ω is simplified into an ellipse and the triangular zoning of Ω1 is hidden. 7 See Fig. 4. 8 See Fig. 5. 9 See Fig. 6. 10 See Fig. 7.

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

simple). Then we can have that:

¯ M Ω1 = Ω|

⇔

∀q1 q2 ∈ Edge(Ω1 ),

q1 q2 q3 ∩ int(Ω) = ∅,

where q2 = q1 → next (in counterclockwise order within Vertex(Ω1 )) and q3 is the one described in the mesh-closure-contracting algorithm. Proof. With the given condition, we will prove the following equivalent proposition instead:

¯ M Ω1 = Ω|

⇔

∃q1 q2 ∈ Edge(Ω1 ),

q1 q2 q3 ∩ int(Ω) = ∅,

where q2 , q3 are as above. (⇐) If there exists a triangle q1 q2 q3 such that q1 q2 ∈ Edge(Ω1 ), q1 q2 q3 ∩ int(Ω) = ∅, then we can always get an-

1365

other simple polygon Ω2 such that Ω1  Ω2 ⊇ Ω and Edge(Ω2 ) ⊆ Edge( M ), using the method described in the cases 1–4 of the mesh-closure-contracting algorithm. Therefore Ω1 is not the smallest one that covers Ω which leads ¯ M. to Ω1 = Ω| ¯ M , we have that Ω1  Ω| ¯ M. (⇒) Since Ω1 ⊇ Ω and Ω1 = Ω| ¯ M )) ⊆ Edge( M ), which means Recall that Edge(Ω1 ) (or Edge(Ω| ¯ M can “inherit” the triangular zoning of M that both Ω1 and Ω| (M is properly defined). Hence there exists at least one triangle  ¯ M ). Actually, let M (Π) = { ∈ M: of M such that  ⊆ Ω1 \ int(Ω|  ⊆ Π }, where Π is the union of some triangles of M, such as ¯ M and Ω1 \ int(Ω| ¯ M ). Obviously M ⊇ M (Ω1 )  M (Ω| ¯ M ) and Ω1 , Ω| ¯ ¯ M (Ω1 \ int(Ω| M )) = M (Ω1 ) \ M (Ω| M ) = ∅.

  ¯ M) , ∀ ∈ M Ω1 \ int(Ω|

¯ M ) = ∅.  ∩ int(Ω|

¯ M , we have int(Ω) ⊆ int(Ω| ¯ M ) which results in Since Ω ⊆ Ω|

  ¯ M) , ∀ ∈ M Ω1 \ int(Ω|

Fig. 6. A special case in the mesh-closure-contracting algorithm where q3 = q2 → next and q1 q2 q3 ∩ int(Ω) = ∅.

Fig. 7. A special case in the mesh-closure-contracting algorithm where q3 = q1 → prev and q1 q2 q3 ∩ int(Ω) = ∅.

 ∩ int(Ω) = ∅.

¯ M )), ∃l ∈ So we only need to prove that: ∃ ∈ M (Ω1 \ int(Ω| Edge(), such that l ∈ Edge(Ω1 ). ¯ M )), ∀l ∈ Edge(), l ∈ / Edge(Ω1 ), Suppose not: ∀ ∈ M (Ω1 \ int(Ω| ¯ M )), ∀l ∈ Edge(), ∃ ∈ M (Ω1 ), such that is, ∀ ∈ M (Ω1 \ int(Ω| that  =  and l ∈ Edge( ). ¯ M )), and let Ω2 = Take one fixed triangle 1 ∈ M (Ω1 \ int(Ω| ¯ M )). 1 . So Ω2 is a simple polygon and M (Ω2 ) ⊆ M (Ω1 \ int(Ω| ¯ M ∪ Ω2 )), such that l ∈ (∗) If ∃l ∈ Edge(Ω2 ), ∃ ∈ M (Ω1 \ int(Ω| Edge(), then let Ω2 take the value of Ω2 ∪ , while we still ¯ M )). have that Ω2 is simple and M (Ω2 ) ⊆ M (Ω1 \ int(Ω| Repeat the process (∗), and due to the finiteness of the set ¯ M )), the above process will definitely terminate. M (Ω1 \ int(Ω| ¯ M∪ When it terminates, we will have that either M (Ω1 \ int(Ω| ¯ M ∪ Ω2 )) ( = ∅), Ω2 )) = ∅ or ∀l ∈ Edge(Ω2 ), ∀ ∈ M (Ω1 \ int(Ω| / Edge(). l∈ Then, ∀l ∈ Edge(Ω2 ), ∃1 ∈ M (Ω2 ), such that l ∈ Edge( ). ¯ M )), according to the assumpSince  ∈ M (Ω2 ) ⊆ M (Ω1 \ int(Ω| tion we can get that ∃ ∈ M (Ω1 ), such that  =  and l ∈ Edge( ).

Fig. 8. The computing process of mesh-closure showed in Fig. 3 using the mesh-closure-contracting algorithm.

1366

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

¯ M ) using the conclusion which is established Then,  ∈ M (Ω| as the process (∗) terminates. ¯ M , ∀l ∈ Edge(Ω2 ). So l ⊆  ⊆ Ω| ¯ M are simple. ¯ M , because both Ω2 and Ω| So 1 ⊆ Ω2 ⊆ Ω| ¯ M ) which contradicts with 1 ∈ M (Ω1 \ int(Ω| ¯ M )). So 1 ∈ M (Ω| Hence the assumption is incorrect. Therefore ∃q1 q2 ∈ Edge(Ω1 ), q1 q2 q3 ∩ int(Ω) = ∅. 2 Theorem 3.2.2. The algorithm MESH-CLOSURE-CONTRACTING is correct, that is, for any input described in the mesh-closure problem, the output is exactly what Definition 3.2.2 defines, or symbolically ¯ M. MESH-CLOSURE-CONTRACTING( M , Ω) = Ω| Proof. The proof comprises two parts: Part 1. The “while” loop in the algorithm will definitely terminate. Part 2. The return value of the algorithm is exactly the meshclosure when the while loop terminates. Part 1. Let nt =  M and k be the number of the execution of the while loop. (∗∗) If q1 q2 q3 ∩ int(Ω) = ∅, where q1 = ptr, q2 = q1 → next (in counterclockwise order within Vertex(Ω1 )) and q3 is the one described in the algorithm, then q1 q2 q3 will be deleted from M (Ω1 ), whichever the case is. Then M (Ω1 ) will have at least one triangular element removed (in case 1, at least two triangles will be eliminated). Whereas  M (Ω1 )  nt , the process (∗∗) will definitely terminate. When it terminates, we will have that either M (Ω1 ) = ∅ or ptr(ptr → next)q3 ∩ int(Ω) = ∅. If M (Ω1 ) = ∅, then ∀ ∈ M \ M (Ω1 ) = M,  ∩ int(Ω) = ∅, which is impossible. Then the situation must be the second circumstance. Let ki be the number of the execution of the while loop when it happens the i-th time that ptr(ptr → next)q3 ∩ int(Ω) = ∅ (the addressing of this case is also counted in ki , so ki  1). For i = 1, k1  ( M \ M (Ω1 ))+ 1  ( M \ M (Ω1 ))+ M (Ω1 ) = nt , where the Ω1 is the new one after k1 loops. After (k1 − 1) loops, ptr(ptr → next)q3 ∩ int(Ω) = ∅, let a1 = ptr and after the successive loop, the pointer ptr would be a1 → next (or (a1 → next) → next which corresponds to the case 3 in the algorithm), as well as n(a1 ) = 1 (or n(a1 ) = n(a1 → next) = 1 or n(a1 ) = n(a1 → prev) = 111 ). For i = 2, after another successive (k2 − k1 − 1) loops, ptr(ptr → next)q3 ∩ int(Ω) = ∅, and during these (k2 − k1 − 1) loops, the pointer ptr remains the same because the situation can only be the cases 1–3 in the algorithm (if it is the case 1, then Ω  = Ωi ). Let a2 = ptr and after the successive loop, the pointer ptr would be a2 → next (or (a2 → next) → next), as well as n(a2 ) = 1. It also holds that k2  ( M \ M (Ω1 )) + 2  ( M \ M (Ω1 )) +  M (Ω1 ) = nt , where the Ω1 is the new one after k2 loops. Investigating the point sequence S = {a1 , a1 → next, . . . , ptr → prev}, where ptr is the new one after k2 loops, we have that ∀q ∈ S , n(q) = 1. Proceeding the argument about i as above, we can get the point sequence S = {a1 , a1 → next, . . . , ptr → prev} where ptr is the new one after ki loops, and we have that ∀q ∈ S , n(q) = 1. Since Vertex( M ) is a finite set, it definitely occurs that the sequence S  = {a1 , a1 → next, . . . , ptr → prev, ptr} (the one after ki loops, S ⊆ S  ) has some repeated points as i increases. Consider the first time it happens, at which we have that {a1 , a1 →

11

See cases 3–4 in the mesh-closure-contracting algorithm.

next, . . . , ptr → prev → prev} does not have any repeated points and that it does when ptr → prev, ptr are added into. If it does when ptr → prev is added into, then ptr → prev = a1 because Ω1 is simple. Then for Vertex(Ω1 ) = {a1 , a1 → next, . . . , ptr → prev → prev} ⊆ S , it holds that ∀q ∈ Vertex(Ω1 ), n(q) = 1, which terminates the while loop. Otherwise, it does when ptr is added into, then ptr = a1 . It also holds that ∀q ∈ Vertex(Ω1 ) = {a1 , . . . , ptr → pev} = S , n(q) = 1, which terminates the while loop. Therefore the “while” loop in the algorithm will definitely terminate. It also holds by mathematical induction that k  nt . Part 2. When the while loop terminates, it holds that the return value Ω1 is a simple polygon of M that contains Ω by mathematical induction of the two obvious points as below: (i) before the first while loop, Ω1 is a simple polygon of M that contains Ω ; (ii) and if before each while loop, Ω1 is a simple polygon of M that contains Ω , then after each while loop, so is Ω1 . Yet we can tell that ∀q ∈ Vertex(Ω1 ), n(q) = 1 ⇔ ∀q1 q2 ∈ Edge(Ω1 ), q1 q2 q3 ∩ int(Ω) = ∅, where q2 = q1 → next (in counterclockwise order within Vertex(Ω1 )) and q3 is the one described in the meshclosure-contracting algorithm. Applying Lemma 3.2.1 to the above conclusions, it holds that ¯ M when the while loop terminates. 2 Ω1 = Ω| 3.2.1. Technicalities Although the best data structure for implementing the above algorithm is doubly linked circular list equipped with pointers and objects, it can be executed with several one-dimensional arrays as we did in our program. Besides the mesh-closure-contracting algorithm for determining the mesh-closure of a given simple polygon within a triangular mesh, there can be two more algorithms: the expanding method in which a selected bad triangular element is expanded gradually to approach the mesh-closure, and the bounding method in which the edges of the mesh-closure are first computed and then the vertices are picked up one by one to form the mesh-closure. All the three methods can be proved correct and have almost linear asymptotic time complexity. So the mesh-closure-contracting algorithm is enough for the application. Also we can similarly consider the mesh-interior of a given simple polygon within a triangular mesh. 3.3. The mesh-union of the bad triangles within a triangular mesh This subsection introduces the method for the specifying of worse areas using the CFL-condition and shape-size metric criterion. Although the criteria can be various as we said before, our mesh-union method can yet be suitable for any criterion because only the labels of triangles resulting from the criterion matter. Definition 3.3.1. The triangular metric related to a triangle is 1 if the triangle is good according to some criterion and 0 if the triangle is bad. Let trim( p 1 p 2 p 3 ) be the triangular metric related to  p 1 p 2 p 3 . Definition 3.3.2. A finite polygonal set Λ = {Ωi } is the mesh-union j j j of bad triangles’ set T = { p 1 p 2 p 3 } (or worse areas) within a triangular mesh M (T ⊆ M) if it satisfies:

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

Ω is a simple polygon, ∀Ω ∈ Λ; Edge(Ω) ⊆ Edge

( M ), ∀Ω ∈ Λ;  p 1 p 2 p 3 ⊆ Ω∈Λ Ω , ∀ p 1 p 2 p 3 ∈ T ; Either Ω1 ∩ Ω2 = ∅, or Ω1 ∩ Ω2 has only one element which is

a common vertex of Ω1 and Ω2 , ∀Ω1 = Ω2 ∈ Λ; (5) Ω∈Λ Ω is the minimum (in the sense of inclusion or size) among the unions of polygons in Λ meeting the above four conditions.

(1) (2) (3) (4)

The mesh-union exists uniquely, denoted by Λ T | M . Fig. 9 provides an example for mesh-union. Note the difference between the j j j mesh-union and the union

of bad triangles’ set T = { p 1 p 2 p 3 }, whereas the latter Λ T is  p 1 p 2 p 3 ∈ T  p 1 p 2 p 3 . The algorithm for computing the mesh-union of bad triangles within a triangular mesh is proposed in pseudo-code. Mesh-union problem Input: A triangular mesh M and trim( p 1 p 2 p 3 ), ∀ p 1 p 2 p 3 ∈ M; Output: The mesh-union Λ T | M where T = { p 1 p 2 p 3 ∈ M: trim( p 1 p 2 p 3 ) = 0}. Before the introduction of the MESH-UNION-CONTRACTING algorithm a new kind of data structure is illustrated and termed the group of doubly linked circular lists. See Fig. 10 for an illustration. MESH-UNION-CONTRACTING( M , trim) (1) Let Ω1 = ∂ M (in counterclockwise order) and q1 be the Y – X ordering minimum of Vertex(Ω1 ); (2) build a group node gr pointing to Ω1 and let head, tail and ptr point to gr, ptr → lptr = q1 and ptr → next = NULL;

1367

(3) let n(q) = 0, ∀q ∈ Vertex( M ); (4) while (ptr = NULL) (5) do if n(q) = 1, ∀q ∈ Vertex(Ω(ptr → lptr)) (Ω(ptr → lptr) is the simple polygon in the group that the pointer ptr → lptr points to), then ptr = ptr → next; (6) else q1 = ptr → lptr, q2 = q1 → next (in counterclockwise order within Vertex(Ω(q1 ))), and q3 takes the unique point q such that q ∈ Vertex( M ), Edge(q1 q2 q) ⊆ Edge( M ) and q1 q2 q ⊆ Ω(q1 ); (there are four circumstances discussed below) (7) case 1: if q3 = q1 → prev and q3 = q2 → next and q3 ∈ Vertex(Ω(q1 )), then: if trim(q1 q2 q3 ) = 1, then Ω(q1 ) can be partitioned into two sub-polygons Ωi , Ωii besides q1 q2 q3 , where Ωi = [q1 , q3 , . . . , q1 → prev], Ωii = [q3 , q2 , q2 → next, . . .], then: when both Ωi , Ωii contain triangles with 0 triangular metric, let Ω(q1 ) = Ωi and build another group node gr pointing to Ωii , with letting sequentially tail → next = gr, tail point to gr, and tail → next = NULL; otherwise only one of them contains triangles with 0 triangular metric denoted by Ω  , therefore let Ω(q1 ) = Ω  , and ptr → lptr = q1 , if Ω  = Ωi , or ptr → lptr = q2 , if Ω  = Ωii ; else (which means trim(q1 q2 q3 ) = 0) let n(q1 ) = 1, ptr → lptr = q2 ; / Vertex(Ω(q1 )), then: (8) case 2: q3 ∈ if trim(q1 q2 q3 ) = 1, modify Ω(q1 ) by inserting q3 into between q1 and q2 ; else let n(q1 ) = 1, ptr → lptr = q2 ; (9) case 3: q3 = q2 → next, then: if trim(q1 q2 q3 ) = 1, modify Ω(q1 ) by deleting q2 from between q1 and q3 ; else let n(q1 ) = n(q2 ) = 1, ptr → lptr = q3 ; (10) case 4: q3 = q1 → prev, then: if trim(q1 q2 q3 ) = 1, modify Ω(q1 ) by deleting q1 from between q3 and q2 , with letting ptr → lptr = q3 ; else let n(q3 ) = n(q1 ) = 1, ptr → lptr = q2 ; (11) end; (the end of while loop) (12) return head. Theorem 3.3.1. The algorithm MESH-UNION-CONTRACTING is correct, that is, for any input described in the mesh-union problem, the output is exactly what Definition 3.3.2 defines, or symbolically MESH-UNION-CONTRACTING( M , trim) = Λ T | M .

Fig. 9. An example for the mesh-union.

The proof is completely similar to the one for mesh-closurecontracting algorithm. The only difference is that there is only one simple polygon to be dealt with during and after the while

Fig. 10. An illustration for the group of doubly linked circular lists.

1368

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

loop for mesh-closure-contracting algorithm, whereas there can be several simple polygons for mesh-union-contracting algorithm. We will not repeat it here, but it should be pointed out that: given a simple polygon Ω and a triangular mesh M such that Ω ⊆ ∂ M and ∂ M is also simple, if the triangular metric is taken as trim( p 1 p 2 p 3 ) = 1, if  p 1 p 2 p 3 ∩ int(Ω) = ∅, and ¯ M = ΛT |M , trim( p 1 p 2 p 3 ) = 0 otherwise, ∀ p 1 p 2 p 3 ∈ M, then Ω| where T = { p 1 p 2 p 3 ∈ M: trim( p 1 p 2 p 3 ) = 0}. See Fig. 3 for an interpretation. 3.3.1. Technicalities Although the best data structure for implementing the above algorithm is the group of doubly linked circular lists, it can be executed with several one-dimensional arrays as we did in our program. Besides the mesh-union-contracting algorithm for determining the mesh-union of bad triangles within a triangular mesh, there can be two more algorithms: the expanding method in which a selected bad triangular element is expanded gradually to approach the mesh-union, and the bounding method in which the edges of the mesh-union are first computed and then the vertices are picked up one by one to form the mesh-union. All the three methods have almost linear asymptotic time complexity but the last two cannot be proved. Sometimes for the application purpose the mesh-union-contracting algorithm can be slightly modified to assure the number of vertices of any polygon among the union12 greater than or equal to some integer (e.g., 20), which can be proved using the loop invariant method. 4. The rezoning of worse areas After the specifying of worse areas with the polygonal appearances, the vertices of obtained polygons should be adjusted as follows: (1) if a line segment on the symmetric axis contains only vertices of the polygon and is also a maximal (i.e., the two external points next to the two end points of line segment are neither the vertices of the polygon), then only the end points are kept, the internal points are deleted, and the line segment is equally divided into n subsegments (n can be taken as length(segment)/¯l, where ¯l is the average length over the Edge( M )); (2) if there are edges of the polygon lying on the free boundary, then usually they are lengthened and so will be equally divided into n subsegments (n can be taken as the above one) or just remain unchanged. After the modifying of the vertices of obtained polygons, the worse areas can be rezoned as follows: (1) we can rezone a simple polygon just by calling the routine triangle.c designed by Jonathan R. Shewchuk [26]; (2) the adaptive mesh methods including R-adaptive and H -adaptive can also be adopted; (3) although Patrick M. Knupp’s RJM (Reference Jacobian Matrix) method [4] is mainly for the continuous rezoning, it is acceptable for the rezoning of the obtained polygonal areas. For convenience the routine triangle.c is directly called without adding vertices into the boundary of the input simple polygon in our program. After the rezoning of worse areas, they are united with the remaining areas to compose the new mesh. The algorithm

12

This kind of union can be termed weak mesh-union.

for determining the data of new mesh in the format described in the previous section is provided. New-mesh problem Input: A triangular mesh M = (x, y , tn) and some local triangular meshes LM i = (lxi , ly i , ltni ), 1  i  k, such that Edge( LM i ) ⊆ Edge( M ) (except the edges on the symmetric axes or the free boundary), ∀1  i  k, and LM i ∩ LM j = ∅, ∀i = j; Output: The new mesh N M = (nx, ny , ntn) composed by the triangles in LM i (∀1  i  k) and the triangles of M which stay outside ∂( LM i ) (∀1  i  k). NEW-MESH( M , LM i |ki=1 )13 (1) (2) (3) (4) (5) (6) (7)

(8) (9) (10) (11) (12) (13) (14)

Let (nx, ny , ntn) = (lx1 , ly 1 , ltn1 ); for i = 2 : 1 : k (if k = 1, then this loop can be omitted) for j = 1 : 1 : length(lxi ) if the point (lxi ( j ), ly i ( j )) is not in (nx, ny ), then nx = [nx, lxi ( j )], ny = [ny , ly i ( j )]; end; (the end of j loop) for j = 1 : 1 : (length(ltni )/3) find the indices of the three points (lxi (ltni (3 ∗ j − 2)), ly i (ltni (3 ∗ j − 2))), (lxi (ltni (3 ∗ j − 1)), ly i (ltni (3 ∗ j − 1))) and (lxi (ltni (3 ∗ j )), ly i (ltni (3 ∗ j ))) into (nx, ny ), say id1 , id2 , id3 , then ntn = [ntn, id1 , id2 , id3 ]; end; (the end of j loop) end; (the end of i loop) for s = 1 : 1 : length(x) if the point (x(s), y (s)) lies strictly outside ∂( LM i ) (∀1  i  k), then nx = [nx, x(s)], ny = [ny , y (s)]; end; (the end of s loop) for t = 1 : 1 : (length(tn)/3) if the point ( 13 (x(tn(3 ∗ t − 2)) + x(tn(3 ∗ t − 1)) + x(tn(3 ∗ t ))),

1 ( y (tn(3 ∗ t − 2)) + y (tn(3 ∗ t − 1)) + y (tn(3 ∗ t )))) lies strictly 3 outside ∂( LM i ) (∀1  i  k), then find the indices of the three points (x(tn(3 ∗ t − 2)), y (tn(3 ∗ t − 2))), (x(tn(3 ∗ t − 1)), y (tn(3 ∗ t − 1))) and (x(tn(3 ∗ t )), y (tn(3 ∗ t ))) into (nx, ny ), say id1 , id2 , id3 , then ntn = [ntn, id1 , id2 , id3 ]; (15) end; (the end of t loop) (16) return (nx, ny , ntn).

4.1. Technicalities In practice we need to obtain some other attributes about the new mesh such as node-flag and triangle-flag with which the remapping of physical variables can be performed locally only on the worse areas, which can be obtained by the slight modification of the NEW-MESH algorithm. In some unstructured mesh the ct data14 is added into the triangular mesh data structure [x, y , tn], while in some circumstances it is necessary to insert vertices into the boundary of the input simple polygon. The NEW-MESH algorithm with some simple adaptation is fulfilling for both cases. We will demonstrate some figures for the above methods as follows. All figures in this paper use CGS units (i.e., cm, g , s). Fig. 11 shows the initial Lagrangian mesh. Figs. 12 and 13 display the old mesh and new one at 1 ns respectively from which we can tell that the majority of the old mesh remains unchanged. Fig. 14 illustrates the worse area in Fig. 12 specified by the MESH-CLOSURE-CONTRACTING algorithm and Fig. 15 is the local new mesh after the rezoning of worse area in Fig. 14 by calling triangle.c. 13 In the following parts i : j : k means the array [i , i + j , i + 2 j , i + 3 j , . . . , k] where k = i + nj. 14 The ct is short for cell type, i.e., the mesh triangular element type.

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

1369

Fig. 11. The Lagrangian mesh at initial time.

Fig. 12. The heavily distorted mesh at 1 ns.

Fig. 13. The new mesh after the rezoning of old one in Fig. 12.

The triangle-node connectivities in Figs. 12 and 13 are locally different15 due to the fact that triangle.c cannot maintain the mesh connectivity. If we want to keep the mesh connectivity, calling triangle.c will not be appropriate any more and other methods such as RJM should be considered. Since our rezoning method is occasional, it is difficult to have the same mesh connectivity before

15

This problem makes no trouble to our remapping methods yet.

and after rezoning, especially when the old mesh is heavily distorted just like the situation in Fig. 12. The triangular elements of Fig. 12 which are constrained within the simple polygon of Fig. 14 are continuously compressed by the others, and the corresponding mesh nodes become increasingly dense so that we have to delete some of them if we want to proceed the simulation. Certainly the number of mesh elements will be much larger in the practical simulation.

1370

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

Fig. 14. The polygon of worse area in Fig. 12.

Fig. 15. The local new mesh after the rezoning of worse area in Fig. 14.



5. The remapping of worse areas

1  ρ  j · V ( j )), where V (Ω) means the rotating volume of the planar region Ω , i.e., V (Ω) = 2π · S (Ω) · y c , where S (Ω)

V (i ))/(

The remapping algorithms of physical variables from the old mesh onto the new one have been well developed. The least standard is the conservation of physical variables.16 Also we care the efficiency and accuracy of remapping algorithms. In the following parts, all the quantities of the new mesh are distinguished with a  , the new x-node‘tilde’, e.g., the new mesh triangular element  ˜ the new y-node-velocity v˜ , etc. velocity u,

j

means the size of Ω and y c is the y-coordinate of the gravity center of Ω (e.g., y c () = 13 ( y 1 + y 2 + y 3 )). Here we consider the cylindrical coordinates (i.e., x stands for z, and y for r) not the Cartesian coordinates. Finally the new mass-density can be defined as 1 ρ  j = c ρ  .

(5)

j

5.1. The linear interpolation

Obviously

The linear interpolation interpolates the physical variables from the old mesh onto the new one linearly. The polynomial interpolation can be similarly discussed. Here we only study the linear one of the mass-density which is cell-centered and the node-velocity which is node-centered.  j , find the triangle Let ρ (x, y ) = ρi , ∀(x, y ) ∈ i . Then for   j is located, say of the old mesh in which the gravity center of  i , and let 1 ρ  = ρi . j

16 The rigorous definition can be somehow various just as what we will show in the following subsections.

i

ρi · V (i ) =

 j

 j ), i.e., mtotal = m ˜ total . So ρ  j · V (

the total mass is conserved. Let u (x, y ) = ai · x + bi · y + c i , ∀(x, y ) ∈ i , where ai , bi , c i are constants and can be determined by the following three linear equations: u (x j , y j ) = u j , where (x j , y j ) are the three vertices of i and u j are the three x-node-velocities respectively, j = 1, 2, 3. Then for the new node (˜x, y˜ ), find the triangle of the old mesh in which the point (˜x, y˜ ) is located, say i , and let u  (˜x, y˜ ) = ai · x˜ + bi · y˜ + c i . The computing result of the above procedures is that

(4)

In order to maintain the conservation of total mass, we  need to get ρ 1 multiplied by some constant c, where c = ( i ρi ·



u  (˜q) =

u 1 S (q2 q3 q˜ ) + u 2 S (q1 q3 q˜ ) + u 3 S (q1 q2 q˜ ) S (q1 q2 q3 )

,

(6)

where q˜ = (˜x, y˜ ) and q1 q2 q3 = i . In order to conserve the total momentum, we need to get δ u [17], where δ u meets the u  added  constant  condi by some ˜ j (u j + δ u ), which leads to δ u = ( i mi u i − tion i mi u i = jm

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376







˜ j u j )/( jm

˜ j ), where mi , m ˜ j are the masses of the old mesh jm nodes and new ones respectively. So let u  (˜x, y˜ ) = u  (˜x, y˜ ) + δ u .

(7)

We also need to address u  as follows for the conservation of total kinetic energy:







¯ 

u˜ = u + a u − u , (8)     ˜ j u )/( j m ˜ j ), and the constant a holds for where u¯ = ( j m  1  1 2  + a(u  − u ¯  ))2 . So the new x-node-velocity ˜ m u = m ( u i 2 i i j 2 j u˜ (˜x, y˜ ) is determined. The new y-node-velocity v˜ (˜x, y˜ ) can be similarly calculated. The above algorithm can guarantee the conservation of total mass, total momentum and total kinetic energy, i.e., the conservation is global. Another kind of conservation of physical variables which is termed local conservation is more accurate and popular, and will be presented in the next subsection. 5.2. The particle remapping In this subsection we will discuss the particle remapping method. But at first the integral remapping is defined which corresponds to the local conservation of physical variables. Definition 5.2.1. Given two sets of triangular meshes M = {i }, ˜ = {  j } and the mass-density distribution ρ (x, y ) on M, the M mass-density remapping is integral if the average mass-density of  j satisfies: 

ρ j =

 i



ρ (x, y ) dV



 j ), V (

(9)

 j) V (i ∩

Again here we consider the cylindrical coordinates (i.e., x stands for z, and y for r) not the Cartesian coordinates. The integral remapping maintains the local conservation of the mass in the following sense:

˜ ˜  i =  ˜ 0  j (M 0 ⊆ M , M 0 ⊆ M), then  j ∈M   17 ˜ ρ (x, y ) dV = V (Ω) ρ (x, y ) dV . The proof is quite straightV (Ω)

i ∈ M 0

forward. Apparently the local conservation of mass implies the global conservation of mass. And since we do not know the actual mass-density distribution ρ (x, y ) on M in practice (we only know the average massdensity ρi ), we have to make some approximations about ρ (x, y ) in actual calculations. The integral remapping of node-velocity can be a little different because we don’t know whether the momentum or the kinetic energy should be conserved, e.g., in the above definition ρ (x, y ) will be substituted by ρ (x, y )u (x, y ) (corresponding to momentum) or 1 ρ (x, y )u 2 (x, y ) (corresponding to kinetic energy), the denomina2 tor should be the mass, and the integral regions should be the overlaps (intersection regions) of control volumes18 of the old and new mesh nodes. Usually the local conservation of momentum and the local conservation of kinetic energy are incompatible,19 and

17

 V (i )

18

ρ (x, y ) dV = ρi V (i ),



j) V (

 j ). ρ˜ (x, y ) dV = ρ  j V (

A polygonal region which surrounds the mesh node. E.g., for the case in which two triangles with the same mass and the same speed of opposite directions are merged, the new velocity should be zero according to the conservation of momentum but nonzero according to the conservation of kinetic energy. 19

ρ j = (

i

 j ))/ V (  j ). ρi V (i ∩ 

The particle remapping method [18] is developed in order to  j . The basic idea is to avoid calculating the overlaps of i and  decompose each old mesh element into many small triangles with the same shape, to let their gravity centers possess their masses and to compute the mass of each new mesh element by accumulating the mass of gravity center which lies inside the new mesh element.20 The new node-velocity can be similarly determined using momentum conservation by the decomposition of control volumes of old mesh nodes. The particle remapping maintains almost local conservation (not exactly). For application purpose, the coordinates of the gravity centers of sub-triangles are given using the so-called Gravity Center Matrix. Let qi = (xi , y i ), i = 1, 2, 3, each edge of q1 q2 q3 be equally divided into n sub-segments, and the coordinates of gravity centers of n2 sub-triangles of q1 q2 q3 be denoted by xc = (xc1 , xc2 , xc3 , . . . , xnc 2 )T , y c = ( y c1 , y c2 , y c3 , . . . , ync 2 )T . Then we can have

1

xc = yc =

3n 1 3n

· A · (x1 , x2 , x3 )T ,

(10)

· A · ( y 1 , y 2 , y 3 )T ,

(11)

where A = (ai j )n2 ×3 is the Gravity Center Matrix provided as follows. GRAVITY-CENTER-MATRIX(n)

where dV = 2π y dx d y.

if Ω =

the former is preferred while the loss of kinetic energy is added to the internal energy. It is the integral over the overlap that the new physical variable by integral remapping heavily depends on. Assuming ρi (x, y ) = ρi , ∀(x, y ) ∈ i , the new average mass-density of  j denoted by ρ  j , i.e.,   j relies on the overlaps of i and 





1371

(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)

if n = 1, then A = [1, 1, 1]; else A = []; (empty matrix) for i = 1 : 1 : (n − 1) a1 = ((1 : 1 : (n − i + 1)) ∗ 0 + 3 ∗ i − 2)T , a2 = (1 : 3 : (3 ∗ n − 3 ∗ i + 1))T , a3 = ((3 ∗ n − 3 ∗ i + 1) : (−3) : 1)T , A 1 = [a1 , a2 , a3 ]; a4 = ((1 : 1 : (n − i )) ∗ 0 + 3 ∗ i − 1)T , a5 = (2 : 3 : (3 ∗ n − 3 ∗ i − 1))T , a6 = ((3 ∗ n − 3 ∗ i − 1) : (−3) : 2)T , A 2 = [a4 , a5 , a6 ]; A 3 = [ A 1 ; A 2 ], A = [ A ; A 3 ]; ( A 3 = [ A 1 ; A 2 ] means that A 2 constitutes other rows under A 1 ) end; (the end of i loop) A = [ A ; 3 ∗ n − 2, 1, 1]; return A.

E.g., for n = 3 (see Fig. 16),

⎡

1 ⎢1 ⎢ ⎢1 ⎢ ⎢2 ⎢ A=⎢2 ⎢ ⎢4 ⎢4 ⎢ ⎣5 7

1 4 7 2 5 1 4 2 1

⎤

7 4⎥ ⎥ 1⎥ ⎥ 5⎥ ⎥ 2 ⎥. ⎥ 4⎥ ⎥ 1⎥ 2⎦ 1

There are some interesting properties about the gravity center matrix, e.g., the sum of the elements of any row of A is 3n, and for any row ai = (ai1 ai2 ai3 ) of A, the permutation of ai is still a row of A. 20

The boundary of it doesn’t matter.

1372

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

Fig. 16. The elements of gravity center matrix for n = 3.



Finally it can be proved that limn→+∞ ρ  j (n) = ρ j = (

i

ρ i ×

 j ))/ V (  j ), where n means that each edge of the old V (i ∩  mesh elements is equally divided into n sub-segments. 5.3. The first order integral remapping

The reason why the integral remapping is more accurate than the linear interpolation is that the latter uses only the position information (i.e., x– y-coordinates) of the mesh nodes other than the triangle-node connectivity (i.e., tn), which exactly represents the essence of triangular mesh and is effectively employed by the former. The first order integral remapping of average mass-density computes the new one directly as follows:

ρ j =



   j ), ρi V (i ∩  j ) V (

(12)

i

assuming that ρ (x, y ) = constant = ρi , ∀(x, y ) ∈ i . The crucial algorithm is the OVERLAP algorithm which determines the overlap of any two convex polygons. The doubly linked circular lists are used in our algorithm OVERLAP(Ω1 , Ω2 ) where Ω1 , Ω2 and the return are all convex polygons. We refer to our programs for more detailed information. The first order integral remapping of node-velocity maintains the local conservation of momentum as follows:

u˜ j =



m( I i ∩ ˜I j )u i

i

v˜ j =



  i

m( I i ∩ ˜I j ) v i

 

i

m( I i ∩ ˜I j ) ,

 ˜ m( I i ∩ I j ) ,

(14)

Fig. 17 displays, then the new node-velocity would be

u˜ j =

m( I ik ∩ ˜I j )u ik

1k3

v˜ j =



1k3

m( I ik ∩ ˜I j ) v ik

 

1k3



i

i 2

1

i1

ρ (x, y ) dx d y = ρi2 S (i2 ). The above three equations are lin-

ear and the coefficients can be determined by using the fact that

 x dx d y =

1 3

(x1 + x2 + x3 ) S (q1 q2 q3 ),

(16)

q1 q2 q3



y dx d y =

1 3

( y 1 + y 2 + y 3 ) S (q1 q2 q3 ).

(17)

q1 q2 q3

i

 

where ai , bi , c i are constants and can be determined by the average mass-densities of i and two of its adjacent triangles such that (|ai | + |bi |) is the least one, say i 1 , i 2 , i.e.,  ρ (x, y ) dx d y = ρi S (i ),  ρ (x, y ) dx d y = ρi S (i 1 ), and

(13)

where I i , ˜I j are the control volumes of the old mesh nodes and new ones respectively. And the I i ∩ ˜I j is calculated by decomposing them into convex polygons and calling OVERLAP(Ω1 , Ω2 ). If the control volume ˜I j of the new mesh node j is as what



Fig. 17. An illustration for the control volumes where I ik (k = 1, 2, 3) and ˜I j are the control volumes of the old mesh nodes ik (k = 1, 2, 3) and new one j respectively.

After 

the

computing

of

a i ,

b i ,

and

c i

the

integral

  j ρ (x, y ) dx d y is determined by decomposing i ∩  j into i ∩

triangles, then the new average mass-density can be specified as

ρ j =

  i

ρ (x, y ) dx d y



 j ). S (

(18)

j i ∩

m( I ik ∩ ˜I j ) ,

As for the second order ENO remapping of node-velocity it is assumed that

˜ m( I i k ∩ I j ) .

P x (x, y ) = a I i · x + b I i · y + c I i ,

(19)

where P x = ρ (x, y )u (x, y ) is the x-momentum-density and the I i is the control volume of old mesh nodes. The constants a I i , b I i and c I i can be determined by the x-node-velocity of I i and two of its  neighbors such that (|a I i | + |b I i |) is the least one, say I i 1 , I i 2 , i.e., I P x (x, y ) dx d y = mI u ( I i ), I P x (x, y ) dx d y = mI u ( I i 1 ), and

1k3

5.4. The high order ENO remapping The high order ENO remapping can be regarded as the natural generalization of the first order integral remapping [19,20]. And in this subsection we consider the Cartesian coordinates. The second order ENO remapping assumes that

ρ (x, y ) = ai · x + bi · y + ci , ∀(x, y ) ∈ i ,

∀(x, y ) ∈ I i ,

(15)



i

i

I i2

i1

i1

P x (x, y ) dx d y = mI u ( I i 2 ). The new x-momentum-density can i2

be computed as

P x ( ˜I j ) =

  i

I i ∩ ˜I j

 P x (x, y ) dx d y

S ( ˜I j ).

(20)

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

grams are called to perform the rezoning and remapping phases. So the whole process of ICF can be quantitatively understood now.

Table 1 Efficiency and accuracy. Algorithm

Efficiency

Accuracy

Linear interpolation Particle remapping First order remapping k-th order ENO remapping

O (nt · n˜ t ) O (n2 · nt · n˜ t ) O (nt · n˜ t ) O (nt · n˜ t )

0 1 1 k

Then the linear functions ρ˜ (x, y ) and P˜ x (x, y ) can be reconstructed based on the ρ  j and P x ( ˜I k ) of the new mesh just as the linear functions

ρ (x, y ) and P x (x, y ) of the old mesh. We simply define:

u (˜x j , y˜ j ) = P˜ x j , y˜ j )/ρ¯ (˜x j , y˜ j ), x (˜

ρ¯ (˜x j , y˜ j ) =

1373

x j , y˜ j ) k (˜ k ) ρ˜ (˜x j , y˜ j )∈Vertex(

k : (˜x j , y˜ j ) ∈ Vertex( k )} {

(21)

.

(22)

The reason why we assume that the x-momentum-density P x (x, y ) is linear rather than the x-node-velocity u (x, y ) is that the inte gral function ρ (x, y )u (x, y ) in I ∩ ˜I ρ (x, y )u (x, y ) dx d y will be a i

j

quadratic function for which the order of accuracy will be three when the u (x, y ) is linear. The ENO remapping method with higher order of accuracy can be similarly discussed and all high order ENO remapping methods maintain the local conservation of physical variables. 5.5. The efficiency and accuracy of remapping algorithms This subsection provides the preliminary estimation of efficiency and accuracy of above four remapping methods. The precise analysis can be extremely complicated, especially the asymptotic time complexity analysis. We first introduce the notations which will appear. In the following table, nt and n˜ t stand for the numbers of old mesh triangles and new ones respectively, n means that each edge of the old mesh elements is equally divided into n sub-segments, and (k − 1) represents the order of the polynomial arising in the k-th order ENO remapping method. Table 1 contains the efficiency and accuracy of above four remapping methods. 5.5.1. Technicalities Although two algorithms can have the same asymptotic time complexity, the practical performances of them may be of great differences (e.g., given f 1 (n) = n2 , f 2 (n) = 100n2 , then f 1 (n) = f 2 (n) = O (n2 ) but f 2 is one hundred times of f 1 ). We refer to [19,20] for detailed information on accuracy test. In the following parts we will focus on the accuracy of the first order integral remapping. 6. The implementation of our algorithms This section discusses the issues on the implementation of our algorithms. Since our program is mainly designed for ICF simulation, we will start from introducing some background on ICF and MULTI2D. Indirectly irradiated targets for Inertial Confinement Fusion (ICF) rely on radiative energy transport to induce the implosion of thermonuclear fuel to ignition conditions [27–30]. Intense laser or ion beam power is converted into thermal radiation and is confined in a cavity (hohlraum) having high-Z walls. The almost isotropic radiation field produced in this way uniformly heats the external layers of a spherical capsule, filled with a deuterium– tritium mixture, and induces its hydrodynamic implosion to thermonuclear conditions. MULTI2D [31] created by R. Ramis is used to fulfill the Lagrangian phase of the simulation of the ICF process, while our pro-

Our programs comprise the following components: (1) the routines which import and export data composed by the mesh structure and physical variables, e.g., m2dppinit( ), filewrite( ); m2dppinit( ) is used to provide the input data to Matlab from the computing results of MULTI2D, and filewrite( ) is called to write a file consisting of the rezoning and remapping results; it is required for non-MULTI2D users that these routines are reprogrammed; (2) the routines which calculate the CFL-condition and shape-size metric, e.g., cfl_time( ), f _ref _shape( ); (3) the routines which specify the mesh-closure and mesh-union, e.g., closure(), weak_mesh_union( ); (4) the routines which rezone the worse areas, e.g., localrezone( ); (5) the routines which perform the remapping of physical variables, e.g., remap( ); we only remap the primary physical variables of MULTI2D as follows:

[nmeshwa, mass] = remap(meshwa, lx, ly, ltn); remap( ) uses the first order integral remapping method while the other remapping methods are not included for the moment; remap( ) can be slightly modified to remap the nonMULTI2D primary physical variables; (6) the subroutines which are called by the above ones, e.g., cellremap( ), noderemap( ), overlap( ).21 An example is offered here: (1) ggg = m2dppinit( ); %import data from MULTI2D (2) time = 1.07378e–09; (3) meshwa = initial(ggg, time); %meshwa is a structure constituted by the MULTI2D primary physical variables (4) dt = cfl_time(meshwa. x, meshwa. y, meshwa. tn, meshwa. rho, meshwa. e, p, mid); %p and mid are included in ggg (5) vida1 = closure(ggg, 0); (6) vida2 = closure(ggg, time); (7) vida = [vida1(1), vida2(142 : 178), vida1(367 : end)]; %modify vida1 and vida2 by their figures to obtain vida (8) [lx, ly, ltn] = localrezone([meshwa.x(vida)’, meshwa.y(vida) ]); %the angle and size of new mesh element can be controlled (9) [nmeshwa, mass] = remap(meshwa, lx, ly, ltn); %nmeshwa means the new mesh with argument (10) filewrite(‘. . /1’, time, nmeshwa, mass). 7. The application in ICF simulation This section introduces the application of our method in the ICF simulation. We will also compare our results to the ones from MULTI2D. We will display some figures of the ICF simulation as follows in which the different colors represent different substances: red for Au, green for C–H and blue for D–T. Fig. 18 shows the ICF Lagrangian mesh with 16 433 nodes and 32 263 triangles at initial time. Then MULTI2D is called to simulate the Lagrangian motion of the mesh until (3.08e–10) s when the mesh is heavily distorted as in Fig. 19 and is rezoned and remapped22 as in Fig. 20. After that the simulation continues with several other rezonings and remappings until (8e–10) s when the spherical capsule filled

21 The usage of our programs can be just in the above order and we refer to our programs for more detailed comments. 22 Here we use the mesh-closure-contracting algorithm and the first order integral remapping.

1374

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

Fig. 18. The ICF Lagrangian mesh with 16 433 nodes and 32 263 triangles at initial time. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 19. The old ICF Lagrangian mesh at (3.08e–10) s. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 20. The new ICF Lagrangian mesh at (3.08e–10) s. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 21. The ICF Lagrangian mesh at (1.07e–9) s. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

with a deuterium–tritium mixture is about to be fast compressed. Fig. 21 depicts the ICF Lagrangian meshes at (1.07e–9) s and then the simulation halts. Fig. 22 zooms in Fig. 21 to expose the symmetry property23 of the capsule while it is awful.

In order to investigate the accuracy of the first order integral remapping, we can do some ideal tests (i.e., assuming that we know the actual mass-density distribution ρ (x, y ) on M). We can find such tests in [20] which also show that the first order integral remapping is indeed first order accurate.24 The ideal tests for

23 The symmetry property is extremely important in ICF research and we will not focus on it in this paper.

24 In [20] it is termed CIB/DC which has the same order of accuracy as SFB/DC or “da”.

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

high order ENO remapping can also be found in [19]. Instead of these tests we compare the results with remapping performed to the ones without remapping performed to show the accuracy of the first order integral remapping. Fig. 23 displays the x-node-velocity distribution u (x, y ) at (5e–10) s without remapping performed. The whole process only uses MULTI2D which took about four days to advance to (5e–10) s on my computer, and at (5e–10) s the computing time step is about (1e–18) s due to the great distortion of Lagrangian mesh. Fig. 24 illustrates the x-node-velocity distribution u˜ (x, y ) at (5e–10) s with 4 remapping phases (the first order integral remapping) performed at (3.08e–10), (3.96e–10), (4.34e–10) and (4.75e–10) s respectively. The whole process uses MULTI2D to per-

1375

form the Lagrangian phase and totally cost my computer about 4 hours to advance to (5e–10) s. We can tell from Figs. 23 and 24 that the remapping causes little change over the simulation of MULTI2D, while considerably improves the efficiency. The only deficiency is that the details about the heavily distorted area were kind of lost. 8. Concluding remarks and future work The local rezoning and remapping method for unstructured mesh is explicitly proposed in this paper including the criteria for worse areas, the specifying, the rezoning, and the remapping of worse areas. Four remapping methods are discussed and compared. The application in ICF simulation suggests the success of our methods. The future work will be the remapping error estimation and control. References

Fig. 22. The capsule within the ICF Lagrangian mesh at (1.07e–9) s. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

[1] C.W. Hirt, A.A. Amsden, J. Cook, J. Comput. Phys. 14 (1974) 227; reprinted in C.W. Hirt, A.A. Amsden, J. Cook, J. Comput. Phys. 135 (1997) 203. [2] P. Knupp, SIAM J. Sci. Comput. 23 (1) (2001) 193. [3] P. Knupp, Finite Elem. Anal. Des. 39 (2003) 217. [4] P. Knupp, L.G. Margolin, M. Shashkov, J. Comput. Phys. 176 (2002) 93. [5] D.J. Benson, Comput. Methods Appl. Mech. Engrg. 72 (1989) 305. [6] S. Giuliani, Nucl. Eng. Des. 72 (1982) 205. [7] J.S. Peery, D.E. Carroll, Comput. Methods Appl. Mech. Engrg. 187 (2000) 591. [8] M. Ainsworth, J.T. Oden, Comput. Methods Appl. Mech. Engrg. 142 (1997) 1. [9] J.U. Brackbill, J.S. Saltzman, J. Comput. Phys. 46 (1982) 342. [10] J.K. Dukowicz, J. Comput. Phys. 56 (1984) 324. [11] A. Anderson, X.M. Zheng, V. Cristini, J. Comput. Phys. 208 (2005) 616. [12] X.M. Zheng, J. Lowengrub, A. Anderson, V. Cristini, J. Comput. Phys. 208 (2005) 626. [13] V. Cristini, J. Blawzdziewicz, M. Loewenberg, J. Comput. Phys. 168 (2001) 445. [14] K.R. Moyle, Y. Ventikos, J. Comput. Phys. 227 (2008) 2781. [15] W.B. Castelló, F.G. Flores, Comput. Methods Appl. Mech. Engrg. 198 (2008) 332. [16] Q.D. Cai, H.S. Shui, S.W. Fu Chin, J. Comput. Phys. 1 (2001) 17. [17] J. Wu, Z.X. Liang, Y.F. Li, Chin. Numer. Comput. Comput. Appl. 4 (1994) 261. [18] H.G. Horak, E.M. Jones, J.W. Kodis, et al., J. Comput. Phys. 26 (1977) 277.

Fig. 23. The x-node-velocity distribution at (5e–10) s without remapping performed. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 24. The x-node-velocity distribution at (5e–10) s with several remapping performed. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

1376

Z. Lin et al. / Computer Physics Communications 182 (2011) 1361–1376

[19] J. Cheng, C.W. Shu, Appl. Numer. Math. 58 (2008) 1042. [20] L.G. Margolin, M. Shashkov, J. Comput. Phys. 184 (2003) 266. [21] L.G. Margolin, M. Shashkov, Comput. Methods Appl. Mech. Engrg. 193 (2004) 4139. [22] R. Loubère, M. Shashkov, J. Comput. Phys. 209 (2005) 105. [23] J. Grandy, J. Comput. Phys. 148 (1999) 433. [24] J.K. Dukowicz, J.W. Kodis, SIAM J. Sci. Stat. Comput. 8 (1987) 305.

[25] [26] [27] [28] [29] [30] [31]

L.G. Margolin, J. Comput. Phys. 135 (1997) 198. J.R. Shewchuk, Comput. Geom. 22 (2002) 21. F.P. Mattar, M.C. Newstein, Comput. Phys. Comm. 20 (1980) 139. F.L. Addessio, M. Cline, J.K. Dukowicz, Comput. Phys. Comm. 48 (1988) 65. R.A. Clark, Comput. Phys. Comm. 48 (1988) 61. S. Atzeni, A. Schiavi, F. Califano, et al., Comput. Phys. Comm. 169 (2005) 153. R. Ramis, J. Meyer-ter-Vehn, J. Ramírez, Comput. Phys. Comm. 180 (2009) 977.