From f0a5a35e097d47e06e38157fa7ac29868aef7a2a Mon Sep 17 00:00:00 2001
From: Benjamin Bolm <benjamin.bolm@gmx.de>
Date: Fri, 26 Jan 2024 13:24:49 +0100
Subject: [PATCH 1/9] Add `delaunay_compute_neighbors` routine

---
 src/Smesh.jl | 26 +++++++++++++++++++-------
 1 file changed, 19 insertions(+), 7 deletions(-)

diff --git a/src/Smesh.jl b/src/Smesh.jl
index 1731af7..8363939 100644
--- a/src/Smesh.jl
+++ b/src/Smesh.jl
@@ -1,8 +1,8 @@
 module Smesh
 
-using Preferences: @load_preference, @has_preference 
+using Preferences: @load_preference, @has_preference
 
-export build_delaunay_triangulation
+export build_delaunay_triangulation, delaunay_compute_neighbors
 
 if !@has_preference("libsmesh")
     error("""
@@ -23,27 +23,39 @@ if !@has_preference("libsmesh")
 end
 const libsmesh = @load_preference("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
+    # TODO: Converting to Julia `Int`s for convenience.
+    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
 end # module Smesh

From 7377f2686ab2873c611effdcb3506228e86197e9 Mon Sep 17 00:00:00 2001
From: Benjamin Bolm <benjamin.bolm@gmx.de>
Date: Fri, 26 Jan 2024 14:05:50 +0100
Subject: [PATCH 2/9] Add test for `delaunay_compute_neighbors`

---
 examples/build_delaunay_triangulation.jl |  2 ++
 test/test_unit.jl                        | 10 ++++++++++
 2 files changed, 12 insertions(+)

diff --git a/examples/build_delaunay_triangulation.jl b/examples/build_delaunay_triangulation.jl
index bda809e..df5e1cf 100644
--- a/examples/build_delaunay_triangulation.jl
+++ b/examples/build_delaunay_triangulation.jl
@@ -10,3 +10,5 @@ data_points = collect([0.0 0.0
 
 # Create triangulation
 ve = build_delaunay_triangulation(data_points; verbose = true)
+
+ne = delaunay_compute_neighbors(data_points, ve)
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 021b37b..c19276d 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -14,6 +14,16 @@ 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
+
 end # @testset "test_unit.jl"
 
 end # module

From 6ee15ea09b419ba636b59b234de71bf1868dc38a Mon Sep 17 00:00:00 2001
From: Benjamin Bolm <benjamin.bolm@gmx.de>
Date: Tue, 6 Feb 2024 12:35:46 +0100
Subject: [PATCH 3/9] Add routines for polygon mesh

---
 examples/build_delaunay_triangulation.jl | 26 ++++++++-
 examples/build_polygon_mesh.jl           | 69 +++++++++++++++++++++++
 src/Smesh.jl                             | 72 +++++++++++++++++++++++-
 3 files changed, 164 insertions(+), 3 deletions(-)
 create mode 100644 examples/build_polygon_mesh.jl

diff --git a/examples/build_delaunay_triangulation.jl b/examples/build_delaunay_triangulation.jl
index df5e1cf..96f81cd 100644
--- a/examples/build_delaunay_triangulation.jl
+++ b/examples/build_delaunay_triangulation.jl
@@ -9,6 +9,28 @@ data_points = collect([0.0 0.0
                        0.0 1.0]')
 
 # Create triangulation
-ve = build_delaunay_triangulation(data_points; verbose = true)
+vertices = build_delaunay_triangulation(data_points; verbose = true)
 
-ne = delaunay_compute_neighbors(data_points, ve)
+neighbors = delaunay_compute_neighbors(data_points, vertices)
+
+### Plotting ###
+using Plots; gr()
+
+# Vertices
+p = scatter(data_points[1, :], data_points[2, :], legend=false, title="vertices")
+display(p)
+
+# Triangles
+for element in axes(vertices, 2)
+    v1 = vertices[1, element]
+    v2 = vertices[2, element]
+    v3 = vertices[3, element]
+    x1 = data_points[:, v1]
+    x2 = data_points[:, v2]
+    x3 = data_points[:, v3]
+    plot!(p, [x1[1], x2[1]], [x1[2], x2[2]])
+    plot!(p, [x2[1], x3[1]], [x2[2], x3[2]])
+    plot!(p, [x3[1], x1[1]], [x3[2], x1[2]])
+end
+plot!(title = "triangles")
+display(p)
diff --git a/examples/build_polygon_mesh.jl b/examples/build_polygon_mesh.jl
new file mode 100644
index 0000000..b252c4e
--- /dev/null
+++ b/examples/build_polygon_mesh.jl
@@ -0,0 +1,69 @@
+using Smesh
+
+# Create data points
+corners = collect([0.0 0.0
+                   1.0 0.0
+                   1.0 1.0
+                   0.0 1.0]')
+# inner points (randomly generated)
+# n_points = 10
+# data_points = rand(Float64, 2, n_points)
+data_points = [0.110127  0.995047  0.636537  0.942174   0.22912   0.162025  0.616885  0.376891  0.475242  0.448486;
+               0.554234  0.431985  0.540326  0.0252587  0.702442  0.379256  0.80191   0.237447  0.745391  0.868326]
+data_points = hcat(data_points, corners)
+
+# 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
+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)
+
+##### Plotting
+using Plots; gr()
+
+# Vertices
+p = scatter(data_points[1, :], data_points[2, :], legend=false)
+
+# Triangles
+# p = scatter(data_points, legend=false)
+for element in axes(vertices, 2)
+    v1 = vertices[1, element]
+    v2 = vertices[2, element]
+    v3 = vertices[3, element]
+    x1 = data_points[:, v1]
+    x2 = data_points[:, v2]
+    x3 = data_points[:, v3]
+    plot!(p, [x1[1], x2[1]], [x1[2], x2[2]])
+    plot!(p, [x2[1], x3[1]], [x2[2], x3[2]])
+    plot!(p, [x3[1], x1[1]], [x3[2], x1[2]])
+end
+plot!(p, title="vertices and triangles")
+display(p)
+
+# Polygon mesh
+p1 = scatter(data_points[1, :], data_points[2, :], legend=false)
+for element_vor in axes(voronoi_vertices_interval, 2)
+    vertex_first = voronoi_vertices_interval[1, element_vor]
+    vertex_last = voronoi_vertices_interval[2, element_vor]
+    for i in vertex_first:(vertex_last-1)
+        v1 = voronoi_vertices[i]
+        v2 = voronoi_vertices[i + 1]
+        x1 = voronoi_vertices_coordinates[:, v1]
+        x2 = voronoi_vertices_coordinates[:, v2]
+        plot!(p1, [x1[1], x2[1]], [x1[2], x2[2]])
+    end
+    v1 = voronoi_vertices[vertex_last]
+    v2 = voronoi_vertices[vertex_first]
+    x1 = voronoi_vertices_coordinates[:, v1]
+    x2 = voronoi_vertices_coordinates[:, v2]
+    plot!(p1, [x1[1], x2[1]], [x1[2], x2[2]], title="mesh type: $mesh_type")
+end
+display(p1)
diff --git a/src/Smesh.jl b/src/Smesh.jl
index 67ca629..5caa2e1 100644
--- a/src/Smesh.jl
+++ b/src/Smesh.jl
@@ -3,7 +3,7 @@ module Smesh
 using Preferences: @load_preference
 using smesh_jll: smesh_jll
 
-export build_delaunay_triangulation, delaunay_compute_neighbors
+export build_delaunay_triangulation, delaunay_compute_neighbors, build_polygon_mesh, voronoi_compute_neighbors
 
 const libsmesh = @load_preference("libsmesh", smesh_jll.libsmesh)
 
@@ -42,4 +42,74 @@ function delaunay_compute_neighbors(data_points, vertices)
 
     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
+
+"""
+function build_polygon_mesh(data_points, triangulation_vertices; mesh_type=:standard_voronoi, orthogonal_boundary_edges=true)
+    mesh_type_dict = Dict(: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
 end # module Smesh

From c9757a9277fc676254373338ed8a3480aef271d6 Mon Sep 17 00:00:00 2001
From: Benjamin Bolm <benjamin.bolm@gmx.de>
Date: Tue, 6 Feb 2024 12:37:26 +0100
Subject: [PATCH 4/9] Remove comment

---
 src/Smesh.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/Smesh.jl b/src/Smesh.jl
index 5caa2e1..ec2136d 100644
--- a/src/Smesh.jl
+++ b/src/Smesh.jl
@@ -24,7 +24,6 @@ function build_delaunay_triangulation(data_points; shuffle = false, verbose = fa
                                                                 verbose::Cint)::Cint
 
     # Resize array to appropriate size
-    # TODO: Converting to Julia `Int`s for convenience.
     ve_out = ve_out[:, 1:ntriangles]
 
     return ve_out

From 7b0a05af2e5ee76cddbebd7bad8155d96610f32a Mon Sep 17 00:00:00 2001
From: Benjamin Bolm <benjamin.bolm@gmx.de>
Date: Tue, 6 Feb 2024 12:48:25 +0100
Subject: [PATCH 5/9] Add tests

---
 test/test_examples.jl |  4 ++++
 test/test_unit.jl     | 25 +++++++++++++++++++++++++
 2 files changed, 29 insertions(+)

diff --git a/test/test_examples.jl b/test/test_examples.jl
index 4ad85fe..429d346 100644
--- a/test/test_examples.jl
+++ b/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
diff --git a/test/test_unit.jl b/test/test_unit.jl
index c19276d..15f1428 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -24,6 +24,31 @@ end
     @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

From 5526dd3edcd8eb3cecf6f9f7c7acd05624bc93aa Mon Sep 17 00:00:00 2001
From: Benjamin Bolm <benjamin.bolm@gmx.de>
Date: Tue, 6 Feb 2024 12:50:12 +0100
Subject: [PATCH 6/9] Remove plotting

---
 examples/build_delaunay_triangulation.jl | 22 -------------
 examples/build_polygon_mesh.jl           | 42 ------------------------
 2 files changed, 64 deletions(-)

diff --git a/examples/build_delaunay_triangulation.jl b/examples/build_delaunay_triangulation.jl
index 96f81cd..54e53d4 100644
--- a/examples/build_delaunay_triangulation.jl
+++ b/examples/build_delaunay_triangulation.jl
@@ -12,25 +12,3 @@ data_points = collect([0.0 0.0
 vertices = build_delaunay_triangulation(data_points; verbose = true)
 
 neighbors = delaunay_compute_neighbors(data_points, vertices)
-
-### Plotting ###
-using Plots; gr()
-
-# Vertices
-p = scatter(data_points[1, :], data_points[2, :], legend=false, title="vertices")
-display(p)
-
-# Triangles
-for element in axes(vertices, 2)
-    v1 = vertices[1, element]
-    v2 = vertices[2, element]
-    v3 = vertices[3, element]
-    x1 = data_points[:, v1]
-    x2 = data_points[:, v2]
-    x3 = data_points[:, v3]
-    plot!(p, [x1[1], x2[1]], [x1[2], x2[2]])
-    plot!(p, [x2[1], x3[1]], [x2[2], x3[2]])
-    plot!(p, [x3[1], x1[1]], [x3[2], x1[2]])
-end
-plot!(title = "triangles")
-display(p)
diff --git a/examples/build_polygon_mesh.jl b/examples/build_polygon_mesh.jl
index b252c4e..86677ee 100644
--- a/examples/build_polygon_mesh.jl
+++ b/examples/build_polygon_mesh.jl
@@ -25,45 +25,3 @@ 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)
-
-##### Plotting
-using Plots; gr()
-
-# Vertices
-p = scatter(data_points[1, :], data_points[2, :], legend=false)
-
-# Triangles
-# p = scatter(data_points, legend=false)
-for element in axes(vertices, 2)
-    v1 = vertices[1, element]
-    v2 = vertices[2, element]
-    v3 = vertices[3, element]
-    x1 = data_points[:, v1]
-    x2 = data_points[:, v2]
-    x3 = data_points[:, v3]
-    plot!(p, [x1[1], x2[1]], [x1[2], x2[2]])
-    plot!(p, [x2[1], x3[1]], [x2[2], x3[2]])
-    plot!(p, [x3[1], x1[1]], [x3[2], x1[2]])
-end
-plot!(p, title="vertices and triangles")
-display(p)
-
-# Polygon mesh
-p1 = scatter(data_points[1, :], data_points[2, :], legend=false)
-for element_vor in axes(voronoi_vertices_interval, 2)
-    vertex_first = voronoi_vertices_interval[1, element_vor]
-    vertex_last = voronoi_vertices_interval[2, element_vor]
-    for i in vertex_first:(vertex_last-1)
-        v1 = voronoi_vertices[i]
-        v2 = voronoi_vertices[i + 1]
-        x1 = voronoi_vertices_coordinates[:, v1]
-        x2 = voronoi_vertices_coordinates[:, v2]
-        plot!(p1, [x1[1], x2[1]], [x1[2], x2[2]])
-    end
-    v1 = voronoi_vertices[vertex_last]
-    v2 = voronoi_vertices[vertex_first]
-    x1 = voronoi_vertices_coordinates[:, v1]
-    x2 = voronoi_vertices_coordinates[:, v2]
-    plot!(p1, [x1[1], x2[1]], [x1[2], x2[2]], title="mesh type: $mesh_type")
-end
-display(p1)

From 2b0e6e71965bdd3dca326839567ad491279feeae Mon Sep 17 00:00:00 2001
From: Benjamin Bolm <benjamin.bolm@gmx.de>
Date: Tue, 6 Feb 2024 13:55:57 +0100
Subject: [PATCH 7/9] Remove empty line

---
 src/Smesh.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/Smesh.jl b/src/Smesh.jl
index ec2136d..ca69eb0 100644
--- a/src/Smesh.jl
+++ b/src/Smesh.jl
@@ -93,8 +93,9 @@ function build_polygon_mesh(data_points, triangulation_vertices; mesh_type=:stan
     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)

From 1160e6d2896103684faa7b44747ffe932c3f237a Mon Sep 17 00:00:00 2001
From: Benjamin Bolm <benjamin.bolm@gmx.de>
Date: Tue, 6 Feb 2024 15:57:34 +0100
Subject: [PATCH 8/9] Add basic mesh and use it in examples

---
 examples/build_delaunay_triangulation.jl | 11 +++++------
 examples/build_polygon_mesh.jl           | 15 +++++---------
 src/Smesh.jl                             |  3 +++
 src/standard_meshes.jl                   | 25 ++++++++++++++++++++++++
 4 files changed, 38 insertions(+), 16 deletions(-)
 create mode 100644 src/standard_meshes.jl

diff --git a/examples/build_delaunay_triangulation.jl b/examples/build_delaunay_triangulation.jl
index 54e53d4..568ca90 100644
--- a/examples/build_delaunay_triangulation.jl
+++ b/examples/build_delaunay_triangulation.jl
@@ -1,12 +1,11 @@
 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
 vertices = build_delaunay_triangulation(data_points; verbose = true)
diff --git a/examples/build_polygon_mesh.jl b/examples/build_polygon_mesh.jl
index 86677ee..0762187 100644
--- a/examples/build_polygon_mesh.jl
+++ b/examples/build_polygon_mesh.jl
@@ -1,16 +1,11 @@
 using Smesh
 
 # Create data points
-corners = collect([0.0 0.0
-                   1.0 0.0
-                   1.0 1.0
-                   0.0 1.0]')
-# inner points (randomly generated)
-# n_points = 10
-# data_points = rand(Float64, 2, n_points)
-data_points = [0.110127  0.995047  0.636537  0.942174   0.22912   0.162025  0.616885  0.376891  0.475242  0.448486;
-               0.554234  0.431985  0.540326  0.0252587  0.702442  0.379256  0.80191   0.237447  0.745391  0.868326]
-data_points = hcat(data_points, corners)
+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)
diff --git a/src/Smesh.jl b/src/Smesh.jl
index ca69eb0..8bbfb2d 100644
--- a/src/Smesh.jl
+++ b/src/Smesh.jl
@@ -4,6 +4,7 @@ using Preferences: @load_preference
 using smesh_jll: smesh_jll
 
 export build_delaunay_triangulation, delaunay_compute_neighbors, build_polygon_mesh, voronoi_compute_neighbors
+export mesh_basic
 
 const libsmesh = @load_preference("libsmesh", smesh_jll.libsmesh)
 
@@ -112,4 +113,6 @@ function voronoi_compute_neighbors(vertices, voronoi_vertices, voronoi_vertices_
 
     return voronoi_neighbors
 end
+
+include("standard_meshes.jl")
 end # module Smesh
diff --git a/src/standard_meshes.jl b/src/standard_meshes.jl
new file mode 100644
index 0000000..26f65e5
--- /dev/null
+++ b/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

From 83b4faec705d6d2afc41f911632a40708e7ba42e Mon Sep 17 00:00:00 2001
From: Benjamin Bolm <benjamin.bolm@gmx.de>
Date: Tue, 6 Feb 2024 16:27:28 +0100
Subject: [PATCH 9/9] Add pure Voronoi mesh

---
 examples/build_polygon_mesh.jl | 7 ++++---
 src/Smesh.jl                   | 3 ++-
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/examples/build_polygon_mesh.jl b/examples/build_polygon_mesh.jl
index 0762187..4a82e7c 100644
--- a/examples/build_polygon_mesh.jl
+++ b/examples/build_polygon_mesh.jl
@@ -13,9 +13,10 @@ 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
+# :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)
 
diff --git a/src/Smesh.jl b/src/Smesh.jl
index 8bbfb2d..da23e96 100644
--- a/src/Smesh.jl
+++ b/src/Smesh.jl
@@ -50,10 +50,11 @@ 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(:standard_voronoi => Cint(0), :centroids => Cint(1), :incenters => Cint(2))
+    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