CarpetRegrid2: Rely on CoordBase information

This thorn works now.

darcs-hash:20060420025155-dae7b-01881883f067668faf7686ed39bbd1e6af077c79.gz

3 files changed, 189 insertions, 31 deletions

diff --git a/Carpet/CarpetRegrid2/interface.ccl b/Carpet/CarpetRegrid2/interface.ccl
index 2ffb68d22..a624bc1a1 100644
--- a/Carpet/CarpetRegrid2/interface.ccl
+++ b/Carpet/CarpetRegrid2/interface.ccl
@@ -12,15 +12,51 @@ USES INCLUDE HEADER: carpet.hh
 
 
 
+# The location of the boundary points
+CCTK_INT FUNCTION GetBoundarySpecification \
+    (CCTK_INT IN        size,              \
+     CCTK_INT OUT ARRAY nboundaryzones,    \
+     CCTK_INT OUT ARRAY is_internal,       \
+     CCTK_INT OUT ARRAY is_staggered,      \
+     CCTK_INT OUT ARRAY shiftout)
+USES FUNCTION GetBoundarySpecification
+
+# The overall size of the domain
+CCTK_INT FUNCTION GetDomainSpecification \
+    (CCTK_INT  IN        size,           \
+     CCTK_REAL OUT ARRAY physical_min,   \
+     CCTK_REAL OUT ARRAY physical_max,   \
+     CCTK_REAL OUT ARRAY interior_min,   \
+     CCTK_REAL OUT ARRAY interior_max,   \
+     CCTK_REAL OUT ARRAY exterior_min,   \
+     CCTK_REAL OUT ARRAY exterior_max,   \
+     CCTK_REAL OUT ARRAY spacing)
+USES FUNCTION GetDomainSpecification
+
+# Conversion between boundary types
+CCTK_INT FUNCTION ConvertFromPhysicalBoundary \
+    (CCTK_INT  IN        size,                \
+     CCTK_REAL IN  ARRAY physical_min,        \
+     CCTK_REAL IN  ARRAY physical_max,        \
+     CCTK_REAL OUT ARRAY interior_min,        \
+     CCTK_REAL OUT ARRAY interior_max,        \
+     CCTK_REAL OUT ARRAY exterior_min,        \
+     CCTK_REAL OUT ARRAY exterior_max,        \
+     CCTK_REAL IN  ARRAY spacing)
+USES FUNCTION ConvertFromPhysicalBoundary
+
+
+
 # The true prototype of the routine below:
 #   int Carpet_Regrid (const cGH * cctkGH,
 #                      gh::mexts  * bbsss,
 #                      gh::rbnds  * obss,
 #                      gh::rprocs * pss,
 #                      int          force);
-CCTK_INT FUNCTION Carpet_Regrid (CCTK_POINTER_TO_CONST IN cctkGH, \
-                                 CCTK_POINTER IN bbsss,           \
-                                 CCTK_POINTER IN obss,            \
-                                 CCTK_POINTER IN pss,		  \
-				 CCTK_INT IN force)
+CCTK_INT FUNCTION Carpet_Regrid       \
+    (CCTK_POINTER_TO_CONST IN cctkGH, \
+     CCTK_POINTER          IN bbsss,  \
+     CCTK_POINTER          IN obss,   \
+     CCTK_POINTER          IN pss,    \
+     CCTK_INT IN force)
 PROVIDES FUNCTION Carpet_Regrid WITH CarpetRegrid2_Regrid LANGUAGE C
diff --git a/Carpet/CarpetRegrid2/par/qc0.par b/Carpet/CarpetRegrid2/par/qc0.par
index bf57c5259..3cdd69eff 100644
--- a/Carpet/CarpetRegrid2/par/qc0.par
+++ b/Carpet/CarpetRegrid2/par/qc0.par
@@ -164,6 +164,8 @@ TwoPunctures::par_m_minus    =  0.45
 TwoPunctures::par_P_plus [1] =  0.333
 TwoPunctures::par_P_minus[1] = -0.333
 
+TwoPunctures::grid_setup_method = "evaluation"
+
 TwoPunctures::TP_epsilon = 1e-4
 
 TwoPunctures::verbose = yes
diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index 014101706..29cfa0999 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -1,6 +1,7 @@
 #include <cassert>
 #include <cmath>
 #include <cstdlib>
+#include <iostream>
 #include <vector>
 
 #include "cctk.h"
@@ -24,6 +25,7 @@ namespace CarpetRegrid2 {
   
   typedef bboxset <int, dim> ibboxset;
   typedef vect <CCTK_REAL, dim> rvect;
+  typedef bbox <CCTK_REAL, dim> rbbox;
   
   
   
@@ -90,6 +92,40 @@ namespace CarpetRegrid2 {
   
   
   
+  ivect
+  rpos2ipos (rvect const & rpos,
+             rvect const & origin, rvect const & scale,
+             gh const & hh, int const rl)
+  {
+    ivect const levfac = hh.reffacts.at(rl);
+    assert (all (hh.baseextent.stride() % levfac == 0));
+    ivect const istride = hh.baseextent.stride() / levfac;
+    
+    return (ivect (floor ((rpos - origin) * scale / rvect(istride) + 0.5)) *
+            istride);
+  }
+  
+  rvect
+  ipos2rpos (ivect const & ipos,
+             rvect const & origin, rvect const & scale,
+             gh const & hh, int const rl)
+  {
+    return rvect(ipos) / scale + origin;
+  }
+  
+  rbbox
+  ibbox2rbbox (ibbox const & ib,
+               rvect const & origin, rvect const & scale,
+               gh const & hh, int const rl)
+  {
+    rvect const zero (0);
+    return rbbox (ipos2rpos (ib.lower() , origin, scale, hh, rl),
+                  ipos2rpos (ib.upper() , origin, scale, hh, rl),
+                  ipos2rpos (ib.stride(), zero  , scale, hh, rl));
+  }
+  
+  
+  
   extern "C" {
     CCTK_INT
     CarpetRegrid2_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
@@ -139,18 +175,59 @@ namespace CarpetRegrid2 {
     if (do_recompose) {
       
       // Find extent of domain
-      rvect global_lower, global_upper;
-      for (int d = 0; d < dim; ++ d) {
-        int const ierr = CCTK_CoordRange
-          (cctkGH, & global_lower[d], & global_upper[d], d + 1, 0, "cart3d");
-        if (ierr < 0) {
-          global_lower[d] = 0.0;
-          global_upper[d] = 1.0;
-        }
+      
+      // This requires that CoordBase is used (but this is not
+      // checked)
+      
+      iivect nboundaryzones;
+      iivect is_internal;
+      iivect is_staggered;
+      iivect shiftout;
+      {
+        CCTK_INT const ierr = GetBoundarySpecification
+          (2*dim,
+           & nboundaryzones[0][0],
+           & is_internal[0][0],
+           & is_staggered[0][0],
+           & shiftout[0][0]);
+        assert (not ierr);
+      }
+      
+      rvect physical_lower, physical_upper;
+      rvect interior_lower, interior_upper;
+      rvect exterior_lower, exterior_upper;
+      rvect spacing;
+      {
+        CCTK_INT const ierr = GetDomainSpecification
+          (dim,
+           & physical_lower[0], & physical_upper[0],
+           & interior_lower[0], & interior_upper[0],
+           & exterior_lower[0], & exterior_upper[0],
+           & spacing[0]);
+        assert (not ierr);
+      }
+      
+      // Adapt spacing for convergence level
+      spacing *= ipow ((CCTK_REAL) Carpet::mgfact, Carpet::basemglevel);
+      
+      {
+        CCTK_INT const ierr =
+          ConvertFromPhysicalBoundary
+          (dim,
+           & physical_lower[0], & physical_upper[0],
+           & interior_lower[0], & interior_upper[0],
+           & exterior_lower[0], & exterior_upper[0],
+           & spacing[0]);
+        assert (not ierr);
       }
       
-      ivect const global_extent (hh.baseextent.upper() - hh.baseextent.lower());
-      rvect const scale (rvect (global_extent) / (global_upper - global_lower));
+      rvect const origin (exterior_lower);
+      rvect const scale (rvect (hh.baseextent.stride()) / spacing);
+      
+      ivect const physical_ilower =
+        rpos2ipos (physical_lower, origin, scale, hh, 0);
+      ivect const physical_iupper =
+        rpos2ipos (physical_upper, origin, scale, hh, 0);
       
       // The set of refined regions
       vector <ibboxset> regions;
@@ -167,17 +244,15 @@ namespace CarpetRegrid2 {
           rvect const rmax = centre.position + centre.radius.at(rl);
           
           // Convert to an integer bbox
+          ivect const imin =
+            rpos2ipos (rmin, origin, scale, hh, rl);
+          ivect const imax =
+            rpos2ipos (rmax, origin, scale, hh, rl);
+          
           ivect const levfac = hh.reffacts.at(rl);
           assert (all (hh.baseextent.stride() % levfac == 0));
           ivect const istride = hh.baseextent.stride() / levfac;
           
-          ivect const imin =
-            (ivect (floor ((rmin - global_lower) * scale / rvect(istride)))
-             * istride);
-          ivect const imax =
-            (ivect (ceil  ((rmax - global_lower) * scale / rvect(istride)))
-             * istride);
-          
           ibbox const region (imin, imax, istride);
           
           // Add this region to the list of regions
@@ -192,7 +267,7 @@ namespace CarpetRegrid2 {
       for (size_t rl = regions.size() - 1; rl >= 2; -- rl) {
         ibbox coarse = * regions.at(rl-1).begin();
         
-        regions.at(rl).normalize();        
+        regions.at(rl).normalize();
         ibboxset coarsified;
         for (ibboxset::const_iterator ibb = regions.at(rl).begin();
              ibb != regions.at(rl).end();
@@ -201,17 +276,14 @@ namespace CarpetRegrid2 {
           ivect const distance (min_distance);
           ibbox const fbb = * ibb;
           ibbox const cbb = fbb.expanded_for(coarse);
-          ibbox const ebb = cbb.expand(distance, distance);
+          ibbox const ebb = cbb.expand (distance, distance);
           coarsified |= ebb;
         }
         
         regions.at(rl-1) |= coarsified;
       }
       
-      // Simplify the set of refined regions
-      for (size_t rl = 1; rl < regions.size(); ++ rl) {
-        regions.at(rl).normalize();
-      }
+      
       
       // Convert to (bbsss, obss, pss) triplet
       vector <vector <ibbox> > bbss;
@@ -224,6 +296,51 @@ namespace CarpetRegrid2 {
       
       for (size_t rl = 1; rl < regions.size(); ++ rl) {
         
+        rvect const level_physical_lower = physical_lower;
+        rvect const level_physical_upper = physical_upper;
+        rvect const level_spacing = spacing / hh.reffacts.at(rl);
+        rvect level_interior_lower, level_interior_upper;
+        rvect level_exterior_lower, level_exterior_upper;
+        {
+          CCTK_INT const ierr =
+            ConvertFromPhysicalBoundary
+            (dim,
+             & level_physical_lower[0], & level_physical_upper[0],
+             & level_interior_lower[0], & level_interior_upper[0],
+             & level_exterior_lower[0], & level_exterior_upper[0],
+             & level_spacing[0]);
+          assert (not ierr);
+        }
+        
+        ivect const level_exterior_ilower =
+          rpos2ipos (level_exterior_lower, origin, scale, hh, rl);
+        ivect const level_exterior_iupper =
+          rpos2ipos (level_exterior_upper, origin, scale, hh, rl);
+        
+        // Clip at the outer boundary
+        regions.at(rl).normalize();
+        ibboxset clipped;
+        for (ibboxset::const_iterator ibb = regions.at(rl).begin();
+             ibb != regions.at(rl).end();
+             ++ ibb)
+        {
+          ibbox const bb = * ibb;
+          
+          bvect const lower_is_outer = bb.lower() <= physical_ilower;
+          bvect const upper_is_outer = bb.upper() >= physical_iupper;
+          
+          ibbox const clipped_bb (max (bb.lower(), level_exterior_ilower),
+                                  min (bb.upper(), level_exterior_iupper),
+                                  bb.stride());
+          
+          clipped += clipped_bb;
+        }
+        
+        regions.at(rl) = clipped;
+        
+        // Simplify the set of refined regions
+        regions.at(rl).normalize();
+        
         // Create a vector of bboxes for this level
         gh::cexts bbs;
         gh::cbnds obs;
@@ -234,12 +351,15 @@ namespace CarpetRegrid2 {
              ++ ibb)
         {
           ibbox bb = * ibb;
+          
+          bvect const lower_is_outer = bb.lower() <= physical_ilower;
+          bvect const upper_is_outer = bb.upper() >= physical_iupper;
+          
           bbvect ob;
           for (int d = 0; d < dim; ++ d) {
-            ob[d][0] = bb.lower()[d] <= hh.baseextent.lower()[d];
-            ob[d][1] = bb.upper()[d] >= hh.baseextent.upper()[d];
+            ob[d][0] = lower_is_outer[d];
+            ob[d][1] = upper_is_outer[d];
           }
-          bb = bb & hh.baseextent.expanded_for (bb);
           bbs.push_back (bb);
           obs.push_back (ob);
         }