SciML · DanielVandH · Aug 8, 2023 · Aug 8, 2023 · Aug 8, 2023 · Aug 8, 2023
diff --git a/.gitignore b/.gitignore
@@ -1,5 +1,4 @@
 /Manifest.toml
-/test/Manifest.toml
 test/meshgeometry.geo
 gmsh-4.9.4-Windows64/
 test/meshgeometry.msh

diff --git a/Project.toml b/Project.toml
@@ -15,14 +15,14 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-ChunkSplitters = "^0.1, 1.0, 2.0"
-DelaunayTriangulation = "0.7"
-DiffEqBase = "^6.0"
-FunctionWrappersWrappers = "^0.1"
-MuladdMacro = "^0.2"
-PreallocationTools = "^0.4"
-SciMLBase = "^1.77"
-julia = "^1"
+ChunkSplitters = "0.1, 1.0, 2.0"
+DelaunayTriangulation = "0.7, 0.8"
+DiffEqBase = "6"
+FunctionWrappersWrappers = "0.1"
+MuladdMacro = "0.2"
+PreallocationTools = "0.4"
+SciMLBase = "1"
+julia = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

diff --git a/docs/src/annulus.md b/docs/src/annulus.md
@@ -69,15 +69,13 @@ alg = TRBDF2(linsolve=KLUFactorization())
 sol = solve(prob, alg, parallel=false, saveat=0.2)
 
 ## Visualise 
-pt_mat = Matrix(points')
-T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
 fig = Figure(resolution=(2068.72f0, 686.64f0), fontsize=38)
-ax = Axis(fig[1, 1], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[1], colorrange=(-10, 20), colormap=:viridis)
-ax = Axis(fig[1, 2], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[6], colorrange=(-10, 20), colormap=:viridis)
-ax = Axis(fig[1, 3], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[11], colorrange=(-10, 20), colormap=:viridis)
+for (i, j) in zip(1:3, (1, 6, 11))
+    ax = Axis(fig[1, i], width=600, height=600)
+    tricontourf!(ax, tri, sol.u[j], levels=-10:2:40, colormap=:viridis)
+    tightlimits!(ax)
+end
+fig
 ```
 
 ```@raw html

diff --git a/docs/src/diffusion_equation.md b/docs/src/diffusion_equation.md
@@ -101,20 +101,13 @@ u: 11-element Vector{Vector{Float64}}:
 You can use `sol` as you would any other solution from `DifferentialEquations.jl` (e.g. `sol(t)` returns the solution at time `t`). To visualise the solution at the times `t = 0.0`, `t = 0.25`, and `t = 0.5`, the following code can be used:
 ```julia
 using CairoMakie
-T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
 fig = Figure(resolution=(2068.72f0, 686.64f0), fontsize=38)
-ax = Axis(fig[1, 1], width=600, height=600)
-xlims!(ax, a, b)
-ylims!(ax, c, d)
-mesh!(ax, points, T_mat, color=sol.u[1], colorrange=(0, 50), colormap=:matter)
-ax = Axis(fig[1, 2], width=600, height=600)
-xlims!(ax, a, b)
-ylims!(ax, c, d)
-mesh!(ax, points, T_mat, color=sol.u[6], colorrange=(0, 50), colormap=:matter)
-ax = Axis(fig[1, 3], width=600, height=600)
-xlims!(ax, a, b)
-ylims!(ax, c, d)
-mesh!(ax, points, T_mat, color=sol.u[11], colorrange=(0, 50), colormap=:matter)
+for (i, j) in zip(1:3, (1, 6, 11))
+    ax = Axis(fig[1, i], width=600, height=600)
+    tricontourf!(ax, tri, sol.u[j], levels=0:5:50, colormap=:matter)
+    tightlimits!(ax)
+end
+fig
 ```
 
 ```@raw html

diff --git a/docs/src/diffusion_equation_on_a_wedge.md b/docs/src/diffusion_equation_on_a_wedge.md
@@ -100,16 +100,13 @@ sol = solve(prob, alg; saveat=0.025)
 
 ```julia
 using CairoMakie 
-
-pt_mat = Matrix(points')
-T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
-fig = Figure(resolution=(2131.8438f0, 684.27f0), fontsize=38)
-ax = Axis(fig[1, 1], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[1], colorrange=(0, 0.5), colormap=:matter)
-ax = Axis(fig[1, 2], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[3], colorrange=(0, 0.5), colormap=:matter)
-ax = Axis(fig[1, 3], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[5], colorrange=(0, 0.5), colormap=:matter)
+fig = Figure(resolution=(2068.72f0, 686.64f0), fontsize=38)
+for (i, j) in zip(1:3, (1, 3, 5))
+    ax = Axis(fig[1, i], width=600, height=600)
+    tricontourf!(ax, tri, sol.u[j], levels=0:0.01:1, colormap=:matter)
+    tightlimits!(ax)
+end
+fig
 ```
 
 ```@raw html

diff --git a/docs/src/figures/annulus_mean_exit_time.png b/docs/src/figures/annulus_mean_exit_time.png
diff --git a/docs/src/figures/annulus_test.png b/docs/src/figures/annulus_test.png
diff --git a/docs/src/figures/circle_mean_exit_time.png b/docs/src/figures/circle_mean_exit_time.png
diff --git a/docs/src/figures/diffusion_equation_wedge_test.png b/docs/src/figures/diffusion_equation_wedge_test.png
diff --git a/docs/src/figures/ellipse_mean_exit_time.png b/docs/src/figures/ellipse_mean_exit_time.png
diff --git a/docs/src/figures/heat_equation_test.png b/docs/src/figures/heat_equation_test.png
diff --git a/docs/src/figures/heat_equation_test_error.png b/docs/src/figures/heat_equation_test_error.png
diff --git a/docs/src/figures/heat_equation_wedge_test_error.png b/docs/src/figures/heat_equation_wedge_test_error.png
diff --git a/docs/src/figures/laplace_equation.png b/docs/src/figures/laplace_equation.png
diff --git a/docs/src/figures/perturbed_annulus_mean_exit_time.png b/docs/src/figures/perturbed_annulus_mean_exit_time.png
diff --git a/docs/src/figures/perturbed_circle_mean_exit_time.png b/docs/src/figures/perturbed_circle_mean_exit_time.png
diff --git a/docs/src/figures/perturbed_ellipse_mean_exit_time.png b/docs/src/figures/perturbed_ellipse_mean_exit_time.png
diff --git a/docs/src/figures/porous_medium_linear_source_test.png b/docs/src/figures/porous_medium_linear_source_test.png
diff --git a/docs/src/figures/porous_medium_linear_source_test_error.png b/docs/src/figures/porous_medium_linear_source_test_error.png
diff --git a/docs/src/figures/porous_medium_test.png b/docs/src/figures/porous_medium_test.png
diff --git a/docs/src/figures/porous_medium_test_error.png b/docs/src/figures/porous_medium_test_error.png
diff --git a/docs/src/figures/reaction_diffusion_equation_test_error.png b/docs/src/figures/reaction_diffusion_equation_test_error.png
diff --git a/docs/src/figures/reaction_diffusion_test.png b/docs/src/figures/reaction_diffusion_test.png
diff --git a/docs/src/figures/surface_plots_travelling_wave.png b/docs/src/figures/surface_plots_travelling_wave.png
diff --git a/docs/src/figures/travelling_wave_problem_test.png b/docs/src/figures/travelling_wave_problem_test.png
diff --git a/docs/src/laplace.md b/docs/src/laplace.md
@@ -83,10 +83,10 @@ sol = solve(prob, alg, parallel=true)
 ## Plot 
 fig = Figure(fontsize=38)
 ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y")
-pt_mat = get_points(tri)
-T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
-msh = mesh!(ax, pt_mat, T_mat, color=sol, colorrange=(-15, 15))
+msh = tricontourf!(tri, sol.u, levels = -15:0.5:15)
+tightlimits!(ax)
 Colorbar(fig[1, 2], msh)
+fig
 ```
 
 ```@raw html

diff --git a/docs/src/mean_exit_time.md b/docs/src/mean_exit_time.md
@@ -82,10 +82,11 @@ sol = solve(prob, alg)
 # Visualise 
 fig = Figure(fontsize=38)
 ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y")
-pt_mat = get_points(tri)
-T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
-msh = mesh!(ax, pt_mat, T_mat, color=sol, colorrange=(0, 10000))
+msh = tricontourf!(tri, sol.u, levels = 0:500:10000)
+tightlimits!(ax)
 Colorbar(fig[1, 2], msh)
+resize_to_layout!(fig)
+fig
 
 # Test 
 exact = [(R^2 - norm(p)^2)/(4D) for p in each_point(tri)]
@@ -145,10 +146,11 @@ alg = DynamicSS(TRBDF2(linsolve=FastLUFactorization()))
 sol = solve(prob, alg, parallel = true)
 fig = Figure(fontsize=38)
 ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y")
-pt_mat = get_points(tri)
-T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
-msh = mesh!(ax, pt_mat, T_mat, color=sol, colorrange=(0, 10000))
+msh = tricontourf!(tri, sol.u, levels = 0:500:10000)
+tightlimits!(ax)
 Colorbar(fig[1, 2], msh)
+fig
+resize_to_layout!(fig)
 ```
 
 ```@raw html
@@ -210,11 +212,11 @@ prob = FVMProblem(mesh, BCs;
     steady=true)
 sol = solve(prob, DynamicSS(TRBDF2(linsolve=KrylovJL_GMRES())))
 fig = Figure(fontsize=38)
-ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y", width=600, height=300)
-pt_mat = get_points(tri)
-T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
-msh = mesh!(ax, pt_mat, T_mat, color=sol, colorrange=(0, 16000))
+ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y", width = 900, height = 600)
+msh = tricontourf!(tri, sol.u, levels = 0:500:16000)
+tightlimits!(ax)
 Colorbar(fig[1, 2], msh)
+fig
 resize_to_layout!(fig)
 exact = [a^2 * b^2 / (2 * D * (a^2 + b^2)) * (1 - x^2 / a^2 - y^2 / b^2) for (x, y) in each_point(tri)]
 error_fnc = (fvm_sol, exact_sol) -> 100abs(fvm_sol - exact_sol) / maximum(exact)
@@ -272,12 +274,12 @@ prob = FVMProblem(mesh, BCs;
     steady=true)
 sol = solve(prob, DynamicSS(TRBDF2(linsolve=KrylovJL_GMRES())))
 fig = Figure(fontsize=38)
-ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y", width=600, height=300)
-pt_mat = get_points(tri)
-T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
-msh = mesh!(ax, pt_mat, T_mat, color=sol, colorrange=(0, 16000))
+ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y", width = 900, height = 600)
+msh = tricontourf!(tri, sol.u, levels = 0:500:16000)
+tightlimits!(ax)
 Colorbar(fig[1, 2], msh)
 resize_to_layout!(fig)
+fig
 ```
 
 ```@raw html
@@ -338,11 +340,11 @@ prob = FVMProblem(mesh, BCs;
     steady=true)
 sol = solve(prob, DynamicSS(TRBDF2(linsolve=KrylovJL_GMRES())))
 fig = Figure(fontsize=38)
-ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y", width=600, height=300)
-pt_mat = get_points(tri)
-T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
-msh = mesh!(ax, pt_mat, T_mat, color=sol, colorrange=(0, 900))
+ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y")
+msh = tricontourf!(tri, sol.u, levels = 0:50:900)
+tightlimits!(ax)
 Colorbar(fig[1, 2], msh)
+fig
 resize_to_layout!(fig)
 exact = [(R₂^2 - norm(p)^2) / (4D) + R₁^2 * log(norm(p) / R₂) / (2D) for p in each_point(tri)]
 error_fnc = (fvm_sol, exact_sol) -> 100abs(fvm_sol - exact_sol) / maximum(exact)
@@ -408,11 +410,11 @@ prob = FVMProblem(mesh, BCs;
     steady=true)
 sol = solve(prob, DynamicSS(TRBDF2(linsolve=KrylovJL_GMRES())), parallel = true)
 fig = Figure(fontsize=38)
-ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y", width=600, height=300)
-pt_mat = get_points(tri)
-T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
-msh = mesh!(ax, pt_mat, T_mat, color=sol, colorrange=(0, 900))
+ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y")
+msh = tricontourf!(tri, sol.u, levels = 0:50:900)
+tightlimits!(ax)
 Colorbar(fig[1, 2], msh)
+fig
 resize_to_layout!(fig)
 ```
 

diff --git a/docs/src/porous_fisher.md b/docs/src/porous_fisher.md
@@ -121,16 +121,13 @@ The results we obtain are shown below, with the exact travelling wave from the o
 using CairoMakie 
 
 # The solution 
-pt_mat = points
-T_mat = [[T...] for T in each_solid_triangle(tri)] # each_solid since tri contains some ghost triangles
-T_mat = reduce(hcat, T_mat)'
-fig = Figure(resolution=(3023.5881f0, 684.27f0), fontsize=38)
-ax = Axis(fig[1, 1], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[1], colorrange=(0.0, 1.0), colormap=:matter)
-ax = Axis(fig[1, 2], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[51], colorrange=(0.0, 1.0), colormap=:matter)
-ax = Axis(fig[1, 3], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[101], colorrange=(0.0, 1.0), colormap=:matter)
+fig = Figure(resolution=(3024.72f0, 686.64f0), fontsize=38)
+for (i, j) in zip(1:3, (1, 51, 101))
+    ax = Axis(fig[1, i], width=600, height=600)
+    tricontourf!(ax, tri, sol.u[j], levels=0:0.05:1, colormap=:matter)
+    tightlimits!(ax)
+end
+fig
 
 # The wave comparisons 
 ax = Axis(fig[1, 4], width=900, height=600)

diff --git a/docs/src/porous_medium.md b/docs/src/porous_medium.md
@@ -67,16 +67,13 @@ sol = solve(prob, alg; saveat=3.0)
 
 ## Step 5: Visualisation
 using CairoMakie
-
-pt_mat = points
-T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
-fig = Figure(resolution=(2131.8438f0, 684.27f0), fontsize=38)
-ax = Axis(fig[1, 1], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[1], colorrange=(0.0, 0.05), colormap=:matter)
-ax = Axis(fig[1, 2], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[3], colorrange=(0.0, 0.05), colormap=:matter)
-ax = Axis(fig[1, 3], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[5], colorrange=(0.0, 0.05), colormap=:matter)
+fig = Figure(resolution=(2068.72f0, 686.64f0), fontsize=38)
+for (i, j) in zip(1:3, (1, 3, 5))
+    ax = Axis(fig[1, i], width=600, height=600)
+    tricontourf!(ax, tri, sol.u[j], levels=0:0.005:0.05, colormap=:matter, extendhigh=:auto)
+    tightlimits!(ax)
+end
+fig
 ```
 
 ```@raw html
@@ -93,7 +90,6 @@ Here is a detailed comparison with the exact solution.
 </figure>
 ```
 
-
 ## Linear source
 
 We can continue this example with the Porous medium equation by considering the same equation except with a linear source:
@@ -156,15 +152,13 @@ alg = TRBDF2(linsolve=KLUFactorization())
 sol = solve(prob, alg; saveat=2.5)
 
 ## Step 5: Visualisation 
-pt_mat = points
-T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
-fig = Figure(resolution=(2131.8438f0, 684.27f0), fontsize=38)
-ax = Axis(fig[1, 1], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[1], colorrange=(0.0, 0.5), colormap=:matter)
-ax = Axis(fig[1, 2], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[3], colorrange=(0.0, 0.5), colormap=:matter)
-ax = Axis(fig[1, 3], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[5], colorrange=(0.0, 0.5), colormap=:matter)
+fig = Figure(resolution=(2132.72f0, 686.64f0), fontsize=38)
+for (i, j) in zip(1:3, (1, 3, 5))
+    ax = Axis(fig[1, i], width=600, height=600)
+    tricontourf!(ax, tri, sol.u[j], levels=0:0.05:1, colormap=:matter, extendhigh=:auto)
+    tightlimits!(ax)
+end
+fig
 ```
 
 ```@raw html

diff --git a/docs/src/reaction_diffusion.md b/docs/src/reaction_diffusion.md
@@ -61,16 +61,13 @@ sol = solve(prob, alg; saveat=0.025)
 
 ## Step 5: Visualisation 
 using CairoMakie 
-
-pt_mat = Matrix(points')
-T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
-fig = Figure(resolution=(2131.8438f0, 684.27f0), fontsize=38)
-ax = Axis(fig[1, 1], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[1], colorrange=(1, 1.1), colormap=:matter)
-ax = Axis(fig[1, 2], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[3], colorrange=(1, 1.1), colormap=:matter)
-ax = Axis(fig[1, 3], width=600, height=600)
-mesh!(ax, pt_mat, T_mat, color=sol.u[5], colorrange=(1, 1.1), colormap=:matter)
+fig = Figure(resolution=(2068.72f0, 686.64f0), fontsize=38)
+for (i, j) in zip(1:3, (1, 3, 5))
+    ax = Axis(fig[1, i], width=600, height=600)
+    tricontourf!(ax, tri, sol.u[j], levels=1:0.01:1.4, colormap=:matter)
+    tightlimits!(ax)
+end
+fig
 ```
 
 ```@raw html