From 6f26a6357c55c94885d7093784abcb74c360acc7 Mon Sep 17 00:00:00 2001
From: "Vasco C. B. Ferreira" <vascocbf@gmail.com>
Date: Wed, 8 Jul 2026 12:07:51 +0200
Subject: [PATCH] upgrade to cells; cells were not sorted by lower[0], now yes

---
 nufi/cells.h | 29 +++++++++++++++++++----------
 1 file changed, 19 insertions(+), 10 deletions(-)

diff --git a/nufi/cells.h b/nufi/cells.h
index a6f852c..3369adf 100644
--- a/nufi/cells.h
+++ b/nufi/cells.h
@@ -15,8 +15,8 @@ template <int dim> struct CellInfo {
   typename DoFHandler<dim>::active_cell_iterator cell;
   // usefull for locator
   Point<dim> lower;
-  // Point<dim> upper;
-  // double h;
+  Point<dim> upper;
+  double h;
 };
 
 template <int dim> class CellLocator {
@@ -44,11 +44,20 @@ void CellLocator<dim>::rebuild(const DoFHandler<dim> &dof_handler,
 
     info.cell = cell;
     info.lower = cell->vertex(0);
-    // info.upper = cell->vertex(GeometryInfo<dim>::vertices_per_cell - 1);
-    // info.h = info.upper[0] - info.lower[0];
+    info.upper = cell->vertex(GeometryInfo<dim>::vertices_per_cell - 1);
+
+    if (info.lower[0] > info.upper[0])
+      std::swap(info.lower, info.upper);
+
+    info.h = info.upper[0] - info.lower[0];
 
     cells.push_back(info);
   }
+
+  std::sort(cells.begin(), cells.end(),
+            [](const CellInfo<dim> &a, const CellInfo<dim> &b) {
+              return a.lower[0] < b.lower[0];
+            });
 }
 
 template <int dim>
@@ -62,20 +71,20 @@ CellLocator<dim>::locate(const Point<dim> &p) const {
 
   const double x = p[0];
 
-  // const double x_min = cells.front().lower[0];
-  // const double x_max = cells.back().uppper[0];
-  // const double L = x_max - x_min;
-
-  // binary search loop
   auto it = std::upper_bound(cells.begin(), cells.end(), x,
                              [](double value, const CellInfo<dim> &cell) {
                                return value < cell.lower[0];
-                             });
+                             }); // returns cell to the right of cell with x
+
   if (it == cells.begin())
     it = cells.begin();
   else
     --it;
 
+  // Safety check: make sure the point is really inside this cell
+  AssertThrow(x >= it->lower[0] - 1e-12 && x <= it->upper[0] + 1e-12,
+              ExcMessage("CellLocator failed to find containing cell."));
+
   return it->cell;
 }