# Robust Three-Dimensional Level-Set Method for Evolving Fronts on Complex Unstructured Meshes

Authors: Ran Wei; Futing Bao; Yang Liu; Weihua Hui
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

With a purpose to evolve the surfaces of complex geometries in their normal direction at arbitrarily defined velocities, we have developed a robust level-set approach which runs on three-dimensional unstructured meshes. The approach is built on the basis of an innovative spatial discretization and corresponding gradient-estimating approach. The numerical consistency of the estimating method is mathematically proven. A correction technology is utilized to improve accuracy near sharp geometric features. Validation tests show that the proposed approach is able to accurately handle geometries containing sharp features, computation regions having irregular shapes, discontinuous speed fields, and topological changes. Results of the test problems fit well with the reference results produced by analytical or other numerical methods and converge to reference results as the meshes refine. Compared to level-set method implementations on Cartesian meshes, the proposed approach makes it easier to describe jump boundary conditions and to perform coupling simulations.


---

## Body

## 1. Introduction

We are interested in evolving fronts on complex unstructured meshes in their normal direction using the level-set approach [1]. The level-set approach has been applied to various evolving-front problems. In particular, for evolving in the normal direction, the examples include image segmentation [2] and phase transition of solid materials, such as burning [3, 4] and exploding [5]. The level-set approach uses a higher-dimensional function ϕ to model the evolving front, where the x→∣ϕx→=0 isosurface denotes the front. Since the implicit representation automatically handles topological changes, the level-set approach is capable of processing arbitrary geometry and evolving speed.The level-set method is usually applied on uniform Cartesian structured meshes. Advantages of Cartesian meshes include that they are easy to generate and that numerical technologies to solve Hamilton-Jacobi equations on Cartesian meshes are mature. When applied to image segmentation, the Cartesian mesh is naturally composed of the pixels.However, when solid materials are involved in the simulation, there are cases in which non-Cartesian meshes outperform Cartesian ones, such as the following: (a) Describing complex boundary conditions. Jump conditions across the solid boundaries and different surface patches can be critical [6]. In this situation, tetrahedral or hexahedral meshes can be generated to fit the boundaries directly, while Cartesian meshes must be treated specially to handle the complex boundary conditions. One example of this issue is [7], where additional level-set functions are applied to address different boundaries, thus introducing more complexity. (b) Multi-phase coupling simulation. Stress, fluid, and thermal coupling simulations are often required in addition to evolving-front ones. Many such simulations are easier to perform on non-Cartesian meshes. If an evolving-front simulation is run on Cartesian meshes, interpolating among different meshes would be inevitable. (c) For “sparse” solid objects. Ordinary Cartesian meshes fill a rectangular bounding box of the evolving front. For “sparse” solid parts (e.g., the “L” or ring-shaped ones), a number of Cartesian cells, which are outside the object, would never be swept by the front of concern. However, storing, operating, or skipping these cells costs extra computational resources and consequently lowers the operation efficiency.It is clear from the above discussion that tetrahedral and hexahedral meshes outperform Cartesian ones in evolving fronts with complex boundary conditions or in coupling simulations. Using the level-set method to evolve fronts essentially requires solving the Hamilton-Jacobi equation formed like∂ϕ/∂t+H∇ϕ,x→=0. Previous efforts to solve such equations on unstructured meshes have been focused on two-dimensional (2D) triangular meshes. In [6, 8, 9], various numerical methods were proposed to solve the equation. In [10], the authors have studied handling boundary conditions for Hamilton–Jacobi equations on triangular meshes. There have also been studies on high-order numerical methods for the problem: In [11–13], WENO and Hermite WENO reconstruction methodologies were studied separately on 2D unstructured meshes; Augoula [14] developed several high-order numerical approximations for first-order Hamilton Jacobi equations on triangular meshes on the base of the Lax-Friedrichs scheme. Zhang [15] and Rokicki [16] then extended WENO reconstruction to 3D tetrahedral meshes.Efforts to extend the level-set method to 3D tetrahedral meshes include [17, 18]. Morgan [17] developed 3D level-set methods for evolving fronts on tetrahedral meshes with adaptive mesh refinement, which are suitable for both advecting the level-set field and for evolving a front in the surface normal direction. Fu [18] mainly worked on applying parallel computing technology to speed up level-set solving. We have noticed that both of the aforementioned studies estimate the gradient of ϕ within a tetrahedral cell on the basis of the divergence theorem, i.e., ∇ϕ=(1/Vcell)∮Scellϕn→ds. The method, though appearing to be stable and convergent within smooth regions, produces an O(1) error near corners and edges. The error produced is so large that the sharp features would be smeared out during evolving and reinitializing. Another problem caused by sharp features is that the up-wind discretization [17] fails near the angle bisector of edges or corners. On different sides of the angle bisector, ∇ϕ is discontinuous. The discontinuity may lead ϕx- and ϕx+ to differ in sign, and eventually causes simple up-wind gradient estimation to fail to obtain the correct domain of dependence [19] p57. The problem of the up-wind gradient is not significant when the level-set method is applied to flow simulation, since fluids can hardly develop edges or corners during their evolution. However, if we need to evolve fronts with sharp features, the broken up-wind gradient can be a serious problem.This article addresses development of a robust level-set approach to evolve geometries having sharp features on tetrahedral meshes, but the proposed approach can also be easily extended to hexahedral meshes because each hexahedral cell can be decomposed into five or six tetrahedrons without introducing new vertices [20]. In order to solve the two above-mentioned problems introduced by sharp features, we have developed an innovative spatial discretization scheme and a gradient-estimating technology matching the scheme. The new spatial discretization imitates the majority of behaviors of the Godunov’s scheme [21], which ensures numerical stability in discontinuous regions without introducing artificial viscosity like the Lax-Friedrichs or the Roe-Fix schemes. The gradient-estimating technology is more flexible than the divergence theorem-based method, thus allowing us to estimate the gradient of ϕ on any specific side of a vertex. This property is important to the proposed spatial discretization, which requires to perform gradient estimation on specifiable regions.Compared to previous works, our level-set approach has the following advantages: (a) It is robust on body-fitted tetrahedral meshes of complex geometries, including those composed of multiple materials; (b) Sharp features are preserved during evolving; (c) No ghost cell or similar tricks are required, meaning that the level-set evolving and coupled physical simulations could be run on the same mesh; (d) The∇ϕ estimation used for evolving and reinitializing are obtained through unified flow, so the dual grid computation [17] is evaded. These properties imply that our approach has utility in a wide range of applications.


## 2. Nomenclature

The nomenclature used in this paper is illustrated by Figure1. Within each cell, the four vertices are denoted by Vi, where i∈1,2,3,4 should be attached to the vertices such that V2V1→·V2V4→×V2V3→>0. In Figure 1, V1 is the vertex for which we need to estimate ∇ϕV1. Such a vertex is referred astarget vertex later. The midpoint between Vi and Vj is referred to as Eij. The barycenter of triangle facet VjVkVl is referred to as Fi, where i∪j,k,l=1,2,3,4. A vector which goes from the target vertex to any of the neighboring vertices, midpoints, or barycenters P is referred as anoutgoing vector to P. S(V1) is the set composed of all outgoing vectors that start from V1.Figure 1
Tetrahedral cell nomenclature demonstration.The signed distance field function of the evolving front is defined on each of the vertices and referred to asϕ(V→i). The linear estimation of the directional derivative of ϕ along the outgoing vector to P is referred to as Rc(P) or Rc(V1P→), e.g., Rc(V2)=(ϕV2-ϕV1)/V1V2, Rc(V1F1→)=ϕF1-ϕV1/|V1F1|. Rc denotes the changing rate of ϕ.


## 3. Governing Equations

The level-set equation Eq.(1) describes the movement of the isosurface x→∣ϕx→=0, where un(x→) is the normal speed field at which the front evolves.Usually,un(x→) is not constant across the entire field, so ϕ would no longer be an exact signed distance field after each time-forward iteration step, and we have to reinitialize ϕ regularly. The reinitialization is performed via Eq.(2) [22], where τ is artificial time and S(ϕ)=ϕ/ϕ2+Δx2 (Δx is the cubic root of the volume of the smallest tetrahedral cell).(1)∂ϕ∂t+unx→∇ϕ=0(2)∂ϕ∂τ+Sϕ∇ϕ-1=0Equations (1) and (2) can be solved with respect to time via the Euler method or the total variation diminishing Runge-Kutta (TVD-RK) method. Both of the two time-discretization schemes are mature and straightforward, and therefore are not discussed in detail here. The core problem to solve in (1) and (2) is spatial discretization. Since handling sharp features and discontinuity is required, we have not used the simple up-wind scheme, but developed an innovative spatial discretization, which will be discussed in Section 4.1. The newly constructed estimating technologies will be presented in Sections 4.2–4.4.


## 4. Numerical Method

### 4.1. Spatial Discretization

The Godunov’s spatial discretization scheme is widely used in numerical calculation studies to handle discontinuity. It has different implementations in different situations. The dimension-by-dimension form, which is widely used in level set studies, can be summarized by the following Eq.(3) [19] p58:(3)C=unDuringevolvingϕx→Duringreinitializingϕa2≈max⁡max⁡ϕa-,02,min⁡ϕa+,02C>0max⁡min⁡ϕa-,02,max⁡ϕa+,02C<0where a is any of axes, and C is the coefficient of the Hamiltonians. |∇ϕ| is then computed via:(4)∇ϕ¯=∑a∈x,y,zϕa2While being widely used in previous studies, the above scheme triggers two problems when applied to our study: (a) With non-Cartesian meshes, Eq.(3) cannot be directly applied due to that the cell edges are not axis-aligned. (b) In the vicinity of sharp features (e.g. 2D corners, 3D edges or corners), Eq.(3) may produce ϕx2≈ϕy2≈ϕz2≈1 on however fine grid, leading Eq.(4) to output 2 or 3 and smearing sharp geometry features out during reinitialization.Problem (a) can be settled by a newly developed edge-based estimation method, which will be presented later (Section4.3). The method works on any specifiable side of a vertex and directly estimate the gradient, instead of computing several directional derivatives like ϕx- and ϕx+. However, in order to solve problem (b) elegantly, we need a new spatial discretization scheme, which has to meet the following criteria:(C.1)
In the vicinity of sharp geometry, the scheme shall not introduce theO(1) error like Eq.(4);(C.2)
At peaks or valleys of the field, the result shall tend to that of Eq.(3);(C.3)
In regions whereϕa- and ϕa+ have the same signs, the result shall tend to that of Eq.(3)-(4);(C.4)
In other regions where the situations on the axes are different (e.g. the field is smooth along one axis while there are peaks or valleys along other axes), the scheme shall perform the estimation within reasonable domains of dependence.The newly developed scheme is presented below. We will prove that the presented scheme fits the above criteria later.The proposed scheme defines the directionm→, where ϕ has the maximum directional derivative along m→, as the “main axis”. Beside each vertex, two estimating operations are performed on direction m→ and -m→ respectively. Several main axes and corresponding domains of dependence are demonstrated in Figure 2, where thin blue lines are isolines of ϕ, red arrows are main axes, and gray areas denote the domains of dependence where ϕm- and ϕm+ are estimated. Once the two gradient estimation are obtained, the following Eq.(5)-(6) are integrated to finally get ∇ϕ¯.(5)t-=Cm→·∇ϕm-t+=Cm→·∇ϕm+(6)∇ϕ¯=∇ϕm+t-<0andt+<0∇ϕm-t->0andt+>00→t-<0andt+>0∇ϕm+t->0andt+<0and∇ϕm+>∇ϕm-∇ϕm-t->0andt+<0and∇ϕm+<∇ϕm-Figure 2
Main axes and domains of dependence of several target vertices.For (C.1), it can be seen that it is satisfied only if the gradient estimation does not contain O(1) error. In Section 4.4, we have deduced the root cause of the O(1) error and proposed the “explicit correction” technology to ensure (C.1) to be satisfied.For (C.2), peaks and valleys can be discussed separately. At peaks of the field, both ∇ϕm- and ∇ϕm+ point to the target vertex from their respective dependent domains. Therefore, we always have ∇ϕm-·m→>0 and ∇ϕm+·m→<0. Substituting the values into Eq.(5)-(6), it can be found that the behavior of the proposed scheme is in accordance with that of the Godunov-type scheme: when C>0, either ∇ϕm- or ∇ϕm+ is chosen to depend on which one has the maximum norm; when C<0, 0→ is outputted as ∇ϕ¯. Similarly, at valleys of the field, we have ∇ϕm-·m→<0 and ∇ϕm+·m→>0. When C>0, 0→ is outputted; when C<0, the one having larger norm is chosen. It is evident that in both cases concerning (C.2), the proposed scheme behaves in the same way as the Godunov-type scheme, and (C.2) is satisfied.In cases concerning (C.3), ∀a∈{x,y,z}, ϕa- and ϕa+ have the same signs. As a result, m→ would be appropriately parallel to ∇ϕm- and ∇ϕm+, and m→·∇ϕm->0, m→·∇ϕm+>0. According to Eq.(5)-(6), the proposed scheme outputs ∇ϕm- when C>0 and ∇ϕm+ when C<0. Meanwhile, if ϕa- and ϕa+ are both positive, m→ points to the positive side of axis a, thus the component of ∇ϕm- on a is approximately equal to ϕa-, and that of ∇ϕm+ is approximately ϕa+. Considering Eq.(3), which takes ϕa- when C>0 and ϕa+ when C<0, it can be seen that the proposed scheme takes the same domain of dependence and gives close result compared to the Godunov-type approach. When ϕa- and ϕa+ are both negative, the case can be proved similarly. The requirement of (C.3) is satisfied.In cases concerning (C.4), the proposed scheme operates in a slightly different manner with the Godunov-type scheme, but the behavior can be explained reasonably. Since that there are too many cases to analysis one by one, and that all cases can be explained via similar patterns, we only provide the analysis for two typical 2D cases, which are demonstrated in Figure 3. For both vertex A and B, it is clear that ϕx-<0, ϕx+>0, ϕy->0, ϕy+>0. According to Eq.(3)-(4), when C>0, 02+(ϕy-)2 will be taken as the estimation of ∇ϕ¯. For A, ϕy equals to the norm of local gradient, so the estimation is reasonable. However, ϕy estimated from beside B is significantly smaller than the norm of actual gradient, and the above estimation introduces error.Figure 3
A sample case concerning (C.4): semi-concave region.On the other hand, the proposed scheme processes the two vertices in different manners. ForA, m→ is parallel to the corner bisector, so both m→·∇ϕm- and m→·∇ϕm+ are positive. Eq.(6) will choose ∇ϕm- when C>0 and ∇ϕm+ when C<0, meaning that the proposed scheme is still in accordance with the Godunov-type scheme. However, for vertex B, m→ is perpendicular to one of the adjacent edges, depending on the numerical disturbance near B. In this situation, m→·∇ϕm+ is always positive, but the sign of m→·∇ϕm- depends on the sharpness of the corner. If the corner is relatively sharp, ∇ϕm- would diverge so far from m→ that m→·∇ϕm-<0, and the proposed scheme will treat this situation as a valley. Conversely, if the corner is so smooth that the direction of ∇ϕm- is close to m→, both m→·∇ϕm- and   m→·∇ϕm+ would be positive. and the situation will be treated as a smooth region.To sum up, the proposed scheme treats the vertex as a valley when “valley” is the dominant property, and otherwise chooses either side of the bisector as the information source to evolve the target vertex. Considering that the target vertex lies extremely close to the angle bisector, it may be influenced by the information propagated from either side. In some ways, choosing one side can be viewed as shifting the vertex towards the chosen side by a tiny distance, and thus ensuring the chosen side to be the theoretically correct one. Although the behavior introduces error for potentially incorrect choices, when the scheme is used for solving Eq.(1)-(2), the introduced error is small due to the symmetry across both sides of the bisector. Numerical experiments and validations show that the produced result is stable and accurate enough, implying the proposed scheme to satisfy (C.4) within the scope of this article.A problem which remains unsolved is about establishing the main axis. Theoretically,m→ shall be established according to the central-difference estimation of the gradient. However, such operation requires building dual grid structure [17] and costs extra computational resource. In this paper, we establish m→ such that m→ is parallel to the outgoing vector M→ which has the maximum |Rc|:(7)M→=argmaxVP→∈SV⁡RcVP→For vertices adjacent to theϕ=0 isosurface, since the computing is aimed at evolving the ϕ=0 isosurface, information is supposed to propagate from the isosurface to other regions of the mesh. Therefore, the correct domains of dependence are always on the nearer sides to the ϕ=0 isosurface, and the “explicit correction” technology is applied to obtain proper estimations for such vertices.


### 4.2. Cell-Based Gradient Estimation

Here, we briefly illustrate the problem of the cell-based approach before introducing our approach. The core of the cell-based estimating approach can be summarized by Eq. (8). The cell-centered estimations of ∇ϕ obtained from (8) can then be integrated into various difference schemes, such as the up-wind approach [17] and the WENO approach [15, 16]:(8)∇ϕ=1Vcell∮Scellϕn→ds.Such a cell-based approach, though used in a large portion of previous research, shows two major disadvantages when applied to feature-rich geometries: (a) The orientations of cells with respect to the target vertex cannot be well defined, making it difficult to apply the above-presented spatial discretization. The problem is illustrated in Figure4. (b) As we have proved in the Appendix, cell-based estimation methods rely on a pre-condition that the gradient within one cell tends to be uniform as the mesh grows finer. The pre-condition is true only if the zero level set {x→∣ϕ(x→)=0} is c1 continuous or smoother. Near sharp features (where the zero level set is c0 continuous), the cell-based approach produces an O(1) error.Figure 4
Domain of dependence of the edge- and cell-based approaches. The gray region in the figure shows the ideal domain of dependence with which to estimate∇m-ϕ. △OV1V2, △OV1V6, and △OV5V6 overlap with the gray region. Including cells △OV1V2 and △OV5V6 in the calculation introduces information from V2 and V5, which are not supposed to be included since they lie outside the gray region. Excluding △OV1V2 and △OV5V6, however, results in a lack of robustness.In this article, we have developed the edge-based gradient estimation approach to solve problem (a), and the “explicit correction” technology to settle problem (b). The two methods is described in the following Sections.


### 4.3. Edge-Based Gradient Estimation

Assuming that∇ϕ is uniform within one differential volume, we know from the Appendix that ∀VP→ within this region, Rc(VP→)≡∇ϕ·VP→/|VP→|. Conversely, in N-dimensional space, we can use N non-linear correlative vectors and corresponding changing-rate values to uniquely identify ∇ϕ, as shown in Figure 5. This approach, which uses outgoing vectors instead of cells as sampling primitives to fetch gradient information, is referred to as the edge-based approach.Figure 5
Edge-based∇ϕ estimating.
(a)
2D (b)
3DSince outgoing vectors are smaller primitives than cells, and all outgoing vectors pass the target vertex, the edge-based approach allows finer control over the domain of dependence compared to the cell-based approach. In Figure4, the red vectors demonstrate acceptable outgoing vectors to estimate ∇m-ϕ. Obviously, for each outgoing vector, it is easy to unambiguously tell whether the outgoing vector belongs to the desired domain of dependence or not.Lete→ be the vector that points to the desired domain of dependence. Linear equation set (9) estimates ∇ϕ at point O→ upon the desired domain of dependence, where TG is the threshold value that controls how close the select outgoing vectors should be to e→. For the outgoing vectors that end at midpoints or barycenters, the values of ϕ on their ends can be computed via interpolation.(9)V→·∇ϕ=V→RcV→∀V→∈SO,andV→·e→<TGIn order to solve (9) in N-dimensional space, N outgoing vectors meeting the threshold condition are required. Practically, there are often more outgoing vectors than required, especially in 3D meshes. In this case, an overdetermined linear system can be established to include all the available outgoing vectors. The system can then be solved via the least-squares method. For overdetermined linear systems, multiplying a constant value on both sides on an equation is equivalent to setting the weight of this equation in the system. We assign larger weights to the outgoing vectors that are closer to e→ by using V→·e→ as the weight. The complete form of the estimating equation set can be described by (10):(10)w=V→·e→wV→·∇ϕ=wV→RcV→∀V→∈SO,andw<TGNumerical consistency is a required property for discrete Hamiltonians [14]. The analysis below will demonstrate that Eq. (10) provides an exact estimation for linear fields. Letting ϕL be a linear field, we have(11)ϕLx→=ϕLz→+x→-z→·g→,where z→ is an arbitrary constant point and g→ is the gradient of field ϕ. We rewrite Eq. (10) as the matrix form Ar→=b, where r→ is the unknown gradient estimation obtained from (10). Let V0 be the target vertex and V1,V2,…,Vi be the ends of all acceptable outgoing vectors; then, A and b can be expressed by the following partitioned matrix:(12a)A=w1V0V1→w2V0V2→⋯wiV0Vi→(12b)b=w1V0V1RcV0V1→w2V0V2RcV0V2→⋯wiV0ViRcV0Vi→=w1ϕLV1→-ϕLV0→w2ϕLV2→-ϕLV0→⋯w1ϕLVi→-ϕLV0→Substituting (11) into (12b), we have(13)b=w1V1→-V0→·g→w2V2→-V0→·g→⋯w1Vi→-V0→·g→=AgAs discussed above, we solve (10) via the least-squares method,(14)r→=ATA-1ATbSubstituting (13) into (14) we have(15)r→=ATA-1ATAg=gFrom (15), it is obvious that the proposed edge-based approach gives an exact estimation result for linear fields.


### 4.4. Gradient Estimation Near Sharp Features

While the edge-based method settles the problem of estimating∇ϕ upon the desired domain of dependence, the problem of the O(1) error near sharp features remains unresolved (although using smaller TG eases the problem). To achieve higher accuracy, it is obvious from the above discussion that the uniform gradient assumption must be discarded.The root reason of the problem is that decreasing the cell size near sharp features does not improve the uniformity of the gradient within single cells. Since only the cells near sharp features (or more broadly, near the zero level set) are affected, a reasonable solution is integrating another estimating approach that does not rely on the uniform gradient assumption for these cells. One useful property of signed distance fields is that the zero level set is perpendicular to the local gradient. In discontinuous regions such as the outer sector of the corner in Figure17(a), the direction of gradient starts from the corner and ends at the target vertex. Generally, for a vertex V that locates near the ϕ=0 isosurface, letting the nearest point to V on the isosurface be N, the corresponding gradient ϕV is parallel to NV→ or VN→. Assuming that the gradient distributes uniformly on line segment NV, the norm of the gradient equals to the changing rate on NV. The above process can be summarized by Eq. (16):(16)N→=argmin∀N→∣ϕN→=0⁡NV,∇ϕ¯∇ϕ¯=NV→NV→ϕV→>0,-NV→NV→ϕV→<0,∇ϕ¯=ϕV→NV.The problem is now reduced to locating the nearest pointN→ for each target vertex near the level set. In [23], a brute-force approach is proposed to search N→ in all adjacent cells in an explicit manner. We have employed this approach, together with a minor improvement, to search N→ more reliably. The approach used in this article can be summarized by Procedure 1.Procedure 1:Search nearest point to V→ on {x→∣ϕ(x→)=0}.
Input: V →, ϕ 1: Cell setC={Ci∣V→∈Ci} 2: Vertex setS={Vi→∣Vi→∈C} 3:if∀Vi→∈S, sign(ϕ(Vi→))=sign(ϕ(V→))then 4:V is not adjacent to the level set. Quit procedure and launch edge-based approach. 5:end if 6:for all Cells Ci∈Cdo 7:  Find the line segment (2D) or facet (3D)Zi in Ci belonging to {x→∣ϕ(x→)=0} via linear interpolation. 8:  Record the minimum distanceD and corresponding nearest point N 9:end for 10:if∃Vi→∈S that SN<Dthen 11:  Add neighboring vertices ofVi to S 12:  Repeat step 6 - 11 until step 10 returnsfalse. 13:end if 14:returnNV→, DCompared to its original form, in Procedure1 we have added steps 10 - 13 to ensure that the outputted nearest point is the globally nearest one. Other details about the procedure can be found in [23] and are not discussed here. The method proposed here is referred to as “explicit correction (EC)” below.


### 4.5. Overall Flow

The three core modules of our evolving approach, namely the spatial discretization, the edge-based method, and the explicit correction approach, were discussed above. Combining them, we have a robust approach to solve Eq.(1)-(2). Other aspects of our evolving flow are straightforward. The entire flow is demonstrated in Figure 6.Figure 6
Overall flow.The time integration we have used to solve Eq.(1)-(2) is the two-stage Runge-Kutta method, which can be summarized by equation (17):(17)ϕVn+1/2=ϕVn-12ΔtHVn,ϕVn+1=ϕVn-ΔtHVn+1/2.


## 4.1. Spatial Discretization

The Godunov’s spatial discretization scheme is widely used in numerical calculation studies to handle discontinuity. It has different implementations in different situations. The dimension-by-dimension form, which is widely used in level set studies, can be summarized by the following Eq.(3) [19] p58:(3)C=unDuringevolvingϕx→Duringreinitializingϕa2≈max⁡max⁡ϕa-,02,min⁡ϕa+,02C>0max⁡min⁡ϕa-,02,max⁡ϕa+,02C<0where a is any of axes, and C is the coefficient of the Hamiltonians. |∇ϕ| is then computed via:(4)∇ϕ¯=∑a∈x,y,zϕa2While being widely used in previous studies, the above scheme triggers two problems when applied to our study: (a) With non-Cartesian meshes, Eq.(3) cannot be directly applied due to that the cell edges are not axis-aligned. (b) In the vicinity of sharp features (e.g. 2D corners, 3D edges or corners), Eq.(3) may produce ϕx2≈ϕy2≈ϕz2≈1 on however fine grid, leading Eq.(4) to output 2 or 3 and smearing sharp geometry features out during reinitialization.Problem (a) can be settled by a newly developed edge-based estimation method, which will be presented later (Section4.3). The method works on any specifiable side of a vertex and directly estimate the gradient, instead of computing several directional derivatives like ϕx- and ϕx+. However, in order to solve problem (b) elegantly, we need a new spatial discretization scheme, which has to meet the following criteria:(C.1)
In the vicinity of sharp geometry, the scheme shall not introduce theO(1) error like Eq.(4);(C.2)
At peaks or valleys of the field, the result shall tend to that of Eq.(3);(C.3)
In regions whereϕa- and ϕa+ have the same signs, the result shall tend to that of Eq.(3)-(4);(C.4)
In other regions where the situations on the axes are different (e.g. the field is smooth along one axis while there are peaks or valleys along other axes), the scheme shall perform the estimation within reasonable domains of dependence.The newly developed scheme is presented below. We will prove that the presented scheme fits the above criteria later.The proposed scheme defines the directionm→, where ϕ has the maximum directional derivative along m→, as the “main axis”. Beside each vertex, two estimating operations are performed on direction m→ and -m→ respectively. Several main axes and corresponding domains of dependence are demonstrated in Figure 2, where thin blue lines are isolines of ϕ, red arrows are main axes, and gray areas denote the domains of dependence where ϕm- and ϕm+ are estimated. Once the two gradient estimation are obtained, the following Eq.(5)-(6) are integrated to finally get ∇ϕ¯.(5)t-=Cm→·∇ϕm-t+=Cm→·∇ϕm+(6)∇ϕ¯=∇ϕm+t-<0andt+<0∇ϕm-t->0andt+>00→t-<0andt+>0∇ϕm+t->0andt+<0and∇ϕm+>∇ϕm-∇ϕm-t->0andt+<0and∇ϕm+<∇ϕm-Figure 2
Main axes and domains of dependence of several target vertices.For (C.1), it can be seen that it is satisfied only if the gradient estimation does not contain O(1) error. In Section 4.4, we have deduced the root cause of the O(1) error and proposed the “explicit correction” technology to ensure (C.1) to be satisfied.For (C.2), peaks and valleys can be discussed separately. At peaks of the field, both ∇ϕm- and ∇ϕm+ point to the target vertex from their respective dependent domains. Therefore, we always have ∇ϕm-·m→>0 and ∇ϕm+·m→<0. Substituting the values into Eq.(5)-(6), it can be found that the behavior of the proposed scheme is in accordance with that of the Godunov-type scheme: when C>0, either ∇ϕm- or ∇ϕm+ is chosen to depend on which one has the maximum norm; when C<0, 0→ is outputted as ∇ϕ¯. Similarly, at valleys of the field, we have ∇ϕm-·m→<0 and ∇ϕm+·m→>0. When C>0, 0→ is outputted; when C<0, the one having larger norm is chosen. It is evident that in both cases concerning (C.2), the proposed scheme behaves in the same way as the Godunov-type scheme, and (C.2) is satisfied.In cases concerning (C.3), ∀a∈{x,y,z}, ϕa- and ϕa+ have the same signs. As a result, m→ would be appropriately parallel to ∇ϕm- and ∇ϕm+, and m→·∇ϕm->0, m→·∇ϕm+>0. According to Eq.(5)-(6), the proposed scheme outputs ∇ϕm- when C>0 and ∇ϕm+ when C<0. Meanwhile, if ϕa- and ϕa+ are both positive, m→ points to the positive side of axis a, thus the component of ∇ϕm- on a is approximately equal to ϕa-, and that of ∇ϕm+ is approximately ϕa+. Considering Eq.(3), which takes ϕa- when C>0 and ϕa+ when C<0, it can be seen that the proposed scheme takes the same domain of dependence and gives close result compared to the Godunov-type approach. When ϕa- and ϕa+ are both negative, the case can be proved similarly. The requirement of (C.3) is satisfied.In cases concerning (C.4), the proposed scheme operates in a slightly different manner with the Godunov-type scheme, but the behavior can be explained reasonably. Since that there are too many cases to analysis one by one, and that all cases can be explained via similar patterns, we only provide the analysis for two typical 2D cases, which are demonstrated in Figure 3. For both vertex A and B, it is clear that ϕx-<0, ϕx+>0, ϕy->0, ϕy+>0. According to Eq.(3)-(4), when C>0, 02+(ϕy-)2 will be taken as the estimation of ∇ϕ¯. For A, ϕy equals to the norm of local gradient, so the estimation is reasonable. However, ϕy estimated from beside B is significantly smaller than the norm of actual gradient, and the above estimation introduces error.Figure 3
A sample case concerning (C.4): semi-concave region.On the other hand, the proposed scheme processes the two vertices in different manners. ForA, m→ is parallel to the corner bisector, so both m→·∇ϕm- and m→·∇ϕm+ are positive. Eq.(6) will choose ∇ϕm- when C>0 and ∇ϕm+ when C<0, meaning that the proposed scheme is still in accordance with the Godunov-type scheme. However, for vertex B, m→ is perpendicular to one of the adjacent edges, depending on the numerical disturbance near B. In this situation, m→·∇ϕm+ is always positive, but the sign of m→·∇ϕm- depends on the sharpness of the corner. If the corner is relatively sharp, ∇ϕm- would diverge so far from m→ that m→·∇ϕm-<0, and the proposed scheme will treat this situation as a valley. Conversely, if the corner is so smooth that the direction of ∇ϕm- is close to m→, both m→·∇ϕm- and   m→·∇ϕm+ would be positive. and the situation will be treated as a smooth region.To sum up, the proposed scheme treats the vertex as a valley when “valley” is the dominant property, and otherwise chooses either side of the bisector as the information source to evolve the target vertex. Considering that the target vertex lies extremely close to the angle bisector, it may be influenced by the information propagated from either side. In some ways, choosing one side can be viewed as shifting the vertex towards the chosen side by a tiny distance, and thus ensuring the chosen side to be the theoretically correct one. Although the behavior introduces error for potentially incorrect choices, when the scheme is used for solving Eq.(1)-(2), the introduced error is small due to the symmetry across both sides of the bisector. Numerical experiments and validations show that the produced result is stable and accurate enough, implying the proposed scheme to satisfy (C.4) within the scope of this article.A problem which remains unsolved is about establishing the main axis. Theoretically,m→ shall be established according to the central-difference estimation of the gradient. However, such operation requires building dual grid structure [17] and costs extra computational resource. In this paper, we establish m→ such that m→ is parallel to the outgoing vector M→ which has the maximum |Rc|:(7)M→=argmaxVP→∈SV⁡RcVP→For vertices adjacent to theϕ=0 isosurface, since the computing is aimed at evolving the ϕ=0 isosurface, information is supposed to propagate from the isosurface to other regions of the mesh. Therefore, the correct domains of dependence are always on the nearer sides to the ϕ=0 isosurface, and the “explicit correction” technology is applied to obtain proper estimations for such vertices.


## 4.2. Cell-Based Gradient Estimation

Here, we briefly illustrate the problem of the cell-based approach before introducing our approach. The core of the cell-based estimating approach can be summarized by Eq. (8). The cell-centered estimations of ∇ϕ obtained from (8) can then be integrated into various difference schemes, such as the up-wind approach [17] and the WENO approach [15, 16]:(8)∇ϕ=1Vcell∮Scellϕn→ds.Such a cell-based approach, though used in a large portion of previous research, shows two major disadvantages when applied to feature-rich geometries: (a) The orientations of cells with respect to the target vertex cannot be well defined, making it difficult to apply the above-presented spatial discretization. The problem is illustrated in Figure4. (b) As we have proved in the Appendix, cell-based estimation methods rely on a pre-condition that the gradient within one cell tends to be uniform as the mesh grows finer. The pre-condition is true only if the zero level set {x→∣ϕ(x→)=0} is c1 continuous or smoother. Near sharp features (where the zero level set is c0 continuous), the cell-based approach produces an O(1) error.Figure 4
Domain of dependence of the edge- and cell-based approaches. The gray region in the figure shows the ideal domain of dependence with which to estimate∇m-ϕ. △OV1V2, △OV1V6, and △OV5V6 overlap with the gray region. Including cells △OV1V2 and △OV5V6 in the calculation introduces information from V2 and V5, which are not supposed to be included since they lie outside the gray region. Excluding △OV1V2 and △OV5V6, however, results in a lack of robustness.In this article, we have developed the edge-based gradient estimation approach to solve problem (a), and the “explicit correction” technology to settle problem (b). The two methods is described in the following Sections.


## 4.3. Edge-Based Gradient Estimation

Assuming that∇ϕ is uniform within one differential volume, we know from the Appendix that ∀VP→ within this region, Rc(VP→)≡∇ϕ·VP→/|VP→|. Conversely, in N-dimensional space, we can use N non-linear correlative vectors and corresponding changing-rate values to uniquely identify ∇ϕ, as shown in Figure 5. This approach, which uses outgoing vectors instead of cells as sampling primitives to fetch gradient information, is referred to as the edge-based approach.Figure 5
Edge-based∇ϕ estimating.
(a)
2D (b)
3DSince outgoing vectors are smaller primitives than cells, and all outgoing vectors pass the target vertex, the edge-based approach allows finer control over the domain of dependence compared to the cell-based approach. In Figure4, the red vectors demonstrate acceptable outgoing vectors to estimate ∇m-ϕ. Obviously, for each outgoing vector, it is easy to unambiguously tell whether the outgoing vector belongs to the desired domain of dependence or not.Lete→ be the vector that points to the desired domain of dependence. Linear equation set (9) estimates ∇ϕ at point O→ upon the desired domain of dependence, where TG is the threshold value that controls how close the select outgoing vectors should be to e→. For the outgoing vectors that end at midpoints or barycenters, the values of ϕ on their ends can be computed via interpolation.(9)V→·∇ϕ=V→RcV→∀V→∈SO,andV→·e→<TGIn order to solve (9) in N-dimensional space, N outgoing vectors meeting the threshold condition are required. Practically, there are often more outgoing vectors than required, especially in 3D meshes. In this case, an overdetermined linear system can be established to include all the available outgoing vectors. The system can then be solved via the least-squares method. For overdetermined linear systems, multiplying a constant value on both sides on an equation is equivalent to setting the weight of this equation in the system. We assign larger weights to the outgoing vectors that are closer to e→ by using V→·e→ as the weight. The complete form of the estimating equation set can be described by (10):(10)w=V→·e→wV→·∇ϕ=wV→RcV→∀V→∈SO,andw<TGNumerical consistency is a required property for discrete Hamiltonians [14]. The analysis below will demonstrate that Eq. (10) provides an exact estimation for linear fields. Letting ϕL be a linear field, we have(11)ϕLx→=ϕLz→+x→-z→·g→,where z→ is an arbitrary constant point and g→ is the gradient of field ϕ. We rewrite Eq. (10) as the matrix form Ar→=b, where r→ is the unknown gradient estimation obtained from (10). Let V0 be the target vertex and V1,V2,…,Vi be the ends of all acceptable outgoing vectors; then, A and b can be expressed by the following partitioned matrix:(12a)A=w1V0V1→w2V0V2→⋯wiV0Vi→(12b)b=w1V0V1RcV0V1→w2V0V2RcV0V2→⋯wiV0ViRcV0Vi→=w1ϕLV1→-ϕLV0→w2ϕLV2→-ϕLV0→⋯w1ϕLVi→-ϕLV0→Substituting (11) into (12b), we have(13)b=w1V1→-V0→·g→w2V2→-V0→·g→⋯w1Vi→-V0→·g→=AgAs discussed above, we solve (10) via the least-squares method,(14)r→=ATA-1ATbSubstituting (13) into (14) we have(15)r→=ATA-1ATAg=gFrom (15), it is obvious that the proposed edge-based approach gives an exact estimation result for linear fields.


## 4.4. Gradient Estimation Near Sharp Features

While the edge-based method settles the problem of estimating∇ϕ upon the desired domain of dependence, the problem of the O(1) error near sharp features remains unresolved (although using smaller TG eases the problem). To achieve higher accuracy, it is obvious from the above discussion that the uniform gradient assumption must be discarded.The root reason of the problem is that decreasing the cell size near sharp features does not improve the uniformity of the gradient within single cells. Since only the cells near sharp features (or more broadly, near the zero level set) are affected, a reasonable solution is integrating another estimating approach that does not rely on the uniform gradient assumption for these cells. One useful property of signed distance fields is that the zero level set is perpendicular to the local gradient. In discontinuous regions such as the outer sector of the corner in Figure17(a), the direction of gradient starts from the corner and ends at the target vertex. Generally, for a vertex V that locates near the ϕ=0 isosurface, letting the nearest point to V on the isosurface be N, the corresponding gradient ϕV is parallel to NV→ or VN→. Assuming that the gradient distributes uniformly on line segment NV, the norm of the gradient equals to the changing rate on NV. The above process can be summarized by Eq. (16):(16)N→=argmin∀N→∣ϕN→=0⁡NV,∇ϕ¯∇ϕ¯=NV→NV→ϕV→>0,-NV→NV→ϕV→<0,∇ϕ¯=ϕV→NV.The problem is now reduced to locating the nearest pointN→ for each target vertex near the level set. In [23], a brute-force approach is proposed to search N→ in all adjacent cells in an explicit manner. We have employed this approach, together with a minor improvement, to search N→ more reliably. The approach used in this article can be summarized by Procedure 1.Procedure 1:Search nearest point to V→ on {x→∣ϕ(x→)=0}.
Input: V →, ϕ 1: Cell setC={Ci∣V→∈Ci} 2: Vertex setS={Vi→∣Vi→∈C} 3:if∀Vi→∈S, sign(ϕ(Vi→))=sign(ϕ(V→))then 4:V is not adjacent to the level set. Quit procedure and launch edge-based approach. 5:end if 6:for all Cells Ci∈Cdo 7:  Find the line segment (2D) or facet (3D)Zi in Ci belonging to {x→∣ϕ(x→)=0} via linear interpolation. 8:  Record the minimum distanceD and corresponding nearest point N 9:end for 10:if∃Vi→∈S that SN<Dthen 11:  Add neighboring vertices ofVi to S 12:  Repeat step 6 - 11 until step 10 returnsfalse. 13:end if 14:returnNV→, DCompared to its original form, in Procedure1 we have added steps 10 - 13 to ensure that the outputted nearest point is the globally nearest one. Other details about the procedure can be found in [23] and are not discussed here. The method proposed here is referred to as “explicit correction (EC)” below.


## 4.5. Overall Flow

The three core modules of our evolving approach, namely the spatial discretization, the edge-based method, and the explicit correction approach, were discussed above. Combining them, we have a robust approach to solve Eq.(1)-(2). Other aspects of our evolving flow are straightforward. The entire flow is demonstrated in Figure 6.Figure 6
Overall flow.The time integration we have used to solve Eq.(1)-(2) is the two-stage Runge-Kutta method, which can be summarized by equation (17):(17)ϕVn+1/2=ϕVn-12ΔtHVn,ϕVn+1=ϕVn-ΔtHVn+1/2.


## 5. Validation

The following problems are tested to validate the accuracy and robustness of the proposed level-set method:(i)
(2D) Diffusion into a notched square. This example demonstrates the ability to evolve fronts at non-uniform evolving velocity on irregular regions.(ii)
(2D) Three merging circles. This example demonstrates the ability to handle topological changes.(iii)
(2D) Rate stick [5]. This example demonstrates the compatibility with slender and long regions.(iv)
(3D) Burning and erosion in a solid rocket motor. This example demonstrates the ability to evolve fronts at non-uniform evolving velocity on 3D irregular regions.(v)
(3D) Reconstructing a signed distance function [24]. This example demonstrates the convergence of the reinitializing operation.For the merging-circles example, the reinitializing operation is not integrated since the evolving speed distributes uniformly. Other details of each test problem will be discussed in the corresponding subsections. We have already published the simulation code of all 2D validating problems [25].The error metric we have used to measure the accuracy is the averaged and maximum errors ofϕ near the zero level set, which are defined as equation (18),(18)S=Vi→∣ϕexactVi→<5Δx,Ea=1nS∑i=1nSϕsimVi→-ϕexactVi→,Em=maxVi→∈S⁡ϕsimVi→-ϕexactVi→,where Δx is the scale of cells. Unless otherwise specified, Δx=0.01 in the following examples. The 2D meshes used for simulation are generated by the mesh2d [26] toolkit, while the 3D meshes are generated by Gmsh [27].


### 5.1. Diffusion into a Notched Square

This problem describes the diffusion starting at the upper edge of a U-notched square. Figure7 shows the computational mesh of the notched square, where the red line is the initial location of the field being diffused. Evolving velocity is defined as(19)unV→=0.250≤Vx→<0.15,0.50.15≤Vx→<0.6,10.6≤Vx→≤1,where Vx→ denotes the x coordinate of V→.Figure 7
Computational mesh of the notched square.In order to demonstrate the improvement introduced by explicit correction, in Figure8(a), we show two sets of transient results that are generated by the algorithm with and without explicit correction. The analytical solution is obtained via artificial geometry analysis. It can be seen that sharp geometric features are successfully preserved. The smearing-out effects near sharp features are not completely eliminated, but the explicit correction largely helps in relieving the effect. The error data, plotted in Figure 8(b), imply that the averaged error near the zero level set is kept below 0.7Δx at the end of the simulation. The peaks in the two Em curves at 100~115 steps are the product of the merging procedure of two sharp features near [0.15,0.8].Figure 8
Result of notched square example.
(a)
Evolving result (b)
Evolving error


### 5.2. Three Merging Circles

This example includes three circles, the centers of which are located at[0.3,0.3], [0.7,0.3], and [0.5,0.7], respectively. The radiuses of the three circles are initially 0.2 and grow uniformly at speed 1. We introduce this example to verify the ability of the proposed method to smoothly handle topological changes, especially the peak of ϕ located in the middle of the three circle centers.The result of this example is plotted in Figure9. It can be seen that the proposed method is able to solve the merging-circles problem with high accuracy. The averaged error is kept below 0.06Δx at the end of the simulation.Figure 9
Result of merging-circles example.
(a)
Evolving result (b)
Evolving error


### 5.3. Rate Stick

In the rate-stick problem, a front that starts from a single point and propagates through a slender and long tube is the object of study. The case is similar to the high-explosives (HE) rate-stick experiment described in [5]. Figure 10(a) demonstrates the mesh and initial status. The red point x0→, which denotes a zero-radius circle such that ϕ(x0→)=0, is the starting point of the front. un≡1 across the entire region. The mesh scale used in the example is Δx=0.005.Figure 10
Rate-stick problem.
(a)
Computational mesh (b)
Evolving result (c)
Evolving errorThe slender computational region means that the boundaries must be handled smoothly to evade distortion. The zero-radius initial region requires the proposed method to handle the case where the scale of the initial feature is smaller thanΔx. Figures 10(b) and 10(c) show the result and error, respectively. The result shows that the proposed method works well on the rate-stick problem. The peak in the Em curve is produced when the evolving circle is tangential to the mesh boundary.


### 5.4. Burning and Erosion in a Solid Rocket Motor

This problem contains real engineering components: the solid-rocket-motor grain and the corresponding thermal insulation layer. During the working procedure of the motor, the grain, which is made of solid propellant, burns and transforms into high-temperature gas. As the grain retreats, the thermal insulation layer is exposed to the gas and is slowly eroded. To summarize, the gas propagates into the grain and the insulation at different velocities. Figure11 demonstrates the working procedure on one section view of the motor, where the red lines outline the propagation of the front of the gas, and the gray and blue hatching mark the sections of the thermal insulation layer and the grain, respectively.Figure 11
Working procedure of a solid rocket motor.The computational mesh for this problem is shown in Figure12. In actual solid rocket motors the propagation speed is determined by local environmental parameters such as temperature and pressure. In this article, however, we set the speed as shown by Eq. (20) to ease the verification:(20)unx→=1x→∈grain,0.1x→∈insulation.Accordingly, the mesh for the insulation is finer than that for the grain to better capture the relatively slower front movement inside the insulation. Letting dt(x→) be the minimum distance from the insulation to x→, the mesh scale Δx is defined by Eq. (21):(21)Δxx→=2.4dtx→<2,1.8dtx→-1.22≤dtx→≤6,66<dtx→.Figure 12
Computational mesh of a solid rocket motor.Since the structure of the motor is so complex that an analytical solution is unavailable, we used the level-set method implementation described in [28] to produce a reference result and further compute the error via tri-linear interpolation across meshes. The mesh scale used to generate the reference result is Δx=Δy=Δz=0.6. Figure 13(a) shows two patches of the resulting isosurfaces (the coordinate systems to plot the two patches are aligned): the green piece is the right half of the zero level set at the 25th step and the red piece is the left half of the zero level set at the 80th step. The demonstrated patches are in line with the predictions made in Figure 11. The error data shown in Figure 13(b) reveal that the proposed method is consistent with the level-set method on Cartesian grids.Figure 13
Result of solid-rocket-motor example.
(a)
Evolving result (b)
Evolving error


### 5.5. Reconstructing a Signed Distance Function

This test is aimed at constructing an exact signed distance field from a distorted one. In this article, the source geometry of the sign distance fields is a cube. The initial distorted field is generated via Eq. (22a),(22b), and (22c),(22a)h=0.5,(22b)ϕ′V→=-min⁡h+Vx,h-Vx,h+Vy,h-Vy,h+Vz,h-Vz,(22c)ϕV→=ϕ′V→·Vx-0.32+Vy-0.22+Vz-0.12+0.05,where h is the half-side-length of the cube. Equation (22b) generates a signed field that contains discontinuity and its zero level set denotes the cube. Equation (22c) is the equation used by [24] to further distort the field. One section of the finally generated field is shown in Figure 14, where the thick red line marks its ϕ=0 isosurface. It can be seen that the output field contains two desired properties, namely the sharp features and the non-uniform distortion, which are generated on purpose to verify the ability of the proposed approach to reconstruct an exact signed distance function. The mesh scale used to generate the 3D mesh is Δx=0.05.Figure 14
Contours of initialϕ at section z=0.The result after 80 reinitializing passes is shown in Figure15. It can be seen that the signed distance field is correctly reconstructed and that the initial ϕ=0 isosurface is left unmoved. The error computed from the resulting ϕ field is Ea=0.082Δx, Em=0.54Δx. Such a result implies that the proposed approach can be used as a robust and accurate method to construct signed distance fields on unstructured meshes.Figure 15
Result of reconstructing the signed distance field.


### 5.6. Error Convergence

Here, we demonstrate the convergence of error as the mesh refines. The overall error across the evolving procedure is defined as Eq. (23),(23)Eo=1nT∑i=1nTEai,where n(T) is the number of evolving steps and Ea(i) denotes the Ea computed at step i. For the example that reconstructs a signed distance function, since the only error of concern is that of the last step, we use the corresponding Ea(n(T)) to measure its error.The error-convergence data is shown in Figure16. Since the above examples have different mesh scales, we have unified the mesh-scale axis in Figure 16 using the above-used Δx in corresponding examples. For the solid-rocket-motor example, the Δx used for unifying is 6. The error data of the solid rocket motor example corresponding to 4 are missing because Δx=4×6=24 is so large a scale that the unstructured mesh cannot be properly generated. It can be seen from Figure 16 that, as the mesh refines, the absolute Eo decreases. However, Eo/Δx has slightly increased at the same time. This implies that the error decreases more slowly than the mesh scale. Since the edge-based gradient-estimating approach is proved to have first-order accuracy, there must be other error sources that produce lower-order error and consequently cause Eo/Δx to increase as the mesh refines. Examples of such error sources include floating-point error, the remaining smearing-out effect near sharp features that cannot be completely eliminated by the explicit correction, or the error produced when converging explicit analytical results to implicit ones.Figure 16
Error convergence data.Figure 17
Cell-based∇ϕ estimation.
(a) (b)In order to quantify the convergence, we have performed Richardson’s analysis [29] modeled by Eq. (24), where C is a constant and p is the unknown order of convergence. By substituting the error data into (24), the order of convergence p of each example problem can be solved via the least-squares method. The results are collected in Table 1. The solid-rocket-motor example is excluded from the table since computing the order of convergence form one numerical method to another makes no sense. It can be seen from Table 1 that the implementation’s order of convergence is near 1. Theoretically, arbitrarily small error can be achieved by continually refining the mesh:(24)Eo=CΔxp+oΔx.Table 1
Error-converging orders.
Test example Converging order Notched square 0.8357 Three merging circles 0.9528 Rate stick 0.9121 Reconstructing signed distance fields 0.7586


## 5.1. Diffusion into a Notched Square

This problem describes the diffusion starting at the upper edge of a U-notched square. Figure7 shows the computational mesh of the notched square, where the red line is the initial location of the field being diffused. Evolving velocity is defined as(19)unV→=0.250≤Vx→<0.15,0.50.15≤Vx→<0.6,10.6≤Vx→≤1,where Vx→ denotes the x coordinate of V→.Figure 7
Computational mesh of the notched square.In order to demonstrate the improvement introduced by explicit correction, in Figure8(a), we show two sets of transient results that are generated by the algorithm with and without explicit correction. The analytical solution is obtained via artificial geometry analysis. It can be seen that sharp geometric features are successfully preserved. The smearing-out effects near sharp features are not completely eliminated, but the explicit correction largely helps in relieving the effect. The error data, plotted in Figure 8(b), imply that the averaged error near the zero level set is kept below 0.7Δx at the end of the simulation. The peaks in the two Em curves at 100~115 steps are the product of the merging procedure of two sharp features near [0.15,0.8].Figure 8
Result of notched square example.
(a)
Evolving result (b)
Evolving error


## 5.2. Three Merging Circles

This example includes three circles, the centers of which are located at[0.3,0.3], [0.7,0.3], and [0.5,0.7], respectively. The radiuses of the three circles are initially 0.2 and grow uniformly at speed 1. We introduce this example to verify the ability of the proposed method to smoothly handle topological changes, especially the peak of ϕ located in the middle of the three circle centers.The result of this example is plotted in Figure9. It can be seen that the proposed method is able to solve the merging-circles problem with high accuracy. The averaged error is kept below 0.06Δx at the end of the simulation.Figure 9
Result of merging-circles example.
(a)
Evolving result (b)
Evolving error


## 5.3. Rate Stick

In the rate-stick problem, a front that starts from a single point and propagates through a slender and long tube is the object of study. The case is similar to the high-explosives (HE) rate-stick experiment described in [5]. Figure 10(a) demonstrates the mesh and initial status. The red point x0→, which denotes a zero-radius circle such that ϕ(x0→)=0, is the starting point of the front. un≡1 across the entire region. The mesh scale used in the example is Δx=0.005.Figure 10
Rate-stick problem.
(a)
Computational mesh (b)
Evolving result (c)
Evolving errorThe slender computational region means that the boundaries must be handled smoothly to evade distortion. The zero-radius initial region requires the proposed method to handle the case where the scale of the initial feature is smaller thanΔx. Figures 10(b) and 10(c) show the result and error, respectively. The result shows that the proposed method works well on the rate-stick problem. The peak in the Em curve is produced when the evolving circle is tangential to the mesh boundary.


## 5.4. Burning and Erosion in a Solid Rocket Motor

This problem contains real engineering components: the solid-rocket-motor grain and the corresponding thermal insulation layer. During the working procedure of the motor, the grain, which is made of solid propellant, burns and transforms into high-temperature gas. As the grain retreats, the thermal insulation layer is exposed to the gas and is slowly eroded. To summarize, the gas propagates into the grain and the insulation at different velocities. Figure11 demonstrates the working procedure on one section view of the motor, where the red lines outline the propagation of the front of the gas, and the gray and blue hatching mark the sections of the thermal insulation layer and the grain, respectively.Figure 11
Working procedure of a solid rocket motor.The computational mesh for this problem is shown in Figure12. In actual solid rocket motors the propagation speed is determined by local environmental parameters such as temperature and pressure. In this article, however, we set the speed as shown by Eq. (20) to ease the verification:(20)unx→=1x→∈grain,0.1x→∈insulation.Accordingly, the mesh for the insulation is finer than that for the grain to better capture the relatively slower front movement inside the insulation. Letting dt(x→) be the minimum distance from the insulation to x→, the mesh scale Δx is defined by Eq. (21):(21)Δxx→=2.4dtx→<2,1.8dtx→-1.22≤dtx→≤6,66<dtx→.Figure 12
Computational mesh of a solid rocket motor.Since the structure of the motor is so complex that an analytical solution is unavailable, we used the level-set method implementation described in [28] to produce a reference result and further compute the error via tri-linear interpolation across meshes. The mesh scale used to generate the reference result is Δx=Δy=Δz=0.6. Figure 13(a) shows two patches of the resulting isosurfaces (the coordinate systems to plot the two patches are aligned): the green piece is the right half of the zero level set at the 25th step and the red piece is the left half of the zero level set at the 80th step. The demonstrated patches are in line with the predictions made in Figure 11. The error data shown in Figure 13(b) reveal that the proposed method is consistent with the level-set method on Cartesian grids.Figure 13
Result of solid-rocket-motor example.
(a)
Evolving result (b)
Evolving error


## 5.5. Reconstructing a Signed Distance Function

This test is aimed at constructing an exact signed distance field from a distorted one. In this article, the source geometry of the sign distance fields is a cube. The initial distorted field is generated via Eq. (22a),(22b), and (22c),(22a)h=0.5,(22b)ϕ′V→=-min⁡h+Vx,h-Vx,h+Vy,h-Vy,h+Vz,h-Vz,(22c)ϕV→=ϕ′V→·Vx-0.32+Vy-0.22+Vz-0.12+0.05,where h is the half-side-length of the cube. Equation (22b) generates a signed field that contains discontinuity and its zero level set denotes the cube. Equation (22c) is the equation used by [24] to further distort the field. One section of the finally generated field is shown in Figure 14, where the thick red line marks its ϕ=0 isosurface. It can be seen that the output field contains two desired properties, namely the sharp features and the non-uniform distortion, which are generated on purpose to verify the ability of the proposed approach to reconstruct an exact signed distance function. The mesh scale used to generate the 3D mesh is Δx=0.05.Figure 14
Contours of initialϕ at section z=0.The result after 80 reinitializing passes is shown in Figure15. It can be seen that the signed distance field is correctly reconstructed and that the initial ϕ=0 isosurface is left unmoved. The error computed from the resulting ϕ field is Ea=0.082Δx, Em=0.54Δx. Such a result implies that the proposed approach can be used as a robust and accurate method to construct signed distance fields on unstructured meshes.Figure 15
Result of reconstructing the signed distance field.


## 5.6. Error Convergence

Here, we demonstrate the convergence of error as the mesh refines. The overall error across the evolving procedure is defined as Eq. (23),(23)Eo=1nT∑i=1nTEai,where n(T) is the number of evolving steps and Ea(i) denotes the Ea computed at step i. For the example that reconstructs a signed distance function, since the only error of concern is that of the last step, we use the corresponding Ea(n(T)) to measure its error.The error-convergence data is shown in Figure16. Since the above examples have different mesh scales, we have unified the mesh-scale axis in Figure 16 using the above-used Δx in corresponding examples. For the solid-rocket-motor example, the Δx used for unifying is 6. The error data of the solid rocket motor example corresponding to 4 are missing because Δx=4×6=24 is so large a scale that the unstructured mesh cannot be properly generated. It can be seen from Figure 16 that, as the mesh refines, the absolute Eo decreases. However, Eo/Δx has slightly increased at the same time. This implies that the error decreases more slowly than the mesh scale. Since the edge-based gradient-estimating approach is proved to have first-order accuracy, there must be other error sources that produce lower-order error and consequently cause Eo/Δx to increase as the mesh refines. Examples of such error sources include floating-point error, the remaining smearing-out effect near sharp features that cannot be completely eliminated by the explicit correction, or the error produced when converging explicit analytical results to implicit ones.Figure 16
Error convergence data.Figure 17
Cell-based∇ϕ estimation.
(a) (b)In order to quantify the convergence, we have performed Richardson’s analysis [29] modeled by Eq. (24), where C is a constant and p is the unknown order of convergence. By substituting the error data into (24), the order of convergence p of each example problem can be solved via the least-squares method. The results are collected in Table 1. The solid-rocket-motor example is excluded from the table since computing the order of convergence form one numerical method to another makes no sense. It can be seen from Table 1 that the implementation’s order of convergence is near 1. Theoretically, arbitrarily small error can be achieved by continually refining the mesh:(24)Eo=CΔxp+oΔx.Table 1
Error-converging orders.
Test example Converging order Notched square 0.8357 Three merging circles 0.9528 Rate stick 0.9121 Reconstructing signed distance fields 0.7586


## 6. Conclusions

In this article, we developed a robust approach to evolve fronts in their normal direction on unstructured meshes. The approach is built on an innovative spatial discretization scheme and a gradient-estimating technology matching the scheme. An explicit correction approach is introduced to improve the accuracy near the zero level set. Results of the validating examples show that the proposed method handles sharp geometric features, discontinuous velocity field, and topological changes smoothly and accurately, on both 2D and 3D unstructured meshes. The method theoretically also applies to structured or hybrid meshes, not only due to that hexahedral cells can be decomposed into tetrahedrons, but also because the proposed gradient-estimating technology only relies on the outgoing vectors, which are available in all kinds of meshes.In practice, there are also cases where advecting fronts using externally generated velocity fields is required. We have also studied the advecting problem during this research. Unfortunately, it turned out that solving some of the advecting problems with acceptable accuracy requires higher-order spatial discretization schemes, which is beyond the scope of this article. In future research, we would study possible routes to integrate the proposed method with higher-order schemes and further expand the range of application of the method.



---

*Source: 2730829-2018-09-25.xml*