Initialize at openlb-1-3

1 files changed, 1229 insertions, 0 deletions

diff --git a/src/functors/lattice/blockLatticeLocalF2D.hh b/src/functors/lattice/blockLatticeLocalF2D.hh
new file mode 100644
index 0000000..2dbcfed
--- /dev/null
+++ b/src/functors/lattice/blockLatticeLocalF2D.hh
@@ -0,0 +1,1229 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2012 Lukas Baron, Mathias J. Krause, Albert Mink
+ *  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.
+*/
+
+#ifndef BLOCK_LATTICE_LOCAL_F_2D_HH
+#define BLOCK_LATTICE_LOCAL_F_2D_HH
+
+#include<vector>
+#include<cmath>
+
+#include "blockLatticeLocalF2D.h"
+#include "blockBaseF2D.h"
+#include "functors/genericF.h"
+#include "functors/analytical/analyticalF.h"
+#include "functors/analytical/indicator/indicatorF2D.h"
+#include "core/blockLattice2D.h"
+#include "communication/mpiManager.h"
+#include "core/blockLatticeStructure2D.h"
+#include "dynamics/lbHelpers.h"  // for computation of lattice rho and velocity
+
+namespace olb {
+
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticeDissipation2D<T,DESCRIPTOR>::BlockLatticeDissipation2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticeF2D<T,DESCRIPTOR>(blockLattice,1), _converter(converter)
+{
+  this->getName() = "dissipation";
+}
+
+// todo: get functor working
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticeDissipation2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  //  int globIC = input[0];
+  //  int locix= input[1];
+  //  int lociy= input[2];
+  //  int lociz= input[3];
+  //  if ( this->_blockLattice.get_load().rank(globIC) == singleton::mpi().getRank() ) {
+  //    // local coords are given, fetch local cell and compute value(s)
+  //    T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR >::n];
+  //    int overlap = this->_blockLattice.getOverlap();
+  //    this->_blockLattice.getExtendedBlockLattice(this->_blockLattice.get_load().loc(globIC)).get(locix+overlap, lociy+overlap, lociz+overlap).computeAllMomenta(rho, uTemp, pi);
+
+  //    T PiNeqNormSqr = pi[0]*pi[0] + 2.*pi[1]*pi[1] + pi[2]*pi[2];
+  //    if (util::TensorVal<DESCRIPTOR >::n == 6)
+  //      PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2.*pi[4]*pi[4] +pi[5]*pi[5];
+
+  //    T nuLattice = converter.getLatticeNu();
+  //    T omega = converter.getOmega();
+  //    T finalResult = PiNeqNormSqr*nuLattice*pow(omega*descriptors::invCs2<T,DESCRIPTOR>(),2)/rho/2.;
+
+  //    return std::vector<T>(1,finalResult);
+  //  } else {
+  //    return std::vector<T>(); // empty vector
+  //  }
+
+  return false;
+}
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysDissipation2D<T,DESCRIPTOR>::BlockLatticePhysDissipation2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice,converter,1)
+{
+  this->getName() = "physDissipation";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysDissipation2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR >::n];
+  this->_blockLattice.get( input[0], input[1] ).computeAllMomenta(rho, uTemp, pi);
+
+  T PiNeqNormSqr = pi[0]*pi[0] + 2.*pi[1]*pi[1] + pi[2]*pi[2];
+  if (util::TensorVal<DESCRIPTOR >::n == 6) {
+    PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2.*pi[4]*pi[4] + pi[5]*pi[5];
+  }
+
+  T nuLattice = this->_converter.getLatticeViscosity();
+  T omega = 1. / this->_converter.getLatticeRelaxationTime();
+  T dt = this->_converter.getConversionFactorTime();
+  output[0] = PiNeqNormSqr*nuLattice*pow(omega*descriptors::invCs2<T,DESCRIPTOR>()/rho,2)/2.*this->_converter.getPhysViscosity()/this->_converter.getLatticeViscosity()/dt/dt;
+
+  return true;
+}
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticeDensity2D<T,DESCRIPTOR>::BlockLatticeDensity2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice)
+  : BlockLatticeF2D<T,DESCRIPTOR>(blockLattice,1)
+{
+  this->getName() = "density";
+}
+
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticeDensity2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  output[0] = this->_blockLattice.get( input[0], input[1] ).computeRho();
+  return true;
+}
+
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticeVelocity2D<T,DESCRIPTOR>::BlockLatticeVelocity2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice)
+  : BlockLatticeF2D<T,DESCRIPTOR>(blockLattice,2)
+{
+  this->getName() = "velocity";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticeVelocity2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  T rho;
+  this->_blockLattice.get(input[0], input[1]).computeRhoU(rho,output);
+  return true;
+}
+
+/*template <typename T, typename DESCRIPTOR>
+BlockLatticeStrainRate2D<T,DESCRIPTOR>::BlockLatticeStrainRate2D
+  (BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice, const UnitConverter<T>& converter)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice,converter,4)
+{ this->getName() = "strainRate"; }
+
+template <typename T, typename DESCRIPTOR>
+std::vector<T> BlockLatticeStrainRate2D<T,DESCRIPTOR>::operator() (std::vector<int> input) {
+
+  T strainRate[4];
+  T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR >::n];
+  this->_blockLattice.get( input[0], input[1] ).computeAllMomenta(rho, uTemp, pi);
+
+  T omega = this->_converter.getOmega();
+
+  strainRate[0] = -pi[0]*omega*descriptors::invCs2<T,DESCRIPTOR>()/rho/2.;
+  strainRate[1] = -pi[1]*omega*descriptors::invCs2<T,DESCRIPTOR>()/rho/2.;
+  strainRate[2] = -pi[1]*omega*descriptors::invCs2<T,DESCRIPTOR>()/rho/2.;
+  strainRate[3] = -pi[2]*omega*descriptors::invCs2<T,DESCRIPTOR>()/rho/2.;
+
+  //cout << pi[0] << " " << pi[1] << " " << pi[2] << " " << descriptors::invCs2<T,DESCRIPTOR>() << endl;
+
+  std::vector<T> output(strainRate, strainRate+4); // first adress, last adress
+  return output;
+}*/
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysStrainRate2D<T,DESCRIPTOR>::BlockLatticePhysStrainRate2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice,converter,4)
+{
+  this->getName() = "strainRate";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysStrainRate2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR >::n];
+  this->_blockLattice.get( input[0], input[1] ).computeAllMomenta(rho, uTemp, pi);
+
+  T omega = 1. / this->_converter.getLatticeRelaxationTime();
+  T dt = this->_converter.getConversionFactorTime();
+
+  output[0] = -pi[0]*omega*descriptors::invCs2<T,DESCRIPTOR>()/rho/2./dt;
+  output[1] = -pi[1]*omega*descriptors::invCs2<T,DESCRIPTOR>()/rho/2./dt;
+  output[2] = -pi[1]*omega*descriptors::invCs2<T,DESCRIPTOR>()/rho/2./dt;
+  output[3] = -pi[2]*omega*descriptors::invCs2<T,DESCRIPTOR>()/rho/2./dt;
+
+  return true;
+}
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticeGeometry2D<T,DESCRIPTOR>::BlockLatticeGeometry2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice, BlockGeometryStructure2D<T>& blockGeometry, int material)
+  : BlockLatticeF2D<T,DESCRIPTOR>(blockLattice,1), _blockGeometry(blockGeometry), _material(material)
+{
+  this->getName() = "geometry";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticeGeometry2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  int materialTmp = _blockGeometry.getMaterial( input[0], input[1] );
+
+  if (_material != -1) {
+    if (_material == materialTmp) {
+      output[0] = T(1);
+      return true;
+    }
+    else {
+      output[0] = T();
+      return true;
+    }
+  }
+  output[0]=T(materialTmp);
+  return false;
+}
+
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticeRank2D<T,DESCRIPTOR>::BlockLatticeRank2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice)
+  : BlockLatticeF2D<T,DESCRIPTOR>(blockLattice,1)
+{
+  this->getName() = "rank";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticeRank2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  output[0] = singleton::mpi().getRank() + 1;
+  return false;
+}
+
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticeCuboid2D<T,DESCRIPTOR>::BlockLatticeCuboid2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice, const int iC)
+  : BlockLatticeF2D<T,DESCRIPTOR>(blockLattice,1), _iC(iC)
+{
+  this->getName() = "cuboid";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticeCuboid2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  output[0] = _iC + 1;
+  return false;
+}
+
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysPressure2D<T,DESCRIPTOR>::BlockLatticePhysPressure2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice,converter,1)
+{
+  this->getName() = "physPressure";
+}
+
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysPressure2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  // lattice pressure = c_s^2 ( rho -1 )
+  T latticePressure = ( this->_blockLattice.get( input[0], input[1] ).computeRho() - 1.0 ) / descriptors::invCs2<T,DESCRIPTOR>();
+  output[0] = this->_converter.getPhysPressure(latticePressure);
+
+  return true;
+}
+
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysVelocity2D<T,DESCRIPTOR>::BlockLatticePhysVelocity2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice,converter,2)
+{
+  this->getName() = "physVelocity";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysVelocity2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  T rho;
+  this->_blockLattice.get( input[0], input[1] ).computeRhoU(rho,output);
+  output[0]=this->_converter.getPhysVelocity(output[0]);
+  output[1]=this->_converter.getPhysVelocity(output[1]);
+  return true;
+}
+
+
+template<typename T, typename DESCRIPTOR, typename FIELD>
+BlockLatticeField2D<T,DESCRIPTOR,FIELD>::BlockLatticeField2D(
+  BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice)
+  : BlockLatticeF2D<T, DESCRIPTOR>(blockLattice, DESCRIPTOR::template size<FIELD>())
+{
+  this->getName() = "extField";
+}
+
+template<typename T, typename DESCRIPTOR, typename FIELD>
+bool BlockLatticeField2D<T,DESCRIPTOR,FIELD>::operator()(
+  T output[], const int input[])
+{
+  this->_blockLattice.get(input[0], input[1]).template computeField<FIELD>(output);
+  return true;
+}
+
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysExternalPorosity2D<T,DESCRIPTOR>::BlockLatticePhysExternalPorosity2D(
+  BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice,
+  int overlap,
+  const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice, converter, 2),
+    _overlap(overlap)
+{
+  this->getName() = "ExtPorosityField";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysExternalPorosity2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  this->_blockLattice.get(
+    input[0]+_overlap, input[1]+_overlap
+  ).template computeField<descriptors::POROSITY>(output);
+  return true;
+}
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysExternalVelocity2D<T,DESCRIPTOR>::BlockLatticePhysExternalVelocity2D(
+  BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice,
+  int overlap,
+  const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice, converter, 2),
+    _overlap(overlap)
+{
+  this->getName() = "ExtVelocityField";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysExternalVelocity2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  this->_blockLattice.get(
+    input[0]+_overlap, input[1]+_overlap
+  ).template computeField<descriptors::VELOCITY>(output);
+  output[0] = this->_converter.getPhysVelocity(output[0]);
+  output[1] = this->_converter.getPhysVelocity(output[1]);
+  return true;
+}
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysExternalParticleVelocity2D<T,DESCRIPTOR>::BlockLatticePhysExternalParticleVelocity2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice,converter,2)
+{
+  this->getName() = "ExtParticleVelocityField";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysExternalParticleVelocity2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  const T* velocity_numerator   = this->blockLattice.get(input).template getFieldPointer<descriptors::VELOCITY_NUMERATOR>();
+  const T* velocity_denominator = this->blockLattice.get(input).template getFieldPointer<descriptors::VELOCITY_DENOMINATOR>();
+
+  if (velocity_denominator[0] > std::numeric_limits<T>::epsilon()) {
+    output[0]=this->_converter.getPhysVelocity(velocity_numerator[0]/velocity_denominator[0]);
+    output[1]=this->_converter.getPhysVelocity(velocity_numerator[1]/velocity_denominator[0]);
+    return true;
+  }
+  output[0]=this->_converter.getPhysVelocity(velocity_numerator[0]);
+  output[1]=this->_converter.getPhysVelocity(velocity_numerator[1]);
+  return true;
+}
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysWallShearStress2D<T,DESCRIPTOR>::BlockLatticePhysWallShearStress2D(
+    BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice,
+    BlockGeometryStructure2D<T>& blockGeometry,
+    int overlap,
+    int material,
+    const UnitConverter<T,DESCRIPTOR>& converter,
+    IndicatorF2D<T>& indicator)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice,converter,1),
+    _blockGeometry(blockGeometry), _overlap(overlap), _material(material)
+{
+  this->getName() = "physWallShearStress";
+  const T scaling = this->_converter.getConversionFactorLength() * 0.1;
+  const T omega = 1. / this->_converter.getLatticeRelaxationTime();
+  const T dt = this->_converter.getConversionFactorTime();
+  _physFactor = -omega * descriptors::invCs2<T,DESCRIPTOR>() / dt * this->_converter.getPhysDensity() * this->_converter.getPhysViscosity();
+  std::vector<int> discreteNormalOutwards(3, 0);
+
+  for (int iX = 1; iX < _blockGeometry.getNx() - 1; iX++) {
+    _discreteNormal.resize(_blockGeometry.getNx() - 2);
+    _normal.resize(_blockGeometry.getNx() - 2);
+
+    for (int iY = 1; iY < _blockGeometry.getNy() - 1; iY++) {
+      _discreteNormal[iX-1].resize(_blockGeometry.getNy() - 2);
+      _normal[iX-1].resize(_blockGeometry.getNy() - 2);
+
+      if (_blockGeometry.get(iX, iY) == _material) {
+        discreteNormalOutwards = _blockGeometry.getStatistics().getType(iX, iY);
+        _discreteNormal[iX-1][iY-1].resize(2);
+        _normal[iX-1][iY-1].resize(2);
+
+        _discreteNormal[iX-1][iY-1][0] = -discreteNormalOutwards[1];
+        _discreteNormal[iX-1][iY-1][1] = -discreteNormalOutwards[2];
+
+        T physR[2];
+        _blockGeometry.getPhysR(physR, iX, iY);
+        Vector<T,2> origin(physR[0],physR[1]);
+        Vector<T,2> direction(-_discreteNormal[iX-1][iY-1][0] * scaling,
+                              -_discreteNormal[iX-1][iY-1][1] * scaling);
+        Vector<T,2> normal(0.,0.);
+        origin[0] = physR[0];
+        origin[1] = physR[1];
+
+        indicator.normal(normal, origin, direction);
+        normal.normalize();
+
+        _normal[iX-1][iY-1][0] = normal[0];
+        _normal[iX-1][iY-1][1] = normal[0];
+      }
+    }
+  }
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysWallShearStress2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  output[0] = T();
+
+  if (input[0] + _overlap < 1 ||
+      input[1] + _overlap < 1 ||
+      input[0] + _overlap >= _blockGeometry.getNx()-1 ||
+      input[1] + _overlap >= _blockGeometry.getNy()-1){
+#ifdef OLB_DEBUG
+    std::cout << "Input address not mapped by _discreteNormal, overlap too small" << std::endl;
+#endif
+    return true;
+  }
+
+  if (this->_blockGeometry.get(input[0]+_overlap,input[1]+_overlap)==_material) {
+
+    T traction[2];
+    T stress[3];
+    T rho = this->_blockLattice.get(input[0] + _overlap + _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][0],
+                                    input[1] + _overlap + _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][1]).computeRho();
+
+    this->_blockLattice.get(input[0] + _overlap +   _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][0],
+                            input[1] + _overlap +   _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][1]).computeStress(stress);
+
+    traction[0] = stress[0]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][0] +
+                  stress[1]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][1];
+    traction[1] = stress[1]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][0] +
+                  stress[2]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][1];
+
+    T traction_normal_SP;
+    T tractionNormalComponent[2];
+    // scalar product of traction and normal vector
+    traction_normal_SP = traction[0] * _normal[input[0]+_overlap-1][input[1]+_overlap-1][0] +
+                           traction[1] * _normal[input[0]+_overlap-1][input[1]+_overlap-1][1];
+    tractionNormalComponent[0] = traction_normal_SP * _normal[input[0]+_overlap-1][input[1]+_overlap-1][0];
+    tractionNormalComponent[1] = traction_normal_SP * _normal[input[0]+_overlap-1][input[1]+_overlap-1][1];
+
+    T WSS[2];
+    WSS[0] = traction[0] - tractionNormalComponent[0];
+    WSS[1] = traction[1] - tractionNormalComponent[1];
+    // magnitude of the wall shear stress vector
+    output[0] = sqrt(WSS[0]*WSS[0] + WSS[1]*WSS[1]);
+    return true;
+  }
+  else {
+    return true;
+  }
+}
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysBoundaryForce2D<T,DESCRIPTOR>::BlockLatticePhysBoundaryForce2D(
+  BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice,
+  BlockIndicatorF2D<T>&                  indicatorF,
+  const UnitConverter<T,DESCRIPTOR>&     converter)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice, converter, 2),
+    _indicatorF(indicatorF),
+    _blockGeometry(indicatorF.getBlockGeometryStructure())
+{
+  this->getName() = "physBoundaryForce";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysBoundaryForce2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  for (int i = 0; i < this->getTargetDim(); ++i) {
+    output[i] = T();
+  }
+
+  if (_indicatorF(input)) {
+    for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) {
+      // Get direction
+      // Get next cell located in the current direction
+      // Check if the next cell is a fluid node
+      if (_blockGeometry.get(input[0] + descriptors::c<DESCRIPTOR >(iPop,0), input[1] + descriptors::c<DESCRIPTOR >(iPop,1)) == 1) {
+        // Get f_q of next fluid cell where l = opposite(q)
+        T f = this->_blockLattice.get(input[0] + descriptors::c<DESCRIPTOR >(iPop,0), input[1] + descriptors::c<DESCRIPTOR >(iPop,1))[iPop];
+        // Get f_l of the boundary cell
+        // Add f_q and f_opp
+        f += this->_blockLattice.get(input[0], input[1])[util::opposite<DESCRIPTOR >(iPop)];
+        // Update force
+        for (int i = 0; i < this->getTargetDim(); ++i) {
+          output[i] -= descriptors::c<DESCRIPTOR >(iPop,i) * f;
+        }
+      }
+    }
+    for (int i = 0; i < this->getTargetDim(); ++i) {
+      output[i] = this->_converter.getPhysForce(output[i]);
+    }
+  }
+  return true;
+}
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePSMPhysForce2D<T,DESCRIPTOR>::BlockLatticePSMPhysForce2D(
+  BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice,
+  const UnitConverter<T,DESCRIPTOR>&     converter,
+  int mode_)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice, converter, 2)
+{
+  this->getName() = "physPSMForce";
+  mode = (Mode) mode_;
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticePSMPhysForce2D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  for (int i = 0; i < this->getTargetDim(); ++i) {
+    output[i] = T();
+  }
+
+  T epsilon = 1. - *(this->_blockLattice.get(input[0], input[1]).template getFieldPointer<descriptors::POROSITY>());
+
+  //if ((epsilon > 1e-5 && epsilon < 1 - 1e-5)) {
+  if ((epsilon > 1e-5)) {
+    T rho, u[DESCRIPTOR::d], u_s[DESCRIPTOR::d];
+
+    for (int i = 0; i < DESCRIPTOR::d; i++) {
+      u_s[i] = (this->_blockLattice.get(input[0], input[1]).template getFieldPointer<descriptors::VELOCITY_SOLID>())[i];
+    }
+    T paramA = this->_converter.getLatticeRelaxationTime() - 0.5;
+    // speed up paramB
+    T paramB = (epsilon * paramA) / ((1. - epsilon) + paramA);
+
+    T omega_s;
+    T omega = 1. / this->_converter.getLatticeRelaxationTime();
+
+    this->_blockLattice.get(input[0], input[1]).computeRhoU(rho, u);
+
+    const T uSqr_s = util::normSqr<T,DESCRIPTOR::d>(u_s);
+    T uSqr = util::normSqr<T,DESCRIPTOR::d>(u);
+    for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+      switch (mode) {
+        case M2:
+          omega_s = (lbDynamicsHelpers<T, DESCRIPTOR>::equilibrium(iPop, rho, u_s, uSqr_s)
+              - this->_blockLattice.get(input[0], input[1])[iPop])
+              + (1 - omega)
+                  * (this->_blockLattice.get(input[0], input[1])[iPop]
+                      - lbDynamicsHelpers<T, DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr));
+          break;
+        case M3:
+          omega_s =
+              (this->_blockLattice.get(input[0], input[1])[descriptors::opposite<DESCRIPTOR>(iPop)]
+                  - lbDynamicsHelpers<T, DESCRIPTOR>::equilibrium(
+                      descriptors::opposite<DESCRIPTOR>(iPop), rho, u_s, uSqr_s))
+                  - (this->_blockLattice.get(input[0], input[1])[iPop]
+                      - lbDynamicsHelpers<T, DESCRIPTOR>::equilibrium(iPop, rho, u_s, uSqr_s));
+      }
+
+      for (int i = 0; i < this->getTargetDim(); ++i) {
+        output[i] -= descriptors::c<DESCRIPTOR>(iPop,i) * omega_s;
+      }
+    }
+
+    for (int i = 0; i < this->getTargetDim(); ++i) {
+      output[i] = this->_converter.getPhysForce(output[i] * paramB);
+    }
+  }
+  return true;
+}
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysCorrBoundaryForce2D<T,DESCRIPTOR>::BlockLatticePhysCorrBoundaryForce2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice,BlockGeometry2D<T>& blockGeometry,
+ int material, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice,converter,2),
+    _blockGeometry(blockGeometry), _material(material)
+{
+  this->getName() = "physCorrBoundaryForce";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysCorrBoundaryForce2D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  //  int globIC = input[0];
+  //  int locix= input[1];
+  //  int lociy= input[2];
+  //  int lociz= input[3];
+
+  //  std::vector<T> force(3, T());
+  //  if ( this->_blockLattice.get_load().rank(globIC) == singleton::mpi().getRank() ) {
+  //    int globX = (int)this->_blockLattice.get_cGeometry().get(globIC).get_globPosX() + locix;
+  //    int globY = (int)this->_blockLattice.get_cGeometry().get(globIC).get_globPosY() + lociy;
+  //    int globZ = (int)this->_blockLattice.get_cGeometry().get(globIC).get_globPosZ() + lociz;
+  //    if(BlockGeometry.getMaterial(globX, globY, globZ) == material) {
+  //      for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) {
+  //        // Get direction
+  //        const int* c = DESCRIPTOR::c[iPop];
+  //        // Get next cell located in the current direction
+  //        // Check if the next cell is a fluid node
+  //        if (BlockGeometry.getMaterial(globX + c[0], globY + c[1], globZ + c[2]) == 1) {
+  //          int overlap = this->_blockLattice.getOverlap();
+  //          // Get f_q of next fluid cell where l = opposite(q)
+  //          T f = this->_blockLattice.getExtendedBlockLattice(this->_blockLattice.get_load().loc(globIC)).get(locix+overlap + c[0], lociy+overlap + c[1], lociz+overlap + c[2])[iPop];
+  //          // Get f_l of the boundary cell
+  //          // Add f_q and f_opp
+  //          f += this->_blockLattice.getExtendedBlockLattice(this->_blockLattice.get_load().loc(globIC)).get(locix+overlap, lociy+overlap, lociz+overlap)[util::opposite<DESCRIPTOR >(iPop)];
+  //          // Update force
+  //          force[0] -= c[0]*(f-2.*descriptors::t<T,DESCRIPTOR>(iPop));
+  //          force[1] -= c[1]*(f-2.*descriptors::t<T,DESCRIPTOR>(iPop));
+  //          force[2] -= c[2]*(f-2.*descriptors::t<T,DESCRIPTOR>(iPop));
+  //        }
+  //      }
+  //      force[0]=this->_converter.physForce(force[0]);
+  //      force[1]=this->_converter.physForce(force[1]);
+  //      force[2]=this->_converter.physForce(force[2]);
+  //      return force;
+  //    }
+  //    else {
+  //      return force;
+  //    }
+  //  }
+  //  else {
+  //    return force;
+  //  }
+  return false;
+}
+
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePorosity2D<T,DESCRIPTOR>::BlockLatticePorosity2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice, BlockGeometryStructure2D<T>& blockGeometry,
+ int material, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice,converter,1), _blockGeometry(blockGeometry), _material(material)
+{
+  this->getName() = "porosity";
+}
+
+// under construction
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePorosity2D<T,DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  return false;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeVolumeFractionApproximation2D<T, DESCRIPTOR>::BlockLatticeVolumeFractionApproximation2D(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice,
+    BlockGeometryStructure2D<T>& blockGeometry,
+    IndicatorF2D<T>& indicator,
+    int refinementLevel,
+    const UnitConverter<T,DESCRIPTOR>& converter,
+    bool insideOut)
+  : BlockLatticeF2D<T, DESCRIPTOR>(blockLattice, 1),
+    _blockGeometry(blockGeometry), _indicator(indicator), _refinementLevel(refinementLevel), _converter(converter), _insideOut(insideOut),
+    _physSubGridMinPhysRshift((1./ T(_refinementLevel * 2.) - 0.5) * _converter.getPhysDeltaX()),
+    _physSubGridDeltaX(1. / T(_refinementLevel) * _converter.getPhysDeltaX()),
+    _latticeSubGridVolume(1. / (_refinementLevel * _refinementLevel))
+{
+  this->getName() = "volumeFractionApproximation";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeVolumeFractionApproximation2D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  output[0] = 0.;
+  T physR[2];
+  bool inside[1];
+  _blockGeometry.getPhysR(physR, input[0], input[1]);
+
+  T subGridMinPhysR[2];
+  subGridMinPhysR[0] = physR[0] + _physSubGridMinPhysRshift;
+  subGridMinPhysR[1] = physR[1] + _physSubGridMinPhysRshift;
+
+  for (int x = 0; x < _refinementLevel; x++) {
+    for (int y = 0; y < _refinementLevel; y++) {
+      physR[0] = subGridMinPhysR[0] + x * _physSubGridDeltaX;
+      physR[1] = subGridMinPhysR[1] + y * _physSubGridDeltaX;
+      _indicator(inside, physR);
+      if (!_insideOut) {
+        if (inside[0]) {
+          output[0] += _latticeSubGridVolume;
+        }
+      }
+      else {
+        if (!inside[0]) {
+          output[0] += _latticeSubGridVolume;
+        }
+      }
+    }
+  }
+  return true;
+}
+
+
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeVolumeFractionPolygonApproximation2D<T, DESCRIPTOR>::BlockLatticeVolumeFractionPolygonApproximation2D(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice,
+    BlockGeometryStructure2D<T>& blockGeometry,
+    IndicatorF2D<T>& indicator,
+    const UnitConverter<T,DESCRIPTOR>& converter,
+    bool insideOut)
+  : BlockLatticeF2D<T, DESCRIPTOR>(blockLattice, 1),
+    _blockGeometry(blockGeometry), _indicator(indicator), _converter(converter), _insideOut(insideOut)
+{
+  this->getName() = "volumeFractionPolygonApproximation";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeVolumeFractionPolygonApproximation2D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  output[0] = 0.;
+  T physR[2];
+
+  bool cornerXMYM_inside[1];
+  bool cornerXMYP_inside[1];
+  bool cornerXPYM_inside[1];
+  bool cornerXPYP_inside[1];
+  _blockGeometry.getPhysR(physR, input[0], input[1]);
+
+  T cornerXMYM[2];
+  T cornerXMYP[2];
+  T cornerXPYM[2];
+  T cornerXPYP[2];
+
+  cornerXMYM[0] = physR[0] - 0.5*_converter.getPhysDeltaX();
+  cornerXMYM[1] = physR[1] - 0.5*_converter.getPhysDeltaX();
+
+  cornerXMYP[0] = physR[0] - 0.5*_converter.getPhysDeltaX();
+  cornerXMYP[1] = physR[1] + 0.5*_converter.getPhysDeltaX();
+
+  cornerXPYM[0] = physR[0] + 0.5*_converter.getPhysDeltaX();
+  cornerXPYM[1] = physR[1] - 0.5*_converter.getPhysDeltaX();
+
+  cornerXPYP[0] = physR[0] + 0.5*_converter.getPhysDeltaX();
+  cornerXPYP[1] = physR[1] + 0.5*_converter.getPhysDeltaX();
+
+  _indicator(cornerXMYM_inside, cornerXMYM);
+  _indicator(cornerXMYP_inside, cornerXMYP);
+  _indicator(cornerXPYM_inside, cornerXPYM);
+  _indicator(cornerXPYP_inside, cornerXPYP);
+
+  if ((cornerXMYM_inside[0] && cornerXMYP_inside[0] && cornerXPYM_inside[0] && cornerXPYP_inside[0]) ||
+      (!cornerXMYM_inside[0] && !cornerXMYP_inside[0] && !cornerXPYM_inside[0] && !cornerXPYP_inside[0]) ) {
+    if (!_insideOut) {
+      if (cornerXMYM_inside[0]) {
+        output[0] = 1.0;
+      }
+    }
+    else {
+      if (!cornerXMYM_inside[0]) {
+        output[0] = 1.0;
+      }
+    }
+  } else {
+    Vector<T,2> cornerXMYM_vec(physR[0] - 0.5*_converter.getPhysDeltaX(), physR[1] - 0.5*_converter.getPhysDeltaX());
+    Vector<T,2> cornerXPYP_vec(physR[0] + 0.5*_converter.getPhysDeltaX(), physR[1] + 0.5*_converter.getPhysDeltaX());
+
+    Vector<T,2> directionXP(_converter.getPhysDeltaX(), 0);
+    Vector<T,2> directionXM(-_converter.getPhysDeltaX(), 0);
+    Vector<T,2> directionYP(0, _converter.getPhysDeltaX());
+    Vector<T,2> directionYM(0, -_converter.getPhysDeltaX());
+
+    T distanceXP = 1.01 * _converter.getPhysDeltaX();
+    T distanceXM = 1.01 * _converter.getPhysDeltaX();
+    T distanceYP = 1.01 * _converter.getPhysDeltaX();
+    T distanceYM = 1.01 * _converter.getPhysDeltaX();
+
+    if ((cornerXMYM_inside[0] && !cornerXMYP_inside[0]) ||
+        (!cornerXMYM_inside[0] && cornerXMYP_inside[0]) ) {
+      _indicator.distance(distanceYP, cornerXMYM, directionYP);
+    }
+
+    if ((cornerXMYM_inside[0] && !cornerXPYM_inside[0]) ||
+        (!cornerXMYM_inside[0] && cornerXPYM_inside[0]) ) {
+      _indicator.distance(distanceXP, cornerXMYM, directionXP);
+    }
+
+    if ((cornerXPYP_inside[0] && !cornerXMYP_inside[0]) ||
+        (!cornerXPYP_inside[0] && cornerXMYP_inside[0]) ) {
+      _indicator.distance(distanceXM, cornerXPYP, directionXM);
+    }
+
+    if ((cornerXPYP_inside[0] && !cornerXPYM_inside[0]) ||
+        (!cornerXPYP_inside[0] && cornerXPYM_inside[0]) ) {
+      _indicator.distance(distanceYM, cornerXPYP, directionYM);
+    }
+
+    T volumeFraction = 0.0;
+
+    if (distanceXP < _converter.getPhysDeltaX() && distanceXM < _converter.getPhysDeltaX()) {
+      volumeFraction = distanceXP * _converter.getPhysDeltaX();
+      volumeFraction += 0.5 * (_converter.getPhysDeltaX() - distanceXM - distanceXP) * _converter.getPhysDeltaX();
+      volumeFraction /= _converter.getPhysDeltaX() * _converter.getPhysDeltaX();
+      if (!cornerXMYM_inside[0]) {
+        volumeFraction = 1.0 - volumeFraction;
+      }
+    }
+
+    if (distanceYP < _converter.getPhysDeltaX() && distanceYM < _converter.getPhysDeltaX()) {
+      volumeFraction = distanceYP * _converter.getPhysDeltaX();
+      volumeFraction += 0.5 * (_converter.getPhysDeltaX() - distanceYM - distanceYP) * _converter.getPhysDeltaX();
+      volumeFraction /= _converter.getPhysDeltaX() * _converter.getPhysDeltaX();
+      if (!cornerXMYM_inside[0]) {
+        volumeFraction = 1.0 - volumeFraction;
+      }
+    }
+
+    if (distanceXP < _converter.getPhysDeltaX() && distanceYP < _converter.getPhysDeltaX()) {
+      volumeFraction = 0.5 * distanceXP * distanceYP;
+      volumeFraction /= _converter.getPhysDeltaX() * _converter.getPhysDeltaX();
+      if (!cornerXMYM_inside[0]) {
+        volumeFraction = 1.0 - volumeFraction;
+      }
+    }
+
+    if (distanceXM < _converter.getPhysDeltaX() && distanceYM < _converter.getPhysDeltaX()) {
+      volumeFraction = 0.5 * distanceXM * distanceYM;
+      volumeFraction /= _converter.getPhysDeltaX() * _converter.getPhysDeltaX();
+      if (!cornerXPYP_inside[0]) {
+        volumeFraction = 1.0 - volumeFraction;
+      }
+    }
+
+    if (distanceXM < _converter.getPhysDeltaX() && distanceYP < _converter.getPhysDeltaX()) {
+      volumeFraction = 0.5 * (_converter.getPhysDeltaX() - distanceXM) * (_converter.getPhysDeltaX() - distanceYP);
+      volumeFraction /= _converter.getPhysDeltaX() * _converter.getPhysDeltaX();
+      if (!cornerXMYP_inside[0]) {
+        volumeFraction = 1.0 - volumeFraction;
+      }
+    }
+
+    if (distanceYM < _converter.getPhysDeltaX() && distanceXP < _converter.getPhysDeltaX()) {
+      volumeFraction = 0.5 * (_converter.getPhysDeltaX() - distanceXP) * (_converter.getPhysDeltaX() - distanceYM);
+      volumeFraction /= _converter.getPhysDeltaX() * _converter.getPhysDeltaX();
+      if (!cornerXPYM_inside[0]) {
+        volumeFraction = 1.0 - volumeFraction;
+      }
+    }
+
+    if (!_insideOut) {
+        output[0] = volumeFraction;
+    }
+    else {
+        output[0] = 1.0 - volumeFraction;
+    }
+
+  }
+  return true;
+}
+
+
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysPermeability2D<T,DESCRIPTOR>::BlockLatticePhysPermeability2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice, BlockGeometryStructure2D<T>& blockGeometry,
+ int material, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF2D<T,DESCRIPTOR>(blockLattice,converter,1), _blockGeometry(blockGeometry), _material(material)
+{
+  this->getName() = "permeability";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysPermeability2D<T,DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  return false;
+}
+
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysDarcyForce2D<T,DESCRIPTOR>::BlockLatticePhysDarcyForce2D
+(BlockLatticeStructure2D<T,DESCRIPTOR>& blockLattice, BlockGeometry2D<T>& blockGeometry,
+ int material, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF2D&