Make max_points_per_cell work with templates (#76)

* Make `max_points_per_cell` work with templates

* Reformat code

---------

Co-authored-by: Niklas Neher <[email protected]>

Author: efaulhaber

src/cell_lists/full_grid.jl

@@ -48,12 +48,12 @@ function FullGridCellList(; min_corner, max_corner, search_radius = 0.0,
     # Pad domain to avoid 0 in cell indices due to rounding errors.
     # We can't just use `eps()`, as one might use lower precision types.
     # This padding is safe, and will give us one more layer of cells in the worst case.
-    min_corner = SVector(Tuple(min_corner .- 1e-3 * search_radius))
-    max_corner = SVector(Tuple(max_corner .+ 1e-3 * search_radius))
+    min_corner = SVector(Tuple(min_corner .- 1.0f-3 * search_radius))
+    max_corner = SVector(Tuple(max_corner .+ 1.0f-3 * search_radius))
 
     if search_radius < eps()
         # Create an empty "template" cell list to be used with `copy_cell_list`
-        cells = construct_backend(backend, 0, 0)
+        cells = construct_backend(backend, 0, max_points_per_cell)
         linear_indices = nothing
     else
         # Note that we don't shift everything so that the first cell starts at `min_corner`.
@@ -180,5 +180,15 @@ function copy_cell_list(cell_list::FullGridCellList, search_radius, periodic_box
     (; min_corner, max_corner) = cell_list
 
     return FullGridCellList(; min_corner, max_corner, search_radius,
-                            backend = typeof(cell_list.cells))
+                            backend = typeof(cell_list.cells),
+                            max_points_per_cell = max_points_per_cell(cell_list.cells))
+end
+
+function max_points_per_cell(cells::DynamicVectorOfVectors)
+    return size(cells.backend, 1)
+end
+
+# Fallback when backend is a `Vector{Vector{T}}`
+function max_points_per_cell(cells)
+    return 100
 end

src/nhs_grid.jl

@@ -167,6 +167,12 @@ function initialize_grid!(neighborhood_search::GridNeighborhoodSearch, y::Abstra
 
     empty!(cell_list)
 
+    if neighborhood_search.search_radius < eps()
+        # Cannot initialize with zero search radius.
+        # This is used in TrixiParticles when a neighborhood search is not used.
+        return neighborhood_search
+    end
+
     for point in axes(y, 2)
         # Get cell index of the point's cell
         point_coords = extract_svector(y, Val(ndims(neighborhood_search)), point)

test/nhs_grid.jl

@@ -44,11 +44,13 @@
 
         # Check that the update strategy is preserved
         nhs = GridNeighborhoodSearch{2}(cell_list = FullGridCellList(; min_corner,
-                                                                     max_corner),
+                                                                     max_corner,
+                                                                     max_points_per_cell = 101),
                                         update_strategy = SerialUpdate())
         copy = copy_neighborhood_search(nhs, 1.0, 10)
 
         @test copy.update_strategy == SerialUpdate()
+        @test size(copy.cell_list.cells.backend, 1) == 101
     end
 
     @testset "Cells at Coordinate Limits" begin