From 486cf3beac45910bde47c7fbceee7c24cde527eb Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Mon, 26 Feb 2024 17:15:15 +0100 Subject: [PATCH 1/6] Adapt error messages --- src/Smesh.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Smesh.jl b/src/Smesh.jl index bf351c2..8b06ca1 100644 --- a/src/Smesh.jl +++ b/src/Smesh.jl @@ -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) From 6769eafa97790b5aa92ac84349ea0cfb961fbf6a Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Mon, 26 Feb 2024 19:13:18 +0100 Subject: [PATCH 2/6] Add check --- src/standard_meshes.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/standard_meshes.jl b/src/standard_meshes.jl index ac3dc37..3be5e34 100644 --- a/src/standard_meshes.jl +++ b/src/standard_meshes.jl @@ -4,6 +4,9 @@ Create points in a regular grid. """ 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) From db4205d0f7ac84dced9ac201c9a30cb6c92d18e8 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Mon, 26 Feb 2024 19:20:22 +0100 Subject: [PATCH 3/6] Add new mesh --- examples/build_bisected_rectangle.jl | 22 ++++++++ examples/build_delaunay_triangulation.jl | 2 +- src/Smesh.jl | 2 +- src/standard_meshes.jl | 66 +++++++++++++++++++++++- test/test_examples.jl | 4 ++ test/test_unit.jl | 13 +++++ 6 files changed, 106 insertions(+), 3 deletions(-) create mode 100644 examples/build_bisected_rectangle.jl diff --git a/examples/build_bisected_rectangle.jl b/examples/build_bisected_rectangle.jl new file mode 100644 index 0000000..2e82a5f --- /dev/null +++ b/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_points_x = 4 +n_points_y = 4 +data_points = mesh_bisected_rectangle(coordinates_min, coordinates_max, n_points_x, n_points_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 = (false, false)) diff --git a/examples/build_delaunay_triangulation.jl b/examples/build_delaunay_triangulation.jl index 568ca90..95f1e06 100644 --- a/examples/build_delaunay_triangulation.jl +++ b/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)) diff --git a/src/Smesh.jl b/src/Smesh.jl index 8b06ca1..3f45a9f 100644 --- a/src/Smesh.jl +++ b/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) diff --git a/src/standard_meshes.jl b/src/standard_meshes.jl index 3be5e34..0f33622 100644 --- a/src/standard_meshes.jl +++ b/src/standard_meshes.jl @@ -1,7 +1,8 @@ """ 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." @@ -35,3 +36,66 @@ 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_points_x, n_points_y; + symmetric_shift = false) + # n_points_x -= 1 + # n_points_y -= 1 + @assert n_points_x > 0 "n_points_x has to be at least 1." + @assert n_points_y > 0 "n_points_y has to be at least 1." + + dx = (coordinates_max[1] - coordinates_min[1]) / n_points_x + dy = (coordinates_max[2] - coordinates_min[2]) / n_points_y + + n_points = (n_points_x + 1) * (n_points_y + 1) + points = Matrix{eltype(coordinates_min)}(undef, 2, n_points) + for j in 0:n_points_y + for i = 0:n_points_x + k = j * (n_points_x + 1) + i + 1 + points[:, k] = [coordinates_min[1] + i * dx, coordinates_min[2] + j * dy] + end + end + + 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) + d = sqrt(sum((domain_center .- points[:,i]).^2)) + u = max(sign(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])) - 1.0e-8), + 0.0) + points[:, i] .= points[:, i] .+ 1.0e-6 * u * d * s .* (domain_center .- points[:, i]) + end + end + + # This directly creates the connectivity of a triangulation. Every rectangle is bisected + # in the same direction. + # n_triangles = 2 * n_points_x * n_points_y + # vertices = Matrix{Cint}(undef, 3, n_triangles) + # k = 0 + # for j in 1:n_points_y + # for i in 1:n_points_x + # k = k + 1 + # vertices[:, k] .= [(j - 1) * (n_points_x + 1) + i, + # (j - 1) * (n_points_x + 1) + i + 1, + # j * (n_points_x + 1) + i] + # k = k + 1 + # vertices[:, k] .= [(j - 1) * (n_points_x + 1) + i + 1, + # j * (n_points_x + 1) + i + 1, + # j * (n_points_x + 1) + i] + # end + # end + + return points +end \ No newline at end of file diff --git a/test/test_examples.jl b/test/test_examples.jl index 429d346..c7f26c1 100644 --- a/test/test_examples.jl +++ b/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 diff --git a/test/test_unit.jl b/test/test_unit.jl index a68074e..7c71e64 100644 --- a/test/test_unit.jl +++ b/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 From c6585e00e2f6b21e17e17f7379d5b2a7868b51da Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Wed, 28 Feb 2024 14:21:47 +0100 Subject: [PATCH 4/6] Rename variable --- examples/build_bisected_rectangle.jl | 6 ++--- src/standard_meshes.jl | 40 +++++++++++++--------------- 2 files changed, 22 insertions(+), 24 deletions(-) diff --git a/examples/build_bisected_rectangle.jl b/examples/build_bisected_rectangle.jl index 2e82a5f..143115c 100644 --- a/examples/build_bisected_rectangle.jl +++ b/examples/build_bisected_rectangle.jl @@ -3,9 +3,9 @@ using Smesh # Create data points coordinates_min = [0.0, 0.0] coordinates_max = [1.0, 1.0] -n_points_x = 4 -n_points_y = 4 -data_points = mesh_bisected_rectangle(coordinates_min, coordinates_max, n_points_x, n_points_y, +n_elements_x = 4 +n_elements_y = 4 +data_points = mesh_bisected_rectangle(coordinates_min, coordinates_max, n_elements_x, n_elements_y, symmetric_shift = true) # Create triangulation diff --git a/src/standard_meshes.jl b/src/standard_meshes.jl index 0f33622..e94ffab 100644 --- a/src/standard_meshes.jl +++ b/src/standard_meshes.jl @@ -45,28 +45,26 @@ Creates points in a regular manner. The resulting non-unique triangulation consi 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_points_x, n_points_y; +function mesh_bisected_rectangle(coordinates_min, coordinates_max, n_elements_x, n_elements_y; symmetric_shift = false) - # n_points_x -= 1 - # n_points_y -= 1 - @assert n_points_x > 0 "n_points_x has to be at least 1." - @assert n_points_y > 0 "n_points_y has to be at least 1." + @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_points_x - dy = (coordinates_max[2] - coordinates_min[2]) / n_points_y + dx = (coordinates_max[1] - coordinates_min[1]) / n_elements_x + dy = (coordinates_max[2] - coordinates_min[2]) / n_elements_y - n_points = (n_points_x + 1) * (n_points_y + 1) + n_points = (n_elements_x + 1) * (n_elements_y + 1) points = Matrix{eltype(coordinates_min)}(undef, 2, n_points) - for j in 0:n_points_y - for i = 0:n_points_x - k = j * (n_points_x + 1) + i + 1 + 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 if symmetric_shift domain_center = 0.5 * [coordinates_min[1] + coordinates_max[1], - coordinates_min[2] + coordinates_max[2]] + coordinates_min[2] + coordinates_max[2]] s = [dx, dy] for i in axes(points, 2) d = sqrt(sum((domain_center .- points[:,i]).^2)) @@ -81,19 +79,19 @@ function mesh_bisected_rectangle(coordinates_min, coordinates_max, n_points_x, n # This directly creates the connectivity of a triangulation. Every rectangle is bisected # in the same direction. - # n_triangles = 2 * n_points_x * n_points_y + # n_triangles = 2 * n_elements_x * n_elements_y # vertices = Matrix{Cint}(undef, 3, n_triangles) # k = 0 - # for j in 1:n_points_y - # for i in 1:n_points_x + # for j in 1:n_elements_y + # for i in 1:n_elements_x # k = k + 1 - # vertices[:, k] .= [(j - 1) * (n_points_x + 1) + i, - # (j - 1) * (n_points_x + 1) + i + 1, - # j * (n_points_x + 1) + i] + # 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_points_x + 1) + i + 1, - # j * (n_points_x + 1) + i + 1, - # j * (n_points_x + 1) + i] + # 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 From 662a714302f550e51cf1d4512e2d5568c291dfd5 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Wed, 28 Feb 2024 15:30:00 +0100 Subject: [PATCH 5/6] Fix symmteric shift --- src/standard_meshes.jl | 51 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 44 insertions(+), 7 deletions(-) diff --git a/src/standard_meshes.jl b/src/standard_meshes.jl index e94ffab..879187a 100644 --- a/src/standard_meshes.jl +++ b/src/standard_meshes.jl @@ -62,18 +62,55 @@ function mesh_bisected_rectangle(coordinates_min, coordinates_max, n_elements_x, 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) - d = sqrt(sum((domain_center .- points[:,i]).^2)) - u = max(sign(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])) - 1.0e-8), - 0.0) - points[:, i] .= points[:, i] .+ 1.0e-6 * u * d * s .* (domain_center .- points[:, i]) + # 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 From a0e7e039b5fedf1565c3728b9a30de03a2b55262 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Wed, 28 Feb 2024 15:34:00 +0100 Subject: [PATCH 6/6] Cover symmetric shift for odd number of elements --- examples/build_bisected_rectangle.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/build_bisected_rectangle.jl b/examples/build_bisected_rectangle.jl index 143115c..5479287 100644 --- a/examples/build_bisected_rectangle.jl +++ b/examples/build_bisected_rectangle.jl @@ -3,8 +3,8 @@ using Smesh # Create data points coordinates_min = [0.0, 0.0] coordinates_max = [1.0, 1.0] -n_elements_x = 4 -n_elements_y = 4 +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) @@ -19,4 +19,4 @@ voronoi_vertices_coordinates, voronoi_vertices, voronoi_neighbors = voronoi_compute_neighbors(vertices, voronoi_vertices_coordinates, voronoi_vertices, voronoi_vertices_interval, - neighbors, periodicity = (false, false)) + neighbors, periodicity = (true, true))