Merge pull request #50 from maltezfaria/interpolations-ext

Extension for `Interpolants`

Author: maltezfaria

Project.toml

@@ -8,16 +8,19 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [weakdeps]
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 MMG_jll = "86086c02-e288-5929-a127-40944b0018b7"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 MarchingCubes = "299715c1-40a9-479a-aaf9-4a633d36f717"
 
 [extensions]
+InterpolationsExt = "Interpolations"
 MMGSurfaceExt = ["MMG_jll", "MarchingCubes"]
 MMGVolumeExt = "MMG_jll"
 MakieExt = "Makie"
 
 [compat]
+Interpolations = "0.15"
 LinearAlgebra = "1"
 MMG_jll = "5"
 Makie = "0.21, 0.22"

docs/Project.toml

@@ -3,6 +3,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LevelSetMethods = "ca3df947-0575-4ae6-b9a3-78df393a7451"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 MMG_jll = "86086c02-e288-5929-a127-40944b0018b7"

docs/make.jl

@@ -3,6 +3,7 @@ using Documenter
 using GLMakie
 using MMG_jll
 using MarchingCubes
+using Interpolations
 using StaticArrays
 using DocumenterCitations
 
@@ -22,7 +23,7 @@ const numbered_pages = [
 ]
 
 modules = [LevelSetMethods]
-for extension in [:MakieExt, :MMGSurfaceExt, :MMGVolumeExt]
+for extension in [:MakieExt, :MMGSurfaceExt, :MMGVolumeExt, :InterpolationsExt]
     ext = Base.get_extension(LevelSetMethods, extension)
     isnothing(ext) && "error loading $ext"
     push!(modules, ext)
@@ -42,7 +43,8 @@ makedocs(;
         "terms.md",
         "time-integrators.md",
         "boundary-conditions.md",
-        "Extensions" => ["extension-makie.md", "extension-mmg.md"],
+        "Extensions" =>
+            ["extension-makie.md", "extension-mmg.md", "extension-interpolations.md"],
         "Examples" => ["example-zalesak.md", "example-shape-optim.md"],
         numbered_pages,
     ),

docs/src/extension-interpolations.md

@@ -0,0 +1,56 @@
+# [Interpolations extension](@id extension-interpolations)
+
+This extension overloads the `interpolate` function from
+[`Interpolations.jl`](https://juliamath.github.io/Interpolations.jl/latest/) to provide a
+way to construct a global interpolant from the discrete data in a
+[`LevelSet`](@ref) or [`LevelSetEquation`](@ref). This can be useful in situations where you want
+to evaluate the approximate underlying functions at points that are not on the grid.
+
+Here is an example of how to construct such an interpolant:
+
+```@example interpolations
+using LevelSetMethods, Interpolations
+LevelSetMethods.set_makie_theme!()
+a, b = (-2, -2), (2, 2)
+ϕ   = LevelSetMethods.star(CartesianGrid(a, b, (50, 50)))
+itp = interpolate(ϕ, BSpline(Cubic())) # create the interpolant
+```
+
+Once constructed, the interpolant can be used to evaluate the level-set function anywhere
+inside the grid:
+
+```@example interpolations
+itp(0.5, 0.5)
+```
+
+This can be used e.g. to plot the level-set function using `Makie`:
+
+```@example interpolations
+using CairoMakie
+xx = yy = -2:0.01:2
+contour(xx, yy, [itp(x,y) for x in xx, y in yy]; levels = [0], linewidth = 2)
+```
+
+Trying to evaluate it outside the domain will throw an error:
+
+```@example interpolations
+try
+  itp(3, 0.1)
+catch e
+    println("Error caught")
+end
+```
+
+Using it on three-dimensional level sets is similar:
+
+```@example interpolations
+using 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
+f = (x) -> norm(x .- P1)*norm(x .- P2) - b^2
+ϕ = LevelSet(f, grid)
+itp = interpolate(ϕ) # cubic spline by default
+println("ϕ(0.5, 0.5, 0.5)   = ", f((0.5, 0.5, 0.5)))
+println("itp(0.5, 0.5, 0.5) = ", itp(0.5, 0.5, 0.5))
+```

docs/src/extension-makie.md

@@ -1,8 +1,8 @@
 # [Makie extension](@id extension-makie)
 
-[Makie](https://docs.makie.org/v0.21/) can be used to visualize levelset functions in both 2
-and 3 dimensions. After loading one of the `Makie` backends (we recomment `GLMakie` for 3D),
-you can simply call `plot` on the levelset function to visualize it. For example:
+[Makie](https://docs.makie.org/v0.21/) can be used to visualize level set functions in both 2
+and 3 dimensions. After loading one of the `Makie` backends (we recommend `GLMakie` for 3D),
+you can simply call `plot` on the level set function to visualize it. For example:
 
 ```@example contour2D
 using LevelSetMethods, GLMakie
@@ -11,15 +11,15 @@ grid = CartesianGrid((-2, -2), (2, 2), (100, 100))
 plot(ϕ)
 ```
 
-By default, only the zero level-set is plotted as a contour line. For more control, simply
+By default, only the zero level set is plotted as a contour line. For more control, simply
 call the `contour` (or `contourf`) function from `Makie` directly. For example:
 
 ```@example contour2D
 contour(ϕ; levels = [-0.5, 0, 0.5], labels = true)
 ```
 
 Although you can manually customize the `Axis` attributes for the plot, `LevelSetMethods`
-provides a `Theme` with some reasonable defaults for plotting levelset functions:
+provides a `Theme` with some reasonable defaults for plotting level set functions:
 
 ```@example contour2D
 theme = LevelSetMethods.makie_theme()
@@ -28,7 +28,7 @@ with_theme(theme) do
 end
 ```
 
-In `3D`, the `plot` function will plot the zero level-set as an isosurface. For example:
+In `3D`, the `plot` function will plot the zero level set as an isosurface. For example:
 
 ```@example volume3D
 using LevelSetMethods, GLMakie, LinearAlgebra

docs/src/extension-mmg.md

@@ -1,7 +1,7 @@
 # [MMG extension](@id extension-mmg)
 
-This extension provides functions to generate meshes of levelset functions using [MMG](https://www.mmgtools.org/).
-It define two methods: `export_volume_mesh` and `export_surface_mesh`.
+This extension provides functions to generate meshes of level-set functions using [MMG](https://www.mmgtools.org/).
+It defines two methods: `export_volume_mesh` and `export_surface_mesh`.
 For both of them, it is possible to control the size of the generated mesh using the following optional parameters:
 
 - `hgrad` control the growth ratio between two adjacent edges.
@@ -10,7 +10,7 @@ For both of them, it is possible to control the size of the generated mesh using
 
 ## Generation of 2D and 3D mesh from a level-set
 
-For 2 and 3 dimensional Cartesian levelset, one can use the `export_volume_mesh` function to generate meshes.
+For 2 and 3-dimensional Cartesian level set, one can use the `export_volume_mesh` function to generate meshes.
 This method relies on the `mmg2d_O3` and `mmg3d_O3` utilities.
 Example in 2D:
 

ext/InterpolationsExt.jl

@@ -0,0 +1,22 @@
+module InterpolationsExt
+
+import Interpolations as Itp
+import LevelSetMethods as LSM
+
+function __init__()
+    @info "Loading Interpolations extension for LevelSetMethods.jl"
+end
+
+Itp.interpolate(ϕ::LSM.LevelSet) = Itp.interpolate(ϕ, Itp.BSpline(Itp.Cubic()))
+
+function Itp.interpolate(ϕ::LSM.LevelSet, mode)
+    itp = Itp.interpolate(ϕ.vals, mode)
+    grids = LSM.grid1d(LSM.mesh(ϕ))
+    return Itp.scale(itp, grids...)
+end
+
+function Itp.interpolate(eq::LSM.LevelSetEquation, args...; kwargs...)
+    return Itp.interpolate(LSM.current_state(eq), args...; kwargs...)
+end
+
+end # module

ext/MMGVolumeExt.jl

@@ -78,7 +78,7 @@ function LSM.export_volume_mesh(
                 num_triangles = 2(nx - 1) * (ny - 1)
                 write(file, "\nTriangles\n")
                 write(file, "$num_triangles\n")
-                for x_id in 1:nx-1, y_id in 1:ny-1
+                for x_id in 1:(nx-1), y_id in 1:(ny-1)
                     c00 = (y_id - 1) * nx + x_id
                     c10 = c00 + 1
                     c01 = c00 + nx
@@ -96,7 +96,7 @@ function LSM.export_volume_mesh(
                 num_tetrahedrons = 6(nx - 1) * (ny - 1) * (nz - 1)
                 write(file, "\nTetrahedra\n")
                 write(file, "$num_tetrahedrons\n")
-                for x_id in 1:nx-1, y_id in 1:ny-1, z_id in 1:nz-1
+                for x_id in 1:(nx-1), y_id in 1:(ny-1), z_id in 1:(nz-1)
                     c000 = (z_id - 1) * nx * ny + (y_id - 1) * nx + x_id
                     c100 = c000 + 1
                     c010 = c000 + nx

src/levelsetequation.jl

@@ -102,7 +102,7 @@ function Base.show(io::IO, eq::LevelSetEquation)
     print(io, "Level-set equation given by\n")
     print(io, "\n \t ϕₜ + ")
     terms = eq.terms
-    for term in terms[1:end-1]
+    for term in terms[1:(end-1)]
         print(io, term)
         print(io, " + ")
     end
@@ -138,7 +138,7 @@ function integrate!(ls::LevelSetEquation, tf, Δt = Inf)
            the level-set equation cannot be solved back in time"
     @assert tf >= tc msg
     # append boundary conditions for integration
-    ϕ          = current_state(ls)
+    ϕ         = current_state(ls)
     buf        = buffers(ls)
     integrator = time_integrator(ls)
     # dynamic dispatch. Should not be a problem provided enough computation is
@@ -153,7 +153,7 @@ end
 number_of_buffers(fe::ForwardEuler) = 1
 
 @noinline function _integrate!(ϕ, buffers, integrator::ForwardEuler, terms, tc, tf, Δt)
-    buffer = buffers[1]
+    buffer  = buffers[1]
     α      = cfl(integrator)
     Δt_cfl = α * compute_cfl(terms, ϕ, tc)
     Δt     = min(Δt, Δt_cfl)

src/meshes.jl

@@ -79,7 +79,7 @@ function interior_indices(g::CartesianGrid, P)
     N  = dimension(g)
     sz = size(g)
     I  = ntuple(N) do dim
-        return P+1:sz[dim]-P
+        return (P+1):(sz[dim]-P)
     end
     return CartesianIndices(I)
 end