1 files changed, 100 insertions, 0 deletions
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 9b4a5f443..9a04af1df 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -272,6 +272,106 @@ local_components (int const rl)
 
 
 
+// Find the refinement level and component to which a grid point
+// belongs.  This uses a tree search over the superregions in the grid
+// struction, which should scale reasonably (O(n log n)) better with
+// the number of componets components.
+void
+gh::
+locate_position (rvect const & rpos,
+                 int const ml,
+                 int const minrl, int const maxrl,
+                 int & rl, int & c, ivect & aligned_ipos) const
+{
+  assert (ml>=0 and ml<mglevels());
+  assert (minrl>=0 and minrl<=maxrl and maxrl<=reflevels());
+  
+  // Try finer levels first
+  for (rl = maxrl-1; rl >= minrl; --rl) {
+    
+    // Align (round) the position to the nearest existing grid point
+    // on this refinement level
+    ivect const str = baseextent(ml,rl).stride();
+    aligned_ipos = ivect(floor(rpos / rvect(str) + rvect(0.5))) * str;
+    
+    gh::cregs const & regs = superregions.AT(rl);
+    
+    // Search all superregions linearly.  Each superregion corresponds
+    // to a "refined region", and the number of superregions is thus
+    // presumably independent of the number of processors.
+    for (size_t r = 0; r < regs.size(); ++r) {
+      region_t const & reg = regs.AT(r);
+      if (reg.extent.contains(aligned_ipos)) {
+        // We found the superregion to which this grid point belongs
+        
+        // Search the superregion hierarchically
+        pseudoregion_t const * const preg =
+          reg.processors->search(aligned_ipos);
+        assert (preg);
+        
+        // We now know the refinement level, component, and index to
+        // which this grid point belongs
+        c = preg->component;
+        return;
+      }
+    }
+  } // for rl
+  
+  // The point does not belong to any component on any refinement
+  // level
+  rl = -1;
+  c = -1;
+}
+
+void
+gh::
+locate_position (ivect const & ipos,
+                 int const ml,
+                 int const minrl, int const maxrl,
+                 int & rl, int & c, ivect & aligned_ipos) const
+{
+  assert (ml>=0 and ml<mglevels());
+  assert (minrl>=0 and minrl<=maxrl and maxrl<=reflevels());
+  
+  // Try finer levels first
+  for (rl = maxrl-1; rl >= minrl; --rl) {
+    
+    // Align (round) the position to the nearest existing grid point
+    // on this refinement level
+    ivect const str = baseextent(ml, rl).stride();
+    aligned_ipos = ivect(floor(rvect(ipos) / rvect(str) + rvect(0.5))) * str;
+    
+    gh::cregs const & regs = superregions.AT(rl);
+    
+    // Search all superregions linearly.  Each superregion corresponds
+    // to a "refined region", and the number of superregions is thus
+    // presumably independent of the number of processors.
+    for (size_t r = 0; r < regs.size(); ++r) {
+      region_t const & reg = regs.AT(r);
+      if (reg.extent.contains(aligned_ipos)) {
+        // We found the superregion to which this grid point belongs
+        
+        // Search the superregion hierarchically
+        pseudoregion_t const * const preg =
+          reg.processors->search(aligned_ipos);
+        assert (preg);
+        
+        // We now know the refinement level, component, and index to
+        // which this grid point belongs
+        c = preg->component;
+        return;
+      }
+    }
+  } // for rl
+  
+  // The point does not belong to any component on any refinement
+  // level
+  rl = -1;
+  c = -1;
+}
+
+
+
 // Time hierarchy management
 
 void