Implicit integration extension

Project.toml

@@ -12,18 +12,21 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [weakdeps]
+ImplicitIntegration = "bc256489-3a69-4a66-afc4-127cc87e6182"
 MMG_jll = "86086c02-e288-5929-a127-40944b0018b7"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 MarchingCubes = "299715c1-40a9-479a-aaf9-4a633d36f717"
 
 [extensions]
+ImplicitIntegrationExt = "ImplicitIntegration"
 MMGSurfaceExt = ["MMG_jll", "MarchingCubes"]
 MMGVolumeExt = "MMG_jll"
 MakieExt = "Makie"
 
 [compat]
 DiffResults = "1.1.0"
 ForwardDiff = "1.3.2"
+ImplicitIntegration = "0.2"
 LinearAlgebra = "1"
 MMG_jll = "5"
 Makie = "0.21, 0.22, 0.23, 0.24"

docs/src/boundary-conditions.md

@@ -30,7 +30,7 @@ Here is how it looks in practice:
 ```@example boundary-conditions
 using LevelSetMethods, GLMakie
 grid = CartesianGrid((-1,-1), (1,1), (100, 100))
-ϕ₀    = LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid)
+ϕ₀    = MeshField(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid)
 bc   = PeriodicBC()
 eq   = LevelSetEquation(; ic = deepcopy(ϕ₀), bc, terms = AdvectionTerm((x,t) -> (1,0)))
 fig = Figure(; size = (1200, 300))

docs/src/extension-makie.md

@@ -35,7 +35,7 @@ using LevelSetMethods, GLMakie, LinearAlgebra
 grid = CartesianGrid((-1.5, -1.5, -1.5), (1.5, 1.5, 1.5), (50, 50, 50))
 P1, P2 = (-1, 0, 0), (1, 0, 0)
 b = 1.05
-ϕ = LevelSet(grid) do x
+ϕ = MeshField(grid) do x
   norm(x .- P1)*norm(x .- P2) - b^2
 end
 theme = LevelSetMethods.makie_theme()
@@ -54,5 +54,5 @@ end
 ```
 
 !!! tip "Plotting a `LevelSetEquation`"
-    Calling `plot` on a [`LevelSetEquation`](@ref) defaults to plotting the `LevelSet` given by its
+    Calling `plot` on a [`LevelSetEquation`](@ref) defaults to plotting the `MeshField` given by its
     `current_state`; exactly the same as calling `plot(current_state(equation))`.

docs/src/index.md

@@ -39,7 +39,7 @@ Here is how it looks in practice to create a simple `LevelSetEquation`:
 ```@example ls-intro
 using LevelSetMethods
 grid = CartesianGrid((-1, -1), (1, 1), (100, 100))
-# ϕ    = LevelSet(x -> sqrt(2*x[1]^2 + x[2]^2) - 1/2, grid) # a disk
+# ϕ    = MeshField(x -> sqrt(2*x[1]^2 + x[2]^2) - 1/2, grid) # a disk
 ϕ    = LevelSetMethods.dumbbell(grid) # a predefined shape
 𝐮    = (x,t) -> (-x[2], x[1])
 eq   = LevelSetEquation(;

docs/src/interpolation.md

@@ -2,31 +2,25 @@
 
 LevelSetMethods.jl provides a built-in high-performance piecewise polynomial
 interpolation scheme. This allows you to construct a continuous interpolant from the
-discrete data in a [`LevelSet`](@ref) or [`LevelSetEquation`](@ref). This is useful for
+discrete data in a [`MeshField`](@ref) or [`LevelSetEquation`](@ref). This is useful for
 evaluating the level-set function at arbitrary coordinates, computing analytical
 gradients and Hessians, or for visualization.
 
 ## Basic Usage
 
-To construct an interpolant, use the [`interpolate`](@ref) function:
+Interpolation is enabled by passing `interp_order` when constructing a `MeshField`.
+The field itself becomes callable and evaluates the piecewise polynomial interpolant:
 
 ```@example interpolation
 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(_ -> (LinearExtrapolationBC(), LinearExtrapolationBC()), 2)
-ϕ = LevelSetMethods.add_boundary_conditions(ϕ, bc)
-
-itp = interpolate(ϕ) # cubic interpolation by default (order=3)
+ϕ = LevelSetMethods.star(CartesianGrid(a, b, (50, 50)); interp_order = 3)
 ```
 
-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:
+Once constructed, `ϕ` can be evaluated at any point inside the grid:
 
 ```@example interpolation
-itp(0.5, 0.5)
+ϕ(0.5, 0.5)
 ```
 
 ## Plotting
@@ -39,7 +33,7 @@ using GLMakie
 LevelSetMethods.set_makie_theme!()
 xx = yy = -2:0.01:2
 # Evaluate on a fine grid for plotting
-contour(xx, yy, [itp(x, y) for x in xx, y in yy]; levels = [0], linewidth = 2)
+contour(xx, yy, [ϕ(x, y) for x in xx, y in yy]; levels = [0], linewidth = 2)
 ```
 
 ## Derivatives
@@ -50,10 +44,13 @@ These are zero-allocation and computed using Horner's method on the underlying
 Lagrange polynomials.
 
 ```@example interpolation
-x = (0.1, 0.2)
-val  = itp(x)
-grad = LevelSetMethods.gradient(itp, x)
-hess = LevelSetMethods.hessian(itp, x)
+using StaticArrays
+x = SVector(0.1, 0.2)
+I = LevelSetMethods.compute_index(ϕ, x)
+p = LevelSetMethods.make_interpolant(ϕ, I)
+val  = p(x)
+grad = LevelSetMethods.gradient(p, x)
+hess = LevelSetMethods.hessian(p, x)
 println("Value:    ", val)
 println("Gradient: ", grad)
 ```
@@ -64,15 +61,11 @@ Using it on three-dimensional level sets is identical:
 
 ```@example interpolation
 using LinearAlgebra, StaticArrays
-grid = CartesianGrid((-1.5, -1.5, -1.5), (1.5, 1.5, 1.5), (32, 32, 32))
+grid3 = CartesianGrid((-1.5, -1.5, -1.5), (1.5, 1.5, 1.5), (32, 32, 32))
 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(_ -> (LinearExtrapolationBC(), LinearExtrapolationBC()), 3)
-ϕ3 = LevelSetMethods.add_boundary_conditions(ϕ3, bc3)
-
-itp3 = interpolate(ϕ3)
-println("ϕ(0.5, 0.5, 0.5)   = ", f(SVector(0.5, 0.5, 0.5)))
-println("itp(0.5, 0.5, 0.5) = ", itp3(0.5, 0.5, 0.5))
+f3 = (x) -> norm(x .- P1)*norm(x .- P2) - b^2
+ϕ3 = MeshField(f3, grid3; interp_order = 3)
+println("ϕ(0.5, 0.5, 0.5)   = ", f3(SVector(0.5, 0.5, 0.5)))
+println("ϕ(0.5, 0.5, 0.5) = ", ϕ3(0.5, 0.5, 0.5))
 ```

docs/src/reinitialization.md

@@ -17,15 +17,15 @@ Two approaches are available:
 
 ## Usage
 
-Call `reinitialize!` on a `LevelSet` to reinitialize it in place:
+Call `reinitialize!` on a `MeshField` to reinitialize it in place:
 
 ```@example reinit
 using LevelSetMethods
 using GLMakie
 
 grid = CartesianGrid((-1, -1), (1, 1), (100, 100))
-sdf = LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid)
-ϕ = LevelSet(x -> x[1]^2 + x[2]^2 - 0.5^2, grid)
+sdf = MeshField(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid)
+ϕ = MeshField(x -> x[1]^2 + x[2]^2 - 0.5^2, grid)
 LevelSetMethods.set_makie_theme!()
 fig = Figure(; size = (800, 300))
 ax1 = Axis(fig[1, 1]; title = "Signed Distance")
@@ -107,7 +107,7 @@ When the underlying level set changes, use `LevelSetMethods.update!` to rebuild
 KD-tree, reusing the same interpolation order and tolerances:
 
 ```@example reinit
-ϕ2 = LevelSet(x -> x[1]^2 + x[2]^2 - 0.3^2, grid)  # smaller circle
+ϕ2 = MeshField(x -> x[1]^2 + x[2]^2 - 0.3^2, grid)  # smaller circle
 LevelSetMethods.update!(sdf_obj, ϕ2)
 sdf_obj(SVector(0.3, 0.0))   # should be ≈ 0
 ```

docs/src/terms.md

@@ -43,7 +43,7 @@ point. This is useful if your velocity field is time-independent, or if you only
 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)
+ϕ₀ = MeshField(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid)
 eq = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = ϕ₀, bc = NeumannBC())
 ```
 
@@ -70,7 +70,7 @@ have instead a time-dependent velocity field, we could pass a function to the
 `AdvectionTerm`, and the velocity field would be computed at each time step. For example:
 
 ```@example advection-term
-ϕ₀ = LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid)
+ϕ₀ = MeshField(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid)
 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
@@ -164,11 +164,11 @@ v = zeros(Float64, size(grid)...)
 frozen = abs.(values(ϕext)) .<= 1.5Δ
 for I in CartesianIndices(v)
     frozen[I] || continue
-    x = grid[I]
+    x = getnode(grid, I)
     v[I] = 0.2 + 0.1 * cos(2π * atan(x[2], x[1]))
 end
 extend_along_normals!(v, ϕext; frozen, nb_iters = 80)
-term = NormalMotionTerm(MeshField(v, grid, nothing))
+term = NormalMotionTerm(MeshField(v, grid))
 term
 ```
 
@@ -197,7 +197,7 @@ r0 = 0.5
 α = π / 100.0
 R = [cos(α) -sin(α); sin(α) cos(α)]
 M = R * [1/0.06^2 0; 0 1/(4π^2)] * R'
-ϕ = LevelSet(grid) do (x, y)
+ϕ = MeshField(grid) do (x, y)
     r = sqrt(x^2 + y^2)
     θ = atan(y, x)
     result = 1e30
@@ -239,8 +239,8 @@ and its signed distance function:
 ```@example reinitialization-term
 using LevelSetMethods, GLMakie
 grid = CartesianGrid((-1,-1), (1,1), (100, 100))
-ϕ = LevelSet(x -> x[1]^2 + x[2]^2 - 0.5^2, grid) # circle level-set, but not a signed distance function
-sdf = LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) # signed distance function
+ϕ = MeshField(x -> x[1]^2 + x[2]^2 - 0.5^2, grid) # circle level-set, but not a signed distance function
+sdf = MeshField(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) # signed distance function
 LevelSetMethods.set_makie_theme!()
 fig = Figure(; size = (800, 400))
 ax = Axis(fig[1,1], title = "Signed distance function")
@@ -271,7 +271,7 @@ Alternatively, you can use a modified reinitialization term that applies the sig
   \phi_t + \text{sign}(\phi_0) \left( |\nabla \phi| - 1 \right) = 0
 ```
 
-To enable this behavior, simply pass a `LevelSet` object to the `EikonalReinitializationTerm`:
+To enable this behavior, simply pass a `MeshField` object to the `EikonalReinitializationTerm`:
 
 ```@example reinitialization-term
 eq = LevelSetEquation(; terms = (EikonalReinitializationTerm(ϕ),), ic = deepcopy(ϕ), bc = NeumannBC())

ext/ImplicitIntegrationExt.jl

@@ -0,0 +1,34 @@
+module ImplicitIntegrationExt
+
+using LevelSetMethods
+import ImplicitIntegration as II
+using StaticArrays
+
+_single_cell(I::CartesianIndex{N}) where {N} = CartesianIndices(ntuple(d -> I[d]:I[d], Val(N)))
+
+function LevelSetMethods.quadrature(ϕ::LevelSetMethods.AbstractMeshField; order, surface = false)
+    LevelSetMethods.has_interpolation(ϕ) ||
+        error(
+        "quadrature requires a MeshField constructed with `interp_order`. " *
+            "Pass `interp_order=k` to the MeshField or NarrowBandMeshField constructor."
+    )
+
+    grid = LevelSetMethods.mesh(ϕ)
+    N = ndims(grid)
+
+    results = Tuple{CartesianIndices{N, NTuple{N, UnitRange{Int}}}, II.Quadrature}[]
+
+    for I in LevelSetMethods.cellindices(ϕ)
+        LevelSetMethods.proven_empty(ϕ, I; surface) && continue
+        bp_lsm = LevelSetMethods.make_interpolant(ϕ, I)
+        bp = II.BernsteinPolynomial(copy(bp_lsm.coeffs), bp_lsm.low_corner, bp_lsm.high_corner)
+        out = II.quadgen(bp, bp.low_corner, bp.high_corner; order, surface)
+        if !isempty(out.quad.coords)
+            push!(results, (_single_cell(I), out.quad))
+        end
+    end
+
+    return results
+end
+
+end # module

ext/MMGSurfaceExt.jl

@@ -11,18 +11,18 @@ end
 """
     export_surface_mesh(eq::LevelSetEquation, args...; kwargs...)
 
-Call [`export_surface_mesh`](@ref LSM.export_surface_mesh(::LSM.LevelSet)) on
+Call [`export_surface_mesh`](@ref LSM.export_surface_mesh(::LSM.MeshField)) on
 `current_state(eq)`.
 """
 function LSM.export_surface_mesh(eq::LSM.LevelSetEquation, args...; kwargs...)
     return LSM.export_surface_mesh(LSM.current_state(eq), output; hgrad, hmin, hmax, hausd)
 end
 
 """
-    export_surface_mesh(ϕ::LevelSet, output::String;
+    export_surface_mesh(ϕ::MeshField, output::String;
         hgrad = nothing, hmin = nothing, hmax = nothing, hausd = nothing)
 
-Compute a mesh of the [`LevelSet`](@ref LSM.LevelSet) `ϕ` zero contour using MMGs_O3.
+Compute a mesh of the [`MeshField`](@ref LSM.MeshField) `ϕ` zero contour using MMGs_O3.
 
 `hgrad` control the growth ratio between two adjacent edges
 
@@ -37,7 +37,7 @@ boundary and the reconstructed ideal boundary
     Only works for 3 dimensional level-set.
 """
 function LSM.export_surface_mesh(
-        ϕ::LSM.LevelSet,
+        ϕ::LSM.MeshField,
         output::String;
         hgrad = nothing,
         hmin = nothing,

ext/MMGVolumeExt.jl

@@ -11,17 +11,17 @@ end
 """
     export_volume_mesh(eq::LevelSetEquation, output; kwargs...)
 
-Call [`export_volume_mesh`](@ref LSM.export_volume_mesh(::LSM.LevelSet)) on `current_state(eq)`.
+Call [`export_volume_mesh`](@ref LSM.export_volume_mesh(::LSM.MeshField)) on `current_state(eq)`.
 """
 function LSM.export_volume_mesh(eq::LSM.LevelSetEquation, args...; kwargs...)
     return LSM.export_volume_mesh(LSM.current_state(eq), args...; kwargs...)
 end
 
 """
-    export_volume_mesh(ϕ::LevelSet, output::String;
+    export_volume_mesh(ϕ::MeshField, output::String;
         hgrad = nothing, hmin = nothing, hmax = nothing, hausd = nothing)
 
-Compute a mesh of the domains associated with [`LevelSet`](@ref LSM.LevelSet) `eq` using
+Compute a mesh of the domains associated with [`MeshField`](@ref LSM.MeshField) `eq` using
 either MMG2d_O3 or MMG3d_O3.
 
 `hgrad` control the growth ratio between two adjacent edges.
@@ -39,7 +39,7 @@ For more information, see the official [MMG documentation](http://www.mmgtools.o
     Only works for 2 and 3 dimensional level-set.
 """
 function LSM.export_volume_mesh(
-        ϕ::LSM.LevelSet,
+        ϕ::LSM.MeshField,
         output::String;
         hgrad = nothing,
         hmin = nothing,