From f1a3dd4d93a6412bb7234c87bfb11095c9b3b39f Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 12 Nov 2024 16:13:12 +0100
Subject: [PATCH] Make `max_points_per_cell` work with templates

---
 src/cell_lists/full_grid.jl | 18 ++++++++++++++----
 src/nhs_grid.jl             |  6 ++++++
 test/nhs_grid.jl            |  4 +++-
 3 files changed, 23 insertions(+), 5 deletions(-)

diff --git a/src/cell_lists/full_grid.jl b/src/cell_lists/full_grid.jl
index 4e8603d..dc8f77f 100644
--- a/src/cell_lists/full_grid.jl
+++ b/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 .- 1f-3 * search_radius))
+    max_corner = SVector(Tuple(max_corner .+ 1f-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
diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl
index fef556b..8bb1440 100644
--- a/src/nhs_grid.jl
+++ b/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)
diff --git a/test/nhs_grid.jl b/test/nhs_grid.jl
index 90a7a69..61b4c8a 100644
--- a/test/nhs_grid.jl
+++ b/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