1 files changed, 45 insertions, 32 deletions

diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index 9f92c2a..9ef1ce8 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -313,45 +313,58 @@ public:
     }
   }
 
-  Grid<T,DESCRIPTOR>* refine(Vector<T,2> origin, Vector<T,2> extend)
+  FineCoupler<T,DESCRIPTOR>& addFineCoupling(
+    Grid<T,DESCRIPTOR>& fineGrid, Vector<T,2> origin, Vector<T,2> extend)
+  {
+    _fineCouplers.emplace_back(
+      new FineCoupler<T,DESCRIPTOR>(
+        *this, fineGrid, origin, extend));
+    return *_fineCouplers.back();
+  }
+
+  CoarseCoupler<T,DESCRIPTOR>& addCoarseCoupling(
+    Grid<T,DESCRIPTOR>& fineGrid, Vector<T,2> origin, Vector<T,2> extend)
+  {
+    _coarseCouplers.emplace_back(
+      new CoarseCoupler<T,DESCRIPTOR>(
+        *this, fineGrid, origin, extend));
+    return *_coarseCouplers.back();
+  }
+
+  Grid<T,DESCRIPTOR>& refine(IndicatorF2D<T>& domainF)
   {
-    IndicatorCuboid2D<T> fineCuboid(extend, origin);
     _fineGrids.emplace_back(
       new Grid<T,DESCRIPTOR>(
-        fineCuboid,
+        domainF,
         2*getConverter().getResolution(),
         2.0*getConverter().getLatticeRelaxationTime() - 0.5,
         getConverter().getReynoldsNumber()
       ));
-    Grid<T,DESCRIPTOR>* const fineGrid = _fineGrids.back().get();
+    return *_fineGrids.back();
+  }
 
-    const T coarseDeltaX = getConverter().getPhysDeltaX();
+  Grid<T,DESCRIPTOR>& refine(Vector<T,2> origin, Vector<T,2> extend)
+  {
+    IndicatorCuboid2D<T> fineCuboid(extend, origin);
+    auto& fineGrid = refine(fineCuboid);
 
-    _fineCouplers.emplace_back(
-      new FineCoupler<T,DESCRIPTOR>(
-        *this, *fineGrid, origin, Vector<T,2> {0,extend[1]}));
-    _fineCouplers.emplace_back(
-      new FineCoupler<T,DESCRIPTOR>(
-        *this, *fineGrid, origin + Vector<T,2> {extend[0],0}, Vector<T,2> {0,extend[1]}));
-    _fineCouplers.emplace_back(
-      new FineCoupler<T,DESCRIPTOR>(
-        *this, *fineGrid, origin + Vector<T,2>(0,extend[1]), Vector<T,2> {extend[0],0}));
-    _fineCouplers.emplace_back(
-      new FineCoupler<T,DESCRIPTOR>(
-        *this, *fineGrid, origin, Vector<T,2> {extend[0],0}));
+    const Vector<T,2> extendX = {extend[0],0};
+    const Vector<T,2> extendY = {0,extend[1]};
 
-    _coarseCouplers.emplace_back(
-      new CoarseCoupler<T,DESCRIPTOR>(
-        *this, *fineGrid, origin + coarseDeltaX, Vector<T,2> {0,extend[1]-2*coarseDeltaX}));
-    _coarseCouplers.emplace_back(
-      new CoarseCoupler<T,DESCRIPTOR>(
-        *this, *fineGrid, origin + Vector<T,2>(extend[0]-coarseDeltaX,coarseDeltaX), Vector<T,2> {0,extend[1]-2*coarseDeltaX}));
-    _coarseCouplers.emplace_back(
-      new CoarseCoupler<T,DESCRIPTOR>(
-        *this, *fineGrid, origin + Vector<T,2> {coarseDeltaX,extend[1]-coarseDeltaX}, Vector<T,2> {extend[0]-2*coarseDeltaX,0}));
-    _coarseCouplers.emplace_back(
-      new CoarseCoupler<T,DESCRIPTOR>(
-        *this, *fineGrid, origin + coarseDeltaX, Vector<T,2> {extend[0]-2*coarseDeltaX,0}));
+    addFineCoupling(fineGrid, origin,           extendY);
+    addFineCoupling(fineGrid, origin + extendX, extendY);
+    addFineCoupling(fineGrid, origin + extendY, extendX);
+    addFineCoupling(fineGrid, origin,           extendX);
+
+    const T coarseDeltaX = getConverter().getPhysDeltaX();
+    const Vector<T,2> innerOrigin = origin + coarseDeltaX;
+    const Vector<T,2> innerExtendX = extendX - Vector<T,2> {2*coarseDeltaX,0};
+    const Vector<T,2> innerExtendY = extendY - Vector<T,2> {0,2*coarseDeltaX};
+
+    addCoarseCoupling(fineGrid, innerOrigin,                innerExtendY);
+    addCoarseCoupling(fineGrid, innerOrigin + innerExtendX, innerExtendY);
+    addCoarseCoupling(fineGrid, innerOrigin + innerExtendY, innerExtendX);
+    addCoarseCoupling(fineGrid, innerOrigin,                innerExtendX);
 
     return fineGrid;
   }
@@ -684,8 +697,6 @@ int main(int argc, char* argv[])
   singleton::directories().setOutputDir("./tmp/");
   OstreamManager clout(std::cout,"main");
 
-  const T coarseDeltaX = 1./N;
-
   Vector<T,2> coarseOrigin { 0.0, 0.0 };
   Vector<T,2> coarseExtend { lx, ly  };
   IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin);
@@ -694,11 +705,13 @@ int main(int argc, char* argv[])
   Vector<T,2> fineExtend   { 0.6  * lx, 0.6 };
 
   auto coarseGrid = Grid<T,DESCRIPTOR>::make(coarseCuboid, N, 0.8, Re);
-  auto fineGrid   = coarseGrid->refine(fineOrigin, fineExtend);
+  auto fineGrid   = &coarseGrid->refine(fineOrigin, fineExtend);
 
   prepareGeometry(coarseGrid->getConverter(), coarseGrid->getSuperGeometry());
   prepareGeometry(fineGrid->getConverter(), fineGrid->getSuperGeometry());
 
+  const T coarseDeltaX = coarseGrid->getConverter().getPhysDeltaX();
+
   fineGrid->getSuperGeometry().rename(2,1);
   fineGrid->getSuperGeometry().rename(5,2);
   IndicatorCuboid2D<T> overlapCuboid(fineExtend - 4*coarseDeltaX, fineOrigin + 2*coarseDeltaX);