Cleaned up the code to make it work with multiple multigrid levels

Cleaned up the code to make it work with multiple multigrid levels
(aka shadow hierarchy).  The shadow logic is not yet in place.

Added simple recombining to the clusterer.  This should lead to fewer
grid components.  Not very tested.

darcs-hash:20020109125624-07bb3-f2d22fa4583bf562101ab521606e6142585622a7.gz

1 files changed, 137 insertions, 13 deletions

diff --git a/Carpet/CarpetRegrid/src/regrid.cc b/Carpet/CarpetRegrid/src/regrid.cc
index 5e7e28588..1f3b96b16 100644
--- a/Carpet/CarpetRegrid/src/regrid.cc
+++ b/Carpet/CarpetRegrid/src/regrid.cc
@@ -1,5 +1,4 @@
 #include <assert.h>
-#include <math.h>
 #include <stdarg.h>
 #include <stdlib.h>
 
@@ -12,6 +11,7 @@
 
 #include "Carpet/CarpetLib/src/bbox.hh"
 #include "Carpet/CarpetLib/src/bboxset.hh"
+#include "Carpet/CarpetLib/src/defs.hh"
 #include "Carpet/CarpetLib/src/gh.hh"
 #include "Carpet/CarpetLib/src/gf.hh"
 #include "Carpet/CarpetLib/src/vect.hh"
@@ -19,7 +19,7 @@
 #include "carpet.hh"
 #include "regrid.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.2 2001/12/17 13:32:43 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.3 2002/01/09 13:56:28 schnetter Exp $";
 
 
 
@@ -112,7 +112,9 @@ namespace CarpetRegrid {
     
     vector<bbox<int,dim> > bbs;
     bbs.reserve (bbl.size());
-    for (list<bbox<int,dim> >::const_iterator it = bbl.begin(); it != bbl.end(); ++it) {
+    for (list<bbox<int,dim> >::const_iterator it = bbl.begin();
+	 it != bbl.end();
+	 ++it) {
       bbs.push_back (*it);
     }
     
@@ -121,7 +123,6 @@ namespace CarpetRegrid {
     SplitRegions (cctkGH, bbs);
     
     gh<dim>::cexts bbss;
-    const int mglevels = 1;	// arbitrary value
     bbss = hh->make_reflevel_multigrid_boxes(bbs, mglevels);
     
     gh<dim>::rexts bbsss = hh->extents;
@@ -184,8 +185,13 @@ namespace CarpetRegrid {
       rstr /= hh->reffact;
       // refine (arbitrarily) around the center only
       rlb = rcentre - (newrextent/2 / rstr) * rstr;
-      // // refine (arbitrarily) around the lower boundary only
-      // rlb = rlb;
+#if 0
+      // refine (arbitrarily) around the lower boundary only
+      rlb = rlb;
+#endif
+      // honour multigrid factors
+      const int mgstr = ipow(hh->mgfact, mglevels);
+      rlb = (rlb / mgstr) * mgstr;
       rub = rlb + newrextent;
       // require rub<oldrub because we really want rub-rstr<=oldrub-oldstr
       assert (all(rlb >= oldrlb && rub < oldrub));
@@ -211,7 +217,7 @@ namespace CarpetRegrid {
     
     const int rl = reflevel+1;
     
-    const int levfac = floor(pow((double)hh->reffact, rl) + 0.5);
+    const int levfac = ipow(hh->reffact, rl);
     assert (all (rstr % levfac == 0));
     const vect<int,dim> str (rstr / levfac);
     const vect<int,dim> lb  (lower[rl-1]);
@@ -238,22 +244,135 @@ namespace CarpetRegrid {
   
   
   static void
+  MakeRegions_Adaptively_Recombine (list<bbox<int,dim> >& bbl1,
+				    list<bbox<int,dim> >& bbl2,
+				    list<bbox<int,dim> >& bbl,
+				    const bbox<int,dim>& iface)
+  {
+#if 0
+    // remember old bounding boxes
+    bboxset<int,dim> oldboxes;
+    for (list<bbox<int,dim> >::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) {
+      oldboxes += *ibb1;
+    }
+    for (list<bbox<int,dim> >::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) {
+      oldboxes += *ibb2;
+    }
+#endif
+#if 0
+    cout << endl;
+    cout << "MakeRegions_Adaptively_Recombine: initial list:" << endl;
+    cout << bbl << endl;
+    cout << "MakeRegions_Adaptively_Recombine: initial list 1:" << endl;
+    cout << bbl1 << endl;
+    cout << "MakeRegions_Adaptively_Recombine: initial list 2:" << endl;
+    cout << bbl2 << endl;
+#endif
+    
+    {
+      // prune boxes on the left
+      list<bbox<int,dim> >::iterator ibb1 = bbl1.begin();
+      while (ibb1 != bbl1.end()) {
+	// is this bbox just to the left of the interface?
+	const bbox<int,dim> iface1 (ibb1->upper()+ibb1->stride(), ibb1->upper()+ibb1->stride(), ibb1->stride());
+	assert (!iface1.empty());
+	if (!iface1.is_contained_in(iface)) {
+	  // no: forget it
+	  bbl.push_back (*ibb1);
+	  ibb1 = bbl1.erase(ibb1);
+	  continue;
+	}
+	++ibb1;
+      }	// while
+    }
+    
+    {
+      // prune boxes on the right
+      list<bbox<int,dim> >::iterator ibb2 = bbl2.begin();
+      while (ibb2 != bbl2.end()) {
+	// is this bbox just to the right of the interface?
+	const bbox<int,dim> iface2 (ibb2->lower(), ibb2->lower(), ibb2->stride());
+	assert (!iface2.empty());
+	if (! iface2.is_contained_in(iface)) {
+	  // no: forget it
+	  bbl.push_back (*ibb2);
+	  ibb2 = bbl2.erase(ibb2);
+	  continue;
+	}
+	++ibb2;
+      }	// while
+    }
+    
+    {
+      // walk all boxes on the left
+      list<bbox<int,dim> >::iterator ibb1 = bbl1.begin();
+      while (ibb1 != bbl1.end()) {
+	const bbox<int,dim> iface1 (ibb1->upper()+ibb1->stride(), ibb1->upper()+ibb1->stride(), ibb1->stride());
+	assert (iface1.is_contained_in(iface));
+	
+	{
+	  // walk all boxes on the right
+	  list<bbox<int,dim> >::iterator ibb2 = bbl2.begin();
+	  while (ibb2 != bbl2.end()) {
+	    const bbox<int,dim> iface2 (ibb2->lower(), ibb2->lower(), ibb2->stride());
+	    assert (iface2.is_contained_in(iface));
+	    
+	    // check for a match
+	    if (iface1 == iface2) {
+	      assert (all(ibb1->stride() == ibb2->stride()));
+	      const bbox<int,dim> combined (ibb1->lower(), ibb2->upper(), ibb1->stride());
+	      bbl.push_back (combined);
+	      ibb1 = bbl1.erase(ibb1);
+	      ibb2 = bbl2.erase(ibb2);
+	      goto continue_search;
+	    }
+	    
+	    ++ibb2;
+	  } // while
+	}
+	
+	++ibb1;
+      continue_search:;
+      }	// while
+    }
+    
+    bbl.splice (bbl.end(), bbl1);
+    bbl.splice (bbl.end(), bbl2);
+    
+#if 0
+    // find new bounding boxes
+    bboxset<int,dim> newboxes;
+    for (list<bbox<int,dim> >::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
+      newboxes += *ibb;
+    }
+    
+    // Check that they are equal
+    assert (newboxes.size() <= oldboxes.size());
+    assert ((newboxes.size()==0) == (oldboxes.size()==0));
+    assert (oldboxes == newboxes);
+#endif
+#if 0
+    cout << "MakeRegions_Adaptively_Recombine: final list:" << endl;
+    cout << bbl << endl;
+    cout << endl;
+#endif
+  }
+  
+  static void
   MakeRegions_Adaptively (const cGH* const cctkGH,
 			  const int minwidth, const double minfraction,
 			  const CCTK_REAL maxerror,
 			  const data<CCTK_REAL,dim>& error,
 			  list<bbox<int,dim> >& bbl,
-			  const bbox<int,dim> region)
+			  const bbox<int,dim>& region)
   {
     // Just to be sure
     assert (! region.empty());
     
     // Count grid points that need to be refined
     int cnt = 0;
-    CCTK_REAL me=-999;
     for (bbox<int,dim>::iterator it=region.begin(); it!=region.end(); ++it) {
       if (error[*it] > maxerror) ++cnt;
-      me=max(me,error[*it]);
     }
     const CCTK_REAL fraction = (CCTK_REAL)cnt / region.num_points();
     const int width = minval(region.shape() / region.stride());
@@ -274,7 +393,8 @@ namespace CarpetRegrid {
       const vect<int,dim> str(region.stride());
       vect<int,dim> lo1(lo), lo2(lo);
       vect<int,dim> up1(up), up2(up);
-      up1[d] = (up[d]+lo[d])/2/str[d]*str[d];
+      const int mgstr = ipow(hh->mgfact, mglevels); // honour multigrid factors
+      up1[d] = ((up[d]+lo[d])/2/str[d]/mgstr)*str[d]*mgstr;
       lo2[d] = up1[d]+str[d];
       const bbox<int,dim> region1(lo1,up1,str);
       const bbox<int,dim> region2(lo2,up2,str);
@@ -282,10 +402,14 @@ namespace CarpetRegrid {
       assert (region2.is_contained_in(region));
       assert ((region1 & region2).empty());
       assert (region1 + region2 == region);
+      list<bbox<int,dim> > bbl1, bbl2;
       MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror,
-			      error, bbl, region1);
+			      error, bbl1, region1);
       MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror,
-			      error, bbl, region2);
+			      error, bbl2, region2);
+      // Combine regions if possible
+      const bbox<int,dim> iface(lo2,lo2,str);
+      MakeRegions_Adaptively_Recombine (bbl1, bbl2, bbl, iface);
     }
     
   }