LINEBURG

φo units from the interface. In t units of time the interface will approach

interface as it moves normal to itself in the outward direction.

a t spatial units closer, changing the value of this point to φo ’a t, which

When φ is a signed distance function, equation (6.1) reduces to φt = ’a,

is exactly the forward Euler time update of this point. The exact interface

and the values of φ either increase or decrease, depending on the sign of a.

crossing time can be identi¬ed for all points by solving φo ’ at = 0 to get

Forward Euler time discretization of this equation gives φn+1 = φn ’ a t.

t = φo /a. (Similar arguments hold when a < 0, except that the interface

When a > 0, the φ = 0 isocontour becomes the φ = ’a t isocontour

moves in the opposite direction.)

after one time step. Similarly, the φ = a t isocontour becomes the φ = 0

Here, we see the power of signed distance functions. When φo is a signed

isocontour. That is, the φ = 0 isocontour moves a t units forward in

distance function, we can write down the exact solution of equation (6.1)

the normal direction to the old position of the old φ = a t isocontour.

as φ(t) = φo ’ at. On the other hand, when φo is not a signed distance

The interface is moving in the normal direction with speed a. Taking the

function, equation (6.1) needs to be solved numerically by treating it as a

gradient of this forward Euler time stepping gives φn+1 = φn ’ (a t).

Hamilton-Jacobi equation, as discussed in the last chapter.

Since a t is spatially constant, (a t) = 0, implying that φn+1 = φn .

6.2. Numerical Discretization 57 58 6. Motion in the Normal Direction

6.2 Numerical Discretization viscosity solution. The RF scheme treats the ambiguity associated with up-

winding near sonic points by using central di¬erencing plus some arti¬cial

For instructive purposes, suppose we plug V = aN into equation (4.2) and viscosity.

try a simple upwind di¬erencing approach. That is, we will attempt to Recall that numerical dissipation can smear out ¬ne solution details

discretize on coarse grids. In order to avoid as much numerical smearing as pos-

sible, we have proposed ¬ve di¬erent versions (LF, SLF, LLF, SLLF, and

aφx aφy aφz

·

φt + , , φ=0 (6.2) LLLF) of the central di¬erencing plus arti¬cial viscosity approach to solving

| φ| | φ| | φ|

Hamilton-Jacobi problems. While the RF method is a better alternative,

with simple upwinding. Consider the ¬rst spatial term aφx | φ|’1 φx , where it too resorts to arti¬cial dissipation in the vicinity of sonic points where

aφx | φ|’1 is the “velocity” in the x-direction. Since upwinding is based ambiguities occur. In order to avoid the addition of arti¬cial dissipation,

only on the sign of the velocity, we can ignore the positive | φ| denomina- one must resort to the Godunov scheme.

tor, assuming temporarily that it is nonzero. Then the sign of aφx can be Let us examine the Godunov method in detail. Again, assume a > 0. If

used to decide whether φ’ or φ+ should be used to approximate φx . When φx and φ+ are both positive, then extx minimizes H when φ’ < φ+ and

’

x x x x x

φ’ and φ+ have the same sign, it does not matter which of these is plugged maximizes H when φ’ > φ+ . In either case, we choose φ’ consistent with

x x x x x

into aφx , since only the sign of this term determines whether we use φ’ or upwinding. Similarly, when φ’ and φ+ are both negative, extx minimizes

x x x

φ+ . For example, suppose a > 0. Then when φ’ > 0 and φ+ > 0, aφx > 0 H when φ’ < φ+ and maximizes H when φ’ > φ+ . Again, φ+ is chosen

x x x x x x x x

and φ’ should be used in equation (6.2) everywhere φx appears, including in both instances consistent with upwinding. Now consider the “V”-shaped

x

the velocity term. On the other hand, when φ’ < 0 and φ+ < 0, aφx < 0 case where φ’ < 0 and φ+ > 0, indicating an expansion. Here φ’ < φ+ ,

x x x x x x

and φ+ should be used to approximate φx . so Godunov™s method minimizes H, achieving this minimum by setting

x

This simple upwinding approach works well as long as φ’ and φ+ have φx = 0. This implies that a region of expansion should have a locally ¬‚at φ

x x

the same sign, but consider what happens when they have di¬erent signs. with φx = 0. Instead of adding numerical dissipation to hide the problem,

For example, when φ’ < 0 and φ+ > 0, aφ’ < 0 (still assuming a > 0), Godunov™s method chooses the most meaningful solution. Next, consider

x x x

indicating that φ+ should be used, while aφ+ > 0, indicating that φ’ the case where φ’ > 0 and φ+ < 0, indicating coalescing characteristics.

x x x x x

’ +

should be used. This situation corresponds to a “V”-shaped region where Here φx > φx , so Godunov™s method maximizes H, achieving this max-

imum by setting φx equal to the larger in magnitude of φ’ and φ+ . In

each side of the “V” should move outward. The di¬culty in approximating x x

φx arises because we are in the vicinity of a sonic point, where φx = 0. The this shock case, information is coming from both directions, and the grid

LLLF interval de¬ned by φ’ and φ+ includes this sonic point since φ’ and point feels the e¬ects of the information that gets there ¬rst. The velocities

x x x

+

are characterized by H1 = aφx | φ|’1 , and the side with the fastest speed

φx di¬er in sign. We have to take care to ensure that the expansion takes

place properly. A similar problem occurs when φ’ > 0 and φ+ < 0. Here arrives ¬rst. This is determined by taking the larger in magnitude of φ’

x x x

’ ’ +

and φ+ . Again, Godunov™s method chooses the most meaningful solution,

aφx > 0, indicating that φx should be used, while aφx < 0, indicating that x

φ+ should be used. This upside-down “V” is shaped like a carrot (or hat) avoiding arti¬cial dissipation.

x

and represents the coalescing of information similar to a shock wave. Once Godunov™s method for equation (6.1) can be summarized as follows for

both positive and negative a. If aφ’ and aφ+ are both positive, use φx =

again caution is needed to ensure that the correct solution is obtained. x x

Simple upwinding breaks down when φ’ and φ+ di¬er in sign. Let φ’ . If aφ’ and aφ+ are both negative, use φx = φ+ . If aφ’ ¤ 0 and

x x x x x x x

aφ+ ≥ 0, treat the expansion by setting φx = 0. If aφ’ ≥ 0 and aφ+ ¤ 0,

us examine how the Roe-Fix method works in this case. In order to do x x x

treat the shock by setting φx to either φ’ or φ+ , depending on which gives

this, we need to consider the Hamilton-Jacobi form of the equation, i.e., x x

equation (6.1). Here H1 = aφx | φ|’1 , implying that the simple velocity the largest magnitude for aφx . Note that when φ’ = φ+ = 0 both of the

x x

V = aN we used in equation (6.2) was the correct expression to look at last two cases are activated, and both consistently give φx = 0. We also

for upwinding. The sign of H1 is independent of the y and z directions, have the following elegant formula due to Rouy and Tourin [139]:

depending only on aφx . If both φ’ and φ+ have the same sign, we choose

φ2 ≈ max max(φ’ , 0)2 , min(φ+ , 0)2

x x

(6.3)

x x x

one or the other depending on the sign of H1 as in the usual upwinding.

However, unlike simple upwinding, Roe-Fix gives a consistent method for when a > 0, and

treating the case where φ’ and φ+ di¬er in sign. In that instance we are in

φ2 ≈ max min(φ’ , 0)2 , max(φ+ , 0)2

x x

(6.4)

x x x

the vicinity of a sonic point and switch to the LLF method, adding some

numerical dissipation to the scheme in order to obtain the correct vanishing when a < 0.

6.3. Adding a Curvature-Dependent Term 59 60 6. Motion in the Normal Direction

6.3 Adding a Curvature-Dependent Term 1 1 1

0.5

0.5 0.5

Most ¬‚ames burn with a speed in the normal direction plus extra heating

0

0 0

and cooling e¬ects due to the curvature of the front. This velocity ¬eld can

be modeled by setting Vn = a ’ bκ in the level set equation (4.4) to obtain ’0.5 ’0.5 ’0.5

’1 ’1 ’1

φt + a| φ| = bκ| φ|; (6.5) 1

0

’1 ’1 1

0 ’1 0 1

which has both hyperbolic and parabolic terms. The hyperbolic a| φ| term 1 1 1

can be discretized as outlined above using Hamilton-Jacobi techniques, 0.5 0.5 0.5

while the parabolic bκ| φ| term can be independently discretized using

0 0 0

central di¬erencing as described in Chapter 4.

’0.5 ’0.5

Once both terms have been discretized, either forward Euler or RK ’0.5

time discretization can be used to advance the front forward in time. ’1 ’1 ’1

’1 1

0 ’1 0 1 0

’1 1

The combined CFL condition for equations that contain both hyperbolic

Hamilton-Jacobi terms and parabolic terms is given by 1 1 1

0.5 0.5 0.5

|H1 | |H2 | |H3 | 2b 2b 2b

t + + + 2+ 2+ < 1, (6.6) 0 0

2 0

x y z ( x) ( y) ( z)

’0.5 ’0.5 ’0.5

as one might have guessed from equation (4.12).

’1 ’1 ’1

0 1

’1 1

’1 0 0

’1 1

Figure 6.2. A star-shaped interface being advected by a rigid body rotation as it

6.4 Adding an External Velocity Field moves outward normal to itself.

Equation (6.5) models a ¬‚ame front burning through a material at rest, but

does not account for the velocity of the unburnt material. A more general background velocity and the self-generated velocity in the normal direction

equation is are both moving the front in the same direction, and the upwind direction is

easy to determine. For example, if a, φ’ , and φ+ are all positive, the second

φt + V · φ + a| φ| = bκ| φ|, (6.7) x x

term in H1 indicates that the self-generated normal velocity is moving the

since it includes the velocity V of the unburnt material. This equation front to the right. Additionally, when u > 0 the external velocity is also

combines an external velocity ¬eld with motion normal to the interface and moving the front to the right. In this case, both the RF scheme and the

motion by mean curvature. It is the most general form of the G-equation Godunov scheme set φx = φ’ . x

for burning ¬‚ame fronts; see Markstein [110]. As in equation (6.5), the It is more di¬cult to determine what is happening when u and aφx | φ|’1

parabolic term on the right-hand side can be independently discretized with disagree in sign. In this case the background velocity is moving the front in

central di¬erencing. The hyperbolic Hamilton-Jacobi part of this equation one direction while the self-generated normal velocity is moving the front

consists of two terms, V · φ and a| φ|. Figure 6.2 shows the evolution in the opposite direction. In order to upwind, we must determine which of

of a star-shaped interface under the in¬‚uence of both an externally given these two terms dominates. It helps if φ is a signed distance function, since

rigid-body rotation (a V · φ term) and a self-generated motion outward we obtain the simpli¬ed H1 = u + aφx . If H1 is positive for both φ’ and x

φ+ , then both RF and Godunov set φx = φ’ . If H1 is negative for both

normal to the interface (an a| φ| term). x x

φ’ and φ+ , then both RF and Godunov set φx = φ+ . If H1 is negative

In order to discretize the Hamilton-Jacobi part of equation (6.7), we ¬rst x x x

identify the partial derivatives of H, i.e., H1 = u + aφx | φ|’1 and H2 = for φ’ and positive for φ+ , we have an expansion. If φ’ < φ+ , Godunov™s

x x x x

v + aφy | φ|’1 . The ¬rst term in H1 represents motion of the interface as method chooses the minimum value for H. This relative extremum occurs

when H1 = 0 implying that we set φx = ’u/a. If φ’ > φ+ Godunov™s

it is passively advected in the external velocity ¬eld, while the second term x x

represents the self-generated velocity of the interface as it moves normal method chooses the maximum value for H, which is again obtained by

to itself. If u and aφx have the same sign for both φ’ and φ+ , then the setting φx = ’u/a. If H1 is positive for φ’ and negative for φ+ , we have a

x x x x

6.4. Adding an External Velocity Field 61

shock. Godunov™s method dictates setting φx equal to the value of φ’ or x

φ+ that gives the value of H1 with the largest magnitude.

x

When φ is not a signed distance function, the above simpli¬cations can-

not be made. In general, H1 = u + aφx | φ|’1 , and we need to consider not

only I x , but also the values of φy and φz in I y and I z , respectively. This

can become rather complicated quite quickly. In fact, even the RF method

can become quite complicated in this case, since it is hard to tell when sonic

7

points are nearby and when they are not. In situations like this, the LLF

scheme is ideal, since one merely uses both the values of φ’ and φ+ along

x x

with some arti¬cial dissipation setting ± as dictated by equation (5.12).

Constructing Signed Distance

At ¬rst glance, equation (5.12) might seem complicated to evaluate; e.g.,

Functions

one has to determine the maximum value of |H1 |. However, since ± is just

a dissipation coe¬cient, we can safely overestimate ± and pay the price

of slightly more arti¬cial dissipation. In contrast, it is hard to predict how

certain approximations will a¬ect the Godunov scheme. One way to approx-

imate ± is to separate H1 into parts, i.e., using |H1 | < |u| + |aφx || φ|’1

to treat the ¬rst and second terms independently. Also, when φ is approxi-

mately a signed distance, we can look at |H1 | = |u + aφx |. This function is

easy to maximize, since the maximum occurs at either φ’ or φ+ and the y

x x

and z spatial directions play no role.

7.1 Introduction

As we have seen, a number of simpli¬cations can be made when φ is a signed

distance function. For this reason, we dedicate this chapter to numerical

techniques for constructing approximate signed distance functions. These

techniques can be applied to the initial data in order to initialize φ to a

signed distance function.

As the interface evolves, φ will generally drift away from its initialized

value as signed distance. Thus, the techniques presented in this chapter

need to be applied periodically in order to keep φ approximately equal to

signed distance. For a particular application, one has to decide how sensi-

tive the relevant techniques are to φ™s approximation of a signed distance

function. If they are very sensitive, φ needs to be reinitialized to signed

distance both accurately and often. If they are not sensitive, one can reini-

tialize with a lower-order accurate method on an occasional basis. However,

even if a particular numerical approach doesn™t seem to depend on how ac-

curately φ approximates a signed distance function, one needs to remember

that φ can develop noisy features and steep gradients that are not amenable

to ¬nite-di¬erence approximations. For this reason, it is always advisable

to reinitialize occasionally so that φ stays smooth enough to approximate

its spatial derivatives with some degree of accuracy.

64 7. Constructing Signed Distance Functions 7.3. Crossing Times 65

7.2 Reinitialization creases ¬‚exibility in algorithm design, since di¬culties away from the φ = 0

isocontour can be ignored.

In their seminal level set paper, Osher and Sethian [126] initialized their

numerical calculations with φ = 1±d2 , where d is the distance function and

the “±” sign is negative in „¦’ and positive in „¦+ . Later, it became clear

7.3 Crossing Times

that the signed distance function φ = ±d, was a better choice for initializ-

ing φ. Mulder, Osher, and Sethian [115] demonstrated that initializing φ to

One of the di¬culties associated with the direct computation of signed

a signed distance function results in more accurate numerical solutions than

distance functions is locating and discretizing the interface. This can be

initializing φ to a Heaviside function. While it is obvious that better results

avoided in the following fashion. Consider a point x ∈ „¦+ . If x does not lie

can be obtained with smooth functions than nonsmooth functions, there

on the interface, we wish to know how far from the interface it is so that

are those who insist on using (slightly smeared out) Heaviside functions,

we can set φ(x) = +d. If we move the interface in the normal direction

or color functions, to track interfaces.

using equation (6.1) with a = 1, the interface eventually crosses over x,

In [48], Chopp considered an application where certain regions of the

changing the local value of φ from positive to negative. If we keep a time

¬‚ow had level sets piling up on each other, increasing the local gradient,

history of the local values of φ at x, we can ¬nd the exact time when φ was

while other regions of the ¬‚ow had level sets separating from each other,

equal to zero using interpolation in time. This is the time it takes the zero

¬‚attening out φ. In order to reduce the numerical errors caused by both

level set to reach the point x, and we call that time to the crossing time.

steepening and ¬‚attening e¬ects, [48] introduced the notion that one should

Since equation (6.1) moves the level set normal to itself with speed a = 1,

reinitialize the level set function periodically throughout the calculation.

the time it takes for the zero level set to reach a point x is equal to the

Since only the φ = 0 isocontour has any meaning, one can stop the cal-

distance the interface is from x. That is, the crossing time to is equal to the

culation at any point in time and reset the other isocontours so that φ is

distance d. For points x ∈ „¦’ , the crossing time is similarly determined

again initialized to a signed distance function. The most straightforward

using a = ’1 in equation (6.1).

way of implementing this is to use a contour plotting algorithm to locate

In a series of papers, [20], [97], and [100], Kimmel and Bruckstein intro-

and discretize the φ = 0 isocontour and then explicitly measure distances

duced the notion of using crossing times in image-processing applications.

from it. Unfortunately, this straightforward reinitialization routine can be

For example, [100] used equation (6.1) with a = 1 to create shape o¬sets,

slow, especially if it needs to be done at every time step. In order to ob-

which are distance functions with distance measured from the boundary of

tain reasonable run times, [48] restricted the calculations of the interface

an image. The idea of using crossing times to solve some general Hamilton-

motion and the reinitialization to a small band of points near the φ = 0

Jacobi equations with Dirichlet boundary conditions was later generalized

isocontour, producing the ¬rst version of the local level set method. We

and rigorized by Osher [123].

refer those interested in local level set methods to the more recent works

of Adalsteinsson and Sethian [2] and Peng, Merriman, Osher, Zhao, and

Kang [130].

7.4 The Reinitialization Equation

The concept of frequent reinitialization is a powerful numerical tool. In

a standard numerical method, one starts with initial data and proceeds

In [139], Rouy and Tourin proposed a numerical method for solving | φ| =

forward in time, assuming that the numerical solution stays well behaved

f (x) for a spatially varying function f derived from the intensity of an

until the ¬nal solution is computed. With reinitialization, we have a less-

image. In the trivial case of f (x) = 1, the solution φ is a signed distance

stringent assumption, since only our φ = 0 isocontour needs to stay well

function. They added f (x) to the right-hand side of equation (6.1) as a

behaved. Any problems that creep up elsewhere are wiped out when the

source term to obtain

level set is reinitialized. For example, Merriman, Bence, and Osher [114]

proposed numerical techniques that destroy the nice properties of the level

φt + | φ| = f (x), (7.1)

set function and show that poor numerical solutions are obtained using

these degraded level set functions. Then they show that periodic reini- which is evolved in time until a steady state is reached. At steady state,

tialization to a signed distance function repairs the damage, producing the values of φ cease to change, implying that φt = 0. Then equation (7.1)

reduces to | φ| = f (x), as desired. Since only the steady-state solution

high-quality numerical results. Numerical techniques need to be e¬ective

only for the φ = 0 isocontour, since the rest of the implicit surface can be is desired, [139] used an accelerated iteration method instead of directly

repaired by reinitializing φ to a signed distance function. This greatly in- evolving equation (7.1) forward in time.

66 7. Constructing Signed Distance Functions 7.4. The Reinitialization Equation 67

Equation (7.1) propagates information in the normal direction, so infor- does not change if the interface incorrectly crosses over a grid point. This

mation ¬‚ows from smaller values of φ to larger values of φ. This equation was addressed directly by Fedkiw, Aslam, Merriman, and Osher in the

is of little use in reinitializing the level set function, since the interface lo- appendix of [63], where incorrect interface crossings were identi¬ed as sign

cation will be in¬‚uenced by the negative values of φ. That is, the φ = 0 changes in the nodal values of φ. These incorrect interface crossings were

isocontour is not guaranteed to stay ¬xed, but will instead move around recti¬ed by putting the interface back on the correct side of a grid point x,

setting φ(x) = ± , where ± is a small number with the appropriate sign.

as it is in¬‚uenced by the information ¬‚owing in from the negative values

of φ. One way to avoid this is to compute the signed distance function for In discretizing equation (7.4), the S(φo )| φ| term is treated as motion in

all the grid points adjacent to the interface by hand. Then the normal direction as described in Chapter 6. Here S(φo ) is constant for

all time and can be thought of as a spatially varying “a” term. Numerical

φt + | φ| = 1 (7.2)

tests indicate that better results are obtained when S(φo ) is numerically

can be solved in „¦+ to update φ based on those grid points adjacent to smeared out, so [160] used

the interface. That is, the grid points adjacent to the interface are not

φo

updated, but instead used as boundary conditions. Since there is only a S(φo ) = (7.5)

φ2 + ( x)2

single band of initialized grid cells on each side of the interface, one cannot o

apply higher-order accurate methods such as HJ WENO. However, if a

as a numerical approximation. Later, Peng, Merriman, Osher, Zhao, and

two-grid-cell-thick band is initialized in „¦’ (in addition to the one-grid-

Kang [130] suggested that

cell-thick band in „¦+ ), the total size if the band consists of three grid

cells and the HJ WENO scheme can then be used. Alternatively, one could φ

S(φ) = (7.6)

initialize a three-grid-cell-thick band of boundary conditions in „¦’ and φ2 + | φ|2 ( x)2

use these to update every point in „¦+ including those adjacent to the

interface. Similarly, a three grid cell thick band of boundary conditions can was a better choice, especially when the initial φo was a poor estimate of

be initialized in „¦+ and used to update the values of φ in „¦’ by solving signed distance, i.e., when | φo | was far from 1. In equation (7.6), it is im-

portant to update S(φ) continually as the calculation progresses so that the

φt ’ | φ| = ’1 (7.3)

| φ| term has the intended e¬ect. In contrast, equation (7.5) is evaluated

to steady state. Equations (7.2) and (7.3) reach steady state rather quickly, only once using the initial data. Numerical smearing of the sign function

since they propagate information at speed 1 in the direction normal to the decreases its magnitude, slowing the propagation speed of information near

interface. For example, if t = 0.5 x, it takes only about 10 time steps to the interface. This probably aids the balancing out of the circular depen-

move information from the interface to 5 grid cells away from the interface. dence on the initial data as well, since it produces characteristics that do

In [160], Sussman, Smereka, and Osher put all this together into a not look as far across the interface for their information. We recommend

reinitialization equation using Godunov™s method for discretizing the hyperbolic S(φo )| φ| term.

After ¬nding a numerical approximation to S(φo )| φ|, we combine it with

φt + S(φo )(| φ| ’ 1) = 0, (7.4)

the remaining S(φo ) source term at each grid point and update the resulting

where S(φo ) is a sign function taken as 1 in „¦+ , ’1 in „¦’ , and 0 on the quantity in time with a Runge-Kutta method.

interface, where we want φ to stay identically equal to zero. Using this Ideally, the interface remains stationary during the reinitialization pro-

equation, there is no need to initialize any points near the interface for use cedure, but numerical errors will tend to move it to some degree. In [158],

as boundary conditions. The points near the interface in „¦+ use the points Sussman and Fatemi suggested an improvement to the standard reini-

in „¦’ as boundary conditions, while the points in „¦’ conversely look at tialization procedure. Since their application of interest was two-phase

those in „¦+ . This circular loop of dependencies eventually balances out, incompressible ¬‚ow, they focused on preserving the amount of material

and a steady-state signed distance function is obtained. As long as φ is in each cell, i.e., preserving the area (volume) in two (three) spatial di-

relatively smooth and the initial data are somewhat balanced across the mensions. If the interface does not move during reinitialization, the area is

interface, this method works rather well. Unfortunately, if φ is not smooth preserved. On the other hand, one can preserve the area while allowing the

or φ is much steeper on one side of the interface than the other, circular interface to move, implying that their proposed constraint is weaker than

dependencies on initial data can cause the interface to move incorrectly it should be. In [158] this constraint was applied locally, requiring that the

from its initial starting position. For this reason, [160] de¬ned S(φo ) using area be preserved individually in each cell. Instead of using the exact area,

the initial values of φ (denoted by φo ) so that the domain of dependence the authors used equation (1.15) with f (x) = 1 to approximate the area in

68 7. Constructing Signed Distance Functions 7.5. The Fast Marching Method 69

each cell as procedure is very similar to that used by Rudin, Osher, and Fatemi [142]

as a continuous in time gradient projection method. In [142] a set of con-

Ai,j = H(φ) dx, (7.7) straints needs to be preserved under evolution, while in [158] the evolution

„¦i,j

is not inherited from gradient descent on a functional to be optimized.

where „¦i,j is an individual grid cell and H is the smeared-out Heaviside Reinitialization is still an active area of research. Recently, Russo and

function in equation (1.22). In both [158] and the related [159] by Sussman, Smereka [143] introduced yet another method for computing the signed

Fatemi, Smereka, and Osher this local constraint was shown to signi¬cantly distance function. This method was designed to keep the stencil from in-

improve the results obtained with the HJ ENO scheme. However, this local correctly looking across the interface at values that should not in¬‚uence it,

constraint method has not yet been shown to improve upon the results ob- essentially removing the balancing act between the interdependent initial

tained with the signi¬cantly more accurate HJ WENO scheme. The concern data across the interface. Their idea replaces equation (7.4) with a com-

is that the HJ WENO scheme might be so accurate that the approximations bination of equations (7.2) and (7.3) along with interpolation to ¬nd the

made by [158] could lower the accuracy of the method. interface location. In [143] marked improvement was shown in using low-

This local constraint is implemented in [158] by the addition of a order HJ ENO schemes, but the authors did not address whether they can

correction term to the right-hand side of equation (7.4), obtain improved results over the recommended HJ WENO discretization

of equation (7.4). Moreover, implementing a high-order accurate version of

φt + S(φo )(| φ| ’ 1) = »δ(φ)| φ|, (7.8)

the scheme in [143] requires a number of ghost cells, as discussed above.

ˆ

where the multidimensional delta function δ = δ(φ)| φ| from equa-

tion (1.19) is used, since the modi¬cations are needed only near the interface

where Ai,j is not trivially equal to either zero or the volume of „¦i,j . The

7.5 The Fast Marching Method

constraint that Ai,j in each cell not change, i.e., (Ai,j )t = 0, is equivalent

to

In the crossing-time approach to constructing signed distance functions, the

H (φ)φt dx = 0, (7.9)

zero isocontour moves in the normal direction, crossing over grid points at

„¦i,j

times equal to their distance from the interface. In this fashion, each grid

Copyright Design by: Sunlight webdesign