Research Article  Open Access
Robust ThreeDimensional LevelSet Method for Evolving Fronts on Complex Unstructured Meshes
Abstract
With a purpose to evolve the surfaces of complex geometries in their normal direction at arbitrarily defined velocities, we have developed a robust levelset approach which runs on threedimensional unstructured meshes. The approach is built on the basis of an innovative spatial discretization and corresponding gradientestimating 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 levelset method implementations on Cartesian meshes, the proposed approach makes it easier to describe jump boundary conditions and to perform coupling simulations.
1. Introduction
We are interested in evolving fronts on complex unstructured meshes in their normal direction using the levelset approach [1]. The levelset approach has been applied to various evolvingfront 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 levelset approach uses a higherdimensional function to model the evolving front, where the isosurface denotes the front. Since the implicit representation automatically handles topological changes, the levelset approach is capable of processing arbitrary geometry and evolving speed.
The levelset 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 HamiltonJacobi 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 nonCartesian 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 levelset functions are applied to address different boundaries, thus introducing more complexity. (b) Multiphase coupling simulation. Stress, fluid, and thermal coupling simulations are often required in addition to evolvingfront ones. Many such simulations are easier to perform on nonCartesian 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 ringshaped 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 levelset method to evolve fronts essentially requires solving the HamiltonJacobi equation formed like . Previous efforts to solve such equations on unstructured meshes have been focused on twodimensional (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 highorder 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 highorder numerical approximations for firstorder Hamilton Jacobi equations on triangular meshes on the base of the LaxFriedrichs scheme. Zhang [15] and Rokicki [16] then extended WENO reconstruction to 3D tetrahedral meshes.
Efforts to extend the levelset 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 levelset field and for evolving a front in the surface normal direction. Fu [18] mainly worked on applying parallel computing technology to speed up levelset 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., . The method, though appearing to be stable and convergent within smooth regions, produces an 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 upwind gradient estimation to fail to obtain the correct domain of dependence [19] p57. The problem of the upwind gradient is not significant when the levelset 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 upwind 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 abovementioned problems introduced by sharp features, we have developed an innovative spatial discretization scheme and a gradientestimating 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 LaxFriedrichs or the RoeFix schemes. The gradientestimating technology is more flexible than the divergence theorembased 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 levelset 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 levelset 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 Figure 1. Within each cell, the four vertices are denoted by , where should be attached to the vertices such that . In Figure 1, is the vertex for which we need to estimate . 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 . A vector which goes from the target vertex to any of the neighboring vertices, midpoints, or barycenters is referred as an outgoing vector to . is the set composed of all outgoing vectors that start from .
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 or , e.g., , . denotes the changing rate of .
3. Governing Equations
The levelset equation Eq.(1) describes the movement of the isosurface , where is the normal speed field at which the front evolves.
Usually, is not constant across the entire field, so would no longer be an exact signed distance field after each timeforward iteration step, and we have to reinitialize regularly. The reinitialization is performed via Eq.(2) [22], where is artificial time and ( is the cubic root of the volume of the smallest tetrahedral cell).
Equations (1) and (2) can be solved with respect to time via the Euler method or the total variation diminishing RungeKutta (TVDRK) 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 upwind 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 dimensionbydimension 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 nonCartesian meshes, Eq.(3) cannot be directly applied due to that the cell edges are not axisaligned. (b) In the vicinity of sharp features (e.g. 2D corners, 3D edges or corners), Eq.(3) may produce on however fine grid, leading Eq.(4) to output or 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 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 error. In Section 4.4, we have deduced the root cause of the 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 and . 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 Godunovtype scheme: when , either or is chosen to depend on which one has the maximum norm; when , is outputted as . Similarly, at valleys of the field, we have and . When , is outputted; when , 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), , and have the same signs. As a result, would be appropriately parallel to and , and , . According to Eq.(5)(6), the proposed scheme outputs when and when . Meanwhile, if and are both positive, points to the positive side of axis , thus the component of on is approximately equal to , and that of is approximately . Considering Eq.(3), which takes when and when , it can be seen that the proposed scheme takes the same domain of dependence and gives close result compared to the Godunovtype 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 Godunovtype 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 , , , . According to Eq.(3)(4), when , 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 and when , meaning that the proposed scheme is still in accordance with the Godunovtype 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 , 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 centraldifference 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 isosurface, since the computing is aimed at evolving the 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 isosurface, and the “explicit correction” technology is applied to obtain proper estimations for such vertices.
4.2. CellBased Gradient Estimation
Here, we briefly illustrate the problem of the cellbased approach before introducing our approach. The core of the cellbased estimating approach can be summarized by Eq. (8). The cellcentered estimations of obtained from (8) can then be integrated into various difference schemes, such as the upwind approach [17] and the WENO approach [15, 16]:
Such a cellbased approach, though used in a large portion of previous research, shows two major disadvantages when applied to featurerich geometries: (a) The orientations of cells with respect to the target vertex cannot be well defined, making it difficult to apply the abovepresented spatial discretization. The problem is illustrated in Figure 4. (b) As we have proved in the Appendix, cellbased estimation methods rely on a precondition 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 is continuous or smoother. Near sharp features (where the zero level set is continuous), the cellbased approach produces an error.
In this article, we have developed the edgebased 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. EdgeBased Gradient Estimation
Assuming that is uniform within one differential volume, we know from the Appendix that within this region, . Conversely, in dimensional space, we can use nonlinear correlative vectors and corresponding changingrate 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 edgebased approach.
(a) 2D
(b) 3D
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 cellbased 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 leastsquares 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 havewhere 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 be the target vertex and 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 leastsquares method,Substituting (13) into (14) we haveFrom (15), it is obvious that the proposed edgebased approach gives an exact estimation result for linear fields.
4.4. Gradient Estimation Near Sharp Features
While the edgebased method settles the problem of estimating upon the desired domain of dependence, the problem of the 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 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.
4.5. Overall Flow
The three core modules of our evolving approach, namely the spatial discretization, the edgebased 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 twostage RungeKutta method, which can be summarized by equation (17):
5. Validation
The following problems are tested to validate the accuracy and robustness of the proposed levelset 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.(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 nonuniform 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 mergingcircles 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, 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 Unotched square. 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 aswhere 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 smearingout 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 at the end of the simulation. The peaks in the two curves at steps are the product of the merging procedure of two sharp features near .
(a) Evolving result
(b) Evolving error
5.2. Three Merging Circles
This example includes three circles, the centers of which are located at , , and , respectively. The radiuses of the three circles are initially 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 at the end of the simulation.
(a) Evolving result
(b) Evolving error
5.3. Rate Stick
In the ratestick 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 highexplosives (HE) ratestick experiment described in [5]. Figure 10(a) demonstrates the mesh and initial status. The red point , which denotes a zeroradius circle such that , is the starting point of the front. across the entire region. The mesh scale used in the example is .
(a) Computational mesh
(b) Evolving result
(c) Evolving error
The slender computational region means that the boundaries must be handled smoothly to evade distortion. The zeroradius 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 ratestick problem. The peak in the 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 solidrocketmotor 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 hightemperature 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 trilinear interpolation across meshes. The mesh scale used to generate the reference result is . 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 levelset method on Cartesian grids.
(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),where is the halfsidelength 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 isosurface. It can be seen that the output field contains two desired properties, namely the sharp features and the nonuniform 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 .
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 isosurface is left unmoved. The error computed from the resulting field is , . 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.
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),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 errorconvergence data is shown in Figure 16. Since the above examples have different mesh scales, we have unified the meshscale axis in Figure 16 using the aboveused in corresponding examples. For the solidrocketmotor example, the used for unifying is 6. The error data of the solid rocket motor example corresponding to 4 are missing because 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 edgebased gradientestimating approach is proved to have firstorder accuracy, there must be other error sources that produce lowerorder error and consequently cause to increase as the mesh refines. Examples of such error sources include floatingpoint error, the remaining smearingout 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.
(a)
(b)
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 solidrocketmotor 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:

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 gradientestimating 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 gradientestimating 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 higherorder 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 higherorder schemes and further expand the range of application of the method.
Appendix
Error Analysis of CellBased Gradient Estimation
Taking the 2D cell in Figure 17(a) for example, the target vertex lies on a rectangular corner of the level set (denoted by blue). All edges of are within the outer sector of the corner, so the point on the front that is closest to is , , , . Vectors are outer normals of , where and . According to (8), we have
Letting , , is to the lefthand side of ; therefore, we have , and that the outer normals , , and . We also know that , . The triangle area . Substituting the above formulas into Eq. (A.1), we have The projection lengths of to and (as shown in Figure 17(b)) are
Taking equilateral triangles as an example, the magnitude of the resulting is , which deviates from the true value (i.e., 1) by no matter how small the scale of is. Although the above discussion is based on the 2D case, the 3D cellbased estimation method shows a similar problem near sharp features.
Data Availability
The source code and validation data used to support the findings of this study have been deposited in the GitHub repository https://github.com/metorm/UnstructuredLevelSetMethod.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
References
 S. Osher and J. A. Sethian, “Fronts propagating with curvaturedependent speed: algorithms based on HamiltonJacobi formulations,” Journal of Computational Physics, vol. 79, no. 1, pp. 12–49, 1988. View at: Publisher Site  Google Scholar  MathSciNet
 L. A. Vese and T. F. Chan, “A multiphase level set framework for image segmentation using the Mumford and Shah model,” International Journal of Computer Vision, vol. 50, no. 3, pp. 271–293, 2002. View at: Publisher Site  Google Scholar
 D.H. Wang, Y. Fei, F. Hu, and W.H. Zhang, “An integrated framework for solid rocket motor grain design optimization,” Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, vol. 228, no. 7, pp. 1156–1170, 2014. View at: Publisher Site  Google Scholar
 K. Albarado, A. Shelton, and R. J. Hartfield, “SRM simulation using the level set method and higher order integration schemes,” in Proceedings of the 48th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit 2012, USA, August 2012. View at: Google Scholar
 J. A. Zukas and W. P. Walters, Explosive Effects and Applications. Shock Wave and High Pressure Phenomena, Springer New York, New York, NY, 1998. View at: Publisher Site
 T. J. Barth and J. A. Sethian, “Numerical schemes for the HamiltonJacobi and level set equations on triangulated domains,” Journal of Computational Physics, vol. 145, no. 1, pp. 1–40, 1998. View at: Publisher Site  Google Scholar  MathSciNet
 T. D. Aslam, J. B. Bdzil, and D. S. Stewart, “Level set methods applied to modeling detonation shock dynamics,” Journal of Computational Physics, vol. 126, no. 2, pp. 390–409, 1996. View at: Publisher Site  Google Scholar  MathSciNet
 R. Abgrall, “Numerical discretization of the firstorder HamiltonJacobi equation on triangular meshes,” Communications on Pure and Applied Mathematics, vol. 49, no. 12, pp. 1339–1373, 1996. View at: Publisher Site  Google Scholar  MathSciNet
 X. Li, W. Yan, and C. K. Chan, “Numerical schemes for HamiltonJacobi equations on unstructured meshes,” Numerische Mathematik, vol. 94, no. 2, pp. 315–331, 2003. View at: Publisher Site  Google Scholar  MathSciNet
 R. Abgrall, “Numerical discretization of boundary conditions for first order HamiltonJacobi equations,” SIAM Journal on Numerical Analysis, vol. 41, no. 6, pp. 2233–2261, 2003. View at: Publisher Site  Google Scholar  MathSciNet
 Y.T. Zhang and C.W. Shu, “Highorder {WENO} schemes for HamiltonJacobi equations on triangular meshes,” SIAM Journal on Scientific Computing, vol. 24, no. 3, pp. 1005–1030, 2002. View at: Publisher Site  Google Scholar  MathSciNet
 C. Hu and C.W. Shu, “Weighted essentially nonoscillatory schemes on triangular meshes,” Journal of Computational Physics, vol. 150, no. 1, pp. 97–127, 1999. View at: Publisher Site  Google Scholar  MathSciNet
 J. Zhu and J. Qiu, “Hermite WENO schemes for HamiltonJacobi equations on unstructured meshes,” Journal of Computational Physics, vol. 254, pp. 76–92, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 S. Augoula and R. Abgrall, “High order numerical discretization for HamiltonJacobi equations on triangular meshes,” Journal of Scientific Computing, vol. 15, no. 2, pp. 197–229, 2000. View at: Publisher Site  Google Scholar  MathSciNet
 Y.T. Zhang and C.W. Shu, “Third order WENO scheme on three dimensional tetrahedral meshes,” Communications in Computational Physics, vol. 5, no. 24, pp. 836–848, 2009. View at: Google Scholar  MathSciNet
 J. Rokicki and R. Wieteska, “High Order Weno Schemes on Unstructured Tetrahedral Meshes,” in Proceedings of the In ECCOMAS CFD 2006: Proceedings of the European Conference on Computational Fluid Dynamics, p. 12, 2006. View at: Google Scholar
 N. R. Morgan and J. I. Waltz, “3D level set methods for evolving fronts on tetrahedral meshes with adaptive mesh refinement,” Journal of Computational Physics, vol. 336, pp. 492–512, 2017. View at: Publisher Site  Google Scholar  MathSciNet
 Z. Fu, S. Yakovlev, R. M. Kirby, and R. T. Whitaker, “Fast parallel solver for the levelset equations on unstructured meshes,” Concurrency Computation, vol. 27, no. 7, pp. 1639–1657, 2015. View at: Publisher Site  Google Scholar
 S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, vol. 153, Springer, New York, NY, USA, 2003. View at: Publisher Site  MathSciNet
 D. Hacon and C. Tomei, “Tetrahedral decompositions of hexahedral meshes,” European Journal of Combinatorics, vol. 10, no. 5, pp. 435–443, 1989. View at: Publisher Site  Google Scholar  MathSciNet
 S. K. Godunov, “A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics,” Matematicheskii Sbornik, vol. 47 (89), pp. 271–306, 1959. View at: Google Scholar  MathSciNet
 M. Sussman, P. Smereka, and S. Osher, “A level set approach for computing solutions to incompressible twophase flow,” Journal of Computational Physics, vol. 114, no. 1, pp. 146–159, 1994. View at: Publisher Site  Google Scholar
 O. Fortmeier and H. M. B\"Ucker, “Parallel reinitialization of level set functions on distributed unstructured tetrahedral grids,” Journal of Computational Physics, vol. 230, no. 12, pp. 4437–4453, 2011. View at: Publisher Site  Google Scholar  MathSciNet
 C. Min, “On reinitializing level set functions,” Journal of Computational Physics, vol. 229, no. 8, pp. 2764–2772, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 R. Wei, Unstructured Level Set Method, May 2018, https://github.com/metorm/UnstructuredLevelSetMethod.
 D. Engwirda, Locally Optimal DelaunayRefinement and OptimisationBased Mesh Generation [Ph.D. thesis], University of Sydney, November 2014.
 C. Geuzaine and J. F. Remacle, “Gmsh: A 3D finite element mesh generator with builtin pre and postprocessing facilities,” International Journal for Numerical Methods in Engineering, vol. 79, no. 11, pp. 1309–1331, 2009. View at: Publisher Site  Google Scholar  MathSciNet
 R. Wei, F. Bao, Y. Liu, and W. Hui, “Combined Acceleration Methods for Solid Rocket Motor Grain Burnback Simulation Based on the Level Set Method,” International Journal of Aerospace Engineering, vol. 2018, Article ID 4827810, 12 pages, 2018. View at: Publisher Site  Google Scholar
 E. Christiansen and H. G. Petersen, “Estimation of convergence orders in repeated Richardson extrapolation,” BIT Numerical Mathematics, vol. 29, no. 1, pp. 48–59, 1989. View at: Publisher Site  Google Scholar  MathSciNet
Copyright
Copyright © 2018 Ran Wei et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.