Handle empty grids.

Handle empty grids.
Fix global and level mode handling of the interpolator.

Allow local mode calls for the interpolator, reduction operator, and
hyperslabber, assuming that people know what they are doing.

darcs-hash:20030513101400-07bb3-42af7544c2ea98dbef5a7ad6d5112bea2dcc404f.gz

6 files changed, 55 insertions, 48 deletions

diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index 72b04c264..ed89334e0 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -19,7 +19,7 @@
 #include "carpet.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.43 2003/05/12 16:24:25 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.44 2003/05/13 12:14:00 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_Carpet_SetupGH_cc);
 }
 
@@ -370,12 +370,13 @@ namespace Carpet {
             ((int*)arrdata[group].info.bbox)[2*d+f] = obnds[d][f];
           }
           
-          assert (arrdata[group].info.lsh[d]>=0
-                  && arrdata[group].info.lsh[d]<=arrdata[group].info.gsh[d]);
-          assert (arrdata[group].info.lbnd[d]>=0
-                  && arrdata[group].info.lbnd[d]<=arrdata[group].info.ubnd[d]
-                  && arrdata[group].info.ubnd[d]<arrdata[group].info.gsh[d]);
-          assert (arrdata[group].info.ubnd[d]-arrdata[group].info.lbnd[d]+1 == arrdata[group].info.lsh[d]);
+          assert (arrdata[group].info.lsh[d]>=0);
+          assert (arrdata[group].info.lsh[d]<=arrdata[group].info.gsh[d]);
+          assert (arrdata[group].info.lbnd[d]>=0);
+          assert (arrdata[group].info.lbnd[d]<=arrdata[group].info.ubnd[d]+1);
+          assert (arrdata[group].info.ubnd[d]<arrdata[group].info.gsh[d]);
+          assert (arrdata[group].info.lbnd[d] + arrdata[group].info.lsh[d] - 1
+                  == arrdata[group].info.ubnd[d]);
           assert (arrdata[group].info.lbnd[d]<=arrdata[group].info.ubnd[d]+1);
         }
         
diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc
index 574213b1b..8c8313f1c 100644
--- a/Carpet/CarpetInterp/src/interp.cc
+++ b/Carpet/CarpetInterp/src/interp.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.4 2003/05/12 16:25:40 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.5 2003/05/13 12:14:00 schnetter Exp $
 
 #include <assert.h>
 #include <math.h>
@@ -19,7 +19,7 @@
 #include "interp.hh"
 
 extern "C" {
-  static char const * const rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.4 2003/05/12 16:25:40 schnetter Exp $";
+  static char const * const rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.5 2003/05/13 12:14:00 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetInterp_interp_cc);
 }
 
@@ -66,14 +66,6 @@ namespace CarpetInterp {
     
     
     
-    // We want to be in global or in level mode
-    if (hh->local_components(reflevel) != 1 && component != -1) {
-      CCTK_WARN (0, "Cannot interpolate in local mode");
-    }
-    if (hh->local_components(reflevel) != 1) assert (component == -1);
-    
-    
-    
     // Find out about the coordinates
     const char * coord_system_name
       = CCTK_CoordSystemName (coord_system_handle);
@@ -99,8 +91,8 @@ namespace CarpetInterp {
     int const nprocs = CCTK_nProcs (cgh);
     assert (myproc>=0 && myproc<nprocs);
     
-    int const minrl = reflevel==-1 ? 0 : reflevel;
-    int const maxrl = reflevel==-1 ? hh->reflevels()-1 : reflevel;
+    int const minrl = reflevel==-1 ? 0               : reflevel;
+    int const maxrl = reflevel==-1 ? hh->reflevels() : reflevel+1;
     int maxncomps = 0;
     for (int rl=minrl; rl<maxrl; ++rl) {
       maxncomps = max(maxncomps, hh->components(rl));
@@ -140,15 +132,14 @@ namespace CarpetInterp {
           if (hh->extents[reflevel][c][mglevel].contains(ipos)) {
             rlev[n] = rl;
             home[n] = c;
-            break;
+            goto found;
           }
         }
       }
-      if (! (rlev[n]>=minrl && rlev[n]<maxrl && home[n]>=0 && home[n]<hh->components(rlev[n]))) {
-        CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
-                    "Interpolation point #%d at [%g,%g,%g] is not on any grid patch",
-                    n, pos[0], pos[1], pos[2]);
-      }
+      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Interpolation point #%d at [%g,%g,%g] is not on any grid patch",
+                  n, pos[0], pos[1], pos[2]);
+    found:
       assert (rlev[n]>=minrl && rlev[n]<maxrl);
       assert (home[n]>=0 && home[n]<hh->components(rlev[n]));
       ++ homecnts [home[n] + (rlev[n]-minrl) * maxncomps];
@@ -236,9 +227,15 @@ namespace CarpetInterp {
     //
     int const saved_reflevel = reflevel;
     int const saved_mglevel = mglevel;
+    int const saved_component = component;
+    if (component!=-1) {
+      set_component (cgh, -1);
+    }
+    if (mglevel!=-1) {
+      set_mglevel (cgh, -1);
+    }
     if (reflevel!=-1) {
-      set_mglevel ((cGH*)(cgh), -1);
-      set_reflevel ((cGH*)(cgh), -1);
+      set_reflevel (cgh, -1);
     }
     BEGIN_REFLEVEL_LOOP(cgh) {
       if (reflevel>=minrl && reflevel<maxrl) {
@@ -334,8 +331,13 @@ namespace CarpetInterp {
       }
     } END_REFLEVEL_LOOP(cgh);
     if (saved_reflevel!=-1) {
-      set_reflevel ((cGH*)(cgh), saved_reflevel);
-      set_mglevel ((cGH*)(cgh), saved_mglevel);
+      set_reflevel (cgh, saved_reflevel);
+    }
+    if (saved_mglevel!=-1) {
+      set_mglevel (cgh, saved_mglevel);
+    }
+    if (saved_component!=-1) {
+      set_component (cgh, saved_component);
     }
     
     
diff --git a/Carpet/CarpetLib/src/bbox.cc b/Carpet/CarpetLib/src/bbox.cc
index f208c690d..0e11b9bfc 100644
--- a/Carpet/CarpetLib/src/bbox.cc
+++ b/Carpet/CarpetLib/src/bbox.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.14 2003/03/17 10:24:18 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.15 2003/05/13 12:14:00 schnetter Exp $
 
 #include <assert.h>
 
@@ -131,7 +131,8 @@ bool bbox<T,D>::is_aligned_with (const bbox& b) const {
 // Expand the bbox a little by multiples of the stride
 template<class T, int D>
 bbox<T,D> bbox<T,D>::expand (const vect<T,D>& lo, const vect<T,D>& hi) const {
-  assert (! empty());
+  // Allow expansion only into directions where the extent is not negative
+  assert (all(lower()<=upper() || (lo==0 && hi==0)));
   const vect<T,D> str = stride();
   const vect<T,D> lb = lower() - lo * str;
   const vect<T,D> ub = upper() + hi * str;
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index eb5bdd62f..d914b7206 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.33 2003/05/08 15:35:49 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.34 2003/05/13 12:14:00 schnetter Exp $
 
 #include <assert.h>
 
@@ -76,11 +76,7 @@ void dh<D>::recompose (const int initialise_upto) {
 	  if (h.outer_boundaries[rl][c][d][0]) ldist[d] = 0;
 	  if (h.outer_boundaries[rl][c][d][1]) udist[d] = 0;
 	}
-        if (! intr.empty()) {
-          // If a processor has zero grid points, then it gets an
-          // empty bbox, i.e. no ghost zones are added
-          boxes[rl][c][ml].exterior = intr.expand(ldist, udist);
-        }
+        boxes[rl][c][ml].exterior = intr.expand(ldist, udist);
 	
        	// Boundaries (ghost zones only)
        	// (interior + boundaries = exterior)
diff --git a/Carpet/CarpetReduce/src/reduce.cc b/Carpet/CarpetReduce/src/reduce.cc
index e05bf56d4..118ecf5be 100644
--- a/Carpet/CarpetReduce/src/reduce.cc
+++ b/Carpet/CarpetReduce/src/reduce.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.19 2003/05/12 16:25:28 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.20 2003/05/13 12:14:00 schnetter Exp $
 
 #include <assert.h>
 #include <float.h>
@@ -22,7 +22,7 @@
 #include "reduce.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.19 2003/05/12 16:25:28 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.20 2003/05/13 12:14:00 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetReduce_reduce_cc);
 }
 
@@ -667,7 +667,7 @@ namespace CarpetReduce {
     // TODO: allow all modes for grid scalars and grid arrays, and
     // restrict usage only for grid functions.
     for (int n=0; n<num_invars; ++n) {
-      if (CCTK_GroupTypeFromVarI(invars[n])) {
+      if (CCTK_GroupTypeFromVarI(invars[n]) != CCTK_GF) {
         CCTK_WARN (0, "Reduction operators for grid scalars and grid arrays are not yet implemented");
       }
     }
@@ -676,10 +676,6 @@ namespace CarpetReduce {
     if (reflevel == -1) {
       CCTK_WARN (0, "Reduction operators in global mode are not yet implemented");
     }
-    if (hh->local_components(reflevel) != 1 && component != -1) {
-      CCTK_WARN (0, "It is not possible to use a grid variable reduction operator in local mode");
-    }
-    if (hh->local_components(reflevel) != 1) assert (component == -1);
     
     assert (cgh);
     
@@ -715,6 +711,10 @@ namespace CarpetReduce {
     Initialise (cgh, proc, num_outvals, &myoutvals[0], outtype, &mycounts[0],
                 red);
     
+    int const saved_component = component;
+    if (component!=-1) {
+      set_component (cgh, -1);
+    }
     BEGIN_LOCAL_COMPONENT_LOOP(cgh) {
       
       int lsh[grpdim], bbox[2*grpdim], nghostzones[grpdim];
@@ -761,6 +761,9 @@ namespace CarpetReduce {
               &mycounts[0], red);
       
     } END_LOCAL_COMPONENT_LOOP(cgh);
+    if (saved_component!=-1) {
+      set_component (cgh, saved_component);
+    }
     
     Finalise (cgh, proc, num_outvals, outvals, outtype,
               &myoutvals[0], &mycounts[0],
diff --git a/Carpet/CarpetSlab/src/slab.cc b/Carpet/CarpetSlab/src/slab.cc
index 8f1091a8d..440f27f42 100644
--- a/Carpet/CarpetSlab/src/slab.cc
+++ b/Carpet/CarpetSlab/src/slab.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.cc,v 1.5 2003/05/12 16:25:11 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.cc,v 1.6 2003/05/13 12:14:00 schnetter Exp $
 
 #include <assert.h>
 #include <stdlib.h>
@@ -21,7 +21,7 @@
 #include "slab.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.cc,v 1.5 2003/05/12 16:25:11 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.cc,v 1.6 2003/05/13 12:14:00 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetSlab_slab_cc);
 }
 
@@ -53,7 +53,9 @@ namespace CarpetSlab {
     
     // Save global state
     int saved_component = component;
-    component = -1;
+    if (component!=-1) {
+      set_component (cgh, -1);
+    }
     
     // Check Cactus grid hierarchy
     assert (cgh);
@@ -235,7 +237,9 @@ namespace CarpetSlab {
     } // if components>0
     
     // Restore global state
-    component = saved_component;
+    if (saved_component!=-1) {
+      set_component (cgh, saved_component);
+    }
     
     // Success
     return hdata;