Improve grid refinement interface

By deriving RefiningGrid2D from Grid2D to provide methods specific to refining subgrids.

3 files changed, 98 insertions, 39 deletions

diff --git a/src/refinement/grid2D.cpp b/src/refinement/grid2D.cpp
index b371d0c..86ca284 100644
--- a/src/refinement/grid2D.cpp
+++ b/src/refinement/grid2D.cpp
@@ -30,5 +30,6 @@
 namespace olb {
 
 template class Grid2D<double,descriptors::D2Q9Descriptor>;
+template class RefiningGrid2D<double,descriptors::D2Q9Descriptor>;
 
 }
diff --git a/src/refinement/grid2D.h b/src/refinement/grid2D.h
index 1d15170..d782651 100644
--- a/src/refinement/grid2D.h
+++ b/src/refinement/grid2D.h
@@ -34,23 +34,27 @@
 #include "geometry/superGeometry2D.h"
 #include "communication/loadBalancer.h"
 #include "functors/analytical/indicator/indicatorF2D.h"
+#include "utilities/functorPtr.h"
 
 namespace olb {
 
 
 template <typename T, template<typename> class DESCRIPTOR> class FineCoupler2D;
 template <typename T, template<typename> class DESCRIPTOR> class CoarseCoupler2D;
+template <typename T, template<typename> class DESCRIPTOR> class RefiningGrid2D;
 
 template <typename T, template<typename> class DESCRIPTOR>
 class Grid2D {
-private:
+protected:
+  FunctorPtr<IndicatorF2D<T>> _domainF;
+
   std::unique_ptr<UnitConverter<T,DESCRIPTOR>>  _converter;
   std::unique_ptr<CuboidGeometry2D<T>>          _cuboids;
   std::unique_ptr<LoadBalancer<T>>              _balancer;
   std::unique_ptr<SuperGeometry2D<T>>           _geometry;
   std::unique_ptr<SuperLattice2D<T,DESCRIPTOR>> _lattice;
 
-  std::vector<std::unique_ptr<Grid2D<T,DESCRIPTOR>>> _fineGrids;
+  std::vector<std::unique_ptr<RefiningGrid2D<T,DESCRIPTOR>>> _fineGrids;
 
   std::vector<std::unique_ptr<FineCoupler2D<T,DESCRIPTOR>>>   _fineCouplers;
   std::vector<std::unique_ptr<CoarseCoupler2D<T,DESCRIPTOR>>> _coarseCouplers;
@@ -58,7 +62,7 @@ private:
 public:
   static std::unique_ptr<Grid2D<T,DESCRIPTOR>> make(IndicatorF2D<T>& domainF, int resolution, T tau, int re);
 
-  Grid2D(IndicatorF2D<T>& domainF, int resolution, T tau, int re);
+  Grid2D(FunctorPtr<IndicatorF2D<T>>&& domainF, int resolution, T tau, int re);
 
   UnitConverter<T,DESCRIPTOR>& getConverter();
   CuboidGeometry2D<T>& getCuboidGeometry();
@@ -76,8 +80,26 @@ public:
   Vector<T,2> alignOriginToGrid(Vector<T,2> physR) const;
   Vector<T,2> alignExtendToGrid(Vector<T,2> physR) const;
 
-  Grid2D<T,DESCRIPTOR>& refine(IndicatorF2D<T>& domainF);
-  Grid2D<T,DESCRIPTOR>& refine(Vector<T,2> origin, Vector<T,2> extend);
+  RefiningGrid2D<T,DESCRIPTOR>& refine(Vector<T,2> origin, Vector<T,2> extend, bool addCouplers=true);
+
+};
+
+template <typename T, template<typename> class DESCRIPTOR>
+class RefiningGrid2D : public Grid2D<T,DESCRIPTOR> {
+private:
+  const Vector<T,2> _origin;
+  const Vector<T,2> _extend;
+
+  Grid2D<T,DESCRIPTOR>& _parentGrid;
+
+public:
+  RefiningGrid2D(Grid2D<T,DESCRIPTOR>& parentGrid, Vector<T,2> origin, Vector<T,2> extend);
+
+  Vector<T,2> getOrigin() const;
+  Vector<T,2> getExtend() const;
+
+  /// Indicates the subdomain of the coarse grid rendered moot by refinement
+  std::unique_ptr<IndicatorF2D<T>> getRefinedOverlap() const;
 
 };
 
diff --git a/src/refinement/grid2D.hh b/src/refinement/grid2D.hh
index 99a2563..aa82240 100644
--- a/src/refinement/grid2D.hh
+++ b/src/refinement/grid2D.hh
@@ -27,6 +27,7 @@
 #include "grid2D.h"
 #include "coupler2D.hh"
 
+#include "utilities/vectorHelpers.h"
 #include "communication/heuristicLoadBalancer.h"
 
 namespace olb {
@@ -43,7 +44,8 @@ std::unique_ptr<Grid2D<T,DESCRIPTOR>> Grid2D<T,DESCRIPTOR>::make(
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
-Grid2D<T,DESCRIPTOR>::Grid2D(IndicatorF2D<T>& domainF, int resolution, T tau, int re):
+Grid2D<T,DESCRIPTOR>::Grid2D(FunctorPtr<IndicatorF2D<T>>&& domainF, int resolution, T tau, int re):
+  _domainF(std::move(domainF)),
   _converter(new UnitConverterFromResolutionAndRelaxationTime<T,DESCRIPTOR>(
                resolution, // resolution: number of voxels per charPhysL
                tau,        // latticeRelaxationTime: relaxation time, has to be greater than 0.5!
@@ -53,7 +55,7 @@ Grid2D<T,DESCRIPTOR>::Grid2D(IndicatorF2D<T>& domainF, int resolution, T tau, in
                T{1},       // physDensity: physical density in __kg / m^3__
                T{1})),
   _cuboids(new CuboidGeometry2D<T>(
-             domainF,
+             *_domainF,
              _converter->getConversionFactorLength(),
 #ifdef PARALLEL_MODE_MPI
              singleton::mpi().getSize()
@@ -156,19 +158,6 @@ CoarseCoupler2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::addCoarseCoupling(
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
-Grid2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::refine(IndicatorF2D<T>& domainF)
-{
-  _fineGrids.emplace_back(
-    new Grid2D<T,DESCRIPTOR>(
-      domainF,
-      2*getConverter().getResolution(),
-      2*getConverter().getLatticeRelaxationTime() - 0.5,
-      getConverter().getReynoldsNumber()
-    ));
-  return *_fineGrids.back();
-}
-
-template <typename T, template<typename> class DESCRIPTOR>
 Vector<T,2> Grid2D<T,DESCRIPTOR>::alignOriginToGrid(Vector<T,2> physR) const
 {
   Vector<int,3> latticeR{};
@@ -180,37 +169,84 @@ template <typename T, template<typename> class DESCRIPTOR>
 Vector<T,2> Grid2D<T,DESCRIPTOR>::alignExtendToGrid(Vector<T,2> extend) const
 {
   const T deltaX = _converter->getPhysDeltaX();
-  return floor(extend / deltaX) * deltaX;
+  return util::floor(extend / deltaX) * deltaX;
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
-Grid2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::refine(Vector<T,2> origin, Vector<T,2> extend)
+RefiningGrid2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::refine(
+  Vector<T,2> wantedOrigin, Vector<T,2> wantedExtend, bool addCouplers)
 {
-  IndicatorCuboid2D<T> fineCuboid(extend, origin);
-  auto& fineGrid = refine(fineCuboid);
+  if (addCouplers) {
+    auto& fineGrid = refine(wantedOrigin, wantedExtend, false);
+
+    const Vector<T,2> origin = fineGrid.getOrigin();
+    const Vector<T,2> extend = fineGrid.getExtend();
+
+    const Vector<T,2> extendX = {extend[0],0};
+    const Vector<T,2> extendY = {0,extend[1]};
 
-  const Vector<T,2> extendX = {extend[0],0};
-  const Vector<T,2> extendY = {0,extend[1]};
+    addFineCoupling(fineGrid, origin,           extendY);
+    addFineCoupling(fineGrid, origin + extendX, extendY);
+    addFineCoupling(fineGrid, origin + extendY, extendX);
+    addFineCoupling(fineGrid, origin,           extendX);
 
-  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};
 
-  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);
 
-  addCoarseCoupling(fineGrid, innerOrigin,                innerExtendY);
-  addCoarseCoupling(fineGrid, innerOrigin + innerExtendX, innerExtendY);
-  addCoarseCoupling(fineGrid, innerOrigin + innerExtendY, innerExtendX);
-  addCoarseCoupling(fineGrid, innerOrigin,                innerExtendX);
+    return fineGrid;
+  }
+  else {
+    const Vector<T,2> origin = alignOriginToGrid(wantedOrigin);
+    const Vector<T,2> extend = alignExtendToGrid(wantedExtend);
 
-  return fineGrid;
+    _fineGrids.emplace_back(
+      new RefiningGrid2D<T,DESCRIPTOR>(*this, origin, extend));
+    return *_fineGrids.back();
+  }
 }
 
 
+template <typename T, template<typename> class DESCRIPTOR>
+RefiningGrid2D<T,DESCRIPTOR>::RefiningGrid2D(
+  Grid2D<T,DESCRIPTOR>& parentGrid, Vector<T,2> origin, Vector<T,2> extend):
+  Grid2D<T,DESCRIPTOR>(
+    std::unique_ptr<IndicatorF2D<T>>(new IndicatorCuboid2D<T>(extend, origin)),
+    2*parentGrid.getConverter().getResolution(),
+    2*parentGrid.getConverter().getLatticeRelaxationTime() - 0.5,
+    parentGrid.getConverter().getReynoldsNumber()),
+  _origin(origin),
+  _extend(extend),
+  _parentGrid(parentGrid)
+{ }
+
+template <typename T, template<typename> class DESCRIPTOR>
+Vector<T,2> RefiningGrid2D<T,DESCRIPTOR>::getOrigin() const
+{
+  return _origin;
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+Vector<T,2> RefiningGrid2D<T,DESCRIPTOR>::getExtend() const
+{
+  return _extend;
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+std::unique_ptr<IndicatorF2D<T>> RefiningGrid2D<T,DESCRIPTOR>::getRefinedOverlap() const
+{
+  const T coarseDeltaX = _parentGrid.getConverter().getPhysDeltaX();
+
+  return std::unique_ptr<IndicatorF2D<T>>(
+    new IndicatorCuboid2D<T>(_extend - 4*coarseDeltaX, _origin + 2*coarseDeltaX));
+}
+
 }
 
 #endif