# Efficient Local Level Set Method without Reinitialization and Its Appliance to Topology Optimization

Authors: Wenhui Zhang; Yaoting Zhang
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

The local level set method (LLSM) is higher than the LSMs with global models in computational efficiency, because of the use of narrow-band model. The computational efficiency of the LLSM can be further increased by avoiding the reinitialization procedure by introducing a distance regularized equation (DRE). The numerical stability of the DRE can be ensured by a proposed conditionally stable difference scheme under reverse diffusion constraints. Nevertheless, the proposed method possesses no mechanism to nucleate new holes in the material domain for two-dimensional structures, so that a bidirectional evolutionary algorithm based on discrete level set functions is combined with the LLSM to replace the numerical process of hole nucleation. Numerical examples are given to show high computational efficiency and numerical stability of this algorithm for topology optimization.


---

## Body

## 1. Introduction

Topology optimization is a numerical iterative procedure for making an optimal layout of a structure or the best distribution of material in the conceptual design stage [1]. The level set method (LSM) is a recently developed approach to topology optimization that uses a flexible implicit description of the material domain [2]. The central idea of the LSM is to employ an implicit boundary describing model to parameterize the geometric model, and the boundary of a structure is embedded in a high-dimensional level set function that is called its zero level set [3]. The level set-based method is able to not only fundamentally avoid checkerboards and mesh-dependence, but also maintain smooth boundaries and distinct material interfaces during the topological design process [4]. Hence, many level set-based methods [5] have been developed for topology optimization since the LSM was first introduced into structure optimization.With an implicit local level set model, the computational efficiency of the local level set method (LLSM) [6] is much higher than that of the global level set methods, especially for shape optimization. However, the main shortcoming of the conventional LSM is that it possesses no mechanism to nucleate new holes in the material domain for two-dimensional structures, resulting in the final design heavily dependent on the initial guess [4]. A mechanism named the bubble-method [7] was first proposed to create new holes inside the structures in topology and shape optimization. This idea has been further developed into the mathematical concept of topological derivatives [8]. In the shape-sensitivity-based level set approaches, topological derivatives are incorporated to indicate the best place for introducing a new hole in a separate step of the optimization process [9] or as an additional term in the Hamilton-Jacobi equation [10]. The globally supported radial basis function (RBF) [11] and compactly supported RBF (CSRBF) [12] are typically used to discretize the original time-dependent initial value problem into an interpolation problem. The CSRBF brings about the strictly positive definiteness and sparseness properties of matrices under certain conditions. Hence the CSRBF has generalized the practical applications of RBFs to a larger set of scattered data [12].In the conventional LSM [3], a reinitialization procedure usually needs to reshape the level set function (LSF) to a signed distance function (SDF) periodically. However, the zero level set may drift away from its initial position by iteratively solving a classical reinitialization equation [13]. To suppress this drift, an interface preserving level set redistancing algorithm is proposed by Sussman and Fatemi [14]. Nevertheless, it has been proved that the SDF is not a feasible solution to the H-J equation [15]. In practice, it not only raises serious problems as when and how it should be performed, but also affects numerical accuracy in an undesirable way and thus should be avoided as much as possible [16]. The need for reinitialization was originally eliminated by introducing a penalty term [17] into a variational level set approach [18]. The undesirable boundary effect of the penalty term can be eliminated by taking a distance regularized equation (DRE) instead of this term. Hence, a so-called distance regularized level set evolution (DRLSE) [16] is realized based on the variational approach. As an unnecessary diffusion effect of the DRLSE was found in some locations where the surface is too flat, the DRE was recently modified with a new and balanced formulation to eliminate this effect [19]. Although parts of the diffusion rates in the DRE are negative, the numerical stability can still be maintained by incorporating reverse diffusion constraints in the difference schemes of the DRE, as can the reverse diffusion equations with all negative diffusion rates [20].The aim of this work is to solve the aforementioned numerical issues that still exist in the LLSM for topology optimization of two-dimensional structures. A bidirectional evolutionary algorithm based on the discrete level set functions (DLSFs) is proposed to find a stable topological solution first and then combined with the LLSM to further evolve the local details of the topology and shape of the structure. Transforming the DLSFs into the local level set function of the LLSM is achieved by iteratively solving the DRE. After that, the DRE is incorporated into the LLSM to avoid the reinitialization procedure. A difference scheme under reverse diffusion constraints is formulated for the DRE to improve its numerical stability. Typical examples are given to show the effectiveness of the proposed algorithm in terms of convergence, computational efficiency, and numerical stability.


## 2. Optimization Algorithm

### 2.1. Local Level Set Method Using Narrow-Band Model

In the local level set method (LLSM) [6], the local level set equation is defined as(1)∂ϕ∂t+CϕVn∇ϕ=0,where ϕx,t is defined as the local level set function (LLSF) and Vn is the normal velocity in normal direction n=∇ϕ/∇ϕ; the truncation function Cϕ is (2)Cϕ=1ifϕ≤Δϕ-γ22ϕ+γ-3Δγ-Δ3ifΔ<ϕ≤γ0ifϕ>γ,with T0=x:ϕx<γ being a narrow band with the half-band width γ. The narrow-band model and the corresponding LLSF are described as shown in Figure 1.Figure 1
Narrow-band model and local level set function in LLSM.It can be seen from Figure1 that only the LLSF within the narrow-band T0 needs to be updated during each iteration. Hence, the LLSM is higher in computational efficiency than the LSMs based on global level set models.


### 2.2. Bidirectional Evolutionary Algorithm with Discrete Level Set Functions

A two-dimensional structural model is built in the work regionD⊂R2. And the set S=S1∪S2∪S3 represents the finite elements in the domain D. It can be divided into three parts, S1 that consists of the solid elements with a full-material density; S2 that covers the elements with intermediate material densities; S3 that involves the void elements with a weak-material density. Accordingly, the nodal sets corresponding to the elemental sets S1,S2, and S3 are defined as S1n, S2n, and S3n, respectively. If it is assumed that S2n⊂S3n, S1n∪S3n consist of all the nodes within the region D, then a discrete level set function (DLSF) for node j can be defined as (3)ϕj=-c0j∈S1nc0j∈S3n,where c0 is a predefined constant set as 1 in this study.The values of elemental densities can be derived from the DLSFs. If theith element belongs to S1, that is, i∈S1, then the element density ρi=1; if i∈S3, then ρi=ρmin, where ρmin is a small value 0.001; if i∈S2, then ρi∈ρmin,1, where ρi is calculated in terms of the interpolation criterion given in the code manual [21]. In the structural model, the rectangular element is split into four triangles first, and the value of the DLSF at the common point of the triangles is then given by the average of the values of the four points. After that, each triangle is examined separately in the same logic. Finally, the elemental density is found to be the average of the contributions of the triangles.The structural stiffness design has been widely investigated in numerous literatures for topological sensitivity analysis. The standard notion [1] of minimum compliance design problems under a global volume constraint can be mathematically defined as follows:(4)minimizeϕJΩ=12uTKu,Subject  to∑i=1NViρiϕ=V∗,ρi∈ρmin,1,where J is known as the mean compliance, the open set Ω represents all admissible shapes in the design region D, uϕ is the nodal displacement vector, and K denotes the global stiffness matrix; Vi is the volume of an individual element, V∗ is the prescribed total volume, and N is the number of elements.Similar to the bubble-method [7] and the level set-based optimization methods [4, 10], topological derivatives are taken as topological sensitivities in this study. The topological derivative for node j is given by [9](5)αjn=DTJϕj=πλ+2μ4μλ+μ4μEjεu:εu+λ-μtrEjεutrεu,where Ej denotes the material elasticity tensor for node j, ε is the strain tensor, and the lame constants λ=E0ν/(1-ν2), μ=E0/21+ν with the Poisson ratio ν and Young’s modulus of solid materials E0. trA denotes the trace of a matrix A.Based on an interpolation function proposed by Shepard [22], a filter scheme of mesh independence is proposed to avoid the checkerboard patterns and mesh dependencies. A circular domain Ωw is first defined as the influence region centered round point x with cut-off radius rw, and NΩ denotes the number of points located inside the influence domain. The sensitivity filtering using the Shepard method with scattered points is then defined by(6)α~nx=∑i=1NΩWix·αnx,Wix=Dix∑j=1NDjxi=1,…,N,where Wix is the Shepard interpolation with the basis function Dix=1/dix2+c2, in which dix=x-xi denotes the radial distance from point x to xi, and if dix≥rw, then Dix is set to zero. c is a positive constant and chosen as a onefold mesh size in terms of numerical experiences.Over the last two decades, many topology description models have been developed for topology optimization of structures, which can roughly be classified into two categories, the material distribution model and the boundary description model [23]. Based on the material distribution model, the ESO (Evolutionary Structural Optimization) method has won a great deal of popularity in recent years [24]. The bidirectional ESO (BESO) method [25], as an extension of the ESO method, allows efficient material to be added to the structure while the inefficient one is removed simultaneously. So a bidirectional evolutionary algorithm is developed by integrating both the DLSFs and topological derivatives into the optimization criteria of the BESO method [25]. Note that the design variables and topological sensitivities in the BESO method are based on the elemental pseudo densities while those in the proposed algorithm are based on the discrete level set functions.It is assumed that the volume in thekth iteration Vk is known, and k≥0. The target volume Vk+1 in the next iteration is then updated by(7)Vk+1=min⁡Vk1+ER,V∗ifVk≤V∗max⁡Vk1-ER,V∗ifVk>V∗with the evolutionary volume ratio ER and the volume limit V∗ defined in (4).A parameterARn is defined as the adding number of nodes in the set S1n divided by the total numbers of nodes and ARn≤ARmaxn, where ARmaxn is a predefined positive constant. The definition of ARn is different from that of AR in the original BESO method, since the former parameter corresponds to nodal sensitivity while the latter one corresponds to elemental sensitivity.It is assumed that the DLSFϕjk of node j is known in the kth iteration. Then the DLSF ϕjk+1 in the next iteration is updated by (8)ϕjk+1=-c0αjn≥αthadd,ϕjk=c0c0αjn≥αthdel,ϕjk=-c0,where the threshold sensitivity numbers αthdel and αthadd are determined as the number of nodes decreased from the set S1n and that increased from the set S3n, respectively. These thresholds are similar to those based on elemental sensitivities given in the original BESO method. Full details of determining these thresholds are described in [25].Finally, a stable topological solution is obtained when the following convergence criterion is satisfied:(9)∑k-9k-5JΩ-∑k-4kJΩ∑k-4kJΩ≤τ,where τ is an allowable convergence error with typical values ranging from 0.001 to 0.01.The majority of logical steps of the bidirectional evolutionary algorithm are presented in Figure2.Figure 2
Flow chart depicting logical steps of the bidirectional evolutionary algorithm.


### 2.3. Distance Regularized Equation (DRE) and Its Improvement

In the distance regularized level set evolution (DRLSE) [16], the DRE can retain the signed distance feature ∇ϕ=1 at least within the narrow-band region near boundaries without reinitialization, whose formula is expressed in the standard form of the diffusion equation as(10)∂ϕ∂t=μdiv⁡α1ϕ∇ϕ,with the diffusion rate α1ϕ=μdp1∇ϕ, where the diffusion function is set to dp1s=p1′s/s with s=∇ϕ.In the original DRLSE, the energy densityp1s was defined as (11)p1s=1-cos2πs2π2ifs≤1s-122ifs>1,which is a double-well potential function because there are two minimum points of p1s at s=1 and s=0. So the diffusion function dp1s is given by (12)dp1s=sin2πs2πsifs≤11-1sifs>1.It is easy to verify the boundedness of the diffusion rateα1ϕ=μdp1∇ϕ and α1ϕ≤μ. It can be seen from (12) that, for ∇ϕ>1, dp1∇ϕ is positive, and ∇ϕ will decrease and approach 1; for 0.5<∇ϕ≤1, dp1∇ϕ is negative, and ∇ϕ will increase and approach 1; for ∇ϕ≤0.5, dp1∇ϕ is positive, and ∇ϕ will decrease and approach 0.If∇ϕ0≤0.5 is satisfied for all the initial values ϕ0, the diffusion effect of (10) will make ∇ϕ approach 0. So it loses the ability to regularize ∇ϕ to 1. An improved diffusion rate α2ϕ with a diffusion function like the following is proposed in [19]:(13)dp2∇ϕ=2πarctan⁡∇ϕ-1/σϕ≤ε2πarctan⁡∇ϕ-1/σϕ>ε,where σ is a positive constant and is chosen as fourfold mesh sizes in terms of numerical experiences. T1=x:ϕx<ε is a narrow band with a half-band width ε.The diffusion effect can be divided into two parts: the forward diffusion for∇ϕ≥1 and the backward diffusion for ∇ϕ<1. It will make ∇ϕ approach one within T1 but zero outside T1. However, the two parts are balanced within T1 but unbalanced outside T1 so that multiple iterations are required to retain a flatter level set surface outside T1. In this paper, the diffusion function dp2∇ϕ is further localized by introducing the half-band width γ of the narrow-band T0 in LLSM, thereby resulting in an improved diffusion rate α3ϕ using the diffusion function(14)dp3∇ϕ=2πarctan⁡∇ϕ-1/σϕ<γ0ϕ≥γ.With the diffusion rateα3ϕ, the two parts are balanced within T0 without influencing the level set surface outside T0.


### 2.4. A Conditionally Stable Difference Scheme for DRE

It is noted that the common difference schemes for the DRE with parts of the negative diffusion rates are incapable of remaining stable during an iterative process, according to the stability definition of the difference equation. In our numerical experiments,ϕ is apt to gradually become divergent along with the process of iterations. To enhance the numerical stability of the DRE, a difference scheme similar to that of the mean curvature given in [18] is developed and described as(15)ϕi,jk+1-ϕi,jkΔt=αϕki+1/2,jδ+xϕi,j-αϕki-1/2,jδ-xϕi,jΔx+αϕki,j+1/2δ+yϕi,j-αϕki,j-1/2δ-yϕi,jΔy,where(16)δ+xϕi,j=ϕi+1,jk-ϕi,jkΔx,δcxϕi,j=ϕi+1,jk-ϕi-1,jk2Δx,δ-xϕi,j=ϕi,jk-ϕi-1,jkΔx,δ+yϕi,j=ϕi,j+1k-ϕi,jkΔy,δcyϕi,j=ϕi,j+1k-ϕi,j-1k2Δy,δ-yϕi,j=ϕi,jk-ϕi,j-1kΔy,and αϕ≔α3ϕ, in which the difference schemes of ∇ϕi±1/2,j,∇ϕi,j±1/2 are given by(17)∇ϕi±1/2,j=δ±xϕi,j2+δcyϕi,j+δcyϕi±1,j22,∇ϕi,j±1/2=δ±yϕi,j2+δcxϕi,j+δcxϕi±1,j22.It has been verified by our numerical experiments that the evolution of level set can remain bounded stability even after a large number of iterations for solving (15). However, the maximum of ϕ often exceeds the initial value c0 defined by (1), which leads to the level set surface unsmoothed near the edges of the narrow-band T0, thereby reducing the computational accuracy of (15). Furthermore, multiple iterations are required to find a suitable Courant-Friedrichs-Lewy (CFL) condition to ensure the stability of (15). The issues related to the numerical instability of the DRE can be resolved by imposing reverse diffusion constraints on the difference scheme (see (15)), since it has been proved that the constraints can ensure the numerical stability of the diffusion equations with all negative diffusion rates [20].First, (15) can be subdivided along the direction of x and y into(18)ϕi,jk+1/2-ϕi,jk=ΔtΔxϕi+1,jk-ϕi,jkα¯i+1/2,j+ϕi-1,jk-ϕi,jkα¯i-1/2,j,ϕi,jk+1-ϕi,jk+1/2=ΔtΔyϕi,j+1k-ϕi,jkα¯i,j+1/2+ϕi,j-1k-ϕi,jkα¯i,j-1/2,with α¯i±1/2,j=1/Δxαϕki+1/2,j and α¯i,j±1/2=1/Δyαϕki,j±1/2.Then the four flow functions are defined as(19)Fi,jx=ϕi+1,jk-ϕi,jkα¯i+1/2,j,Fi-1,jx=ϕi-1,jk-ϕi,jkα¯i-1/2,j,Fi,jy=ϕi,j+1k-ϕi,jkα¯i,j+1/2,Fi-1,jy=ϕi,j-1k-ϕi,jkα¯i,j-1/2,where Fi,jx denotes the change from ϕi,jk to ϕi+1,jk in one time step Δt and in the x direction, and the definitions of Fi,jx, Fi,jy, and Fi,j-1y are similar to that of Fi,jx. The lowest and highest limit values of these flow functions are defined as(20)Flow=minp,q=-1,0,1⁡ϕi+p,j+qk-ϕi,jk,Fhigh=maxp,q=-1,0,1⁡ϕi+p,j+qk-ϕi,jk.It can be seen that the four diffusion rates in (15) satisfy the boundedness αϕk≤μ. If μ≤Δx, the reverse diffusion constraints can be defined by (21)Flow≤Fi,jx,Fi-1,jx,Fi,jy,Fi,j-1y≤Fhigh.IfΔt/Δx≤1/4, first substituting (20) into inequalities (21) and then substituting the result into (18), one can obtain a solution using the inequalities(22)minp,q=-1,0,1⁡ϕi+p,j+qk≤ϕi,jk+1≤maxp,q=-1,0,1⁡ϕi+p,j+qk.It can be seen that the absolute values ofϕi+p,j+qk for p,q=-1,0,1 are lower than their initial value c0. That means that all the absolute values ϕ≤c0 if the CFL conditions are satisfied:(23)μΔtΔx2≤14,μΔtΔy2≤14.


### 2.5. Flow Chart for Difference Schemes to LLSM with DRE

The procedure for the LLSM with the DRE consists of two main parts, transforming the models of discrete level set functions into the local level set function and solving the difference schemes of the LLSE associated with the DRE. The final DLSFs corresponding to the stable topological solution can be transformed into the LLSF within the initial narrow-bandT0 by iteratively solving the DRE under reverse diffusion constraints. The LLSE can be solved by difference schemes using the third-order Runge-Kutta (R-K) scheme for temporal discretization and the fifth-order weighted essentially nonoscillatory (WENO) scheme for spatial discretization. The reader is referred to [26] for more numerical details. The logical steps of the two parts can be described by a flow chart given in Figure 3.Figure 3
Flow chart depicting logical steps of the LLSM.


### 2.6. Shape Derivatives and Normal Extension Velocities

With the classical level set model [4], the minimum compliance design problem given in (4) can be converted to an unconstrained problem with the Lagrangian method:(24)MinimizeϕLu,ϕ=Ju,ϕ+λ+∫DHϕdΩ-V∗,where the Lagrangian function Lu,ϕ is the objective functional. λ+ is the Lagrangian multiplier of the volume constraint. Hϕ is the Heaviside function. The mean compliance Ju,ϕ is reformulated as(25)Ju,ϕ=12∫DEεu:εuHϕdΩ.For a number of level set-based approaches [4, 9, 11, 12], the steepest descent method is used to ensure the decrease of the objective function by directly setting the normal velocity field Vn as the negative shape derivative of Lu,ϕ. For the particular case of a 2-D model of a linear elastic structure, the boundary traction is fixed and remains unchanged, the displacement constraint is fixed, and the body force is set to zero; thus Vn can be given by(26)Vn=0.5Eεu:εu-λ+.The reader is referred to the article [4] for more detailed theoretical proofs. In addition, a bisectioning algorithm is used to find the Lagrangian multiplier λ+ to guarantee that the volume constraint be exactly satisfied during each iteration.The normal velocity field can be naturally extended to the whole domain using the so-called “ersatz material” approach, which fills the void areas with a weak material and then the material density is assumed to be piecewise constant in each element and is adequately interpolated in those elements cut by the zero level set (the shape boundary) [4]. In the LLSM, the extension velocity field is localized within the narrow-band T0. By iteratively solving the difference scheme for the DRE (see (15)), one can obtain a smooth velocity field in the region near the edges of the narrow-band T0 to improve the computational accuracy of the extension velocity.


## 2.1. Local Level Set Method Using Narrow-Band Model

In the local level set method (LLSM) [6], the local level set equation is defined as(1)∂ϕ∂t+CϕVn∇ϕ=0,where ϕx,t is defined as the local level set function (LLSF) and Vn is the normal velocity in normal direction n=∇ϕ/∇ϕ; the truncation function Cϕ is (2)Cϕ=1ifϕ≤Δϕ-γ22ϕ+γ-3Δγ-Δ3ifΔ<ϕ≤γ0ifϕ>γ,with T0=x:ϕx<γ being a narrow band with the half-band width γ. The narrow-band model and the corresponding LLSF are described as shown in Figure 1.Figure 1
Narrow-band model and local level set function in LLSM.It can be seen from Figure1 that only the LLSF within the narrow-band T0 needs to be updated during each iteration. Hence, the LLSM is higher in computational efficiency than the LSMs based on global level set models.


## 2.2. Bidirectional Evolutionary Algorithm with Discrete Level Set Functions

A two-dimensional structural model is built in the work regionD⊂R2. And the set S=S1∪S2∪S3 represents the finite elements in the domain D. It can be divided into three parts, S1 that consists of the solid elements with a full-material density; S2 that covers the elements with intermediate material densities; S3 that involves the void elements with a weak-material density. Accordingly, the nodal sets corresponding to the elemental sets S1,S2, and S3 are defined as S1n, S2n, and S3n, respectively. If it is assumed that S2n⊂S3n, S1n∪S3n consist of all the nodes within the region D, then a discrete level set function (DLSF) for node j can be defined as (3)ϕj=-c0j∈S1nc0j∈S3n,where c0 is a predefined constant set as 1 in this study.The values of elemental densities can be derived from the DLSFs. If theith element belongs to S1, that is, i∈S1, then the element density ρi=1; if i∈S3, then ρi=ρmin, where ρmin is a small value 0.001; if i∈S2, then ρi∈ρmin,1, where ρi is calculated in terms of the interpolation criterion given in the code manual [21]. In the structural model, the rectangular element is split into four triangles first, and the value of the DLSF at the common point of the triangles is then given by the average of the values of the four points. After that, each triangle is examined separately in the same logic. Finally, the elemental density is found to be the average of the contributions of the triangles.The structural stiffness design has been widely investigated in numerous literatures for topological sensitivity analysis. The standard notion [1] of minimum compliance design problems under a global volume constraint can be mathematically defined as follows:(4)minimizeϕJΩ=12uTKu,Subject  to∑i=1NViρiϕ=V∗,ρi∈ρmin,1,where J is known as the mean compliance, the open set Ω represents all admissible shapes in the design region D, uϕ is the nodal displacement vector, and K denotes the global stiffness matrix; Vi is the volume of an individual element, V∗ is the prescribed total volume, and N is the number of elements.Similar to the bubble-method [7] and the level set-based optimization methods [4, 10], topological derivatives are taken as topological sensitivities in this study. The topological derivative for node j is given by [9](5)αjn=DTJϕj=πλ+2μ4μλ+μ4μEjεu:εu+λ-μtrEjεutrεu,where Ej denotes the material elasticity tensor for node j, ε is the strain tensor, and the lame constants λ=E0ν/(1-ν2), μ=E0/21+ν with the Poisson ratio ν and Young’s modulus of solid materials E0. trA denotes the trace of a matrix A.Based on an interpolation function proposed by Shepard [22], a filter scheme of mesh independence is proposed to avoid the checkerboard patterns and mesh dependencies. A circular domain Ωw is first defined as the influence region centered round point x with cut-off radius rw, and NΩ denotes the number of points located inside the influence domain. The sensitivity filtering using the Shepard method with scattered points is then defined by(6)α~nx=∑i=1NΩWix·αnx,Wix=Dix∑j=1NDjxi=1,…,N,where Wix is the Shepard interpolation with the basis function Dix=1/dix2+c2, in which dix=x-xi denotes the radial distance from point x to xi, and if dix≥rw, then Dix is set to zero. c is a positive constant and chosen as a onefold mesh size in terms of numerical experiences.Over the last two decades, many topology description models have been developed for topology optimization of structures, which can roughly be classified into two categories, the material distribution model and the boundary description model [23]. Based on the material distribution model, the ESO (Evolutionary Structural Optimization) method has won a great deal of popularity in recent years [24]. The bidirectional ESO (BESO) method [25], as an extension of the ESO method, allows efficient material to be added to the structure while the inefficient one is removed simultaneously. So a bidirectional evolutionary algorithm is developed by integrating both the DLSFs and topological derivatives into the optimization criteria of the BESO method [25]. Note that the design variables and topological sensitivities in the BESO method are based on the elemental pseudo densities while those in the proposed algorithm are based on the discrete level set functions.It is assumed that the volume in thekth iteration Vk is known, and k≥0. The target volume Vk+1 in the next iteration is then updated by(7)Vk+1=min⁡Vk1+ER,V∗ifVk≤V∗max⁡Vk1-ER,V∗ifVk>V∗with the evolutionary volume ratio ER and the volume limit V∗ defined in (4).A parameterARn is defined as the adding number of nodes in the set S1n divided by the total numbers of nodes and ARn≤ARmaxn, where ARmaxn is a predefined positive constant. The definition of ARn is different from that of AR in the original BESO method, since the former parameter corresponds to nodal sensitivity while the latter one corresponds to elemental sensitivity.It is assumed that the DLSFϕjk of node j is known in the kth iteration. Then the DLSF ϕjk+1 in the next iteration is updated by (8)ϕjk+1=-c0αjn≥αthadd,ϕjk=c0c0αjn≥αthdel,ϕjk=-c0,where the threshold sensitivity numbers αthdel and αthadd are determined as the number of nodes decreased from the set S1n and that increased from the set S3n, respectively. These thresholds are similar to those based on elemental sensitivities given in the original BESO method. Full details of determining these thresholds are described in [25].Finally, a stable topological solution is obtained when the following convergence criterion is satisfied:(9)∑k-9k-5JΩ-∑k-4kJΩ∑k-4kJΩ≤τ,where τ is an allowable convergence error with typical values ranging from 0.001 to 0.01.The majority of logical steps of the bidirectional evolutionary algorithm are presented in Figure2.Figure 2
Flow chart depicting logical steps of the bidirectional evolutionary algorithm.


## 2.3. Distance Regularized Equation (DRE) and Its Improvement

In the distance regularized level set evolution (DRLSE) [16], the DRE can retain the signed distance feature ∇ϕ=1 at least within the narrow-band region near boundaries without reinitialization, whose formula is expressed in the standard form of the diffusion equation as(10)∂ϕ∂t=μdiv⁡α1ϕ∇ϕ,with the diffusion rate α1ϕ=μdp1∇ϕ, where the diffusion function is set to dp1s=p1′s/s with s=∇ϕ.In the original DRLSE, the energy densityp1s was defined as (11)p1s=1-cos2πs2π2ifs≤1s-122ifs>1,which is a double-well potential function because there are two minimum points of p1s at s=1 and s=0. So the diffusion function dp1s is given by (12)dp1s=sin2πs2πsifs≤11-1sifs>1.It is easy to verify the boundedness of the diffusion rateα1ϕ=μdp1∇ϕ and α1ϕ≤μ. It can be seen from (12) that, for ∇ϕ>1, dp1∇ϕ is positive, and ∇ϕ will decrease and approach 1; for 0.5<∇ϕ≤1, dp1∇ϕ is negative, and ∇ϕ will increase and approach 1; for ∇ϕ≤0.5, dp1∇ϕ is positive, and ∇ϕ will decrease and approach 0.If∇ϕ0≤0.5 is satisfied for all the initial values ϕ0, the diffusion effect of (10) will make ∇ϕ approach 0. So it loses the ability to regularize ∇ϕ to 1. An improved diffusion rate α2ϕ with a diffusion function like the following is proposed in [19]:(13)dp2∇ϕ=2πarctan⁡∇ϕ-1/σϕ≤ε2πarctan⁡∇ϕ-1/σϕ>ε,where σ is a positive constant and is chosen as fourfold mesh sizes in terms of numerical experiences. T1=x:ϕx<ε is a narrow band with a half-band width ε.The diffusion effect can be divided into two parts: the forward diffusion for∇ϕ≥1 and the backward diffusion for ∇ϕ<1. It will make ∇ϕ approach one within T1 but zero outside T1. However, the two parts are balanced within T1 but unbalanced outside T1 so that multiple iterations are required to retain a flatter level set surface outside T1. In this paper, the diffusion function dp2∇ϕ is further localized by introducing the half-band width γ of the narrow-band T0 in LLSM, thereby resulting in an improved diffusion rate α3ϕ using the diffusion function(14)dp3∇ϕ=2πarctan⁡∇ϕ-1/σϕ<γ0ϕ≥γ.With the diffusion rateα3ϕ, the two parts are balanced within T0 without influencing the level set surface outside T0.


## 2.4. A Conditionally Stable Difference Scheme for DRE

It is noted that the common difference schemes for the DRE with parts of the negative diffusion rates are incapable of remaining stable during an iterative process, according to the stability definition of the difference equation. In our numerical experiments,ϕ is apt to gradually become divergent along with the process of iterations. To enhance the numerical stability of the DRE, a difference scheme similar to that of the mean curvature given in [18] is developed and described as(15)ϕi,jk+1-ϕi,jkΔt=αϕki+1/2,jδ+xϕi,j-αϕki-1/2,jδ-xϕi,jΔx+αϕki,j+1/2δ+yϕi,j-αϕki,j-1/2δ-yϕi,jΔy,where(16)δ+xϕi,j=ϕi+1,jk-ϕi,jkΔx,δcxϕi,j=ϕi+1,jk-ϕi-1,jk2Δx,δ-xϕi,j=ϕi,jk-ϕi-1,jkΔx,δ+yϕi,j=ϕi,j+1k-ϕi,jkΔy,δcyϕi,j=ϕi,j+1k-ϕi,j-1k2Δy,δ-yϕi,j=ϕi,jk-ϕi,j-1kΔy,and αϕ≔α3ϕ, in which the difference schemes of ∇ϕi±1/2,j,∇ϕi,j±1/2 are given by(17)∇ϕi±1/2,j=δ±xϕi,j2+δcyϕi,j+δcyϕi±1,j22,∇ϕi,j±1/2=δ±yϕi,j2+δcxϕi,j+δcxϕi±1,j22.It has been verified by our numerical experiments that the evolution of level set can remain bounded stability even after a large number of iterations for solving (15). However, the maximum of ϕ often exceeds the initial value c0 defined by (1), which leads to the level set surface unsmoothed near the edges of the narrow-band T0, thereby reducing the computational accuracy of (15). Furthermore, multiple iterations are required to find a suitable Courant-Friedrichs-Lewy (CFL) condition to ensure the stability of (15). The issues related to the numerical instability of the DRE can be resolved by imposing reverse diffusion constraints on the difference scheme (see (15)), since it has been proved that the constraints can ensure the numerical stability of the diffusion equations with all negative diffusion rates [20].First, (15) can be subdivided along the direction of x and y into(18)ϕi,jk+1/2-ϕi,jk=ΔtΔxϕi+1,jk-ϕi,jkα¯i+1/2,j+ϕi-1,jk-ϕi,jkα¯i-1/2,j,ϕi,jk+1-ϕi,jk+1/2=ΔtΔyϕi,j+1k-ϕi,jkα¯i,j+1/2+ϕi,j-1k-ϕi,jkα¯i,j-1/2,with α¯i±1/2,j=1/Δxαϕki+1/2,j and α¯i,j±1/2=1/Δyαϕki,j±1/2.Then the four flow functions are defined as(19)Fi,jx=ϕi+1,jk-ϕi,jkα¯i+1/2,j,Fi-1,jx=ϕi-1,jk-ϕi,jkα¯i-1/2,j,Fi,jy=ϕi,j+1k-ϕi,jkα¯i,j+1/2,Fi-1,jy=ϕi,j-1k-ϕi,jkα¯i,j-1/2,where Fi,jx denotes the change from ϕi,jk to ϕi+1,jk in one time step Δt and in the x direction, and the definitions of Fi,jx, Fi,jy, and Fi,j-1y are similar to that of Fi,jx. The lowest and highest limit values of these flow functions are defined as(20)Flow=minp,q=-1,0,1⁡ϕi+p,j+qk-ϕi,jk,Fhigh=maxp,q=-1,0,1⁡ϕi+p,j+qk-ϕi,jk.It can be seen that the four diffusion rates in (15) satisfy the boundedness αϕk≤μ. If μ≤Δx, the reverse diffusion constraints can be defined by (21)Flow≤Fi,jx,Fi-1,jx,Fi,jy,Fi,j-1y≤Fhigh.IfΔt/Δx≤1/4, first substituting (20) into inequalities (21) and then substituting the result into (18), one can obtain a solution using the inequalities(22)minp,q=-1,0,1⁡ϕi+p,j+qk≤ϕi,jk+1≤maxp,q=-1,0,1⁡ϕi+p,j+qk.It can be seen that the absolute values ofϕi+p,j+qk for p,q=-1,0,1 are lower than their initial value c0. That means that all the absolute values ϕ≤c0 if the CFL conditions are satisfied:(23)μΔtΔx2≤14,μΔtΔy2≤14.


## 2.5. Flow Chart for Difference Schemes to LLSM with DRE

The procedure for the LLSM with the DRE consists of two main parts, transforming the models of discrete level set functions into the local level set function and solving the difference schemes of the LLSE associated with the DRE. The final DLSFs corresponding to the stable topological solution can be transformed into the LLSF within the initial narrow-bandT0 by iteratively solving the DRE under reverse diffusion constraints. The LLSE can be solved by difference schemes using the third-order Runge-Kutta (R-K) scheme for temporal discretization and the fifth-order weighted essentially nonoscillatory (WENO) scheme for spatial discretization. The reader is referred to [26] for more numerical details. The logical steps of the two parts can be described by a flow chart given in Figure 3.Figure 3
Flow chart depicting logical steps of the LLSM.


## 2.6. Shape Derivatives and Normal Extension Velocities

With the classical level set model [4], the minimum compliance design problem given in (4) can be converted to an unconstrained problem with the Lagrangian method:(24)MinimizeϕLu,ϕ=Ju,ϕ+λ+∫DHϕdΩ-V∗,where the Lagrangian function Lu,ϕ is the objective functional. λ+ is the Lagrangian multiplier of the volume constraint. Hϕ is the Heaviside function. The mean compliance Ju,ϕ is reformulated as(25)Ju,ϕ=12∫DEεu:εuHϕdΩ.For a number of level set-based approaches [4, 9, 11, 12], the steepest descent method is used to ensure the decrease of the objective function by directly setting the normal velocity field Vn as the negative shape derivative of Lu,ϕ. For the particular case of a 2-D model of a linear elastic structure, the boundary traction is fixed and remains unchanged, the displacement constraint is fixed, and the body force is set to zero; thus Vn can be given by(26)Vn=0.5Eεu:εu-λ+.The reader is referred to the article [4] for more detailed theoretical proofs. In addition, a bisectioning algorithm is used to find the Lagrangian multiplier λ+ to guarantee that the volume constraint be exactly satisfied during each iteration.The normal velocity field can be naturally extended to the whole domain using the so-called “ersatz material” approach, which fills the void areas with a weak material and then the material density is assumed to be piecewise constant in each element and is adequately interpolated in those elements cut by the zero level set (the shape boundary) [4]. In the LLSM, the extension velocity field is localized within the narrow-band T0. By iteratively solving the difference scheme for the DRE (see (15)), one can obtain a smooth velocity field in the region near the edges of the narrow-band T0 to improve the computational accuracy of the extension velocity.


## 3. Numerical Examples

In this section, two widely researched examples, the cantilever beam and the arch bridge, are presented in the context of structural minimum compliance design to demonstrate the characteristics of the proposed method. Some of the system parameters using the same values are defined as follows.Young’s elasticity modulus for the solid material isE0=200GPa and for weak material is Emin=10-3Pa, and the Poisson ratio is 0.3. The volume constraint V∗=0.5Vall, where Vall is the total volume in the design region D. The convergence tolerance τ is set as 0.01 in the bidirectional evolutionary algorithm and 0.001 in the LLSM.


### 3.1. First Cantilever-Beam Model

Shown in Figure4 is the design domain of a cantilever beam with a size 40mm×25mm. The left side of the domain is fixed as the Dirichlet boundary, and a concentrated force P=100N is vertically applied at the central point of the right side as a nonhomogeneous Neumann boundary. In the bidirectional evolutionary algorithm, ER=2%, ARmaxn=5%, and rw=2mm. In the difference schemes for the LLSE, Δt=0.1, Δ=0.3, and γ=0.999. In the difference schemes for the DRE, μ=0.5 and dt=0.001. The design domain is discretized with a mesh of 80×50 quadrilateral elements.Figure 4
Design domain of the first cantilever beam and its boundary conditions.In the design domain as shown in Figure4, the initial volume V0 of the solid region Ω is set to V0=Vall. The structural topologies corresponding to the zero level set and related level set surfaces are shown in Figures 5 and 6, respectively. The convergence histories of the mean compliance and volume fraction are depicted in Figure 7. The result in Figure 5(e) stands for a stable topological solution obtained from the bidirectional evolutionary algorithm. Topological results given by this algorithm are characterized by a smooth boundary attributed to the structural model described by the DLSFs. By comparing Figures 6(e) and 6(f), the sharp level set surface corresponding to the DLSFs has been successfully converted into a smooth one related to a local level set function by iteratively solving the DRE at the initial stage of the LLSM. Then the shape of the boundary is further improved by iteratively solving the LLSE. In addition, all the absolute values of the level set function are less than the initial value c0. Therefore it verifies the effectiveness of reverse diffusion constraints on the numerical stability of the difference scheme for the DRE.Figure 5
Topologies of zero level set at different steps: (a) Step1; (b) Step 10; (c) Step 20; (d) Step 30; (e) Step 43; (f) Step 57.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
The corresponding level set surface at different steps: (a) Step1; (b) Step 10; (c) Step 20; (d) Step 30; (e) Step 43; (f) Step 57.
(a)
(b)
(c)
(d)
(e)
(f)Figure 7
Convergent histories of the objective function and the constraint.


### 3.2. Second Cantilever-Beam Model

This study has also investigated the influence of different initial models of structure on the final design. Figure8 depicts the design domain with a size 4.0mm×2.5mm. The left side of the domain is fixed and a concentrated force P=1N is vertically applied at the bottom of the right side. All the parameters but rw=0.2mm and ARmaxn=1% remain unchanged as those of the first cantilever-beam model. Figure 9 shows two cases of the initial configurations with full materials and the least but essential materials and their topologies during the process of optimization. The two final designs are made with the same topology and almost similar shape of the structure, which shows the complexity of the final topology is not changing appreciably with different initial structures. Therefore, the numerical process of the bidirectional evolutionary algorithm can be used to replace the numerical process of hole nucleation in the LLSM to avoid the final design heavily dependent on the initial guess.Figure 8
Design domain of the second cantilever beam and its boundary conditions.Figure 9
Topologies of zero level set in the two cases of the second cantilever beam.
(a)
The case starts from the full-material initial configurations
(b)
The case starts from the least-material initial configurationsFigure10 shows the topological topologies for almost rw=0.2mm using several mesh sizes. It can be seen that the optimal topology does not depend on the discretization in terms of layout and number of bars.Figure 10
Mesh-independent solutions of the second cantilever beam: (a)40×25; (b) 96×60; (c) 128×80; (d) 160×100.
(a)
(b)
(c)
(d)


### 3.3. Arch-Bridge Model

The design domain of an arch-bridge model with a size2.0mm×1.2mm is shown in Figure 11. Both the bottom corners of the domain are the fixed support. A uniform static pressure is vertically applied on the upper side, and the sum of the pressure is 5 N. In accordance with the same definitions of the cantilever-beam parameters, the parameters are set to ER=3%, ARmaxn=5%, rw=0.1mm, Δt=0.02, Δ=0.3, γ=0.99, dt=0.001, and μ=0.0167. The design domain is discretized with a mesh of 120×72 quadrilateral elements.Figure 11
Design domain of the arch bridge.This example focuses on the new characteristic of the proposed algorithm for improving the convergence of the bidirectional evolutionary algorithm using the LLSM. The structural topologies corresponding to the zero level set are depicted in Figure12. Note that it starts from the initial model with the volume V=0.5Vall and remains unchanged, so as to maintain the stability of the evolution process of the object function. The evolutionary histories of the objective and the volume constraint starting from the initial models are plotted in Figure 13. A design of the structure shown in Figure 12(b) corresponding to a stable topological solution is also the final design of the topology (not shape) made merely using the bidirectional evolutionary algorithm in our numerical experiments. The subsequent topologies given in Figures 12(c)–12(f) show that the LLSM can further optimize the topology of the structure to improve convergence. Moreover, the LLSM can also improve the shape of the boundary to obtain a smoother shape design till it reaches the convergence tolerance in the 38th iteration. It is worth noticing that the final topology obtained by the LLSM is just a local optimal design because of the use of the steepest descent method. Despite the optimal solution of this arch-bridge model obtained by an element-wise ESO method [27], in this case the optimized topology obtained by the proposed bidirectional evolutionary algorithm is not reasonable compared with the final optimal solution shown in Figure 12(f).Figure 12
Topologies of zero level set at different steps starting from the case ofV=0.5Vall: (a) Initial; (b) Step 19; (c) Step 20; (d) Step 21; (e) Step 22; (f) Step 38.
(a)
(b)
(c)
(d)
(e)
(f)Figure 13
Convergent histories of the objective function and the constraint.Starting from different initial models, the final models obtained by the bidirectional evolutionary algorithm and the LLSM, respectively, are shown in Figure14.Figure 14
Topologies of zero level set for the initial model, and the final models obtained by the bidirectional evolutionary algorithm and the LLSM, respectively: (a) case 1,V=0.4Vall; (b) case 2, V=0.6Vall; (c) case 3, V=0.8Vall.
(a)
(b)
(c)It can be seen from Figures12 and 14 that the final topologies obtained by the bidirectional evolutionary algorithm are inconsistent. In contrast, the final optimized results found using the LLSM subsequently are of the same topology and similar shape. Although the local optimal solution of the arch-bridge model can also be obtained by using the ESO/BESO methods with elemental variables, numerical instabilities and zigzag boundaries can result in these elemental variables-based methods. Hence, nodal variables are needed to take the place of the elemental variables in these methods. The combined algorithm with the bidirectional evolutionary algorithm and the LLSM can also resolve this problem and achieve at least the consistent local optimal solution for the different cases of initial models.


## 3.1. First Cantilever-Beam Model

Shown in Figure4 is the design domain of a cantilever beam with a size 40mm×25mm. The left side of the domain is fixed as the Dirichlet boundary, and a concentrated force P=100N is vertically applied at the central point of the right side as a nonhomogeneous Neumann boundary. In the bidirectional evolutionary algorithm, ER=2%, ARmaxn=5%, and rw=2mm. In the difference schemes for the LLSE, Δt=0.1, Δ=0.3, and γ=0.999. In the difference schemes for the DRE, μ=0.5 and dt=0.001. The design domain is discretized with a mesh of 80×50 quadrilateral elements.Figure 4
Design domain of the first cantilever beam and its boundary conditions.In the design domain as shown in Figure4, the initial volume V0 of the solid region Ω is set to V0=Vall. The structural topologies corresponding to the zero level set and related level set surfaces are shown in Figures 5 and 6, respectively. The convergence histories of the mean compliance and volume fraction are depicted in Figure 7. The result in Figure 5(e) stands for a stable topological solution obtained from the bidirectional evolutionary algorithm. Topological results given by this algorithm are characterized by a smooth boundary attributed to the structural model described by the DLSFs. By comparing Figures 6(e) and 6(f), the sharp level set surface corresponding to the DLSFs has been successfully converted into a smooth one related to a local level set function by iteratively solving the DRE at the initial stage of the LLSM. Then the shape of the boundary is further improved by iteratively solving the LLSE. In addition, all the absolute values of the level set function are less than the initial value c0. Therefore it verifies the effectiveness of reverse diffusion constraints on the numerical stability of the difference scheme for the DRE.Figure 5
Topologies of zero level set at different steps: (a) Step1; (b) Step 10; (c) Step 20; (d) Step 30; (e) Step 43; (f) Step 57.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
The corresponding level set surface at different steps: (a) Step1; (b) Step 10; (c) Step 20; (d) Step 30; (e) Step 43; (f) Step 57.
(a)
(b)
(c)
(d)
(e)
(f)Figure 7
Convergent histories of the objective function and the constraint.


## 3.2. Second Cantilever-Beam Model

This study has also investigated the influence of different initial models of structure on the final design. Figure8 depicts the design domain with a size 4.0mm×2.5mm. The left side of the domain is fixed and a concentrated force P=1N is vertically applied at the bottom of the right side. All the parameters but rw=0.2mm and ARmaxn=1% remain unchanged as those of the first cantilever-beam model. Figure 9 shows two cases of the initial configurations with full materials and the least but essential materials and their topologies during the process of optimization. The two final designs are made with the same topology and almost similar shape of the structure, which shows the complexity of the final topology is not changing appreciably with different initial structures. Therefore, the numerical process of the bidirectional evolutionary algorithm can be used to replace the numerical process of hole nucleation in the LLSM to avoid the final design heavily dependent on the initial guess.Figure 8
Design domain of the second cantilever beam and its boundary conditions.Figure 9
Topologies of zero level set in the two cases of the second cantilever beam.
(a)
The case starts from the full-material initial configurations
(b)
The case starts from the least-material initial configurationsFigure10 shows the topological topologies for almost rw=0.2mm using several mesh sizes. It can be seen that the optimal topology does not depend on the discretization in terms of layout and number of bars.Figure 10
Mesh-independent solutions of the second cantilever beam: (a)40×25; (b) 96×60; (c) 128×80; (d) 160×100.
(a)
(b)
(c)
(d)


## 3.3. Arch-Bridge Model

The design domain of an arch-bridge model with a size2.0mm×1.2mm is shown in Figure 11. Both the bottom corners of the domain are the fixed support. A uniform static pressure is vertically applied on the upper side, and the sum of the pressure is 5 N. In accordance with the same definitions of the cantilever-beam parameters, the parameters are set to ER=3%, ARmaxn=5%, rw=0.1mm, Δt=0.02, Δ=0.3, γ=0.99, dt=0.001, and μ=0.0167. The design domain is discretized with a mesh of 120×72 quadrilateral elements.Figure 11
Design domain of the arch bridge.This example focuses on the new characteristic of the proposed algorithm for improving the convergence of the bidirectional evolutionary algorithm using the LLSM. The structural topologies corresponding to the zero level set are depicted in Figure12. Note that it starts from the initial model with the volume V=0.5Vall and remains unchanged, so as to maintain the stability of the evolution process of the object function. The evolutionary histories of the objective and the volume constraint starting from the initial models are plotted in Figure 13. A design of the structure shown in Figure 12(b) corresponding to a stable topological solution is also the final design of the topology (not shape) made merely using the bidirectional evolutionary algorithm in our numerical experiments. The subsequent topologies given in Figures 12(c)–12(f) show that the LLSM can further optimize the topology of the structure to improve convergence. Moreover, the LLSM can also improve the shape of the boundary to obtain a smoother shape design till it reaches the convergence tolerance in the 38th iteration. It is worth noticing that the final topology obtained by the LLSM is just a local optimal design because of the use of the steepest descent method. Despite the optimal solution of this arch-bridge model obtained by an element-wise ESO method [27], in this case the optimized topology obtained by the proposed bidirectional evolutionary algorithm is not reasonable compared with the final optimal solution shown in Figure 12(f).Figure 12
Topologies of zero level set at different steps starting from the case ofV=0.5Vall: (a) Initial; (b) Step 19; (c) Step 20; (d) Step 21; (e) Step 22; (f) Step 38.
(a)
(b)
(c)
(d)
(e)
(f)Figure 13
Convergent histories of the objective function and the constraint.Starting from different initial models, the final models obtained by the bidirectional evolutionary algorithm and the LLSM, respectively, are shown in Figure14.Figure 14
Topologies of zero level set for the initial model, and the final models obtained by the bidirectional evolutionary algorithm and the LLSM, respectively: (a) case 1,V=0.4Vall; (b) case 2, V=0.6Vall; (c) case 3, V=0.8Vall.
(a)
(b)
(c)It can be seen from Figures12 and 14 that the final topologies obtained by the bidirectional evolutionary algorithm are inconsistent. In contrast, the final optimized results found using the LLSM subsequently are of the same topology and similar shape. Although the local optimal solution of the arch-bridge model can also be obtained by using the ESO/BESO methods with elemental variables, numerical instabilities and zigzag boundaries can result in these elemental variables-based methods. Hence, nodal variables are needed to take the place of the elemental variables in these methods. The combined algorithm with the bidirectional evolutionary algorithm and the LLSM can also resolve this problem and achieve at least the consistent local optimal solution for the different cases of initial models.


## 4. Conclusions

The LLSM is intended to remarkably increase the computational efficiency of the conventional LSMs using global models. To overcome the issue of hole nucleation of the LLSM, a bidirectional evolutionary algorithm is combined with the LLSM. This proposed algorithm has been used successfully in topology optimization of two-dimensional (2-D) structures, and it is easy to be extended to 3-D structures. The main features of this algorithm unknown to the conventional LSMs and the LLSM can be summarized as follows:(a)
The discrete level set functions can be efficiently transformed into the local level set function by iteratively solving the distance regularized equation (DRE).(b)
The DRE can be used instead of the reinitialization equation to further increase the computational efficiency of the LLSM.(c)
A conditionally stable difference scheme under reverse diffusion constraints is formulated to ensure the numerical stability of the DRE.(d)
If the stable topological solutions of the bidirectional evolutionary algorithm are inconsistent, the LLSM can achieve at least the consistent local optimal solution for the different cases of initial models.High computational efficiency and numerical stability of the proposed algorithm have been verified by three typical numerical examples.



---

*Source: 6392901-2016-01-14.xml*