Carpet: Initialise outer_boundaries during recursive load balancing

2 files changed, 49 insertions, 29 deletions

diff --git a/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.F90 b/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.F90
index bdd2f9601..441bf2d78 100644
--- a/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.F90
+++ b/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.F90
@@ -41,11 +41,12 @@ subroutine splitregions_recursively ( &
        CCTK_POINTER, intent(out) :: cxx_superreg
      end subroutine carpet_get_region
      
-     subroutine carpet_get_bbox (cxx_superreg, box)
+     subroutine carpet_get_bbox (cxx_superreg, box, obound)
        use carpet_boxtypes
        implicit none
-       CCTK_POINTER, intent(in)  :: cxx_superreg
-       type(bbox),   intent(out) :: box
+       CCTK_POINTER,   intent(in)  :: cxx_superreg
+       type(bbox),     intent(out) :: box
+       type(boundary), intent(out) :: obound
      end subroutine carpet_get_bbox
      
      subroutine carpet_insert_region (cxx_regs, reg)
@@ -87,11 +88,10 @@ subroutine splitregions_recursively ( &
   
   
   
-  outbound%obound(:,:) = 1
   allocate (sregions(nsuperregs))
   do i=1, nsuperregs
      call carpet_get_region (cxx_superregs, i-1, cxx_superreg)
-     call carpet_get_bbox (cxx_superreg, box)
+     call carpet_get_bbox (cxx_superreg, box, outbound)
      call create_sregion (box, outbound, i-1, sregions(i)%point)
   end do
   
@@ -119,6 +119,7 @@ contains
     
     integer                   :: nch, ich
     integer                   :: dir
+    integer                   :: mydir
     integer,      allocatable :: bounds(:)
     CCTK_POINTER, allocatable :: cxx_subtrees(:)
     type(superregion2slim)    :: sregslim
@@ -133,31 +134,38 @@ contains
        if (nch /= 2) then
           call CCTK_WARN (CCTK_WARN_ABORT, "number of children is not 2")
        end if
+       mydir = -1
        do dir=1, 3
           if (sreg%children(1)%point%extent%dim(dir)%upper + &
                sreg%children(2)%point%extent%dim(dir)%stride == &
                sreg%children(2)%point%extent%dim(dir)%lower) then
-             goto 100
-          end if
-          if (sreg%children(1)%point%extent%dim(dir)%lower /= &
-               sreg%children(2)%point%extent%dim(dir)%lower .or. &
-               sreg%children(1)%point%extent%dim(dir)%upper /= &
-               sreg%children(2)%point%extent%dim(dir)%upper) then
-             call CCTK_WARN (CCTK_WARN_ABORT, "children differ in unexpected ways")
+             ! Found direction
+             if (mydir > 0) then
+                call CCTK_WARN (CCTK_WARN_ABORT, "could not determine direction")
+             end if
+             mydir = dir
+          else
+             if (sreg%children(1)%point%extent%dim(dir)%lower /= &
+                  sreg%children(2)%point%extent%dim(dir)%lower .or. &
+                  sreg%children(1)%point%extent%dim(dir)%upper /= &
+                  sreg%children(2)%point%extent%dim(dir)%upper) then
+                call CCTK_WARN (CCTK_WARN_ABORT, "children differ in unexpected ways")
+             end if
           end if
        end do
-       call CCTK_WARN (CCTK_WARN_ABORT, "could not determine direction")
-100    continue
-       bounds(1) = sreg%children(1)%point%extent%dim(dir)%lower
-       bounds(2) = sreg%children(2)%point%extent%dim(dir)%lower
-       bounds(3) = sreg%children(2)%point%extent%dim(dir)%upper + &
-            sreg%children(2)%point%extent%dim(dir)%stride
+       if (mydir < 0) then
+          call CCTK_WARN (CCTK_WARN_ABORT, "could not determine direction")
+       end if
+       bounds(1) = sreg%children(1)%point%extent%dim(mydir)%lower
+       bounds(2) = sreg%children(2)%point%extent%dim(mydir)%lower
+       bounds(3) = sreg%children(2)%point%extent%dim(mydir)%upper + &
+            sreg%children(2)%point%extent%dim(mydir)%stride
        do ich=1, nch
           call insert_region &
                (sreg%children(ich)%point, cxx_subtrees(ich), cxx_regs)
        end do
        call carpet_create_tree_branch &
-            (nch, dir-1, bounds, cxx_subtrees, cxx_tree)
+            (nch, mydir-1, bounds, cxx_subtrees, cxx_tree)
     else
        ! The region is a leaf: insert it
        sregslim%extent           = sreg%extent
diff --git a/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.cc b/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.cc
index 9359e109d..893545f32 100644
--- a/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.cc
+++ b/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.cc
@@ -76,7 +76,7 @@ namespace Carpet {
   struct f_bbox {
     f_range dim[3];
     
-    f_bbox () { }
+    f_bbox () {}
     f_bbox (ibbox const& box)
     {
       assert (::dim == 3);
@@ -102,6 +102,16 @@ namespace Carpet {
   struct f_boundary {
     int obound[2][3];
     
+    f_boundary () {}
+    f_boundary (b2vect const& ob)
+    {
+      assert (::dim == 3);
+      for (int d=0; d<3; ++d) {
+        for (int f=0; f<2; ++f) {
+          obound[f][d] = ob[f][d];
+        }
+      }
+    }
     /*explicit*/ operator b2vect () const
     {
       b2vect ob;
@@ -433,7 +443,8 @@ namespace Carpet {
                                  int const& i,
                                  CCTK_POINTER& cxx_superreg)
   {
-    vector<region_t>& superregs = *(vector<region_t>*)cxx_superregs;
+    vector<region_t>& superregs =
+      *static_cast<vector<region_t>*>(cxx_superregs);
     region_t& superreg = superregs.AT(i);
     cxx_superreg = &superreg;
   }
@@ -441,10 +452,11 @@ namespace Carpet {
   extern "C"
   CCTK_FCALL void
   CCTK_FNAME(carpet_get_bbox) (CCTK_POINTER& cxx_superreg,
-                               f_bbox& box)
+                               f_bbox& box, f_boundary& obound)
   {
-    region_t& superreg = *(region_t*)cxx_superreg;
+    region_t& superreg = *static_cast<region_t*>(cxx_superreg);
     box = f_bbox(superreg.extent);
+    obound = f_boundary(superreg.outer_boundaries);
   }
   
   extern "C"
@@ -452,8 +464,8 @@ namespace Carpet {
   CCTK_FNAME(carpet_insert_region) (CCTK_POINTER& cxx_regs,
                                     f_superregion2slim const& reg)
   {
-    vector<region_t>& regs = *(vector<region_t>*)cxx_regs;
-    regs.push_back (region_t (reg));
+    vector<region_t>& regs = *static_cast<vector<region_t>*>(cxx_regs);
+    regs.push_back(region_t(reg));
   }
   
   extern "C"
@@ -470,7 +482,7 @@ namespace Carpet {
       bounds.AT(i) = fbounds[i];
     }
     for (int i=0; i<nch; ++i) {
-      ipfulltree* const tree = (ipfulltree*)cxx_subtrees[i];
+      ipfulltree* const tree = static_cast<ipfulltree*>(cxx_subtrees[i]);
       assert (tree->invariant());
       subtrees.AT(i) = tree;
     }
@@ -482,7 +494,7 @@ namespace Carpet {
   CCTK_FNAME(carpet_create_tree_leaf) (f_superregion2slim const& sreg,
                                        CCTK_POINTER& cxx_tree)
   {
-    cxx_tree = new ipfulltree (pseudoregion_t (sreg));
+    cxx_tree = new ipfulltree(pseudoregion_t(sreg));
   }
   
   extern "C"
@@ -490,8 +502,8 @@ namespace Carpet {
   CCTK_FNAME(carpet_set_tree) (CCTK_POINTER& cxx_superreg,
                                CCTK_POINTER& cxx_tree)
   {
-    region_t& superreg = *(region_t*)cxx_superreg;
-    ipfulltree* tree = (ipfulltree*)cxx_tree;
+    region_t& superreg = *static_cast<region_t*>(cxx_superreg);
+    ipfulltree* tree = static_cast<ipfulltree*>(cxx_tree);
     assert (not superreg.processors);
     superreg.processors = tree;
   }