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

,


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 higherdimensional function  to model the evolving front, where the {  →  | (  →  ) = 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 evolvingfront 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 / + (∇,  →  ) = 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][12][13], WENO and Hermite WENO reconstruction methodologies were studied separately on 2D unstructured meshes; Augoula [14] developed several highorder 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 levelset 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/  )∮   (  →  )d.The method, though appearing to be stable and convergent within smooth regions, produces an (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 upwind discretization [17] fails near the angle bisector of edges or corners.On different sides of the angle bisector, ∇ is discontinuous.The discontinuity may lead  −  and  +  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 levelset 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 bodyfitted 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.

Nomenclature
The nomenclature used in this paper is illustrated by Figure 1.Within each cell, the four vertices are denoted by   , where  ∈ {1, 2, 3, 4} should be attached to the vertices such that In Figure 1,  1 is the vertex for which we need to estimate ∇  1 .Such a vertex is referred as target vertex later.The midpoint between   and   is referred to as   .The barycenter of triangle facet       is referred to as   , where {} ∪ {, , } = {1, 2, 3, 4}.A vector which goes from the target vertex to any of the neighboring vertices, midpoints, or barycenters  is referred as an outgoing vector to . ( 1 ) is the set composed of all outgoing vectors that start from  1 .
The signed distance field function of the evolving front is defined on each of the vertices and referred to as (  →   ).The linear estimation of the directional derivative of  along the outgoing vector to  is referred to as   (𝑃)

Governing Equations
The level-set equation Eq. (1) Equations ( 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 timediscretization 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.

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: where  is any of axes, and  is the coefficient of the Hamiltonians.|∇| is then computed via: While 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  2  ≈  2  ≈  2  ≈ 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 edgebased estimation method, which will be presented later (Section 4.3).The method works on any specifiable side of a vertex and directly estimate the gradient, instead of computing several directional derivatives like  −  and  +  .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 the (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  −  and  +  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 direction  → , where  has the maximum directional derivative along  → , as the "main axis".Beside each vertex, two estimating operations are performed on direction  →  and −  →  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  −  and  +  are estimated.Once the two gradient estimation are obtained, the following Eq.( 5)-( 6) are integrated to finally get ∇.
For (C.1), it can be seen that it is satisfied only if the gradient estimation does not contain (1) error.In Section 4.4, we have deduced the root cause of the (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 ∇ −  and ∇ +  point to the target vertex from their respective dependent domains.Therefore, we always have ∇ −  ⋅  →  > 0 and ∇ +  ⋅  →  < 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  > 0, either ∇ −  or ∇ +  is chosen to depend on which one has the maximum norm; when  < 0,  → 0 is outputted as ∇.Similarly, at valleys of the field, we have ∇ −  ⋅  →  < 0 and ∇ +  ⋅  →  > 0. When  > 0,  → 0 is outputted; when  < 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 Godunovtype scheme, and (C.2) is satisfied.
In cases concerning (C.3), which takes  −  when  > 0 and  +  when  < 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  −  and  +  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  and , it is clear that  −  < 0,  +  > 0,  −  > 0,  +  > 0. According to Eq.( 3)-( 4), when  > 0, √0 2 + ( −  ) 2 will be taken as the estimation of |∇|.For ,   equals to the norm of local gradient, so the estimation is reasonable.However,   estimated from beside  is significantly smaller than the norm of actual gradient, and the above estimation introduces error.
On the other hand, the proposed scheme processes the two vertices in different manners.For ,  →  is parallel to the corner bisector, so both  →  ⋅ ∇ −  and  →  ⋅ ∇ +  are positive.Eq.( 6) will choose ∇ −  when  > 0 and ∇ +  when  < 0, meaning that the proposed scheme is still in accordance with the Godunov-type scheme.However, for vertex ,  →  is perpendicular to one of the adjacent edges, depending on the numerical disturbance near .In this situation,  →  ⋅ ∇ +  is always positive, but the sign of  →  ⋅ ∇ −  depends on the sharpness of the corner.If the corner is relatively sharp, ∇ −  would diverge so far from  →  that  →  ⋅ ∇ −  < 0, and the proposed scheme will treat this situation as a valley.Conversely, if the corner is so smooth that the direction of ∇ −  is close to  → , both  → ⋅∇ −  and  → ⋅∇ +  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,  →  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  →  such that  →  is parallel to the outgoing vector  →  which has the maximum |  |: 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.

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]: 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 Figure 4. (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 precondition is true only if the zero level set {  →  | (  →  ) = 0} is 1 continuous or smoother.Near sharp features (where the zero level set is 0 continuous), the cell-based approach produces an (1) error.
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.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.

Edge-Based Gradient
Since outgoing vectors are smaller primitives than cells, and all outgoing vectors pass the target vertex, the edgebased approach allows finer control over the domain of dependence compared to the cell-based approach.In Figure 4, the red vectors demonstrate acceptable outgoing vectors to estimate ∇  − .Obviously, for each outgoing vector, it is easy to unambiguously tell whether the outgoing vector belongs to the desired domain of dependence or not.
Let  →  be the vector that points to the desired domain of dependence.Linear equation set ( 9) estimates ∇ at point  →  upon the desired domain of dependence, where   is the threshold value that controls how close the select outgoing vectors should be to  →  .For the outgoing vectors that end at midpoints or barycenters, the values of  on their ends can be computed via interpolation.
In order to solve (9) in -dimensional space,  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  →  by using  →  ⋅  →  as the weight.The complete form of the estimating equation set can be described by (10): Numerical 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   be a linear field, we have where  →  is an arbitrary constant point and  →  is the gradient of field .We rewrite Eq. ( 10) as the matrix form   →  = , where  →  is the unknown gradient estimation obtained from (10).
Let  0 be the target vertex and  1 ,  2 , . . .,   be the ends of all acceptable outgoing vectors; then,  and  can be expressed by the following partitioned matrix: Substituting ( 11) into (12b), we have As discussed above, we solve (10) via the least-squares method,  →  = (  ) Substituting ( 13) into ( 14) we have From (15), it is obvious that the proposed edge-based approach gives an exact estimation result for linear fields.

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 (1) error near sharp features remains unresolved (although using smaller   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 Figure 17(a), the direction of gradient starts from the corner and ends at the target vertex.Generally, for a vertex  that locates near the  = 0 isosurface, letting the nearest point to  on the isosurface be , the corresponding gradient   is parallel to   →  or   → .Assuming that the gradient distributes uniformly on line segment , the norm of the gradient equals to the changing rate on .The above process can be summarized by Eq. ( 16): The problem is now reduced to locating the nearest point  →  for each target vertex near the level set.In [23], a bruteforce approach is proposed to search  →  in all adjacent cells in an explicit manner.We have employed this approach, together with a minor improvement, to search  →  more reliably.The approach used in this article can be summarized by Procedure 1.
Compared to its original form, in Procedure 1 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.

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.
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):

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 nonuniform evolving velocity on irregular regions.
(ii) (2D) Three merging circles.This example demonstrates the ability to handle topological changes.
(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.
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), where Δ is the scale of cells.Unless otherwise specified, Δ = 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].

Diffusion into a Notched
Square.This problem describes the diffusion starting at the upper edge of a U-notched square.
Initialize 
Estimate ∇ via Procedure 1 Apply spatial discretization according to equations ( 5)-( 6 Figure 7 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 where  →   denotes the  coordinate of  → .In order to demonstrate the improvement introduced by explicit correction, in Figure 8(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Δ at the end of the simulation.The peaks in the two   curves at 100 ∼ 115 steps are the product of the merging procedure of two sharp features near [0.15, 0.8].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 Figure 9.It can be seen that the proposed method is able to solve the mergingcircles problem with high accuracy.The averaged error is kept below 0.06Δ at the end of the simulation.

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  →  0 , which denotes a zero-radius circle such that (  →  0 ) = 0, is the starting point of the front.  ≡ 1 across the entire region.The mesh scale used in the example is Δ = 0.005.The 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 Δ.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   curve is produced when the evolving circle is tangential to the mesh boundary.

Burning and Erosion in a Solid Rocket
Motor.This problem contains real engineering components: the solidrocket-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.Figure 11 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.
The computational mesh for this problem is shown in Figure 12.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: 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   (  →  ) be the minimum distance from the insulation to  →  , the mesh scale Δ is defined by Eq. ( 21): Since the structure of the motor is so complex that an analytical solution is unavailable, we used the levelset 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 Δ = Δ = Δ = 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.

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 where ℎ 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 Δ = 0.05.The result after 80 reinitializing passes is shown in Figure 15.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   = 0.082Δ,   = 0.54Δ.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.

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), where () is the number of evolving steps and   () denotes the   computed at step .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   (()) to measure its error.
The error-convergence data is shown in Figure 16.Since the above examples have different mesh scales, we have unified the mesh-scale axis in Figure 16 using the aboveused Δ in corresponding examples.For the solid-rocketmotor example, the Δ used for unifying is 6.The error data of the solid rocket motor example corresponding to 4 are missing because Δ = 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   decreases.However,   /Δ 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   /Δ 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.In order to quantify the convergence, we have performed Richardson's analysis [29] modeled by Eq. ( 24), where  is a constant and  is the unknown order of convergence.By substituting the error data into (24), the order of convergence  of each example problem can be solved via the leastsquares 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)

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.
( Taking equilateral triangles as an example, the magnitude of the resulting ∇ is 2 √ 3/3 ≈ 1.1547, which deviates from the true value (i.e., 1) by 15.47% no matter how small the scale of △ is.Although the above discussion is based on the 2D case, the 3D cell-based estimation method shows a similar problem near sharp features.

Figure 2 :
Figure 2: Main axes and domains of dependence of several target vertices.

Figure 4 :
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 ∇ − .△ 1  2 , △ 1  6 , and △ 5  6 overlap with the gray region.Including cells △ 1  2 and △ 5  6 in the calculation introduces information from  2 and  5 , which are not supposed to be included since they lie outside the gray region.Excluding △ 1  2 and △ 5  6 , however, results in a lack of robustness.

Figure 7 :
Figure 7: Computational mesh of the notched square.

Figure 8 :
Figure 8: Result of notched square example.

Figure 11 :
Figure 11: Working procedure of a solid rocket motor.

Figure 12 :
Figure 12: Computational mesh of a solid rocket motor.

Figure 15 :
Figure 15: Result of reconstructing the signed distance field.
1 |.  denotes the changing rate of .
Estimation.Assuming that ∇ is uniform within one differential volume, we know from the Conversely, in -dimensional space, we can use  non-linear correlative vectors and corresponding changingrate values to uniquely identify ∇, as shown in Figure