Add mesh with bisected rectangles

examples/build_bisected_rectangle.jl

@@ -0,0 +1,22 @@
+using Smesh
+
+# Create data points
+coordinates_min = [0.0, 0.0]
+coordinates_max = [1.0, 1.0]
+n_elements_x = 5
+n_elements_y = 5
+data_points = mesh_bisected_rectangle(coordinates_min, coordinates_max, n_elements_x, n_elements_y,
+ symmetric_shift = true)
+
+# Create triangulation
+vertices = build_delaunay_triangulation(data_points; verbose = false, shuffle = false)
+
+neighbors = delaunay_compute_neighbors(data_points, vertices)
+
+mesh_type = :centroids
+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_coordinates,
+ voronoi_vertices, voronoi_vertices_interval,
+ neighbors, periodicity = (true, true))

examples/build_delaunay_triangulation.jl

@@ -10,4 +10,4 @@ data_points = mesh_basic(coordinates_min, coordinates_max, n_points_x, n_points_
 # Create triangulation
 vertices = build_delaunay_triangulation(data_points; verbose = true)
 
-neighbors = delaunay_compute_neighbors(data_points, vertices)
+neighbors = delaunay_compute_neighbors(data_points, vertices, periodicity=(true, true))

src/Smesh.jl

@@ -7,7 +7,7 @@ using LinearAlgebra: normalize
 
 export build_delaunay_triangulation, delaunay_compute_neighbors
 export build_polygon_mesh, voronoi_compute_neighbors
-export mesh_basic
+export mesh_basic, mesh_bisected_rectangle
 
 const libsmesh = @load_preference("libsmesh", smesh_jll.libsmesh)
 
@@ -101,7 +101,7 @@ function delaunay_compute_periodic_neighbors!(neighbors, periodicity, data_point
  end
  end
  # Check whether there are the same number of elements on both sides
- @assert length(boundary_elements_left) == length(boundary_elements_right) "Different number of elements at boundaries!"
+ @assert length(boundary_elements_left) == length(boundary_elements_right) "Different number of elements at boundaries in $dim-th direction!"
  @assert length(boundary_elements_left) != 0 "No detected boundary edge in $dim-th direction!"
  # Get coordinates for sorting
  # Note: In vertices the points are ordered counterclockwise:
@@ -124,12 +124,12 @@ function delaunay_compute_periodic_neighbors!(neighbors, periodicity, data_point
  for i in 1:(length(boundary_elements_left) - 1)
  face_length_left = abs(coord_elements_left[i] - coord_elements_left[i + 1])
  face_length_right = abs(coord_elements_right[i] - coord_elements_right[i + 1])
- @assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces do not match!"
+ @assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces in $dim-th direction do not match!"
  end
  # Check length of last boundary face
  face_length_left = abs(coord_elements_left[end] - data_points[dim % 2 + 1, vertices[boundary_faces_left[end] % 3 + 1, boundary_elements_left[end]]])
  face_length_right = abs(coord_elements_right[end] - data_points[dim % 2 + 1, vertices[(boundary_faces_right[end] + 1) % 3 + 1, boundary_elements_right[end]]])
- @assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces do not match!"
+ @assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces in $dim-th direction do not match!"
 
  # Add neighboring elements to neighbor data structure
  for i in eachindex(boundary_elements_left)
@@ -293,7 +293,7 @@ function voronoi_compute_periodic_neighbors!(voronoi_neighbors, periodicity,
  end
  end
  # Check whether there are the same number of elements on both sides
- @assert length(boundary_elements_left) == length(boundary_elements_right) "Different number of elements at boundaries!"
+ @assert length(boundary_elements_left) == length(boundary_elements_right) "Different number of elements at boundaries in $dim-th direction!"
  @assert length(boundary_elements_left) != 0 "No detected boundary edge in $dim-th direction!"
  # Get coordinates for sorting
  # Note: In voronoi_vertices the points are ordered counterclockwise:
@@ -318,12 +318,12 @@ function voronoi_compute_periodic_neighbors!(voronoi_neighbors, periodicity,
  for i in 1:(length(boundary_elements_left) - 1)
  face_length_left = abs(coord_elements_left[i] - coord_elements_left[i + 1])
  face_length_right = abs(coord_elements_right[i] - coord_elements_right[i + 1])
- @assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces do not match!"
+ @assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces in $dim-th direction do not match!"
  end
  # Check length of last boundary face.
  face_length_left = abs(coord_elements_left[end] - voronoi_vertices_coordinates[dim % 2 + 1, voronoi_vertices[boundary_faces_left[end]]])
  face_length_right = abs(coord_elements_right[end] - voronoi_vertices_coordinates[dim % 2 + 1, voronoi_vertices[boundary_faces_right[end] + 1]])
- @assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces do not match!"
+ @assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces in $dim-th direction do not match!"
 
  # Add neighboring elements to neighbor data structure
  for i in eachindex(boundary_elements_left)

src/standard_meshes.jl

@@ -1,9 +1,13 @@
 """
  mesh_basic(coordinates_min, coordinates_max, n_points_x, n_points_y)
 
-Create points in a regular grid.
+Creates points for a regular grid. Shifting every second column of points to avoid a simple
+mesh with bisected rectangles. This results in a unique triangulation.
 """
 function mesh_basic(coordinates_min, coordinates_max, n_points_x, n_points_y)
+ @assert n_points_x > 1 "n_points_x has to be at least 2."
+ @assert n_points_y > 1 "n_points_y has to be at least 2."
+
  dx = (coordinates_max[1] - coordinates_min[1]) / (n_points_x - 1)
  dy = (coordinates_max[2] - coordinates_min[2]) / (n_points_y - 1)
 
@@ -32,3 +36,101 @@ function mesh_basic(coordinates_min, coordinates_max, n_points_x, n_points_y)
 
  return points
 end
+
+"""
+ mesh_bisected_rectangle(coordinates_min, coordinates_max, n_points_x, n_points_y;
+ symmetric_shift = false)
+
+Creates points in a regular manner. The resulting non-unique triangulation consists of bisected
+rectangles. To allow periodic boundaries for the resulting polygon mesh, it is possible to enable
+a symmetric shift.
+"""
+function mesh_bisected_rectangle(coordinates_min, coordinates_max, n_elements_x, n_elements_y;
+ symmetric_shift = false)
+ @assert n_elements_x > 0 "n_elements_x has to be at least 1."
+ @assert n_elements_y > 0 "n_elements_y has to be at least 1."
+
+ dx = (coordinates_max[1] - coordinates_min[1]) / n_elements_x
+ dy = (coordinates_max[2] - coordinates_min[2]) / n_elements_y
+
+ n_points = (n_elements_x + 1) * (n_elements_y + 1)
+ points = Matrix{eltype(coordinates_min)}(undef, 2, n_points)
+ for j in 0:n_elements_y
+ for i = 0:n_elements_x
+ k = j * (n_elements_x + 1) + i + 1
+ points[:, k] = [coordinates_min[1] + i * dx, coordinates_min[2] + j * dy]
+ end
+ end
+
+ # Symmetric shift to get unique triangulation and therefore possible periodic boundaries in
+ # the polygon mesh
+ if symmetric_shift
+ domain_center = 0.5 * [coordinates_min[1] + coordinates_max[1],
+ coordinates_min[2] + coordinates_max[2]]
+ s = [dx, dy]
+ for i in axes(points, 2)
+ # Do not move boundary points with boundary_distance <= 10^-6
+ boundary_distance = min(abs(coordinates_min[1] - points[1, i]),
+ abs(coordinates_max[1] - points[1, i]),
+ abs(coordinates_min[2] - points[2, i]),
+ abs(coordinates_max[2] - points[2, i]))
+ if boundary_distance > 1.0e-8 # inner point
+ d = sqrt(sum((domain_center .- points[:,i]).^2))
+ points[:, i] .+= 1.0e-6 * d * s .* (domain_center .- points[:, i])
+ end
+ end
+
+ if isodd(n_elements_x)
+ for i in axes(points, 2)
+ # Do not move boundary points with boundary_distance <= 10^-6
+ boundary_distance = min(abs(coordinates_min[1] - points[1, i]),
+ abs(coordinates_max[1] - points[1, i]),
+ abs(coordinates_min[2] - points[2, i]),
+ abs(coordinates_max[2] - points[2, i]))
+ if boundary_distance > 1.0e-8 # inner point
+ # Only move the two most inner points columns
+ distance_center_x = abs(domain_center[1] - points[1, i])
+ if distance_center_x <= dx
+ points[1, i] += 1.0e-6 * dx
+ end
+ end
+ end
+ end
+ if isodd(n_elements_y)
+ for i in axes(points, 2)
+ # Do not move boundary points with boundary_distance <= 10^-6
+ boundary_distance = min(abs(coordinates_min[1] - points[1, i]),
+ abs(coordinates_max[1] - points[1, i]),
+ abs(coordinates_min[2] - points[2, i]),
+ abs(coordinates_max[2] - points[2, i]))
+ if boundary_distance > 1.0e-8 # inner point
+ # Only move the two most inner points rows
+ distance_center_y = abs(domain_center[2] - points[2, i])
+ if distance_center_y <= dy
+ points[2, i] += 1.0e-6 * dy
+ end
+ end
+ end
+ end
+ end
+
+ # This directly creates the connectivity of a triangulation. Every rectangle is bisected
+ # in the same direction.
+ # n_triangles = 2 * n_elements_x * n_elements_y
+ # vertices = Matrix{Cint}(undef, 3, n_triangles)
+ # k = 0
+ # for j in 1:n_elements_y
+ # for i in 1:n_elements_x
+ # k = k + 1
+ # vertices[:, k] .= [(j - 1) * (n_elements_x + 1) + i,
+ # (j - 1) * (n_elements_x + 1) + i + 1,
+ # j * (n_elements_x + 1) + i]
+ # k = k + 1
+ # vertices[:, k] .= [(j - 1) * (n_elements_x + 1) + i + 1,
+ # j * (n_elements_x + 1) + i + 1,
+ # j * (n_elements_x + 1) + i]
+ # end
+ # end
+
+ return points
+end

test/test_examples.jl

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

test/test_unit.jl

@@ -5,6 +5,19 @@ using Smesh
 
 @testset verbose=true showtiming=true "test_unit.jl" begin
 
+@testset verbose=true showtiming=true "meshes" begin
+ coordinates_min = [0.0, 0.0]
+ coordinates_max = [1.0, 1.0]
+ @testset verbose=true showtiming=true "mesh_basic" begin
+ points = mesh_basic(coordinates_min, coordinates_max, 2, 2)
+ @test points == [0.0 1.0 0.0 0.5 1.0; 0.0 0.0 1.0 1.0 1.0]
+ end
+ @testset verbose=true showtiming=true "mesh_bisected_rectangle" begin
+ points = mesh_bisected_rectangle(coordinates_min, coordinates_max, 1, 1)
+ @test points == [0.0 1.0 0.0 1.0; 0.0 0.0 1.0 1.0]
+ end
+end
+
 @testset verbose=true showtiming=true "build_delaunay_triangulation" begin
  data_points = collect([0.0 0.0
  1.0 0.0