From fec75e20b0964d901b5c09b124d54cd5e823c93e Mon Sep 17 00:00:00 2001
From: Erik Schnetter <schnetter@cct.lsu.edu>
Date: Mon, 14 Jul 2008 12:38:57 -0500
Subject: CarpetLib: Implement function gh::locate_position

Implement a new function gh::locate_position which finds the refinement
level and component which owns a certain location.
This is e.g. useful for interpolation.
---
 Carpet/CarpetLib/src/gh.cc | 100 +++++++++++++++++++++++++++++++++++++++++++++
 Carpet/CarpetLib/src/gh.hh |  10 +++++
 2 files changed, 110 insertions(+)

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
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index ffac38804..b80d71ca3 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -128,6 +128,16 @@ public:
   
   int local_components (const int rl) const;
   
+  void locate_position (rvect const & rpos,
+                        int const ml,
+                        int const minrl, int const maxrl,
+                        int & rl, int & c, ivect & aligned_ipos) const;
+  
+  void locate_position (ivect const & ipos,
+                        int const ml,
+                        int const minrl, int const maxrl,
+                        int & rl, int & c, ivect & aligned_ipos) const;
+  
   // Time hierarchy management
   void add (th * t);
   void remove (th * t);
-- 
cgit v1.2.3