avoids initial GLOBAL_REFINMENT iterations of while(cell->has_children())

README.md

@@ -11,7 +11,9 @@ dimensions: 1x1v
 
 notes:
 
-- about 3 times slower than previous locator
+- Locator not optimized for 1d.
+- works for higher dimensions
+-
 
 status: Working, to be re-reviewed
 

nufi/cells.h

@@ -1,11 +1,15 @@
 #ifndef CELLS_H
 #define CELLS_H
 
+#include <array>
 #include <cmath>
 #include <deal.II/base/geometry_info.h>
 #include <deal.II/base/point.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/grid/tria.h>
+#include <vector>
+
+#include "nufi/parameters.h"
 
 using namespace dealii;
 
@@ -23,9 +27,16 @@ public:
 
 private:
   const DoFHandler<dim> *dof_handler_ptr = nullptr;
-  typename Triangulation<dim>::cell_iterator root;
   Point<dim> lower;
   Point<dim> upper;
+
+  // Cached base level = Parameters::GLOBAL_REFINEMENT.
+  unsigned int base_level = 0;
+  unsigned int base_n_per_axis = 1; // 2^base_level
+  // Flat lookup, indexed in the same bit-interleaved order that
+  // GeometryInfo<dim>::child_cell_from_point produces at each level,
+  // so base_cells[idx] can be found with pure bit math, no tree walk.
+  std::vector<typename Triangulation<dim>::cell_iterator> base_cells;
 };
 
 template <int dim>
@@ -33,14 +44,59 @@ void CellLocator<dim>::rebuild(const DoFHandler<dim> &dof_handler,
                                const Triangulation<dim> &triangulation) {
   dof_handler_ptr = &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."));
 
-  root = triangulation.begin(0);
+  typename Triangulation<dim>::cell_iterator root = triangulation.begin(0);
   lower = root->vertex(0);
   upper = root->vertex(GeometryInfo<dim>::vertices_per_cell - 1);
+
+  base_level = Parameters::GLOBAL_REFINEMENT;
+  base_n_per_axis = 1u << base_level;
+
+  // Walk down exactly base_level times ONCE per base cell to build the
+  // flat table. This costs O(2^(dim*base_level)) total at rebuild time
+  // (proportional to the number of base cells), not per locate() call.
+  const unsigned int n_base_cells = 1u << (dim * base_level);
+  base_cells.assign(n_base_cells, typename Triangulation<dim>::cell_iterator());
+
+  // Recursive-free BFS/DFS: descend from root, tracking the accumulated
+  // child-index bits per level to know where to store each level-G cell.
+  std::vector<typename Triangulation<dim>::cell_iterator> stack;
+  std::vector<unsigned int> index_stack;
+  std::vector<unsigned int> depth_stack;
+  stack.push_back(root);
+  index_stack.push_back(0);
+  depth_stack.push_back(0);
+
+  while (!stack.empty()) {
+    auto cell = stack.back();
+    unsigned int idx = index_stack.back();
+    unsigned int depth = depth_stack.back();
+    stack.pop_back();
+    index_stack.pop_back();
+    depth_stack.pop_back();
+
+    if (depth == base_level) {
+      base_cells[idx] = cell;
+      continue;
+    }
+
+    // At this point cell must have children, since refine_global(base_level)
+    // guarantees a fully uniform tree down to base_level.
+    AssertThrow(cell->has_children(),
+                ExcMessage("CellLocator: mesh is not uniformly refined to "
+                           "Parameters::GLOBAL_REFINEMENT; base-level cache "
+                           "cannot be built. Did you coarsen below the "
+                           "global refinement level?"));
+
+    const unsigned int n_children = GeometryInfo<dim>::max_children_per_cell;
+    for (unsigned int c = 0; c < n_children; ++c) {
+      stack.push_back(cell->child(c));
+      index_stack.push_back(idx * n_children + c);
+      depth_stack.push_back(depth + 1);
+    }
+  }
 }
 
 template <int dim>
@@ -65,9 +121,27 @@ CellLocation<dim> CellLocator<dim>::locate(const Point<dim> &p) const {
     xi[d] = std::min(std::max(xi[d], 0.0), 1.0);
   }
 
-  // Steps 3-5: descend the refinement tree using deal.II's own
+  // Step 3: O(1) jump to the base-level (GLOBAL_REFINEMENT) cell via index
-  // reference-cell child logic (branch-free, handles dim=1,2,3 uniformly).
+  // math, replacing what used to be `base_level` iterations of the
-  typename Triangulation<dim>::cell_iterator cell = root;
+  // has_children() loop. Must mirror the same bit convention used when
+  // building base_cells in rebuild() (child index accumulated as
+  // idx = idx*n_children + child_cell_from_point(xi) at each level).
+  unsigned int idx = 0;
+  Point<dim> xi_local = xi;
+  for (unsigned int l = 0; l < base_level; ++l) {
+    const unsigned int child_index =
+        GeometryInfo<dim>::child_cell_from_point(xi_local);
+    xi_local =
+        GeometryInfo<dim>::cell_to_child_coordinates(xi_local, child_index);
+    idx = idx * GeometryInfo<dim>::max_children_per_cell + child_index;
+  }
+
+  typename Triangulation<dim>::cell_iterator cell = base_cells[idx];
+  xi = xi_local;
+
+  // Step 4: continue descending only through ADAPTIVE refinement beyond
+  // the base level -- this loop now only runs `depth - base_level` times
+  // instead of `depth` times.
   while (cell->has_children()) {
     const unsigned int child_index =
         GeometryInfo<dim>::child_cell_from_point(xi);
@@ -75,7 +149,6 @@ CellLocation<dim> CellLocator<dim>::locate(const Point<dim> &p) const {
     cell = cell->child(child_index);
   }
 
-  // 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);

src/main.cc

@@ -16,7 +16,7 @@ void clear_results_directory(const std::string &dir) {
 }
 
 int main() {
-  omp_set_max_active_levels(1);
+  // omp_set_max_active_levels(1);
   std::cout << "Threads: " << omp_get_max_threads() << "\n";
   try {
     clear_results_directory("results");