Further indicatorize geometry setup

Grid refinement is at its easiest when discrete materials are set using
analytical - i.e. resolution independent - indicators.

3 files changed, 54 insertions, 33 deletions

diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index f30e9b4..c8a4659 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -35,8 +35,8 @@ typedef double T;
 
 #define DESCRIPTOR D2Q9Descriptor
 
-const T lx  = 4.;             // length of the channel
-const T ly  = 1.;             // height of the channel
+const T lx  = 4.0;            // length of the channel
+const T ly  = 1.0;            // height of the channel
 const int N = 30;             // resolution of the model
 const T Re  = 100.;           // Reynolds number
 const T maxPhysT = 60.;       // max. simulation time in s, SI unit
@@ -51,32 +51,45 @@ void prepareGeometry(UnitConverter<T,DESCRIPTOR> const& converter,
   OstreamManager clout(std::cout,"prepareGeometry");
   clout << "Prepare Geometry ..." << std::endl;
 
-  superGeometry.rename(0,2);
+  superGeometry.rename(0,1);
+
+  const T physSpacing = converter.getPhysDeltaX();
+
   // Set material number for bounce back boundaries
-  superGeometry.rename(2,1,0,1);
+  {
+    const Vector<T,2> wallExtend {lx+physSpacing,  physSpacing/2};
+    const Vector<T,2> wallOrigin {-physSpacing/4, -physSpacing/4};
+
+    IndicatorCuboid2D<T> lowerWall(wallExtend, wallOrigin);
+    superGeometry.rename(1,2,lowerWall);
 
-  T physSpacing = converter.getPhysDeltaX();
-  Vector<T,2> extend{  physSpacing / 2, ly };
-  Vector<T,2> origin{ -physSpacing / 4, 0  };
+    IndicatorCuboid2D<T> upperWall(wallExtend, wallOrigin + Vector<T,2> {0,ly});
+    superGeometry.rename(1,2,upperWall);
+  }
+
+  // Set material number for inflow and outflow
+  {
+    const Vector<T,2> extend { physSpacing/2,  ly};
+    const Vector<T,2> origin {-physSpacing/4, -physSpacing/4};
 
-  // Set material number for inflow
-  IndicatorCuboid2D<T> inflow(extend, origin);
-  superGeometry.rename(1,3,1,inflow);
+    IndicatorCuboid2D<T> inflow(extend, origin);
+    superGeometry.rename(1,3,inflow);
 
-  // Set material number for outflow
-  origin[0] = lx - physSpacing / 4;
-  IndicatorCuboid2D<T> outflow(extend, origin);
-  superGeometry.rename(1,4,1,outflow);
+    IndicatorCuboid2D<T> outflow(extend, origin + Vector<T,2> {lx,0});
+    superGeometry.rename(1,4,outflow);
+  }
 
-  IndicatorCuboid2D<T> obstacle(Vector<T,2> {0.2,0.1}, Vector<T,2> {1.25,0.45});
-  superGeometry.rename(1,5,obstacle);
+  // Set material number for vertically centered obstacle
+  {
+    const Vector<T,2> extend {0.1, 0.1};
+    const Vector<T,2> origin {1.25, (ly-extend[1])/2};
+    IndicatorCuboid2D<T> obstacle(extend, origin);
+    superGeometry.rename(1,2,obstacle);
+  }
 
-  // Removes all not needed boundary voxels outside the surface
   superGeometry.clean();
-  // Removes all not needed boundary voxels inside the surface
   superGeometry.innerClean();
   superGeometry.checkForErrors();
-
   superGeometry.print();
 
   clout << "Prepare Geometry ... OK" << std::endl;
@@ -203,27 +216,25 @@ int main(int argc, char* argv[])
   singleton::directories().setOutputDir("./tmp/");
   OstreamManager clout(std::cout,"main");
 
-  Vector<T,2> coarseOrigin { 0.0, 0.0 };
-  Vector<T,2> coarseExtend { lx, ly  };
+  const Vector<T,2> coarseOrigin {0.0, 0.0};
+  const Vector<T,2> coarseExtend {lx, ly};
   IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin);
 
-  Vector<T,2> fineOrigin   { 0.25 * lx, 0.2 };
-  Vector<T,2> fineExtend   { 0.6  * lx, 0.6 };
-
   auto coarseGrid = Grid2D<T,DESCRIPTOR>::make(coarseCuboid, N, 0.8, Re);
-  auto fineGrid   = &coarseGrid->refine(fineOrigin, fineExtend);
-
   prepareGeometry(coarseGrid->getConverter(), coarseGrid->getSuperGeometry());
+
+  const Vector<T,2> wantedFineExtend {2.0, 0.75};
+  const Vector<T,2> fineOrigin = coarseGrid->alignToGrid({0.8, (ly-wantedFineExtend[1])/2});
+  const Vector<T,2> fineExtend = coarseGrid->alignToGrid(fineOrigin + wantedFineExtend) - fineOrigin;
+
+  auto fineGrid = &coarseGrid->refine(fineOrigin, fineExtend);
   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);
-  coarseGrid->getSuperGeometry().rename(1,0,overlapCuboid);
-  coarseGrid->getSuperGeometry().rename(2,0,overlapCuboid);
-  coarseGrid->getSuperGeometry().rename(5,0,overlapCuboid);
+  IndicatorCuboid2D<T> refinedOverlap(fineExtend - 4*coarseDeltaX, fineOrigin + 2*coarseDeltaX);
+  coarseGrid->getSuperGeometry().rename(1,0,refinedOverlap);
+  coarseGrid->getSuperGeometry().rename(2,0,refinedOverlap);
 
   Dynamics<T, DESCRIPTOR>* coarseBulkDynamics;
   coarseBulkDynamics = new BGKdynamics<T, DESCRIPTOR>(
diff --git a/src/refinement/grid2D.h b/src/refinement/grid2D.h
index fda38ed..6d3579c 100644
--- a/src/refinement/grid2D.h
+++ b/src/refinement/grid2D.h
@@ -76,6 +76,8 @@ public:
   CoarseCoupler2D<T,DESCRIPTOR>& addCoarseCoupling(
     Grid2D<T,DESCRIPTOR>& fineGrid, Vector<T,2> origin, Vector<T,2> extend);
 
+  Vector<T,2> alignToGrid(Vector<T,2> physR) const;
+
   Grid2D<T,DESCRIPTOR>& refine(IndicatorF2D<T>& domainF);
   Grid2D<T,DESCRIPTOR>& refine(Vector<T,2> origin, Vector<T,2> extend);
 
diff --git a/src/refinement/grid2D.hh b/src/refinement/grid2D.hh
index 86c3102..d609746 100644
--- a/src/refinement/grid2D.hh
+++ b/src/refinement/grid2D.hh
@@ -169,13 +169,21 @@ Grid2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::refine(IndicatorF2D<T>& domainF)
     new Grid2D<T,DESCRIPTOR>(
       domainF,
       2*getConverter().getResolution(),
-      2.0*getConverter().getLatticeRelaxationTime() - 0.5,
+      2*getConverter().getLatticeRelaxationTime() - 0.5,
       getConverter().getReynoldsNumber()
     ));
   return *_fineGrids.back();
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
+Vector<T,2> Grid2D<T,DESCRIPTOR>::alignToGrid(Vector<T,2> physR) const
+{
+  Vector<int,3> latticeR{};
+  _cuboids->getLatticeR(physR, latticeR);
+  return _cuboids->getPhysR(latticeR.toStdVector());
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
 Grid2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::refine(Vector<T,2> origin, Vector<T,2> extend)
 {
   IndicatorCuboid2D<T> fineCuboid(extend, origin);