Added a bit of convenient functionality.

Added a bit of convenient functionality.
Fixed a few bugs that were uncovered by dynamic regridding.

darcs-hash:20011214153941-07bb3-c31c15cec651b5c3c992f529b7409c081fda5ba8.gz

9 files changed, 204 insertions, 139 deletions

diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc
index 5713f2d35..0d7522940 100644
--- a/Carpet/CarpetLib/src/bboxset.cc
+++ b/Carpet/CarpetLib/src/bboxset.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.8 2001/07/04 12:29:51 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.9 2001/12/14 16:39:41 schnetter Exp $
 
  ***************************************************************************/
 
@@ -95,7 +95,6 @@ T bboxset<T,D>::size () const {
 }
 
 
-
 // Add (bboxes that don't overlap)
 template<class T, int D>
 bboxset<T,D>& bboxset<T,D>::operator+= (const box& b) {
@@ -134,6 +133,16 @@ bboxset<T,D> bboxset<T,D>::operator+ (const bboxset& s) const {
   return r;
 }
 
+template<class T, int D>
+bboxset<T,D> bboxset<T,D>::plus (const bbox<T,D>& b1, const bbox<T,D>& b2) {
+  return bboxset(b1) + b2;
+}
+
+template<class T, int D>
+bboxset<T,D> bboxset<T,D>::plus (const bbox<T,D>& b, const bboxset<T,D>& s) {
+  return s + b;
+}
+
 
 
 // Union
@@ -293,6 +302,39 @@ bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b, const bboxset<T,D>& s) {
 
 
 
+// Equality
+template<class T, int D>
+bool bboxset<T,D>::operator<= (const bboxset<T,D>& s) const {
+  return (*this - s).empty();
+}
+
+template<class T, int D>
+bool bboxset<T,D>::operator< (const bboxset<T,D>& s) const {
+  return (*this - s).empty() && ! (s - *this).empty();
+}
+
+template<class T, int D>
+bool bboxset<T,D>::operator>= (const bboxset<T,D>& s) const {
+  return s <= *this;
+}
+
+template<class T, int D>
+bool bboxset<T,D>::operator> (const bboxset<T,D>& s) const {
+  return s < *this;
+}
+
+template<class T, int D>
+bool bboxset<T,D>::operator== (const bboxset<T,D>& s) const {
+  return (*this <= s) && (*this >= s);
+}
+
+template<class T, int D>
+bool bboxset<T,D>::operator!= (const bboxset<T,D>& s) const {
+  return ! (*this == s);
+}
+
+
+
 // Output
 template<class T,int D>
 void bboxset<T,D>::output (ostream& os) const {
diff --git a/Carpet/CarpetLib/src/bboxset.hh b/Carpet/CarpetLib/src/bboxset.hh
index 6a832ae82..500f5821e 100644
--- a/Carpet/CarpetLib/src/bboxset.hh
+++ b/Carpet/CarpetLib/src/bboxset.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.hh,v 1.5 2001/03/27 22:26:31 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.hh,v 1.6 2001/12/14 16:39:41 schnetter Exp $
 
  ***************************************************************************/
 
@@ -37,10 +37,15 @@ using namespace std;
 // Forward declaration
 template<class T, int D> class bboxset;
 
-template<class T,int D>
-bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2);
-template<class T,int D>
-bboxset<T,D> operator- (const bbox<T,D>& b, const bboxset<T,D>& s);
+// template<class T,int D>
+// bboxset<T,D> operator+ (const bbox<T,D>& b1, const bbox<T,D>& b2);
+// template<class T,int D>
+// bboxset<T,D> operator+ (const bbox<T,D>& b, const bboxset<T,D>& s);
+
+// template<class T,int D>
+// bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2);
+// template<class T,int D>
+// bboxset<T,D> operator- (const bbox<T,D>& b, const bboxset<T,D>& s);
 
 // Output
 template<class T,int D>
@@ -86,6 +91,8 @@ public:
   bboxset& operator+= (const bboxset& s);
   bboxset operator+ (const box& b) const;
   bboxset operator+ (const bboxset& s) const;
+  static bboxset plus (const box& b1, const box& b2);
+  static bboxset plus (const box& b, const bboxset& s);
   
   // Union
   bboxset& operator|= (const box& b);
@@ -109,6 +116,14 @@ public:
   // friend bboxset operator- <T,D>(const box& b, const bboxset& s);
   static bboxset minus (const box& b, const bboxset& s);
   
+  // Equality
+  bool operator== (const bboxset& s) const;
+  bool operator!= (const bboxset& s) const;
+  bool operator< (const bboxset& s) const;
+  bool operator<= (const bboxset& s) const;
+  bool operator> (const bboxset& s) const;
+  bool operator>= (const bboxset& s) const;
+  
   // Iterators
   typedef typename bset::const_iterator const_iterator;
   typedef typename bset::iterator       iterator;
@@ -123,6 +138,16 @@ public:
 
 
 template<class T,int D>
+inline bboxset<T,D> operator+ (const bbox<T,D>& b1, const bbox<T,D>& b2) {
+  return bboxset<T,D>::plus(b1,b2);
+}
+
+template<class T,int D>
+inline bboxset<T,D> operator+ (const bbox<T,D>& b, const bboxset<T,D>& s) {
+  return bboxset<T,D>::plus(b,s);
+}
+
+template<class T,int D>
 inline bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2) {
   return bboxset<T,D>::minus(b1,b2);
 }
@@ -132,6 +157,8 @@ inline bboxset<T,D> operator- (const bbox<T,D>& b, const bboxset<T,D>& s) {
   return bboxset<T,D>::minus(b,s);
 }
 
+
+
 // Output
 template<class T,int D>
 inline ostream& operator<< (ostream& os, const bboxset<T,D>& s) {
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index efca6c699..e2a7a8dae 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.14 2001/12/09 16:43:09 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.15 2001/12/14 16:39:41 schnetter Exp $
 
  ***************************************************************************/
 
@@ -20,7 +20,6 @@
 
 #include <assert.h>
 
-#include <fstream>
 #include <iostream>
 #include <string>
 
@@ -226,7 +225,7 @@ void data<T,D>
   
   assert (order_space <= 1);
   
-  assert (srcs.size() > order_time);
+  assert ((int)srcs.size() > order_time);
   vector<T> src_fac(order_time+1);
   for (int t=0; t<(int)src_fac.size(); ++t) {
     src_fac[t] = 1;
diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc
index 41e265238..1d8f3cfd7 100644
--- a/Carpet/CarpetLib/src/gdata.cc
+++ b/Carpet/CarpetLib/src/gdata.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.15 2001/12/09 16:43:10 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.16 2001/12/14 16:39:42 schnetter Exp $
 
  ***************************************************************************/
 
@@ -20,8 +20,7 @@
 
 #include <assert.h>
 
-#include <fstream>
-#include <iomanip>
+#include <iostream>
 
 #include "bbox.hh"
 #include "defs.hh"
@@ -132,7 +131,7 @@ void generic_data<D>
 // Output
 template<int D>
 template<int DD>
-void generic_data<D>::write_ascii (const string name, const int time,
+void generic_data<D>::write_ascii (ostream& os, const int time,
 				   const vect<int,D>& org,
 				   const vect<int,DD>& dirs,
 				   const int tl, const int rl,
@@ -151,19 +150,16 @@ void generic_data<D>::write_ascii (const string name, const int time,
     MPI_Comm_rank (dist::comm, &rank);
     if (rank == 0) {
       
-      ofstream file(name.c_str(), ios::app);
-      assert (file.good());
+      assert (os.good());
       
-      file << setprecision(15);
-      
-      file << "# iteration " << time << endl
-	   << "# time level " << tl << "   refinement level " << rl
-	   << "   component " << c << "   multigrid level " << ml << endl
-	   << "# column format: it tl rl c ml";
+      os << "# iteration " << time << endl
+	 << "# time level " << tl << "   refinement level " << rl
+	 << "   component " << c << "   multigrid level " << ml << endl
+	 << "# column format: it tl rl c ml";
       assert (D>=1 && D<=3);
       const char* const coords = "xyz";
-      for (int d=0; d<D; ++d) file << " " << coords[d];
-      file << " data" << endl;
+      for (int d=0; d<D; ++d) os << " " << coords[d];
+      os << " data" << endl;
       
       const vect<int,DD> lo = extent().lower()[dirs];
       const vect<int,DD> up = extent().upper()[dirs];
@@ -179,25 +175,24 @@ void generic_data<D>::write_ascii (const string name, const int time,
 	for (bbox<int,DD>::iterator it=ext.begin(); it!=ext.end(); ++it) {
 	  ivect index(org);
 	  for (int d=0; d<DD; ++d) index[dirs[d]] = (*it)[d];
-	  file << time << " " << tl << " " << rl << " " << c << " " << ml
-	       << " ";
-	  for (int d=0; d<D; ++d) file << index[d] << " ";
-	  write_ascii_output_element (file, index);
-	  file << endl;
+	  os << time << " " << tl << " " << rl << " " << c << " " << ml
+	     << " ";
+	  for (int d=0; d<D; ++d) os << index[d] << " ";
+	  write_ascii_output_element (os, index);
+	  os << endl;
 	  for (int d=0; d<DD; ++d) {
 	    if (index[dirs[d]]!=extent().upper()[dirs[d]]) break;
-	    file << endl;
+	    os << endl;
 	  }
 	}
 	
       }	else {
 	
-	file << "#" << endl;
+	os << "#" << endl;
 	
       }	// if ! ext contains org
       
-      file.close();
-      assert (file.good());
+      assert (os.good());
       
     }
     
@@ -207,7 +202,7 @@ void generic_data<D>::write_ascii (const string name, const int time,
     generic_data* const tmp = make_typed();
     tmp->allocate(extent(), 0);
     tmp->copy_from (this, extent());
-    tmp->write_ascii (name, time, org, dirs, tl, rl, c, ml);
+    tmp->write_ascii (os, time, org, dirs, tl, rl, c, ml);
     delete tmp;
     
   }
@@ -220,15 +215,15 @@ void generic_data<D>::write_ascii (const string name, const int time,
 template class generic_data<3>;
 
 template void generic_data<3>
-::write_ascii (const string name, const int time,
+::write_ascii (ostream& os, const int time,
 	       const vect<int,3>& org, const vect<int,1>& dirs,
 	       const int tl, const int rl, const int c, const int ml) const;
 template void generic_data<3>
-::write_ascii (const string name, const int time,
+::write_ascii (ostream& os, const int time,
 	       const vect<int,3>& org, const vect<int,2>& dirs,
 	       const int tl, const int rl, const int c, const int ml) const;
 template void generic_data<3>
-::write_ascii (const string name, const int time,
+::write_ascii (ostream& os, const int time,
 	       const vect<int,3>& org, const vect<int,3>& dirs,
 	       const int tl, const int rl, const int c, const int ml) const;
 
diff --git a/Carpet/CarpetLib/src/gdata.hh b/Carpet/CarpetLib/src/gdata.hh
index e827737f8..587613f4e 100644
--- a/Carpet/CarpetLib/src/gdata.hh
+++ b/Carpet/CarpetLib/src/gdata.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.10 2001/12/09 16:43:10 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.11 2001/12/14 16:39:42 schnetter Exp $
 
  ***************************************************************************/
 
@@ -24,7 +24,6 @@
 #include <assert.h>
 #include <stdlib.h>
 
-#include <fstream>
 #include <iostream>
 #include <string>
 
@@ -118,7 +117,7 @@ public:
   
   // Output
   template<int DD>
-  void write_ascii (const string name, const int time,
+  void write_ascii (ostream& os, const int time,
                     const vect<int,D>& org, const vect<int,DD>& dirs,
 		    const int tl, const int rl,
                     const int c, const int ml)
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 0408417d4..4cc2ccfae 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -6,7 +6,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.13 2001/12/09 16:43:10 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.14 2001/12/14 16:39:42 schnetter Exp $
 
  ***************************************************************************/
 
@@ -70,26 +70,9 @@ bool generic_gf<D>::operator== (const generic_gf<D>& f) const {
 // Modifiers
 template<int D>
 void generic_gf<D>::recompose () {
-  // Retain storage that might be needed
-  fdata oldstorage;
   
-  oldstorage.resize(tmax-tmin+1);
-  for (int tl=tmin; tl<=tmax; ++tl) {
-    oldstorage[tl-tmin].resize
-      (min(h.reflevels(), (int)storage[tl-tmin].size()));
-    for (int rl=0; rl<(int)storage[tl-tmin].size(); ++rl) {
-      oldstorage[tl-tmin][rl].resize
-        (min(h.components(rl), (int)storage[tl-tmin][rl].size()));
-      for (int c=0; c<(int)storage[tl-tmin][rl].size(); ++c) {
-        oldstorage[tl-tmin][rl][c].resize
-	  (min(h.mglevels(rl,c), (int)storage[tl-tmin][rl][c].size()));
-        for (int ml=0; ml<(int)storage[tl-tmin][rl][ml].size(); ++ml) {
-          oldstorage[tl-tmin][rl][c][ml]->transfer_from
-            (storage[tl-tmin][tl][c][ml]);
-        } // for ml
-      } // for c
-    } // for rl
-  } // for tl
+  // Retain storage that might be needed
+  fdata oldstorage = storage;
   
   // Resize structure
   storage.resize(tmax-tmin+1);
@@ -99,10 +82,14 @@ void generic_gf<D>::recompose () {
       storage[tl-tmin][rl].resize(h.components(rl));
       for (int c=0; c<h.components(rl); ++c) {
       	storage[tl-tmin][rl][c].resize(h.mglevels(rl,c));
+	for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
+	  storage[tl-tmin][rl][c][ml] = 0;
+	} // for ml
       } // for c
     } // for rl
   } // for tl
   
+  // Initialise the new storage
   for (int rl=0; rl<h.reflevels(); ++rl) {
     for (int tl=tmin; tl<=tmax; ++tl) {
       for (int c=0; c<h.components(rl); ++c) {
@@ -136,16 +123,22 @@ void generic_gf<D>::recompose () {
 	} // for ml
       } // for c
       
-      // Delete old storage
-      if (rl<(int)oldstorage[tl-tmin].size()) {
-      	oldstorage[tl-tmin][rl].clear();
-      }
-      
+    } // for tl
+  } // for rl
+  
+  // Delete old storage
+  for (int tl=tmin; tl<=tmax; ++tl) {
+    for (int rl=0; rl<(int)oldstorage[tl-tmin].size(); ++rl) {
+      for (int c=0; c<(int)oldstorage[tl-tmin][rl].size(); ++c) {
+        for (int ml=0; ml<(int)oldstorage[tl-tmin][rl][c].size(); ++ml) {
+	  delete oldstorage[tl-tmin][rl][c][ml];
+        } // for ml
+      } // for c
     } // for rl
   } // for tl
   
-  for (int tl=tmin; tl<=tmax; ++tl) {
-    for (int rl=0; rl<h.reflevels(); ++rl) {
+  for (int rl=0; rl<h.reflevels(); ++rl) {
+    for (int tl=tmin; tl<=tmax; ++tl) {
       
       // Set boundaries
       for (int c=0; c<h.components(rl); ++c) {
@@ -158,8 +151,8 @@ void generic_gf<D>::recompose () {
       	} // for ml
       } // for c
       
-    } // for rl
-  } // for tl
+    } // for tl
+  } // for rl
 }
 
 // Cycle the time levels by rotating the data sets
@@ -411,20 +404,11 @@ template<int D>
 void generic_gf<D>::ref_bnd_prolongate (int tl, int rl, int c, int ml) {
   // Interpolate
   assert (rl>=1);
-//   const int tmod
-//     = ((t.time(tl,rl,ml) - t.get_time(rl-1,ml)) % t.get_delta(rl-1,ml)
-//        + t.get_delta(rl-1,ml)) % t.get_delta(rl-1,ml);
   vector<int> tl2s;
-//   if (tmod == 0) {
-//     // No interpolation in time
-//     tl2s.resize(1);
-//     tl2s[0] = tl;
-//   } else {
-    // Interpolation in time
-    assert (tmax-tmin+1 >= prolongation_order_time+1);
-    tl2s.resize(prolongation_order_time+1);
-    for (int i=0; i<=prolongation_order_time; ++i) tl2s[i] = tmax - i;
-//   }
+  // Interpolation in time
+  assert (tmax-tmin+1 >= prolongation_order_time+1);
+  tl2s.resize(prolongation_order_time+1);
+  for (int i=0; i<=prolongation_order_time; ++i) tl2s[i] = tmax - i;
   intercat (tl  ,rl  ,c,ml, &dh<D>::dboxes::recv_ref_bnd_coarse,
 	    tl2s,rl-1,  ml, &dh<D>::dboxes::send_ref_bnd_fine);
 }
@@ -466,23 +450,12 @@ void generic_gf<D>::ref_restrict (int tl, int rl, int c, int ml)
 template<int D>
 void generic_gf<D>::ref_prolongate (int tl, int rl, int c, int ml)
 {
-  // Require same times
-  assert (t.get_time(rl,ml) == t.get_time(rl-1,ml));
   assert (rl>=1);
-//   const int tmod
-//     = ((t.time(tl,rl,ml) - t.get_time(rl-1,ml)) % t.get_delta(rl-1,ml)
-//        + t.get_delta(rl-1,ml)) % t.get_delta(rl-1,ml);
   vector<int> tl2s;
-//   if (tmod == 0) {
-//     // No interpolation in time
-//     tl2s.resize(1);
-//     tl2s[0] = tl;
-//   } else {
-    // Interpolation in time
-    assert (tmax-tmin+1 >= prolongation_order_time+1);
-    tl2s.resize(prolongation_order_time+1);
-    for (int i=0; i<=prolongation_order_time; ++i) tl2s[i] = tmax - i;
-//   }
+  // Interpolation in time
+  assert (tmax-tmin+1 >= prolongation_order_time+1);
+  tl2s.resize(prolongation_order_time+1);
+  for (int i=0; i<=prolongation_order_time; ++i) tl2s[i] = tmax - i;
   intercat (tl  ,rl  ,c,ml, &dh<D>::dboxes::recv_ref_coarse,
 	    tl2s,rl-1,  ml, &dh<D>::dboxes::send_ref_fine);
 }
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 6cb6fcdda..7c307860b 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -7,7 +7,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.9 2001/07/04 12:29:52 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.10 2001/12/14 16:39:42 schnetter Exp $
 
  ***************************************************************************/
 
@@ -151,50 +151,61 @@ void gh<D>::recompose (const rexts& exts, const rprocs& procs) {
 
 // Helpers
 template<int D>
+gh<D>::cexts gh<D>::make_reflevel_multigrid_boxes (const vector<ibbox>& exts,
+						   const int mglevels)
+  const
+{
+  assert (mglevels>0);
+  
+  cexts mexts (exts.size());
+  for (int c=0; c<(int)exts.size(); ++c) {
+    
+    mexts[c].resize(mglevels);
+    
+    ibbox ext = exts[c];
+    for (int ml=0; ml<mglevels; ++ml) {
+      
+      mexts[c][ml] = ext;
+      
+      if (ml == mglevels-1) break;
+      
+      // This level's characteristics
+      ivect str = ext.stride();
+      ivect lo = ext.lower();
+      ivect up = ext.upper();
+      
+      // Transform to next (coarser) level
+      switch (mgcent) {
+      case vertex_centered:
+	break;
+      case cell_centered:
+	for (int d=0; d<D; ++d) assert (str[d]%2 == 0);
+	lo += str/2;
+	break;
+      default:
+	abort();
+      }
+      str *= mgfact;
+      up = up - (up - lo) % str;
+      
+      ext = ibbox(lo,up,str);
+    }	// for ml
+  } // for c
+  
+  return mexts;
+}
+
+template<int D>
 gh<D>::rexts gh<D>::make_multigrid_boxes (const vector<vector<ibbox> >& exts,
 					  const int mglevels)
   const
 {
   assert (mglevels>0);
   
-  rexts mexts;
-  mexts.resize(exts.size());
+  rexts mexts (exts.size());
   for (int rl=0; rl<(int)exts.size(); ++rl) {
-    mexts[rl].resize(exts[rl].size());
-    for (int c=0; c<(int)exts[rl].size(); ++c) {
-      
-      mexts[rl][c].resize(mglevels);
-      
-      ibbox ext = exts[rl][c];
-      for (int ml=0; ml<mglevels; ++ml) {
-	
-	mexts[rl][c][ml] = ext;
-	
-	if (ml == mglevels-1) break;
-	
-	// This level's characteristics
-	ivect str = ext.stride();
-	ivect lo = ext.lower();
-	ivect up = ext.upper();
-	
-	// Transform to next (coarser) level
-	switch (mgcent) {
-	case vertex_centered:
-	  break;
-	case cell_centered:
-	  for (int d=0; d<D; ++d) assert (str[d]%2 == 0);
-	  lo += str/2;
-	  break;
-	default:
-	  abort();
-	}
-	str *= mgfact;
-	up = up - (up - lo) % str;
-	
-	ext = ibbox(lo,up,str);
-      }	// for ml
-    } // for c
-  } // for rl
+    mexts[rl] = make_reflevel_multigrid_boxes (exts[rl], mglevels);
+  }
   
   return mexts;
 }
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index 0b26b48e9..4748fb1de 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -7,7 +7,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.7 2001/09/06 11:18:28 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.8 2001/12/14 16:39:43 schnetter Exp $
 
  ***************************************************************************/
 
@@ -87,6 +87,9 @@ public:
   void recompose (const rexts& exts, const rprocs& procs);
   
   // Helpers
+  cexts make_reflevel_multigrid_boxes (const vector<ibbox>& exts,
+				       const int mglevels)
+    const;
   rexts make_multigrid_boxes (const vector<vector<ibbox> >& exts,
 			      const int mglevels)
     const;
diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh
index 5905d14b9..00a31da04 100644
--- a/Carpet/CarpetLib/src/vect.hh
+++ b/Carpet/CarpetLib/src/vect.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.5 2001/12/05 03:31:57 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.6 2001/12/14 16:39:43 schnetter Exp $
 
  ***************************************************************************/
 
@@ -505,7 +505,7 @@ template<class T,int D>
 inline T maxval (const vect<T,D>& a) {
   assert (D>0);
   T r(a[0]);
-  for (int d=0; d<D; ++d) r=max(r,a[d]);
+  for (int d=1; d<D; ++d) r=max(r,a[d]);
   return r;
 }
 
@@ -513,7 +513,23 @@ template<class T,int D>
 inline T minval (const vect<T,D>& a) {
   assert (D>0);
   T r(a[0]);
-  for (int d=0; d<D; ++d) r=min(r,a[d]);
+  for (int d=1; d<D; ++d) r=min(r,a[d]);
+  return r;
+}
+
+template<class T,int D>
+inline int maxloc (const vect<T,D>& a) {
+  assert (D>0);
+  int r(0);
+  for (int d=1; d<D; ++d) if (a[d]>a[r]) r=d;
+  return r;
+}
+
+template<class T,int D>
+inline int minloc (const vect<T,D>& a) {
+  assert (D>0);
+  int r(0);
+  for (int d=1; d<D; ++d) if (a[d]<a[r]) r=d;
   return r;
 }