dh_decompose

  Broke up dh<D>::recompose( bool ) into bite-sized chunks.

  Made use of member function pointers to cut down as many as three
  levels of blocks.  See 'foreach_reflevel_component_mglevel'.

  Advantages:

  1) the code of the main function 'recompose' now presents an outline of
     its functionality, rather than a post-ice-age deluge.

  2) component functions are named according to some comment found in the
     code.  For those who want to minimise comments.

  3) most of the functions conform to my "rule of thumb" (comprehensible
     chunks of code are no longer than my thumb)

  4) nesting substantially reduced (The functions that still have more
     than 5 levels of nested blocks are mostly those that are also
     unthumbworthy.  More room on the screen, but still room for
     improvement.)

darcs-hash:20041122173415-32473-616506e73f5c93be610bc773adaa7867663cb726.gz

1 files changed, 543 insertions, 542 deletions

diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index d946abdc1..23e1f2923 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -1,4 +1,4 @@
-// $Header:$
+// $Header: /home/cvs/carpet/Carpet/CarpetLib/src/dh.cc,v 1.57 2004/09/17 16:37:26 schnetter Exp $
 
 #include <assert.h>
 
@@ -20,7 +20,7 @@ using namespace std;
 template<int D>
 dh<D>::dh (gh<D>& h,
            const ivect& lghosts, const ivect& ughosts,
-	   const int prolongation_order_space, const int buffer_width)
+           const int prolongation_order_space, const int buffer_width)
   : h(h),
     lghosts(lghosts), ughosts(ughosts),
     prolongation_order_space(prolongation_order_space),
@@ -57,603 +57,563 @@ void dh<D>::recompose (const bool do_prolongate) {
   CHECKPOINT;
   
   boxes.clear();
+
+  allocate_bboxes();
   
+  foreach_reflevel_component_mglevel (&dh<D>::setup_sync_and_refine_boxes );
+  foreach_reflevel_component_mglevel (&dh<D>::intersect_sync_with_interior );
+  foreach_reflevel_component_mglevel (&dh<D>::setup_multigrid_boxes );
+  foreach_reflevel_component_mglevel (&dh<D>::setup_refinement_interior_boxes );
+  foreach_reflevel_component_mglevel (&dh<D>::setup_refinement_exterior_boxes );
+  foreach_reflevel_component_mglevel (&dh<D>::setup_restrict_interior_boxes );
+  foreach_reflevel_component_mglevel (&dh<D>::trim_unsynced_boundaries );
+
+  calculate_bases();
+
+  if (output_bboxes) {
+    cout << endl << h << endl;
+    foreach_reflevel_component_mglevel (&dh<D>::output_bboxes );
+    output_bases();
+  }
+  
+  foreach_reflevel_component_mglevel (&dh<D>::assert_assert_assert );
+
+  if (! save_memory_during_regridding) {
+    save_time( do_prolongate);
+  } else {
+    save_memory( do_prolongate);
+  }
+}
+
+template<int D>
+void dh<D>::allocate_bboxes() {
   boxes.resize(h.reflevels());
   for (int rl=0; rl<h.reflevels(); ++rl) {
     boxes.at(rl).resize(h.components(rl));
     for (int c=0; c<h.components(rl); ++c) {
       boxes.at(rl).at(c).resize(h.mglevels(rl,c));
       for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-       	const ibbox intr = h.extents.at(rl).at(c).at(ml);
-	
-       	// Interior
-       	// (the interior of the grid has the extent as specified by
-       	// the user)
-       	boxes.at(rl).at(c).at(ml).interior = intr;
-	
-       	// Exterior (add ghost zones)
-       	// (the content of the exterior is completely determined by
-       	// the interior of this or other components; the content of
-       	// the exterior is redundant)
-	ivect ldist(lghosts), udist(ughosts);
-	for (int d=0; d<D; ++d) {
-	  if (h.outer_boundaries.at(rl).at(c)[d][0]) ldist[d] = 0;
-	  if (h.outer_boundaries.at(rl).at(c)[d][1]) udist[d] = 0;
-	}
-        boxes.at(rl).at(c).at(ml).exterior = intr.expand(ldist, udist);
-	
-       	// Boundaries (ghost zones only)
-       	// (interior + boundaries = exterior)
-       	boxes.at(rl).at(c).at(ml).boundaries
-          = boxes.at(rl).at(c).at(ml).exterior - intr;
-	
-      } // for ml
-    } // for c
-  } // for rl
-  
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    for (int c=0; c<h.components(rl); ++c) {
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-	
-      	// Sync boxes
-      	const int cs = h.components(rl);
-      	boxes.at(rl).at(c).at(ml).send_sync.resize(cs);
-      	boxes.at(rl).at(c).at(ml).recv_sync.resize(cs);
-      	
-      	// Refinement boxes
-      	if (rl>0) {
-      	  const int csm1 = h.components(rl-1);
-      	  boxes.at(rl).at(c).at(ml).send_ref_coarse.resize(csm1);
-      	  boxes.at(rl).at(c).at(ml).recv_ref_coarse.resize(csm1);
-      	  boxes.at(rl).at(c).at(ml).recv_ref_bnd_coarse.resize(csm1);
-      	}
-      	if (rl<h.reflevels()-1) {
-      	  const int csp1 = h.components(rl+1);
-      	  boxes.at(rl).at(c).at(ml).recv_ref_fine.resize(csp1);
-      	  boxes.at(rl).at(c).at(ml).send_ref_fine.resize(csp1);
-      	  boxes.at(rl).at(c).at(ml).send_ref_bnd_fine.resize(csp1);
-      	}
-	
-      } // for ml
-    } // for c
-  } // for rl
-  
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    for (int c=0; c<h.components(rl); ++c) {
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-      	const ibset& bnds = boxes.at(rl).at(c).at(ml).boundaries;
-	
-      	// Sync boxes
-      	for (int cc=0; cc<h.components(rl); ++cc) {
-      	  assert (ml<h.mglevels(rl,cc));
-      	  // intersect boundaries with interior of that component
-          ibset ovlp = bnds & boxes.at(rl).at(cc).at(ml).interior;
-          ovlp.normalize();
-      	  for (typename ibset::const_iterator b=ovlp.begin();
-	       b!=ovlp.end(); ++b) {
-	    boxes.at(rl).at(c ).at(ml).recv_sync.at(cc).push_back(*b);
-	    boxes.at(rl).at(cc).at(ml).send_sync.at(c ).push_back(*b);
-	  }
-      	} // for cc
-      	
+        const ibbox intr = h.extents.at(rl).at(c).at(ml);
+        dboxes & b = boxes.at(rl).at(c).at(ml);
+
+        // Interior
+        // (the interior of the grid has the extent as specified by
+        // the user)
+        b.interior = intr;
+
+        // Exterior (add ghost zones)
+        // (the content of the exterior is completely determined by
+        // the interior of this or other components; the content of
+        // the exterior is redundant)
+        ivect ldist(lghosts), udist(ughosts);
+        for (int d=0; d<D; ++d) {
+          if (h.outer_boundaries.at(rl).at(c)[d][0]) ldist[d] = 0;
+          if (h.outer_boundaries.at(rl).at(c)[d][1]) udist[d] = 0;
+        }
+        b.exterior = intr.expand(ldist, udist);
+
+        // Boundaries (ghost zones only)
+        // (interior + boundaries = exterior)
+        b.boundaries = b.exterior - intr;
+
       } // for ml
     } // for c
   } // for rl
-  
+}
+
+// Loops over each refinement level, each component, and each multigrid
+// level, executing the "boxesop" member function argument on the corresponding
+// element of the "boxes" member 
+template<int D>
+void dh<D>::foreach_reflevel_component_mglevel (dh<D>::boxesop op) {
+
   for (int rl=0; rl<h.reflevels(); ++rl) {
     for (int c=0; c<h.components(rl); ++c) {
       for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-      	const ibbox& intr = boxes.at(rl).at(c).at(ml).interior;
-      	const ibbox& extr = boxes.at(rl).at(c).at(ml).exterior;
-	
-      	// Multigrid boxes
-      	if (ml>0) {
-      	  const ibbox intrf = boxes.at(rl).at(c).at(ml-1).interior;
-      	  const ibbox extrf = boxes.at(rl).at(c).at(ml-1).exterior;
-      	  // Restriction (interior)
-      	  {
-      	    // (the restriction must fill all of the interior of the
-      	    // coarse grid, and may use the exterior of the fine grid)
-      	    const ibbox recv = intr;
+        dboxes & box = boxes.at(rl).at(c).at(ml);
+        (this->*op)( box, rl, c, ml ); // evaluate member function
+      }
+    }
+  }
+}
+
+template<int D>
+void dh<D>::setup_multigrid_boxes( dh<D>::dboxes & box, int rl, int c, int ml )
+{
+  const ibbox& intr = box.interior;
+  const ibbox& extr = box.exterior;
+
+  // Multigrid boxes
+  if (ml>0) {
+    dboxes & bbox = boxes.at(rl).at(c).at(ml-1);
+    const ibbox intrf = bbox.interior;
+    const ibbox extrf = bbox.exterior;
+    // Restriction (interior)
+    {
+      // (the restriction must fill all of the interior of the
+      // coarse grid, and may use the exterior of the fine grid)
+      const ibbox recv = intr;
             assert (intr.empty() || ! recv.empty());
-      	    const ibbox send = recv.expanded_for(extrf);
+      const ibbox send = recv.expanded_for(extrf);
             assert (intr.empty() || ! send.empty());
-            // TODO: put the check back in, taking outer boundaries
-            // into account
+      // TODO: put the check back in, taking outer boundaries
+      // into account
 #if 0
-      	    assert (send.is_contained_in(extrf));
+      assert (send.is_contained_in(extrf));
 #endif
-      	    boxes.at(rl).at(c).at(ml-1).send_mg_coarse.push_back(send);
-      	    boxes.at(rl).at(c).at(ml  ).recv_mg_fine  .push_back(recv);
-      	  }
-      	  // Prolongation (interior)
-      	  {
-      	    // (the prolongation may use the exterior of the coarse
-      	    // grid, and may fill only the interior of the fine grid,
-      	    // and the bbox must be as large as possible)
-      	    const ibbox recv = extr.contracted_for(intrf) & intrf;
+      bbox.send_mg_coarse.push_back(send);
+      box .recv_mg_fine  .push_back(recv);
+    }
+    // Prolongation (interior)
+    {
+      // (the prolongation may use the exterior of the coarse
+      // grid, and may fill only the interior of the fine grid,
+      // and the bbox must be as large as possible)
+      const ibbox recv = extr.contracted_for(intrf) & intrf;
             assert (intr.empty() || ! recv.empty());
-      	    const ibbox send = recv.expanded_for(extr);
+      const ibbox send = recv.expanded_for(extr);
             assert (intr.empty() || ! send.empty());
-      	    boxes.at(rl).at(c).at(ml-1).recv_mg_coarse.push_back(recv);
-      	    boxes.at(rl).at(c).at(ml  ).send_mg_fine  .push_back(send);
-      	  }
-      	} // if not finest multigrid level
-        
-      } // for ml
-    } // for c
-  } // for rl
-  
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    for (int c=0; c<h.components(rl); ++c) {
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-      	const ibbox& intr = boxes.at(rl).at(c).at(ml).interior;
-      	const ibbox& extr = boxes.at(rl).at(c).at(ml).exterior;
-	
-      	// Refinement boxes
-      	if (rl<h.reflevels()-1) {
-      	  for (int cc=0; cc<h.components(rl+1); ++cc) {
-      	    const ibbox intrf = boxes.at(rl+1).at(cc).at(ml).interior;
-      	    // Prolongation (interior)
-            // TODO: prefer boxes from the same processor
-      	    {
-      	      // (the prolongation may use the exterior of the coarse
-      	      // grid, and must fill all of the interior of the fine
-      	      // grid)
-              const int pss = prolongation_stencil_size();
-              ibset recvs
-                = extr.expand(-pss,-pss).contracted_for(intrf) & intrf;
-              const iblistvect& rrc
-                = boxes.at(rl+1).at(cc).at(ml).recv_ref_coarse;
-              for (typename iblistvect::const_iterator lvi=rrc.begin();
-                   lvi!=rrc.end(); ++lvi) {
-                for (typename iblist::const_iterator li=lvi->begin();
-                     li!=lvi->end(); ++li) {
-                  recvs -= *li;
-                }
-              }
-              recvs.normalize();
-              assert (recvs.setsize() <= 1);
-              if (recvs.setsize() == 1) {
-                const ibbox recv = *recvs.begin();
-                const ibbox send = recv.expanded_for(extr);
-                assert (! send.empty());
-                assert (send.is_contained_in(extr));
-                boxes.at(rl+1).at(cc).at(ml).recv_ref_coarse.at(c )
-                  .push_back(recv);
-                boxes.at(rl  ).at(c ).at(ml).send_ref_fine  .at(cc)
-                  .push_back(send);
-              }
-      	    }
-            
-      	  } // for cc
-      	} // if not finest refinement level
-	
-      } // for ml
-    } // for c
-  } // for rl
-  
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    for (int c=0; c<h.components(rl); ++c) {
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-      	const ibbox& intr = boxes.at(rl).at(c).at(ml).interior;
-      	const ibbox& extr = boxes.at(rl).at(c).at(ml).exterior;
-	
-      	// Refinement boxes
-      	if (rl<h.reflevels()-1) {
-      	  for (int cc=0; cc<h.components(rl+1); ++cc) {
-      	    const ibbox intrf = boxes.at(rl+1).at(cc).at(ml).interior;
-            const ibbox& extrf = boxes.at(rl+1).at(cc).at(ml).exterior;
-            const ibset& bndsf = boxes.at(rl+1).at(cc).at(ml).boundaries;
-      	    // Prolongation (boundaries)
-            // TODO: prefer boxes from the same processor
-      	    {
-              // (the prolongation may use the exterior of the coarse
-              // grid, and must fill all of the boundary of the fine
-              // grid)
-              const int pss = prolongation_stencil_size();
-              // Prolongation boundaries
-              ibset pbndsf = bndsf;
-              {
-                // Do not count what is synced
-                const iblistvect& rs
-                  = boxes.at(rl+1).at(cc).at(ml).recv_sync;
-                for (typename iblistvect::const_iterator lvi=rs.begin();
-                     lvi!=rs.end(); ++lvi) {
-                  for (typename iblist::const_iterator li=lvi->begin();
-                       li!=lvi->end(); ++li) {
-                    pbndsf -= *li;
-                  }
-                }
-                pbndsf.normalize();
-              }
-              // Buffer zones
-              ibset buffers;
-              {
-                for (typename ibset::const_iterator pbi=pbndsf.begin();
-                     pbi!=pbndsf.end(); ++pbi) {
-                  buffers |= (*pbi).expand(buffer_width, buffer_width) & extrf;
-                }
-                buffers.normalize();
-              }
-              // Add boundaries
-              const ibbox maxrecvs
-                = extr.expand(-pss,-pss).contracted_for(extrf);
-              ibset recvs = buffers & maxrecvs;
-              recvs.normalize();
-              {
-                // Do not prolongate what is already prolongated
-                const iblistvect& rrbc
-                  = boxes.at(rl+1).at(cc).at(ml).recv_ref_bnd_coarse;
-                for (typename iblistvect::const_iterator lvi=rrbc.begin();
-                     lvi!=rrbc.end(); ++lvi) {
-                  for (typename iblist::const_iterator li=lvi->begin();
-                       li!=lvi->end(); ++li) {
-                    recvs -= *li;
-                  }
-                }
-                recvs.normalize();
-              }
-              {
-                for (typename ibset::const_iterator ri = recvs.begin();
-                     ri != recvs.end(); ++ri) {
-                  const ibbox & recv = *ri;
-                  const ibbox send = recv.expanded_for(extr);
-                  assert (! send.empty());
-                  assert (send.is_contained_in(extr));
-                  assert (send.is_contained_in(extr.expand(-pss,-pss)));
-                  boxes.at(rl+1).at(cc).at(ml).recv_ref_bnd_coarse.at(c )
-                    .push_back(recv);
-                  boxes.at(rl  ).at(c ).at(ml).send_ref_bnd_fine  .at(cc)
-                    .push_back(send);
-                }
-      	      }
-      	    }
-            
-      	  } // for cc
-      	} // if not finest refinement level
-	
-      } // for ml
-    } // for c
-  } // for rl
-  
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    for (int c=0; c<h.components(rl); ++c) {
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-      	const ibbox& intr = boxes.at(rl).at(c).at(ml).interior;
-      	const ibbox& extr = boxes.at(rl).at(c).at(ml).exterior;
-	
-      	// Refinement boxes
-      	if (rl<h.reflevels()-1) {
-      	  for (int cc=0; cc<h.components(rl+1); ++cc) {
-      	    const ibbox intrf = boxes.at(rl+1).at(cc).at(ml).interior;
-      	    // Restriction (interior)
-      	    {
-      	      // (the restriction may fill the interior of the of the
-      	      // coarse grid, and may use the interior of the fine
-      	      // grid, and the bbox must be as large as possible)
-              // (the restriction must not use points that are filled
-              // by boundary prolongation)
-              ibset sends = intrf & intr.expanded_for(intrf);
-              // remove what is received during boundary prolongation
-              for (int ccc=0; ccc<h.components(rl); ++ccc) {
-                const iblist& sendlist
-                  = boxes.at(rl+1).at(cc).at(ml).recv_ref_bnd_coarse.at(ccc);
-                for (typename iblist::const_iterator sli = sendlist.begin();
-                     sli != sendlist.end(); ++sli) {
-                  sends -= *sli;
-                }
-              }
-              sends.normalize();
-              for (typename ibset::const_iterator si = sends.begin();
-                   si != sends.end(); ++si) {
-                const ibbox recv = (*si).contracted_for(intr);
-                if (! recv.empty()) {
-                  const ibbox & send = recv.expanded_for(intrf);
-                  assert (! send.empty());
-                  boxes.at(rl+1).at(cc).at(ml).send_ref_coarse.at(c )
-                    .push_back(send);
-                  boxes.at(rl  ).at(c ).at(ml).recv_ref_fine  .at(cc)
-                    .push_back(recv);
-                }
-              }
-      	    }
-            
-      	  } // for cc
-      	} // if not finest refinement level
-	
-      } // for ml
-    } // for c
-  } // for rl
-  
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    for (int c=0; c<h.components(rl); ++c) {
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-	
-	// Boundaries that are not synced, or are neither synced nor
-	// prolonged to from coarser grids (outer boundaries)
-        ibset& sync_not = boxes.at(rl).at(c).at(ml).sync_not;
-        ibset& recv_not = boxes.at(rl).at(c).at(ml).recv_not;
-        
-        // The whole boundary
-        sync_not = boxes.at(rl).at(c).at(ml).boundaries;
-        recv_not = boxes.at(rl).at(c).at(ml).boundaries;
-        
-        // Subtract boxes received during synchronisation
-        const iblistvect& recv_sync = boxes.at(rl).at(c).at(ml).recv_sync;
-        for (typename iblistvect::const_iterator lvi=recv_sync.begin();
-             lvi!=recv_sync.end(); ++lvi) {
+      bbox.recv_mg_coarse.push_back(recv);
+      box .send_mg_fine  .push_back(send);
+    }
+  } // if not finest multigrid level
+}
+
+template<int D>
+void dh<D>::setup_sync_and_refine_boxes( dh<D>::dboxes & b, int rl, int c, int ml )
+{
+  // Sync boxes
+  const int cs = h.components(rl);
+  b.send_sync.resize(cs);
+  b.recv_sync.resize(cs);
+
+  // Refinement boxes
+  if (rl>0) {
+    const int csm1 = h.components(rl-1);
+    b.send_ref_coarse.resize(csm1);
+    b.recv_ref_coarse.resize(csm1);
+    b.recv_ref_bnd_coarse.resize(csm1);
+  }
+  if (rl<h.reflevels()-1) {
+    const int csp1 = h.components(rl+1);
+    b.recv_ref_fine.resize(csp1);
+    b.send_ref_fine.resize(csp1);
+    b.send_ref_bnd_fine.resize(csp1);
+  }
+}
+
+template<int D>
+void dh<D>::intersect_sync_with_interior( dh<D>::dboxes & box, int rl, int c, int ml )
+{
+  const ibset& bnds = box.boundaries;
+
+  // Sync boxes
+  for (int cc=0; cc<h.components(rl); ++cc) {
+    dboxes & box1 = boxes.at(rl).at(cc).at(ml);
+    assert (ml<h.mglevels(rl,cc));
+    // intersect boundaries with interior of that component
+    ibset ovlp = bnds & box1.interior;
+    ovlp.normalize();
+    for (typename ibset::const_iterator b=ovlp.begin();b!=ovlp.end(); ++b) {
+      box .recv_sync.at(cc).push_back(*b);
+      box1.send_sync.at(c ).push_back(*b);
+    }
+  }
+}
+
+template<int D>
+void dh<D>::setup_refinement_interior_boxes( dh<D>::dboxes & box, int rl, int c, int ml )
+{
+  const ibbox& intr = box.interior;
+  const ibbox& extr = box.exterior;
+
+  // Refinement boxes
+  if (rl<h.reflevels()-1) {
+    for (int cc=0; cc<h.components(rl+1); ++cc) {
+      dboxes & box1 = boxes.at(rl).at(cc).at(ml);
+      const ibbox intrf = box1.interior;
+      // Prolongation (interior)
+      // TODO: prefer boxes from the same processor
+      {
+        // (the prolongation may use the exterior of the coarse
+        // grid, and must fill all of the interior of the fine
+        // grid)
+        const int pss = prolongation_stencil_size();
+        ibset recvs = extr.expand(-pss,-pss).contracted_for(intrf) & intrf;
+        const iblistvect& rrc = box1.recv_ref_coarse;
+        for (typename iblistvect::const_iterator lvi=rrc.begin();
+             lvi!=rrc.end(); ++lvi) {
           for (typename iblist::const_iterator li=lvi->begin();
                li!=lvi->end(); ++li) {
-            sync_not -= *li;
-            recv_not -= *li;
+            recvs -= *li;
           }
         }
-        
-        // Subtract boxes received during prolongation
-        const iblistvect& recv_ref_bnd_coarse
-          = boxes.at(rl).at(c).at(ml).recv_ref_bnd_coarse;
-        for (typename iblistvect::const_iterator
-               lvi=recv_ref_bnd_coarse.begin();
-             lvi!=recv_ref_bnd_coarse.end(); ++lvi) {
-          for (typename iblist::const_iterator li=lvi->begin();
-               li!=lvi->end(); ++li) {
-            recv_not -= *li;
+        recvs.normalize();
+        assert (recvs.setsize() <= 1);
+        if (recvs.setsize() == 1) {
+          const ibbox recv = *recvs.begin();
+          const ibbox send = recv.expanded_for(extr);
+          assert (! send.empty());
+          assert (send.is_contained_in(extr));
+          box1.recv_ref_coarse.at(c ).push_back(recv);
+          box. send_ref_fine  .at(cc).push_back(send);
+        }
+      }
+    } // for cc
+  } // if not finest refinement level
+}
+
+template<int D>
+void dh<D>::setup_refinement_exterior_boxes( dh<D>::dboxes & box, int rl, int c, int ml )
+{
+  const ibbox& intr = box.interior;
+  const ibbox& extr = box.exterior;
+
+  // Refinement boxes
+  if (rl<h.reflevels()-1) {
+    for (int cc=0; cc<h.components(rl+1); ++cc) {
+      dboxes & box1 = boxes.at(rl+1).at(cc).at(ml);
+      const ibbox intrf = box1.interior;
+            const ibbox& extrf = box1.exterior;
+            const ibset& bndsf = box1.boundaries;
+      // Prolongation (boundaries)
+            // TODO: prefer boxes from the same processor
+      {
+        // (the prolongation may use the exterior of the coarse
+        // grid, and must fill all of the boundary of the fine
+        // grid)
+        const int pss = prolongation_stencil_size();
+        // Prolongation boundaries
+        ibset pbndsf = bndsf;
+        {
+          // Do not count what is synced
+          const iblistvect& rs = box1.recv_sync;
+          for (typename iblistvect::const_iterator lvi=rs.begin();
+               lvi!=rs.end(); ++lvi) {
+            for (typename iblist::const_iterator li=lvi->begin();
+                 li!=lvi->end(); ++li) {
+              pbndsf -= *li;
+            }
+          }
+          pbndsf.normalize();
+        }
+        // Buffer zones
+        ibset buffers;
+        {
+          for (typename ibset::const_iterator pbi=pbndsf.begin();
+               pbi!=pbndsf.end(); ++pbi) {
+            buffers |= (*pbi).expand(buffer_width, buffer_width) & extrf;
+          }
+          buffers.normalize();
+        }
+        // Add boundaries
+        const ibbox maxrecvs = extr.expand(-pss,-pss).contracted_for(extrf);
+        ibset recvs = buffers & maxrecvs;
+        recvs.normalize();
+        {
+          // Do not prolongate what is already prolongated
+          const iblistvect& rrbc = box1.recv_ref_bnd_coarse;
+          for (typename iblistvect::const_iterator lvi=rrbc.begin();
+               lvi!=rrbc.end(); ++lvi) {
+            for (typename iblist::const_iterator li=lvi->begin();
+                 li!=lvi->end(); ++li) {
+              recvs -= *li;
+            }
+          }
+          recvs.normalize();
+        }
+        {
+          for (typename ibset::const_iterator ri = recvs.begin();
+               ri != recvs.end(); ++ri) {
+            const ibbox & recv = *ri;
+            const ibbox send = recv.expanded_for(extr);
+            assert (! send.empty());
+            assert (send.is_contained_in(extr));
+            assert (send.is_contained_in(extr.expand(-pss,-pss)));
+            box1.recv_ref_bnd_coarse.at(c ).push_back(recv);
+            box .send_ref_bnd_fine  .at(cc).push_back(send);
           }
         }
-	
-      } // for ml
-    } // for c
-  } // for rl
-  
-  // Calculate bases
-  bases.resize(h.reflevels());
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    if (h.components(rl)==0) {
-      bases.at(rl).resize(0);
-    } else {
-      bases.at(rl).resize(h.mglevels(rl,0));
-      for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
-	bases.at(rl).at(ml).exterior = ibbox();
-	bases.at(rl).at(ml).interior = ibbox();
-	for (int c=0; c<h.components(rl); ++c) {
-	  bases.at(rl).at(ml).exterior
-	    = (bases.at(rl).at(ml).exterior
-	       .expanded_containing(boxes.at(rl).at(c).at(ml).exterior));
-	  bases.at(rl).at(ml).interior
-	    = (bases.at(rl).at(ml).interior
-	       .expanded_containing(boxes.at(rl).at(c).at(ml).interior));
-	}
-	bases.at(rl).at(ml).boundaries
-	  = bases.at(rl).at(ml).exterior - bases.at(rl).at(ml).interior;
       }
+            
+    } // for cc
+  } // if not finest refinement level
+}
+
+template<int D>
+void dh<D>::setup_restrict_interior_boxes( dh<D>::dboxes & box, int rl, int c, int ml )
+{
+  const ibbox& intr = box.interior;
+  const ibbox& extr = box.exterior;
+
+  // Refinement boxes
+  if (rl<h.reflevels()-1) {
+    for (int cc=0; cc<h.components(rl+1); ++cc) {
+      dboxes & box1 = boxes.at(rl+1).at(cc).at(ml);
+      const ibbox intrf = box1.interior;
+      // Restriction (interior)
+      {
+        // (the restriction may fill the interior of the of the
+        // coarse grid, and may use the interior of the fine
+        // grid, and the bbox must be as large as possible)
+        // (the restriction must not use points that are filled
+        // by boundary prolongation)
+        ibset sends = intrf & intr.expanded_for(intrf);
+        // remove what is received during boundary prolongation
+        for (int ccc=0; ccc<h.components(rl); ++ccc) {
+          const iblist& sendlist = box1.recv_ref_bnd_coarse.at(ccc);
+          for (typename iblist::const_iterator sli = sendlist.begin();
+               sli != sendlist.end(); ++sli) {
+            sends -= *sli;
+          }
+        }
+        sends.normalize();
+        for (typename ibset::const_iterator si = sends.begin();
+             si != sends.end(); ++si) {
+          const ibbox recv = (*si).contracted_for(intr);
+          if (! recv.empty()) {
+            const ibbox & send = recv.expanded_for(intrf);
+            assert (! send.empty());
+            box1.send_ref_coarse.at(c ).push_back(send);
+            box .recv_ref_fine  .at(cc).push_back(recv);
+          }
+        }
+      }
+            
+    } // for cc
+  } // if not finest refinement level
+}
+
+template<int D>
+void dh<D>::trim_unsynced_boundaries( dh<D>::dboxes & box, int rl, int c, int ml )
+{
+  // Boundaries that are not synced, or are neither synced nor
+  // prolonged to from coarser grids (outer boundaries)
+  ibset& sync_not = box.sync_not;
+  ibset& recv_not = box.recv_not;
+  
+  // The whole boundary
+  sync_not = box.boundaries;
+  recv_not = box.boundaries;
+  
+  // Subtract boxes received during synchronisation
+  const iblistvect& recv_sync = box.recv_sync;
+  for (typename iblistvect::const_iterator lvi=recv_sync.begin();
+       lvi!=recv_sync.end(); ++lvi) {
+    for (typename iblist::const_iterator li=lvi->begin();
+         li!=lvi->end(); ++li) {
+      sync_not -= *li;
+      recv_not -= *li;
     }
   }
   
-  if (output_bboxes) {
-    cout << endl << h << endl;
-    for (int rl=0; rl<h.reflevels(); ++rl) {
-      for (int c=0; c<h.components(rl); ++c) {
-	for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-	  cout << endl;
-	  cout << "dh bboxes:" << endl;
-	  cout << "rl=" << rl << " c=" << c << " ml=" << ml << endl;
-	  cout << "exterior=" << boxes.at(rl).at(c).at(ml).exterior << endl;
-	  cout << "interior=" << boxes.at(rl).at(c).at(ml).interior << endl;
-	  cout << "send_mg_fine=" << boxes.at(rl).at(c).at(ml).send_mg_fine << endl;
-	  cout << "send_mg_coarse=" << boxes.at(rl).at(c).at(ml).send_mg_coarse << endl;
-	  cout << "recv_mg_fine=" << boxes.at(rl).at(c).at(ml).recv_mg_fine << endl;
-	  cout << "recv_mg_coarse=" << boxes.at(rl).at(c).at(ml).recv_mg_coarse << endl;
-	  cout << "send_ref_fine=" << boxes.at(rl).at(c).at(ml).send_ref_fine << endl;
-	  cout << "send_ref_coarse=" << boxes.at(rl).at(c).at(ml).send_ref_coarse << endl;
-	  cout << "recv_ref_fine=" << boxes.at(rl).at(c).at(ml).recv_ref_fine << endl;
-	  cout << "recv_ref_coarse=" << boxes.at(rl).at(c).at(ml).recv_ref_coarse << endl;
-	  cout << "send_sync=" << boxes.at(rl).at(c).at(ml).send_sync << endl;
-	  cout << "send_ref_bnd_fine=" << boxes.at(rl).at(c).at(ml).send_ref_bnd_fine << endl;
-	  cout << "boundaries=" << boxes.at(rl).at(c).at(ml).boundaries << endl;
-	  cout << "recv_sync=" << boxes.at(rl).at(c).at(ml).recv_sync << endl;
-	  cout << "recv_ref_bnd_coarse=" << boxes.at(rl).at(c).at(ml).recv_ref_bnd_coarse << endl;
-	  cout << "sync_not=" << boxes.at(rl).at(c).at(ml).sync_not << endl;
-	  cout << "recv_not=" << boxes.at(rl).at(c).at(ml).recv_not << endl;
-	}
-      }
-    }
-    for (int rl=0; rl<h.reflevels(); ++rl) {
-      if (h.components(rl)>0) {
-	for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
-	  cout << endl;
-	  cout << "dh bases:" << endl;
-	  cout << "rl=" << rl << " ml=" << ml << endl;
-	  cout << "exterior=" << bases.at(rl).at(ml).exterior << endl;
-	  cout << "interior=" << bases.at(rl).at(ml).interior << endl;
-	  cout << "boundaries=" << bases.at(rl).at(ml).boundaries << endl;
-	}
-      }
+  // Subtract boxes received during prolongation
+  const iblistvect& recv_ref_bnd_coarse = box.recv_ref_bnd_coarse;
+  for (typename iblistvect::const_iterator lvi=recv_ref_bnd_coarse.begin();
+       lvi!=recv_ref_bnd_coarse.end(); ++lvi) {
+    for (typename iblist::const_iterator li=lvi->begin();
+         li!=lvi->end(); ++li) {
+      recv_not -= *li;
     }
-  } // if output_bboxes
-  
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    for (int c=0; c<h.components(rl); ++c) {
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-	
-	// Assert that all boundaries are synced or received
-        {
-          const ibset& sync_not = boxes.at(rl).at(c).at(ml).sync_not;
+  }
+}
+
+template<int D>
+void dh<D>::assert_assert_assert( dh<D>::dboxes & box, int rl, int c, int ml )
+{
+// Assert that all boundaries are synced or received
+  {
+    const ibset& sync_not = box.sync_not;
 #if 0
-          const ibset& recv_not = boxes.at(rl).at(c).at(ml).recv_not;
+    const ibset& recv_not = boxes.at(rl).at(c).at(ml).recv_not;
 #endif
           
-          // Check that no boundaries are left over
-          if (rl==0) assert (sync_not.empty());
+    // Check that no boundaries are left over
+    if (rl==0) assert (sync_not.empty());
 #if 0
-          assert (recv_not.empty());
+    assert (recv_not.empty());
 #endif
+  }
+
+  // Assert that the interior is received exactly once during
+  // prolongation, and that nothing else is received
+  {
+    if (rl==0) {
+      const iblistvect& recv_ref_coarse = box.recv_ref_coarse;
+      assert (recv_ref_coarse.empty());
+    } else {              // rl!=0
+      const iblistvect& recv_ref_coarse = box.recv_ref_coarse;
+      ibset intr = box.interior;
+      for (typename iblistvect::const_iterator
+             lvi=recv_ref_coarse.begin();
+           lvi!=recv_ref_coarse.end(); ++lvi) {
+        for (typename iblist::const_iterator li=lvi->begin();
+             li!=lvi->end(); ++li) {
+          const int old_sz = intr.size();
+          const int this_sz = li->size();
+          intr -= *li;
+          const int new_sz = intr.size();
+          // TODO
+          assert (new_sz + this_sz == old_sz);
         }
-	
-        // Assert that the interior is received exactly once during
-        // prolongation, and that nothing else is received
-        {
-          if (rl==0) {
-            const iblistvect& recv_ref_coarse
-              = boxes.at(rl).at(c).at(ml).recv_ref_coarse;
-            assert (recv_ref_coarse.empty());
-          } else {              // rl!=0
-            const iblistvect& recv_ref_coarse
-              = boxes.at(rl).at(c).at(ml).recv_ref_coarse;
-            ibset intr = boxes.at(rl).at(c).at(ml).interior;
-            for (typename iblistvect::const_iterator
-                   lvi=recv_ref_coarse.begin();
-                 lvi!=recv_ref_coarse.end(); ++lvi) {
-              for (typename iblist::const_iterator li=lvi->begin();
-                   li!=lvi->end(); ++li) {
-                const int old_sz = intr.size();
-                const int this_sz = li->size();
-                intr -= *li;
-                const int new_sz = intr.size();
-                // TODO
-                assert (new_sz + this_sz == old_sz);
-              }
-            }
-            // TODO
-            // This need not be empty at outer boundaries.  Check that
-            // those are indeed outer boundaries!  But what size of the
-            // boundary region should be used for that?
+      }
+      // TODO
+      // This need not be empty at outer boundaries.  Check that
+      // those are indeed outer boundaries!  But what size of the
+      // boundary region should be used for that?
 #if 0
-            assert (intr.empty());
+      assert (intr.empty());
 #endif
-          }
-        }
-        
-        // Assert that the boundaries are received at most once during
-        // prolongation and synchronisation, and that nothing else is
-        // received
-        {
-          const iblistvect& recv_sync = boxes.at(rl).at(c).at(ml).recv_sync;
-          const iblistvect& recv_ref_bnd_coarse
-            = boxes.at(rl).at(c).at(ml).recv_ref_bnd_coarse;
-          ibset bnds = boxes.at(rl).at(c).at(ml).boundaries;
-          for (typename iblistvect::const_iterator lvi=recv_sync.begin();
-               lvi!=recv_sync.end(); ++lvi) {
-            for (typename iblist::const_iterator li=lvi->begin();
-                 li!=lvi->end(); ++li) {
-              const int old_sz = bnds.size();
-              const int this_sz = li->size();
-              bnds -= *li;
-              const int new_sz = bnds.size();
-              assert (new_sz + this_sz == old_sz);
-            }
-          }
-          for (typename iblistvect::const_iterator
-                 lvi=recv_ref_bnd_coarse.begin();
-               lvi!=recv_ref_bnd_coarse.end(); ++lvi) {
-            for (typename iblist::const_iterator li=lvi->begin();
-                 li!=lvi->end(); ++li) {
-              const int old_sz = bnds.size();
-              const int this_sz = li->size();
-              bnds -= *li;
-              const int new_sz = bnds.size();
-              // TODO
-              // The new size can be larger if part of the
-              // prolongation went into the buffer zone.
-//               assert (new_sz + this_sz == old_sz);
-              assert (new_sz + this_sz >= old_sz);
-            }
-          }
-          // TODO
-          // This need not be empty at outer boundaries.  Check that
-          // those are indeed outer boundaries!  But what size of the
-          // boundary region should be used for that?
+    }
+  }
+  
+  // Assert that the boundaries are received at most once during
+  // prolongation and synchronisation, and that nothing else is
+  // received
+  {
+    const iblistvect& recv_sync = box.recv_sync;
+    const iblistvect& recv_ref_bnd_coarse = box.recv_ref_bnd_coarse;
+    ibset bnds = box.boundaries;
+    for (typename iblistvect::const_iterator lvi=recv_sync.begin();
+         lvi!=recv_sync.end(); ++lvi) {
+      for (typename iblist::const_iterator li=lvi->begin();
+           li!=lvi->end(); ++li) {
+        const int old_sz = bnds.size();
+        const int this_sz = li->size();
+        bnds -= *li;
+        const int new_sz = bnds.size();
+        assert (new_sz + this_sz == old_sz);
+      }
+    }
+    for (typename iblistvect::const_iterator
+           lvi=recv_ref_bnd_coarse.begin();
+         lvi!=recv_ref_bnd_coarse.end(); ++lvi) {
+      for (typename iblist::const_iterator li=lvi->begin();
+           li!=lvi->end(); ++li) {
+        const int old_sz = bnds.size();
+        const int this_sz = li->size();
+        bnds -= *li;
+        const int new_sz = bnds.size();
+        // TODO
+        // The new size can be larger if part of the
+        // prolongation went into the buffer zone.
+//      assert (new_sz + this_sz == old_sz);
+        assert (new_sz + this_sz >= old_sz);
+      }
+    }
+    // TODO
+    // This need not be empty at outer boundaries.  Check that
+    // those are indeed outer boundaries!  But what size of the
+    // boundary region should be used for that?
 #if 0
-          assert (bnds.empty());
+    assert (bnds.empty());
 #endif
+  }
+}
+
+template<int D>
+void dh<D>::calculate_bases () {
+  // Calculate bases
+  bases.resize(h.reflevels());
+  for (int rl=0; rl<h.reflevels(); ++rl) {
+    if (h.components(rl)==0) {
+      bases.at(rl).resize(0);
+    } else {
+      bases.at(rl).resize(h.mglevels(rl,0));
+      for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
+        dbases & base = bases.at(rl).at(ml);
+        base.exterior = ibbox();
+        base.interior = ibbox();
+        for (int c=0; c<h.components(rl); ++c) {
+          dboxes & box = boxes.at(rl).at(c).at(ml);
+          base.exterior = (base.exterior.expanded_containing(box.exterior));
+          base.interior = (base.interior.expanded_containing(box.interior));
         }
-        
-      } // for ml
-    } // for c
-  } // for rl
-  
-  if (! save_memory_during_regridding) {
-    
+        base.boundaries = base.exterior - base.interior;
+      }
+    }
+  }
+}
+
+template<int D>
+void dh<D>::save_time ( bool do_prolongate ) {
+  for (typename list<ggf<D>*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+    (*f)->recompose_crop ();
+  }
+  for (int rl=0; rl<h.reflevels(); ++rl) {
     for (typename list<ggf<D>*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-      (*f)->recompose_crop ();
+      (*f)->recompose_allocate (rl);
     }
-    for (int rl=0; rl<h.reflevels(); ++rl) {
+    for (comm_state<D> state; !state.done(); state.step()) {
       for (typename list<ggf<D>*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-        (*f)->recompose_allocate (rl);
+        (*f)->recompose_fill (state, rl, do_prolongate);
       }
+    }
+    for (typename list<ggf<D>*>::reverse_iterator f=gfs.rbegin(); f!=gfs.rend(); ++f) {
+      (*f)->recompose_free (rl);
+    }
+    for (comm_state<D> state; !state.done(); state.step()) {
+      for (typename list<ggf<D>*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+        (*f)->recompose_bnd_prolongate (state, rl, do_prolongate);
+      }
+    }
+    for (comm_state<D> state; !state.done(); state.step()) {
+      for (typename list<ggf<D>*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+        (*f)->recompose_sync (state, rl, do_prolongate);
+      }
+    }
+  } // for rl
+}
+
+template<int D>
+void dh<D>::save_memory ( bool do_prolongate ) {
+  ggf<D>* vectorleader = NULL;
+  for (typename list<ggf<D>*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+    
+    (*f)->recompose_crop ();
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      (*f)->recompose_allocate (rl);
       for (comm_state<D> state; !state.done(); state.step()) {
-        for (typename list<ggf<D>*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-          (*f)->recompose_fill (state, rl, do_prolongate);
-        }
+        (*f)->recompose_fill (state, rl, do_prolongate);
       }
-      for (typename list<ggf<D>*>::reverse_iterator f=gfs.rbegin(); f!=gfs.rend(); ++f) {
+      assert ((*f)->vectorlength >= 1);
+      if ((*f)->vectorlength == 1) {
+        assert (! vectorleader);
         (*f)->recompose_free (rl);
+      } else {
+        // treat vector groups specially: delete the leader only
+        // after all other elements have been deleted
+        if ((*f)->vectorindex == 0) {
+          // first element: delete nothing
+          if (rl == 0) {
+            assert (! vectorleader);
+            vectorleader = *f;
+          }
+          assert (vectorleader);
+        } else {
+          assert (vectorleader);
+          (*f)->recompose_free (rl);
+          if ((*f)->vectorindex == (*f)->vectorlength-1) {
+            // this was the last element: delete the leader as well
+            vectorleader->recompose_free (rl);
+            if (rl == h.reflevels()-1) {
+              vectorleader = NULL;
+            }
+          }
+        }
       }
       for (comm_state<D> state; !state.done(); state.step()) {
-        for (typename list<ggf<D>*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-          (*f)->recompose_bnd_prolongate (state, rl, do_prolongate);
-        }
+        (*f)->recompose_bnd_prolongate (state, rl, do_prolongate);
       }
       for (comm_state<D> state; !state.done(); state.step()) {
-        for (typename list<ggf<D>*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-          (*f)->recompose_sync (state, rl, do_prolongate);
-        }
+        (*f)->recompose_sync (state, rl, do_prolongate);
       }
     } // for rl
     
-  } else {                      // save memory
-    
-    ggf<D>* vectorleader = NULL;
-    for (typename list<ggf<D>*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-      
-      (*f)->recompose_crop ();
-      for (int rl=0; rl<h.reflevels(); ++rl) {
-        (*f)->recompose_allocate (rl);
-        for (comm_state<D> state; !state.done(); state.step()) {
-          (*f)->recompose_fill (state, rl, do_prolongate);
-        }
-        assert ((*f)->vectorlength >= 1);
-        if ((*f)->vectorlength == 1) {
-          assert (! vectorleader);
-          (*f)->recompose_free (rl);
-        } else {
-          // treat vector groups specially: delete the leader only
-          // after all other elements have been deleted
-          if ((*f)->vectorindex == 0) {
-            // first element: delete nothing
-            if (rl == 0) {
-              assert (! vectorleader);
-              vectorleader = *f;
-            }
-            assert (vectorleader);
-          } else {
-            assert (vectorleader);
-            (*f)->recompose_free (rl);
-            if ((*f)->vectorindex == (*f)->vectorlength-1) {
-              // this was the last element: delete the leader as well
-              vectorleader->recompose_free (rl);
-              if (rl == h.reflevels()-1) {
-                vectorleader = NULL;
-              }
-            }
-          }
-        }
-        for (comm_state<D> state; !state.done(); state.step()) {
-          (*f)->recompose_bnd_prolongate (state, rl, do_prolongate);
-        }
-        for (comm_state<D> state; !state.done(); state.step()) {
-          (*f)->recompose_sync (state, rl, do_prolongate);
-        }
-      } // for rl
-      
-    } // for gf
-    assert (! vectorleader);
-    
-  } // save memory
+  } // for gf
+  assert (! vectorleader);
 }
 
-
-
 // Grid function management
 template<int D>
 void dh<D>::add (ggf<D>* f) {
@@ -684,6 +644,47 @@ void dh<D>::output (ostream& os) const {
   os << "}";
 }
 
+template<int D>
+void dh<D>::output_bboxes( dh<D>::dboxes & box, int rl, int c, int ml )
+{
+  cout << endl;
+  cout << "dh bboxes:" << endl;
+  cout << "rl=" << rl << " c=" << c << " ml=" << ml << endl;
+  cout << "exterior=" << box.exterior << endl;
+  cout << "interior=" << box.interior << endl;
+  cout << "send_mg_fine=" << box.send_mg_fine << endl;
+  cout << "send_mg_coarse=" << box.send_mg_coarse << endl;
+  cout << "recv_mg_fine=" << box.recv_mg_fine << endl;
+  cout << "recv_mg_coarse=" << box.recv_mg_coarse << endl;
+  cout << "send_ref_fine=" << box.send_ref_fine << endl;
+  cout << "send_ref_coarse=" << box.send_ref_coarse << endl;
+  cout << "recv_ref_fine=" << box.recv_ref_fine << endl;
+  cout << "recv_ref_coarse=" << box.recv_ref_coarse << endl;
+  cout << "send_sync=" << box.send_sync << endl;
+  cout << "send_ref_bnd_fine=" << box.send_ref_bnd_fine << endl;
+  cout << "boundaries=" << box.boundaries << endl;
+  cout << "recv_sync=" << box.recv_sync << endl;
+  cout << "recv_ref_bnd_coarse=" << box.recv_ref_bnd_coarse << endl;
+  cout << "sync_not=" << box.sync_not << endl;
+  cout << "recv_not=" << box.recv_not << endl;
+}
+
+template<int D>
+void dh<D>::output_bases () {
+  for (int rl=0; rl<h.reflevels(); ++rl) {
+    if (h.components(rl)>0) {
+      for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
+        dbases & base = bases.at(rl).at(ml);
+        cout << endl;
+        cout << "dh bases:" << endl;
+        cout << "rl=" << rl << " ml=" << ml << endl;
+        cout << "exterior=" << base.exterior << endl;
+        cout << "interior=" << base.interior << endl;
+        cout << "boundaries=" << base.boundaries << endl;
+      }
+    }
+  }
+}
 
 
 template class dh<3>;