Narrow band implementation and major refactoring

.gitignore

@@ -18,3 +18,4 @@ sandbox
 *.msh
 *.sol
 *.mesh
+*.markdown-preview.html

docs/src/boundary-conditions.md

@@ -7,10 +7,10 @@ The following boundary conditions are available:
 | [`PeriodicBC`](@ref) | Periodic (wrap-around) |
 | [`DirichletBC`](@ref) | Prescribed boundary value |
 | [`ExtrapolationBC{P}`](@ref ExtrapolationBC) | P-th order one-sided polynomial extrapolation |
-| `NeumannBC` | Alias for `ExtrapolationBC{1}` (constant extension, ∂ϕ/∂n = 0) |
-| `NeumannGradientBC` | Alias for `ExtrapolationBC{2}` (linear extrapolation, ∂²ϕ/∂n² = 0) |
+| `NeumannBC` | Alias for `ExtrapolationBC{0}` (constant extension, ∂ϕ/∂n = 0) |
+| `LinearExtrapolationBC` | Alias for `ExtrapolationBC{1}` (linear extrapolation, ∂²ϕ/∂n² = 0) |
 
-`ExtrapolationBC{P}` uses the `P` nearest interior cells to build a degree `P-1` polynomial
+`ExtrapolationBC{P}` uses the `P+1` nearest interior cells to build a degree-`P` polynomial
 and extrapolates it into the ghost region. Higher `P` gives smoother outflow at the cost of
 a wider stencil.
 
@@ -32,7 +32,7 @@ using LevelSetMethods, GLMakie
 grid = CartesianGrid((-1,-1), (1,1), (100, 100))
 ϕ₀    = LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid)
 bc   = PeriodicBC()
-eq   = LevelSetEquation(; levelset = deepcopy(ϕ₀), bc, terms = AdvectionTerm((x,t) -> (1,0)))
+eq   = LevelSetEquation(; ic = deepcopy(ϕ₀), bc, terms = AdvectionTerm((x,t) -> (1,0)))
 fig = Figure(; size = (1200, 300))
 for (n,t) in enumerate([0.0, 0.5, 0.75, 1.0])
     integrate!(eq, t)
@@ -42,10 +42,10 @@ end
 fig
 ```
 
-Changing `PeriodicBC()` to `NeumannBC()` gives allows for the level-set to "leak" out of the domain:
+Changing `PeriodicBC()` to `NeumannBC()` allows the level-set to "leak" out of the domain:
 
 ```@example boundary-conditions
-eq   = LevelSetEquation(; levelset = deepcopy(ϕ₀), bc = NeumannBC(), terms = AdvectionTerm((x,t) -> (1,0)))
+eq   = LevelSetEquation(; ic = deepcopy(ϕ₀), bc = NeumannBC(), terms = AdvectionTerm((x,t) -> (1,0)))
 fig = Figure(; size = (1200, 300))
 for (n,t) in enumerate([0.0, 0.5, 0.75, 1.0])
     integrate!(eq, t)
@@ -59,7 +59,7 @@ To combine both boundary conditions you can use
 
 ```@example boundary-conditions
 bc = (NeumannBC(), PeriodicBC()) # Neumann in x, periodic in y
-eq   = LevelSetEquation(; levelset = deepcopy(ϕ₀), bc, terms = AdvectionTerm((x,t) -> (1,1)))
+eq   = LevelSetEquation(; ic = deepcopy(ϕ₀), bc, terms = AdvectionTerm((x,t) -> (1,1)))
 fig = Figure(; size = (1200, 300))
 for (n,t) in enumerate([0.0, 0.5, 0.75, 1.0])
     integrate!(eq, t)
@@ -70,11 +70,11 @@ fig
 ```
 
 For higher-order outflow you can use `ExtrapolationBC{P}` directly. For example,
-`ExtrapolationBC{5}` fits a degree-4 polynomial through the 5 nearest interior cells:
+`ExtrapolationBC{5}` fits a degree-5 polynomial through the 6 nearest interior cells:
 
 ```@example boundary-conditions
 bc = ExtrapolationBC(5)
-eq   = LevelSetEquation(; levelset = deepcopy(ϕ₀), bc, terms = AdvectionTerm((x,t) -> (1,0)))
+eq   = LevelSetEquation(; ic = deepcopy(ϕ₀), bc, terms = AdvectionTerm((x,t) -> (1,0)))
 fig = Figure(; size = (1200, 300))
 for (n,t) in enumerate([0.0, 0.5, 0.75, 1.0])
     integrate!(eq, t)

docs/src/example-mass-conservation.md

@@ -39,7 +39,7 @@ timestamps = range(0, tf; length = nframes + 1)
 function track_volume(integrator)
     ϕ = deepcopy(ϕ₀)
     eq = LevelSetEquation(;
-        levelset = ϕ,
+        ic = ϕ,
         terms = AdvectionTerm((x, t) -> (-x[2], x[1])),
         bc = NeumannBC(),
         integrator,

docs/src/example-shape-optim.md

@@ -93,9 +93,9 @@ term1 = NormalMotionTerm(MeshField(X -> 0.0, grid))
 term2 = CurvatureTerm(MeshField(X -> -1.0, grid))
 terms = (term1, term2)
 
-bc = NeumannGradientBC()
+bc = LinearExtrapolationBC()
 integrator = ForwardEuler(0.5)
-eq = LevelSetEquation(; terms, integrator, levelset = ϕ, t = 0, bc)
+eq = LevelSetEquation(; terms, integrator, ic = ϕ, t = 0, bc)
 
 using GLMakie
 

docs/src/example-zalesak.md

@@ -40,7 +40,7 @@ simulates rotation:
 ```@example zalesak_disk_example
 # Define the advection equation with a rotational velocity field
 eq = LevelSetEquation(;
-    levelset = ϕ,
+    ic = ϕ,
     terms = AdvectionTerm((x, t) -> (-x[2], x[1])),
     bc = NeumannBC(),
 )
@@ -97,7 +97,7 @@ rec = LevelSetMethods.rectangle(
 )
 ϕ = setdiff(disk, rec)
 eq = LevelSetEquation(;
-    levelset = ϕ,
+    ic = ϕ,
     terms = AdvectionTerm((x, t) -> π * SVector(x[2], -x[1], 0)),
     bc = NeumannBC(),
 )

docs/src/extension-mmg.md

@@ -1,6 +1,6 @@
 # [MMG extension](@id extension-mmg)
 
-This extension provides functions to generate meshes of level-set functions using [MMG](https://www.mmgtools.org/).
+This extension provides functions to generate meshes of level-set functions using [MMG](https://github.com/MmgTools/mmg).
 It defines two methods: `export_volume_mesh` and `export_surface_mesh`.
 For both of them, it is possible to control the size of the generated mesh using the following optional parameters:
 

docs/src/index.md

@@ -44,8 +44,8 @@ grid = CartesianGrid((-1, -1), (1, 1), (100, 100))
 𝐮    = (x,t) -> (-x[2], x[1])
 eq   = LevelSetEquation(;
   terms = (AdvectionTerm(𝐮),),
-  levelset = ϕ,
-  bc = PeriodicBC()
+  ic = ϕ,
+  bc = NeumannBC()
 )
 ```
 
@@ -99,10 +99,10 @@ Here is what the `.gif` file looks like:
 For more interesting applications and advanced usage, see the examples section!
 
 !!! note "Other resources"
-    There is an almost one-to-one correspondance between each of the [`LevelSetTerm`](@ref)s
-    described above and individual chapters of the book by Osher and Fedwick on level set
+    There is an almost one-to-one correspondence between each of the [`LevelSetTerm`](@ref)s
+    described above and individual chapters of the book by Osher and Fedkiw on level set
     methods [osher2003level](@cite), so users interested in digging deeper into the
-    theory/algorithms are encourage to consult that refenrence. We also drew some
+    theory/algorithms are encouraged to consult that reference. We also drew some
     inspiration from the great Matlab library `ToolboxLS` by Ian Mitchell
     [mitchell2007toolbox](@cite).
 

docs/src/interpolation.md

@@ -15,13 +15,13 @@ using LevelSetMethods
 a, b = (-2.0, -2.0), (2.0, 2.0)
 ϕ   = LevelSetMethods.star(CartesianGrid(a, b, (50, 50)))
 # Add boundary conditions for safe evaluation near edges
-bc = ntuple(_ -> (NeumannGradientBC(), NeumannGradientBC()), 2)
+bc = ntuple(_ -> (LinearExtrapolationBC(), LinearExtrapolationBC()), 2)
 ϕ = LevelSetMethods.add_boundary_conditions(ϕ, bc)
 
 itp = interpolate(ϕ) # cubic interpolation by default (order=3)
 ```
 
-The returned object is a [`PiecewisePolynomialInterpolation`](@ref LevelSetMethods.PiecewisePolynomialInterpolation), which is callable and
+The returned object is a [`PiecewisePolynomialInterpolant`](@ref LevelSetMethods.PiecewisePolynomialInterpolant), which is callable and
 efficient. Once constructed, the interpolant can be used to evaluate the level-set function
 anywhere inside (and even slightly outside, using boundary conditions) the grid:
 
@@ -69,7 +69,7 @@ P1, P2 = (-1.0, 0.0, 0.0), (1.0, 0.0, 0.0)
 b = 1.05
 f = (x) -> norm(x .- P1)*norm(x .- P2) - b^2
 ϕ3 = LevelSet(f, grid)
-bc3 = ntuple(_ -> (NeumannGradientBC(), NeumannGradientBC()), 3)
+bc3 = ntuple(_ -> (LinearExtrapolationBC(), LinearExtrapolationBC()), 3)
 ϕ3 = LevelSetMethods.add_boundary_conditions(ϕ3, bc3)
 
 itp3 = interpolate(ϕ3)

docs/src/reinitialization.md

@@ -56,8 +56,8 @@ Pass `reinit` to [`LevelSetEquation`](@ref) to reinitialize automatically every
 ```julia
 eq = LevelSetEquation(;
     terms  = (AdvectionTerm(𝐮),),
-    levelset = ϕ,
-    bc     = PeriodicBC(),
+    ic = ϕ,
+    bc     = NeumannBC(),
     reinit = 5,          # reinitialize every 5 time steps
 )
 ```
@@ -69,8 +69,8 @@ directly:
 ```julia
 eq = LevelSetEquation(;
     terms    = (AdvectionTerm(𝐮),),
-    levelset = ϕ,
-    bc       = PeriodicBC(),
+    ic = ϕ,
+    bc       = NeumannBC(),
     reinit   = NewtonReinitializer(; reinit_freq = 5, upsample = 4),
 )
 ```

docs/src/terms.md

@@ -44,7 +44,7 @@ grid points. Lets construct a level-set equation with an advection term:
 
 ```@example advection-term
 ϕ₀ = LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid)
-eq = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), levelset = ϕ₀, bc = PeriodicBC())
+eq = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = ϕ₀, bc = NeumannBC())
 ```
 
 To see how the advection term affects the level-set, we can solve the equation for a few
@@ -71,7 +71,7 @@ have instead a time-dependent velocity field, we could pass a function to the
 
 ```@example advection-term
 ϕ₀ = LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid)
-eq = LevelSetEquation(; terms = (AdvectionTerm((x,t) -> SVector(x[1]^2, 0)),), levelset = ϕ₀, bc = PeriodicBC())
+eq = LevelSetEquation(; terms = (AdvectionTerm((x,t) -> SVector(x[1]^2, 0)),), ic = ϕ₀, bc = NeumannBC())
 fig = Figure(; size = (1200, 300))
 # create a 2 x 2 figure
 for (n,t) in enumerate([0.0, 0.5, 0.75, 1.0])
@@ -102,8 +102,8 @@ let us compare both schemes for a purely rotational velocity field:
 𝐮  = MeshField(grid) do (x,y)
     SVector(-y, x)
 end
-eq_upwind = LevelSetEquation(; terms = AdvectionTerm(𝐮, Upwind()), levelset = deepcopy(ϕ₀), bc = PeriodicBC())
-eq_weno   = LevelSetEquation(; terms = AdvectionTerm(𝐮), levelset  = deepcopy(ϕ₀), bc = PeriodicBC())
+eq_upwind = LevelSetEquation(; terms = AdvectionTerm(𝐮, Upwind()), ic = deepcopy(ϕ₀), bc = NeumannBC())
+eq_weno   = LevelSetEquation(; terms = AdvectionTerm(𝐮), ic = deepcopy(ϕ₀), bc = NeumannBC())
 fig = Figure(size = (1000, 400))
 ax = Axis(fig[1,1], title = "Initial")
 plot!(ax, eq_upwind)
@@ -134,7 +134,7 @@ using LevelSetMethods
 using GLMakie
 grid = CartesianGrid((-2,-2), (2,2), (100, 100))
 ϕ = LevelSetMethods.star(grid)
-eq = LevelSetEquation(; terms = (NormalMotionTerm((x,t) -> 0.5),), levelset = ϕ, bc = PeriodicBC())
+eq = LevelSetEquation(; terms = (NormalMotionTerm((x,t) -> 0.5),), ic = ϕ, bc = NeumannBC())
 fig = Figure(; size = (1200, 300))
 for (n,t) in enumerate([0.0, 0.5, 0.75, 1.0])
     integrate!(eq, t)
@@ -185,7 +185,7 @@ where ``\kappa = \nabla \cdot (\nabla \phi / |\nabla \phi|)`` is the mean curvat
 coefficient ``b`` should be negative; a positive value of ``b`` would yield an ill-posed
 evolution problem (akin to a negative diffusion coefficient).
 
-Here is the classic example of motion by mean curavature for a spiral-like level-set:
+Here is the classic example of motion by mean curvature for a spiral-like level-set:
 
 ```@example curvature-term
 using LevelSetMethods, GLMakie
@@ -208,7 +208,7 @@ M = R * [1/0.06^2 0; 0 1/(4π^2)] * R'
     end
     return result
 end
-eq = LevelSetEquation(; terms = (CurvatureTerm((x,t) -> -0.1),), levelset = ϕ, bc = PeriodicBC())
+eq = LevelSetEquation(; terms = (CurvatureTerm((x,t) -> -0.1),), ic = ϕ, bc = NeumannBC())
 fig = Figure(; size = (1200, 300))
 for (n,t) in enumerate([0.0, 0.1, 0.2, 0.3])
     integrate!(eq, t)
@@ -253,7 +253,7 @@ fig
 We will now evolve the level-set using the reinitialization term:
 
 ```@example reinitialization-term
-eq = LevelSetEquation(; terms = (EikonalReinitializationTerm(),), levelset = deepcopy(ϕ), bc = PeriodicBC())
+eq = LevelSetEquation(; terms = (EikonalReinitializationTerm(),), ic = deepcopy(ϕ), bc = NeumannBC())
 fig = Figure(; size = (1200, 300))
 for (n,t) in enumerate([0.0, 0.25, 0.5, 0.75])
     integrate!(eq, t)
@@ -274,7 +274,7 @@ Alternatively, you can use a modified reinitialization term that applies the sig
 To enable this behavior, simply pass a `LevelSet` object to the `EikonalReinitializationTerm`:
 
 ```@example reinitialization-term
-eq = LevelSetEquation(; terms = (EikonalReinitializationTerm(ϕ),), levelset = deepcopy(ϕ), bc = PeriodicBC())
+eq = LevelSetEquation(; terms = (EikonalReinitializationTerm(ϕ),), ic = deepcopy(ϕ), bc = NeumannBC())
 fig = Figure(; size = (1200, 300))
 for (n,t) in enumerate([0.0, 0.25, 0.5, 0.75])
     integrate!(eq, t)

ext/MMGSurfaceExt.jl

@@ -44,7 +44,7 @@ function LSM.export_surface_mesh(
         hmax = nothing,
         hausd = nothing,
     )
-    N = LSM.dimension(ϕ)
+    N = ndims(ϕ)
     if N != 3
         throw(ArgumentError("export_mesh of $N dimensional level-set not supported."))
     end

ext/MMGVolumeExt.jl

@@ -46,7 +46,7 @@ function LSM.export_volume_mesh(
         hmax = nothing,
         hausd = nothing,
     )
-    N = LSM.dimension(ϕ)
+    N = ndims(ϕ)
     if N != 2 && N != 3
         throw(ArgumentError("export_mesh of $N dimensional level-set not supported."))
     end