Initialize at openlb-1-3

1 files changed, 677 insertions, 0 deletions

diff --git a/src/geometry/cuboidGeometry2D.hh b/src/geometry/cuboidGeometry2D.hh
new file mode 100644
index 0000000..f7a68e9
--- /dev/null
+++ b/src/geometry/cuboidGeometry2D.hh
@@ -0,0 +1,677 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2007, 2014 Mathias J. Krause
+ *  E-mail contact: info@openlb.net
+ *  The most recent release of OpenLB can be downloaded at
+ *  <http://www.openlb.net/>
+ *
+ *  This program is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU General Public License
+ *  as published by the Free Software Foundation; either version 2
+ *  of the License, or (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public
+ *  License along with this program; if not, write to the Free
+ *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ *  Boston, MA  02110-1301, USA.
+*/
+
+/** \file
+ * The description of a vector of 2D cuboid  -- generic implementation.
+ */
+
+
+#ifndef CUBOID_GEOMETRY_2D_HH
+#define CUBOID_GEOMETRY_2D_HH
+
+
+#include <iostream>
+#include <math.h>
+#include "geometry/cuboidGeometry2D.h"
+#include "functors/analytical/indicator/indicatorF2D.h"
+
+namespace olb {
+
+
+////////////////////// Class CuboidGeometry2D /////////////////////////
+
+template<typename T>
+CuboidGeometry2D<T>::CuboidGeometry2D()
+  : _motherCuboid(0,0,0,0,0), _periodicityOn(3, bool(false)), clout(std::cout, "CuboidGeometry3D")
+{
+  add(_motherCuboid);
+  split(0, 1);
+}
+
+template<typename T>
+CuboidGeometry2D<T>::CuboidGeometry2D(T originX, T originY, T deltaR, int nX, int nY, int nC)
+  : _motherCuboid(originX, originY, deltaR, nX, nY), _periodicityOn(2, bool(false)),
+    clout(std::cout, "CuboidGeometry2D")
+{
+  add(_motherCuboid);
+  split(0, nC);
+}
+
+template<typename T>
+CuboidGeometry2D<T>::CuboidGeometry2D(IndicatorF2D<T>& indicatorF, T voxelSize, int nC)
+  : _motherCuboid(indicatorF.getMin()[0],  indicatorF.getMin()[1], voxelSize,
+                  (int)((indicatorF.getMax()[0] - indicatorF.getMin()[0]) / voxelSize + 1.5),
+                  (int)((indicatorF.getMax()[1] - indicatorF.getMin()[1]) / voxelSize + 1.5)),
+    _periodicityOn(2, bool(false)), clout(std::cout, "CuboidGeometry2D")
+{
+
+  add(_motherCuboid);
+  split(0, nC);
+  shrink(indicatorF);
+}
+
+
+template<typename T>
+void CuboidGeometry2D<T>::reInit(T globPosX, T globPosY, T delta, int nX, int nY, int nC)
+{
+  _cuboids.clear();
+  _motherCuboid = Cuboid2D<T>(globPosX, globPosY, delta, nX, nY);
+  Cuboid2D<T> cuboid(globPosX, globPosY, delta, nX, nY);
+  if (_oldApproach) {
+    cuboid.init(0, 0, 1, nX, nY);
+  }
+
+  add(cuboid);
+  split(0, nC);
+}
+
+template<typename T>
+Cuboid2D<T>& CuboidGeometry2D<T>::get(int i)
+{
+  return _cuboids[i];
+}
+
+template<typename T>
+Cuboid2D<T> const& CuboidGeometry2D<T>::get(int i) const
+{
+  return _cuboids[i];
+}
+
+template<typename T>
+void CuboidGeometry2D<T>::setPeriodicity(bool periodicityX, bool periodicityY)
+{
+  _periodicityOn.resize(2);
+  _periodicityOn[0] = periodicityX;
+  _periodicityOn[1] = periodicityY;
+}
+
+
+template<typename T>
+bool CuboidGeometry2D<T>::getC(std::vector<T> physR, int& iC) const
+{
+  int iCtmp = get_iC(physR[0], physR[1]);
+  if (iCtmp < getNc()) {
+    iC = iCtmp;
+    return true;
+  }
+  else {
+    return false;
+  }
+}
+
+template<typename T>
+bool CuboidGeometry2D<T>::getC(const Vector<T,2>& physR, int& iC) const
+{
+  int iCtmp = get_iC(physR[0], physR[1]);
+  if (iCtmp < getNc()) {
+    iC = iCtmp;
+    return true;
+  }
+  else {
+    return false;
+  }
+}
+
+template<typename T>
+bool CuboidGeometry2D<T>::getLatticeR(std::vector<T> physR, std::vector<int>& latticeR) const
+{
+  return getLatticeR(&latticeR[0], &physR[0]);
+  /*  int iCtmp = get_iC(physR[0], physR[1]);
+    if (iCtmp < getNc()) {
+      latticeR[0] = iCtmp;
+      latticeR[1] = (int)floor( (physR[0] - _cuboids[latticeR[0]].getOrigin()[0] ) / _cuboids[latticeR[0]].getDeltaR() + .5);
+      latticeR[2] = (int)floor( (physR[1] - _cuboids[latticeR[0]].getOrigin()[1] ) / _cuboids[latticeR[0]].getDeltaR() + .5);
+      return true;
+    } else {
+      return false;
+    }*/
+}
+
+template<typename T>
+bool CuboidGeometry2D<T>::getLatticeR(int latticeR[], const T physR[]) const
+{
+  int iCtmp = get_iC(physR[0], physR[1]);
+  if (iCtmp < getNc()) {
+    latticeR[0] = iCtmp;
+    latticeR[1] = (int)floor( (physR[0] - _cuboids[latticeR[0]].getOrigin()[0] ) / _cuboids[latticeR[0]].getDeltaR() + .5);
+    latticeR[2] = (int)floor( (physR[1] - _cuboids[latticeR[0]].getOrigin()[1] ) / _cuboids[latticeR[0]].getDeltaR() + .5);
+    return true;
+  }
+  else {
+    return false;
+  }
+}
+
+template<typename T>
+bool CuboidGeometry2D<T>::getLatticeR(
+  const Vector<T,2>& physR, Vector<int,3>& latticeR) const
+{
+  int iCtmp = get_iC(physR[0], physR[1]);
+  if (iCtmp < getNc()) {
+    latticeR[0] = iCtmp;
+    latticeR[1] = (int)floor( (physR[0] - _cuboids[latticeR[0]].getOrigin()[0] ) / _cuboids[latticeR[0]].getDeltaR() + .5);
+    latticeR[2] = (int)floor( (physR[1] - _cuboids[latticeR[0]].getOrigin()[1] ) / _cuboids[latticeR[0]].getDeltaR() + .5);
+    return true;
+  }
+  else {
+    return false;
+  }
+}
+
+template<typename T>
+bool CuboidGeometry2D<T>::getFloorLatticeR(std::vector<T> physR, std::vector<int>& latticeR) const
+{
+  int iCtmp = get_iC(physR[0], physR[1]);
+  if (iCtmp < getNc()) {
+    latticeR[0] = iCtmp;
+    latticeR[1] = (int)floor( (physR[0] - _cuboids[latticeR[0]].getOrigin()[0] ) / _cuboids[latticeR[0]].getDeltaR() );
+    latticeR[2] = (int)floor( (physR[1] - _cuboids[latticeR[0]].getOrigin()[1] ) / _cuboids[latticeR[0]].getDeltaR() );
+    return true;
+  }
+  else {
+    return false;
+  }
+}
+
+template<typename T>
+bool CuboidGeometry2D<T>::getFloorLatticeR(
+  const Vector<T,2>& physR, Vector<int,3>& latticeR) const
+{
+  int iCtmp = get_iC(physR[0], physR[1]);
+  if (iCtmp < getNc()) {
+    latticeR[0] = iCtmp;
+    latticeR[1] = (int)floor( (physR[0] - _cuboids[latticeR[0]].getOrigin()[0] ) / _cuboids[latticeR[0]].getDeltaR() );
+    latticeR[2] = (int)floor( (physR[1] - _cuboids[latticeR[0]].getOrigin()[1] ) / _cuboids[latticeR[0]].getDeltaR() );
+    return true;
+  }
+  else {
+    return false;
+  }
+}
+
+template<typename T>
+std::vector<T> CuboidGeometry2D<T>::getPhysR(int iCglob, int iX, int iY) const
+{
+  std::vector<T> physR(2,T());
+  _cuboids[iCglob].getPhysR(&(physR[0]), iX, iY);
+  for (int iDim = 0; iDim < 2; iDim++) {
+    if (_periodicityOn[iDim]) {
+      //std::cout << iDim << _periodicityOn[iDim] <<":"<< _motherCuboid.getDeltaR()*(_motherCuboid.getExtend()[iDim]) << std::endl;
+      physR[iDim] = remainder( physR[iDim] - _motherCuboid.getOrigin()[iDim]
+                               + _motherCuboid.getDeltaR() * (_motherCuboid.getExtend()[iDim]),
+                               _motherCuboid.getDeltaR() * (_motherCuboid.getExtend()[iDim]));
+      // solving the rounding error problem for double
+      if ( physR[iDim]*physR[iDim] < 0.001 * _motherCuboid.getDeltaR()*_motherCuboid.getDeltaR() ) {
+        physR[iDim] = T();
+      }
+      // make it to mod instead remainer
+      if ( physR[iDim] < 0 ) {
+        physR[iDim] += _motherCuboid.getDeltaR() * _motherCuboid.getExtend()[iDim];
+      }
+      // add origin
+      physR[iDim] += _motherCuboid.getOrigin()[iDim];
+    }
+  }
+  return physR;
+}
+
+template<typename T>
+std::vector<T> CuboidGeometry2D<T>::getPhysR(std::vector<int> latticeR) const
+{
+  return getPhysR(latticeR[0], latticeR[1], latticeR[2]);
+}
+
+template<typename T>
+void CuboidGeometry2D<T>::getPhysR(T output[2], const int latticeR[3]) const
+{
+  getPhysR(output, latticeR[0], latticeR[1], latticeR[2]);
+}
+
+template<typename T>
+void CuboidGeometry2D<T>::getPhysR(T output[2], const int iCglob, const int iX, const int iY) const
+{
+  _cuboids[iCglob].getPhysR(output, iX, iY);
+  for (int iDim = 0; iDim < 2; iDim++) {
+    if (_periodicityOn[iDim]) {
+      //std::cout << iDim << _periodicityOn[iDim] <<":"<< _motherCuboid.getDeltaR()*(_motherCuboid.getExtend()[iDim]) << std::endl;
+      output[iDim] = remainder( output[iDim] - _motherCuboid.getOrigin()[iDim]
+                                + _motherCuboid.getDeltaR() * (_motherCuboid.getExtend()[iDim]),
+                                _motherCuboid.getDeltaR() * (_motherCuboid.getExtend()[iDim]));
+      // solving the rounding error problem for double
+      if ( output[iDim]*output[iDim] < 0.001 * _motherCuboid.getDeltaR()*_motherCuboid.getDeltaR() ) {
+        if ( output[iDim] > 0 ) {
+          output[iDim] = _motherCuboid.getDeltaR() * _motherCuboid.getExtend()[iDim];
+        }
+        else {
+          output[iDim] = T();
+        }
+      }
+      // make it to mod instead remainer
+      if ( output[iDim] < 0 ) {
+        output[iDim] += _motherCuboid.getDeltaR() * _motherCuboid.getExtend()[iDim];
+      }
+      // add origin
+      output[iDim] += _motherCuboid.getOrigin()[iDim];
+    }
+  }
+}
+
+template<typename T>
+int CuboidGeometry2D<T>::getNc() const
+{
+  return _cuboids.size();
+}
+
+template<typename T>
+T CuboidGeometry2D<T>::getMinRatio() const
+{
+  T minRatio = 1.;
+  for ( auto& cuboid : _cuboids ) {
+    if ((T)cuboid.getNx() / (T)cuboid.getNy() < minRatio) {
+      minRatio = (T)cuboid.getNx() / (T)cuboid.getNy();
+    }
+  }
+  return minRatio;
+}
+
+template<typename T>
+T CuboidGeometry2D<T>::getMaxRatio() const
+{
+  T maxRatio = 1.;
+  for (unsigned i = 0; i < _cuboids.size(); i++) {
+    if ((T)_cuboids[i].getNx() / (T)_cuboids[i].getNy() > maxRatio) {
+      maxRatio = (T)_cuboids[i].getNx() / (T)_cuboids[i].getNy();
+    }
+  }
+  return maxRatio;
+}
+
+template<typename T>
+std::vector<T> CuboidGeometry2D<T>::getMinPhysR() const
+{
+  Vector<T,2> output(_cuboids[0].getOrigin());
+  for (unsigned i = 0; i < _cuboids.size(); i++) {
+    if (_cuboids[i].getOrigin()[0] < output[0]) {
+      output[0] = _cuboids[i].getOrigin()[0];
+    }
+    if (_cuboids[i].getOrigin()[1] < output[1]) {
+      output[1] = _cuboids[i].getOrigin()[1];
+    }
+  }
+  return std::vector<T> {output[0],output[1]};
+}
+
+template<typename T>
+std::vector<T> CuboidGeometry2D<T>::getMaxPhysR() const
+{
+  Vector<T,2> output(_cuboids[0].getOrigin());
+  output[0] += _cuboids[0].getNx()*_cuboids[0].getDeltaR();
+  output[1] += _cuboids[0].getNy()*_cuboids[0].getDeltaR();
+  for (unsigned i = 0; i < _cuboids.size(); i++) {
+    if (_cuboids[i].getOrigin()[0] + _cuboids[i].getNx()*_cuboids[i].getDeltaR() > output[0]) {
+      output[0] = _cuboids[i].getOrigin()[0] + _cuboids[i].getNx()*_cuboids[i].getDeltaR();
+    }
+    if (_cuboids[i].getOrigin()[1] + _cuboids[i].getNy()*_cuboids[i].getDeltaR() > output[1]) {
+      output[1] = _cuboids[i].getOrigin()[1] + _cuboids[i].getNy()*_cuboids[i].getDeltaR();
+    }
+  }
+  return std::vector<T> {output[0],output[1]};
+}
+
+template<typename T>
+T CuboidGeometry2D<T>::getMinPhysVolume() const
+{
+  T minVolume = _cuboids[0].getPhysVolume();
+  for (unsigned i = 0; i < _cuboids.size(); i++) {
+    if (_cuboids[i].getPhysVolume() < minVolume) {
+      minVolume = _cuboids[i].getPhysVolume();
+    }
+  }
+  return minVolume;
+}
+
+template<typename T>
+T CuboidGeometry2D<T>::getMaxPhysVolume() const
+{
+  T maxVolume = _cuboids[0].getPhysVolume();
+  for (unsigned i = 0; i < _cuboids.size(); i++) {
+    if (_cuboids[i].getPhysVolume() > maxVolume) {
+      maxVolume = _cuboids[i].getPhysVolume();
+    }
+  }
+  return maxVolume;
+}
+
+template<typename T>
+size_t CuboidGeometry2D<T>::getMinLatticeVolume() const
+{
+  size_t minNodes = _cuboids[0].getLatticeVolume();
+  for (unsigned i = 0; i < _cuboids.size(); i++) {
+    if (_cuboids[i].getLatticeVolume() < minNodes) {
+      minNodes = _cuboids[i].getLatticeVolume();
+    }
+  }
+  return minNodes;
+}
+
+template<typename T>
+size_t CuboidGeometry2D<T>::getMaxLatticeVolume() const
+{
+  size_t maxNodes = _cuboids[0].getLatticeVolume();
+  for (unsigned i = 0; i < _cuboids.size(); i++) {
+    if (_cuboids[i].getLatticeVolume() > maxNodes) {
+      maxNodes = _cuboids[i].getLatticeVolume();
+    }
+  }
+  return maxNodes;
+}
+
+template<typename T>
+T CuboidGeometry2D<T>::getMinDeltaR() const
+{
+  T minDelta = _cuboids[0].getDeltaR();
+  for (unsigned i = 0; i < _cuboids.size(); i++) {
+    if (_cuboids[i].getDeltaR() < minDelta) {
+      minDelta = _cuboids[i].getDeltaR();
+    }
+  }
+  return minDelta;
+}
+
+template<typename T>
+T CuboidGeometry2D<T>::getMaxDeltaR() const
+{
+  T maxDelta = _cuboids[0].getDeltaR();
+  for (unsigned i = 0; i < _cuboids.size(); i++) {
+    if (_cuboids[i].getDeltaR() > maxDelta) {
+      maxDelta = _cuboids[i].getDeltaR();
+    }
+  }
+  return maxDelta;
+}
+
+template<typename T>
+Cuboid2D<T> CuboidGeometry2D<T>::getMotherCuboid() const
+{
+
+  /*Cuboid2D<T> found;
+  if(_cuboids.size()==0) {
+      found.init(0, 0, 0, 0, 0);
+      return found;
+  }
+
+  T delta = _cuboids[0].getDeltaR();
+  T globPosXmin = _cuboids[0].get_globPosX();
+  T globPosYmin = _cuboids[0].get_globPosY();
+  T globPosXmax = _cuboids[0].get_globPosX() + delta*(_cuboids[0].getNx()-1);
+  T globPosYmax = _cuboids[0].get_globPosY() + delta*(_cuboids[0].getNy()-1);
+
+  for (unsigned i=1; i<_cuboids.size(); i++) {
+      if(delta > _cuboids[i].getDeltaR() ) {
+          delta = _cuboids[i].getDeltaR();
+      }
+      if(globPosXmin > _cuboids[i].get_globPosX() ) {
+          globPosXmin = _cuboids[i].get_globPosX();
+      }
+      if(globPosYmin > _cuboids[i].get_globPosY() ) {
+          globPosYmin = _cuboids[i].get_globPosY();
+      }
+      if(globPosXmax < _cuboids[i].get_globPosX()
+                          + delta*(_cuboids[i].getNx()-1)) {
+          globPosXmax = _cuboids[i].get_globPosX()
+                          + delta*(_cuboids[i].getNx()-1);
+      }
+      if(globPosYmax < _cuboids[i].get_globPosY()
+                          + delta*(_cuboids[i].getNy()-1)) {
+          globPosYmax = _cuboids[i].get_globPosY()
+                          + delta*(_cuboids[i].getNy()-1);
+      }
+  }
+  int nX = int(ceil((globPosXmax - globPosXmin)/delta))+1;
+  int nY = int(ceil((globPosYmax - globPosYmin)/delta))+1;
+
+  found.init(globPosXmin, globPosYmin, delta, nX, nY);
+
+  return found;*/
+  return _motherCuboid;
+}
+
+template<typename T>
+int CuboidGeometry2D<T>::get_iC(T globX, T globY, int offset) const
+{
+  unsigned i;
+  for (i = 0; i < _cuboids.size(); i++) {
+    if (_cuboids[i].checkPoint(globX, globY, offset)) {
+      return (int)i;
+    }
+  }
+  return (int)i;
+}
+
+template<typename T>
+int CuboidGeometry2D<T>::get_iC(T globX, T globY, int orientationX, int orientationY) const
+{
+  unsigned i;
+  for (i = 0; i < _cuboids.size(); i++) {
+    if (_cuboids[i].checkPoint(globX, globY) &&
+        _cuboids[i].checkPoint(globX + orientationX / _cuboids[i].getDeltaR(),
+                               globY + orientationY / _cuboids[i].getDeltaR())) {
+      return (int)i;
+    }
+  }
+  return (int)i;
+}
+
+
+template<typename T>
+void CuboidGeometry2D<T>::add(Cuboid2D<T> cuboid)
+{
+
+  _cuboids.push_back(cuboid);
+}
+
+template<typename T>
+void CuboidGeometry2D<T>::remove(int iC)
+{
+
+  _cuboids.erase(_cuboids.begin() + iC);
+}
+/*
+template<typename T>
+void CuboidGeometry2D<T>::remove(olb::ScalarField2D<int>* geometryData)
+{
+
+  std::vector<Cuboid2D<T> > cuboids;
+  unsigned size = _cuboids.size();
+
+  std::vector<bool> allZero;
+  for (unsigned i = 0; i < size; i++) {
+    allZero.push_back(1);
+    for (int iX = 0; iX < _cuboids[i].getNx(); iX++) {
+      for (int iY = 0; iY < _cuboids[i].getNy(); iY++) {
+        if (geometryData->get(_cuboids[i].get_globPosX() + iX,
+                              _cuboids[i].get_globPosY() + iY) != 0 ) {
+          allZero[i] = 0;
+        }
+      }
+    }
+  }
+  for (unsigned i = 0; i < size; i++) {
+    if (!allZero[i] ) {
+      cuboids.push_back(_cuboids[i]);
+    }
+  }
+  _cuboids.clear();
+  for (unsigned i = 0; i < cuboids.size(); i++) {
+    _cuboids.push_back(cuboids[i]);
+  }
+}*/
+
+template<typename T>
+void CuboidGeometry2D<T>::shrink(IndicatorF2D<T>& indicatorF)
+{
+  //IndicatorIdentity3D<T> tmpIndicatorF(indicatorF);
+  int newX, newY, maxX, maxY;
+  int nC = getNc();
+  std::vector<int> latticeR(3, 0);
+  std::vector<T> physR(2, T());
+  bool inside[1];
+  for (int iC = nC - 1; iC >= 0; iC--) {
+    latticeR[0] = iC;
+    int fullCells = 0;
+    int xN = get(iC).getNx();
+    int yN = get(iC).getNy();
+    maxX = 0;
+    maxY = 0;
+    newX = xN - 1;
+    newY = yN - 1;
+    for (int iX = 0; iX < xN; iX++) {
+      for (int iY = 0; iY < yN; iY++) {
+        latticeR[1] = iX;
+        latticeR[2] = iY;
+        physR = getPhysR(latticeR);
+        indicatorF(inside,&physR[0]);
+        if (inside[0]) {
+          fullCells++;
+          maxX = std::max(maxX, iX);
+          maxY = std::max(maxY, iY);
+          newX = std::min(newX, iX);
+          newY = std::min(newY, iY);
+        }
+      }
+    }
+    if (fullCells > 0) {
+      get(iC).setWeight(fullCells);
+      _cuboids[iC].resize(newX, newY, maxX - newX + 1, maxY - newY + 1);
+    }
+    else {
+      remove(iC);
+    }
+  }
+  // shrink mother cuboid
+  std::vector<T> minPhysR = getMinPhysR();
+  std::vector<T> maxPhysR = getMaxPhysR();
+  T minDelataR = getMinDeltaR();
+  _motherCuboid = Cuboid2D<T>(minPhysR[0], minPhysR[1], minDelataR, (int)((maxPhysR[0]-minPhysR[0])/minDelataR + 0.5), (int)((maxPhysR[1]-minPhysR[1])/minDelataR + 0.5));
+}
+
+template<typename T>
+void CuboidGeometry2D<T>::split(int iC, int p)
+{
+
+  Cuboid2D<T> temp(_cuboids[iC].get_globPosX(), _cuboids[iC].get_globPosY(),
+                   _cuboids[iC].getDeltaR(), _cuboids[iC].getNx(), _cuboids[iC].getNy());
+  temp.divide(p, _cuboids);
+  remove(iC);
+}
+
+template<typename T>
+void CuboidGeometry2D<T>::getNeighbourhood(int cuboid, std::vector<int> neighbours, int offset)
+{
+  for (int iC = 0; iC < getNc(); iC++) {
+    if (cuboid == iC) {
+      continue;
+    }
+    T globX = get(iC).get_globPosX();
+    T globY = get(iC).get_globPosY();
+    T nX = get(iC).getNx();
+    T nY = get(iC).getNy();
+    if (get(cuboid).checkInters(globX, globX + nX, globY, globY + nY, offset)) {
+      neighbours.push_back(iC);
+    }
+  }
+}
+
+template<typename T>
+void CuboidGeometry2D<T>::refineArea(T x0, T x1, T y0, T y1, int coarse_level)
+{
+
+  for (int iC = 0; iC < getNc(); iC++) {
+    if (get(iC).get_refinementLevel() != coarse_level) {
+      continue;
+    }
+    int locX0, locX1, locY0, locY1;
+    bool inter = get(iC).checkInters(x0, y1, y0, y1, locX0, locX1, locY0, locY1, 0);
+    if (!inter) {
+      continue;
+    }
+
+    T globX = get(iC).get_globPosX();
+    T globY = get(iC).get_globPosY();
+    T delta = get(iC).getDeltaR();
+    int nx = get(iC).getNx();
+    int ny = get(iC).getNy();
+
+    if (locX0 != 0) {
+      Cuboid2D<T> right_side(globX, globY, delta, locX0, ny, coarse_level);
+      add(right_side);
+    }
+
+    if (locY0 != 0) {
+      Cuboid2D<T> down_side(globX + locX0 * delta, globY, delta, nx - locX0, locY0, coarse_level);
+      add(down_side);
+    }
+
+    if (locX1 != get(iC).getNx() - 1) {
+      Cuboid2D<T> left_side(globX + (locX1 + 1)*delta, globY + locY0 * delta, delta, nx - locX1 - 1, ny - locY0, coarse_level);
+      add(left_side);
+    }
+
+    if (locY1 != get(iC).getNy() - 1) {
+      Cuboid2D<T> top_side(globX + locX0 * delta, globY + (locY1 + 1)*delta, delta, locX1 - locX0 + 1, ny - locY1 - 1, coarse_level);
+      add(top_side);
+    }
+    get(iC).init(globX + locX0 * delta, globY + locY0 * delta, delta, locX1 - locX0 + 1, locY1 - locY0 + 1, coarse_level);
+    get(iC).refineIncrease();
+  }
+}
+
+template<typename T>
+void CuboidGeometry2D<T>::print() const
+{
+  clout << "---Cuboid Stucture Statistics---" << std::endl;
+  clout << " Number of Cuboids: " << "\t" << getNc() << std::endl;
+  clout << " Delta (min): " << "\t" << "\t" << getMinDeltaR() << std::endl;
+  clout << "       (max): " << "\t" << "\t" << getMaxDeltaR() << std::endl;
+  clout << " Ratio (min): " << "\t" << "\t" << getMinRatio() << std::endl;
+  clout << "       (max): " << "\t" << "\t" << getMaxRatio() << std::endl;
+  clout << " Nodes (min): " << "\t" << "\t" << getMinLatticeVolume() << std::endl;
+  clout << "       (max): " << "\t" << "\t" << getMaxLatticeVolume() << std::endl;
+  clout << "--------------------------------" << std::endl;
+}
+
+template<typename T>
+void CuboidGeometry2D<T>::printExtended()
+{
+  clout << "Mothercuboid :" << std::endl;
+  getMotherCuboid().print();
+
+  for (int iC = 0; iC < getNc(); iC++) {
+    clout << "Cuboid #" << iC << ": " << std::endl;
+    get(iC).print();
+  }
+}
+
+}  // namespace olb
+
+#endif