Continue cylinder2d refinement setup encapsulation

4 files changed, 52 insertions, 25 deletions

diff --git a/apps/adrian/cylinder2d/optimized_grid/cylinder2d.cpp b/apps/adrian/cylinder2d/optimized_grid/cylinder2d.cpp
index e5a934c..d516394 100644
--- a/apps/adrian/cylinder2d/optimized_grid/cylinder2d.cpp
+++ b/apps/adrian/cylinder2d/optimized_grid/cylinder2d.cpp
@@ -137,7 +137,6 @@ void prepareLattice(Grid2D<T,DESCRIPTOR>& grid)
           instances::getBulkMomenta<T,DESCRIPTOR>())));
 
   sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition = grid.getOnLatticeBoundaryCondition();
-  //createInterpBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryCondition);
   createLocalBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryCondition);
 
   const T omega = converter.getLatticeRelaxationFrequency();
@@ -251,26 +250,8 @@ void takeMeasurements(Grid2D<T,DESCRIPTOR>& grid, int iT, bool print=true)
   }
 }
 
-int main(int argc, char* argv[])
+void setupRefinement(Grid2D<T,DESCRIPTOR>& coarseGrid, Vector<T,2> domainOrigin, Vector<T,2> domainExtend)
 {
-  olbInit(&argc, &argv);
-  singleton::directories().setOutputDir("./tmp/");
-  OstreamManager clout(std::cout,"main");
-
-  const Vector<T,2> coarseOrigin {0.0, 0.0};
-  const Vector<T,2> coarseExtend {lx, ly};
-  IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin);
-
-  Grid2D<T,DESCRIPTOR> coarseGrid(
-    coarseCuboid,
-    RelaxationTime<T>(tau),
-    N,
-    PhysCharacteristics);
-  const Vector<T,2> domainOrigin = coarseGrid.getSuperGeometry().getStatistics().getMinPhysR(0);
-  const Vector<T,2> domainExtend = coarseGrid.getSuperGeometry().getStatistics().getPhysExtend(0);
-
-  prepareGeometry(coarseGrid, domainOrigin, domainExtend);
-
   const auto coarseDeltaX = coarseGrid.getConverter().getPhysDeltaX();
 
   const Vector<T,2> fineOutflowExtend {1*cylinderD, domainExtend[1]};
@@ -331,6 +312,28 @@ int main(int argc, char* argv[])
   auto& fineGrid3 = fineGrid2.refine(fineOrigin3, fineExtend3);
   prepareGeometry(fineGrid3, domainOrigin, domainExtend);
   disableRefinedArea(fineGrid2, fineGrid3);
+}
+
+int main(int argc, char* argv[])
+{
+  olbInit(&argc, &argv);
+  singleton::directories().setOutputDir("./tmp/");
+  OstreamManager clout(std::cout,"main");
+
+  const Vector<T,2> coarseOrigin {0.0, 0.0};
+  const Vector<T,2> coarseExtend {lx, ly};
+  IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin);
+
+  Grid2D<T,DESCRIPTOR> coarseGrid(
+    coarseCuboid,
+    RelaxationTime<T>(tau),
+    N,
+    PhysCharacteristics);
+  const Vector<T,2> domainOrigin = coarseGrid.getSuperGeometry().getStatistics().getMinPhysR(0);
+  const Vector<T,2> domainExtend = coarseGrid.getSuperGeometry().getStatistics().getPhysExtend(0);
+
+  prepareGeometry(coarseGrid, domainOrigin, domainExtend);
+  setupRefinement(coarseGrid, domainOrigin, domainExtend);
 
   coarseGrid.forEachGrid(prepareLattice);
 
@@ -343,6 +346,8 @@ int main(int argc, char* argv[])
     coarseGrid.getSuperGeometry().getStatistics().getNvoxel());
   timer.start();
 
+  Grid2D<T,DESCRIPTOR>& cylinderGrid = coarseGrid.locate({cylinderX, cylinderY});
+
   for (int iT = 0; iT <= coarseGrid.getConverter().getLatticeTime(maxPhysT); ++iT) {
     setBoundaryValues(coarseGrid, iT);
 
@@ -356,7 +361,7 @@ int main(int argc, char* argv[])
         getResults(grid, id, iT);
       });
 
-      takeMeasurements(fineGrid3, iT);
+      takeMeasurements(cylinderGrid, iT);
     }
   }
 
diff --git a/src/refinement/coupler2D.hh b/src/refinement/coupler2D.hh
index 92c8d26..373a6de 100644
--- a/src/refinement/coupler2D.hh
+++ b/src/refinement/coupler2D.hh
@@ -73,7 +73,8 @@ Coupler2D<T,DESCRIPTOR>::Coupler2D(Grid2D<T,DESCRIPTOR>& coarse, Grid2D<T,DESCRI
   const auto& fineGeometry   = _fine.getCuboidGeometry();
 
   const T deltaX = _fine.getConverter().getPhysDeltaX();
-  const Vector<T,2> stepPhysR = _vertical ? Vector<T,2> {0, deltaX} : Vector<T,2> {deltaX, 0};
+  const Vector<T,2> stepPhysR = _vertical ? Vector<T,2> {0, deltaX} :
+                                Vector<T,2> {deltaX, 0};
 
   for (int i=0; i < _fineSize; ++i) {
     if (i % 2 == 0) {
@@ -212,7 +213,8 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
     fineLattice.get(finePos, fineCell);
 
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr) + this->getScalingFactor() * fneq[iPop];
+      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr)
+                       + this->getScalingFactor() * fneq[iPop];
     }
 
     fineLattice.set(finePos, fineCell);
@@ -230,7 +232,8 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
     fineLattice.get(finePos, fineCell);
 
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr) + this->getScalingFactor() * fneq[iPop];
+      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr)
+                       + this->getScalingFactor() * fneq[iPop];
     }
 
     fineLattice.set(finePos, fineCell);
@@ -248,7 +251,8 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
     fineLattice.get(finePos, fineCell);
 
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr) + this->getScalingFactor() * fneq[iPop];
+      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr)
+                       + this->getScalingFactor() * fneq[iPop];
     }
 
     fineLattice.set(finePos, fineCell);
diff --git a/src/refinement/grid2D.h b/src/refinement/grid2D.h
index 6cbbf6f..f538850 100644
--- a/src/refinement/grid2D.h
+++ b/src/refinement/grid2D.h
@@ -147,6 +147,12 @@ public:
   void forEachGrid(std::function<void(Grid2D<T,DESCRIPTOR>&)>&& f);
   void forEachGrid(const std::string& id, std::function<void(Grid2D<T,DESCRIPTOR>&,const std::string&)>&& f);
 
+  /// Returns the finest grid representing a physical position
+  /**
+   * Only works if pos is actually contained in a node of the refinement tree.
+   **/
+  Grid2D<T,DESCRIPTOR>& locate(Vector<T,2> pos);
+
   std::size_t getActiveVoxelN() const;
 
 };
diff --git a/src/refinement/grid2D.hh b/src/refinement/grid2D.hh
index fdf6461..9ca8aaf 100644
--- a/src/refinement/grid2D.hh
+++ b/src/refinement/grid2D.hh
@@ -281,6 +281,18 @@ void Grid2D<T,DESCRIPTOR>::forEachGrid(
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
+Grid2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::locate(Vector<T,2> pos)
+{
+  int iC;
+  for (auto& grid : _fineGrids) {
+    if (grid->getCuboidGeometry().getC(pos, iC)) {
+      return grid->locate(pos);
+    }
+  }
+  return *this;
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
 std::size_t Grid2D<T,DESCRIPTOR>::getActiveVoxelN() const
 {
   std::size_t n = _geometry->getStatistics().getNvoxel();