|
| 1 | +@testset verbose=true "GridNeighborhoodSearch" begin |
| 2 | + @testset "Coordinate Limits" begin |
| 3 | + # Test the threshold for very large and very small coordinates. |
| 4 | + coords1 = [Inf, -Inf] |
| 5 | + coords2 = [NaN, 0] |
| 6 | + coords3 = [typemax(Int) + 1.0, -typemax(Int) - 1.0] |
| 7 | + |
| 8 | + @test TrixiNeighborhoodSearch.cell_coords(coords1, nothing, (1.0, 1.0)) == |
| 9 | + (typemax(Int), typemin(Int)) |
| 10 | + @test TrixiNeighborhoodSearch.cell_coords(coords2, nothing, (1.0, 1.0)) == |
| 11 | + (typemax(Int), 0) |
| 12 | + @test TrixiNeighborhoodSearch.cell_coords(coords3, nothing, (1.0, 1.0)) == |
| 13 | + (typemax(Int), typemin(Int)) |
| 14 | + end |
| 15 | + |
| 16 | + @testset "Rectangular Point Cloud 2D" begin |
| 17 | + #### Setup |
| 18 | + # Rectangular filled with equidistant spaced particles |
| 19 | + # from (x, y) = (-0.25, -0.25) to (x, y) = (0.35, 0.35) |
| 20 | + range = -0.25:0.1:0.35 |
| 21 | + coordinates1 = hcat(collect.(Iterators.product(range, range))...) |
| 22 | + nparticles = size(coordinates1, 2) |
| 23 | + |
| 24 | + particle_position1 = [0.05, 0.05] |
| 25 | + particle_spacing = 0.1 |
| 26 | + radius = particle_spacing |
| 27 | + |
| 28 | + # Create neighborhood search |
| 29 | + nhs1 = GridNeighborhoodSearch{2}(radius, nparticles) |
| 30 | + |
| 31 | + coords_fun(i) = coordinates1[:, i] |
| 32 | + initialize!(nhs1, coords_fun) |
| 33 | + |
| 34 | + # Get each neighbor for `particle_position1` |
| 35 | + neighbors1 = sort(collect(TrixiNeighborhoodSearch.eachneighbor(particle_position1, |
| 36 | + nhs1))) |
| 37 | + |
| 38 | + # Move particles |
| 39 | + coordinates2 = coordinates1 .+ [1.4, -3.5] |
| 40 | + |
| 41 | + # Update neighborhood search |
| 42 | + coords_fun2(i) = coordinates2[:, i] |
| 43 | + update!(nhs1, coords_fun2) |
| 44 | + |
| 45 | + # Get each neighbor for updated NHS |
| 46 | + neighbors2 = sort(collect(TrixiNeighborhoodSearch.eachneighbor(particle_position1, |
| 47 | + nhs1))) |
| 48 | + |
| 49 | + # Change position |
| 50 | + particle_position2 = particle_position1 .+ [1.4, -3.5] |
| 51 | + |
| 52 | + # Get each neighbor for `particle_position2` |
| 53 | + neighbors3 = sort(collect(TrixiNeighborhoodSearch.eachneighbor(particle_position2, |
| 54 | + nhs1))) |
| 55 | + |
| 56 | + # Double search radius |
| 57 | + nhs2 = GridNeighborhoodSearch{2}(2 * radius, size(coordinates1, 2)) |
| 58 | + initialize!(nhs2, coords_fun) |
| 59 | + |
| 60 | + # Get each neighbor in double search radius |
| 61 | + neighbors4 = sort(collect(TrixiNeighborhoodSearch.eachneighbor(particle_position1, |
| 62 | + nhs2))) |
| 63 | + |
| 64 | + # Move particles |
| 65 | + coordinates2 = coordinates1 .+ [0.4, -0.4] |
| 66 | + |
| 67 | + # Update neighborhood search |
| 68 | + update!(nhs2, coords_fun2) |
| 69 | + |
| 70 | + # Get each neighbor in double search radius |
| 71 | + neighbors5 = sort(collect(TrixiNeighborhoodSearch.eachneighbor(particle_position1, |
| 72 | + nhs2))) |
| 73 | + |
| 74 | + #### Verification |
| 75 | + @test neighbors1 == [17, 18, 19, 24, 25, 26, 31, 32, 33] |
| 76 | + |
| 77 | + @test neighbors2 == Int[] |
| 78 | + |
| 79 | + @test neighbors3 == [17, 18, 19, 24, 25, 26, 31, 32, 33] |
| 80 | + |
| 81 | + @test neighbors4 == [ |
| 82 | + 9, 10, 11, 12, 13, 14, 16, 17, 18, 19, 20, 21, 23, 24, 25, |
| 83 | + 26, 27, 28, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 41, 42, 44, 45, 46, 47, |
| 84 | + 48, 49] |
| 85 | + |
| 86 | + @test neighbors5 == [36, 37, 38, 43, 44, 45] |
| 87 | + end |
| 88 | + |
| 89 | + @testset "Rectangular Point Cloud 3D" begin |
| 90 | + #### Setup |
| 91 | + # Rectangular filled with equidistant spaced particles |
| 92 | + # from (x, y, z) = (-0.25, -0.25, -0.25) to (x, y) = (0.35, 0.35, 0.35) |
| 93 | + range = -0.25:0.1:0.35 |
| 94 | + coordinates1 = hcat(collect.(Iterators.product(range, range, range))...) |
| 95 | + nparticles = size(coordinates1, 2) |
| 96 | + |
| 97 | + particle_position1 = [0.05, 0.05, 0.05] |
| 98 | + particle_spacing = 0.1 |
| 99 | + radius = particle_spacing |
| 100 | + |
| 101 | + # Create neighborhood search |
| 102 | + nhs1 = GridNeighborhoodSearch{3}(radius, nparticles) |
| 103 | + |
| 104 | + coords_fun(i) = coordinates1[:, i] |
| 105 | + initialize!(nhs1, coords_fun) |
| 106 | + |
| 107 | + # Get each neighbor for `particle_position1` |
| 108 | + neighbors1 = sort(collect(TrixiNeighborhoodSearch.eachneighbor(particle_position1, |
| 109 | + nhs1))) |
| 110 | + |
| 111 | + # Move particles |
| 112 | + coordinates2 = coordinates1 .+ [1.4, -3.5, 0.8] |
| 113 | + |
| 114 | + # Update neighborhood search |
| 115 | + coords_fun2(i) = coordinates2[:, i] |
| 116 | + update!(nhs1, coords_fun2) |
| 117 | + |
| 118 | + # Get each neighbor for updated NHS |
| 119 | + neighbors2 = sort(collect(TrixiNeighborhoodSearch.eachneighbor(particle_position1, |
| 120 | + nhs1))) |
| 121 | + |
| 122 | + # Change position |
| 123 | + particle_position2 = particle_position1 .+ [1.4, -3.5, 0.8] |
| 124 | + |
| 125 | + # Get each neighbor for `particle_position2` |
| 126 | + neighbors3 = sort(collect(TrixiNeighborhoodSearch.eachneighbor(particle_position2, |
| 127 | + nhs1))) |
| 128 | + |
| 129 | + #### Verification |
| 130 | + @test neighbors1 == |
| 131 | + [115, 116, 117, 122, 123, 124, 129, 130, 131, 164, 165, 166, 171, 172, |
| 132 | + 173, 178, 179, 180, 213, 214, 215, 220, 221, 222, 227, 228, 229] |
| 133 | + |
| 134 | + @test neighbors2 == Int[] |
| 135 | + |
| 136 | + @test neighbors3 == |
| 137 | + [115, 116, 117, 122, 123, 124, 129, 130, 131, 164, 165, 166, 171, 172, |
| 138 | + 173, 178, 179, 180, 213, 214, 215, 220, 221, 222, 227, 228, 229] |
| 139 | + end |
| 140 | + |
| 141 | + @testset verbose=true "Periodicity 2D" begin |
| 142 | + @testset "Simple Example" begin |
| 143 | + coords = [-0.08 0.0 0.18 0.1 -0.08 |
| 144 | + -0.12 -0.05 -0.09 0.15 0.39] |
| 145 | + |
| 146 | + # 3 x 6 cells |
| 147 | + nhs = GridNeighborhoodSearch{2}(0.1, size(coords, 2), |
| 148 | + periodic_box_min_corner = [-0.1, -0.2], |
| 149 | + periodic_box_max_corner = [0.2, 0.4]) |
| 150 | + |
| 151 | + initialize!(nhs, coords) |
| 152 | + |
| 153 | + neighbors = [sort(collect(TrixiNeighborhoodSearch.eachneighbor(coords[:, i], |
| 154 | + nhs))) |
| 155 | + for i in 1:5] |
| 156 | + |
| 157 | + # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells |
| 158 | + @test neighbors[1] == [1, 2, 3, 5] |
| 159 | + @test neighbors[2] == [1, 2, 3] |
| 160 | + @test neighbors[3] == [1, 2, 3] |
| 161 | + @test neighbors[4] == [4] |
| 162 | + @test neighbors[5] == [1, 5] |
| 163 | + |
| 164 | + neighbors_loop = [Int[] for _ in axes(coords, 2)] |
| 165 | + |
| 166 | + for_particle_neighbor(coords, coords, nhs, |
| 167 | + particles = axes(coords, 2)) do particle, |
| 168 | + neighbor, |
| 169 | + pos_diff, |
| 170 | + distance |
| 171 | + append!(neighbors_loop[particle], neighbor) |
| 172 | + end |
| 173 | + |
| 174 | + @test sort(neighbors_loop[1]) == [1, 3, 5] |
| 175 | + @test sort(neighbors_loop[2]) == [2] |
| 176 | + @test sort(neighbors_loop[3]) == [1, 3] |
| 177 | + @test sort(neighbors_loop[4]) == [4] |
| 178 | + @test sort(neighbors_loop[5]) == [1, 5] |
| 179 | + end |
| 180 | + |
| 181 | + @testset "Rounding Up Cell Sizes" begin |
| 182 | + coords = [-0.08 0.0 0.18 0.1 -0.08 |
| 183 | + -0.12 -0.05 -0.09 0.15 0.42] |
| 184 | + |
| 185 | + # 3 x 6 cells |
| 186 | + nhs = GridNeighborhoodSearch{2}(0.1, size(coords, 2), |
| 187 | + periodic_box_min_corner = [-0.1, -0.2], |
| 188 | + periodic_box_max_corner = [0.205, 0.43]) |
| 189 | + |
| 190 | + initialize!(nhs, coords) |
| 191 | + |
| 192 | + neighbors = [sort(collect(TrixiNeighborhoodSearch.eachneighbor(coords[:, i], |
| 193 | + nhs))) |
| 194 | + for i in 1:5] |
| 195 | + |
| 196 | + # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells |
| 197 | + @test neighbors[1] == [1, 2, 3, 5] |
| 198 | + @test neighbors[2] == [1, 2, 3] |
| 199 | + @test neighbors[3] == [1, 2, 3] |
| 200 | + @test neighbors[4] == [4] |
| 201 | + @test neighbors[5] == [1, 5] |
| 202 | + |
| 203 | + neighbors_loop = [Int[] for _ in axes(coords, 2)] |
| 204 | + |
| 205 | + for_particle_neighbor(coords, coords, nhs, |
| 206 | + particles = axes(coords, 2)) do particle, |
| 207 | + neighbor, |
| 208 | + pos_diff, |
| 209 | + distance |
| 210 | + append!(neighbors_loop[particle], neighbor) |
| 211 | + end |
| 212 | + |
| 213 | + @test sort(neighbors_loop[1]) == [1, 3, 5] |
| 214 | + @test sort(neighbors_loop[2]) == [2] |
| 215 | + @test sort(neighbors_loop[3]) == [1, 3] |
| 216 | + @test sort(neighbors_loop[4]) == [4] |
| 217 | + @test sort(neighbors_loop[5]) == [1, 5] |
| 218 | + end |
| 219 | + |
| 220 | + @testset "Offset Domain Triggering Split Cells" begin |
| 221 | + # This used to trigger a "split cell bug", where the left and right boundary |
| 222 | + # cells were only partially contained in the domain. |
| 223 | + # The left particle was placed inside a ghost cells, which causes it to not |
| 224 | + # see the right particle, even though it is within the search distance. |
| 225 | + # The domain size is an integer multiple of the cell size, but the NHS did not |
| 226 | + # offset the grid based on the domain position. |
| 227 | + # See https://github.com/trixi-framework/TrixiNeighborhoodSearch.jl/pull/211 |
| 228 | + # for a more detailed explanation. |
| 229 | + coords = [-1.4 1.9 |
| 230 | + 0.0 0.0] |
| 231 | + |
| 232 | + # 5 x 1 cells |
| 233 | + nhs = GridNeighborhoodSearch{2}(1.0, size(coords, 2), |
| 234 | + periodic_box_min_corner = [-1.5, 0.0], |
| 235 | + periodic_box_max_corner = [2.5, 3.0]) |
| 236 | + |
| 237 | + initialize!(nhs, coords) |
| 238 | + |
| 239 | + neighbors = [sort(unique(collect(TrixiNeighborhoodSearch.eachneighbor(coords[:, |
| 240 | + i], |
| 241 | + nhs)))) |
| 242 | + for i in 1:2] |
| 243 | + |
| 244 | + @test neighbors[1] == [1, 2] |
| 245 | + @test neighbors[2] == [1, 2] |
| 246 | + end |
| 247 | + end |
| 248 | + |
| 249 | + @testset verbose=true "Periodicity 3D" begin |
| 250 | + coords = [-0.08 0.0 0.18 0.1 -0.08 |
| 251 | + -0.12 -0.05 -0.09 0.15 0.39 |
| 252 | + 0.14 0.34 0.12 0.06 0.13] |
| 253 | + |
| 254 | + # 3 x 6 x 3 cells |
| 255 | + nhs = GridNeighborhoodSearch{3}(0.1, size(coords, 2), |
| 256 | + periodic_box_min_corner = [-0.1, -0.2, 0.05], |
| 257 | + periodic_box_max_corner = [0.2, 0.4, 0.35]) |
| 258 | + |
| 259 | + initialize!(nhs, coords) |
| 260 | + |
| 261 | + neighbors = [sort(collect(TrixiNeighborhoodSearch.eachneighbor(coords[:, i], nhs))) |
| 262 | + for i in 1:5] |
| 263 | + |
| 264 | + # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells |
| 265 | + @test neighbors[1] == [1, 2, 3, 5] |
| 266 | + @test neighbors[2] == [1, 2, 3] |
| 267 | + @test neighbors[3] == [1, 2, 3] |
| 268 | + @test neighbors[4] == [4] |
| 269 | + @test neighbors[5] == [1, 5] |
| 270 | + |
| 271 | + neighbors_loop = [Int[] for _ in axes(coords, 2)] |
| 272 | + |
| 273 | + for_particle_neighbor(coords, coords, |
| 274 | + nhs) do particle, neighbor, pos_diff, distance |
| 275 | + append!(neighbors_loop[particle], neighbor) |
| 276 | + end |
| 277 | + |
| 278 | + @test sort(neighbors_loop[1]) == [1, 3, 5] |
| 279 | + @test sort(neighbors_loop[2]) == [2] |
| 280 | + @test sort(neighbors_loop[3]) == [1, 3] |
| 281 | + @test sort(neighbors_loop[4]) == [4] |
| 282 | + @test sort(neighbors_loop[5]) == [1, 5] |
| 283 | + end |
| 284 | +end |
0 commit comments