diff --git a/Project.toml b/Project.toml index 5868f29..36f40d8 100644 --- a/Project.toml +++ b/Project.toml @@ -12,11 +12,13 @@ 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" @@ -24,6 +26,7 @@ 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" diff --git a/docs/src/boundary-conditions.md b/docs/src/boundary-conditions.md index 40d6fcd..6a6e695 100644 --- a/docs/src/boundary-conditions.md +++ b/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)) diff --git a/docs/src/extension-makie.md b/docs/src/extension-makie.md index e1f23af..103bfda 100644 --- a/docs/src/extension-makie.md +++ b/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))`. diff --git a/docs/src/index.md b/docs/src/index.md index 4e05826..dc1f72d 100644 --- a/docs/src/index.md +++ b/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(; diff --git a/docs/src/interpolation.md b/docs/src/interpolation.md index 9e16154..0bfff87 100644 --- a/docs/src/interpolation.md +++ b/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)) ``` diff --git a/docs/src/reinitialization.md b/docs/src/reinitialization.md index b0ff7e1..aca9654 100644 --- a/docs/src/reinitialization.md +++ b/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 ``` diff --git a/docs/src/terms.md b/docs/src/terms.md index 29ae025..2c7ecfa 100644 --- a/docs/src/terms.md +++ b/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()) diff --git a/ext/ImplicitIntegrationExt.jl b/ext/ImplicitIntegrationExt.jl new file mode 100644 index 0000000..f4f0cf7 --- /dev/null +++ b/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 diff --git a/ext/MMGSurfaceExt.jl b/ext/MMGSurfaceExt.jl index d694547..1cd4f0e 100644 --- a/ext/MMGSurfaceExt.jl +++ b/ext/MMGSurfaceExt.jl @@ -11,7 +11,7 @@ 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...) @@ -19,10 +19,10 @@ function LSM.export_surface_mesh(eq::LSM.LevelSetEquation, args...; kwargs...) 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, diff --git a/ext/MMGVolumeExt.jl b/ext/MMGVolumeExt.jl index 1b6193e..7b46ba1 100644 --- a/ext/MMGVolumeExt.jl +++ b/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, diff --git a/ext/MakieExt.jl b/ext/MakieExt.jl index 24d55af..6f9d286 100644 --- a/ext/MakieExt.jl +++ b/ext/MakieExt.jl @@ -28,7 +28,7 @@ function LSM.set_makie_theme!() return Makie.set_theme!(LSM.makie_theme()) end -# function Makie.plottype(ϕ::LSM.LevelSet) +# function Makie.plottype(ϕ::LSM.MeshField) # N = ndims(ϕ) # if N == 2 # return Contour @@ -41,14 +41,14 @@ end function Makie.convert_arguments( P::Union{Type{<:Contour}, Type{<:Contourf}, Type{<:Heatmap}}, - ϕ::LSM.LevelSet, + ϕ::LSM.MeshField, ) ndims(ϕ) == 2 || throw(ArgumentError("Contour plot only supported for 2D level-sets.")) return Makie.convert_arguments(P, _contour_plot(ϕ)...) end -function Makie.convert_arguments(P::Type{<:Volume}, ϕ::LSM.LevelSet) +function Makie.convert_arguments(P::Type{<:Volume}, ϕ::LSM.MeshField) ndims(ϕ) == 3 || throw(ArgumentError("Volume plot only supported for 3D level-sets.")) x, y, z, v = _volume_plot(ϕ) @@ -57,7 +57,7 @@ end function Makie.convert_arguments( P::Union{Type{<:Contour}, Type{<:Contourf}, Type{<:Heatmap}}, - nb::LSM.NarrowBandLevelSet, + nb::LSM.NarrowBandMeshField, ) ndims(nb) == 2 || throw(ArgumentError("Contour plot only supported for 2D level-sets.")) @@ -67,7 +67,7 @@ function Makie.convert_arguments( return Makie.convert_arguments(P, x, y, v) end -function Makie.convert_arguments(P::Type{<:Volume}, nb::LSM.NarrowBandLevelSet) +function Makie.convert_arguments(P::Type{<:Volume}, nb::LSM.NarrowBandMeshField) ndims(nb) == 3 || throw(ArgumentError("Volume plot only supported for 3D level-sets.")) msh = LSM.mesh(nb) @@ -78,14 +78,14 @@ function Makie.convert_arguments(P::Type{<:Volume}, nb::LSM.NarrowBandLevelSet) return Makie.convert_arguments(P, xlims, ylims, zlims, v) end -function _contour_plot(ϕ::LSM.LevelSet) +function _contour_plot(ϕ::LSM.MeshField) msh = LSM.mesh(ϕ) x, y = LSM.xgrid(msh), LSM.ygrid(msh) v = LSM.values(ϕ) return x, y, v end -function _volume_plot(ϕ::LSM.LevelSet) +function _volume_plot(ϕ::LSM.MeshField) msh = LSM.mesh(ϕ) x, y, z = LSM.xgrid(msh), LSM.ygrid(msh), LSM.zgrid(msh) v = LSM.values(ϕ) @@ -93,8 +93,9 @@ function _volume_plot(ϕ::LSM.LevelSet) return xlims, ylims, zlims, v end -# Build a dense array from a NarrowBandLevelSet, with NaN outside the active band. -function _nb_to_dense(nb::LSM.NarrowBandLevelSet{N, T}) where {N, T} +# Build a dense array from a narrow-band field, with NaN outside the active band. +function _nb_to_dense(nb::LSM.NarrowBandMeshField) + T = eltype(nb) msh = LSM.mesh(nb) arr = fill(T(NaN), size(msh)) for (I, v) in LSM.values(nb) @@ -104,12 +105,12 @@ function _nb_to_dense(nb::LSM.NarrowBandLevelSet{N, T}) where {N, T} end # Collect active cell rectangles for 2D narrow band. -function _active_cell_rects(nb::LSM.NarrowBandLevelSet{2}) +function _active_cell_rects(nb::LSM.NarrowBandMeshField{<:Any, <:LSM.AbstractMesh{2}}) grid = LSM.mesh(nb) h = LSM.meshsize(grid) rects = Rect2f[] - for I in LSM.active_cells(nb) - x, y = grid[I] + for I in LSM.cellindices(nb) + x, y = LSM.getnode(grid, I) push!(rects, Rect2f(x, y, h[1], h[2])) end return rects @@ -132,9 +133,9 @@ function _grid_linesegments(msh::LSM.CartesianGrid{2}) end # Collect active node coordinates as a vector of Points (3D only). -function _active_node_coords(nb::LSM.NarrowBandLevelSet{3}) +function _active_node_coords(nb::LSM.NarrowBandMeshField{<:Any, <:LSM.AbstractMesh{3}}) msh = LSM.mesh(nb) - return [Point3f(msh[I]) for I in LSM.active_indices(nb)] + return [Point3f(LSM.getnode(msh, I)) for I in LSM.nodeindices(nb)] end @@ -146,13 +147,13 @@ function Makie.plot!(p::LevelSetPlot) eq = p.eq ϕ = @lift LSM.current_state($eq) N = @lift ndims($ϕ) - is_nb = to_value(ϕ) isa LSM.NarrowBandLevelSet + is_nb = @lift $ϕ isa LSM.NarrowBandMeshField if to_value(N) == 2 if to_value(p.showgrid) segs = _grid_linesegments(LSM.mesh(to_value(ϕ))) linesegments!(p, segs; color = (:black, 0.15), linewidth = 0.5) end - if is_nb + if to_value(is_nb) rects = @lift _active_cell_rects($ϕ) poly!(p, rects; color = (:steelblue, 0.2), strokewidth = 0.5, strokecolor = (:steelblue, 0.5)) else @@ -160,7 +161,7 @@ function Makie.plot!(p::LevelSetPlot) end contour!(p, ϕ; levels = [0], linewidth = 2, color = :black, overdraw = true) elseif to_value(N) == 3 - if is_nb + if to_value(is_nb) pts = @lift _active_node_coords($ϕ) scatter!(p, pts; color = (:steelblue, 0.3), markersize = 4) end @@ -172,8 +173,7 @@ function Makie.plot!(p::LevelSetPlot) end Makie.plottype(::LSM.LevelSetEquation) = LevelSetPlot -Makie.plottype(::LSM.LevelSet) = LevelSetPlot -Makie.plottype(::LSM.NarrowBandLevelSet) = LevelSetPlot +Makie.plottype(::LSM.AbstractMeshField) = LevelSetPlot ## Pick correct axis type based on dimension of the level-set function Makie.preferred_axis_type(p::LevelSetPlot) diff --git a/src/LevelSetMethods.jl b/src/LevelSetMethods.jl index 49c5dab..5bbadc0 100644 --- a/src/LevelSetMethods.jl +++ b/src/LevelSetMethods.jl @@ -9,6 +9,7 @@ using NearestNeighbors include("meshes.jl") include("boundaryconditions.jl") +include("bernstein.jl") include("meshfield.jl") include("levelset.jl") include("derivatives.jl") @@ -21,17 +22,17 @@ include("timestepping.jl") include("logging.jl") include("levelsetequation.jl") -export AdvectionTerm, +export AbstractMeshField, + AdvectionTerm, CartesianGrid, CurvatureTerm, DirichletBC, EikonalReinitializationTerm, ExtrapolationBC, ForwardEuler, - LevelSet, LevelSetEquation, MeshField, - NarrowBandLevelSet, + NarrowBandMeshField, NeumannBC, LinearExtrapolationBC, NewtonReinitializer, @@ -42,11 +43,16 @@ export AdvectionTerm, SemiImplicitI2OE, Upwind, WENO5, + cellindices, + check_real_valued, current_state, current_time, extend_along_normals!, + getcell, + getnode, integrate!, - interpolate, + nodeindices, + quadrature, reinitialize! @@ -80,4 +86,22 @@ Requires the `MMG` extension to be loaded. """ function export_surface_mesh end +""" + quadrature(ϕ::AbstractMeshField; order, surface=false) + +Generate a quadrature for the implicit domain defined by `ϕ`. +If `surface=true`, generate a quadrature for the interface `ϕ=0`; +otherwise for the interior `ϕ < 0`. + +Returns a `Vector` of `(region, quadrature)` pairs, where `region` is a +single-cell `CartesianIndices` and `quadrature` is an `ImplicitIntegration.Quadrature`. + +!!! note + Requires loading `ImplicitIntegration.jl` to activate the extension. + `ϕ` must be constructed with `interp_order` set. +""" +function quadrature(ϕ; order, surface = false) + error("ImplicitIntegration extension not loaded. Load ImplicitIntegration to use this functionality.") +end + end # module diff --git a/src/bernstein.jl b/src/bernstein.jl new file mode 100644 index 0000000..6ac6b5f --- /dev/null +++ b/src/bernstein.jl @@ -0,0 +1,234 @@ +""" + struct BernsteinPolynomial{N, T} <: Function + +A multidimensional Bernstein polynomial with coefficients of type `T` in `N` dimensions, +defined on a hyperrectangle `[low_corner[i], high_corner[i]]`. +""" +struct BernsteinPolynomial{N, T} <: Function + coeffs::Array{T, N} + low_corner::SVector{N, T} + high_corner::SVector{N, T} +end + +""" + BernsteinPolynomial(c::AbstractArray, lc, hc) + +Create a multidimensional [Bernstein +polynomial](https://en.wikipedia.org/wiki/Bernstein_polynomial#Generalizations_to_higher_dimension) +with coefficients `c` defined on the hyperrectangle `[lc[1], hc[1]] × … × [lc[N], hc[N]]`. + +Calling `p(x)` evaluates the polynomial at the point `x = (x[1], …, x[N])` by the formula + +```math +p(x_1,\\dots,x_D)=\\sum_{i_j=0}^{d_j}c_{i_1\\dots i_D}\\prod_{j=1}^D\\binom{d_j}{i_j}(x_j-l_j)^{i_j}(r_j-x_j)^{d_j-i_j} +``` + +where ``l_j = lc[j]`` and ``r_j = hc[j]`` are the lower and upper bounds of the +hyperrectangle, respectively, and ``d_j = size(c)[j] - 1`` is the degree of the polynomial +in dimension `j`. + +`BernsteinPolynomial`s can be differentiated using ForwardDiff; see [`gradient`](@ref) and +[`hessian`](@ref). +""" +function BernsteinPolynomial(c::AbstractArray, lc, hc) + N = ndims(c) + T = eltype(c) + coeffs = c isa Array{T, N} ? c : Array{T, N}(c) + return BernsteinPolynomial{N, T}(coeffs, SVector{N}(lc), SVector{N}(hc)) +end + +coefficients(p::BernsteinPolynomial) = p.coeffs +low_corner(p::BernsteinPolynomial) = p.low_corner +high_corner(p::BernsteinPolynomial) = p.high_corner +degree(p::BernsteinPolynomial) = size(coefficients(p)) .- 1 + +# evaluation +function (p::BernsteinPolynomial{N})(x) where {N} + x_ = SVector{N}(x) # try conversion to SVector + l = low_corner(p) + r = high_corner(p) + x₀ = (x_ - l) ./ (r - l) + c = coefficients(p) + return _evaluate_bernstein(x₀, c, Val{N}(), 1, length(c)) +end + +@fastmath function _evaluate_bernstein( + x::SVector{N}, + c::AbstractArray, + ::Val{dim}, + i1, + len, + ) where {N, dim} + n = size(c, dim) + @inbounds xd = x[dim] + # inspired by https://personal.math.ubc.ca/~cass/graphics/text/www/pdf/a6.pdf and the + # FastChebInterp.jl package + if dim == 1 + s = 1 - xd + @inbounds P = c[i1] + C = (n - 1) * xd + for k in 1:(n - 1) + @inbounds P = P * s + C * c[i1 + k] + C = C * (n - k - 1) / (k + 1) * xd + end + return P + else + Δi = len ÷ n # column-major stride of current dimension + + # we recurse downward on dim for cache locality, + # since earlier dimensions are contiguous + dim′ = Val{dim - 1}() + + s = 1 - xd + P = _evaluate_bernstein(x, c, dim′, i1, Δi) + C = (n - 1) * xd + for k in 1:(n - 1) + P = P * s + C * _evaluate_bernstein(x, c, dim′, i1 + k * Δi, Δi) + C = C * (n - k - 1) / (k + 1) * xd + end + return P + end +end + +""" + gradient(p, x) + +Compute the gradient of `p` at point `x` using ForwardDiff. +""" +function gradient(p, x) + return ForwardDiff.gradient(p, SVector{length(x)}(x)) +end + +""" + hessian(p, x) + +Compute the hessian of `p` at point `x` using ForwardDiff. +""" +function hessian(p, x) + return ForwardDiff.hessian(p, SVector{length(x)}(x)) +end + +function value_and_gradient(p, x) + res = ForwardDiff.gradient!(DiffResults.GradientResult(x), p, x) + return DiffResults.value(res), DiffResults.gradient(res) +end + +""" + value_gradient_hessian(p, x) + +Fused computation of the value, gradient, and hessian of `p` at point `x` using ForwardDiff. +""" +function value_gradient_hessian(p, x) + res = ForwardDiff.hessian!(DiffResults.HessianResult(x), p, x) + return DiffResults.value(res), DiffResults.gradient(res), DiffResults.hessian(res) +end + +# Batched right-multiply by Aᵀ: for r in 1:n_batch, C[:,:,r] = B[:,:,r] * Aᵀ, +# where C is (m, n) and B is (m, p), all stored flat in column-major order. +# C and B are accessed via linear indexing +function _rmult!(C, A::Matrix{T}, B, m, p, n, n_batch) where {T} + @inbounds for r in 1:n_batch + offset_B = (r - 1) * m * p + offset_C = (r - 1) * m * n + for j in 1:n + for i in 1:m + s = zero(T) + for k in 1:p + s += A[j, k] * B[offset_B + (k - 1) * m + i] + end + C[offset_C + (j - 1) * m + i] = s + end + end + end + return C +end + +# Apply the N-fold Kronecker product (mat ⊗ mat ⊗ … ⊗ mat) to vals using the vec trick, +# writing the result into coeffs. temp1 and temp2 are flat pre-allocated scratch buffers, +# each of length nc * nv^(N-1). +function _apply_kron!( + coeffs::AbstractArray{T, N}, + mat::Matrix{T}, + vals::Array{T, N}, + temp1::Vector{T}, + temp2::Vector{T}, + ) where {T, N} + nc, nv = size(mat) + N == 1 && return mul!(coeffs, mat, vals) + # Apply mat along each dimension sequentially via batched right-multiply by matᵀ. + # After step k, dims 1..k have size nc and dims k+1..N have size nv. + _rmult!(temp1, mat, vals, 1, nv, nc, nv^(N - 1)) + src, dst = temp1, temp2 + for k in 2:N + n_left = nc^(k - 1) + n_right = k < N ? nv^(N - k) : 1 + out = k < N ? dst : coeffs + _rmult!(out, mat, src, n_left, nv, nc, n_right) + src, dst = dst, src + end + return coeffs +end + +""" + mutable struct InterpolationData{N, T} + +Buffers and precomputed data for piecewise polynomial interpolation. Stores the 1D +interpolation matrix and scratch arrays needed to evaluate Bernstein coefficients on a +cell, but holds no reference to the underlying field. + +Construct via `InterpolationData(N, order, T)`. + +!!! warning "Not thread-safe" + The internal buffers (`coeffs`, `vals`, `temp1`, `temp2`) are shared mutable state. + Create one `InterpolationData` per thread if concurrent evaluation is needed. +""" +mutable struct InterpolationData{N, T} + order::Int + mat::Matrix{T} + coeffs::Array{T, N} + vals::Array{T, N} + Ic::CartesianIndex{N} # sentinel CartesianIndex(0,…,0) means "no cached cell" + temp1::Vector{T} + temp2::Vector{T} +end + +""" + InterpolationData(N::Int, order::Int, T::Type) + +Allocate an `InterpolationData` for `N`-dimensional fields with element type `T` and +polynomial interpolation of degree `order`. +""" +function InterpolationData(N::Int, order::Int, ::Type{T}) where {T} + # Build a 1D interpolation matrix mapping `order+1` values at equispaced nodes in [0,1] + # and returning the coefficients of the Bernstein basis defined on the interval + # [floor(order/2)/order, ceil(order/2)/order], which is symmetric around 0.5 and + # contains the central node at 0.5. When order is even, we use a stencil of size + # order+1 and do a least-squares fit. This matrix is used to compute the interpolant in + # a cell given values on a super-cell around it. + stencil_order = isodd(order) ? order : order + 1 + nc, nv = order + 1, stencil_order + 1 + nodes = ntuple(i -> (i - 1) / stencil_order, nv) + a, b = (stencil_order - 1) / (2 * stencil_order), (stencil_order + 1) / (2 * stencil_order) + B = (i, k, x) -> binomial(k, i) * (x^i) * ((1 - x)^(k - i)) + V = [B(j - 1, order, (nodes[i] - a) / (b - a)) for i in 1:nv, j in 1:nc] + mat = pinv(V) + coeffs = zeros(T, ntuple(_ -> nc, N)) + vals = zeros(T, ntuple(_ -> nv, N)) + buf_size = N > 1 ? nc * nv^(N - 1) : 0 + temp1 = zeros(T, buf_size) + temp2 = zeros(T, buf_size) + Ic = CartesianIndex(ntuple(_ -> 0, N)) + return InterpolationData{N, T}(order, mat, coeffs, vals, Ic, temp1, temp2) +end + +function Base.copy(itp::InterpolationData{N, T}) where {N, T} + return InterpolationData{N, T}( + itp.order, + copy(itp.mat), + copy(itp.coeffs), + copy(itp.vals), + CartesianIndex(ntuple(_ -> 0, Val(N))), # reset sentinel — new owner must rebuild + copy(itp.temp1), + copy(itp.temp2), + ) +end diff --git a/src/boundaryconditions.jl b/src/boundaryconditions.jl index 41daa34..cd8bbda 100644 --- a/src/boundaryconditions.jl +++ b/src/boundaryconditions.jl @@ -124,7 +124,7 @@ Base.show(io::IO, ::DirichletBC) = print(io, "Dirichlet") _normalize_bc(bc, dim) Normalize the `bc` argument into a `dim`-tuple of `(left, right)` pairs, one per spatial -dimension, as expected by [`add_boundary_conditions`](@ref). +dimension, as expected by [`_add_boundary_conditions`](@ref). - A single `BoundaryCondition` is applied to all sides. - A length-`dim` collection applies each entry to both sides of the corresponding dimension. diff --git a/src/interpolation.jl b/src/interpolation.jl index 3c028d9..7b58df0 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -1,204 +1,12 @@ -struct BernsteinPolynomial{N, T} <: Function - coeffs::Array{T, N} - low_corner::SVector{N, T} - high_corner::SVector{N, T} -end - -""" - BernsteinPolynomial(c::AbstractArray, lc, hc) - -Create a multidimensional [Bernstein -polynomial](https://en.wikipedia.org/wiki/Bernstein_polynomial#Generalizations_to_higher_dimension) -with coefficients `c` defined on the hyperrectangle `[lc[1], hc[1]] × … × [lc[N], hc[N]]`. - -Calling `p(x)` evaluates the polynomial at the point `x = (x[1], …, x[N])` by the formula - -```math -p(x_1,\\dots,x_D)=\\sum_{i_j=0}^{d_j}c_{i_1\\dots i_D}\\prod_{j=1}^D\\binom{d_j}{i_j}(x_j-l_j)^{i_j}(r_j-x_j)^{d_j-i_j} -``` - -where ``l_j = lc[j]`` and ``r_j = hc[j]`` are the lower and upper bounds of the -hyperrectangle, respectively, and ``d_j = size(c)[j] - 1`` is the degree of the polynomial -in dimension `j`. - -`BernsteinPolynomial`s can be differentiated using ForwardDiff; see [`gradient`](@ref) and -[`hessian`](@ref). -""" -function BernsteinPolynomial(c::AbstractArray, lc, hc) - N = ndims(c) - T = eltype(c) - coeffs = c isa Array{T, N} ? c : Array{T, N}(c) - return BernsteinPolynomial{N, T}(coeffs, SVector{N}(lc), SVector{N}(hc)) -end - -coefficients(p::BernsteinPolynomial) = p.coeffs -low_corner(p::BernsteinPolynomial) = p.low_corner -high_corner(p::BernsteinPolynomial) = p.high_corner -degree(p::BernsteinPolynomial) = size(coefficients(p)) .- 1 - -# evaluation -function (p::BernsteinPolynomial{N})(x) where {N} - x_ = SVector{N}(x) # try conversion to SVector - l = low_corner(p) - r = high_corner(p) - x₀ = (x_ - l) ./ (r - l) - c = coefficients(p) - return _evaluate_bernstein(x₀, c, Val{N}(), 1, length(c)) -end - -@fastmath function _evaluate_bernstein( - x::SVector{N}, - c::AbstractArray, - ::Val{dim}, - i1, - len, - ) where {N, dim} - n = size(c, dim) - @inbounds xd = x[dim] - # inspired by https://personal.math.ubc.ca/~cass/graphics/text/www/pdf/a6.pdf and the - # FastChebInterp.jl package - if dim == 1 - s = 1 - xd - @inbounds P = c[i1] - C = (n - 1) * xd - for k in 1:(n - 1) - @inbounds P = P * s + C * c[i1 + k] - C = C * (n - k - 1) / (k + 1) * xd - end - return P - else - Δi = len ÷ n # column-major stride of current dimension - - # we recurse downward on dim for cache locality, - # since earlier dimensions are contiguous - dim′ = Val{dim - 1}() - - s = 1 - xd - P = _evaluate_bernstein(x, c, dim′, i1, Δi) - C = (n - 1) * xd - for k in 1:(n - 1) - P = P * s + C * _evaluate_bernstein(x, c, dim′, i1 + k * Δi, Δi) - C = C * (n - k - 1) / (k + 1) * xd - end - return P - end -end - -""" - gradient(p, x) - -Compute the gradient of `p` at point `x` using ForwardDiff. """ -function gradient(p, x) - return ForwardDiff.gradient(p, SVector{length(x)}(x)) -end - -""" - hessian(p, x) - -Compute the hessian of `p` at point `x` using ForwardDiff. -""" -function hessian(p, x) - return ForwardDiff.hessian(p, SVector{length(x)}(x)) -end - - -function value_and_gradient(p, x) - res = ForwardDiff.gradient!(DiffResults.GradientResult(x), p, x) - return DiffResults.value(res), DiffResults.gradient(res) -end - -""" - value_gradient_hessian(p, x) - -Fused computation of the value, gradient, and hessian of `p` at point `x` using ForwardDiff. -""" -function value_gradient_hessian(p, x) - res = ForwardDiff.hessian!(DiffResults.HessianResult(x), p, x) - return DiffResults.value(res), DiffResults.gradient(res), DiffResults.hessian(res) -end - -""" - struct PiecewisePolynomialInterpolant{Φ, N, T} - -A piecewise polynomial interpolant built over a `MeshField`. - -The domain is partitioned into cells of the underlying mesh. On each cell a -`BernsteinPolynomial` is fitted to the surrounding stencil of grid values, so -that evaluating the interpolant at a point `x` amounts to: -1. locating the cell containing `x`, -2. (lazily) computing the local Bernstein coefficients, and -3. evaluating the resulting polynomial. - -Construct via [`interpolate`](@ref). The returned object `itp` supports: -- `itp(x)` — evaluate at point `x` -- [`make_interpolant`](@ref)`(itp, I)` — local `BernsteinPolynomial` for cell `I` -- [`gradient`](@ref)`(itp, x)` / [`hessian`](@ref)`(itp, x)` — via ForwardDiff -- [`cell_extrema`](@ref)`(itp, I)` — certified bounds via convex-hull property -- [`proven_empty`](@ref)`(itp, I)` — certified absence of interface or interior - -The struct is mutable because the local coefficients and stencil values are -cached in place (`coeffs`, `vals`) to avoid allocation on repeated evaluations. + compute_index(ϕ::AbstractMeshField, x) -!!! warning "Not thread-safe" - The internal buffers (`coeffs`, `vals`, `temp1`, `temp2`, `Ic`) are shared - state. Concurrent evaluation from multiple threads will cause data races. - Create one interpolant per thread if parallelism is needed. +Compute the multi-index of the cell containing point `x`. If `x` is outside +the domain, the index of the closest cell is returned (clamped to the boundary). """ -mutable struct PiecewisePolynomialInterpolant{Φ, N, T} - ϕ::Φ # underlying mesh field - mat::Matrix{T} # map from grid to bernstein vals in 1D - coeffs::Array{T, N} # buffer for bernstein coeffs - vals::Array{T, N} # buffer for grid values in the stencil - Ic::CartesianIndex{N} # multi-index of the currently loaded cell - # scratch buffers for applying the N-fold kron mat to vals - temp1::Vector{T} - temp2::Vector{T} -end - -function PiecewisePolynomialInterpolant(ϕ::Φ, order::Int) where {Φ} +function compute_index(ϕ::AbstractMeshField, x) grid = mesh(ϕ) - N = ndims(grid) - T = eltype(ϕ) - - # Build a 1D interpolation matrix mapping `order+1` values at equispaced nodes in [0,1] - # and returning the coefficients of the Bernstein basis defined on the interval - # [floor(order/2)/order, ceil(order/2)/order], which is symmetric around 0.5 and - # contains the central node at 0.5. When order is even, we use a stencil of size - # order+1 and do a least-squares fit. This matrix is used to compute the interpolant in - # a cell given values on a super-cell around it. - stencil_order = isodd(order) ? order : order + 1 - nc, nv = order + 1, stencil_order + 1 - nodes = ntuple(i -> (i - 1) / stencil_order, Val(nv)) - a, b = (stencil_order - 1) / (2 * stencil_order), (stencil_order + 1) / (2 * stencil_order) - B = (i, k, x) -> binomial(k, i) * (x^i) * ((1 - x)^(k - i)) - V = [B(j - 1, order, (nodes[i] - a) / (b - a)) for i in 1:nv, j in 1:nc] - - mat = pinv(V) - coeffs = zeros(T, ntuple(_ -> nc, N)) - vals = zeros(T, ntuple(_ -> nv, N)) - - # Intermediate buffers: each needs nc * nv^(N-1) elements (the largest intermediate size) - buf_size = N > 1 ? nc * nv^(N - 1) : 0 - temp1 = zeros(T, buf_size) - temp2 = zeros(T, buf_size) - - Ic = CartesianIndex(ntuple(_ -> 0, Val(N))) - - return PiecewisePolynomialInterpolant(ϕ, mat, coeffs, vals, Ic, temp1, temp2) -end - -Base.ndims(::PiecewisePolynomialInterpolant{Φ, N}) where {Φ, N} = N - -""" - compute_index(itp::PiecewisePolynomialInterpolant, x) - -Compute the multi-index of the cell containing point `x`. If `x` is outside the domain, the -index of the closest cell is returned (clamping to the boundary). -""" -function compute_index(itp::PiecewisePolynomialInterpolant, x) - grid = mesh(itp.ϕ) - N = ndims(itp) + N = ndims(ϕ) cell_ax = cellindices(grid) return ntuple( d -> clamp( @@ -209,177 +17,90 @@ function compute_index(itp::PiecewisePolynomialInterpolant, x) ) |> CartesianIndex end -# Batched right-multiply by Aᵀ: for r in 1:n_batch, C[:,:,r] = B[:,:,r] * Aᵀ, -# where C is (m, n) and B is (m, p), all stored flat in column-major order. -# C and B are accessed via linear indexing -function _rmult!(C, A::Matrix{T}, B, m, p, n, n_batch) where {T} - @inbounds for r in 1:n_batch - offset_B = (r - 1) * m * p - offset_C = (r - 1) * m * n - for j in 1:n - for i in 1:m - s = zero(T) - for k in 1:p - s += A[j, k] * B[offset_B + (k - 1) * m + i] - end - C[offset_C + (j - 1) * m + i] = s - end - end - end - return C -end - -# Apply the N-fold Kronecker product (mat ⊗ mat ⊗ … ⊗ mat) to vals using the vec trick, -# writing the result into coeffs. temp1 and temp2 are flat pre-allocated scratch buffers, -# each of length nc * nv^(N-1). -function _apply_kron!( - coeffs::AbstractArray{T, N}, - mat::Matrix{T}, - vals::Array{T, N}, - temp1::Vector{T}, - temp2::Vector{T}, - ) where {T, N} - nc, nv = size(mat) - N == 1 && return mul!(coeffs, mat, vals) - # Apply mat along each dimension sequentially via batched right-multiply by matᵀ. - # After step k, dims 1..k have size nc and dims k+1..N have size nv. - _rmult!(temp1, mat, vals, 1, nv, nc, nv^(N - 1)) - src, dst = temp1, temp2 - for k in 2:N - n_left = nc^(k - 1) - n_right = k < N ? nv^(N - k) : 1 - out = k < N ? dst : coeffs - _rmult!(out, mat, src, n_left, nv, nc, n_right) - src, dst = dst, src - end - return coeffs -end - """ - fill_coefficients!(itp::PiecewisePolynomialInterpolant, base_idxs::CartesianIndex) + fill_coefficients!(ϕ::AbstractMeshField, base_idxs::CartesianIndex) -Fill the internal buffer of `itp` with the Bernstein coefficients for the cell at `base_idxs`. +Fill the internal `InterpolationData` buffer of `ϕ` with the Bernstein coefficients +for the cell at `base_idxs`. """ -@inline function fill_coefficients!( - itp::PiecewisePolynomialInterpolant{Φ, N, T}, - base_idxs::CartesianIndex{N}, - ) where {Φ, N, T} - ϕ = itp.ϕ +@inline function fill_coefficients!(ϕ::AbstractMeshField, base_idxs::CartesianIndex{N}) where {N} + itp = interp_data(ϕ) mat = itp.mat nc, nn = size(mat) KS = nn - 1 # order of the interpolation stencil off = -(KS - 1) ÷ 2 - # Gather grid values into vals. Since the field ϕ may generate values on demand via BCs, - # we can't just copy or have a view + # Gather grid values into vals. Since ϕ may generate values on demand via BCs, + # we can't just copy or have a view. for I in CartesianIndices(itp.vals) J = CartesianIndex(ntuple(d -> base_idxs[d] + off + I[d] - 1, N)) @inbounds itp.vals[I] = ϕ[J] end - # The Vandermonde matrix is a Kronecker product V = V₁ ⊗ … ⊗ V₁ _apply_kron!(itp.coeffs, mat, itp.vals, itp.temp1, itp.temp2) itp.Ic = base_idxs - return itp + return ϕ end """ - make_interpolant(itp::PiecewisePolynomialInterpolant, I::CartesianIndex) + make_interpolant(ϕ::AbstractMeshField, I::CartesianIndex) + +Return a `BernsteinPolynomial` for the cell at multi-index `I`, lazily computing +(and caching) its Bernstein coefficients from the surrounding stencil of grid values. -Create a `BernsteinPolynomial` for the cell at multi-index `I`. +!!! warning "Aliased coefficients" + The returned polynomial's `coeffs` array aliases the internal buffer of `ϕ`'s + `InterpolationData`. It remains valid only until the next `make_interpolant` call on + `ϕ` with a different cell index. Copy `coefficients(p)` if the polynomial must outlive + the current cell iteration. """ -function make_interpolant(itp::PiecewisePolynomialInterpolant{Φ, N}, I::CartesianIndex{N}) where {Φ, N} - I == itp.Ic || fill_coefficients!(itp, I) - cell = getcell(mesh(itp.ϕ), I) +function make_interpolant(ϕ::AbstractMeshField, I::CartesianIndex) + itp = interp_data(ϕ) + I == itp.Ic || fill_coefficients!(ϕ, I) + cell = _getcell(mesh(ϕ), I) return BernsteinPolynomial(itp.coeffs, cell.lc, cell.hc) end -function Base.show(io::IO, ::MIME"text/plain", itp::PiecewisePolynomialInterpolant) - order = size(itp.mat, 1) - 1 - N = ndims(itp) - println(io, "PiecewisePolynomialInterpolant") - println(io, " ├─ order: $order") - println(io, " └─ field: MeshField on CartesianGrid in ℝ$(_superscript(N))") - return _show_fields(io, itp.ϕ; prefix = " ") -end - -@inline _evaluate(p::P, x) where {P} = p(x) - -@inline function (itp::PiecewisePolynomialInterpolant)(x) - I = compute_index(itp, x) - p = make_interpolant(itp, I) - return _evaluate(p, x) -end - -@inline (itp::PiecewisePolynomialInterpolant)(x::Vararg{Real}) = itp(SVector(x)) -@inline (itp::PiecewisePolynomialInterpolant)(x::Tuple) = itp(SVector(x)) - """ - interpolate(ϕ::MeshField, order::Int = 3) - -Create a piecewise polynomial interpolant of the given `order` for `ϕ`. - -A deep copy of `ϕ` is made so that the interpolant is independent of future modifications to -`ϕ`. If `ϕ` has no boundary conditions, `ExtrapolationBC{2}` is added automatically on all -sides (this is necessary to evaluate the interpolant near the boundary, where the stencil -may require out-of-bounds values). + cell_extrema(ϕ::AbstractMeshField, I::CartesianIndex) -The returned object `itp` behaves like a function and supports: -- `itp(x)`: evaluate the interpolant at point `x` -- [`make_interpolant`](@ref)`(itp, I)`: return the local interpolant for cell `I` -- [`gradient`](@ref)`(itp, x)`: gradient at `x` (via `make_interpolant` + ForwardDiff) -- [`hessian`](@ref)`(itp, x)`: hessian at `x` (via `make_interpolant` + ForwardDiff) -- [`cell_extrema`](@ref)`(itp, I)`: lower and upper bounds of the interpolant in cell `I` - -To avoid unnecessary copying, call `update!(itp, ϕ)` to update the interpolant's internal -field with new values from `ϕ` while reusing the existing interpolation matrix and scratch -buffers; see [`update!`](@ref) for details. +Compute the minimum and maximum values of the local Bernstein interpolant in cell `I`. """ -function interpolate(ϕ, order::Int = 3) - ϕ_copy = deepcopy(ϕ) - if !has_boundary_conditions(ϕ_copy) - N = ndims(mesh(ϕ_copy)) - bc = ntuple(_ -> (ExtrapolationBC{2}(), ExtrapolationBC{2}()), N) - ϕ_copy = add_boundary_conditions(ϕ_copy, bc) - end - return PiecewisePolynomialInterpolant(ϕ_copy, order) +function cell_extrema(ϕ::AbstractMeshField, I::CartesianIndex) + p = make_interpolant(ϕ, I) + return extrema(coefficients(p)) end """ - update!(itp::PiecewisePolynomialInterpolant, ϕ) + proven_empty(ϕ::AbstractMeshField, I::CartesianIndex; surface=false) -Copy the values of `ϕ` into the interpolant's internal field and invalidate the cell -cache. This is cheaper than calling [`interpolate`](@ref) again because it reuses the -existing interpolation matrix and scratch buffers. +Return `true` if cell `I` is guaranteed to contain no interface (when +`surface=true`) or no interior (when `surface=false`), based on the convex-hull +property of the Bernstein basis. """ -function update!(itp::PiecewisePolynomialInterpolant{Φ, N}, ϕ) where {Φ, N} - copy!(itp.ϕ, ϕ) - itp.Ic = CartesianIndex(ntuple(_ -> 0, Val(N))) - return itp +function proven_empty(ϕ::AbstractMeshField, I::CartesianIndex; surface = false) + m, M = cell_extrema(ϕ, I) + return surface ? (m * M > 0) : (m > 0) end -""" - cell_extrema(itp::PiecewisePolynomialInterpolant, I::CartesianIndex) - -Compute the minimum and maximum values of the interpolant in the cell `I`. -""" -function cell_extrema(itp::PiecewisePolynomialInterpolant{Φ, N}, I::CartesianIndex{N}) where {Φ, N} - p = make_interpolant(itp, I) - return extrema(coefficients(p)) +# Make a MeshField with InterpolationData callable as a continuous function. +@inline function (ϕ::MeshField{V, M, B, <:InterpolationData})(x) where {V, M, B} + I = compute_index(ϕ, x) + p = make_interpolant(ϕ, I) + return p(x) end -""" - proven_empty(itp::PiecewisePolynomialInterpolant, I::CartesianIndex; surface=false) - -Return `true` if the cell `I` is guaranteed to not contain the interface (if `surface=true`) -or to not contain any part of the interior (if `surface=false`). +@inline (ϕ::MeshField{V, M, B, <:InterpolationData})(x::Vararg{Real}) where {V, M, B} = + ϕ(SVector(x)) +@inline (ϕ::MeshField{V, M, B, <:InterpolationData})(x::Tuple) where {V, M, B} = + ϕ(SVector(x)) -Note that `proven_empty` being `false` does not mean the cell is non-empty, but rather that -we can't guarantee emptiness based on the convex hull property of the Bernstein basis. -""" -function proven_empty(itp::PiecewisePolynomialInterpolant, I::CartesianIndex; surface = false) - m, M = cell_extrema(itp, I) - if surface - return m * M > 0 - else - return m > 0 - end +# Make a NarrowBandMeshField with InterpolationData callable as a continuous function. +@inline function (nb::NarrowBandMeshField{V, M, B, T, <:InterpolationData})(x) where {V, M, B, T} + I = compute_index(nb, x) + p = make_interpolant(nb, I) + return p(x) end + +@inline (nb::NarrowBandMeshField{V, M, B, T, <:InterpolationData})(x::Vararg{Real}) where {V, M, B, T} = + nb(SVector(x)) +@inline (nb::NarrowBandMeshField{V, M, B, T, <:InterpolationData})(x::Tuple) where {V, M, B, T} = + nb(SVector(x)) diff --git a/src/levelset.jl b/src/levelset.jl index d9d9b4f..4c3e9b9 100644 --- a/src/levelset.jl +++ b/src/levelset.jl @@ -1,32 +1,14 @@ -""" - const LevelSet{N, T, B} - -Alias for [`MeshField`](@ref) on a `CartesianGrid{N,T}` with values stored as -an `Array{T,N}` and a [`FullDomain`](@ref). `B` is the type of the boundary -conditions. -""" -const LevelSet{N, T, B} = - MeshField{Array{T, N}, CartesianGrid{N, T}, B, FullDomain} - -LevelSet(f::Function, grid, bc = nothing) = MeshField(f, grid, bc) - -active_indices(ϕ::LevelSet) = CartesianIndices(mesh(ϕ)) - -""" - _ensure_boundary_conditions(ϕ) - -Return `ϕ` unchanged if it already carries boundary conditions, otherwise wrap -it with `LinearExtrapolationBC` on every face. -""" +# Return `ϕ` unchanged if it already carries boundary conditions, otherwise wrap +# it with `LinearExtrapolationBC` on every face (needed for finite-difference stencils). function _ensure_boundary_conditions(ϕ) has_boundary_conditions(ϕ) && return ϕ N = ndims(ϕ) bc = _normalize_bc(LinearExtrapolationBC(), N) - return add_boundary_conditions(ϕ, bc) + return _add_boundary_conditions(ϕ, bc) end """ - volume(ϕ::LevelSet) + volume(ϕ::MeshField) Compute the volume of the level-set function. @@ -43,7 +25,7 @@ LevelSetMethods.volume(ϕ), V0 (0.7854362890190668, 0.7853981633974483) ``` """ -function volume(ϕ::LevelSet) +function volume(ϕ::MeshField) δ = meshsize(mesh(ϕ)) δmin = minimum(δ) vol = prod(δ) @@ -52,7 +34,7 @@ function volume(ϕ::LevelSet) end """ - perimeter(ϕ::LevelSet) + perimeter(ϕ::MeshField) Compute the perimeter area of the level-set function. @@ -68,10 +50,10 @@ LevelSetMethods.perimeter(ϕ), S0 # output -(3.1426415491430366, 3.141592653589793) +(3.142641549143036, 3.141592653589793) ``` """ -function perimeter(ϕ::LevelSet) +function perimeter(ϕ::MeshField) ϕ = _ensure_boundary_conditions(ϕ) δ = meshsize(mesh(ϕ)) δmin = minimum(δ) @@ -99,7 +81,7 @@ end # Method used to numerically integrate N-th dimensional matrices. # return a matrix with coefficients equal to 2^(-n) where n is # the number of times this index if on the borders of a dimension. -function trapezoidal_coefficients(ϕ::LevelSet) +function trapezoidal_coefficients(ϕ::MeshField) N = ndims(ϕ) ax = axes(ϕ) coeffs = zeros(size(values(ϕ))) @@ -111,7 +93,7 @@ function trapezoidal_coefficients(ϕ::LevelSet) end """ - curvature(ϕ::LevelSet, I) + curvature(ϕ::MeshField, I) Compute the mean curvature of ϕ at I using κ = ∇ ⋅ (∇ϕ / |∇ϕ|). We use the formula κ = (Δϕ |∇ϕ|^2 - ∇ϕ^T Hϕ ∇ϕ) / |∇ϕ|^3 with @@ -153,10 +135,10 @@ function curvature(ϕ, I) end """ - curvature(ϕ::LevelSet) + curvature(ϕ::MeshField) Compute the mean curvature of ϕ at I using κ = ∇ ⋅ (∇ϕ / |∇ϕ|). -See [`curvature(ϕ::LevelSet, I)`](@ref) for more details. +See [`curvature(ϕ::MeshField, I)`](@ref) for more details. ```julia using LevelSetMethods @@ -173,13 +155,13 @@ Colorbar(fig[:, end+1], hm) contour!(xs, ys, values(ϕ); levels = [0.0]) ``` """ -function curvature(ϕ::LevelSet) - ϕ_ = _ensure_boundary_conditions(ϕ) - return [curvature(ϕ_, I) for I in eachindex(ϕ_)] +function curvature(ϕ::MeshField) + _check_bc(ϕ) + return [curvature(ϕ, I) for I in eachindex(ϕ)] end """ - gradient(ϕ::LevelSet, I::CartesianIndex) + gradient(ϕ::MeshField, I::CartesianIndex) Compute the gradient vector ``∇ϕ`` at grid index `I` using centered finite differences. Returns an `SVector` (or `Vector`) of derivatives. @@ -190,16 +172,16 @@ function gradient(ϕ, I::CartesianIndex) end """ - gradient(ϕ::LevelSet) + gradient(ϕ::MeshField) Compute the gradient vector ``∇ϕ`` for all grid points. """ -function gradient(ϕ::LevelSet) +function gradient(ϕ::MeshField) return [gradient(ϕ, I) for I in eachindex(ϕ)] end """ - normal(ϕ::LevelSet, I::CartesianIndex) + normal(ϕ::MeshField, I::CartesianIndex) Compute the unit exterior normal vector ``\\mathbf{n} = \\frac{∇ϕ}{\\|∇ϕ\\|}`` at grid index `I`. """ @@ -209,7 +191,7 @@ function normal(ϕ, I) end """ - normal(ϕ::LevelSet) + normal(ϕ::MeshField) Compute the unit exterior normal vector ``\\mathbf{n} = \\frac{∇ϕ}{\\|∇ϕ\\|}`` for all grid points. @@ -229,13 +211,13 @@ arrows(xs, ys, us, vs; arrowsize = 10 * vec(coeff), lengthscale = 2.0 / (N - 1)) contour!(xs, ys, values(ϕ); levels = [0.0]) ``` """ -function normal(ϕ::LevelSet) - ϕ_ = _ensure_boundary_conditions(ϕ) - return [normal(ϕ_, I) for I in eachindex(ϕ_)] +function normal(ϕ::MeshField) + _check_bc(ϕ) + return [normal(ϕ, I) for I in eachindex(ϕ)] end """ - hessian(ϕ::LevelSet, I::CartesianIndex) + hessian(ϕ::MeshField, I::CartesianIndex) Compute the Hessian matrix ``\\mathbf{H}ϕ = ∇∇ϕ`` at grid index `I` using second-order finite differences. Returns a `Symmetric` matrix. @@ -246,25 +228,22 @@ function hessian(ϕ, I::CartesianIndex) end """ - hessian(ϕ::LevelSet) + hessian(ϕ::MeshField) Compute the Hessian matrix for all grid points. """ -function hessian(ϕ::LevelSet) +function hessian(ϕ::MeshField) return [hessian(ϕ, I) for I in eachindex(ϕ)] end """ - grad_norm(ϕ::LevelSet) + grad_norm(ϕ::MeshField) Compute the norm of the gradient of ϕ, i.e. `|∇ϕ|`, at all grid points. """ -function grad_norm(ϕ::LevelSet) - msg = """level-set must have boundary conditions to compute gradient. See - `add_boundary_conditions`.""" - has_boundary_conditions(ϕ) || error(msg) - idxs = eachindex(ϕ) - return map(i -> _compute_∇_norm(sign(ϕ[i]), ϕ, i), idxs) +function grad_norm(ϕ::MeshField) + _check_bc(ϕ) + return map(i -> _compute_∇_norm(sign(ϕ[i]), ϕ, i), eachindex(ϕ)) end #= @@ -277,49 +256,49 @@ are optional keyword arguments (e.g. the `center` or `radius` of a circle). =# """ - circle(grid; center = (0, 0), radius = 1) + circle(grid; center = (0, 0), radius = 1, kwargs...) Create a 2D circle with the specified `center` and `radius` on a `grid`. -Returns a [`LevelSet`](@ref) field. +Returns a [`MeshField`](@ref) field. """ -function circle(grid; center = (0, 0), radius = 1) +function circle(grid; center = (0, 0), radius = 1, kwargs...) ndims(grid) == 2 || throw(ArgumentError("circle shape is only available in two dimensions")) - return LevelSet(x -> sqrt(sum((x .- center) .^ 2)) - radius, grid) + return MeshField(x -> sqrt(sum((x .- center) .^ 2)) - radius, grid; kwargs...) end """ - rectangle(grid; center = (0, 0), width = (1, 1)) + rectangle(grid; center = (0, 0), width = (1, 1), kwargs...) Create a rectangle (or N-dimensional box) with the specified `center` and `width` on a `grid`. -Returns a [`LevelSet`](@ref) field. +Returns a [`MeshField`](@ref) field. """ -function rectangle(grid; center = zero(grid.lc), width = (1, 1)) - return LevelSet(x -> maximum(abs.(x .- center) .- width ./ 2), grid) +function rectangle(grid; center = zero(grid.lc), width = (1, 1), kwargs...) + return MeshField(x -> maximum(abs.(x .- center) .- width ./ 2), grid; kwargs...) end """ - sphere(grid; center = (0, 0, 0), radius) + sphere(grid; center = (0, 0, 0), radius, kwargs...) Create a 3D sphere with the specified `center` and `radius` on a `grid`. -Returns a [`LevelSet`](@ref) field. +Returns a [`MeshField`](@ref) field. """ -function sphere(grid; center = (0, 0, 0), radius) +function sphere(grid; center = (0, 0, 0), radius, kwargs...) ndims(grid) == 3 || throw(ArgumentError("sphere shape is only available in three dimensions")) - return LevelSet(x -> sqrt(sum((x .- center) .^ 2)) - radius, grid) + return MeshField(x -> sqrt(sum((x .- center) .^ 2)) - radius, grid; kwargs...) end """ - star(grid; radius = 1, deformation = 0.25, n = 5.0) + star(grid; radius = 1, deformation = 0.25, n = 5.0, kwargs...) Create a 2D star shape defined in polar coordinates by ``r = R(1 + d \\cos(nθ))``. -Returns a [`LevelSet`](@ref) field. +Returns a [`MeshField`](@ref) field. """ -function star(grid; radius = 1, deformation = 0.25, n = 5.0) +function star(grid; radius = 1, deformation = 0.25, n = 5.0, kwargs...) # ndims(grid) == 2 || # throw(ArgumentError("star shape is only available in two dimensions")) - return LevelSet(grid) do x + return MeshField(grid; kwargs...) do x r = norm(x) θ = atan(x[2], x[1]) return r - radius * (1.0 + deformation * cos(n * θ)) @@ -327,29 +306,29 @@ function star(grid; radius = 1, deformation = 0.25, n = 5.0) end """ - dumbbell(grid; width = 1, height = 0.2, radius = 0.25, center = (0, 0)) + dumbbell(grid; width = 1, height = 0.2, radius = 0.25, center = (0, 0), kwargs...) Create a 2D dumbbell shape consisting of two circles connected by a rectangle. -Returns a [`LevelSet`](@ref) field. +Returns a [`MeshField`](@ref) field. """ -function dumbbell(grid; width = 1, height = 1 / 5, radius = 1 / 4, center = (0, 0)) - cl = circle(grid; center = center .- (width / 2, 0), radius) - cr = circle(grid; center = center .+ (width / 2, 0), radius) - rec = rectangle(grid; center, width = (width, height)) +function dumbbell(grid; width = 1, height = 1 / 5, radius = 1 / 4, center = (0, 0), kwargs...) + cl = circle(grid; center = center .- (width / 2, 0), radius, kwargs...) + cr = circle(grid; center = center .+ (width / 2, 0), radius, kwargs...) + rec = rectangle(grid; center, width = (width, height), kwargs...) return cl ∪ cr ∪ rec end """ - zalesak_disk(grid; center = (0, 0), radius = 0.5, width = 0.25, height = 1) + zalesak_disk(grid; center = (0, 0), radius = 0.5, width = 0.25, height = 1, kwargs...) Create a Zalesak disk (a circle with a rectangular slot cut out). -Used for testing advection schemes. Returns a [`LevelSet`](@ref) field. +Used for testing advection schemes. Returns a [`MeshField`](@ref) field. """ -function zalesak_disk(grid; center = (0, 0), radius = 0.5, width = 0.25, height = 1) +function zalesak_disk(grid; center = (0, 0), radius = 0.5, width = 0.25, height = 1, kwargs...) ndims(grid) == 2 || throw(ArgumentError("zalesak disk shape is only available in two dimensions")) - disk = circle(grid; center = center, radius = radius) - rec = rectangle(grid; center = center .- (0, radius), width = (width, height)) + disk = circle(grid; center = center, radius = radius, kwargs...) + rec = rectangle(grid; center = center .- (0, radius), width = (width, height), kwargs...) return setdiff(disk, rec) end @@ -360,73 +339,73 @@ Set operations for level set functions. =# """ - union!(ϕ1::LevelSet, ϕ2::LevelSet) + union!(ϕ1::MeshField, ϕ2::MeshField) In-place union of two level sets: ``ϕ_1 = \\min(ϕ_1, ϕ_2)``. """ -function Base.union!(ϕ1::LevelSet, ϕ2::LevelSet) +function Base.union!(ϕ1::MeshField, ϕ2::MeshField) v1, v2 = values(ϕ1), values(ϕ2) v1 .= min.(v1, v2) return ϕ1 end """ - union(ϕ1::LevelSet, ϕ2::LevelSet) + union(ϕ1::MeshField, ϕ2::MeshField) Return the union of two level sets: ``\\min(ϕ_1, ϕ_2)``. """ -Base.union(ϕ1::LevelSet, ϕ2::LevelSet) = union!(deepcopy(ϕ1), ϕ2) +Base.union(ϕ1::MeshField, ϕ2::MeshField) = union!(deepcopy(ϕ1), ϕ2) """ - intersect!(ϕ1::LevelSet, ϕ2::LevelSet) + intersect!(ϕ1::MeshField, ϕ2::MeshField) In-place intersection of two level sets: ``ϕ_1 = \\max(ϕ_1, ϕ_2)``. """ -function Base.intersect!(ϕ1::LevelSet, ϕ2::LevelSet) +function Base.intersect!(ϕ1::MeshField, ϕ2::MeshField) v1, v2 = values(ϕ1), values(ϕ2) v1 .= max.(v1, v2) return ϕ1 end """ - intersect(ϕ1::LevelSet, ϕ2::LevelSet) + intersect(ϕ1::MeshField, ϕ2::MeshField) Return the intersection of two level sets: ``\\max(ϕ_1, ϕ_2)``. """ -Base.intersect(ϕ1::LevelSet, ϕ2::LevelSet) = intersect!(deepcopy(ϕ1), ϕ2) +Base.intersect(ϕ1::MeshField, ϕ2::MeshField) = intersect!(deepcopy(ϕ1), ϕ2) """ - complement!(ϕ::LevelSet) + complement!(ϕ::MeshField) In-place complement of a level set (negates the values). """ -function complement!(ϕ::LevelSet) +function complement!(ϕ::MeshField) v = values(ϕ) v .= -v return ϕ end """ - complement(ϕ::LevelSet) + complement(ϕ::MeshField) Return the complement of a level set (negates the values). """ -complement(ϕ::LevelSet) = complement!(deepcopy(ϕ)) +complement(ϕ::MeshField) = complement!(deepcopy(ϕ)) """ - setdiff!(ϕ1::LevelSet, ϕ2::LevelSet) + setdiff!(ϕ1::MeshField, ϕ2::MeshField) In-place set difference: ``ϕ_1 = \\max(ϕ_1, -ϕ_2)``. """ -function Base.setdiff!(ϕ1::LevelSet, ϕ2::LevelSet) +function Base.setdiff!(ϕ1::MeshField, ϕ2::MeshField) v1, v2 = values(ϕ1), values(ϕ2) v1 .= max.(v1, -v2) return ϕ1 end """ - setdiff(ϕ1::LevelSet, ϕ2::LevelSet) + setdiff(ϕ1::MeshField, ϕ2::MeshField) Return the set difference: ``\\max(ϕ_1, -ϕ_2)``. """ -Base.setdiff(ϕ1::LevelSet, ϕ2::LevelSet) = setdiff!(deepcopy(ϕ1), ϕ2) +Base.setdiff(ϕ1::MeshField, ϕ2::MeshField) = setdiff!(deepcopy(ϕ1), ϕ2) diff --git a/src/levelsetequation.jl b/src/levelsetequation.jl index d64504e..a9e23c5 100644 --- a/src/levelsetequation.jl +++ b/src/levelsetequation.jl @@ -1,7 +1,7 @@ mutable struct LevelSetEquation - terms::Tuple + terms::Tuple{Vararg{LevelSetTerm}} integrator::TimeIntegrator - state::MeshField + state::AbstractMeshField t::Float64 reinit::Union{Nothing, NewtonReinitializer} log::SimulationLog @@ -11,8 +11,8 @@ end LevelSetEquation(; terms, ic, bc, t = 0, integrator = RK2(), reinit = nothing) Create a level-set equation of the form `ϕₜ + sum(terms) = 0`, where each `t ∈ terms` -is a [`LevelSetTerm`](@ref) and `ic` is the initial condition — either a [`LevelSet`](@ref) -for a full-grid discretization or a [`NarrowBandLevelSet`](@ref) for a narrow-band one. +is a [`LevelSetTerm`](@ref) and `ic` is the initial condition — either a [`MeshField`](@ref) +for a full-grid discretization or a [`NarrowBandMeshField`](@ref) for a narrow-band discretization. Calling [`integrate!(eq, tf)`](@ref) will evolve the equation up to time `tf`, modifying `current_state(eq)` and `current_time(eq)` in place. @@ -36,7 +36,7 @@ Reinitialization is controlled by the `reinit` keyword, which accepts: ```jldoctest; output = true using LevelSetMethods, StaticArrays grid = CartesianGrid((-1, -1), (1, 1), (50, 50)) # define the grid -ϕ = LevelSet(x -> x[1]^2 + x[2]^2 - 0.5^2, grid) # initial shape +ϕ = MeshField(x -> x[1]^2 + x[2]^2 - 0.5^2, grid) # initial shape 𝐮 = MeshField(x -> SVector(1, 0), grid) # advection velocity terms = (AdvectionTerm(𝐮),) # advection term bc = NeumannBC() # zero-gradient boundary conditions @@ -65,7 +65,7 @@ LevelSetEquation function LevelSetEquation(; terms, integrator = RK2(), - ic::MeshField, + ic::AbstractMeshField, bc, t = 0, reinit = nothing, @@ -76,7 +76,7 @@ function LevelSetEquation(; # bc is the authoritative source has_boundary_conditions(ic) && @warn "ic already has boundary conditions; these will be overwritten by bc" - state = add_boundary_conditions(ic, bc) + state = _add_boundary_conditions(ic, bc) log = SimulationLog(t, terms) return LevelSetEquation(terms, integrator, state, t, reinit, log) end @@ -150,12 +150,12 @@ end """ current_state(eq::LevelSetEquation) -Return the current state of the level-set equation (a [`LevelSet`](@ref)). +Return the current state of the level-set equation (a [`MeshField`](@ref)). """ current_state(ls::LevelSetEquation) = ls.state -# Allow current_state on a bare MeshField so that plotting recipes work uniformly. -current_state(ϕ::MeshField) = ϕ +# Allow current_state on a bare AbstractMeshField so that plotting recipes work uniformly. +current_state(ϕ::AbstractMeshField) = ϕ """ current_time(eq::LevelSetEquation) @@ -200,7 +200,6 @@ Return the [`NewtonReinitializer`](@ref) (if any) attached to the equation. reinitializer(ls::LevelSetEquation) = ls.reinit # Convenience delegations to current state -interpolate(eq::LevelSetEquation, order::Int = 3) = interpolate(current_state(eq), order) volume(eq::LevelSetEquation) = volume(current_state(eq)) perimeter(eq::LevelSetEquation) = perimeter(current_state(eq)) diff --git a/src/levelsetterms.jl b/src/levelsetterms.jl index 4f0b8ef..f3bad43 100644 --- a/src/levelsetterms.jl +++ b/src/levelsetterms.jl @@ -14,7 +14,7 @@ called at each stage of the time integration. update_term!(::LevelSetTerm, _, _) = nothing """ - compute_cfl(terms, ϕ::LevelSet, t) + compute_cfl(terms, ϕ::MeshField, t) Compute the maximum stable time-step ``Δt`` for the given `terms` and level set `ϕ` at time `t`, based on the Courant-Friedrichs-Lewy (CFL) condition. @@ -65,13 +65,13 @@ end Base.show(io::IO, _::AdvectionTerm) = print(io, "𝐮 ⋅ ∇ ϕ") -@inline function _compute_term(term::AdvectionTerm{V}, ϕ::MeshField, I, t) where {V} +@inline function _compute_term(term::AdvectionTerm{V}, ϕ::AbstractMeshField, I, t) where {V} sch = scheme(term) N = ndims(ϕ) 𝐮 = if V <: MeshField velocity(term)[I] elseif V <: Function - x = mesh(ϕ)[I] + x = getnode(mesh(ϕ), I) velocity(term)(x, t) else error("velocity field type $V not supported") @@ -105,7 +105,7 @@ function _compute_cfl(term::AdvectionTerm{V}, ϕ, I, t) where {V} 𝐮 = if V <: MeshField velocity(term)[I] elseif V <: Function - x = mesh(ϕ)[I] + x = getnode(mesh(ϕ), I) velocity(term)(x, t) else error("velocity field type $V not supported") @@ -127,14 +127,14 @@ coefficient(cterm::CurvatureTerm) = cterm.b Base.show(io::IO, _::CurvatureTerm) = print(io, "b κ|∇ϕ|") -function _compute_term(term::CurvatureTerm, ϕ::MeshField, I, t) +function _compute_term(term::CurvatureTerm, ϕ::AbstractMeshField, I, t) N = ndims(ϕ) κ = curvature(ϕ, I) b = coefficient(term) bI = if b isa MeshField b[I] elseif b isa Function - x = mesh(ϕ)[I] + x = getnode(mesh(ϕ), I) b(x, t) else error("curvature field type $b not supported") @@ -152,7 +152,7 @@ function _compute_cfl(term::CurvatureTerm, ϕ, I, t) bI = if b isa MeshField b[I] elseif b isa Function - x = mesh(ϕ)[I] + x = getnode(mesh(ϕ), I) b(x, t) else error("curvature field type $b not supported") @@ -188,13 +188,13 @@ end Base.show(io::IO, _::NormalMotionTerm) = print(io, "v|∇ϕ|") -function _compute_term(term::NormalMotionTerm, ϕ::MeshField, I, t) +function _compute_term(term::NormalMotionTerm, ϕ::AbstractMeshField, I, t) N = ndims(ϕ) u = speed(term) v = if u isa MeshField u[I] elseif u isa Function - x = mesh(ϕ)[I] + x = getnode(mesh(ϕ), I) u(x, t) else error("velocity field type $u not supported") @@ -217,7 +217,7 @@ function _compute_cfl(term::NormalMotionTerm, ϕ, I, t) v = if u isa MeshField u[I] elseif u isa Function - x = mesh(ϕ)[I] + x = getnode(mesh(ϕ), I) u(x, t) else error("velocity field type $u not supported") @@ -260,10 +260,10 @@ struct EikonalReinitializationTerm{T} <: LevelSetTerm S₀::T end -function EikonalReinitializationTerm(ϕ₀::LevelSet) +function EikonalReinitializationTerm(ϕ₀::MeshField) Δx = minimum(meshsize(ϕ₀)) # equation 7.5 of Osher and Fedkiw - S₀ = map(CartesianIndices(mesh(ϕ₀))) do I + S₀ = map(nodeindices(mesh(ϕ₀))) do I return ϕ₀[I] / sqrt(ϕ₀[I]^2 + Δx^2) end return EikonalReinitializationTerm(S₀) diff --git a/src/logging.jl b/src/logging.jl index aaa55b4..89e1fb9 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -130,13 +130,10 @@ end """ _level_set_extrema(src) -> (ϕ_min, ϕ_max) -Compute the extrema of the level-set values. Handles both `Array` and `Dict` storage. +Compute the extrema of the level-set values. """ -function _level_set_extrema(src) - v = values(src) - vals_iter = v isa AbstractDict ? Base.values(v) : v - return extrema(vals_iter) -end +_level_set_extrema(src::MeshField) = extrema(values(src)) +_level_set_extrema(src::NarrowBandMeshField) = extrema(Base.values(values(src))) """ _push_record!(log, tc, t_step, reinit_time, did_reinit, update_times, compute_times, ϕ_min, ϕ_max) diff --git a/src/meshes.jl b/src/meshes.jl index eaeded6..81679d4 100644 --- a/src/meshes.jl +++ b/src/meshes.jl @@ -17,6 +17,10 @@ end Create a uniform cartesian grid with lower corner `lc`, upper corner `hc` and `n` nodes in each direction. +Node coordinates are accessed via [`getnode`](@ref); cell geometry via [`getcell`](@ref). +Use [`nodeindices`](@ref) to iterate over all node indices and [`cellindices`](@ref) to +iterate over all cell indices. + # Examples ```jldoctest; output = true @@ -54,6 +58,9 @@ grid1d(g::CartesianGrid{N}) where {N} = ntuple(i -> grid1d(g, i), N) grid1d(g::CartesianGrid, dim) = LinRange(g.lc[dim], g.hc[dim], g.n[dim]) Base.ndims(g::CartesianGrid{N}) where {N} = N +Base.size(g::CartesianGrid) = g.n +Base.length(g::CartesianGrid) = prod(size(g)) +Base.eachindex(g::CartesianGrid) = nodeindices(g) xgrid(g::CartesianGrid) = grid1d(g, 1) ygrid(g::CartesianGrid) = grid1d(g, 2) @@ -68,29 +75,26 @@ If `dim` is not provided, return a `SVector` of spacings for all dimensions. meshsize(g::CartesianGrid) = (g.hc .- g.lc) ./ (g.n .- 1) meshsize(g::CartesianGrid, dim) = (g.hc[dim] - g.lc[dim]) / (g.n[dim] - 1) -Base.size(g::CartesianGrid) = g.n -Base.length(g::CartesianGrid) = prod(size(g)) - -function Base.getindex(g::CartesianGrid{N}, I::CartesianIndex{N}) where {N} - I ∈ CartesianIndices(g) || throw(ArgumentError("index $I is out of bounds")) - return _getindex(g, I) +# Internal unchecked coordinate lookup. Used by boundary-condition code to evaluate +# coordinates at ghost indices that lie outside the grid. +function _getindex(g::CartesianGrid{N, T}, I::CartesianIndex{N}) where {N, T} + h = meshsize(g) + return g.lc .+ (SVector{N, T}(I.I) .- 1) .* h end +_getindex(g::CartesianGrid, I::Int...) = _getindex(g, CartesianIndex(I...)) -Base.getindex(g::CartesianGrid, I::Int...) = g[CartesianIndex(I...)] - -Base.eltype(g::CartesianGrid) = typeof(g.lc) +""" + getnode(g::CartesianGrid, I::CartesianIndex) + getnode(g::CartesianGrid, i1, i2, ...) -function _getindex(g::CartesianGrid, I::CartesianIndex) - N = ndims(g) - @assert N == length(I) - return ntuple(N) do dim - return g.lc[dim] + (I[dim] - 1) / (g.n[dim] - 1) * (g.hc[dim] - g.lc[dim]) - end |> SVector +Return the coordinates of the node at multi-index `I` as an `SVector`. +`I` must be a valid node index (i.e. `I ∈ nodeindices(g)`). +""" +function getnode(g::CartesianGrid{N}, I::CartesianIndex{N}) where {N} + I ∈ nodeindices(g) || throw(ArgumentError("$I is not a valid node index for this grid")) + return _getindex(g, I) end -_getindex(g::CartesianGrid, I::Int...) = _getindex(g, CartesianIndex(I...)) - -Base.CartesianIndices(g::CartesianGrid) = CartesianIndices(size(g)) -Base.eachindex(g::CartesianGrid) = CartesianIndices(g) +getnode(g::CartesianGrid, I::Int...) = getnode(g, CartesianIndex(I...)) """ nodeindices(g::CartesianGrid) @@ -98,7 +102,7 @@ Base.eachindex(g::CartesianGrid) = CartesianIndices(g) Return a `CartesianIndices` ranging over all node indices of `g`. Nodes are indexed `1:n[d]` in each dimension `d`. """ -nodeindices(g::CartesianGrid) = CartesianIndices(g) +nodeindices(g::CartesianGrid) = CartesianIndices(size(g)) """ cellindices(g::CartesianGrid) @@ -120,6 +124,12 @@ struct CartesianCell{N, T} hc::SVector{N, T} end +# Internal unchecked cell lookup. Used when extrapolation or boundary handling is required. +function _getcell(g::CartesianGrid{N}, I::CartesianIndex{N}) where {N} + lc = _getindex(g, I) + return CartesianCell(lc, lc .+ meshsize(g)) +end + """ getcell(g::CartesianGrid, I::CartesianIndex) @@ -127,32 +137,11 @@ Return the `CartesianCell` with lower corner at node `I` and upper corner at nod `I` must be a valid cell index, i.e. `I ∈ cellindices(g)`. """ function getcell(g::CartesianGrid{N}, I::CartesianIndex{N}) where {N} - lc = g[I] - return CartesianCell(lc, lc .+ meshsize(g)) + I ∈ cellindices(g) || throw(ArgumentError("$I is not a valid cell index for this grid")) + return _getcell(g, I) end - -# iterate over all nodes -function Base.iterate(g::CartesianGrid) - i = first(CartesianIndices(g)) - return g[i], i -end - -function Base.iterate(g::CartesianGrid, state) - idxs = CartesianIndices(g) - next = iterate(idxs, state) - if next === nothing - return nothing - else - i, state = next - return g[i], state - end -end - -# Base.IteratorSize(::Type{CartesianGrid{N}}) where {N} = Base.HasShape{N}() -Base.IteratorSize(::CartesianGrid{N}) where {N} = Base.HasShape{N}() - -# --- Display --- +# Display methods """ _superscript(n::Int) -> String diff --git a/src/meshfield.jl b/src/meshfield.jl index 6090bbf..66f3251 100644 --- a/src/meshfield.jl +++ b/src/meshfield.jl @@ -1,65 +1,83 @@ """ - abstract type AbstractDomain end + abstract type AbstractMeshField -Abstract type for the domain of a [`MeshField`](@ref). +Abstract type for node-centered fields on a mesh. Concrete subtypes: +- [`MeshField`](@ref): full-domain field backed by a dense `Array`. +- [`NarrowBandMeshField`](@ref): narrow-band field backed by a sparse `Dict`. """ -abstract type AbstractDomain end +abstract type AbstractMeshField end """ - struct FullDomain <: AbstractDomain + struct MeshField{V,M,B,I} -Represents a field defined on the entire mesh. -""" -struct FullDomain <: AbstractDomain end - - -""" - struct MeshField{V,M,B,D} +A node-centered field defined on the entire mesh, described by its discrete values +at each node. -A field described by its discrete values on a mesh. - -- `vals`: the discrete values of the field. +- `vals`: dense array of values at each node. - `mesh`: the underlying mesh (e.g. [`CartesianGrid`](@ref)). - `bcs`: boundary conditions, used for indexing outside the mesh bounds. -- `domain`: the domain on which the field is defined (e.g. [`FullDomain`](@ref)). +- `itp_data`: optional `InterpolationData` for piecewise polynomial interpolation + (`nothing` if the field was constructed without an interpolation order). + +`Base.getindex(ϕ, I)` returns the value of the field at node index `I`. If `I` lies +outside the mesh bounds, `bcs` are applied to determine the value. -`Base.getindex` of an `MeshField` is overloaded to handle indices that lie outside the -`CartesianIndices` of its `MeshField` by using `bcs`. +Use [`nodeindices`](@ref) and [`cellindices`](@ref) to iterate over active node and cell +indices respectively. """ -struct MeshField{V, M, B, D <: AbstractDomain} +struct MeshField{V, M, B, I} <: AbstractMeshField vals::V mesh::M bcs::B - domain::D + itp_data::I +end + +# getters (on AbstractMeshField where applicable) +mesh(ϕ::AbstractMeshField) = ϕ.mesh +Base.values(ϕ::AbstractMeshField) = ϕ.vals +has_boundary_conditions(ϕ::AbstractMeshField) = !isnothing(ϕ.bcs) +boundary_conditions(ϕ::AbstractMeshField) = ϕ.bcs +interp_data(ϕ::AbstractMeshField) = ϕ.itp_data +has_interpolation(ϕ::AbstractMeshField) = !isnothing(ϕ.itp_data) +Base.valtype(ϕ::AbstractMeshField) = valtype(values(ϕ)) + +function check_real_valued(mf::AbstractMeshField) + valtype(mf) <: Real || throw(ArgumentError("Expected a real-valued MeshField, but got valtype $(valtype(mf))")) + return nothing end -# getters -mesh(ϕ::MeshField) = ϕ.mesh -Base.values(ϕ::MeshField) = ϕ.vals -domain(ϕ::MeshField) = ϕ.domain -has_boundary_conditions(ϕ::MeshField) = !isnothing(ϕ.bcs) -boundary_conditions(ϕ::MeshField) = ϕ.bcs +function _check_bc(ϕ::AbstractMeshField) + has_boundary_conditions(ϕ) || throw( + ArgumentError( + "this operation requires boundary conditions. Pass `bc=...` when constructing the MeshField." + ) + ) + return nothing +end -meshsize(ϕ::MeshField, args...) = meshsize(mesh(ϕ), args...) +meshsize(ϕ::AbstractMeshField, args...) = meshsize(mesh(ϕ), args...) """ - add_boundary_conditions(ϕ::MeshField, bc) + _add_boundary_conditions(ϕ::MeshField, bc) Return a new `MeshField` with `bc` as boundary conditions. All of the underlying data is -aliased (shared) with the original `MeshField`. +aliased (shared) with the original `MeshField`. The `itp_data` is copied with its cache +invalidated so the new field owns independent interpolation state. """ -function add_boundary_conditions(ϕ::MeshField, bc) +function _add_boundary_conditions(ϕ::MeshField, bc) N = ndims(ϕ) - return MeshField(values(ϕ), mesh(ϕ), _normalize_bc(bc, N), domain(ϕ)) + itp = interp_data(ϕ) + new_itp = isnothing(itp) ? nothing : copy(itp) + return MeshField(values(ϕ), mesh(ϕ), _normalize_bc(bc, N), new_itp) end """ - update_bcs!(ϕ::MeshField, t) + update_bcs!(ϕ::AbstractMeshField, t) Update the current time in all [`DirichletBC`](@ref) boundary conditions of `ϕ`. Called automatically by the time-stepper at each stage. """ -function update_bcs!(ϕ::MeshField, t) +function update_bcs!(ϕ::AbstractMeshField, t) has_boundary_conditions(ϕ) || (return ϕ) for bc_pair in boundary_conditions(ϕ) update_bc!(bc_pair[1], t) @@ -69,17 +87,36 @@ function update_bcs!(ϕ::MeshField, t) end """ - MeshField(vals, mesh, bcs) + MeshField(vals, grid::AbstractMesh; bc=nothing, interp_order=nothing) + +Construct a node-centered `MeshField` with explicit values on a `grid`. -Construct a `MeshField` with explicit values, mesh, and boundary conditions. -Defaults to `FullDomain`. +- `bc`: boundary conditions (normalized via [`_normalize_bc`](@ref)). +- `interp_order`: optional polynomial order for piecewise interpolation. + +If `interp_order` is provided but `bc` is `nothing`, `bc` defaults to +`ExtrapolationBC{2}` to ensure valid stencils at the boundary. """ -MeshField(vals, mesh, bcs) = MeshField(vals, mesh, bcs, FullDomain()) +function MeshField(vals, grid::AbstractMesh; bc = nothing, interp_order = nothing) + N = ndims(grid) + T = valtype(vals) + bcs = isnothing(bc) ? nothing : _normalize_bc(bc, N) + + # If interpolation is requested but no BCs are provided, default to 2nd-order + # extrapolation to allow stencil access at boundaries. + if !isnothing(interp_order) && isnothing(bcs) + bcs = ntuple(_ -> (ExtrapolationBC{2}(), ExtrapolationBC{2}()), N) + end + + itp = isnothing(interp_order) ? nothing : InterpolationData(N, interp_order, T) + return MeshField(vals, grid, bcs, itp) +end """ - MeshField(f::Function, m) + MeshField(f::Function, grid; kwargs...) -Create a `MeshField` by evaluating a function `f` on a mesh `m`. +Create a node-centered `MeshField` by evaluating `f` at each node of `grid`. +All keyword arguments are passed to the `vals`-based constructor. # Examples @@ -98,55 +135,33 @@ MeshField on CartesianGrid in ℝ² ├─ eltype: Float64 └─ values: min = -0.25, max = 1.75 ``` - -```jldoctest; output = true -using LevelSetMethods -grid = CartesianGrid((-1, -1), (1, 1), (5, 5)) -# vector-valued field -MeshField(x -> (x[1], x[2]), grid) - -# output - -MeshField on CartesianGrid in ℝ² - ├─ domain: [-1.0, 1.0] × [-1.0, 1.0] - ├─ nodes: 5 × 5 - ├─ spacing: h = (0.5, 0.5) - └─ eltype: Tuple{Float64, Float64} -``` """ -function MeshField(f::Function, m, bc = nothing) - bc_ = isnothing(bc) ? nothing : _normalize_bc(bc, ndims(m)) - vals = map(f, m) - return MeshField(vals, m, bc_, FullDomain()) +function MeshField(f::Function, grid::AbstractMesh; kwargs...) + vals = map(I -> f(getnode(grid, I)), nodeindices(grid)) + return MeshField(vals, grid; kwargs...) end # geometric dimension -Base.ndims(f::MeshField) = ndims(mesh(f)) +Base.ndims(f::AbstractMeshField) = ndims(mesh(f)) # Base.length -Base.length(ϕ::MeshField) = length(eachindex(ϕ)) +Base.length(ϕ::AbstractMeshField) = length(eachindex(ϕ)) # overload base methods for convenience -function Base.getindex(ϕ::MeshField, I::CartesianIndex{N}) where {N} +function Base.getindex(ϕ::AbstractMeshField, I::CartesianIndex{N}) where {N} if has_boundary_conditions(ϕ) return _getindexbc(ϕ, I, N) else return _base_lookup(ϕ, I) end end -function Base.getindex(ϕ::MeshField, I...) +function Base.getindex(ϕ::AbstractMeshField, I...) return ϕ[CartesianIndex(I...)] end -# Recursive getindex with boundary conditions, processing one dimension per -# call (dim = N down to 1). If I[dim] is in-bounds, recurse; if out-of-bounds, -# apply the BC for that dimension then recurse to dim-1. Base case dim=0 does -# the raw array lookup. Corner ghost points (out-of-bounds in multiple -# dimensions) are handled correctly because each dimension's BC is applied in -# turn. _base_lookup(ϕ::MeshField, I) = getindex(values(ϕ), I) -function _getindexbc(ϕ::MeshField, I, dim) +function _getindexbc(ϕ::AbstractMeshField, I, dim) dim == 0 && return _base_lookup(ϕ, I) bcs = boundary_conditions(ϕ)[dim] ax = axes(ϕ)[dim] @@ -203,50 +218,60 @@ function _apply_extrapolation_bc(ϕ, I::CartesianIndex{N}, ::ExtrapolationBC{P}, end -Base.setindex!(ϕ::MeshField, vals, I...) = setindex!(values(ϕ), vals, I...) +function Base.setindex!(ϕ::AbstractMeshField, val, I...) + _invalidate_itp!(ϕ) + setindex!(values(ϕ), val, I...) + return ϕ +end + +_invalidate_itp!(ϕ::MeshField{V, M, B, Nothing}) where {V, M, B} = nothing +function _invalidate_itp!(ϕ::MeshField{V, M, B, <:InterpolationData{N}}) where {V, M, B, N} + ϕ.itp_data.Ic = CartesianIndex(ntuple(_ -> 0, Val(N))) + return nothing +end Base.eltype(ϕ::MeshField) = eltype(values(ϕ)) -function Base.axes(ϕ::MeshField) +function Base.axes(ϕ::AbstractMeshField) sz = size(mesh(ϕ)) return ntuple(d -> Base.OneTo(sz[d]), Val(ndims(ϕ))) end -Base.eachindex(ϕ::MeshField) = _eachindex(domain(ϕ), ϕ) -_eachindex(::FullDomain, ϕ) = eachindex(mesh(ϕ)) +Base.eachindex(ϕ::MeshField) = nodeindices(mesh(ϕ)) + +""" + nodeindices(ϕ::MeshField) +Return the active node indices of `ϕ`. For a `MeshField` this is `nodeindices(mesh(ϕ))`. """ - Base.copy!(dest::MeshField, src::MeshField) +nodeindices(ϕ::MeshField) = nodeindices(mesh(ϕ)) + +""" + cellindices(ϕ::MeshField) + +Return the active cell indices of `ϕ`. For a `MeshField` this is `cellindices(mesh(ϕ))`. +""" +cellindices(ϕ::MeshField) = cellindices(mesh(ϕ)) + +""" + Base.copy!(dest::AbstractMeshField, src::AbstractMeshField) Copy the values from `src` to `dest`. The meshes, boundary conditions, and domains of the -`dest` fields are not modified. +`dest` fields are not modified. The interpolation cache of `dest` is invalidated. """ function Base.copy!(dest::MeshField, src::MeshField) + _invalidate_itp!(dest) copy!(values(dest), values(src)) return dest end -""" - _show_fields(io, ϕ::MeshField; prefix=" ") - -Print the fields of `ϕ` as indented tree lines: grid info (via `_show_fields` for -`CartesianGrid`), boundary conditions, narrow-band info (if applicable), element type, -and value range (for real-valued fields). -""" function _show_fields(io::IO, ϕ::MeshField{<:Any, <:CartesianGrid}; prefix = " ") _show_fields(io, mesh(ϕ); prefix, last = false) if has_boundary_conditions(ϕ) println(io, "$(prefix)├─ bc: $(_bc_str(boundary_conditions(ϕ)))") end - dom = domain(ϕ) - if dom isa NarrowBandDomain - hw = dom.halfwidth - nlayers = round(Int, hw / minimum(meshsize(mesh(ϕ)))) - println(io, "$(prefix)├─ active: $(length(active_indices(ϕ))) nodes ($nlayers layers, halfwidth = $(round(hw; sigdigits = 4)))") - end return if eltype(ϕ) <: Real v = values(ϕ) - vals_iter = v isa AbstractDict ? Base.values(v) : v - vmin, vmax = extrema(vals_iter) + vmin, vmax = extrema(v) println(io, "$(prefix)├─ eltype: $(eltype(ϕ))") print(io, "$(prefix)└─ values: min = $(round(vmin; sigdigits = 4)), max = $(round(vmax; sigdigits = 4))") else @@ -258,3 +283,207 @@ function Base.show(io::IO, ::MIME"text/plain", ϕ::MeshField{<:Any, <:CartesianG println(io, "MeshField on CartesianGrid in ℝ$(_superscript(ndims(ϕ)))") return _show_fields(io, ϕ) end + +# ---- NarrowBandMeshField -------------------------------------------------------- +# Struct and basic methods here; constructors that call NewtonReinitializer are in narrowband.jl. + +""" + struct NarrowBandMeshField{V,M,B,T,I} + +A node-centered field defined on a narrow band around the interface, described by its +discrete values at each *active* (in-band) node. + +- `vals`: sparse `Dict{CartesianIndex{N},T}` mapping active node indices to values. +- `mesh`: the underlying mesh (e.g. [`CartesianGrid`](@ref)). +- `bcs`: boundary conditions, used for indexing outside the mesh bounds. +- `halfwidth`: half-width of the narrow band in physical units (typically a few grid spacings). +- `itp_data`: optional `InterpolationData` for piecewise polynomial interpolation. + +Use [`nodeindices`](@ref) and [`cellindices`](@ref) to iterate over active node and cell +indices respectively. +""" +struct NarrowBandMeshField{V, M, B, T, I} <: AbstractMeshField + vals::V + mesh::M + bcs::B + halfwidth::T + itp_data::I +end + +halfwidth(nb::NarrowBandMeshField) = nb.halfwidth +Base.eltype(nb::NarrowBandMeshField) = valtype(nb) + +""" + nodeindices(nb::NarrowBandMeshField) + +Return the set of active (in-band) node indices of `nb`. +""" +nodeindices(nb::NarrowBandMeshField) = keys(values(nb)) + +""" + cellindices(nb::NarrowBandMeshField) + +Return the set of active cell indices of `nb`: cells whose corners are all in-band nodes. +""" +function cellindices(nb::NarrowBandMeshField{<:Any, <:AbstractMesh{N}}) where {N} + grid = mesh(nb) + cell_axes = cellindices(grid) + active_nodes = nodeindices(nb) + cells = Set{CartesianIndex{N}}() + checked = Set{CartesianIndex{N}}() + offsets = Iterators.product(ntuple(_ -> 0:1, Val(N))...) + for I in active_nodes + for offset in offsets + J = I - CartesianIndex(offset) + J in cell_axes || continue + J in checked && continue + push!(checked, J) + if all(J + CartesianIndex(off) in active_nodes for off in offsets) + push!(cells, J) + end + end + end + return cells +end + +Base.eachindex(nb::NarrowBandMeshField) = nodeindices(nb) + +_invalidate_itp!(nb::NarrowBandMeshField{V, M, B, T, Nothing}) where {V, M, B, T} = nothing +function _invalidate_itp!(nb::NarrowBandMeshField{V, M, B, T, <:InterpolationData{N}}) where {V, M, B, T, N} + nb.itp_data.Ic = CartesianIndex(ntuple(_ -> 0, Val(N))) + return nothing +end + +function Base.copy!(dest::NarrowBandMeshField, src::NarrowBandMeshField) + _invalidate_itp!(dest) + empty!(values(dest)) + merge!(values(dest), values(src)) + return dest +end + +""" + _add_boundary_conditions(nb::NarrowBandMeshField, bc) + +Return a new `NarrowBandMeshField` with `bc` as boundary conditions, preserving +`halfwidth`. All underlying data is aliased with the original. +""" +function _add_boundary_conditions(nb::NarrowBandMeshField, bc) + N = ndims(nb) + itp = interp_data(nb) + new_itp = isnothing(itp) ? nothing : copy(itp) + return NarrowBandMeshField(values(nb), mesh(nb), _normalize_bc(bc, N), nb.halfwidth, new_itp) +end + +""" + _base_lookup(nb::NarrowBandMeshField, I) -> value + +Lookup the value at in-grid index `I`. Tries the dict first; if not stored, falls back +to linear extrapolation from nearby band nodes. +""" +function _base_lookup(nb::NarrowBandMeshField{<:Any, <:AbstractMesh{N}}, I) where {N} + val = get(values(nb), I, nothing) + val !== nothing && return val + val = _extrapolate_nb_rec(nb, I, N) + val !== nothing && return val + error("extrapolation failed at index $I: no resolvable path to stored values") +end + +function _extrapolate_nb_rec(nb::NarrowBandMeshField{<:Any, <:AbstractMesh{N}}, I::CartesianIndex{N}, max_dim) where {N} + haskey(values(nb), I) && return values(nb)[I] + grid_axes = axes(nb) + P = 1 + for dim in 1:max_dim + for k in 1:length(grid_axes[dim]) + for side in (-1, 1) + anchor = I[dim] + side * k + anchor in grid_axes[dim] || continue + val = _lagrange_extrap_from(nb, I, dim, anchor, side, k, P) + val !== nothing && return val + end + end + end + return nothing +end + +function _lagrange_extrap_from(nb::NarrowBandMeshField{<:Any, <:AbstractMesh{N}}, I, dim, anchor, side, k, P) where {N} + T = eltype(nb) + grid_axes = axes(nb) + result = zero(float(T)) + for j in 0:P + pos = anchor + side * j + pos in grid_axes[dim] || return nothing + Ij = CartesianIndex(ntuple(s -> s == dim ? pos : I[s], Val(N))) + Vj = _extrapolate_nb_rec(nb, Ij, dim - 1) + Vj === nothing && return nothing + result += _lagrange_extrap_weight(j, k, P) * Vj + end + return result +end + +""" + _clear_buffer!(nb::NarrowBandMeshField) + +Empty the active-node dict before reuse as a write target in time integration. +""" +_clear_buffer!(::MeshField) = nothing +_clear_buffer!(nb::NarrowBandMeshField) = empty!(values(nb)) + +""" + NarrowBandMeshField(vals, grid; halfwidth, bc=nothing, interp_order=nothing) + +Construct a `NarrowBandMeshField` from a pre-built `vals` dict. `halfwidth` is required. +""" +function NarrowBandMeshField(vals, grid::AbstractMesh; halfwidth, bc = nothing, interp_order = nothing) + N = ndims(grid) + T = valtype(vals) + bcs = isnothing(bc) ? nothing : _normalize_bc(bc, N) + if !isnothing(interp_order) && isnothing(bcs) + bcs = ntuple(_ -> (ExtrapolationBC{2}(), ExtrapolationBC{2}()), N) + end + itp = isnothing(interp_order) ? nothing : InterpolationData(N, interp_order, T) + return NarrowBandMeshField(vals, grid, bcs, halfwidth, itp) +end + +""" + NarrowBandMeshField(f, grid, halfwidth; bc=nothing, interp_order=nothing) + +Construct a narrow-band field by evaluating `f` at each node of `grid` and keeping +only those where `|f(x)| < halfwidth`. No dense array is allocated. + +Pass `interp_order=k` to enable piecewise polynomial interpolation (same semantics +as [`MeshField`](@ref)). + +!!! warning + Since the `halfwidth` threshold is applied to the raw values of `f`, the resulting band + width in physical space will only match `halfwidth` if `f` is already a signed distance + function. +""" +function NarrowBandMeshField(f::Function, grid::AbstractMesh, halfwidth::Real; bc = nothing, interp_order = nothing) + T = float(eltype(grid.lc)) + γ = T(halfwidth) + vals = _nb_dict(I -> f(getnode(grid, I)), grid, γ) + return NarrowBandMeshField(vals, grid; bc = bc, halfwidth = γ, interp_order = interp_order) +end + +function _show_fields(io::IO, nb::NarrowBandMeshField{<:Any, <:CartesianGrid}; prefix = " ") + _show_fields(io, mesh(nb); prefix, last = false) + if has_boundary_conditions(nb) + println(io, "$(prefix)├─ bc: $(_bc_str(boundary_conditions(nb)))") + end + hw = halfwidth(nb) + nlayers = round(Int, hw / minimum(meshsize(mesh(nb)))) + println(io, "$(prefix)├─ active: $(length(nodeindices(nb))) nodes ($nlayers layers, halfwidth = $(round(hw; sigdigits = 4)))") + return if eltype(nb) <: Real + vals_iter = Base.values(values(nb)) + vmin, vmax = extrema(vals_iter) + println(io, "$(prefix)├─ eltype: $(eltype(nb))") + print(io, "$(prefix)└─ values: min = $(round(vmin; sigdigits = 4)), max = $(round(vmax; sigdigits = 4))") + else + print(io, "$(prefix)└─ eltype: $(eltype(nb))") + end +end + +function Base.show(io::IO, ::MIME"text/plain", nb::NarrowBandMeshField{<:Any, <:CartesianGrid}) + println(io, "NarrowBandMeshField on CartesianGrid in ℝ$(_superscript(ndims(nb)))") + return _show_fields(io, nb) +end diff --git a/src/narrowband.jl b/src/narrowband.jl index 68c9bb4..048b29a 100644 --- a/src/narrowband.jl +++ b/src/narrowband.jl @@ -1,38 +1,8 @@ -""" - struct NarrowBandDomain{T} <: AbstractDomain - -Domain for a narrow-band level set. - -- `halfwidth`: half-width of the narrow band, typically on the order of a few grid spacings. - -Active indices are the keys of the associated values dict and need not be stored separately. -""" -struct NarrowBandDomain{T} <: AbstractDomain - halfwidth::T -end - -""" - const NarrowBandLevelSet{N, T, B} - -Alias for [`MeshField`](@ref) on a `CartesianGrid{N,T}` with values stored as -a `Dict{CartesianIndex{N},T}` and a [`NarrowBandDomain{T}`](@ref). `B` is -the type of the boundary conditions. -""" -const NarrowBandLevelSet{N, T, B} = - MeshField{Dict{CartesianIndex{N}, T}, CartesianGrid{N, T}, B, NarrowBandDomain{T}} - -active_indices(nb::NarrowBandLevelSet) = keys(values(nb)) -halfwidth(nb::NarrowBandLevelSet) = domain(nb).halfwidth -Base.eachindex(nb::NarrowBandLevelSet) = active_indices(nb) -_eachindex(::NarrowBandDomain, nb) = active_indices(nb) -Base.eltype(::NarrowBandLevelSet{N, T}) where {N, T} = T - # Build the active-index dict by evaluating `f_at_idx(I)` at every grid node and -# keeping only those where `|v| < γ`. `f_at_idx` is either `I -> ϕ[I]` (LevelSet path) -# or `I -> f(grid[I])` (function path). +# keeping only those where `|v| < γ`. function _nb_dict(f_at_idx, grid::CartesianGrid{N}, γ::T) where {N, T} vals = Dict{CartesianIndex{N}, T}() - for I in CartesianIndices(grid) + for I in nodeindices(grid) v = T(f_at_idx(I)) abs(v) < γ && (vals[I] = v) end @@ -40,21 +10,23 @@ function _nb_dict(f_at_idx, grid::CartesianGrid{N}, γ::T) where {N, T} end """ - NarrowBandLevelSet(ϕ::LevelSet, halfwidth::Real; reinitialize = true) + NarrowBandMeshField(ϕ::MeshField, halfwidth::Real; reinitialize=true) -Construct a `NarrowBandLevelSet` from a full-grid `LevelSet`. Active nodes are those -where `|ϕ[I]| < halfwidth`. Boundary conditions are inherited from `ϕ`. +Construct a narrow-band field from a full-grid [`MeshField`](@ref). +Active nodes are those where `|ϕ[I]| < halfwidth`. Boundary conditions are +inherited from `ϕ`. If `reinitialize` is `true` (the default), `ϕ` is first reinitialized to a signed distance function using [`NewtonReinitializer`](@ref). """ -function NarrowBandLevelSet(ϕ::LevelSet, halfwidth::Real; reinitialize::Bool = true) +function NarrowBandMeshField(ϕ::MeshField, halfwidth::Real; reinitialize::Bool = true) bcs = boundary_conditions(ϕ) # preserve the caller's BCs (may be nothing) + itp = interp_data(ϕ) # preserve interpolation data if reinitialize ϕ = deepcopy(ϕ) # reinit needs BCs for gradient computation; add temporary ones if missing if !has_boundary_conditions(ϕ) - ϕ = add_boundary_conditions(ϕ, ExtrapolationBC(2)) + ϕ = _add_boundary_conditions(ϕ, ExtrapolationBC(2)) end reinitialize!(ϕ, NewtonReinitializer()) end @@ -62,198 +34,50 @@ function NarrowBandLevelSet(ϕ::LevelSet, halfwidth::Real; reinitialize::Bool = T = float(eltype(ϕ)) γ = T(halfwidth) vals = _nb_dict(I -> ϕ[I], grid, γ) - return MeshField(vals, grid, bcs, NarrowBandDomain(γ)) -end - -""" - NarrowBandLevelSet(ϕ::LevelSet; nlayers = 3, reinitialize = true) - -Construct a `NarrowBandLevelSet` with halfwidth automatically computed as -`nlayers * minimum(meshsize(ϕ))`. `nlayers` sets the number of cell layers -on each side of the interface included in the band. -""" -function NarrowBandLevelSet(ϕ::LevelSet; nlayers::Int = 3, reinitialize::Bool = true) - return NarrowBandLevelSet(ϕ, nlayers * minimum(meshsize(ϕ)); reinitialize) -end - -""" - NarrowBandLevelSet(f, grid::CartesianGrid, halfwidth::Real; bc = nothing) - -Construct a `NarrowBandLevelSet` by evaluating `f` at each node of `grid` and -keeping only those where `|f(x)| < halfwidth`. No dense array is allocated. - -!!! warning - Since the `halfwidth` threshold is applied to the raw values of `f`, the resulting band - width in physical space will only match `halfwidth` if `f` is already a signed distance - function. Otherwise the band width will depend on the gradient of `f` near the interface - and may not correspond to a fixed number of cell layers. -""" -function NarrowBandLevelSet(f, grid::CartesianGrid, halfwidth::Real; bc = nothing) - T = float(eltype(eltype(grid))) - γ = T(halfwidth) - vals = _nb_dict(I -> f(grid[I]), grid, γ) - return MeshField(vals, grid, bc, NarrowBandDomain(γ)) + return NarrowBandMeshField(vals, grid, bcs, γ, isnothing(itp) ? nothing : copy(itp)) end """ - NarrowBandLevelSet(f, grid::CartesianGrid; nlayers = 8, bc = nothing) - -Construct a `NarrowBandLevelSet` with halfwidth automatically computed as -`nlayers * minimum(meshsize(grid))`. - -!!! warning - The `nlayers` interpretation is only correct if `f` is already a signed distance - function. Otherwise the band width in cell layers will not match `nlayers`. -""" -function NarrowBandLevelSet(f, grid::CartesianGrid; nlayers::Int = 8, bc = nothing) - return NarrowBandLevelSet(f, grid, nlayers * minimum(meshsize(grid)); bc) -end + NarrowBandMeshField(ϕ::MeshField; nlayers=3, reinitialize=true) +Construct a narrow-band field with halfwidth automatically computed as +`nlayers * minimum(meshsize(ϕ))`. """ - _base_lookup(nb::NarrowBandLevelSet, I) -> value - -Entry point for value lookup on a `NarrowBandLevelSet` at index `I`, which is -assumed to be inside the grid (out-of-grid indices are handled by `_getindexbc` -before reaching this function). - -Tries the dict first; if `I` is not stored (i.e. it is inside the grid but -outside the narrow band), falls back to [`_extrapolate_nb_rec`](@ref) to -approximate the value from nearby band nodes. Throws an error if no path to -stored values can be found. -""" -function _base_lookup(nb::NarrowBandLevelSet{N}, I) where {N} - val = get(values(nb), I, nothing) - val !== nothing && return val - val = _extrapolate_nb_rec(nb, I, N) - val !== nothing && return val - error("extrapolation failed at index $I: no resolvable path to stored values") +function NarrowBandMeshField(ϕ::MeshField; nlayers::Int = 3, reinitialize::Bool = true) + return NarrowBandMeshField(ϕ, nlayers * minimum(meshsize(ϕ)); reinitialize) end -""" - _extrapolate_nb_rec(nb::NarrowBandLevelSet, I, max_dim) -> value or nothing +# ---- Reinitialization ------------------------------------------------------------ -Approximate the value at an in-grid index `I` that is not stored in the band -dict, by bilinear/trilinear (degree-1) Lagrange extrapolation from nearby band values. +reinitialize!(nb::NarrowBandMeshField, ::Nothing, _) = nb -The algorithm processes dimensions `1` through `max_dim` in order. For each -dimension, it searches outward from `I` (nearest first, both sides) for an -anchor point in the dict and uses the two consecutive nodes starting there -(anchor and anchor+step) as a linear stencil. - -Stencil values are resolved by calling `_extrapolate_nb_rec` recursively with -`max_dim = dim - 1`, so each stencil point can itself be extrapolated using -lower dimensions. This produces a tensor-product extrapolation that handles -indices outside the band in multiple dimensions simultaneously. - -Returns `nothing` if no dimension yields a valid stencil. - -# Why degree 1 and no higher? -Outside the band, `|ϕ| ≥ halfwidth`. The only property we need from extrapolated -ghost values is sign correctness (no spurious zeros). Linear extrapolation from a -well-conditioned SDF preserves sign as long as `halfwidth > Δx`, which any -reasonable band satisfies. Quadratic or higher extrapolation introduces polynomial -oscillations that can flip the sign, creating false interface crossings that corrupt -`NewtonSDF` and cause the band to migrate or collapse. -""" -function _extrapolate_nb_rec(nb::NarrowBandLevelSet{N, T}, I::CartesianIndex{N}, max_dim) where {N, T} - haskey(values(nb), I) && return values(nb)[I] - grid_axes = axes(nb) - # Degree-1 (linear) stencil — see docstring for why we never go higher. - P = 1 - for dim in 1:max_dim - for k in 1:length(grid_axes[dim]) - for side in (-1, 1) - anchor = I[dim] + side * k - anchor in grid_axes[dim] || continue - val = _lagrange_extrap_from(nb, I, dim, anchor, side, k, P) - val !== nothing && return val - end - end - end - return nothing -end - -""" - _lagrange_extrap_from(nb, I, dim, anchor, side, k, P) -> value or nothing - -Attempt to evaluate a degree-`P` Lagrange extrapolant at `I[dim]` using `P+1` -consecutive stencil nodes along dimension `dim`: - - anchor, anchor + side, anchor + 2*side, …, anchor + P*side - -The stencil extends from `anchor` deeper into the band (away from `I`). The -target `I[dim]` is at distance `k` from `anchor` in the opposite direction, -corresponding to local coordinate `ξ = -k` relative to stencil nodes at -`ξ = 0, 1, …, P`. This matches the convention of `_lagrange_extrap_weight(j, k, P)`. - -Each stencil value is resolved via `_extrapolate_nb_rec(nb, Ij, dim - 1)`, -using only dimensions lower than `dim`. Returns `nothing` if any stencil point -falls outside the grid or cannot be resolved. -""" -function _lagrange_extrap_from(nb::NarrowBandLevelSet{N, T}, I, dim, anchor, side, k, P) where {N, T} - grid_axes = axes(nb) - result = zero(float(T)) - for j in 0:P - pos = anchor + side * j - pos in grid_axes[dim] || return nothing - Ij = CartesianIndex(ntuple(s -> s == dim ? pos : I[s], Val(N))) - Vj = _extrapolate_nb_rec(nb, Ij, dim - 1) - Vj === nothing && return nothing - result += _lagrange_extrap_weight(j, k, P) * Vj - end - return result -end - -function active_cells(nb::NarrowBandLevelSet) - grid = mesh(nb) - cell_axes = cellindices(grid) - active_nodes = active_indices(nb) - active_cells = Set{CartesianIndex}() - - for I in cell_axes - N = ndims(grid) - all_corners_active = true - for offset in Iterators.product(ntuple(_ -> 0:1, Val(N))...) - corner_idx = I + CartesianIndex(offset) - if !(corner_idx in active_nodes) - all_corners_active = false - break - end - end - all_corners_active && push!(active_cells, I) - end - return active_cells -end - -reinitialize!(nb::NarrowBandLevelSet, ::Nothing, _) = nb - -function reinitialize!(nb::NarrowBandLevelSet, r::NewtonReinitializer) +function reinitialize!(nb::NarrowBandMeshField, r::NewtonReinitializer) sdf = NewtonSDF(nb; order = r.order, upsample = r.upsample, maxiters = r.maxiters, xtol = r.xtol, ftol = r.ftol) rebuild_band!(nb, sdf) return nb end -function reinitialize!(nb::NarrowBandLevelSet, r::NewtonReinitializer, nsteps::Int) +function reinitialize!(nb::NarrowBandMeshField, r::NewtonReinitializer, nsteps::Int) mod(nsteps, r.reinit_freq) == 0 || return nb return reinitialize!(nb, r) end +# ---- Band rebuild ---------------------------------------------------------------- + """ - rebuild_band!(nb::NarrowBandLevelSet, sdf) + rebuild_band!(nb::NarrowBandMeshField, sdf) Rebuild the active node set from the signed distance function `sdf` using a breadth-first search seeded from all previously active nodes. The BFS expands axis-aligned neighbors and adds a node whenever `|sdf(x)| < halfwidth`, stopping a branch when a node falls outside the band. """ -function rebuild_band!(nb::NarrowBandLevelSet{N, T}, sdf) where {N, T} +function rebuild_band!(nb::NarrowBandMeshField{<:Any, <:AbstractMesh{N}}, sdf) where {N} grid = mesh(nb) γ = halfwidth(nb) grid_axes = axes(nb) vals = values(nb) - # Use a vector queue for BFS, keeping track of the head explicitly. Duplicate indices in - # a set to have O(1) membership check. queue = collect(keys(vals)) queue_set = Set{CartesianIndex{N}}(queue) empty!(vals) @@ -262,8 +86,8 @@ function rebuild_band!(nb::NarrowBandLevelSet{N, T}, sdf) where {N, T} while head <= length(queue) I = queue[head] head += 1 - v = sdf(grid[I]) - abs(v) >= γ && continue # outside band — don't expand further + v = sdf(getnode(grid, I)) + abs(v) >= γ && continue vals[I] = v for d in 1:N, s in (-1, 1) J = _increment_index(I, d, s) @@ -275,14 +99,3 @@ function rebuild_band!(nb::NarrowBandLevelSet{N, T}, sdf) where {N, T} end return nb end - -""" - _clear_buffer!(ϕ::MeshField) - -Clear the active entries of a buffer before it is used as a write target in a -time-integration step. For a [`NarrowBandLevelSet`](@ref) this empties the -values dict so that stale entries from a previous band do not survive the -src/dst swap. For a full-domain field this is a no-op. -""" -_clear_buffer!(::MeshField) = nothing -_clear_buffer!(nb::NarrowBandLevelSet) = empty!(values(nb)) diff --git a/src/reinitializer.jl b/src/reinitializer.jl index 249fa8e..67c19e6 100644 --- a/src/reinitializer.jl +++ b/src/reinitializer.jl @@ -9,7 +9,7 @@ from a new level set. abstract type AbstractSignedDistanceFunction <: Function end """ - update!(sdf::AbstractSignedDistanceFunction, ϕ::LevelSet) + update!(sdf::AbstractSignedDistanceFunction, ϕ::MeshField) Rebuild `sdf` in place from the new level set `ϕ`. """ @@ -30,8 +30,8 @@ optional second argument `sdf(x, s)` can be used to supply the sign directly (e. internal mutable buffers. For multi-threaded evaluation, use `deepcopy(sdf)` for each thread. """ -mutable struct NewtonSDF{I, Tr, P} <: AbstractSignedDistanceFunction - itp::I +mutable struct NewtonSDF{Φ, Tr, P} <: AbstractSignedDistanceFunction + meshfield::Φ tree::Tr pts::P upsample::Int @@ -57,45 +57,55 @@ function Base.show(io::IO, ::MIME"text/plain", sdf::NewtonSDF) end """ - NewtonSDF(itp; upsample=8, maxiters=20, xtol=1e-8, ftol=1e-8) + NewtonSDF(ϕ; order=3, upsample=2, maxiters=10, xtol=1e-8, ftol=1e-8) -Construct a [`NewtonSDF`](@ref) from a `PiecewisePolynomialInterpolant`. +Construct a [`NewtonSDF`](@ref) from a level set `ϕ`. If `ϕ` already carries +`InterpolationData` (i.e. was constructed with an `interp_order` keyword), it is used +directly; otherwise a deep copy is made with `ExtrapolationBC{2}` added and +`InterpolationData` of the given `order` attached. -The interface is sampled by projecting uniformly-spaced points in each cell onto the zero -level set. A KDTree is built from these samples for fast nearest-neighbor queries. +The interface is sampled by projecting uniformly-spaced points in each cell onto the +zero level set, building a KDTree for fast nearest-neighbour queries. + +Works for both [`MeshField`](@ref) and narrow-band fields. # Keyword arguments -- `upsample`: sampling density per cell side. Larger values means a denser sampling of the - interface is used to build the `KDTree`, which in turn usually means a better initial - guess for the Newton solver. +- `order`: polynomial interpolation order (only used when `ϕ` has no interpolation) +- `upsample`: sampling density per cell side - `maxiters`: maximum Newton iterations -- `xtol`: tolerance on iterate updates for convergence of the Newton solver -- `ftol`: tolerance on the function value for convergence of the Newton solver +- `xtol`: tolerance on the KKT residual +- `ftol`: tolerance on the function value """ function NewtonSDF( - itp::PiecewisePolynomialInterpolant; + ϕ; + order = 3, upsample = 2, maxiters = 10, xtol = 1.0e-8, - ftol = 1.0e-8 + ftol = 1.0e-8, ) - grid = mesh(itp.ϕ) - pts = _sample_interface(grid, itp, active_cells(itp.ϕ), upsample, maxiters, ftol) + if has_interpolation(ϕ) + mf = ϕ + else + ϕ_copy = deepcopy(ϕ) + if !has_boundary_conditions(ϕ_copy) + N = ndims(mesh(ϕ_copy)) + ϕ_copy = _add_boundary_conditions(ϕ_copy, ntuple(_ -> (ExtrapolationBC{2}(), ExtrapolationBC{2}()), N)) + end + N, T = ndims(mesh(ϕ_copy)), eltype(ϕ_copy) + if ϕ_copy isa NarrowBandMeshField + mf = NarrowBandMeshField( + values(ϕ_copy), mesh(ϕ_copy), boundary_conditions(ϕ_copy), + halfwidth(ϕ_copy), InterpolationData(N, order, T), + ) + else + mf = MeshField(values(ϕ_copy), mesh(ϕ_copy), boundary_conditions(ϕ_copy), InterpolationData(N, order, T)) + end + end + grid = mesh(mf) + pts = _sample_interface(grid, mf, cellindices(mf), upsample, maxiters, ftol) tree = KDTree(pts) - return NewtonSDF(itp, tree, pts, upsample, maxiters, xtol, ftol) -end - -""" - NewtonSDF(ϕ; order=3, kwargs...) - -Construct a [`NewtonSDF`](@ref) from a level set by first creating a piecewise polynomial -interpolant of the given `order`. Additional keyword arguments are forwarded to -`NewtonSDF(itp; ...)`. Works for both [`LevelSet`](@ref) and [`NarrowBandLevelSet`](@ref), -sampling only the relevant candidate cells in each case. -""" -function NewtonSDF(ϕ; order = 3, kwargs...) - itp = interpolate(ϕ, order) - return NewtonSDF(itp; kwargs...) + return NewtonSDF(mf, tree, pts, upsample, maxiters, xtol, ftol) end """ @@ -105,14 +115,14 @@ Rebuild `sdf` in place from the new level set `ϕ`, reusing the existing interpo buffers, upsample density, and solver tolerances. """ function update!(sdf::NewtonSDF, ϕ) - update!(sdf.itp, ϕ) - grid = mesh(sdf.itp.ϕ) - sdf.pts = _sample_interface(grid, sdf.itp, active_cells(sdf.itp.ϕ), sdf.upsample, sdf.maxiters, sdf.ftol) + copy!(sdf.meshfield, ϕ) + grid = mesh(sdf.meshfield) + sdf.pts = _sample_interface(grid, sdf.meshfield, cellindices(sdf.meshfield), sdf.upsample, sdf.maxiters, sdf.ftol) sdf.tree = KDTree(sdf.pts) return sdf end -function (sdf::NewtonSDF)(x, s = sign(sdf.itp(x))) +function (sdf::NewtonSDF)(x, s = sign(sdf.meshfield(x))) cp, _ = _closest_point_on_interface(sdf, x) return s * norm(x - cp) end @@ -128,29 +138,21 @@ neighbouring polynomial patch), a single retry is attempted using the best itera from the failed solve as a new seed on its own patch. """ function _closest_point_on_interface(sdf::NewtonSDF, x, max_retries = 3) - safeguard_dist = 1.5 * maximum(meshsize(mesh(sdf.itp.ϕ))) + safeguard_dist = 1.5 * maximum(meshsize(mesh(sdf.meshfield))) idx, _ = nn(sdf.tree, x) cp = sdf.pts[idx] - cell = compute_index(sdf.itp, cp) + cell = compute_index(sdf.meshfield, cp) converged = false for _ in 1:max_retries - p = make_interpolant(sdf.itp, cell) + p = make_interpolant(sdf.meshfield, cell) cp, converged = _closest_point(p, x, cp, sdf.maxiters, sdf.xtol, sdf.ftol, safeguard_dist) - new_cell = compute_index(sdf.itp, cp) + new_cell = compute_index(sdf.meshfield, cp) (converged || new_cell == cell) && break cell = new_cell end return cp, converged end -""" - active_cells(ϕ) - -Return the cell indices that are active for interface sampling. A cell is active if all -its corners are active indices. For a full-grid level set, this is all cells; for a narrow -band, only cells where all corners are within the band. Dispatch point for NewtonSDF construction. -""" -active_cells(ϕ) = cellindices(mesh(ϕ)) """ _sample_interface(grid, itp, cells, upsample, maxiters, ftol) @@ -326,16 +328,16 @@ function Base.show(io::IO, ::MIME"text/plain", r::NewtonReinitializer) end """ - reinitialize!(ϕ::LevelSet, r::NewtonReinitializer) + reinitialize!(ϕ::MeshField, r::NewtonReinitializer) Reinitialize the level set `ϕ` to a signed distance function in place using the Newton closest-point method. """ -function reinitialize!(ϕ::LevelSet, r::NewtonReinitializer) +function reinitialize!(ϕ::MeshField, r::NewtonReinitializer) sdf = NewtonSDF(ϕ; order = r.order, upsample = r.upsample, maxiters = r.maxiters, xtol = r.xtol, ftol = r.ftol) nfail = 0 for I in eachindex(ϕ) - x = mesh(ϕ)[I] + x = getnode(mesh(ϕ), I) cp, converged = _closest_point_on_interface(sdf, x) converged || (nfail += 1) ϕ[I] = sign(ϕ[I]) * norm(x - cp) @@ -348,11 +350,11 @@ end # Called at each time step with the current step count. Reinitializes only when # nsteps is a multiple of reinit_freq; the nothing method is a no-op. -reinitialize!(ϕ::LevelSet, ::Nothing, _) = ϕ -function reinitialize!(ϕ::LevelSet, r::NewtonReinitializer, nsteps::Int) +reinitialize!(ϕ::MeshField, ::Nothing, _) = ϕ +function reinitialize!(ϕ::MeshField, r::NewtonReinitializer, nsteps::Int) mod(nsteps, r.reinit_freq) == 0 || return ϕ return reinitialize!(ϕ, r) end # tomatoes tomatos ... -reinitialise!(ϕ::LevelSet, r::NewtonReinitializer) = reinitialize!(ϕ, r) +reinitialise!(ϕ::MeshField, r::NewtonReinitializer) = reinitialize!(ϕ, r) diff --git a/src/timestepping.jl b/src/timestepping.jl index 9b176e5..ca2fe3d 100644 --- a/src/timestepping.jl +++ b/src/timestepping.jl @@ -110,7 +110,7 @@ function Base.show(io::IO, ::MIME"text/plain", s::SemiImplicitI2OE) end # common integration logic -@noinline function _integrate!(ϕ::MeshField, integrator::TimeIntegrator, terms, reinit, tc, tf, Δt_max, log) +@noinline function _integrate!(ϕ::AbstractMeshField, integrator::TimeIntegrator, terms, reinit, tc, tf, Δt_max, log) src = ϕ buffers = _alloc_buffers(integrator, ϕ) α = cfl(integrator) @@ -269,8 +269,8 @@ function _integrate!( Δt_max, log, ) - # Check domain compatibility (FullDomain only for now) - domain(ϕ) isa FullDomain || throw(ArgumentError("SemiImplicitI2OE only supports FullDomain")) + # Check domain compatibility (full-domain MeshField only) + ϕ isa MeshField || throw(ArgumentError("SemiImplicitI2OE only supports full-domain MeshField")) _validate_i2oe_setup(ϕ, terms) term = only(terms) @@ -342,7 +342,7 @@ function _fill_advection_velocity_components!(out, term::AdvectionTerm{<:Functio N = ndims(ϕ) g = mesh(ϕ) for I in eachindex(ϕ) - vI = vel(g[I], t) + vI = vel(getnode(g, I), t) for dim in 1:N out[dim][I] = vI[dim] end @@ -511,9 +511,3 @@ function _i2oe_face_measure(Δ, dim) N == 1 && return one(eltype(Δ)) return prod(Δ[d] for d in 1:N if d != dim) end - -function _compute_terms(terms, ϕ, I, t) - return sum(terms) do term - return _compute_term(term, ϕ, I, t) - end -end diff --git a/src/velocityextension.jl b/src/velocityextension.jl index 0179b87..1f3e42a 100644 --- a/src/velocityextension.jl +++ b/src/velocityextension.jl @@ -1,5 +1,5 @@ """ - extend_along_normals!(F, ϕ::LevelSet; + extend_along_normals!(F, ϕ::MeshField; nb_iters = 50, cfl = 0.45, frozen = nothing, @@ -19,7 +19,7 @@ Reference from [Peng et al. 1999] """ function extend_along_normals!( F::AbstractArray{T, N}, - ϕ::LevelSet; + ϕ::MeshField; nb_iters::Integer = 50, cfl::Real = 0.45, frozen = nothing, @@ -40,8 +40,8 @@ function extend_along_normals!( else ntuple(_ -> (LinearExtrapolationBC(), LinearExtrapolationBC()), N) end - ϕw = has_boundary_conditions(ϕ) ? ϕ : add_boundary_conditions(ϕ, bc) - Fw = MeshField(F, mesh(ϕ), bc) + ϕw = has_boundary_conditions(ϕ) ? ϕ : _add_boundary_conditions(ϕ, bc) + Fw = MeshField(F, mesh(ϕ); bc = bc) Δ = minimum(meshsize(ϕ)) τ = cfl * Δ @@ -68,14 +68,14 @@ function extend_along_normals!( return F end -function extend_along_normals!(F::MeshField, ϕ::LevelSet; kwargs...) +function extend_along_normals!(F::MeshField, ϕ::MeshField; kwargs...) mesh(F) == mesh(ϕ) || throw(ArgumentError("F and ϕ must be defined on the same mesh")) extend_along_normals!(values(F), ϕ; kwargs...) return F end -function _normalize_frozen_mask(frozen, ϕ::LevelSet, interface_band, Δ) +function _normalize_frozen_mask(frozen, ϕ::MeshField, interface_band, Δ) vals = values(ϕ) if isnothing(frozen) return abs.(vals) .<= interface_band * Δ @@ -92,7 +92,7 @@ function _normalize_frozen_mask(frozen, ϕ::LevelSet, interface_band, Δ) return BitArray(frozen) end -function _signed_normal_components(ϕ::LevelSet, Δ, min_norm) +function _signed_normal_components(ϕ::MeshField, Δ, min_norm) N = ndims(ϕ) T = float(eltype(values(ϕ))) components = ntuple(_ -> Array{T}(undef, size(values(ϕ))), N) diff --git a/test/Project.toml b/test/Project.toml index 628af5f..b21d4f0 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,6 @@ [deps] DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" +ImplicitIntegration = "bc256489-3a69-4a66-afc4-127cc87e6182" LevelSetMethods = "ca3df947-0575-4ae6-b9a3-78df393a7451" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MMG_jll = "86086c02-e288-5929-a127-40944b0018b7" diff --git a/test/test-derivatives.jl b/test/test-derivatives.jl index 37cbc35..e16617c 100644 --- a/test/test-derivatives.jl +++ b/test/test-derivatives.jl @@ -8,9 +8,9 @@ using LevelSetMethods: D⁺, D⁻, D⁰, D2⁰, D2, weno5⁻, weno5⁺ # Exact derivatives: ∂_x = 3x²+y², ∂_y = 2xy, ∂_xx = 6x, ∂_yy = 2x, ∂_xy = 2y grid = CartesianGrid((-2.0, -2.0), (2.0, 2.0), (100, 50)) h = LevelSetMethods.meshsize(grid) -ϕ = LevelSet(v -> v[1]^3 + v[1] * v[2]^2, grid) +ϕ = MeshField(v -> v[1]^3 + v[1] * v[2]^2, grid) I = CartesianIndex(9, 7) -x, y = grid[I] +x, y = getnode(grid, I) @testset "First derivatives" begin exact = SVector(3x^2 + y^2, 2x * y) diff --git a/test/test-interpolation.jl b/test/test-interpolation.jl index 49cb91c..d920315 100644 --- a/test/test-interpolation.jl +++ b/test/test-interpolation.jl @@ -37,29 +37,27 @@ using Test grad_f(x) = SVector(2 * x[1], 4 * x[2]) hess_f(x) = SMatrix{2, 2}(2.0, 0.0, 0.0, 4.0) - ϕ = LevelSet(f, grid) - itp = interpolate(ϕ, 3) + ϕ = MeshField(f, grid; interp_order = 3) x_test = SVector(0.15, -0.25) - @test itp(x_test) ≈ f(x_test) atol = 1.0e-12 - I = LevelSetMethods.compute_index(itp, x_test) - p = LevelSetMethods.make_interpolant(itp, I) + @test ϕ(x_test) ≈ f(x_test) atol = 1.0e-12 + I = LevelSetMethods.compute_index(ϕ, x_test) + p = LevelSetMethods.make_interpolant(ϕ, I) @test LevelSetMethods.gradient(p, x_test) ≈ grad_f(x_test) atol = 1.0e-12 - @test check_allocs(itp, x_test) == 0 + @test check_allocs(ϕ, x_test) == 0 end @testset "Least Squares Approximation (K=2)" begin grid = CartesianGrid((-1.0, -1.0), (1.0, 1.0), (21, 21)) f(x) = x[1]^2 + 2 * x[2]^2 - 0.5 - ϕ = LevelSet(f, grid) - itp = interpolate(ϕ, 2) + ϕ = MeshField(f, grid; interp_order = 2) x_test = SVector(0.15, -0.25) - @test itp(x_test) ≈ f(x_test) atol = 1.0e-12 - I = LevelSetMethods.compute_index(itp, x_test) - p = LevelSetMethods.make_interpolant(itp, I) + @test ϕ(x_test) ≈ f(x_test) atol = 1.0e-12 + I = LevelSetMethods.compute_index(ϕ, x_test) + p = LevelSetMethods.make_interpolant(ϕ, I) @test LevelSetMethods.gradient(p, x_test) ≈ SVector(2 * 0.15, 4 * (-0.25)) atol = 1.0e-12 - @test check_allocs(itp, x_test) == 0 + @test check_allocs(ϕ, x_test) == 0 end @testset "Mesh Interpolation (3D)" begin @@ -67,41 +65,39 @@ using Test f(x) = x[1]^2 + x[2]^2 + x[3]^2 - 0.5 grad_f(x) = SVector(2 * x[1], 2 * x[2], 2 * x[3]) - ϕ = LevelSet(f, grid) - itp = interpolate(ϕ, 3) + ϕ = MeshField(f, grid; interp_order = 3) x_test = SVector(0.1, -0.2, 0.3) - @test itp(x_test) ≈ f(x_test) atol = 1.0e-12 - I = LevelSetMethods.compute_index(itp, x_test) - p = LevelSetMethods.make_interpolant(itp, I) + @test ϕ(x_test) ≈ f(x_test) atol = 1.0e-12 + I = LevelSetMethods.compute_index(ϕ, x_test) + p = LevelSetMethods.make_interpolant(ϕ, I) @test LevelSetMethods.gradient(p, x_test) ≈ grad_f(x_test) atol = 1.0e-12 - @test check_allocs(itp, x_test) == 0 + @test check_allocs(ϕ, x_test) == 0 end @testset "Convex Hull & proven_empty" begin grid = CartesianGrid((-1.0, -1.0), (1.0, 1.0), (20, 20)) f(x) = x[1] - ϕ = LevelSet(f, grid) - itp = interpolate(ϕ, 3) + ϕ = MeshField(f, grid; interp_order = 3) I_inside = CartesianIndex(1, 1) - m, M = LevelSetMethods.cell_extrema(itp, I_inside) + m, M = LevelSetMethods.cell_extrema(ϕ, I_inside) @test M < 0 - @test LevelSetMethods.proven_empty(itp, I_inside; surface = true) # no surface - @test !LevelSetMethods.proven_empty(itp, I_inside; surface = false) # fully inside + @test LevelSetMethods.proven_empty(ϕ, I_inside; surface = true) # no surface + @test !LevelSetMethods.proven_empty(ϕ, I_inside; surface = false) # fully inside I_outside = CartesianIndex(20, 20) - m, M = LevelSetMethods.cell_extrema(itp, I_outside) + m, M = LevelSetMethods.cell_extrema(ϕ, I_outside) @test m > 0 - @test LevelSetMethods.proven_empty(itp, I_outside; surface = true) - @test LevelSetMethods.proven_empty(itp, I_outside; surface = false) + @test LevelSetMethods.proven_empty(ϕ, I_outside; surface = true) + @test LevelSetMethods.proven_empty(ϕ, I_outside; surface = false) # cell crossing the interface I_interface = CartesianIndex(10, 4) - m, M = LevelSetMethods.cell_extrema(itp, I_interface) + m, M = LevelSetMethods.cell_extrema(ϕ, I_interface) @test m < 0 && M > 0 - @test !LevelSetMethods.proven_empty(itp, I_interface; surface = true) - @test !LevelSetMethods.proven_empty(itp, I_interface; surface = false) + @test !LevelSetMethods.proven_empty(ϕ, I_interface; surface = true) + @test !LevelSetMethods.proven_empty(ϕ, I_interface; surface = false) end end diff --git a/test/test-levelsetequation.jl b/test/test-levelsetequation.jl index 3312f94..1e02bf1 100644 --- a/test/test-levelsetequation.jl +++ b/test/test-levelsetequation.jl @@ -9,11 +9,11 @@ function _convergence_orders(errors, Ns) return [log(errors[i] / errors[i + 1]) / log(Ns[i + 1] / Ns[i]) for i in 1:(length(Ns) - 1)] end -# Max error between a NarrowBandLevelSet and a full-grid field, restricted to +# Max error between a MeshField and a full-grid field, restricted to # active nodes with |ϕ| < halfwidth/2 (near the interface where accuracy matters). function _nb_full_error(nb_s, full_s) γ = LSM.halfwidth(nb_s) - return maximum(LSM.active_indices(nb_s)) do I + return maximum(LSM.nodeindices(nb_s)) do I abs(nb_s[I]) < γ / 2 || return 0.0 abs(nb_s[I] - full_s[I]) end @@ -28,14 +28,14 @@ end Ns = [20, 40, 80] errors = map(Ns) do N grid = LSM.CartesianGrid((-1.0,), (1.0,), (N,)) - ϕ = LSM.LevelSet(x -> ϕ_exact(x, 0.0), grid) + ϕ = LSM.MeshField(x -> ϕ_exact(x, 0.0), grid) eq = LSM.LevelSetEquation(; terms = (LSM.AdvectionTerm((x, t) -> SVector(u)),), ic = ϕ, bc = PeriodicBC(), integrator = RK3(; cfl = 1.0e-2), ) integrate!(eq, tf) ϕ_out = current_state(eq) - maximum(I -> abs(ϕ_out[I] - ϕ_exact(grid[I], tf)), CartesianIndices(LSM.mesh(ϕ_out))) + maximum(I -> abs(ϕ_out[I] - ϕ_exact(getnode(grid, I), tf)), nodeindices(LSM.mesh(ϕ_out))) end @test all(≥(4.5), _convergence_orders(errors, Ns)) end @@ -50,15 +50,15 @@ end Ns = [30, 60, 120] errors = map(Ns) do N grid = LSM.CartesianGrid((-2.0, -2.0), (2.0, 2.0), (N, N)) - ϕ = LSM.LevelSet(x -> norm(x) - r0, grid) + ϕ = LSM.MeshField(x -> norm(x) - r0, grid) eq = LSM.LevelSetEquation(; terms = (LSM.NormalMotionTerm((x, t) -> v),), ic = ϕ, bc = ExtrapolationBC(2), integrator = RK3(), ) integrate!(eq, tf) ϕ_out = current_state(eq) - maximum(CartesianIndices(LSM.mesh(ϕ_out))) do I - x = grid[I] + maximum(nodeindices(LSM.mesh(ϕ_out))) do I + x = getnode(grid, I) (0.5 ≤ norm(x) ≤ 1.5) || return 0.0 abs(ϕ_out[I] - ϕ_exact(x)) end @@ -78,15 +78,15 @@ end Ns = [30, 60, 120] errors = map(Ns) do N grid = LSM.CartesianGrid((-2.0, -2.0), (2.0, 2.0), (N, N)) - ϕ = LSM.LevelSet(x -> norm(x) - r0, grid) + ϕ = LSM.MeshField(x -> norm(x) - r0, grid) eq = LSM.LevelSetEquation(; terms = (LSM.CurvatureTerm((x, t) -> b),), ic = ϕ, bc = ExtrapolationBC(2), integrator = RK3(), ) integrate!(eq, tf) ϕ_out = current_state(eq) - maximum(CartesianIndices(LSM.mesh(ϕ_out))) do I - x = grid[I] + maximum(nodeindices(LSM.mesh(ϕ_out))) do I + x = getnode(grid, I) (0.5 ≤ norm(x) ≤ 1.5) || return 0.0 abs(ϕ_out[I] - ϕ_exact(x)) end @@ -97,7 +97,7 @@ end @testset "Reinitialization inside LevelSetEquation" begin @testset "with RK2" begin grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (33, 33)) - ϕ = LSM.LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) + ϕ = LSM.MeshField(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) eq = LSM.LevelSetEquation(; terms = (LSM.AdvectionTerm((x, t) -> @SVector [1.0, 0.0]),), ic = ϕ, bc = LSM.PeriodicBC(), reinit = 2, @@ -108,7 +108,7 @@ end @testset "with SemiImplicitI2OE" begin grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (33, 33)) - ϕ = LSM.LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) + ϕ = LSM.MeshField(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) eq = LSM.LevelSetEquation(; terms = (LSM.AdvectionTerm((x, t) -> @SVector [1.0, 0.0]),), ic = ϕ, bc = LSM.PeriodicBC(), integrator = LSM.SemiImplicitI2OE(), reinit = 2, @@ -119,10 +119,10 @@ end @testset "NarrowBand integrate! — advection matches full grid" begin grid = LSM.CartesianGrid((-2.0, -2.0), (2.0, 2.0), (60, 60)) - ϕ = LSM.LevelSet(x -> norm(x) - 0.5, grid) + ϕ = LSM.MeshField(x -> norm(x) - 0.5, grid) 𝐮 = (x, t) -> SVector(1.0, 0.0) bc = ExtrapolationBC(2) - eq_nb = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = NarrowBandLevelSet(ϕ; nlayers = 5), bc, reinit = 1) + eq_nb = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = NarrowBandMeshField(ϕ; nlayers = 5), bc, reinit = 1) eq_full = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = deepcopy(ϕ), bc, reinit = 1) integrate!(eq_full, 0.1); integrate!(eq_nb, 0.1) @test _nb_full_error(current_state(eq_nb), current_state(eq_full)) < 1.0e-3 @@ -130,19 +130,19 @@ end @testset "NarrowBand integrate! — advection with reinitialization" begin grid = LSM.CartesianGrid((-2.0, -2.0), (2.0, 2.0), (60, 60)) - ϕ = LSM.LevelSet(x -> norm(x) - 0.5, grid) + ϕ = LSM.MeshField(x -> norm(x) - 0.5, grid) eq_nb = LevelSetEquation(; terms = (AdvectionTerm((x, t) -> SVector(1.0, 0.0)),), - ic = NarrowBandLevelSet(ϕ), bc = ExtrapolationBC(2), reinit = NewtonReinitializer(), + ic = NarrowBandMeshField(ϕ), bc = ExtrapolationBC(2), reinit = NewtonReinitializer(), ) integrate!(eq_nb, 0.1) nb_s = current_state(eq_nb) γ = LSM.halfwidth(nb_s) exact_sdf = x -> norm(x - SVector(0.1, 0.0)) - 0.5 - @test length(LSM.active_indices(nb_s)) > 0 - @test maximum(LSM.active_indices(nb_s)) do I + @test length(LSM.nodeindices(nb_s)) > 0 + @test maximum(LSM.nodeindices(nb_s)) do I abs(nb_s[I]) < γ / 2 || return 0.0 - abs(nb_s[I] - exact_sdf(grid[I])) + abs(nb_s[I] - exact_sdf(getnode(grid, I))) end < 0.01 end @@ -153,7 +153,7 @@ end r0, θ0, α = 0.5, -π / 3, π / 100 R = [cos(α) -sin(α); sin(α) cos(α)] M = R * [1 / 0.06^2 0; 0 1 / (4π^2)] * R' - ϕ = LSM.LevelSet(grid) do (x, y) + ϕ = LSM.MeshField(grid) do (x, y) r, θ = sqrt(x^2 + y^2), atan(y, x) minimum(0:4) do i v = [r - r0; θ + (2i - 4) * π - θ0] @@ -161,7 +161,7 @@ end end end b, reinit, bc = (x, t) -> -0.1, NewtonReinitializer(), ExtrapolationBC(2) - eq_nb = LevelSetEquation(; ic = NarrowBandLevelSet(ϕ), bc, terms = (CurvatureTerm(b),), reinit) + eq_nb = LevelSetEquation(; ic = NarrowBandMeshField(ϕ), bc, terms = (CurvatureTerm(b),), reinit) eq_full = LevelSetEquation(; ic = deepcopy(ϕ), bc, terms = (CurvatureTerm(b),), reinit) integrate!(eq_full, 0.1); integrate!(eq_nb, 0.1) @test _nb_full_error(current_state(eq_nb), current_state(eq_full)) < 0.05 @@ -169,27 +169,27 @@ end @testset "NarrowBand integrate! — full rotation" begin grid = LSM.CartesianGrid((-2.0, -2.0), (2.0, 2.0), (40, 40)) - ϕ = LSM.LevelSet(x -> norm(x - SVector(0.8, 0.0)) - 0.5, grid) + ϕ = LSM.MeshField(x -> norm(x - SVector(0.8, 0.0)) - 0.5, grid) 𝐮, bc = (x, t) -> SVector(-x[2], x[1]), ExtrapolationBC(2) - eq_nb = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = NarrowBandLevelSet(ϕ), bc, reinit = NewtonReinitializer()) + eq_nb = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = NarrowBandMeshField(ϕ), bc, reinit = NewtonReinitializer()) eq_full = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = deepcopy(ϕ), bc) integrate!(eq_full, 2π); integrate!(eq_nb, 2π) nb_s = current_state(eq_nb) - @test length(LSM.active_indices(nb_s)) > 0 + @test length(LSM.nodeindices(nb_s)) > 0 @test _nb_full_error(nb_s, current_state(eq_full)) < 0.01 end @testset "NarrowBand integrate! — star rotation" begin grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (40, 40)) - ϕ₀ = LevelSet(grid) do (x, y) + ϕ₀ = MeshField(grid) do (x, y) r, θ = sqrt(x^2 + y^2), atan(y, x) - π / 2 r - (0.5 + 0.1 * cos(5θ)) end 𝐮, bc = (x, t) -> SVector(-x[2], x[1]), ExtrapolationBC(2) - eq_nb = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = NarrowBandLevelSet(ϕ₀), bc, reinit = NewtonReinitializer()) + eq_nb = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = NarrowBandMeshField(ϕ₀), bc, reinit = NewtonReinitializer()) eq_full = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = deepcopy(ϕ₀), bc) integrate!(eq_full, pi); integrate!(eq_nb, pi) nb_s = current_state(eq_nb) - @test length(LSM.active_indices(nb_s)) > 0 + @test length(LSM.nodeindices(nb_s)) > 0 @test _nb_full_error(nb_s, current_state(eq_full)) < 0.05 end diff --git a/test/test-levelsetterms.jl b/test/test-levelsetterms.jl index bba9ad0..2583d1b 100644 --- a/test/test-levelsetterms.jl +++ b/test/test-levelsetterms.jl @@ -6,7 +6,7 @@ import LevelSetMethods as LSM @testset "AdvectionTerm CFL" begin grid = LSM.CartesianGrid((-1.0,), (1.0,), (100,)) - ϕ = LSM.LevelSet(x -> x[1], grid) + ϕ = LSM.MeshField(x -> x[1], grid) Δx = LSM.meshsize(ϕ, 1) term = LSM.AdvectionTerm((x, t) -> SVector(2.0)) @test LSM.compute_cfl((term,), ϕ, 0.0) ≈ Δx / 2.0 @@ -14,7 +14,7 @@ end @testset "CurvatureTerm CFL" begin grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (50, 50)) - ϕ = LSM.LevelSet(x -> norm(x) - 0.5, grid) + ϕ = LSM.MeshField(x -> norm(x) - 0.5, grid) Δx = minimum(LSM.meshsize(ϕ)) b = 0.5 term = LSM.CurvatureTerm((x, t) -> b) @@ -23,7 +23,7 @@ end @testset "NormalMotionTerm CFL" begin grid = LSM.CartesianGrid((-1.0,), (1.0,), (100,)) - ϕ = LSM.LevelSet(x -> x[1], grid) + ϕ = LSM.MeshField(x -> x[1], grid) Δx = LSM.meshsize(ϕ, 1) v = 3.0 term = LSM.NormalMotionTerm((x, t) -> v) @@ -34,7 +34,7 @@ end # ϕ = 2*(x - 0.3): correct zero set but |∇ϕ| = 2 ≠ 1. # After pseudo-time marching it should converge to x - 0.3. grid = LSM.CartesianGrid((-1.0,), (1.0,), (101,)) - ϕ = LSM.LevelSet(x -> 2 * (x[1] - 0.3), grid) + ϕ = LSM.MeshField(x -> 2 * (x[1] - 0.3), grid) eq = LSM.LevelSetEquation(; terms = (LSM.EikonalReinitializationTerm(ϕ),), ic = deepcopy(ϕ), @@ -42,8 +42,8 @@ end ) integrate!(eq, 2.0) ϕ_out = LSM.current_state(eq) - ϕ_exact = LSM.LevelSet(x -> x[1] - 0.3, grid) - err = maximum(CartesianIndices(LSM.mesh(ϕ_out))) do I + ϕ_exact = LSM.MeshField(x -> x[1] - 0.3, grid) + err = maximum(nodeindices(LSM.mesh(ϕ_out))) do I abs(ϕ_out[I]) > 0.5 && return 0.0 abs(ϕ_out[I] - ϕ_exact[I]) end diff --git a/test/test-meshes.jl b/test/test-meshes.jl index 6ec74f2..8896ca6 100644 --- a/test/test-meshes.jl +++ b/test/test-meshes.jl @@ -7,7 +7,7 @@ using StaticArrays a, b = (-1, 0), (1, 3) grid = CartesianGrid(a, b, (nx, ny)) @test size(grid) === (nx, ny) - @test length(CartesianIndices(grid)) == nx * ny - @test grid[1, 1] == SVector(a[1], a[2]) - @test grid[nx, ny] == SVector(b[1], b[2]) + @test length(nodeindices(grid)) == nx * ny + @test getnode(grid, 1, 1) == SVector(a[1], a[2]) + @test getnode(grid, nx, ny) == SVector(b[1], b[2]) end diff --git a/test/test-meshfield.jl b/test/test-meshfield.jl index 2102b77..dd48752 100644 --- a/test/test-meshfield.jl +++ b/test/test-meshfield.jl @@ -9,27 +9,27 @@ import LevelSetMethods as LSM f = x -> x[1]^2 + x[2]^2 - 0.5 ϕ = MeshField(f, grid) @test LSM.mesh(ϕ) === grid - @test LSM.domain(ϕ) isa LSM.FullDomain + @test ϕ isa MeshField @test !LSM.has_boundary_conditions(ϕ) @test ndims(ϕ) == 2 @test size(values(ϕ)) == (10, 5) - @test ϕ[3, 2] ≈ f(grid[3, 2]) + @test ϕ[3, 2] ≈ f(getnode(grid, 3, 2)) @test LSM.meshsize(ϕ) == LSM.meshsize(grid) end -@testset "add_boundary_conditions" begin +@testset "_add_boundary_conditions" begin grid = CartesianGrid((0.0,), (1.0,), (5,)) ϕ = MeshField(x -> x[1], grid) @test !LSM.has_boundary_conditions(ϕ) - ϕ_bc = LSM.add_boundary_conditions(ϕ, ((NeumannBC(), NeumannBC()),)) + ϕ_bc = LSM._add_boundary_conditions(ϕ, ((NeumannBC(), NeumannBC()),)) @test LSM.has_boundary_conditions(ϕ_bc) @test values(ϕ_bc) === values(ϕ) # underlying data is aliased end @testset "copy!" begin grid = CartesianGrid((0.0, 0.0), (1.0, 1.0), (5, 5)) - ϕ = LevelSet(x -> x[1] + x[2], grid) - ψ = LevelSet(x -> 0.0, grid) + ϕ = MeshField(x -> x[1] + x[2], grid) + ψ = MeshField(x -> 0.0, grid) copy!(ψ, ϕ) @test values(ψ) == values(ϕ) @test values(ψ) !== values(ϕ) # copied, not aliased @@ -38,15 +38,14 @@ end @testset "update_bcs!" begin grid = CartesianGrid((0.0,), (1.0,), (5,)) bc = DirichletBC((x, t) -> t) - ϕ = MeshField(x -> 0.0, grid, ((bc, NeumannBC()),)) + ϕ = MeshField(x -> 0.0, grid; bc = ((bc, NeumannBC()),)) LSM.update_bcs!(ϕ, 3.0) @test bc.t == 3.0 end @testset "Type aliases" begin grid = CartesianGrid((-1.0, -1.0), (1.0, 1.0), (5, 5)) - ϕ = LevelSet(x -> norm(x) - 0.5, grid) - @test ϕ isa LevelSet{2, Float64} + ϕ = MeshField(x -> norm(x) - 0.5, grid) @test ϕ isa MeshField end @@ -57,7 +56,7 @@ end grid = CartesianGrid(a, b, n) vals = rand(n...) bcs = ((PeriodicBC(), PeriodicBC()), (PeriodicBC(), PeriodicBC())) - mf = MeshField(vals, grid, bcs) + mf = MeshField(vals, grid; bc = bcs) @test mf[1, 1] == vals[1, 1] @test mf[1, 0] == vals[1, 4] # wraps around in dim 2 @test mf[11, 5] == mf[2, 5] # wraps around in dim 1 @@ -73,7 +72,7 @@ end bcs = ((ExtrapolationBC(P), ExtrapolationBC(P)),) for k in 0:P f = x -> x[1]^k - ϕ_bc = LSM.add_boundary_conditions(MeshField(f, grid), bcs) + ϕ_bc = LSM._add_boundary_conditions(MeshField(f, grid), bcs) # node i is at x = a + (i-1)*h, so ghost 1-j is at x = a - j*h # and ghost n+j is at x = b + j*h for j in 1:(P + 1) @@ -94,11 +93,12 @@ end bcs2 = ((ExtrapolationBC(P), ExtrapolationBC(P)), (ExtrapolationBC(P), ExtrapolationBC(P))) for j in 0:P, k in 0:P f = x -> x[1]^j * x[2]^k - ϕ_bc = LSM.add_boundary_conditions(MeshField(f, grid2), bcs2) + ϕ_bc = LSM._add_boundary_conditions(MeshField(f, grid2), bcs2) # interior ghost in dim 1, in-bounds in dim 2 - @test ϕ_bc[0, 3] ≈ f((a1 - h1, grid2[1, 3][2])) atol = 1.0e-10 - @test ϕ_bc[n1 + 1, 3] ≈ f((b1 + h1, grid2[1, 3][2])) atol = 1.0e-10 + @test ϕ_bc[0, 3] ≈ f((a1 - h1, getnode(grid2, 1, 3)[2])) atol = 1.0e-10 + @test ϕ_bc[n1 + 1, 3] ≈ f((b1 + h1, getnode(grid2, 1, 3)[2])) atol = 1.0e-10 # corner ghost (out-of-bounds in both dims) + @test ϕ_bc[0, 0] ≈ f((a1 - h1, a2 - h2)) atol = 1.0e-10 @test ϕ_bc[n1 + 1, n2 + 1] ≈ f((b1 + h1, b2 + h2)) atol = 1.0e-10 end @@ -110,7 +110,7 @@ end h = LSM.meshsize(grid)[1] bc = DirichletBC((x, t) -> x[1] + t) bcs = ((bc, NeumannBC()),) - ϕ = MeshField(x -> 0.0, grid, bcs) + ϕ = MeshField(x -> 0.0, grid; bc = bcs) @test ϕ[0] ≈ -h + 0.0 # t = 0 initially LSM.update_bc!(bc, 2.0) @test ϕ[0] ≈ -h + 2.0 # t updated to 2.0 diff --git a/test/test-narrow-band.jl b/test/test-narrow-band.jl index 29c2027..f10037c 100644 --- a/test/test-narrow-band.jl +++ b/test/test-narrow-band.jl @@ -7,25 +7,25 @@ using LevelSetMethods: D⁺, D⁻, D⁰, D2⁰, D2, weno5⁻, weno5⁺ @testset "Construction" begin grid = LSM.CartesianGrid((-2.0, -2.0), (2.0, 2.0), (100, 100)) - ϕ = LSM.LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) + ϕ = LSM.MeshField(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) halfwidth = 0.3 - nb = NarrowBandLevelSet(ϕ, halfwidth; reinitialize = false) + nb = NarrowBandMeshField(ϕ, halfwidth; reinitialize = false) - @test 0 < length(LSM.active_indices(nb)) < length(grid) - @test all(I -> abs(nb[I]) < LSM.halfwidth(nb), LSM.active_indices(nb)) - @test all(I -> nb[I] ≈ ϕ[I], LSM.active_indices(nb)) + @test 0 < length(LSM.nodeindices(nb)) < length(grid) + @test all(I -> abs(nb[I]) < LSM.halfwidth(nb), LSM.nodeindices(nb)) + @test all(I -> nb[I] ≈ ϕ[I], LSM.nodeindices(nb)) # Check that active points satisfy requirements - active_idxs = LSM.active_indices(nb) + active_idxs = LSM.nodeindices(nb) @test all(I -> abs(nb[I]) < halfwidth, active_idxs) - inactive_idxs = setdiff(CartesianIndices(grid), active_idxs) + inactive_idxs = setdiff(nodeindices(grid), active_idxs) @test all(I -> abs(ϕ[I]) >= halfwidth, inactive_idxs) # Automatic halfwidth via nlayers - nb2 = NarrowBandLevelSet(ϕ; nlayers = 8) + nb2 = NarrowBandMeshField(ϕ; nlayers = 8) Δx = minimum(LSM.meshsize(grid)) @test LSM.halfwidth(nb2) ≈ 8 * Δx - @test length(LSM.active_indices(nb2)) > 0 + @test length(LSM.nodeindices(nb2)) > 0 end @@ -33,8 +33,8 @@ end # Bilinear extrapolation is exact for linear functions. grid = LSM.CartesianGrid((-2.0, -2.0), (2.0, 2.0), (100, 100)) f = x -> 3x[1] - 2x[2] + x[1] * x[2] + 1.0 - nb = NarrowBandLevelSet(f, grid, 0.3) - active_idxs = LSM.active_indices(nb) + nb = NarrowBandMeshField(f, grid, 0.3) + active_idxs = LSM.nodeindices(nb) Imin = map(d -> minimum(I[d] for I in active_idxs), (1, 2)) |> CartesianIndex Imax = map(d -> maximum(I[d] for I in active_idxs), (1, 2)) |> CartesianIndex k = 5 @@ -51,12 +51,12 @@ end @testset "Derivatives match full grid" begin grid = LSM.CartesianGrid((-2.0, -2.0), (2.0, 2.0), (100, 100)) - ϕ = LSM.LevelSet(x -> x[1]^2 + x[2]^2 - 1, grid) + ϕ = LSM.MeshField(x -> x[1]^2 + x[2]^2 - 1, grid) bc = LSM._normalize_bc(LSM.ExtrapolationBC{2}(), 2) - ϕ_bc = LSM.add_boundary_conditions(ϕ, bc) - nb = NarrowBandLevelSet(ϕ_bc, 0.5; reinitialize = false) + ϕ_bc = LSM._add_boundary_conditions(ϕ, bc) + nb = NarrowBandMeshField(ϕ_bc, 0.5; reinitialize = false) - best_I = argmin(I -> abs(nb[I]), LSM.active_indices(nb)) + best_I = argmin(I -> abs(nb[I]), LSM.nodeindices(nb)) @testset "First order" begin for op in (D⁺, D⁻, D⁰) @@ -86,16 +86,13 @@ end @testset "Interpolation" begin grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (50, 50)) f(x) = x[1]^2 + 2 * x[2]^2 - 0.5 - ϕ = LSM.LevelSet(f, grid) - bc = LSM._normalize_bc(LSM.ExtrapolationBC{2}(), 2) - ϕ_bc = LSM.add_boundary_conditions(ϕ, bc) - nb = NarrowBandLevelSet(ϕ_bc, 0.4; reinitialize = false) - - itp_nb = interpolate(nb, 3) - itp_full = interpolate(ϕ_bc, 3) + # Full-grid field with interpolation + ϕ_full = LSM.MeshField(f, grid; bc = LSM.ExtrapolationBC{2}(), interp_order = 3) + # Narrow-band field: build from ϕ_full (which already has BCs and interp_order) + nb = NarrowBandMeshField(ϕ_full, 0.4; reinitialize = false) for x in [SVector(0.5, 0.0), SVector(0.0, 0.5), SVector(0.3, 0.3)] - @test itp_nb(x) ≈ itp_full(x) + @test nb(x) ≈ ϕ_full(x) end end @@ -103,10 +100,10 @@ end grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (50, 50)) r = 0.5 exact_sdf(x) = norm(x) - r - ϕ = LSM.LevelSet(exact_sdf, grid) + ϕ = LSM.MeshField(exact_sdf, grid) bc = LSM._normalize_bc(LSM.ExtrapolationBC{2}(), 2) - ϕ_bc = LSM.add_boundary_conditions(ϕ, bc) - nb = NarrowBandLevelSet(ϕ_bc, 0.3) + ϕ_bc = LSM._add_boundary_conditions(ϕ, bc) + nb = NarrowBandMeshField(ϕ_bc, 0.3) sdf = LSM.NewtonSDF(nb; upsample = 4) @test sdf(SVector(r, 0.0)) ≈ 0.0 atol = 2.0e-5 @@ -122,59 +119,59 @@ end grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (50, 50)) r = 0.5 exact_sdf(x) = norm(x) - r - ϕ = LSM.LevelSet(exact_sdf, grid) + ϕ = LSM.MeshField(exact_sdf, grid) bc = LSM._normalize_bc(LSM.ExtrapolationBC{2}(), 2) - ϕ_bc = LSM.add_boundary_conditions(ϕ, bc) - nb = NarrowBandLevelSet(ϕ_bc, 0.3) - n_before = length(LSM.active_indices(nb)) + ϕ_bc = LSM._add_boundary_conditions(ϕ, bc) + nb = NarrowBandMeshField(ϕ_bc, 0.3) + n_before = length(LSM.nodeindices(nb)) sdf = LSM.NewtonSDF(nb; upsample = 4) LSM.rebuild_band!(nb, sdf) - @test length(LSM.active_indices(nb)) > 0 - max_err = maximum(LSM.active_indices(nb)) do I - x = grid[I] + @test length(LSM.nodeindices(nb)) > 0 + max_err = maximum(LSM.nodeindices(nb)) do I + x = getnode(grid, I) abs(nb[I] - exact_sdf(x)) end @test max_err < 1.0e-5 - @test abs(length(LSM.active_indices(nb)) - n_before) < 0.1 * n_before + @test abs(length(LSM.nodeindices(nb)) - n_before) < 0.1 * n_before end @testset "Copy behavior" begin grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (30, 30)) - ϕ = LSM.LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) + ϕ = LSM.MeshField(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) bc = LSM._normalize_bc(LSM.LinearExtrapolationBC(), 2) - ϕ_bc = LSM.add_boundary_conditions(ϕ, bc) + ϕ_bc = LSM._add_boundary_conditions(ϕ, bc) - nb1 = NarrowBandLevelSet(ϕ_bc, 0.2; reinitialize = false) - nb2 = NarrowBandLevelSet(ϕ_bc, 0.4; reinitialize = false) + nb1 = NarrowBandMeshField(ϕ_bc, 0.2; reinitialize = false) + nb2 = NarrowBandMeshField(ϕ_bc, 0.4; reinitialize = false) copy!(nb2, nb1) - @test LSM.active_indices(nb1) == LSM.active_indices(nb2) - @test all(I -> nb1[I] ≈ nb2[I], LSM.active_indices(nb1)) + @test LSM.nodeindices(nb1) == LSM.nodeindices(nb2) + @test all(I -> nb1[I] ≈ nb2[I], LSM.nodeindices(nb1)) end @testset "eachindex iteration" begin grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (30, 30)) - ϕ = LSM.LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) + ϕ = LSM.MeshField(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) bc = LSM._normalize_bc(LSM.LinearExtrapolationBC(), 2) - ϕ_bc = LSM.add_boundary_conditions(ϕ, bc) - nb = NarrowBandLevelSet(ϕ_bc, 0.3; reinitialize = false) + ϕ_bc = LSM._add_boundary_conditions(ϕ, bc) + nb = NarrowBandMeshField(ϕ_bc, 0.3; reinitialize = false) - @test collect(eachindex(nb)) == collect(LSM.active_indices(nb)) + @test collect(eachindex(nb)) == collect(LSM.nodeindices(nb)) end @testset "Extrapolation at band boundary" begin grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (50, 50)) - ϕ = LSM.LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) + ϕ = LSM.MeshField(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) bc = LSM._normalize_bc(LSM.LinearExtrapolationBC(), 2) - ϕ_bc = LSM.add_boundary_conditions(ϕ, bc) - nb = NarrowBandLevelSet(ϕ_bc, 0.3) + ϕ_bc = LSM._add_boundary_conditions(ϕ, bc) + nb = NarrowBandMeshField(ϕ_bc, 0.3) # Find an in-grid node outside the band non_band = findfirst( I -> !haskey(values(nb), I) && all(s -> I[s] in axes(nb)[s], 1:2), - CartesianIndices(grid) + nodeindices(grid) ) @test non_band !== nothing @@ -189,14 +186,14 @@ end # guaranteed for non-linear functions, but sign correctness is (no spurious zeros). grid = LSM.CartesianGrid((-2.0, -2.0), (2.0, 2.0), (40, 40)) f(x) = norm(x) - 0.5 # circle SDF with halfwidth 0.2 → both ± nodes outside band - ϕ = LSM.LevelSet(f, grid) + ϕ = LSM.MeshField(f, grid) bc = LSM._normalize_bc(LSM.ExtrapolationBC{2}(), 2) - ϕ_bc = LSM.add_boundary_conditions(ϕ, bc) - nb = NarrowBandLevelSet(ϕ_bc, 0.2; reinitialize = false) + ϕ_bc = LSM._add_boundary_conditions(ϕ, bc) + nb = NarrowBandMeshField(ϕ_bc, 0.2; reinitialize = false) vals = values(nb) N = ndims(nb) - for I in LSM.active_indices(nb) + for I in LSM.nodeindices(nb) for d in 1:N Im = CartesianIndex(ntuple(i -> i == d ? I[i] - 1 : I[i], N)) Ip = CartesianIndex(ntuple(i -> i == d ? I[i] + 1 : I[i], N)) @@ -208,10 +205,10 @@ end @testset "Periodic BC" begin grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (50, 50)) - ϕ = LSM.LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) + ϕ = LSM.MeshField(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid) bc = LSM._normalize_bc(LSM.PeriodicBC(), 2) - ϕ_bc = LSM.add_boundary_conditions(ϕ, bc) - nb = NarrowBandLevelSet(ϕ_bc, 0.3) + ϕ_bc = LSM._add_boundary_conditions(ϕ, bc) + nb = NarrowBandMeshField(ϕ_bc, 0.3) I = CartesianIndex(0, 25) @test isfinite(nb[I]) @@ -222,14 +219,14 @@ end # preserves sign for all in-grid out-of-band nodes. grid = LSM.CartesianGrid((-2.0, -2.0, -2.0), (2.0, 2.0, 2.0), (20, 20, 20)) f(x) = norm(x) - 0.5 # sphere SDF with halfwidth 0.2 → both ± nodes outside band - ϕ = LSM.LevelSet(f, grid) + ϕ = LSM.MeshField(f, grid) bc = LSM._normalize_bc(LSM.ExtrapolationBC{2}(), 3) - ϕ_bc = LSM.add_boundary_conditions(ϕ, bc) - nb = NarrowBandLevelSet(ϕ_bc, 0.2; reinitialize = false) + ϕ_bc = LSM._add_boundary_conditions(ϕ, bc) + nb = NarrowBandMeshField(ϕ_bc, 0.2; reinitialize = false) vals = values(nb) N = ndims(nb) - for I in LSM.active_indices(nb) + for I in LSM.nodeindices(nb) for d in 1:N Im = CartesianIndex(ntuple(i -> i == d ? I[i] - 1 : I[i], N)) Ip = CartesianIndex(ntuple(i -> i == d ? I[i] + 1 : I[i], N)) @@ -243,7 +240,7 @@ end grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (100, 100)) d = 1; r0 = 0.5; θ0 = -π / 3; α = π / 100.0 R = [cos(α) -sin(α); sin(α) cos(α)]; M = R * [1 / 0.06^2 0; 0 1 / (4π^2)] * R' - ϕ = LSM.LevelSet(grid) do (x, y) + ϕ = LSM.MeshField(grid) do (x, y) r = sqrt(x^2 + y^2); θ = atan(y, x); res = 1.0e30 for i in 0:4 θ1 = θ + (2i - 4) * π; v = [r - r0; θ1 - θ0] @@ -252,24 +249,24 @@ end res end bc = LSM._normalize_bc(LSM.ExtrapolationBC{2}(), 2) - ϕ_bc = LSM.add_boundary_conditions(ϕ, bc) + ϕ_bc = LSM._add_boundary_conditions(ϕ, bc) - nb = NarrowBandLevelSet(ϕ_bc; nlayers = 6) - @test length(LSM.active_indices(nb)) > 3000 + nb = NarrowBandMeshField(ϕ_bc; nlayers = 6) + @test length(LSM.nodeindices(nb)) > 3000 end @testset "3D narrow band" begin grid = LSM.CartesianGrid((-1.0, -1.0, -1.0), (1.0, 1.0, 1.0), (25, 25, 25)) - ϕ = LSM.LevelSet(x -> norm(x) - 0.45, grid) + ϕ = LSM.MeshField(x -> norm(x) - 0.45, grid) bc = LSM._normalize_bc(LSM.ExtrapolationBC{2}(), 3) - ϕ_bc = LSM.add_boundary_conditions(ϕ, bc) - nb = NarrowBandLevelSet(ϕ_bc, 0.3; reinitialize = false) + ϕ_bc = LSM._add_boundary_conditions(ϕ, bc) + nb = NarrowBandMeshField(ϕ_bc, 0.3; reinitialize = false) @test ndims(nb) == 3 - @test length(LSM.active_indices(nb)) > 0 - @test length(LSM.active_indices(nb)) < length(grid) + @test length(LSM.nodeindices(nb)) > 0 + @test length(LSM.nodeindices(nb)) < length(grid) - best_I = argmin(I -> abs(nb[I]), LSM.active_indices(nb)) + best_I = argmin(I -> abs(nb[I]), LSM.nodeindices(nb)) for dim in 1:3 @test D⁰(nb, best_I, dim) ≈ D⁰(ϕ_bc, best_I, dim) end diff --git a/test/test-quadrature.jl b/test/test-quadrature.jl new file mode 100644 index 0000000..8a82d9f --- /dev/null +++ b/test/test-quadrature.jl @@ -0,0 +1,53 @@ +using LevelSetMethods +using ImplicitIntegration +using Test + +_total(quads, f = x -> 1.0) = sum(integrate(f, q) for (_, q) in quads) + +@testset "ImplicitIntegration Quadrature" begin + @testset "2D circle" begin + R = 0.5 + grid = CartesianGrid((-1.0, -1.0), (1.0, 1.0), (21, 21)) + ϕ = MeshField(x -> x[1]^2 + x[2]^2 - R^2, grid; interp_order = 3) + @test _total(quadrature(ϕ; order = 4, surface = false)) ≈ π * R^2 atol = 1.0e-4 + @test _total(quadrature(ϕ; order = 4, surface = true)) ≈ 2π * R atol = 1.0e-3 + end + + @testset "2D ellipse" begin + a, b = 0.6, 0.3 + grid = CartesianGrid((-1.0, -1.0), (1.0, 1.0), (41, 41)) + ϕ = MeshField(x -> (x[1] / a)^2 + (x[2] / b)^2 - 1.0, grid; interp_order = 3) + h = ((a - b) / (a + b))^2 + peri_approx = π * (a + b) * (1 + 3h / (10 + sqrt(4 - 3h))) + @test _total(quadrature(ϕ; order = 4, surface = false)) ≈ π * a * b rtol = 1.0e-3 + @test _total(quadrature(ϕ; order = 4, surface = true)) ≈ peri_approx rtol = 1.0e-3 + end + + @testset "3D sphere" begin + R = 0.5 + grid = CartesianGrid((-1.0, -1.0, -1.0), (1.0, 1.0, 1.0), (11, 11, 11)) + ϕ = MeshField(x -> x[1]^2 + x[2]^2 + x[3]^2 - R^2, grid; interp_order = 3) + @test _total(quadrature(ϕ; order = 2, surface = false)) ≈ 4π / 3 * R^3 atol = 1.0e-3 + @test _total(quadrature(ϕ; order = 2, surface = true)) ≈ 4π * R^2 atol = 1.0e-2 + end + + @testset "3D ellipsoid" begin + a, b, c = 0.61, 0.37, 0.29 + grid = CartesianGrid((-1.0, -1.0, -1.0), (1.0, 1.0, 1.0), (21, 21, 21)) + ϕ = MeshField(grid; interp_order = 3) do x + (x[1] / a)^2 + (x[2] / b)^2 + (x[3] / c)^2 - 1.0 + end + @test _total(quadrature(ϕ; order = 3, surface = false)) ≈ (4 / 3) * π * a * b * c rtol = 1.0e-3 + end + + @testset "Narrow band" begin + R = 0.5 + grid = CartesianGrid((-1.0, -1.0), (1.0, 1.0), (41, 41)) + ϕ_full = MeshField(x -> x[1]^2 + x[2]^2 - R^2, grid; interp_order = 3) + ϕ_nb = NarrowBandMeshField(ϕ_full, 0.3; reinitialize = false) + @test _total(quadrature(ϕ_full; order = 4, surface = false)) ≈ + _total(quadrature(ϕ_nb; order = 4, surface = false)) rtol = 1.0e-10 + @test _total(quadrature(ϕ_full; order = 4, surface = true)) ≈ + _total(quadrature(ϕ_nb; order = 4, surface = true)) rtol = 1.0e-10 + end +end diff --git a/test/test-reinitializer.jl b/test/test-reinitializer.jl index b3a6b14..1321f72 100644 --- a/test/test-reinitializer.jl +++ b/test/test-reinitializer.jl @@ -2,7 +2,7 @@ using Test using LinearAlgebra using StaticArrays import LevelSetMethods as LSM -using LevelSetMethods: NewtonSDF, update!, get_sample_points, interpolate +using LevelSetMethods: NewtonSDF, update!, get_sample_points function check_allocs(f, x) f(x) # warmup @@ -14,7 +14,7 @@ end grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (50, 50)) r = 0.5 exact_sdf(x) = norm(x) - r - ϕ = LSM.LevelSet(exact_sdf, grid) + ϕ = LSM.MeshField(exact_sdf, grid) sdf = NewtonSDF(ϕ; upsample = 4) @test sdf isa LSM.AbstractSignedDistanceFunction @@ -25,9 +25,9 @@ end @test sdf(SVector(1.0, 0.0)) ≈ 1 - r atol = 2.0e-5 # global accuracy over sampled grid points - indices = CartesianIndices(grid) + indices = nodeindices(grid) max_err = maximum(1:10:length(indices)) do k - abs(sdf(grid[indices[k]]) - exact_sdf(grid[indices[k]])) + abs(sdf(getnode(grid, indices[k])) - exact_sdf(getnode(grid, indices[k]))) end @test max_err < 1.0e-5 @@ -38,52 +38,44 @@ end grid = LSM.CartesianGrid((-1.0, -1.0, -1.0), (1.0, 1.0, 1.0), (25, 25, 25)) r = 0.45 exact_sdf(x) = norm(x) - r - ϕ = LSM.LevelSet(exact_sdf, grid) + ϕ = LSM.MeshField(exact_sdf, grid) sdf = NewtonSDF(ϕ; upsample = 3) @test sdf(SVector(r, 0.0, 0.0)) ≈ 0.0 atol = 1.0e-4 @test sdf(SVector(0.0, 0.0, 0.0)) ≈ -r atol = 1.0e-4 - indices = CartesianIndices(grid) + indices = nodeindices(grid) max_err = maximum(1:20:length(indices)) do k - abs(sdf(grid[indices[k]]) - exact_sdf(grid[indices[k]])) + abs(sdf(getnode(grid, indices[k])) - exact_sdf(getnode(grid, indices[k]))) end @test max_err < 5.0e-3 @test check_allocs(sdf, SVector(0.25, 0.25, 0.25)) == 0 end - @testset "from PiecewisePolynomialInterpolant" begin - grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (30, 30)) - ϕ = LSM.LevelSet(x -> norm(x) - 0.5, grid) - itp = interpolate(ϕ, 3) - sdf = NewtonSDF(itp; upsample = 4) - @test sdf(SVector(0.5, 0.0)) ≈ 0.0 atol = 2.0e-5 - end - @testset "get_sample_points" begin grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (20, 20)) - ϕ = LSM.LevelSet(x -> norm(x) - 0.5, grid) + ϕ = LSM.MeshField(x -> norm(x) - 0.5, grid) sdf = NewtonSDF(ϕ; upsample = 3) pts = get_sample_points(sdf) @test length(pts) > 0 - @test maximum(p -> abs(sdf.itp(p)), pts) < 1.0e-6 + @test maximum(p -> abs(sdf.meshfield(p)), pts) < 1.0e-6 end @testset "update!" begin grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (40, 40)) r1, r2 = 0.5, 0.3 - sdf = NewtonSDF(LSM.LevelSet(x -> norm(x) - r1, grid); upsample = 4) + sdf = NewtonSDF(LSM.MeshField(x -> norm(x) - r1, grid); upsample = 4) @test sdf(SVector(r1, 0.0)) ≈ 0.0 atol = 2.0e-4 - update!(sdf, LSM.LevelSet(x -> norm(x) - r2, grid)) + update!(sdf, LSM.MeshField(x -> norm(x) - r2, grid)) @test sdf(SVector(r2, 0.0)) ≈ 0.0 atol = 2.0e-4 @test sdf(SVector(r1, 0.0)) ≈ r1 - r2 atol = 1.0e-4 end @testset "deepcopy" begin grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (30, 30)) - sdf = NewtonSDF(LSM.LevelSet(x -> norm(x) - 0.5, grid); upsample = 4) + sdf = NewtonSDF(LSM.MeshField(x -> norm(x) - 0.5, grid); upsample = 4) sdf2 = deepcopy(sdf) @test sdf2(SVector(0.5, 0.0)) ≈ 0.0 atol = 2.0e-5 @test sdf2(SVector(0.0, 0.0)) ≈ -0.5 atol = 2.0e-5 @@ -94,14 +86,14 @@ end @testset "2D circle" begin grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (100, 100)) exact_sdf(x) = sqrt(x[1]^2 + x[2]^2) - 0.5 - ϕ = LSM.LevelSet(x -> (x[1]^2 + x[2]^2) - 0.25, grid) + ϕ = LSM.MeshField(x -> (x[1]^2 + x[2]^2) - 0.25, grid) @test abs(LSM.volume(ϕ) - π / 4) < 1.0e-2 LSM.reinitialize!(ϕ, LSM.NewtonReinitializer()) - max_err = maximum(eachindex(grid)) do i - abs(ϕ[i] - exact_sdf(grid[i])) + max_err = maximum(nodeindices(grid)) do i + abs(ϕ[i] - exact_sdf(getnode(grid, i))) end @test max_err < 1.0e-8 @test abs(LSM.volume(ϕ) - π / 4) < 1.0e-2 # volume preserved @@ -110,12 +102,12 @@ end @testset "3D sphere" begin grid = LSM.CartesianGrid((-1.0, -1.0, -1.0), (1.0, 1.0, 1.0), (31, 31, 31)) exact_sdf(x) = sqrt(x[1]^2 + x[2]^2 + x[3]^2) - 0.45 - ϕ = LSM.LevelSet(x -> (x[1]^2 + x[2]^2 + x[3]^2) - 0.45^2, grid) + ϕ = LSM.MeshField(x -> (x[1]^2 + x[2]^2 + x[3]^2) - 0.45^2, grid) LSM.reinitialize!(ϕ, LSM.NewtonReinitializer(; upsample = 4)) - max_err = maximum(eachindex(grid)) do i - abs(ϕ[i] - exact_sdf(grid[i])) + max_err = maximum(nodeindices(grid)) do i + abs(ϕ[i] - exact_sdf(getnode(grid, i))) end @test max_err < 5.0e-3 end diff --git a/test/test-semi-implicit.jl b/test/test-semi-implicit.jl index 2d5f0d8..48f21b7 100644 --- a/test/test-semi-implicit.jl +++ b/test/test-semi-implicit.jl @@ -4,7 +4,7 @@ using Test @testset "SemiImplicitI2OE periodic transport" begin grid = CartesianGrid((0.0,), (1.0,), (201,)) - ϕ0 = LevelSet(grid) do (x,) + ϕ0 = MeshField(grid) do (x,) sin(2π * x) + 0.15 * cos(6π * x) end vel = MeshField(x -> SVector(1.0), grid) @@ -18,7 +18,7 @@ using Test tf = 0.35 integrate!(eq, tf) - ϕ_ref = LevelSet(grid) do (x,) + ϕ_ref = MeshField(grid) do (x,) xshift = mod(x - tf, 1.0) sin(2π * xshift) + 0.15 * cos(6π * xshift) end @@ -28,7 +28,7 @@ end @testset "SemiImplicitI2OE periodic transport 2D" begin grid = CartesianGrid((0.0, 0.0), (1.0, 1.0), (121, 111)) - ϕ0 = LevelSet(grid) do (x, y) + ϕ0 = MeshField(grid) do (x, y) sin(2π * x) + 0.4 * cos(2π * y) end vel = MeshField(x -> SVector(0.75, -0.35), grid) @@ -41,7 +41,7 @@ end tf = 0.2 integrate!(eq, tf) - ϕ_ref = LevelSet(grid) do (x, y) + ϕ_ref = MeshField(grid) do (x, y) xshift = mod(x - 0.75 * tf, 1.0) yshift = mod(y + 0.35 * tf, 1.0) sin(2π * xshift) + 0.4 * cos(2π * yshift) @@ -52,7 +52,7 @@ end @testset "SemiImplicitI2OE supports non-periodic BCs" begin grid = CartesianGrid((0.0,), (1.0,), (121,)) - ϕ0 = LevelSet(x -> 0.7, grid) + ϕ0 = MeshField(x -> 0.7, grid) term_neumann = AdvectionTerm((x, t) -> sin(2π * x[1]), Upwind()) eq_neumann = LevelSetEquation(; terms = (term_neumann,), @@ -67,7 +67,7 @@ end eq_dirichlet = LevelSetEquation(; terms = (term_dirichlet,), integrator = SemiImplicitI2OE(cfl = 2.0), - ic = LevelSet(x -> 0.0, grid), + ic = MeshField(x -> 0.0, grid), bc = DirichletBC((x, t) -> 0.0), ) integrate!(eq_dirichlet, 0.4) @@ -78,7 +78,7 @@ end @testset "SemiImplicitI2OE checks invalid setup" begin grid1d = CartesianGrid((0.0,), (1.0,), (41,)) - ϕ1d = LevelSet(x -> x[1], grid1d) + ϕ1d = MeshField(x -> x[1], grid1d) eq_multiterm = LevelSetEquation(; terms = ( AdvectionTerm((x, t) -> 1.0, Upwind()), @@ -91,7 +91,7 @@ end @test_throws ArgumentError integrate!(eq_multiterm, 0.1) grid_small = CartesianGrid((0.0,), (1.0,), (2,)) - ϕ_small = LevelSet(x -> x[1], grid_small) + ϕ_small = MeshField(x -> x[1], grid_small) eq_small = LevelSetEquation(; terms = (AdvectionTerm((x, t) -> 1.0, Upwind()),), integrator = SemiImplicitI2OE(), @@ -103,7 +103,7 @@ end @testset "SemiImplicitI2OE tolerates larger timesteps than explicit FE" begin grid = CartesianGrid((0.0,), (1.0,), (401,)) - ϕ0 = LevelSet(grid) do (x,) + ϕ0 = MeshField(grid) do (x,) sin(2π * x) + 0.2 * cos(4π * x) end vel = MeshField(x -> SVector(1.0), grid) @@ -126,7 +126,7 @@ end integrate!(eq_semi, tf) integrate!(eq_explicit, tf) - ϕ_ref = LevelSet(grid) do (x,) + ϕ_ref = MeshField(grid) do (x,) xshift = mod(x - tf, 1.0) sin(2π * xshift) + 0.2 * cos(4π * xshift) end @@ -142,7 +142,7 @@ end @testset "SemiImplicitI2OE outperforms explicit FE at high CFL in 2D" begin grid = CartesianGrid((0.0, 0.0), (1.0, 1.0), (121, 121)) - ϕ0 = LevelSet(grid) do (x, y) + ϕ0 = MeshField(grid) do (x, y) sin(2π * x) + 0.25 * cos(4π * y) end vel = MeshField(x -> SVector(0.9, -0.55), grid) @@ -165,7 +165,7 @@ end integrate!(eq_semi, tf) integrate!(eq_explicit, tf) - ϕ_ref = LevelSet(grid) do (x, y) + ϕ_ref = MeshField(grid) do (x, y) xshift = mod(x - 0.9 * tf, 1.0) yshift = mod(y + 0.55 * tf, 1.0) sin(2π * xshift) + 0.25 * cos(4π * yshift) diff --git a/test/test-show.jl b/test/test-show.jl index 1adb923..3542179 100644 --- a/test/test-show.jl +++ b/test/test-show.jl @@ -38,7 +38,7 @@ end end @testset "scalar, with bc" begin - ϕ = LevelSet(x -> x[1]^2 + x[2]^2 - 0.5^2, grid) + ϕ = MeshField(x -> x[1]^2 + x[2]^2 - 0.5^2, grid) eq = LevelSetEquation(; terms = (NormalMotionTerm(1.0),), ic = ϕ, bc = NeumannBC()) s = showstr(current_state(eq)) @test occursin("├─ bc: Degree 0 extrapolation (all)", s) @@ -54,6 +54,17 @@ end end end +@testset "NarrowBandMeshField" begin + grid = CartesianGrid((-1, -1), (1, 1), (20, 20)) + ϕ = MeshField(x -> norm(x) - 0.5, grid) + nb = NarrowBandMeshField(ϕ, 0.3; reinitialize = false) + s = showstr(nb) + @test startswith(s, "NarrowBandMeshField on CartesianGrid in ℝ²") + @test occursin("├─ active:", s) + @test occursin("halfwidth", s) + @test occursin("└─ values:", s) +end + @testset "Time integrators" begin @test showstr(ForwardEuler()) == "ForwardEuler (1st order explicit)\n └─ cfl: 0.5" @test showstr(RK2()) == "RK2 (2nd order TVD Runge-Kutta, Heun's method)\n └─ cfl: 0.5" @@ -79,7 +90,7 @@ end @testset "LevelSetEquation" begin grid = CartesianGrid((-1, -1), (1, 1), (20, 20)) - ϕ = LevelSet(x -> x[1]^2 + x[2]^2 - 0.5^2, grid) + ϕ = MeshField(x -> x[1]^2 + x[2]^2 - 0.5^2, grid) 𝐮 = MeshField(x -> SVector(1.0, 0.0), grid) eq = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = ϕ, bc = NeumannBC()) @@ -102,7 +113,7 @@ end @testset "SimulationLog" begin grid = CartesianGrid((-1, -1), (1, 1), (20, 20)) - ϕ = LevelSet(x -> x[1]^2 + x[2]^2 - 0.5^2, grid) + ϕ = MeshField(x -> x[1]^2 + x[2]^2 - 0.5^2, grid) 𝐮 = MeshField(x -> SVector(1.0, 0.0), grid) eq = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = ϕ, bc = NeumannBC()) @@ -125,7 +136,7 @@ end end @testset "with reinit" begin - ϕ2 = LevelSet(x -> x[1]^2 + x[2]^2 - 0.5^2, grid) + ϕ2 = MeshField(x -> x[1]^2 + x[2]^2 - 0.5^2, grid) eq2 = LevelSetEquation(; terms = (AdvectionTerm(𝐮),), ic = ϕ2, diff --git a/test/test-spiral.jl b/test/test-spiral.jl index a82e347..3e532cb 100644 --- a/test/test-spiral.jl +++ b/test/test-spiral.jl @@ -10,7 +10,7 @@ import LevelSetMethods as LSM grid = LSM.CartesianGrid((-1.0, -1.0), (1.0, 1.0), (50, 50)) d = 1; r0 = 0.5; θ0 = -π / 3; α = π / 100.0 R = [cos(α) -sin(α); sin(α) cos(α)]; M = R * [1 / 0.06^2 0; 0 1 / (4π^2)] * R' - ϕ = LSM.LevelSet(grid) do (x, y) + ϕ = LSM.MeshField(grid) do (x, y) r = sqrt(x^2 + y^2); θ = atan(y, x); res = 1.0e30 for i in 0:4 θ1 = θ + (2i - 4) * π; v = [r - r0; θ1 - θ0] @@ -24,7 +24,7 @@ import LevelSetMethods as LSM ic = deepcopy(ϕ), terms = (CurvatureTerm(b),), bc = ExtrapolationBC(2), reinit = LSM.NewtonReinitializer(; reinit_freq = 1), ) - nb = NarrowBandLevelSet(deepcopy(ϕ); nlayers = 3, reinitialize = true) + nb = NarrowBandMeshField(deepcopy(ϕ); nlayers = 3, reinitialize = true) eq_nb = LevelSetEquation(; ic = nb, terms = (CurvatureTerm(b),), bc = ExtrapolationBC(2), reinit = LSM.NewtonReinitializer(; reinit_freq = 1), @@ -36,7 +36,7 @@ import LevelSetMethods as LSM ϕ_full = current_state(eq_full) ϕ_nb = current_state(eq_nb) - max_err = maximum(LSM.active_indices(ϕ_nb)) do I + max_err = maximum(LSM.nodeindices(ϕ_nb)) do I abs(ϕ_nb[I] - ϕ_full[I]) end @test max_err < 0.05 diff --git a/test/test-timestepping.jl b/test/test-timestepping.jl index 50ea5f6..07db288 100644 --- a/test/test-timestepping.jl +++ b/test/test-timestepping.jl @@ -7,7 +7,7 @@ import LevelSetMethods as LSM # at the final time. WENO5 (5th order spatial) is used so that temporal error dominates. function _advection_error_1d(integrator, N; u = 1.0, tf = 0.5) grid = LSM.CartesianGrid((-1.0,), (1.0,), (N,)) - ϕ = LSM.LevelSet(x -> sin(π * x[1]), grid) + ϕ = LSM.MeshField(x -> sin(π * x[1]), grid) eq = LSM.LevelSetEquation(; terms = (LSM.AdvectionTerm((x, t) -> SVector(u)),), ic = ϕ, @@ -16,8 +16,8 @@ function _advection_error_1d(integrator, N; u = 1.0, tf = 0.5) ) integrate!(eq, tf) ϕ_out = LSM.current_state(eq) - return maximum(CartesianIndices(LSM.mesh(ϕ_out))) do I - abs(ϕ_out[I] - sin(π * (grid[I][1] - u * tf))) + return maximum(nodeindices(LSM.mesh(ϕ_out))) do I + abs(ϕ_out[I] - sin(π * (getnode(grid, I)[1] - u * tf))) end end diff --git a/test/test-velocityextension.jl b/test/test-velocityextension.jl index 327d791..573035c 100644 --- a/test/test-velocityextension.jl +++ b/test/test-velocityextension.jl @@ -3,13 +3,13 @@ using LevelSetMethods @testset "Normal Motion update hook" begin grid = CartesianGrid((-1.0, -1.0), (1.0, 1.0), (21, 21)) - ϕ = LevelSet(grid) do (x, y) + ϕ = MeshField(grid) do (x, y) sqrt(x^2 + y^2) - 0.5 end bcs = ((PeriodicBC(), PeriodicBC()), (PeriodicBC(), PeriodicBC())) - ϕ = LevelSetMethods.add_boundary_conditions(ϕ, bcs) + ϕ = LevelSetMethods._add_boundary_conditions(ϕ, bcs) - v = MeshField(zeros(Float64, size(grid)...), grid, nothing) + v = MeshField(zeros(Float64, size(grid)...), grid; bc = nothing) term = NormalMotionTerm(v, (speed, ϕ, t) -> (fill!(values(speed), 2t); nothing)) LevelSetMethods.update_term!(term, ϕ, 0.3) @@ -18,16 +18,16 @@ end @testset "Extend Along Normals" begin grid = CartesianGrid((-1.0, -1.0), (1.0, 1.0), (81, 61)) - ϕ = LevelSet(grid) do (x, y) + ϕ = MeshField(grid) do (x, y) x end F = zeros(Float64, size(grid)...) frozen = falses(size(F)) Δ = minimum(LevelSetMethods.meshsize(grid)) - for I in CartesianIndices(F) + for I in nodeindices(grid) if abs(ϕ[I]) <= Δ - y = grid[I][2] + y = getnode(grid, I)[2] F[I] = sin(π * y) frozen[I] = true end @@ -45,19 +45,19 @@ end @testset "Circle periodic extension" begin grid = CartesianGrid((-1.0, -1.0), (1.0, 1.0), (121, 121)) R = 0.55 - ϕ = LevelSet(grid) do (x, y) + ϕ = MeshField(grid) do (x, y) sqrt(x^2 + y^2) - R end bcs = ((PeriodicBC(), PeriodicBC()), (PeriodicBC(), PeriodicBC())) - ϕ = LevelSetMethods.add_boundary_conditions(ϕ, bcs) + ϕ = LevelSetMethods._add_boundary_conditions(ϕ, bcs) v = zeros(Float64, size(grid)...) Δ = minimum(LevelSetMethods.meshsize(grid)) seed_band = 1.1 frozen = abs.(values(ϕ)) .<= seed_band * Δ - for I in CartesianIndices(v) + for I in nodeindices(grid) frozen[I] || continue - x, y = grid[I] + x, y = getnode(grid, I) r = sqrt(x^2 + y^2) v[I] = y / max(r, eps(Float64)) end @@ -66,7 +66,7 @@ end extend_along_normals!(v, ϕ; nb_iters = 100, frozen, cfl = 0.45) @test v[frozen] == v_seed[frozen] - vf = MeshField(v, grid, bcs) + vf = MeshField(v, grid; bc = bcs) extend_band = 5.0 sum_abs_n_dot_grad = 0.0 nb_samples = 0 @@ -86,11 +86,11 @@ end @testset "MeshField and argument checks" begin grid = CartesianGrid((-1.0, -1.0), (1.0, 1.0), (41, 41)) - ϕ = LevelSet(grid) do (x, y) + ϕ = MeshField(grid) do (x, y) x + y end vals = zeros(Float64, size(grid)...) - F = MeshField(vals, grid, nothing) + F = MeshField(vals, grid; bc = nothing) out = extend_along_normals!(F, ϕ; nb_iters = 5) @test out === F @@ -113,7 +113,7 @@ function _run_curvature_extension_cycle!( ) grid = LevelSetMethods.mesh(ϕ) speed_vals = zeros(Float64, size(grid)...) - speed = MeshField(speed_vals, grid, nothing) + speed = MeshField(speed_vals, grid; bc = nothing) update_speed = function (v, ϕstate, t) vals = values(v) @@ -156,7 +156,7 @@ function _interface_radius_stats(ϕ; band = 1.5) radii = Float64[] for I in eachindex(ϕ) if abs(ϕ[I]) <= band * Δ - x = grid[I] + x = getnode(grid, I) push!(radii, sqrt(sum(abs2, x))) end end @@ -169,7 +169,7 @@ end @testset "Classical circular reconstruction (2D)" begin grid = CartesianGrid((-0.5, -0.5), (0.5, 0.5), (128, 128)) R0 = 0.45 # diameter = 9/10 - ϕ = LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - R0, grid) + ϕ = MeshField(x -> sqrt(x[1]^2 + x[2]^2) - R0, grid) Δ = minimum(LevelSetMethods.meshsize(grid)) ϕf = _run_curvature_extension_cycle!( @@ -189,7 +189,7 @@ end @testset "Classical spherical reconstruction (3D)" begin grid = CartesianGrid((-0.5, -0.5, -0.5), (0.5, 0.5, 0.5), (48, 48, 48)) R0 = 0.45 - ϕ = LevelSet(x -> sqrt(x[1]^2 + x[2]^2 + x[3]^2) - R0, grid) + ϕ = MeshField(x -> sqrt(x[1]^2 + x[2]^2 + x[3]^2) - R0, grid) Δ = minimum(LevelSetMethods.meshsize(grid)) ϕf = _run_curvature_extension_cycle!( @@ -225,7 +225,7 @@ end ϕ = LevelSetMethods.star(grid; radius = R, deformation = deformation, n = nfacets) bcs = ((PeriodicBC(), PeriodicBC()), (PeriodicBC(), PeriodicBC())) - ϕb = LevelSetMethods.add_boundary_conditions(ϕ, bcs) + ϕb = LevelSetMethods._add_boundary_conditions(ϕ, bcs) Δ = minimum(LevelSetMethods.meshsize(grid)) v = zeros(Float64, size(grid)...) @@ -262,7 +262,7 @@ end δ = minimum(LevelSetMethods.meshsize(ϕstate)) for I in eachindex(ϕstate) if abs(ϕstate[I]) <= 1.5 * δ - x = g[I] + x = getnode(g, I) push!(rs, sqrt(x[1]^2 + x[2]^2)) end end @@ -271,7 +271,7 @@ end return std_rs / mean_rs end cv0 = _radius_cv(ϕb) - term = NormalMotionTerm(MeshField(v, grid, nothing)) + term = NormalMotionTerm(MeshField(v, grid; bc = nothing)) eq = LevelSetEquation(; terms = (term,), ic = ϕ,