diff --git a/README.md b/README.md index b5151a7..c72f85f 100644 --- a/README.md +++ b/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 diff --git a/nufi/cells.h b/nufi/cells.h index de4102b..a09bac3 100644 --- a/nufi/cells.h +++ b/nufi/cells.h @@ -1,11 +1,15 @@ #ifndef CELLS_H #define CELLS_H +#include #include #include #include #include #include +#include + +#include "nufi/parameters.h" using namespace dealii; @@ -23,9 +27,16 @@ public: private: const DoFHandler *dof_handler_ptr = nullptr; - typename Triangulation::cell_iterator root; Point lower; Point 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::child_cell_from_point produces at each level, + // so base_cells[idx] can be found with pure bit math, no tree walk. + std::vector::cell_iterator> base_cells; }; template @@ -33,14 +44,59 @@ void CellLocator::rebuild(const DoFHandler &dof_handler, const Triangulation &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::cell_iterator root = triangulation.begin(0); lower = root->vertex(0); upper = root->vertex(GeometryInfo::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::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::cell_iterator> stack; + std::vector index_stack; + std::vector 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::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 @@ -65,9 +121,27 @@ CellLocation CellLocator::locate(const Point &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 - // reference-cell child logic (branch-free, handles dim=1,2,3 uniformly). - typename Triangulation::cell_iterator cell = root; + // Step 3: O(1) jump to the base-level (GLOBAL_REFINEMENT) cell via index + // math, replacing what used to be `base_level` iterations of the + // 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 xi_local = xi; + for (unsigned int l = 0; l < base_level; ++l) { + const unsigned int child_index = + GeometryInfo::child_cell_from_point(xi_local); + xi_local = + GeometryInfo::cell_to_child_coordinates(xi_local, child_index); + idx = idx * GeometryInfo::max_children_per_cell + child_index; + } + + typename Triangulation::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::child_cell_from_point(xi); @@ -75,7 +149,6 @@ CellLocation CellLocator::locate(const Point &p) const { cell = cell->child(child_index); } - // cell is now active (leaf) -> bind it to the DoFHandler. typename DoFHandler::active_cell_iterator dof_cell( &cell->get_triangulation(), cell->level(), cell->index(), dof_handler_ptr); diff --git a/src/main.cc b/src/main.cc index 11bf2a6..b6f9e42 100644 --- a/src/main.cc +++ b/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");