From 4d7078d56da605768e500cf2bdcd4ff939bebaf2 Mon Sep 17 00:00:00 2001
From: "Vasco C. B. Ferreira" <vascocbf@gmail.com>
Date: Fri, 31 Jul 2026 12:47:13 +0200
Subject: [PATCH] testing different locator with while(cell->has_children)

---
 nufi/cells.h           | 165 +++++++++++------------------------------
 nufi/grids.h           |   3 +-
 nufi/poisson_problem.h |   3 +-
 3 files changed, 45 insertions(+), 126 deletions(-)

diff --git a/nufi/cells.h b/nufi/cells.h
index 91d1312..de4102b 100644
--- a/nufi/cells.h
+++ b/nufi/cells.h
@@ -1,13 +1,7 @@
 #ifndef CELLS_H
 #define CELLS_H
 
-#include <algorithm>
-#include <array>
 #include <cmath>
-#include <limits>
-#include <vector>
-
-#include <deal.II/base/exceptions.h>
 #include <deal.II/base/geometry_info.h>
 #include <deal.II/base/point.h>
 #include <deal.II/dofs/dof_handler.h>
@@ -15,8 +9,6 @@
 
 using namespace dealii;
 
-// Current Locator needs Grid to be a single hypercube coarse cell,
-// uniformly refined to `base_level`, with isotropic local refinement on top.
 template <int dim> struct CellLocation {
   typename DoFHandler<dim>::active_cell_iterator cell;
   Point<dim> reference_point;
@@ -25,144 +17,73 @@ template <int dim> struct CellLocation {
 template <int dim> class CellLocator {
 public:
   void rebuild(const DoFHandler<dim> &dof_handler,
-               const Triangulation<dim> &triangulation,
-               unsigned int base_level);
+               const Triangulation<dim> &triangulation);
 
-  CellLocation<dim> locate(const Point<dim> &point) const;
+  CellLocation<dim> locate(const Point<dim> &p) const;
 
 private:
-  const DoFHandler<dim> *dof_handler_ = nullptr;
-
-  Point<dim> domain_lower_;
-  Point<dim> domain_upper_;
-  std::array<unsigned int, dim> base_n_{};
-  std::array<double, dim> base_dx_{};
-
-  // Flat, O(1)-indexable array of the base_level cells.
-  std::vector<typename DoFHandler<dim>::cell_iterator> base_cells_;
+  const DoFHandler<dim> *dof_handler_ptr = nullptr;
+  typename Triangulation<dim>::cell_iterator root;
+  Point<dim> lower;
+  Point<dim> upper;
 };
 
 template <int dim>
 void CellLocator<dim>::rebuild(const DoFHandler<dim> &dof_handler,
-                               const Triangulation<dim> &triangulation,
-                               unsigned int base_level) {
-  AssertThrow(&dof_handler.get_triangulation() == &triangulation,
-              ExcMessage("CellLocator::rebuild(): DoFHandler and Triangulation "
-                         "do not refer to the same mesh."));
+                               const Triangulation<dim> &triangulation) {
+  dof_handler_ptr = &dof_handler;
 
-  dof_handler_ = &dof_handler;
+  // Everything below assumes the mesh is a single hyper_cube coarse cell
+  // (true for your create_mesh(): GridGenerator::hyper_cube + refine_global).
+  AssertThrow(triangulation.n_cells(0) == 1,
+              ExcMessage("CellLocator assumes exactly one coarse/root cell."));
 
-  // --- bounding box of the (single) coarse cell ---
-  auto root = dof_handler.begin(0);
-  AssertThrow(root != dof_handler.end(0),
-              ExcMessage("CellLocator::rebuild(): no level-zero cell."));
-  {
-    auto after_root = root;
-    ++after_root;
-    AssertThrow(after_root == dof_handler.end(0),
-                ExcMessage("CellLocator currently requires exactly one "
-                           "coarse cell."));
-  }
-
-  for (unsigned int d = 0; d < dim; ++d) {
-    domain_lower_[d] = std::numeric_limits<double>::max();
-    domain_upper_[d] = std::numeric_limits<double>::lowest();
-  }
-  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
-    for (unsigned int d = 0; d < dim; ++d) {
-      domain_lower_[d] = std::min(domain_lower_[d], root->vertex(v)[d]);
-      domain_upper_[d] = std::max(domain_upper_[d], root->vertex(v)[d]);
-    }
-
-  // --- build the O(1)-indexable array of base_level cells ---
-  unsigned int total = 1;
-  for (unsigned int d = 0; d < dim; ++d) {
-    const double length = domain_upper_[d] - domain_lower_[d];
-    AssertThrow(length > 0.0,
-                ExcMessage("CellLocator::rebuild(): non-positive domain "
-                           "extent."));
-    base_n_[d] =
-        1u << base_level; // refine_global(base_level) -> 2^level per axis
-    base_dx_[d] = length / base_n_[d];
-    total *= base_n_[d];
-  }
-
-  base_cells_.assign(total, typename DoFHandler<dim>::cell_iterator());
-
-  for (auto cell = dof_handler.begin(base_level);
-       cell != dof_handler.end(base_level); ++cell) {
-    const auto bb = cell->bounding_box();
-    const auto c_lower = bb.get_boundary_points().first;
-
-    unsigned int idx = 0, stride = 1;
-    for (unsigned int d = 0; d < dim; ++d) {
-      unsigned int i = static_cast<unsigned int>(
-          std::round((c_lower[d] - domain_lower_[d]) / base_dx_[d]));
-      i = std::min(i, base_n_[d] - 1);
-      idx += i * stride;
-      stride *= base_n_[d];
-    }
-    base_cells_[idx] = cell;
-  }
-
-  for (const auto &c : base_cells_)
-    AssertThrow(c.state() == IteratorState::valid,
-                ExcMessage("CellLocator::rebuild(): failed to fill the base "
-                           "grid — mesh isn't uniformly refined to "
-                           "'base_level' as expected."));
+  root = triangulation.begin(0);
+  lower = root->vertex(0);
+  upper = root->vertex(GeometryInfo<dim>::vertices_per_cell - 1);
 }
 
 template <int dim>
-CellLocation<dim> CellLocator<dim>::locate(const Point<dim> &point) const {
-  AssertThrow(dof_handler_ != nullptr,
-              ExcMessage("CellLocator::locate(): rebuild() has not been "
-                         "called."));
+CellLocation<dim> CellLocator<dim>::locate(const Point<dim> &p) const {
+  AssertThrow(dof_handler_ptr != nullptr,
+              ExcMessage("CellLocator::rebuild() has not been called."));
 
-  std::array<unsigned int, dim> base_index;
-  Point<dim> reference_point; // local coords, updated at each descent step
-
-  // 1) periodic wrap + O(1) base-cell index per axis
+  // Step 1: periodic wrap into [lower, upper) per axis.
+  Point<dim> p_wrapped;
   for (unsigned int d = 0; d < dim; ++d) {
-    const double length = domain_upper_[d] - domain_lower_[d];
-    double shifted = point[d] - domain_lower_[d];
-    shifted -= length * std::floor(shifted / length);
-    if (shifted >= length)
-      shifted = 0.0;
-
-    const double xi = shifted / base_dx_[d];
-    const unsigned int i =
-        std::min(base_n_[d] - 1, static_cast<unsigned int>(std::floor(xi)));
-    base_index[d] = i;
-    reference_point[d] = xi - i; // fraction within the base cell
+    const double L = upper[d] - lower[d];
+    double x = p[d] - lower[d];
+    x = x - L * std::floor(x / L);
+    p_wrapped[d] = lower[d] + x;
   }
 
-  // 2) O(1) lookup of the base-level cell
-  unsigned int idx = 0, stride = 1;
+  // Step 2: reference coordinates in the root cell, clamped against
+  // floating-point drift at the domain boundary.
+  Point<dim> xi;
   for (unsigned int d = 0; d < dim; ++d) {
-    idx += base_index[d] * stride;
-    stride *= base_n_[d];
+    xi[d] = (p_wrapped[d] - lower[d]) / (upper[d] - lower[d]);
+    xi[d] = std::min(std::max(xi[d], 0.0), 1.0);
   }
-  auto cell = base_cells_[idx];
 
-  // 3) O(R_level_max) descent — only costs anything for cells that are
-  //    actually locally refined beyond base_level
+  // Steps 3-5: descend the refinement tree using deal.II's own
+  // reference-cell child logic (branch-free, handles dim=1,2,3 uniformly).
+  typename Triangulation<dim>::cell_iterator cell = root;
   while (cell->has_children()) {
-    AssertThrow(cell->n_children() == GeometryInfo<dim>::max_children_per_cell,
-                ExcMessage("CellLocator currently supports isotropic "
-                           "refinement only."));
     const unsigned int child_index =
-        GeometryInfo<dim>::child_cell_from_point(reference_point);
-    reference_point = GeometryInfo<dim>::cell_to_child_coordinates(
-        reference_point, child_index);
+        GeometryInfo<dim>::child_cell_from_point(xi);
+    xi = GeometryInfo<dim>::cell_to_child_coordinates(xi, child_index);
     cell = cell->child(child_index);
   }
 
-  AssertThrow(cell->is_active(),
-              ExcMessage("CellLocator tree traversal did not finish on an "
-                         "active cell."));
+  // cell is now active (leaf) -> bind it to the DoFHandler.
+  typename DoFHandler<dim>::active_cell_iterator dof_cell(
+      &cell->get_triangulation(), cell->level(), cell->index(),
+      dof_handler_ptr);
 
-  return CellLocation<dim>{typename DoFHandler<dim>::active_cell_iterator(cell),
-                           reference_point};
+  CellLocation<dim> location;
+  location.cell = dof_cell;
+  location.reference_point = xi;
+  return location;
 }
 
-#endif // CELLS_H
+#endif // !CELLS_H
diff --git a/nufi/grids.h b/nufi/grids.h
index 5540eac..65ce753 100644
--- a/nufi/grids.h
+++ b/nufi/grids.h
@@ -140,8 +140,7 @@ GridStructure<dim> make_grid_snapshot(PoissonProblem<dim> &poisson) {
   grid.constraints->close();
 
   grid.locator = std::make_unique<CellLocator<dim>>();
-  grid.locator->rebuild(*grid.dof_handler, *grid.triangulation,
-                        Parameters::GLOBAL_REFINEMENT);
+  grid.locator->rebuild(*grid.dof_handler, *grid.triangulation);
 
   // START: diagnostics
   AssertThrow(
diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h
index 332c309..0d15e69 100644
--- a/nufi/poisson_problem.h
+++ b/nufi/poisson_problem.h
@@ -338,8 +338,7 @@ template <int dim> void PoissonProblem<dim>::setup_system() {
   system_rhs.reinit(dof_handler.n_dofs());
 
   // used for evaluator to avoid running it anytime there is an eval
-  cell_locator.rebuild(dof_handler, triangulation,
-                       Parameters::GLOBAL_REFINEMENT);
+  cell_locator.rebuild(dof_handler, triangulation);
 }
 
 template <int dim> void PoissonProblem<dim>::assemble_system() {