Update to DelaunayTriangulation 1.0 (#53)

Project.toml

@@ -16,7 +16,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 [compat]
 ChunkSplitters = "0.1, 1.0, 2.0"
 CommonSolve = "0.2"
-DelaunayTriangulation = "0.7, 0.8"
+DelaunayTriangulation = "1.0"
 PreallocationTools = "0.4"
 PrecompileTools = "1.2"
 SciMLBase = "1, 2"

README.md

@@ -26,19 +26,19 @@ If this package doesn't suit what you need, you may like to review some of the o
  As a very quick demonstration, here is how we could solve a diffusion equation with Dirichlet boundary conditions on a square domain using the standard `FVMProblem` formulation; please see the docs for more information.
 
 ```julia
-using FiniteVolumeMethod, DelaunayTriangulation, CairoMakie, DifferentialEquations
+using FiniteVolumeMethod, DelaunayTriangulation, CairoMakie, OrdinaryDiffEq
 a, b, c, d = 0.0, 2.0, 0.0, 2.0
 nx, ny = 50, 50
 tri = triangulate_rectangle(a, b, c, d, nx, ny, single_boundary=true)
 mesh = FVMGeometry(tri)
 bc = (x, y, t, u, p) -> zero(u)
 BCs = BoundaryConditions(mesh, bc, Dirichlet)
 f = (x, y) -> y ≤ 1.0 ? 50.0 : 0.0
-initial_condition = [f(x, y) for (x, y) in each_point(tri)]
+initial_condition = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
 D = (x, y, t, u, p) -> 1 / 9
 final_time = 0.5
 prob = FVMProblem(mesh, BCs; diffusion_function=D, initial_condition, final_time)
-sol = solve(prob, saveat=0.001)
+sol = solve(prob, Tsit5(), saveat=0.001)
 u = Observable(sol.u[1])
 fig, ax, sc = tricontourf(tri, u, levels=0:5:50, colormap=:matter)
 tightlimits!(ax)

docs/liveserver.jl

@@ -6,6 +6,11 @@ import LiveServer
 withenv("LIVESERVER_ACTIVE" => "true") do
     LiveServer.servedocs(;
         launch_browser=true,
+        foldername=joinpath(repo_root, "docs"),
+        include_dirs=[joinpath(repo_root, "src")],
+        skip_dirs=[joinpath(repo_root, "docs/src/tutorials"),
+            joinpath(repo_root, "docs/src/wyos"),
+            joinpath(repo_root, "docs/src/figures"),
+        ],
     )
-end
-
+end

docs/src/interface.md

@@ -37,7 +37,7 @@ TriangleProperties
 
 # `BoundaryConditions`: Defining boundary conditions
 
-Once a mesh is defined, you need to associate each part of the boundary with a set of boundary nodes. Since you have a `Triangulation`, the boundary of the mesh already meets the necessary assumptions made by this package about the boundary; these assumptions are simply that they match the specification of a boundary [here in DelaunayTriangulation.jl's docs](https://SciML.github.io/DelaunayTriangulation.jl/dev/boundary_handling/#Boundary-Specification) (for example, the boundary points connect, the boundary is positively oriented, etc.).
+Once a mesh is defined, you need to associate each part of the boundary with a set of boundary nodes. Since you have a `Triangulation`, the boundary of the mesh already meets the necessary assumptions made by this package about the boundary; these assumptions are simply that they match the specification of a boundary [here in DelaunayTriangulation.jl's docs](https://juliageometry.github.io/DelaunayTriangulation.jl/dev/manual/boundaries/) (for example, the boundary points connect, the boundary is positively oriented, etc.).
 
 You can specify boundary conditions using `BoundaryConditions`, whose docstring is below; the fields of `BoundaryConditions` are not public API, only this wrapper is.
 

docs/src/literate_tutorials/diffusion_equation_in_a_wedge_with_mixed_boundary_conditions.jl

@@ -29,28 +29,20 @@ tc = DisplayAs.withcontext(:displaysize => (15, 80), :limit => true); #hide
 # one part for each boundary condition. This is accomplished 
 # by providing a single vector for each part of the boundary as follows
 # (and as described in DelaunayTriangulation.jl's documentation),
-# where we also `refine!` the mesh to get a better mesh. 
+# where we also `refine!` the mesh to get a better mesh. For the arc, 
+# we use the `CircularArc` so that the mesh knows that it is triangulating 
+# a certain arc in that area.
 using DelaunayTriangulation, FiniteVolumeMethod, ElasticArrays
 using ReferenceTests, Bessels, FastGaussQuadrature, Cubature #src
-n = 50
+
 α = π / 4
-## The bottom edge 
-x₁ = [0.0, 1.0]
-y₁ = [0.0, 0.0]
-## The arc 
-r₂ = fill(1, n)
-θ₂ = LinRange(0, α, n)
-x₂ = @. r₂ * cos(θ₂)
-y₂ = @. r₂ * sin(θ₂)
-## The upper edge 
-x₃ = [cos(α), 0.0]
-y₃ = [sin(α), 0.0]
-## Now combine and create the mesh 
-x = [x₁, x₂, x₃]
-y = [y₁, y₂, y₃]
-boundary_nodes, points = convert_boundary_points_to_indices(x, y; existing_points=ElasticMatrix{Float64}(undef, 2, 0))
+points = [(0.0, 0.0), (1.0, 0.0), (cos(α), sin(α))]
+bottom_edge = [1, 2]
+arc = CircularArc((1.0, 0.0), (cos(α), sin(α)), (0.0, 0.0))
+upper_edge = [3, 1]
+boundary_nodes = [bottom_edge, [arc], upper_edge]
 tri = triangulate(points; boundary_nodes)
-A = get_total_area(tri)
+A = get_area(tri)
 refine!(tri; max_area=1e-4A)
 mesh = FVMGeometry(tri)
 
@@ -74,7 +66,7 @@ BCs = BoundaryConditions(mesh, (lower_bc, arc_bc, upper_bc), types)
 # specifying the diffusion function as a constant. 
 f = (x, y) -> 1 - sqrt(x^2 + y^2)
 D = (x, y, t, u, p) -> one(u)
-initial_condition = [f(x, y) for (x, y) in each_point(tri)]
+initial_condition = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
 final_time = 0.1
 prob = FVMProblem(mesh, BCs; diffusion_function=D, initial_condition, final_time)
 
@@ -89,12 +81,18 @@ flux = (x, y, t, α, β, γ, p) -> (-α, -β)
 # has the best performance for these problems.
 using OrdinaryDiffEq, LinearSolve
 sol = solve(prob, TRBDF2(linsolve=KLUFactorization()), saveat=0.01, parallel=Val(false))
+ind = findall(DelaunayTriangulation.each_point_index(tri)) do i #hide
+    !DelaunayTriangulation.has_vertex(tri, i) #hide
+end #hide
+using Test #hide
+@test sol[ind, :] ≈ reshape(repeat(initial_condition, length(sol)), :, length(sol))[ind, :] # make sure that missing vertices don't change #hide
 sol |> tc #hide
 
 #-
 using CairoMakie
 fig = Figure(fontsize=38)
 for (i, j) in zip(1:3, (1, 6, 11))
+    local ax
     ax = Axis(fig[1, i], width=600, height=600,
         xlabel="x", ylabel="y",
         title="t = $(sol.t[j])",
@@ -145,7 +143,7 @@ function exact_solution(x, y, t, A, ζ, f, α) #src
     return s #src
 end #src
 function compare_solutions(sol, tri, α, f) #src
-    n = DelaunayTriangulation.num_solid_vertices(tri) #src
+    n = DelaunayTriangulation.num_points(tri) #src
     x = zeros(n, length(sol)) #src
     y = zeros(n, length(sol)) #src
     u = zeros(n, length(sol)) #src
@@ -162,6 +160,7 @@ end #src
 x, y, u = compare_solutions(sol, tri, α, f) #src
 fig = Figure(fontsize=64) #src
 for i in eachindex(sol) #src
+    local ax
     ax = Axis(fig[1, i], width=600, height=600) #src
     tricontourf!(ax, tri, sol.u[i], levels=0:0.01:1, colormap=:matter) #src
     ax = Axis(fig[2, i], width=600, height=600) #src

docs/src/literate_tutorials/diffusion_equation_on_a_square_plate.jl

@@ -36,7 +36,7 @@ BCs = BoundaryConditions(mesh, bc, Dirichlet)
 
 # We can now define the actual PDE. We start by defining the initial condition and the diffusion function. 
 f = (x, y) -> y ≤ 1.0 ? 50.0 : 0.0
-initial_condition = [f(x, y) for (x, y) in each_point(tri)]
+initial_condition = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
 D = (x, y, t, u, p) -> 1 / 9
 
 # We can now define the problem:
@@ -50,17 +50,19 @@ prob.flux_function
 # where $(\alpha, \beta, \gamma)$ defines the approximation to $u$ via $u(x, y) = \alpha x + \beta y + \gamma$ so that 
 # $\grad u(\vb x, t) = (\alpha,\beta)^{\mkern-1.5mu\mathsf{T}}$.
 
-# To now solve the problem, we simply use `solve`. When no algorithm 
-# is provided, as long as DifferentialEquations is loaded (instead of e.g. 
-# OrdinaryDiffEq), the algorithm is chosen automatically. Moreover, note that, 
+# To now solve the problem, we simply use `solve`. Note that, 
 # in the `solve` call below, multithreading is enabled by default.
-using DifferentialEquations
-sol = solve(prob, saveat=0.05)
+# (If you don't know what algorithm to consider, do `using DifferentialEquations` instead 
+# and simply call `solve(prob, saveat=0.05)` so that the algorithm is chosen automatically instead 
+# of using `Tsit5()`.)
+using OrdinaryDiffEq
+sol = solve(prob, Tsit5(), saveat=0.05)
 sol |> tc #hide
 
 # To visualise the solution, we can use `tricontourf!` from Makie.jl. 
 fig = Figure(fontsize=38)
 for (i, j) in zip(1:3, (1, 6, 11))
+    local ax
     ax = Axis(fig[1, i], width=600, height=600,
         xlabel="x", ylabel="y",
         title="t = $(sol.t[j])",
@@ -88,7 +90,7 @@ function exact_solution(x, y, t) #src
     end #src
 end #src
 function compare_solutions(sol, tri) #src
-    n = DelaunayTriangulation.num_solid_vertices(tri) #src
+    n = DelaunayTriangulation.num_points(tri) #src
     x = zeros(n, length(sol)) #src
     y = zeros(n, length(sol)) #src
     u = zeros(n, length(sol)) #src
@@ -103,6 +105,7 @@ end #src
 x, y, u = compare_solutions(sol, tri) #src
 fig = Figure(fontsize=64) #src
 for i in eachindex(sol) #src
+    local ax
     ax = Axis(fig[1, i], width=600, height=600) #src
     tricontourf!(ax, tri, sol.u[i], levels=0:5:50, colormap=:matter) #src
     ax = Axis(fig[2, i], width=600, height=600) #src

docs/src/literate_tutorials/diffusion_equation_on_an_annulus.jl

@@ -20,33 +20,19 @@ tc = DisplayAs.withcontext(:displaysize => (15, 80), :limit => true); #hide
 # ```math
 # u_0(x) = 10\mathrm{e}^{-25\left[\left(x+\frac12\right)^2+\left(y+\frac12\right)^2\right]} - 10\mathrm{e}^{-45\left[\left(x-\frac12\right)^2+\left(y-\frac12\right)^2\right]} - 5\mathrm{e}^{-50\left[\left(x+\frac{3}{10}\right)^2+\left(y+\frac12\right)^2\right]}.
 # ```
-# The complicated task for this problem is the definition 
-# of the mesh of the annulus. We need to follow the boundary 
-# specification from DelaunayTriangulation.jl, discussed 
-# [here](https://SciML.github.io/DelaunayTriangulation.jl/dev/boundary_handling/).
-# In particular, the outer boundary must be counter-clockwise, 
-# the inner boundary be clockwise, and we need to provide 
-# the nodes as a `Vector{Vector{Vector{Int}}}`.
-# We define this mesh below. 
+# For the mesh, we use two `CircularArc`s to define the annulus.
 using DelaunayTriangulation, FiniteVolumeMethod, CairoMakie
-R₁ = 0.2
-R₂ = 1.0
-θ = collect(LinRange(0, 2π, 100))
-θ[end] = 0.0 # get the endpoints to match
-x = [
-    [R₂ .* cos.(θ)], # outer first
-    [reverse(R₁ .* cos.(θ))] # then inner - reverse to get clockwise orientation
-]
-y = [
-    [R₂ .* sin.(θ)], # 
-    [reverse(R₁ .* sin.(θ))]
-]
-boundary_nodes, points = convert_boundary_points_to_indices(x, y)
+R₁, R₂ = 0.2, 1.0
+inner = CircularArc((R₁, 0.0), (R₁, 0.0), (0.0, 0.0), positive=false)
+outer = CircularArc((R₂, 0.0), (R₂, 0.0), (0.0, 0.0))
+boundary_nodes = [[[outer]], [[inner]]]
+points = NTuple{2,Float64}[]
 tri = triangulate(points; boundary_nodes)
-A = get_total_area(tri)
+A = get_area(tri)
 refine!(tri; max_area=1e-4A)
 triplot(tri)
 
+
 #-
 mesh = FVMGeometry(tri)
 
@@ -76,7 +62,7 @@ initial_condition_f = (x, y) -> begin
     10 * exp(-25 * ((x + 0.5) * (x + 0.5) + (y + 0.5) * (y + 0.5))) - 5 * exp(-50 * ((x + 0.3) * (x + 0.3) + (y + 0.5) * (y + 0.5))) - 10 * exp(-45 * ((x - 0.5) * (x - 0.5) + (y - 0.5) * (y - 0.5)))
 end
 diffusion_function = (x, y, t, u, p) -> one(u)
-initial_condition = [initial_condition_f(x, y) for (x, y) in each_point(tri)]
+initial_condition = [initial_condition_f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
 final_time = 2.0
 prob = FVMProblem(mesh, BCs;
     diffusion_function,
@@ -91,6 +77,7 @@ sol |> tc #hide
 #-
 fig = Figure(fontsize=38)
 for (i, j) in zip(1:3, (1, 6, 11))
+    local ax
     ax = Axis(fig[1, i], width=600, height=600,
         xlabel="x", ylabel="y",
         title="t = $(sol.t[j])",
@@ -115,6 +102,9 @@ using ReferenceTests #src
 # apply `jump_and_march`. This is done with `add_ghost_triangles!`.
 add_ghost_triangles!(tri)
 
+# (Actually, `tri` already had these ghost triangles,
+# but we are just showing how you would add them back in if needed.)
+
 # Now let's interpolate.
 x = LinRange(-R₂, R₂, 400)
 y = LinRange(-R₂, R₂, 400)
@@ -124,7 +114,7 @@ last_triangle = Ref((1, 1, 1))
 for (j, _y) in enumerate(y)
     for (i, _x) in enumerate(x)
         T = jump_and_march(tri, (_x, _y), try_points=last_triangle[])
-        last_triangle[] = indices(T) # used to accelerate jump_and_march, since the points we're looking for are close to each other
+        last_triangle[] = triangle_vertices(T) # used to accelerate jump_and_march, since the points we're looking for are close to each other
         if DelaunayTriangulation.is_ghost_triangle(T) # don't extrapolate
             interp_vals[i, j] = NaN
         else
@@ -161,6 +151,6 @@ itp_vals |> tc #hide
 
 #-
 fig, ax, sc = contourf(x, y, reshape(itp_vals, length(x), length(y)), colormap=:matter, levels=-10:2:40)
-fig 
+fig
 tricontourf!(Axis(fig[1, 2]), tri, u, levels=-10:2:40, colormap=:matter) #src
 @test_reference joinpath(@__DIR__, "../figures", "diffusion_equation_on_an_annulus_interpolated_with_naturalneighbours.png") fig #src

docs/src/literate_tutorials/equilibrium_temperature_distribution_with_mixed_boundary_conditions_and_using_ensembleproblems.jl

@@ -51,7 +51,7 @@ bn4 = [F, G]
 bn = [bn1, bn2, bn3, bn4]
 boundary_nodes, points = convert_boundary_points_to_indices(bn)
 tri = triangulate(points; boundary_nodes)
-refine!(tri; max_area=1e-4get_total_area(tri))
+refine!(tri; max_area=1e-4get_area(tri))
 triplot(tri)
 
 #-
@@ -78,7 +78,7 @@ BCs = BoundaryConditions(mesh, (bc1, bc2, bc3, bc4),
 # down to $T=40$ at $y=0$.
 diffusion_function = (x, y, t, T, p) -> one(T)
 f = (x, y) -> 500y + 40
-initial_condition = [f(x, y) for (x, y) in each_point(tri)]
+initial_condition = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
 final_time = Inf
 prob = FVMProblem(mesh, BCs;
     diffusion_function,
@@ -98,70 +98,3 @@ fig, ax, sc = tricontourf(tri, sol.u, levels=40:70, axis=(xlabel="x", ylabel="y"
 fig
 using ReferenceTests #src
 @test_reference joinpath(@__DIR__, "../figures", "equilibrium_temperature_distribution_with_mixed_boundary_conditions_and_using_ensembleproblems.png") fig #src
-
-# Let us now suppose we are interested in how the ambient temperature, $T_{\infty}$, 
-# affects the temperature distribution. In particular, let us ask the following question:
-# > What range of $T_{\infty}$ will allow the temperature at $(0.03, 0.03)$ to be between $50$ and $55$?
-# To answer this question, we use an `EnsembleProblem` so that we can solve the problem over many 
-# values of $T_{\infty}$ efficiently. For these new problems, it would be 
-# a good idea to use a new initial condition given by the solution of the previous problem.
-copyto!(prob.initial_condition, sol.u)
-using Accessors
-T∞_range = LinRange(-100, 100, 101)
-ens_prob = EnsembleProblem(steady_prob,
-    prob_func=(prob, i, repeat) -> let T∞_range = T∞_range, h = h, k = k
-        _prob =
-            @set prob.problem.conditions.functions[3].parameters =
-                (h=h, T∞=T∞_range[i], k=k)
-        return _prob
-    end)
-esol = solve(ens_prob, DynamicSS(Rosenbrock23()), EnsembleSerial(); trajectories=length(T∞_range))
-esol |> tc #hide
-
-# From these results, let us now extract the temperature at $(0.03, 0.03)$. We will use 
-# NaturalNeighbours.jl for this.
-using NaturalNeighbours
-itps = [interpolate(tri, esol[i].u) for i in eachindex(esol)];
-itp_vals = [itp(0.03, 0.03; method=Sibson()) for itp in itps]
-## If you want piecewise linear interpolation, use either method=Triangle()
-## or itp_vals = [pl_interpolate(prob, T, sol.u, 0.03, 0.03) for sol in esol], where 
-## T = jump_and_march(tri, (0.03, 0.03)).
-using Test #src
-_T = jump_and_march(tri, (0.03, 0.03)) #src
-_itp_vals = [pl_interpolate(prob, _T, sol.u, 0.03, 0.03) for sol in esol] #src
-@test _itp_vals ≈ itp_vals rtol = 1e-4 #src
-fig = Figure(fontsize=33)
-ax = Axis(fig[1, 1], xlabel=L"T_{\infty}", ylabel=L"T(0.03, 0.03)")
-lines!(ax, T∞_range, itp_vals, linewidth=4)
-fig
-
-# We see that the temperature at this point seems to increase linearly 
-# with $T_{\infty}$. Let us find precisely where this curve 
-# meets $T=50$ and $T=55$.
-using NonlinearSolve, DataInterpolations
-itp = LinearInterpolation(itp_vals, T∞_range)
-rootf = (u, p) -> p.itp(u) - p.τ[]
-Tthresh = Ref(50.0)
-prob = IntervalNonlinearProblem(rootf, (-100.0, 100.0), (itp=itp, τ=Tthresh))
-sol50 = solve(prob, ITP())
-
-#-
-Tthresh[] = 55.0
-sol55 = solve(prob, ITP())
-
-# So, it seems like the answer to our question is $-11.8 \leq T_{\infty} \leq 55$.
-# Here is an an animation of the temperature distribution as $T_{\infty}$ varies.
-fig = Figure(fontsize=33)
-i = Observable(1)
-tt = map(i -> L"T_{\infty} = %$(rpad(round(T∞_range[i], digits=3),5,'0'))", i)
-u = map(i -> esol.u[i], i)
-ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y",
-    title=tt, titlealign=:left)
-tricontourf!(ax, tri, u, levels=40:70, extendlow=:auto, extendhigh=:auto)
-tightlimits!(ax)
-record(fig, joinpath(@__DIR__, "../figures", "temperature_animation.mp4"), eachindex(esol);
-    framerate=12) do _i
-    i[] = _i
-end;
-
-# ![Animation of the temperature distribution](../figures/temperature_animation.mp4)

docs/src/literate_tutorials/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.jl

@@ -46,8 +46,8 @@ u_Sp = b
 v_Sp = d
 u_icf = (x, y) -> 1 - exp(-80 * (x^2 + y^2))
 v_icf = (x, y) -> exp(-80 * (x^ 2 + y^2))
-u_ic = [u_icf(x, y) for (x, y) in each_point(tri)]
-v_ic = [v_icf(x, y) for (x, y) in each_point(tri)]
+u_ic = [u_icf(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
+v_ic = [v_icf(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
 u_prob = FVMProblem(mesh, u_BCs;
     flux_function=u_q, flux_parameters=u_qp,
     source_function=u_S, source_parameters=u_Sp,
@@ -59,8 +59,8 @@ v_prob = FVMProblem(mesh, v_BCs;
 prob = FVMSystem(u_prob, v_prob)
 
 # Now that we have our system, we can solve.
-using OrdinaryDiffEq, Sundials
-sol = solve(prob, CVODE_BDF(linear_solver=:GMRES), saveat=10.0, parallel=Val(false))
+using OrdinaryDiffEq, LinearSolve
+sol = solve(prob, TRBDF2(linsolve=KLUFactorization()), saveat=10.0, parallel=Val(false))
 sol |> tc #hide
 
 # Here is an animation of the solution, looking only at the $v$ variable.

docs/src/literate_tutorials/laplaces_equation_with_internal_dirichlet_conditions.jl

@@ -93,7 +93,7 @@ ICs = InternalConditions((x, y, t, u, p) -> zero(u),
 # one suitable guess is $u(x, y) = 100y$ with $u(1/2, y) = 0$ for $0 \leq y \leq 2/5$;
 # in fact, $u(x, y) = 100y$ is the solution of the problem without the internal condition. 
 # Let us now use this to define our initial condition. 
-initial_condition = zeros(DelaunayTriangulation.num_solid_vertices(tri))
+initial_condition = zeros(DelaunayTriangulation.num_points(tri))
 for i in each_solid_vertex(tri)
     x, y = get_point(tri, i)
     initial_condition[i] = ifelse(x == 1 / 2 && 0 ≤ y ≤ 2 / 5, 0, 100y)