diff --git a/examples/build_bisected_rectangle.jl b/examples/build_bisected_rectangle.jl new file mode 100644 index 0000000..5479287 --- /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_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)) 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 bf351c2..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) @@ -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) diff --git a/src/standard_meshes.jl b/src/standard_meshes.jl index ac3dc37..879187a 100644 --- a/src/standard_meshes.jl +++ b/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 \ 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