WIP: Add all routines to set up a voronoi mesh

examples/build_delaunay_triangulation.jl

@@ -1,12 +1,13 @@
 using Smesh
 
 # Create data points
-# Note: the transpose + collect is just such that we can write the matrix in human readable
-# form here
-data_points = collect([0.0 0.0
-                       1.0 0.0
-                       1.0 1.0
-                       0.0 1.0]')
+coordinates_min = [0.0, 0.0]
+coordinates_max = [1.0, 1.0]
+n_points_x = 4
+n_points_y = 5
+data_points = mesh_basic(coordinates_min, coordinates_max, n_points_x, n_points_y)
 
 # Create triangulation
-ve = build_delaunay_triangulation(data_points; verbose = true)
+vertices = build_delaunay_triangulation(data_points; verbose = true)
+
+neighbors = delaunay_compute_neighbors(data_points, vertices)

examples/build_polygon_mesh.jl

@@ -0,0 +1,23 @@
+using Smesh
+
+# Create data points
+coordinates_min = [0.0, 0.0]
+coordinates_max = [1.0, 1.0]
+n_points_x = 4
+n_points_y = 5
+data_points = mesh_basic(coordinates_min, coordinates_max, n_points_x, n_points_y)
+
+# Create triangulation
+vertices = build_delaunay_triangulation(data_points; verbose = false)
+
+neighbors = delaunay_compute_neighbors(data_points, vertices)
+
+# 3 options for the mesh type
+# :standard_voronoi => standard Voronoi, but use centroid if the circumcenter lies outside the triangle
+# :centroids        => not an actual Voronoi, always use centroids and not circumcenters as vertices for the mesh
+# :incenters        => not an actual Voronoi, always use incenters and not circumcenters as vertices for the mesh
+# :pure_voronoi     => pure Voronoi mesh (just for experiments, should not be used for computation)
+mesh_type = :standard_voronoi
+voronoi_vertices_coordinates, voronoi_vertices, voronoi_vertices_interval = build_polygon_mesh(data_points, vertices, mesh_type=mesh_type)
+
+voronoi_neighbors = voronoi_compute_neighbors(vertices, voronoi_vertices, voronoi_vertices_interval, neighbors)

src/Smesh.jl

@@ -1,33 +1,119 @@
 module Smesh
 
-using Preferences: @load_preference, @has_preference 
+using Preferences: @load_preference
 using smesh_jll: smesh_jll
 
-export build_delaunay_triangulation
+export build_delaunay_triangulation, delaunay_compute_neighbors, build_polygon_mesh, voronoi_compute_neighbors
+export mesh_basic
 
 const libsmesh = @load_preference("libsmesh", smesh_jll.libsmesh)
 
-
 """
 """
 function build_delaunay_triangulation(data_points; shuffle = false, verbose = false)
     # Pre-allocate output array
     npoints = size(data_points, 2)
     ve_max = @ccall libsmesh.delaunay_triangulation_temparray_size_c(npoints::Cint)::Cint
-    ve_internal = Matrix{Cint}(undef, 3, ve_max)
+    ve_out = Matrix{Cint}(undef, 3, ve_max)
 
     # Perform triangulation
-    ntriangles = @ccall libsmesh.build_delaunay_triangulation_c(ve_internal::Ref{Cint},
+    ntriangles = @ccall libsmesh.build_delaunay_triangulation_c(ve_out::Ref{Cint},
                                                                 data_points::Ref{Float64},
                                                                 npoints::Cint,
                                                                 ve_max::Cint,
                                                                 shuffle::Cint,
                                                                 verbose::Cint)::Cint
 
-    # Copy to array of appropriate size and convert to Julia `Int`s for convenience
-    ve_out = convert(Matrix{Int}, ve_internal[:, 1:ntriangles])
+    # Resize array to appropriate size
+    ve_out = ve_out[:, 1:ntriangles]
 
     return ve_out
 end
 
+"""
+"""
+function delaunay_compute_neighbors(data_points, vertices)
+    n_nodes = size(data_points, 2)
+    n_elements = size(vertices, 2)
+    neighbors = Matrix{Cint}(undef, 3, n_elements)
+
+    @ccall libsmesh.delaunay_compute_neighbors_c(neighbors::Ref{Cint}, vertices::Ref{Cint},
+                                                 n_elements::Cint, n_nodes::Cint)::Cvoid
+
+    return neighbors
+end
+
+"""
+    build_polygon_mesh(data_points, triangulation_vertices; mesh_type=:standard_voronoi, orthogonal_boundary_edges=true)
+
+There are three different mesh types:
+- `:standard_voronoi` => standard voronoi, but use centroid if the circumcenter lies outside the triangle
+- `:centroids` => not an actual voronoi, always use centroids and not circumcenters as vertices for the mesh
+- `:incenters` => not an actual voronoi, always use incenters and not circumcenters as vertices for the mesh
+- `:pure_voronoi` => pur Voronoi mesh (just for experiments, should not be used for computation)
+
+"""
+function build_polygon_mesh(data_points, triangulation_vertices; mesh_type=:standard_voronoi, orthogonal_boundary_edges=true)
+    mesh_type_dict = Dict(:pure_voronoi => Cint(-1), :standard_voronoi => Cint(0), :centroids => Cint(1), :incenters => Cint(2))
+
+    array_sizes = Vector{Cint}(undef, 3) # npt_voronoi, nve_voronoi, nelem_voronoi==nnode
+
+    npt_delaunay = size(data_points, 2)
+    nelem_delaunay = size(triangulation_vertices, 2)
+    nnode = npt_delaunay
+
+    orthogonal_boundary_edges_bool = orthogonal_boundary_edges ? 1 : 0
+
+    @ccall libsmesh.polygon_mesh_temparray_size_c(array_sizes::Ref{Cint},
+                                                  triangulation_vertices::Ref{Cint},
+                                                  data_points::Ref{Float64},
+                                                  mesh_type_dict[mesh_type]::Cint,
+                                                  orthogonal_boundary_edges_bool::Cint,
+                                                  npt_delaunay::Cint,
+                                                  nelem_delaunay::Cint,
+                                                  nnode::Cint)::Cvoid
+
+    npt_voronoi, nve_voronoi, nelem_voronoi = array_sizes
+
+    voronoi_vertices_coordinates = Matrix{Cdouble}(undef, 2, nve_voronoi)
+    voronoi_vertices = Array{Cint}(undef, nve_voronoi)
+    voronoi_vertices_interval = Matrix{Cint}(undef, 2, nelem_voronoi)
+
+    @ccall libsmesh.build_polygon_mesh_c(voronoi_vertices_coordinates::Ref{Float64},
+                                         voronoi_vertices::Ref{Cint},
+                                         voronoi_vertices_interval::Ref{Cint},
+                                         triangulation_vertices::Ref{Cint},
+                                         data_points::Ref{Float64},
+                                         mesh_type_dict[mesh_type]::Cint,
+                                         orthogonal_boundary_edges_bool::Cint,
+                                         npt_delaunay::Cint,
+                                         nelem_delaunay::Cint,
+                                         npt_voronoi::Cint,
+                                         nve_voronoi::Cint,
+                                         nelem_voronoi::Cint)::Cvoid
+
+    return voronoi_vertices_coordinates, voronoi_vertices, voronoi_vertices_interval
+end
+
+"""
+"""
+function voronoi_compute_neighbors(vertices, voronoi_vertices, voronoi_vertices_interval, delaunay_neighbors)
+    n_vertices_voronoi = length(voronoi_vertices)
+    n_elements_voronoi = size(voronoi_vertices_interval, 2)
+    n_element_delaunay = size(delaunay_neighbors, 2)
+
+    voronoi_neighbors = Vector{Cint}(undef, n_vertices_voronoi)
+
+    @ccall libsmesh.voronoi_compute_neighbors_c(voronoi_neighbors::Ref{Cint},
+                                                vertices::Ref{Cint},
+                                                voronoi_vertices::Ref{Cint},
+                                                voronoi_vertices_interval::Ref{Cint},
+                                                n_element_delaunay::Cint,
+                                                n_vertices_voronoi::Cint,
+                                                n_elements_voronoi::Cint)::Cvoid
+
+    return voronoi_neighbors
+end
+
+include("standard_meshes.jl")
 end # module Smesh

src/standard_meshes.jl

@@ -0,0 +1,25 @@
+function mesh_basic(coordinates_min, coordinates_max, n_points_x, n_points_y)
+    dx = (coordinates_max[1] - coordinates_min[1]) / n_points_x
+    dy = (coordinates_max[2] - coordinates_min[2]) / n_points_y
+
+    points = Matrix{eltype(coordinates_min)}(undef, 2, n_points_x * n_points_y  +
+                                                       div(n_points_y - n_points_y % 2, 2))
+    count = 1
+    for j in 1:n_points_y
+        for i in 1:n_points_x
+            points[1, count] = coordinates_min[1] + (i - 1) * dx
+            points[2, count] = coordinates_min[2] + (j - 1) * dy
+            if j % 2 == 0 && i != 1
+                points[1, count] -= 0.5dx
+            end
+            count += 1
+            if j % 2 == 0 && i == n_points_x
+                points[1, count] = points[1, count - 1] + 0.5dx
+                points[2, count] = points[2, count - 1]
+                count += 1
+            end
+        end
+    end
+
+    return points
+end

test/test_examples.jl

@@ -9,6 +9,10 @@ using Smesh
     @test_nowarn include("../examples/build_delaunay_triangulation.jl")
 end
 
+@testset verbose=true showtiming=true "examples/build_polygon_mesh.jl" begin
+    @test_nowarn include("../examples/build_polygon_mesh.jl")
+end
+
 end # @testset "test_examples.jl"
 
 end # module

test/test_unit.jl

@@ -14,6 +14,41 @@ using Smesh
     @test build_delaunay_triangulation(data_points) == [3 1; 1 3; 2 4]
 end
 
+@testset verbose=true showtiming=true "delaunay_compute_neighbors" begin
+    data_points = collect([0.0 0.0
+                           1.0 0.0
+                           1.0 1.0
+                           0.0 1.0]')
+    vertices = Cint[3 1; 1 3; 2 4]
+
+    @test delaunay_compute_neighbors(data_points, vertices) == [0 0; 0 0; 2 1]
+end
+
+@testset verbose=true showtiming=true "build_polygon_mesh" begin
+    data_points = collect([0.0 0.0
+                           1.0 0.0
+                           1.0 1.0
+                           0.0 1.0]')
+    vertices = Cint[3 1; 1 3; 2 4]
+    neighbors = Cint[0 0; 0 0; 2 1]
+
+    voronoi_vertices_coordinates, voronoi_vertices, voronoi_vertices_interval = build_polygon_mesh(data_points, vertices)
+    @test voronoi_vertices_interval == [1 7 12 18; 5 10 16 21]
+end
+
+@testset verbose=true showtiming=true "voronoi_compute_neighbors" begin
+    data_points = collect([0.0 0.0
+                           1.0 0.0
+                           1.0 1.0
+                           0.0 1.0]')
+    vertices = Cint[3 1; 1 3; 2 4]
+    neighbors = Cint[0 0; 0 0; 2 1]
+    voronoi_vertices_coordinates, voronoi_vertices, voronoi_vertices_interval = build_polygon_mesh(data_points, vertices)
+
+    voronoi_neighbor = voronoi_compute_neighbors(vertices, voronoi_vertices, voronoi_vertices_interval, neighbors)
+    @test voronoi_neighbor == Cint[3, 4, 0, 0, 2, 0, 1, 0, 0, 3, 0, 1, 2, 0, 0, 4, 0, 3, 0, 0, 1, 0]
+end
+
 end # @testset "test_unit.jl"
 
 end # module