Change to v0.3 for DelaunayTriangulation. Fix missing using DelaunayT…

…riangulation in the examples (#21).

Project.toml

@@ -1,7 +1,7 @@
 name = "FiniteVolumeMethod"
 uuid = "d4f04ab7-4f65-4d72-8a28-7087bc7f46f4"
 authors = ["DanielVandH <danj.vandenheuvel@gmail.com> and contributors"]
-version = "0.3.2"
+version = "0.3.3"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
@@ -16,15 +16,15 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 
 [compat]
-DelaunayTriangulation = "^0.1, 0.2"
+DelaunayTriangulation = "0.3"
+ChunkSplitters = "^0.1"
 DiffEqBase = "^6.0"
 FunctionWrappersWrappers = "^0.1"
 MuladdMacro = "^0.2"
 PreallocationTools = "^0.4"
 SciMLBase = "^1.77"
 StaticArraysCore = "^1.4"
 julia = "^1"
-ChunkSplitters = "^0.1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -193,6 +193,7 @@ GMSH_PATH = "./gmsh-4.9.4-Windows64/gmsh.exe"
 We also have the following packages loaded:
 ```julia
 using FiniteVolumeMethod
+using DelaunayTriangulation
 using OrdinaryDiffEq 
 using LinearSolve 
 using CairoMakie 
@@ -243,7 +244,7 @@ unique!(xy)
 x = getx.(xy)
 y = gety.(xy)
 r = 0.03
-T, adj, adj2v, DG, points, BN = generate_mesh(x, y, r; gmsh_path=GMSH_PATH)
+(T, adj, adj2v, DG, points), BN = generate_mesh(x, y, r; gmsh_path=GMSH_PATH)
 mesh = FVMGeometry(T, adj, adj2v, DG, points, BN)
 ```
 Here I start by defining the square boundary as four segments, but then to have a single boundary segment I combine the segments into a single vector. I then create the mesh using `generate_mesh`, and then put the geometry together using `FVMGeometry`. 
@@ -374,8 +375,11 @@ y₃ = @. r₃ * sin(θ₃)
 x = [x₁, x₂, x₃]
 y = [y₁, y₂, y₃]
 r = 0.01
-T, adj, adj2v, DG, points, BN = generate_mesh(x, y, r; gmsh_path=GMSH_PATH)
+(T, adj, adj2v, DG, points), BN = generate_mesh(x, y, r; gmsh_path=GMSH_PATH)
 mesh = FVMGeometry(T, adj, adj2v, DG, points, BN)
+# You could also do:
+# tri, BN = generate_mesh(x, y, r; gmsh_path=GMSH_PATH)
+# mesh = FVMGeometry(tri, BN) # for DelaunayTriangulation >= v0.3)
 ```
 
 Now we define the boundary conditions.
@@ -438,7 +442,7 @@ r = LinRange(1, 1, 1000)
 x = @. r * cos(θ)
 y = @. r * sin(θ)
 r = 0.05
-T, adj, adj2v, DG, points, BN = generate_mesh(x, y, r; gmsh_path=GMSH_PATH)
+(T, adj, adj2v, DG, points), BN = generate_mesh(x, y, r; gmsh_path=GMSH_PATH)
 mesh = FVMGeometry(T, adj, adj2v, DG, points, BN)
 
 ## Step 2: Define the boundary conditions 
@@ -509,7 +513,7 @@ unique!(xy)
 x = getx.(xy)
 y = gety.(xy)
 r = 0.1
-T, adj, adj2v, DG, points, BN = generate_mesh(x, y, r; gmsh_path=GMSH_PATH)
+(T, adj, adj2v, DG, points), BN = generate_mesh(x, y, r; gmsh_path=GMSH_PATH)
 mesh = FVMGeometry(T, adj, adj2v, DG, points, BN)
 
 ## Step 2: Define the boundary conditions 
@@ -572,7 +576,7 @@ unique!(xy)
 x = getx.(xy)
 y = gety.(xy)
 r = 0.07
-T, adj, adj2v, DG, points, BN = generate_mesh(x, y, r; gmsh_path=GMSH_PATH)
+(T, adj, adj2v, DG, points), BN = generate_mesh(x, y, r; gmsh_path=GMSH_PATH)
 mesh = FVMGeometry(T, adj, adj2v, DG, points, BN)
 
 ## Step 2: Define the boundary conditions 
@@ -634,7 +638,7 @@ Let us now solve the problem. For this problem, rather than using `generate_mesh
 ```julia
 ## Step 1: Define the mesh 
 a, b, c, d, Nx, Ny = 0.0, 3.0, 0.0, 40.0, 60, 80
-T, adj, adj2v, DG, points, BN = triangulate_structured(a, b, c, d, Nx, Ny; return_boundary_types=true)
+(T, adj, adj2v, DG, points), BN = triangulate_structured(a, b, c, d, Nx, Ny; return_boundary_types=true)
 mesh = FVMGeometry(T, adj, adj2v, DG, points, BN)
 
 ## Step 2: Define the boundary conditions 

src/geometry.jl

@@ -228,6 +228,8 @@ Constructor for [`FVMGeometry`](@ref).
 - `pts`: The points for the mesh. 
 - `BNV`: The boundary node vector. This should be a vector of vectors, where each nested vector is a list of indices that define the nodes for the corresponding segment, and `first(BNV[i]) == last(BNV[i-1])`. The nodes must be listed in counter-clockwise order.
 
+There is also a constructor where these arguments are `FVMGeometry(tri::Triangulation, BNV; same kwargs as below)`.
+
 # Keyword Arguments 
 - `coordinate_type=Vector{number_type(pts)}`: How coordinates are represented.
 - `control_volume_storage_type_vector=NTuple{3,coordinate_type}`: How information for triples of coordinates is represented. The element type must be `coordinate_type`.
@@ -292,6 +294,28 @@ function FVMGeometry(T::Ts, adj, adj2v, DG, pts, BNV;
         element_information_list=element_infos,
         volumes=vols)
 end
+function FVMGeometry(tri::Triangulation, BNV;
+    coordinate_type=Vector{number_type(DelaunayTriangulation.get_points(tri))},
+    control_volume_storage_type_vector=NTuple{3,coordinate_type},
+    control_volume_storage_type_scalar=NTuple{3,number_type(DelaunayTriangulation.get_points(tri))},
+    shape_function_coefficient_storage_type=NTuple{9,number_type(DelaunayTriangulation.get_points(tri))},
+    interior_edge_storage_type=NTuple{2,Int64},
+    interior_edge_pair_storage_type=NTuple{2,interior_edge_storage_type})
+    return FVMGeometry(
+        DelaunayTriangulation.get_triangles(tri),
+        DelaunayTriangulation.get_adjacent(tri),
+        DelaunayTriangulation.get_adjacent2vertex(tri),
+        DelaunayTriangulation.get_graph(tri),
+        DelaunayTriangulation.get_points(tri),
+        BNV;
+        coordinate_type,
+        control_volume_storage_type_vector,
+        control_volume_storage_type_scalar,
+        shape_function_coefficient_storage_type,
+        interior_edge_storage_type,
+        interior_edge_pair_storage_type
+    )
+end
 
 function compute_centroid(p₁, p₂, p₃; coordinate_type)
     cx = (getx(p₁) + getx(p₂) + getx(p₃)) / 3
@@ -396,7 +420,7 @@ end
 
 function compute_sub_control_volume_area(p, q)
     F = number_type(p)
-    S = oneunit(F)/2 * abs(getx(p) * gety(q) - gety(p) * getx(q))
+    S = oneunit(F) / 2 * abs(getx(p) * gety(q) - gety(p) * getx(q))
     return S
 end
 function sub_control_volume_areas(p, q)