Allow for a "freezed" signed function in ReinitializationTerm

Project.toml

@@ -1,7 +1,7 @@
 name = "LevelSetMethods"
 uuid = "ca3df947-0575-4ae6-b9a3-78df393a7451"
 authors = ["Luiz M. Faria and Nicolas Lebbe"]
-version = "0.1.3"
+version = "0.1.4"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/src/terms.md

@@ -225,7 +225,7 @@ fig
 We will now evolve the level-set using the reinitialization term:
 
 ```@example reinitialization-term
-eq = LevelSetEquation(; terms = (ReinitializationTerm(),), levelset = ϕ, bc = PeriodicBC())
+eq = LevelSetEquation(; terms = (ReinitializationTerm(),), levelset = deepcopy(ϕ), bc = PeriodicBC())
 fig = Figure(; size = (1200, 300))
 for (n,t) in enumerate([0.0, 0.25, 0.5, 0.75])
     integrate!(eq, t)
@@ -235,4 +235,25 @@ end
 fig
 ```
 
-Note that `ϕ` converges to the signed distance function `sdf` shown in the first figure.
+Observe that as the reinitialization equation evolves, `ϕ` approaches the signed distance function `sdf` depicted in the first figure.
+
+Alternatively, you can use a modified reinitialization term that applies the sign function to the *initial level-set function* only:
+
+```math
+  \phi_t + \text{sign}(\phi_0) \left( |\nabla \phi| - 1 \right) = 0
+```
+
+To enable this behavior, simply pass a `LevelSet` object to the `ReinitializationTerm`:
+
+```@example reinitialization-term
+eq = LevelSetEquation(; terms = (ReinitializationTerm(ϕ),), levelset = deepcopy(ϕ), bc = PeriodicBC())
+fig = Figure(; size = (1200, 300))
+for (n,t) in enumerate([0.0, 0.25, 0.5, 0.75])
+    integrate!(eq, t)
+    ax = Axis(fig[1,n], title = "t = $t")
+    contour!(ax, LevelSetMethods.current_state(eq); levels = [0.25, 0, 0.5], labels = true, labelsize = 14)
+end
+fig
+```
+
+The outcome closely matches that of the previous approach.

src/levelsetequation.jl

@@ -118,6 +118,7 @@ buffers(ls::LevelSetEquation) = ls.buffers
 time_integrator(ls) = ls.integrator
 terms(ls) = ls.terms
 boundary_conditions(ls) = boundary_conditions(ls.state)
+mesh(ls::LevelSetEquation) = mesh(ls.state)
 
 """
     integrate!(ls::LevelSetEquation,tf,Δt=Inf)

src/levelsetterms.jl

@@ -79,6 +79,8 @@ function _compute_cfl(term::AdvectionTerm{V}, ϕ, I, t) where {V}
     elseif V <: Function
         x = mesh(ϕ)[I]
         velocity(term)(x, t)
+    else
+        error("velocity field type $V not supported")
     end
     Δx = meshsize(ϕ)
     return 1 / maximum(abs.(𝐮) ./ Δx)
@@ -208,19 +210,51 @@ end
     struct ReinitializationTerm <: LevelSetTerm
 
 Level-set term representing  `sign(ϕ) (|∇ϕ| - 1)`. This `LevelSetTerm` should be
-used for reinitializing the level set into a signed distance function: for a
-sufficiently large number of time steps this term allows one to solve the
-Eikonal equation |∇ϕ| = 1.
+used for reinitializing the level set into a signed distance function: for a sufficiently
+large number of time steps this term allows one to solve the Eikonal equation |∇ϕ| = 1.
+
+There are two ways of constructing a `ReinitializationTerm`:
+
+  - using `ReinitializationTerm(ϕ₀::LevelSet)` precomputes the `sign` term on the initial
+    level set `ϕ₀`, as in equation 7.5 of Osher and Fedkiw;
+  - using `ReinitializationTerm()` constructs a term that computes the `sign` term
+    on-the-fly at each time step, as in equation 7.6 of Osher and Fedkiw.
 """
-@kwdef struct ReinitializationTerm <: LevelSetTerm end
+struct ReinitializationTerm{T} <: LevelSetTerm
+    S₀::T
+end
+
+function ReinitializationTerm(ϕ₀::LevelSet)
+    Δx = minimum(meshsize(ϕ₀))
+    # equation 7.5 of Osher and Fedkiw
+    S₀ = map(CartesianIndices(mesh(ϕ₀))) do I
+        return ϕ₀[I] / sqrt(ϕ₀[I]^2 + Δx^2)
+    end
+    return ReinitializationTerm(S₀)
+end
+ReinitializationTerm() = ReinitializationTerm(nothing)
 
-Base.show(io::IO, t::ReinitializationTerm) = print(io, "sign(ϕ) (|∇ϕ| - 1)")
+function Base.show(io::IO, t::ReinitializationTerm)
+    S₀ = t.S₀
+    if isnothing(S₀)
+        print(io, "sign(ϕ) (|∇ϕ| - 1)")
+    else
+        print(io, "sign(ϕ₀) (|∇ϕ| - 1)")
+    end
+    return io
+end
 
-function _compute_term(::ReinitializationTerm, ϕ, I, t)
+function _compute_term(term::ReinitializationTerm, ϕ, I, t)
+    S₀ = term.S₀
     norm_∇ϕ = _compute_∇_norm(sign(ϕ[I]), ϕ, I)
-    # equation 7.6 of Osher and Fedkiw
-    Δx = minimum(meshsize(ϕ))
-    S = ϕ[I] / sqrt(ϕ[I]^2 + norm_∇ϕ^2 * Δx^2)
+    if isnothing(S₀)
+        # equation 7.6 of Osher and Fedkiw
+        Δx = minimum(meshsize(ϕ))
+        S = ϕ[I] / sqrt(ϕ[I]^2 + norm_∇ϕ^2 * Δx^2)
+    else
+        # precomputed S₀ term
+        S = S₀[I]
+    end
     return S * (norm_∇ϕ - 1.0)
 end