Introduce a tree data structure to speed up domain decomposition

Introduce a tree data structure "fulltree", which decomposes a single,
rectangular region into a tree of non-overlapping, rectangular sub-regions.

Move the processor decomposition from the regridding thorns into Carpet.
Create such trees during processor decomposition.
Store these trees with the grid hierarchy.

11 files changed, 44 insertions, 93 deletions

diff --git a/Carpet/CarpetRegrid/interface.ccl b/Carpet/CarpetRegrid/interface.ccl
index 1cece9380..ac9d5b112 100644
--- a/Carpet/CarpetRegrid/interface.ccl
+++ b/Carpet/CarpetRegrid/interface.ccl
@@ -52,9 +52,11 @@ USES FUNCTION ConvertFromPhysicalBoundary
 
 # The true prototype of the routine below:
 #   int Carpet_Regrid (cGH const * cctkGH,
+#                      gh::rregs * superregss,
 #                      gh::mregs * regsss,
 #                      int         force);
 CCTK_INT FUNCTION Carpet_Regrid (CCTK_POINTER_TO_CONST IN cctkGH, \
+                                 CCTK_POINTER IN superregss,      \
                                  CCTK_POINTER IN regsss,          \
 				 CCTK_INT IN force)
 PROVIDES FUNCTION Carpet_Regrid WITH CarpetRegrid_Regrid LANGUAGE C
diff --git a/Carpet/CarpetRegrid/src/automatic.cc b/Carpet/CarpetRegrid/src/automatic.cc
index 16c12063c..442224f2c 100644
--- a/Carpet/CarpetRegrid/src/automatic.cc
+++ b/Carpet/CarpetRegrid/src/automatic.cc
@@ -26,15 +26,12 @@ namespace CarpetRegrid {
   
   int Automatic (cGH const * const cctkGH,
                  gh const & hh,
-                 gh::mregs & regsss)
+                 gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     
     assert (refinement_levels >= 1);
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     
     
     const int vi = CCTK_VarIndex (errorvar);
@@ -60,9 +57,6 @@ namespace CarpetRegrid {
       (cctkGH, hh, reflevel, min(reflevels+1, (int)refinement_levels),
        errorgf, regs);
     
-    // make multiprocessor aware
-    SplitRegions (cctkGH, regs);
-    
     if (regs.size() == 0) {
       // remove all finer levels
       regss.resize(reflevel+1);
@@ -77,9 +71,6 @@ namespace CarpetRegrid {
       }
     }
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/baselevel.cc b/Carpet/CarpetRegrid/src/baselevel.cc
index 105b9f72c..a7229fdad 100644
--- a/Carpet/CarpetRegrid/src/baselevel.cc
+++ b/Carpet/CarpetRegrid/src/baselevel.cc
@@ -19,14 +19,12 @@ namespace CarpetRegrid {
   
   int BaseLevel (cGH const * const cctkGH,
                  gh const & hh,
-                 gh::mregs & regsss)
+                 gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     
     assert (refinement_levels == 1);
     
-    assert (regsss.size() == 1);
-    
     return 0;
   }
   
diff --git a/Carpet/CarpetRegrid/src/centre.cc b/Carpet/CarpetRegrid/src/centre.cc
index 5bda85c85..078dd0181 100644
--- a/Carpet/CarpetRegrid/src/centre.cc
+++ b/Carpet/CarpetRegrid/src/centre.cc
@@ -19,7 +19,7 @@ namespace CarpetRegrid {
   
   int Centre (cGH const * const cctkGH,
               gh const & hh,
-              gh::mregs & regsss)
+              gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     
@@ -28,9 +28,6 @@ namespace CarpetRegrid {
     // do nothing if the levels already exist
     if (reflevel == refinement_levels) return 0;
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     regss.resize (refinement_levels);
     
     bvect const symmetric (symmetry_x, symmetry_y, symmetry_z);
@@ -40,7 +37,7 @@ namespace CarpetRegrid {
     ivect rlb  = hh.baseextents.at(0).at(0).lower();
     ivect rub  = hh.baseextents.at(0).at(0).upper();
     
-    assert (! smart_outer_boundaries);
+    assert (not smart_outer_boundaries);
     
     for (size_t rl=1; rl<regss.size(); ++rl) {
       
@@ -69,16 +66,10 @@ namespace CarpetRegrid {
       
       regs.at(0).map = Carpet::map;
       
-      // make multiprocessor aware
-      SplitRegions (cctkGH, regs);
-      
       regss.at(rl) = regs;
       
     } // for rl
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
index 7becee434..f32ec597c 100644
--- a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
+++ b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
@@ -23,7 +23,7 @@ namespace CarpetRegrid {
   
   int ManualCoordinateList (cGH const * const cctkGH,
                             gh const & hh,
-                            gh::mregs & regsss)
+                            gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     int ierr;
@@ -33,9 +33,6 @@ namespace CarpetRegrid {
     // do nothing if the levels already exist
     if (reflevel == refinement_levels && !tracking) return 0;
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     jjvect nboundaryzones, is_internal, is_staggered, shiftout;
     ierr = GetBoundarySpecification
       (2*dim, &nboundaryzones[0][0], &is_internal[0][0],
@@ -116,7 +113,7 @@ namespace CarpetRegrid {
         } // for c
       } // for rl
       
-    } else {                    // if ! smart_outer_boundaries
+    } else {                    // if not smart_outer_boundaries
       
       vector<vector<bbvect> > newobss1;
       if (strcmp(outerbounds, "") != 0) {
@@ -159,7 +156,7 @@ namespace CarpetRegrid {
         }
       }
 
-    } // if ! smart_outer_boundaries
+    } // if not smart_outer_boundaries
     
     for (int rl=1; rl<refinement_levels; ++rl) {
       
@@ -224,17 +221,11 @@ namespace CarpetRegrid {
         }
         
       } // if merge_overlapping_components
-
-      // make multiprocessor aware
-      SplitRegions (cctkGH, regs);
       
       regss.at(rl) = regs;
       
     } // for rl
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/manualcoordinates.cc b/Carpet/CarpetRegrid/src/manualcoordinates.cc
index e6d2d2ae0..1c4cea5b0 100644
--- a/Carpet/CarpetRegrid/src/manualcoordinates.cc
+++ b/Carpet/CarpetRegrid/src/manualcoordinates.cc
@@ -20,7 +20,7 @@ namespace CarpetRegrid {
   
   int ManualCoordinates (cGH const * const cctkGH,
                          gh const & hh,
-                         gh::mregs & regsss)
+                         gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     
@@ -32,9 +32,6 @@ namespace CarpetRegrid {
     // do nothing if the levels already exist
     if (reflevel == refinement_levels) return 0;
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     regss.resize (refinement_levels);
     
     vector<rvect> lower(3), upper(3);
@@ -59,16 +56,10 @@ namespace CarpetRegrid {
         (cctkGH, hh, rl, refinement_levels,
          lower.at(rl-1), upper.at(rl-1), reg, regs);
       
-      // make multiprocessor aware
-      SplitRegions (cctkGH, regs);
-      
       regss.at(rl) = regs;
       
     } // for rl
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/manualgridpointlist.cc b/Carpet/CarpetRegrid/src/manualgridpointlist.cc
index 60968e59d..fe5bdb4ab 100644
--- a/Carpet/CarpetRegrid/src/manualgridpointlist.cc
+++ b/Carpet/CarpetRegrid/src/manualgridpointlist.cc
@@ -23,7 +23,7 @@ namespace CarpetRegrid {
   
   int ManualGridpointList (cGH const * const cctkGH,
                            gh const & hh,
-                           gh::mregs & regsss)
+                           gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     
@@ -32,9 +32,6 @@ namespace CarpetRegrid {
     // do nothing if the levels already exist
     if (reflevel == refinement_levels) return 0;
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     regss.resize (refinement_levels);
     
     vector<vector<ibbox> > newbbss;
@@ -116,16 +113,10 @@ namespace CarpetRegrid {
            ext.lower(), ext.upper(), reg, regs);
       }
       
-      // make multiprocessor aware
-      SplitRegions (cctkGH, regs);
-      
       regss.at(rl) = regs;
       
     } // for rl
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/manualgridpoints.cc b/Carpet/CarpetRegrid/src/manualgridpoints.cc
index 6481247c9..3df1f2119 100644
--- a/Carpet/CarpetRegrid/src/manualgridpoints.cc
+++ b/Carpet/CarpetRegrid/src/manualgridpoints.cc
@@ -22,7 +22,7 @@ namespace CarpetRegrid {
   
   int ManualGridpoints (cGH const * const cctkGH,
                         gh const & hh,
-                        gh::mregs & regsss)
+                        gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     
@@ -34,9 +34,6 @@ namespace CarpetRegrid {
     // do nothing if the levels already exist
     if (reflevel == refinement_levels) return 0;
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     regss.resize (refinement_levels);
     
     vector<ivect> ilower(3), iupper(3);
@@ -47,7 +44,7 @@ namespace CarpetRegrid {
     ilower.at(2) = ivect (l3ixmin, l3iymin, l3izmin);
     iupper.at(2) = ivect (l3ixmax, l3iymax, l3izmax);
     
-    assert (! smart_outer_boundaries);
+    assert (not smart_outer_boundaries);
     
     for (size_t rl=1; rl<regss.size(); ++rl) {
       
@@ -61,16 +58,10 @@ namespace CarpetRegrid {
         (cctkGH, hh, rl,refinement_levels,
          ilower.at(rl-1), iupper.at(rl-1), reg, regs);
       
-      // make multiprocessor aware
-      SplitRegions (cctkGH, regs);
-      
       regss.at(rl) = regs;
       
     } // for rl
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/moving.cc b/Carpet/CarpetRegrid/src/moving.cc
index 41980ea45..075fdbb93 100644
--- a/Carpet/CarpetRegrid/src/moving.cc
+++ b/Carpet/CarpetRegrid/src/moving.cc
@@ -20,16 +20,13 @@ namespace CarpetRegrid {
   
   int Moving (cGH const * const cctkGH,
               gh const & hh,
-              gh::mregs & regsss)
+              gh::rregs & regss)
   {
     DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
     
     assert (refinement_levels >= 1);
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     regss.resize (refinement_levels);
     
     bvect const symmetric (symmetry_x, symmetry_y, symmetry_z);
@@ -60,16 +57,10 @@ namespace CarpetRegrid {
       ManualCoordinates_OneLevel
         (cctkGH, hh, rl, refinement_levels, rlb, rub, reg, regs);
       
-      // make multiprocessor aware
-      SplitRegions (cctkGH, regs);
-      
       regss.at(rl) = regs;
       
     } // for rl
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/regrid.cc b/Carpet/CarpetRegrid/src/regrid.cc
index 2dd2a33a9..0a5567e2b 100644
--- a/Carpet/CarpetRegrid/src/regrid.cc
+++ b/Carpet/CarpetRegrid/src/regrid.cc
@@ -19,6 +19,7 @@ namespace CarpetRegrid {
   using namespace Carpet;
   
   CCTK_INT CarpetRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
+                                CCTK_POINTER const superregss_,
                                 CCTK_POINTER const regsss_,
 				CCTK_INT force)
   {
@@ -26,6 +27,7 @@ namespace CarpetRegrid {
     
     const cGH * const cctkGH = (const cGH *) cctkGH_;
     
+    gh::rregs & superregss = * (gh::rregs *) superregss_;
     gh::mregs & regsss = * (gh::mregs *) regsss_;
     
     gh const & hh = *vhh.at(Carpet::map);
@@ -151,44 +153,54 @@ namespace CarpetRegrid {
     
     if (CCTK_EQUALS(refined_regions, "none")) {
       
-      do_recompose = BaseLevel (cctkGH, hh, regsss);
+      do_recompose = BaseLevel (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "centre")) {
       
-      do_recompose = Centre (cctkGH, hh, regsss);
+      do_recompose = Centre (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "manual-gridpoints")) {
       
       do_recompose
-        = ManualGridpoints (cctkGH, hh, regsss);
+        = ManualGridpoints (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "manual-coordinates")) {
       
       do_recompose
-        = ManualCoordinates (cctkGH, hh, regsss);
+        = ManualCoordinates (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "manual-gridpoint-list")) {
       
       do_recompose
-        = ManualGridpointList (cctkGH, hh, regsss);
+        = ManualGridpointList (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "manual-coordinate-list")) {
       
       do_recompose
-        = ManualCoordinateList (cctkGH, hh, regsss);
+        = ManualCoordinateList (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "moving")) {
       
-      do_recompose = Moving (cctkGH, hh, regsss);
+      do_recompose = Moving (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "automatic")) {
       
-      do_recompose = Automatic (cctkGH, hh, regsss);
+      do_recompose = Automatic (cctkGH, hh, superregss);
                  
     } else {
       assert (0);
     }
     
+    // make multiprocessor aware
+    vector<vector<region_t> > regss(superregss.size());
+    for (size_t rl=0; rl<superregss.size(); ++rl) {
+#warning "TODO: delete .processors"
+      SplitRegions (cctkGH, superregss.at(rl), regss.at(rl));
+    }
+    
+    // make multigrid aware
+    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
+    
     assert (regsss.size() > 0);
     
     return do_recompose;
diff --git a/Carpet/CarpetRegrid/src/regrid.hh b/Carpet/CarpetRegrid/src/regrid.hh
index b7d727317..f37c9cb1d 100644
--- a/Carpet/CarpetRegrid/src/regrid.hh
+++ b/Carpet/CarpetRegrid/src/regrid.hh
@@ -31,8 +31,10 @@ namespace CarpetRegrid {
     
     // Aliased functions
 //     CCTK_INT CarpetRegrid_Regrid (const cGH * const cctkGH,
+//                                   gh<dim>::rregs * superregss,
 //                                   gh<dim>::mregs * regsss);
     CCTK_INT CarpetRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
+                                  CCTK_POINTER const superregss_,
                                   CCTK_POINTER const regsss_,
 				  CCTK_INT force);
   }
@@ -41,15 +43,15 @@ namespace CarpetRegrid {
   
   int BaseLevel (cGH const * const cctkGH,
                  gh const & hh,
-                 gh::mregs & regsss);
+                 gh::rregs & regss);
 
   int Centre (cGH const * const cctkGH,
               gh const & hh,
-              gh::mregs & regsss);
+              gh::rregs & regss);
 
   int ManualGridpoints (cGH const * const cctkGH,
                         gh const & hh,
-                        gh::mregs & regsss);
+                        gh::rregs & regss);
   
   void ManualGridpoints_OneLevel (const cGH * const cctkGH,
                                   const gh & hh,
@@ -62,7 +64,7 @@ namespace CarpetRegrid {
 
   int ManualCoordinates (cGH const * const cctkGH,
                          gh const & hh,
-                         gh::mregs & regsss);
+                         gh::rregs & regss);
   
   void ManualCoordinates_OneLevel (const cGH * const cctkGH,
                                    const gh & hh,
@@ -84,19 +86,19 @@ namespace CarpetRegrid {
 
   int ManualGridpointList (cGH const * const cctkGH,
                            gh const & hh,
-                           gh::mregs & regsss);
+                           gh::rregs & regss);
   
   int ManualCoordinateList (cGH const * const cctkGH,
                             gh const & hh,
-                            gh::mregs & regsss);
+                            gh::rregs & regss);
 
   int Moving (cGH const * const cctkGH,
               gh const & hh,
-              gh::mregs & regsss);
+              gh::rregs & regss);
   
   int Automatic (cGH const * const cctkGH,
                  gh const & hh,
-                 gh::mregs & regsss);
+                 gh::rregs & regss);
   
   void Automatic_OneLevel (const cGH * const cctkGH,
                            const gh & hh,