CarpetLib: Suggest new functions ipos2rpos and rpos2ipos for gh

Implement (but leave commented out) new functions
  rvect ipos2rpos
  ivect rpos2ipos
  ivect rpos2ipos1
in the gh class.

2 files changed, 82 insertions, 0 deletions

diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 52ed1cf60..b281031a4 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -330,6 +330,75 @@ get_local_component (int const rl, int const c)
 {
   return local_components_.AT(rl).AT(c);
 }
+  
+  
+  
+#if 0
+// Convert an index location to a coordinate location
+rvect
+gh::
+ipos2rpos (ivect const & ipos,
+           rvect const & origin, rvect const & scale,
+           int const ml, int const rl)
+  const
+{
+  return rvect(ipos) / scale + origin;
+}
+
+
+
+// Convert a coordinate location to the nearest index location,
+// rounding downwards to break ties.  For cell centring, shift
+// upwards.
+ivect
+gh::
+rpos2ipos (rvect const & rpos,
+           rvect const & origin, rvect const & scale,
+           int const ml, int const rl)
+  const
+{
+  ivect const& istride = baseextents.at(ml).at(rl).stride();
+  
+  if (refcent == cell_centered) {
+    assert (all (istride % 2 == 0));
+  }
+  
+  rvect const spos = (rpos - origin) * scale;
+  ivect const ipos =
+    refcent == vertex_centered
+    ? ivect (floor (spos / rvect(istride) + rvect(0.5))) * istride
+    : ivect (floor (spos / rvect(istride)             )) * istride + istride/2;
+  
+  return ipos;
+}
+
+
+
+// Convert a coordinate location to the nearest index location,
+// rounding upwards to break ties.  For cell centring, shift downwards
+// instead of upwards.
+ivect
+gh::
+rpos2ipos1 (rvect const & rpos,
+            rvect const & origin, rvect const & scale,
+            int const ml, int const rl)
+  const
+{
+  ivect const& istride = baseextents.at(ml).at(rl).stride();
+  
+  if (refcent == cell_centered) {
+    assert (all (istride % 2 == 0));
+  }
+  
+  rvect const spos = (rpos - origin) * scale;
+  ivect const ipos =
+    refcent == vertex_centered
+    ? ivect (ceil (spos / rvect(istride) - rvect(0.5))) * istride
+    : ivect (ceil (spos / rvect(istride)             )) * istride - istride/2;
+  
+  return ipos;
+}
+#endif
 
 
 
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index 1942ae38b..1db0238c9 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -141,6 +141,19 @@ public:
   int get_component (int rl, int lc) const CCTK_ATTRIBUTE_PURE;
   int get_local_component (int rl, int c) const CCTK_ATTRIBUTE_PURE;
   
+#if 0
+  // Convert between index positions and coordinate positions
+  rvect ipos2rpos (ivect const & ipos,
+                   rvect const & origin, rvect const & scale,
+                   int const ml, int const rl) const;
+  ivect rpos2ipos (rvect const & rpos,
+                   rvect const & origin, rvect const & scale,
+                   int const ml, int const rl) const;
+  ivect rpos2ipos1 (rvect const & rpos,
+                    rvect const & origin, rvect const & scale,
+                    int const ml, int const rl) const;
+#endif
+  
   void locate_position (rvect const & rpos,
                         int const ml,
                         int const minrl, int const maxrl,