<< . .

( 23)

. . >>

φ is constant for all time. Figure 6.1 shows the evolution of a star-shaped
φ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
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
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.
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
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
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
Most ¬‚ames burn with a speed in the normal direction plus extra heating
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
’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.
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
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.,
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
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
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
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
times equal to their distance from the interface. In this fashion, each grid

<< . .

( 23)

. . >>

Copyright Design by: Sunlight webdesign