Skip to content

Commit 57d35ef

Browse files
bennibolmbennibolm
authored
Add mesh with bisected rectangles (#18)
* Adapt error messages * Add check * Add new mesh * Rename variable * Fix symmteric shift * Cover symmetric shift for odd number of elements
1 parent 1bd3b26 commit 57d35ef

File tree

6 files changed

+150
-9
lines changed

6 files changed

+150
-9
lines changed

‎examples/build_bisected_rectangle.jl‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
using Smesh
2+
3+
# Create data points
4+
coordinates_min = [0.0, 0.0]
5+
coordinates_max = [1.0, 1.0]
6+
n_elements_x = 5
7+
n_elements_y = 5
8+
data_points = mesh_bisected_rectangle(coordinates_min, coordinates_max, n_elements_x, n_elements_y,
9+
symmetric_shift = true)
10+
11+
# Create triangulation
12+
vertices = build_delaunay_triangulation(data_points; verbose = false, shuffle = false)
13+
14+
neighbors = delaunay_compute_neighbors(data_points, vertices)
15+
16+
mesh_type = :centroids
17+
voronoi_vertices_coordinates, voronoi_vertices,
18+
voronoi_vertices_interval = build_polygon_mesh(data_points, vertices, mesh_type=mesh_type)
19+
20+
voronoi_neighbors = voronoi_compute_neighbors(vertices, voronoi_vertices_coordinates,
21+
voronoi_vertices, voronoi_vertices_interval,
22+
neighbors, periodicity = (true, true))

‎examples/build_delaunay_triangulation.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,4 +10,4 @@ data_points = mesh_basic(coordinates_min, coordinates_max, n_points_x, n_points_
1010
# Create triangulation
1111
vertices = build_delaunay_triangulation(data_points; verbose = true)
1212

13-
neighbors = delaunay_compute_neighbors(data_points, vertices)
13+
neighbors = delaunay_compute_neighbors(data_points, vertices, periodicity=(true, true))

‎src/Smesh.jl‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ using LinearAlgebra: normalize
77

88
export build_delaunay_triangulation, delaunay_compute_neighbors
99
export build_polygon_mesh, voronoi_compute_neighbors
10-
export mesh_basic
10+
export mesh_basic, mesh_bisected_rectangle
1111

1212
const libsmesh = @load_preference("libsmesh", smesh_jll.libsmesh)
1313

@@ -101,7 +101,7 @@ function delaunay_compute_periodic_neighbors!(neighbors, periodicity, data_point
101101
end
102102
end
103103
# Check whether there are the same number of elements on both sides
104-
@assert length(boundary_elements_left) == length(boundary_elements_right) "Different number of elements at boundaries!"
104+
@assert length(boundary_elements_left) == length(boundary_elements_right) "Different number of elements at boundaries in $dim-th direction!"
105105
@assert length(boundary_elements_left) != 0 "No detected boundary edge in $dim-th direction!"
106106
# Get coordinates for sorting
107107
# Note: In vertices the points are ordered counterclockwise:
@@ -124,12 +124,12 @@ function delaunay_compute_periodic_neighbors!(neighbors, periodicity, data_point
124124
for i in 1:(length(boundary_elements_left) - 1)
125125
face_length_left = abs(coord_elements_left[i] - coord_elements_left[i + 1])
126126
face_length_right = abs(coord_elements_right[i] - coord_elements_right[i + 1])
127-
@assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces do not match!"
127+
@assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces in $dim-th direction do not match!"
128128
end
129129
# Check length of last boundary face
130130
face_length_left = abs(coord_elements_left[end] - data_points[dim % 2 + 1, vertices[boundary_faces_left[end] % 3 + 1, boundary_elements_left[end]]])
131131
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]]])
132-
@assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces do not match!"
132+
@assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces in $dim-th direction do not match!"
133133

134134
# Add neighboring elements to neighbor data structure
135135
for i in eachindex(boundary_elements_left)
@@ -293,7 +293,7 @@ function voronoi_compute_periodic_neighbors!(voronoi_neighbors, periodicity,
293293
end
294294
end
295295
# Check whether there are the same number of elements on both sides
296-
@assert length(boundary_elements_left) == length(boundary_elements_right) "Different number of elements at boundaries!"
296+
@assert length(boundary_elements_left) == length(boundary_elements_right) "Different number of elements at boundaries in $dim-th direction!"
297297
@assert length(boundary_elements_left) != 0 "No detected boundary edge in $dim-th direction!"
298298
# Get coordinates for sorting
299299
# Note: In voronoi_vertices the points are ordered counterclockwise:
@@ -318,12 +318,12 @@ function voronoi_compute_periodic_neighbors!(voronoi_neighbors, periodicity,
318318
for i in 1:(length(boundary_elements_left) - 1)
319319
face_length_left = abs(coord_elements_left[i] - coord_elements_left[i + 1])
320320
face_length_right = abs(coord_elements_right[i] - coord_elements_right[i + 1])
321-
@assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces do not match!"
321+
@assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces in $dim-th direction do not match!"
322322
end
323323
# Check length of last boundary face.
324324
face_length_left = abs(coord_elements_left[end] - voronoi_vertices_coordinates[dim % 2 + 1, voronoi_vertices[boundary_faces_left[end]]])
325325
face_length_right = abs(coord_elements_right[end] - voronoi_vertices_coordinates[dim % 2 + 1, voronoi_vertices[boundary_faces_right[end] + 1]])
326-
@assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces do not match!"
326+
@assert isapprox(face_length_left, face_length_right, atol=eps()) "Length of boundary faces in $dim-th direction do not match!"
327327

328328
# Add neighboring elements to neighbor data structure
329329
for i in eachindex(boundary_elements_left)

‎src/standard_meshes.jl‎

Lines changed: 103 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,13 @@
11
"""
22
mesh_basic(coordinates_min, coordinates_max, n_points_x, n_points_y)
33
4-
Create points in a regular grid.
4+
Creates points for a regular grid. Shifting every second column of points to avoid a simple
5+
mesh with bisected rectangles. This results in a unique triangulation.
56
"""
67
function mesh_basic(coordinates_min, coordinates_max, n_points_x, n_points_y)
8+
@assert n_points_x > 1 "n_points_x has to be at least 2."
9+
@assert n_points_y > 1 "n_points_y has to be at least 2."
10+
711
dx = (coordinates_max[1] - coordinates_min[1]) / (n_points_x - 1)
812
dy = (coordinates_max[2] - coordinates_min[2]) / (n_points_y - 1)
913

@@ -32,3 +36,101 @@ function mesh_basic(coordinates_min, coordinates_max, n_points_x, n_points_y)
3236

3337
return points
3438
end
39+
40+
"""
41+
mesh_bisected_rectangle(coordinates_min, coordinates_max, n_points_x, n_points_y;
42+
symmetric_shift = false)
43+
44+
Creates points in a regular manner. The resulting non-unique triangulation consists of bisected
45+
rectangles. To allow periodic boundaries for the resulting polygon mesh, it is possible to enable
46+
a symmetric shift.
47+
"""
48+
function mesh_bisected_rectangle(coordinates_min, coordinates_max, n_elements_x, n_elements_y;
49+
symmetric_shift = false)
50+
@assert n_elements_x > 0 "n_elements_x has to be at least 1."
51+
@assert n_elements_y > 0 "n_elements_y has to be at least 1."
52+
53+
dx = (coordinates_max[1] - coordinates_min[1]) / n_elements_x
54+
dy = (coordinates_max[2] - coordinates_min[2]) / n_elements_y
55+
56+
n_points = (n_elements_x + 1) * (n_elements_y + 1)
57+
points = Matrix{eltype(coordinates_min)}(undef, 2, n_points)
58+
for j in 0:n_elements_y
59+
for i = 0:n_elements_x
60+
k = j * (n_elements_x + 1) + i + 1
61+
points[:, k] = [coordinates_min[1] + i * dx, coordinates_min[2] + j * dy]
62+
end
63+
end
64+
65+
# Symmetric shift to get unique triangulation and therefore possible periodic boundaries in
66+
# the polygon mesh
67+
if symmetric_shift
68+
domain_center = 0.5 * [coordinates_min[1] + coordinates_max[1],
69+
coordinates_min[2] + coordinates_max[2]]
70+
s = [dx, dy]
71+
for i in axes(points, 2)
72+
# Do not move boundary points with boundary_distance <= 10^-6
73+
boundary_distance = min(abs(coordinates_min[1] - points[1, i]),
74+
abs(coordinates_max[1] - points[1, i]),
75+
abs(coordinates_min[2] - points[2, i]),
76+
abs(coordinates_max[2] - points[2, i]))
77+
if boundary_distance > 1.0e-8 # inner point
78+
d = sqrt(sum((domain_center .- points[:,i]).^2))
79+
points[:, i] .+= 1.0e-6 * d * s .* (domain_center .- points[:, i])
80+
end
81+
end
82+
83+
if isodd(n_elements_x)
84+
for i in axes(points, 2)
85+
# Do not move boundary points with boundary_distance <= 10^-6
86+
boundary_distance = min(abs(coordinates_min[1] - points[1, i]),
87+
abs(coordinates_max[1] - points[1, i]),
88+
abs(coordinates_min[2] - points[2, i]),
89+
abs(coordinates_max[2] - points[2, i]))
90+
if boundary_distance > 1.0e-8 # inner point
91+
# Only move the two most inner points columns
92+
distance_center_x = abs(domain_center[1] - points[1, i])
93+
if distance_center_x <= dx
94+
points[1, i] += 1.0e-6 * dx
95+
end
96+
end
97+
end
98+
end
99+
if isodd(n_elements_y)
100+
for i in axes(points, 2)
101+
# Do not move boundary points with boundary_distance <= 10^-6
102+
boundary_distance = min(abs(coordinates_min[1] - points[1, i]),
103+
abs(coordinates_max[1] - points[1, i]),
104+
abs(coordinates_min[2] - points[2, i]),
105+
abs(coordinates_max[2] - points[2, i]))
106+
if boundary_distance > 1.0e-8 # inner point
107+
# Only move the two most inner points rows
108+
distance_center_y = abs(domain_center[2] - points[2, i])
109+
if distance_center_y <= dy
110+
points[2, i] += 1.0e-6 * dy
111+
end
112+
end
113+
end
114+
end
115+
end
116+
117+
# This directly creates the connectivity of a triangulation. Every rectangle is bisected
118+
# in the same direction.
119+
# n_triangles = 2 * n_elements_x * n_elements_y
120+
# vertices = Matrix{Cint}(undef, 3, n_triangles)
121+
# k = 0
122+
# for j in 1:n_elements_y
123+
# for i in 1:n_elements_x
124+
# k = k + 1
125+
# vertices[:, k] .= [(j - 1) * (n_elements_x + 1) + i,
126+
# (j - 1) * (n_elements_x + 1) + i + 1,
127+
# j * (n_elements_x + 1) + i]
128+
# k = k + 1
129+
# vertices[:, k] .= [(j - 1) * (n_elements_x + 1) + i + 1,
130+
# j * (n_elements_x + 1) + i + 1,
131+
# j * (n_elements_x + 1) + i]
132+
# end
133+
# end
134+
135+
return points
136+
end

‎test/test_examples.jl‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,10 @@ end
1313
@test_nowarn include("../examples/build_polygon_mesh.jl")
1414
end
1515

16+
@testset verbose=true showtiming=true "examples/build_bisected_rectangle.jl" begin
17+
@test_nowarn include("../examples/build_bisected_rectangle.jl")
18+
end
19+
1620
end # @testset "test_examples.jl"
1721

1822
end # module

‎test/test_unit.jl‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,19 @@ using Smesh
55

66
@testset verbose=true showtiming=true "test_unit.jl" begin
77

8+
@testset verbose=true showtiming=true "meshes" begin
9+
coordinates_min = [0.0, 0.0]
10+
coordinates_max = [1.0, 1.0]
11+
@testset verbose=true showtiming=true "mesh_basic" begin
12+
points = mesh_basic(coordinates_min, coordinates_max, 2, 2)
13+
@test points == [0.0 1.0 0.0 0.5 1.0; 0.0 0.0 1.0 1.0 1.0]
14+
end
15+
@testset verbose=true showtiming=true "mesh_bisected_rectangle" begin
16+
points = mesh_bisected_rectangle(coordinates_min, coordinates_max, 1, 1)
17+
@test points == [0.0 1.0 0.0 1.0; 0.0 0.0 1.0 1.0]
18+
end
19+
end
20+
821
@testset verbose=true showtiming=true "build_delaunay_triangulation" begin
922
data_points = collect([0.0 0.0
1023
1.0 0.0

0 commit comments

Comments
 (0)