diff --git a/src/cell_lists/full_grid.jl b/src/cell_lists/full_grid.jl index 1edd165e..c1702df6 100644 --- a/src/cell_lists/full_grid.jl +++ b/src/cell_lists/full_grid.jl @@ -81,33 +81,11 @@ function FullGridCellList(; min_corner, max_corner, return FullGridCellList(cells, linear_indices, min_corner, max_corner) end -@inline function cell_coords(coords, periodic_box::Nothing, cell_list::FullGridCellList, - cell_size) - (; min_corner) = cell_list - - # Subtract `min_corner` to offset coordinates so that the min corner of the grid - # corresponds to the (1, 1, 1) cell. - return Tuple(floor_to_int.((coords .- min_corner) ./ cell_size)) .+ 1 -end - -@inline function cell_coords(coords, periodic_box::PeriodicBox, cell_list::FullGridCellList, - cell_size) - # Subtract `periodic_box.min_corner` to offset coordinates so that the min corner - # of the grid corresponds to the (0, 0, 0) cell. - offset_coords = periodic_coords(coords, periodic_box) .- periodic_box.min_corner - - # Add 2, so that the min corner will be the (2, 2, 2)-cell. - # With this, we still have one padding layer in each direction around the periodic box, - # just like without using a periodic box. - # This is not needed for finding neighbor cells, but to make the bounds check - # work the same way as without a periodic box. - return Tuple(floor_to_int.(offset_coords ./ cell_size)) .+ 2 -end - -@inline function periodic_cell_index(cell_index, ::PeriodicBox, n_cells, - cell_list::FullGridCellList) - # 2-based modulo to match the indexing of the periodic box explained above. - return mod.(cell_index .- 2, n_cells) .+ 2 +@inline function nonperiodic_cell_coords(coords, cell_list::FullGridCellList, cell_size) + # Subtract `min_corner` to offset coordinates so that the (padded) min corner + # of the grid corresponds to the (1, 1, 1) cell. + # The unpadded min corner (that is passed by the user) corresponds to the (2, 2, 2) cell. + return Tuple(floor_to_int.((coords .- cell_list.min_corner) ./ cell_size) .+ 1) end function Base.empty!(cell_list::FullGridCellList) diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 8e3c1ff2..4ebcd203 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -586,37 +586,44 @@ end end @inline function periodic_cell_index(cell_index, neighborhood_search) - (; n_cells, periodic_box, cell_list) = neighborhood_search + (; n_cells, periodic_box) = neighborhood_search - periodic_cell_index(cell_index, periodic_box, n_cells, cell_list) + periodic_cell_index(cell_index, periodic_box, n_cells) end -@inline periodic_cell_index(cell_index, ::Nothing, n_cells, cell_list) = cell_index - -@inline function periodic_cell_index(cell_index, ::PeriodicBox, n_cells, cell_list) - # 1-based modulo - return mod1.(cell_index, n_cells) +@inline periodic_cell_index(cell_index, ::Nothing, n_cells) = cell_index + +@inline function periodic_cell_index(cell_index, ::PeriodicBox, n_cells) + # With the `FullGridCellList`, nonperiodic cell coords are in the range 2:n_cells+1, + # so we need an offset of 2 to preserve this, so that the min corner will still + # be the (2, 2, 2)-cell. + # With this, we still have one padding layer in each direction around the periodic box, + # just like without using a periodic box. + # This padding is not needed for finding neighbor cells, but to make the bounds check + # work the same way as without a periodic box (cells in bounds are 2:end-1). + # + # With the `DictionaryCellList`, the grid is conceptually infinite, + # so we can use any offset for the modulo operation. + # The offset only determines to which part of R^n the periodic box is mapped. + return mod.(cell_index .- 2, n_cells) .+ 2 end @inline function cell_coords(coords, neighborhood_search) - (; periodic_box, cell_list, cell_size) = neighborhood_search + # Compute cell coordinates ignoring periodicity at first + nonperiodic_cell_coords_ = nonperiodic_cell_coords(coords, neighborhood_search) - return cell_coords(coords, periodic_box, cell_list, cell_size) + # Now return (without periodicity) or convert to periodic cell index (with periodicity) + return periodic_cell_index(nonperiodic_cell_coords_, neighborhood_search) end -@inline function cell_coords(coords, periodic_box::Nothing, cell_list, cell_size) - return Tuple(floor_to_int.(coords ./ cell_size)) -end +@inline function nonperiodic_cell_coords(coords, neighborhood_search) + (; cell_list, cell_size) = neighborhood_search -@inline function cell_coords(coords, periodic_box::PeriodicBox, cell_list, cell_size) - # Subtract `min_corner` to offset coordinates so that the min corner of the periodic - # box corresponds to the (0, 0, 0) cell of the NHS. - # This way, there are no partial cells in the domain if the domain size is an integer - # multiple of the cell size (which is required, see the constructor). - offset_coords = periodic_coords(coords, periodic_box) .- periodic_box.min_corner + return nonperiodic_cell_coords(coords, cell_list, cell_size) +end - # Add one for 1-based indexing. The min corner will be the (1, 1, 1)-cell. - return Tuple(floor_to_int.(offset_coords ./ cell_size)) .+ 1 +@inline function nonperiodic_cell_coords(coords, cell_list, cell_size) + return Tuple(floor_to_int.(coords ./ cell_size)) end function copy_neighborhood_search(nhs::GridNeighborhoodSearch, search_radius, n_points; diff --git a/test/nhs_grid.jl b/test/nhs_grid.jl index 17ff2fa9..9caabf8d 100644 --- a/test/nhs_grid.jl +++ b/test/nhs_grid.jl @@ -62,12 +62,11 @@ coords2 = [NaN, 0] coords3 = [typemax(Int) + 1.0, -typemax(Int) - 1.0] - @test PointNeighbors.cell_coords(coords1, nothing, nothing, (1.0, 1.0)) == - (typemax(Int), typemin(Int)) - @test PointNeighbors.cell_coords(coords2, nothing, nothing, (1.0, 1.0)) == - (typemax(Int), 0) - @test PointNeighbors.cell_coords(coords3, nothing, nothing, (1.0, 1.0)) == - (typemax(Int), typemin(Int)) + nhs = GridNeighborhoodSearch{2}(search_radius = 1.0, n_points = 1) + + @test PointNeighbors.cell_coords(coords1, nhs) == (typemax(Int), typemin(Int)) + @test PointNeighbors.cell_coords(coords2, nhs) == (typemax(Int), 0) + @test PointNeighbors.cell_coords(coords3, nhs) == (typemax(Int), typemin(Int)) # The full grid cell list adds one to the coordinates to avoid zero-indexing. # This corner case is not relevant, as `typemax` coordinates will always be out of @@ -75,12 +74,12 @@ cell_list = FullGridCellList(min_corner = (0.0, 0.0), max_corner = (1.0, 1.0), search_radius = 1.0) - @test PointNeighbors.cell_coords(coords1, nothing, cell_list, (1.0, 1.0)) == - (typemax(Int), typemin(Int)) .+ 1 - @test PointNeighbors.cell_coords(coords2, nothing, cell_list, (1.0, 1.0)) == - (typemax(Int), 1) .+ 1 - @test PointNeighbors.cell_coords(coords3, nothing, cell_list, (1.0, 1.0)) == - (typemax(Int), typemin(Int)) .+ 1 + nhs = GridNeighborhoodSearch{2}(search_radius = 1.0, n_points = 1, + cell_list = cell_list) + + @test PointNeighbors.cell_coords(coords1, nhs) == (typemax(Int), typemin(Int)) .+ 1 + @test PointNeighbors.cell_coords(coords2, nhs) == (typemax(Int), 1) .+ 1 + @test PointNeighbors.cell_coords(coords3, nhs) == (typemax(Int), typemin(Int)) .+ 1 end @testset "Rectangular Point Cloud 2D" begin @@ -239,7 +238,7 @@ [-0.08 0.0 0.18 0.1 -0.08 -0.12 -0.05 -0.09 0.15 0.39], [-0.08 0.0 0.18 0.1 -0.08 - -0.12 -0.05 -0.09 0.15 0.42], + -0.12 -0.05 -0.09 0.15 0.42] .- [0.0017, 0.005], [-0.08 0.0 0.18 0.1 -0.08 -0.12 -0.05 -0.09 0.15 0.39 0.14 0.34 0.12 0.06 0.13] @@ -249,28 +248,44 @@ PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.2, 0.4]), # The `GridNeighborhoodSearch` is forced to round up the cell sizes in this test # to avoid split cells. - PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.205, 0.43]), + # In order for the points to have the same relative positions inside their cells + # with the two cell lists, we need to offset the coordinates and periodic box + # to make the grid aligned with the origin. + PeriodicBox(min_corner = [-0.1, -0.2] .- [0.0017, 0.005], + max_corner = [0.205, 0.43] .- [0.0017, 0.005]), PeriodicBox(min_corner = [-0.1, -0.2, 0.05], max_corner = [0.2, 0.4, 0.35]) ] @testset verbose=true "$(names[i])" for i in eachindex(names) coords = coordinates[i] - nhs = GridNeighborhoodSearch{size(coords, 1)}(search_radius = 0.1, - n_points = size(coords, 2), - periodic_box = periodic_boxes[i]) - - initialize_grid!(nhs, coords) - - neighbors = [sort(collect(PointNeighbors.eachneighbor(coords[:, i], nhs))) - for i in 1:5] - - # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells - @test neighbors[1] == [1, 2, 3, 5] - @test neighbors[2] == [1, 2, 3] - @test neighbors[3] == [1, 2, 3] - @test neighbors[4] == [4] - @test neighbors[5] == [1, 5] + # `DictionaryCellList` + nhs1 = GridNeighborhoodSearch{size(coords, 1)}(search_radius = 0.1, + n_points = size(coords, 2), + periodic_box = periodic_boxes[i]) + + # `FullGridCellList` + cell_list = FullGridCellList(; min_corner = periodic_boxes[i].min_corner, + max_corner = periodic_boxes[i].max_corner, + search_radius = 0.1) + nhs2 = GridNeighborhoodSearch{size(coords, 1)}(search_radius = 0.1, + n_points = size(coords, 2), + periodic_box = periodic_boxes[i], + cell_list = cell_list) + + for nhs in (nhs1, nhs2) + initialize_grid!(nhs, coords) + + neighbors = [sort(collect(PointNeighbors.eachneighbor(coords[:, i], nhs))) + for i in 1:5] + + # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells + @test neighbors[1] == [1, 2, 3, 5] + @test neighbors[2] == [1, 2, 3] + @test neighbors[3] == [1, 2, 3] + @test neighbors[4] == [4] + @test neighbors[5] == [1, 5] + end end @testset "Offset Domain Triggering Split Cells" begin