Initialize at openlb-1-3

1 files changed, 1556 insertions, 0 deletions

diff --git a/src/functors/lattice/blockLatticeLocalF3D.hh b/src/functors/lattice/blockLatticeLocalF3D.hh
new file mode 100644
index 0000000..6359b55
--- /dev/null
+++ b/src/functors/lattice/blockLatticeLocalF3D.hh
@@ -0,0 +1,1556 @@
+/*  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_3D_HH
+#define BLOCK_LATTICE_LOCAL_F_3D_HH
+
+
+#include<cmath>
+#include<math.h>
+
+#include "blockLatticeLocalF3D.h"
+#include "blockBaseF3D.h"
+#include "core/blockLatticeStructure3D.h"
+#include "dynamics/lbHelpers.h"  // for computation of lattice rho and velocity
+#include "communication/mpiManager.h"
+#include "utilities/vectorHelpers.h"
+
+namespace olb {
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeFpop3D<T, DESCRIPTOR>::BlockLatticeFpop3D(BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, DESCRIPTOR::q)
+{
+  this->getName() = "fPop";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeFpop3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
+    output[iPop] =
+      this->_blockLattice.get(input[0], input[1], input[2])[iPop]
+      + descriptors::t<T,DESCRIPTOR>(iPop);
+  }
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeDissipation3D<T, DESCRIPTOR>::BlockLatticeDissipation3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1), _converter(converter)
+{
+  this->getName() = "dissipation";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeDissipation3D<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], input[2]).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.getLatticeViscosity();
+  T omega = 1. / _converter.getLatticeRelaxationTime();
+  output[0] = PiNeqNormSqr * nuLattice
+              * pow(omega * descriptors::invCs2<T,DESCRIPTOR>(), 2) / rho / 2.;
+
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticePhysDissipation3D<T, DESCRIPTOR>::BlockLatticePhysDissipation3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice,
+  int overlap,
+  const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1),
+    _overlap(overlap),
+    _converter(converter)
+{
+  this->getName() = "physDissipation";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticePhysDissipation3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR >::n];
+  this->_blockLattice.get(
+    input[0]+_overlap, input[1]+_overlap, input[2]+_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 = 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() / nuLattice / dt / dt;
+
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeEffevtiveDissipation3D<T, DESCRIPTOR>::BlockLatticeEffevtiveDissipation3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter, T smagoConst,
+  LESDynamics<T, DESCRIPTOR>& LESdynamics)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1),
+    _converter(converter), _smagoConst(smagoConst), _LESdynamics(LESdynamics)
+{
+  this->getName() = "EffevtiveDissipation";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeEffevtiveDissipation3D<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], input[2]).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 omegaEff = _LESdynamics.getEffectiveOmega(this->_blockLattice.get(input[0], input[1], input[2]));
+  T nuEff = ((1./omegaEff)-0.5)/descriptors::invCs2<T,DESCRIPTOR>();
+
+  output[0] = PiNeqNormSqr * (nuEff)
+              * pow(omegaEff * descriptors::invCs2<T,DESCRIPTOR>() / rho, 2)  / 2.;
+
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticePhysEffevtiveDissipation3D<T, DESCRIPTOR>::BlockLatticePhysEffevtiveDissipation3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter, T smagoConst,
+  LESDynamics<T, DESCRIPTOR>& LESdynamics)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1),
+    _converter(converter), _smagoConst(smagoConst), _LESdynamics(LESdynamics)
+{
+  this->getName() = "physEffevtiveDissipation";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticePhysEffevtiveDissipation3D<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], input[2]).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 dt = 1./ _converter.getConversionFactorTime();
+  T omegaEff = _LESdynamics.getEffectiveOmega(this->_blockLattice.get(input[0], input[1], input[2]));
+  T nuEff = ((1./omegaEff)-0.5)/descriptors::invCs2<T,DESCRIPTOR>();  // BGK shear viscosity definition
+
+  output[0] = PiNeqNormSqr * nuEff
+              * pow(omegaEff * descriptors::invCs2<T,DESCRIPTOR>() / rho, 2) / 2.
+              * _converter.getPhysViscosity() / _converter.getLatticeViscosity() / dt / dt;
+
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeDensity3D<T, DESCRIPTOR>::BlockLatticeDensity3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1)
+{
+  this->getName() = "density";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeDensity3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  output[0] = this->_blockLattice.get(input[0], input[1], input[2]).computeRho();
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeVelocity3D<T, DESCRIPTOR>::BlockLatticeVelocity3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 3)
+{
+  this->getName() = "velocity";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeVelocity3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  T rho;
+  this->_blockLattice.get(input[0], input[1], input[2]).computeRhoU(rho, output);
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeExternalVelocity3D<T, DESCRIPTOR>::BlockLatticeExternalVelocity3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 3)
+{
+  this->getName() = "externalVelocity";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeExternalVelocity3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  T* ExtVel = this->_blockLattice.get(input[0], input[1], input[2]).template getFieldPointer<descriptors::VELOCITY>();
+  for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) {
+    output[iVel] = ExtVel[iVel];
+  }
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeFlux3D<T, DESCRIPTOR>::BlockLatticeFlux3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 3)
+{
+  this->getName() = "flux";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeFlux3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  this->_blockLattice.get(input[0], input[1], input[2]).computeJ(output);
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeStrainRate3D<T, DESCRIPTOR>::BlockLatticeStrainRate3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF3D<T, DESCRIPTOR>(blockLattice, converter, 9)
+{
+  this->getName() = "strainRate";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeStrainRate3D<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], input[2]).computeAllMomenta(rho,
+      uTemp,
+      pi);
+
+  T omega = 1. / this->_converter.getLatticeRelaxationTime();
+
+  output[0] = -pi[0] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.;
+  output[1] = -pi[1] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.;
+  output[2] = -pi[2] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.;
+  output[3] = -pi[1] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.;
+  output[4] = -pi[3] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.;
+  output[5] = -pi[4] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.;
+  output[6] = -pi[2] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.;
+  output[7] = -pi[4] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.;
+  output[8] = -pi[5] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2.;
+
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticePhysStrainRate3D<T, DESCRIPTOR>::BlockLatticePhysStrainRate3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice,
+  int overlap,
+  const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF3D<T, DESCRIPTOR>(blockLattice, converter, 9),
+    _overlap(overlap)
+{
+  this->getName() = "strainRate";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticePhysStrainRate3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal<DESCRIPTOR >::n];
+  this->_blockLattice.get(
+    input[0]+_overlap, input[1]+_overlap, input[2]+_overlap
+  ).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[2] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt;
+  output[3] = -pi[1] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt;
+  output[4] = -pi[3] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt;
+  output[5] = -pi[4] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt;
+  output[6] = -pi[2] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt;
+  output[7] = -pi[4] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt;
+  output[8] = -pi[5] * omega * descriptors::invCs2<T,DESCRIPTOR>() / rho / 2. / dt;
+
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeGeometry3D<T, DESCRIPTOR>::BlockLatticeGeometry3D(BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice,
+    BlockGeometryStructure3D<T>& blockGeometry, int material)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1),
+    _blockGeometry(blockGeometry),
+    _material(material)
+{
+  this->getName() = "geometry";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeGeometry3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  output[0] = _blockGeometry.getMaterial(input[0], input[1], input[2]);
+
+  if (_material != -1) {
+    if ( util::nearZero(_material-output[0]) ) {
+      output[0] = 1.;
+      return true;
+    }
+    else {
+      output[0] = 0.;
+      return true;
+    }
+  }
+  return true;
+}
+
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeRank3D<T,DESCRIPTOR>::BlockLatticeRank3D(
+  BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice)
+  : BlockLatticeF3D<T,DESCRIPTOR>(blockLattice, 1)
+{
+  this->getName() = "rank";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeRank3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  output[0] = singleton::mpi().getRank() + 1;
+  return true;
+}
+
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeCuboid3D<T,DESCRIPTOR>::BlockLatticeCuboid3D(
+  BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, int iC)
+  : BlockLatticeF3D<T,DESCRIPTOR>(blockLattice, 1), _iC(iC)
+{
+  this->getName() = "cuboid";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeCuboid3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  output[0] = _iC + 1;
+  return true;
+}
+
+
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticePhysPressure3D<T, DESCRIPTOR>::BlockLatticePhysPressure3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice,
+  int overlap,
+  const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF3D<T, DESCRIPTOR>(blockLattice, converter, 1),
+    _overlap(overlap)
+{
+  this->getName() = "physPressure";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticePhysPressure3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  // lattice pressure = c_s^2 ( rho -1 )
+  T latticePressure = ( this->_blockLattice.get(input[0]+_overlap, input[1]+_overlap, input[2]+_overlap).computeRho() - 1.0) / descriptors::invCs2<T,DESCRIPTOR>();
+  output[0] = this->_converter.getPhysPressure(latticePressure);
+
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticePhysVelocity3D<T, DESCRIPTOR>::BlockLatticePhysVelocity3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice,
+  int overlap,
+  const UnitConverter<T,DESCRIPTOR>& converter,
+  bool print)
+  : BlockLatticePhysF3D<T, DESCRIPTOR>(blockLattice, converter, 3),
+    _overlap(overlap),
+    _print(print)
+{
+  this->getName() = "physVelocity";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticePhysVelocity3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  if (_print) {
+    std::cout << input[0] << " " << input[1] << " " << input[2] << " | "
+              << singleton::mpi().getRank() << std::endl;
+  }
+
+  T rho;
+  this->_blockLattice.get(
+    input[0]+_overlap, input[1]+_overlap, input[2]+_overlap).computeRhoU(rho, output);
+  output[0] = this->_converter.getPhysVelocity(output[0]);
+  output[1] = this->_converter.getPhysVelocity(output[1]);
+  output[2] = this->_converter.getPhysVelocity(output[2]);
+
+  return true;
+}
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysExternalPorosity3D<T,DESCRIPTOR>::BlockLatticePhysExternalPorosity3D(
+  BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice,
+  int overlap,
+  const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF3D<T,DESCRIPTOR>(blockLattice,converter,2),
+    _overlap(overlap)
+{
+  this->getName() = "ExtPorosityField";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysExternalPorosity3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  this->_blockLattice.get(
+    input[0]+_overlap, input[1]+_overlap, input[2]+_overlap
+  ).template computeField<descriptors::POROSITY>(output);
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticePhysExternalVelocity3D<T, DESCRIPTOR>::BlockLatticePhysExternalVelocity3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF3D<T, DESCRIPTOR>(blockLattice, converter, 3)
+{
+  this->getName() = "physVelExtField";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticePhysExternalVelocity3D<T, DESCRIPTOR>::operator()(
+  T output[], const int input[])
+{
+  this->_blockLattice.get(input[0], input[1], input[2]).template computeField<descriptors::VELOCITY>(output);
+  output[0] = this->_converter.getPhysVelocity(output[0]);
+  output[1] = this->_converter.getPhysVelocity(output[1]);
+  output[2] = this->_converter.getPhysVelocity(output[2]);
+  return true;
+}
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysExternalParticleVelocity3D<T,DESCRIPTOR>::BlockLatticePhysExternalParticleVelocity3D
+(BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticePhysF3D<T,DESCRIPTOR>(blockLattice,converter,2)
+{
+  this->getName() = "ExtParticleVelocityField";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysExternalParticleVelocity3D<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]);
+    output[2]=this->_converter.getPhysVelocity(velocity_numerator[2]/velocity_denominator[0]);
+    return true;
+  }
+  output[0]=this->_converter.getPhysVelocity(velocity_numerator[0]);
+  output[1]=this->_converter.getPhysVelocity(velocity_numerator[1]);
+  output[2]=this->_converter.getPhysVelocity(velocity_numerator[2]);
+  return true;
+}
+
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticePhysExternal3D<T, DESCRIPTOR>::BlockLatticePhysExternal3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice,
+  T convFactorToPhysUnits, int offset, int size)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 3),
+    _convFactorToPhysUnits(convFactorToPhysUnits),
+    _offset(offset), _size(size)
+{
+  this->getName() = "physExtField";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticePhysExternal3D<T, DESCRIPTOR>::operator()(
+  T output[], const int input[])
+{
+  const auto& cell = this->_blockLattice.get(input[0], input[1], input[2]);
+  output[0] = _convFactorToPhysUnits*cell[_offset+0];
+  output[1] = _convFactorToPhysUnits*cell[_offset+1];
+  output[2] = _convFactorToPhysUnits*cell[_offset+2];
+  return true;
+}
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysBoundaryForce3D<T,DESCRIPTOR>::BlockLatticePhysBoundaryForce3D(
+  BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice,
+  BlockIndicatorF3D<T>&                  indicatorF,
+  const UnitConverter<T,DESCRIPTOR>&     converter)
+  : BlockLatticePhysF3D<T,DESCRIPTOR>(blockLattice, converter, 3),
+    _indicatorF(indicatorF),
+    _blockGeometry(indicatorF.getBlockGeometryStructure())
+{
+  this->getName() = "physBoundaryForce";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticePhysBoundaryForce3D<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
+      const Vector<int,3> c = descriptors::c<DESCRIPTOR>(iPop);
+      // Get next cell located in the current direction
+      // Check if the next cell is a fluid node
+      if (_blockGeometry.get(input[0] + c[0], input[1] + c[1], input[2] + c[2]) == 1) {
+        // Get f_q of next fluid cell where l = opposite(q)
+        T f = this->_blockLattice.get(input[0] + c[0], input[1] + c[1], input[2] + c[2])[iPop];
+        // Get f_l of the boundary cell
+        // Add f_q and f_opp
+        f += this->_blockLattice.get(input)[util::opposite<DESCRIPTOR>(iPop)];
+        // Update force
+        for (int i = 0; i < this->getTargetDim(); ++i) {
+          output[i] -= c[i] * f;
+        }
+      }
+    }
+    for (int i = 0; i < this->getTargetDim(); ++i) {
+      output[i] = this->_converter.getPhysForce(output[i]);
+    }
+  }
+  return true;
+}
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePSMPhysForce3D<T,DESCRIPTOR>::BlockLatticePSMPhysForce3D(
+  BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice,
+  const UnitConverter<T,DESCRIPTOR>&     converter,
+  int mode_)
+  : BlockLatticePhysF3D<T,DESCRIPTOR>(blockLattice, converter, 3)
+{
+  this->getName() = "physPSMForce";
+  mode = (Mode) mode_;
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticePSMPhysForce3D<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).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).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).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], input[2])[iPop])
+                  + (1 - omega)
+                  * (this->_blockLattice.get(input[0], input[1], input[2])[iPop]
+                     - lbDynamicsHelpers<T, DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr));
+        break;
+      case M3:
+        omega_s =
+          (this->_blockLattice.get(input[0], input[1], input[2])[descriptors::opposite<DESCRIPTOR>(iPop)]
+           - lbDynamicsHelpers<T, DESCRIPTOR>::equilibrium(
+             descriptors::opposite<DESCRIPTOR>(iPop), rho, u_s, uSqr_s))
+          - (this->_blockLattice.get(input[0], input[1], input[2])[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>
+BlockLatticePhysWallShearStress3D<T,DESCRIPTOR>::BlockLatticePhysWallShearStress3D(
+    BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice,
+    BlockGeometryStructure3D<T>& blockGeometry,
+    int overlap,
+    int material,
+    const UnitConverter<T,DESCRIPTOR>& converter,
+    IndicatorF3D<T>& indicator)
+    : BlockLatticePhysF3D<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(4, 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);
+
+      for (int iZ = 1; iZ < _blockGeometry.getNz() - 1; iZ++) {
+        _discreteNormal[iX-1][iY-1].resize(_blockGeometry.getNz() - 2);
+        _normal[iX-1][iY-1].resize(_blockGeometry.getNz() - 2);
+
+        if (_blockGeometry.get(iX, iY, iZ) == _material) {
+          discreteNormalOutwards = _blockGeometry.getStatistics().getType(iX, iY, iZ);
+          _discreteNormal[iX - 1][iY - 1][iZ - 1].resize(3);
+          _normal[iX - 1][iY - 1][iZ - 1].resize(3);
+
+          _discreteNormal[iX - 1][iY- 1][iZ- 1][0] = -discreteNormalOutwards[1];
+          _discreteNormal[iX- 1][iY- 1][iZ- 1][1] = -discreteNormalOutwards[2];
+          _discreteNormal[iX- 1][iY- 1][iZ- 1][2] = -discreteNormalOutwards[3];
+
+          T physR[3];
+          _blockGeometry.getPhysR(physR,iX, iY, iZ);
+          Vector<T,3> origin(physR[0],physR[1],physR[2]);
+          Vector<T,3> direction(-_discreteNormal[iX- 1][iY- 1][iZ- 1][0] * scaling,
+                                -_discreteNormal[iX- 1][iY- 1][iZ- 1][1] * scaling,
+                                -_discreteNormal[iX- 1][iY- 1][iZ- 1][2] * scaling);
+          Vector<T,3> normal(0.,0.,0.);
+          origin[0] = physR[0];
+          origin[1] = physR[1];
+          origin[2] = physR[2];
+
+          indicator.normal(normal, origin, direction);
+          normal.normalize();
+
+          _normal[iX- 1][iY- 1][iZ- 1][0] = normal[0];
+          _normal[iX- 1][iY- 1][iZ- 1][1] = normal[1];
+          _normal[iX- 1][iY- 1][iZ- 1][2] = normal[2];
+        }
+      }
+    }
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticePhysWallShearStress3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  output[0] = T();
+
+  if (input[0] + _overlap < 1 ||
+      input[1] + _overlap < 1 ||
+      input[2] + _overlap < 1 ||
+      input[0] + _overlap >= _blockGeometry.getNx()-1 ||
+      input[1] + _overlap >= _blockGeometry.getNy()-1 ||
+      input[2] + _overlap >= _blockGeometry.getNz()-1 ){
+
+#ifdef OLB_DEBUG
+    std::cout << "Input address not mapped by _discreteNormal, overlap too small" << std::endl;
+#endif
+    return true;
+  }
+
+  if (_blockGeometry.get(input[0]+_overlap,input[1]+_overlap,input[2]+_overlap) == _material) {
+
+    T traction[3];
+    T stress[6];
+    T rho = this->_blockLattice.get(input[0] + _overlap + _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0],
+                                    input[1] + _overlap + _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1],
+                                    input[2] + _overlap + _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2]).computeRho();
+     this->_blockLattice.get(input[0] + _overlap +   _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0],
+                      input[1] + _overlap +   _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1],
+                      input[2] + _overlap +   _discreteNormal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2]).computeStress(stress);
+
+    traction[0] = stress[0]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0] +
+                  stress[1]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1] +
+                  stress[2]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2];
+    traction[1] = stress[1]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0] +
+                  stress[3]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1] +
+                  stress[4]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2];
+    traction[2] = stress[2]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0] +
+                  stress[4]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1] +
+                  stress[5]*_physFactor/rho*_normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2];
+
+    T traction_normal_SP;
+    T tractionNormalComponent[3];
+    // scalar product of traction and normal vector
+    traction_normal_SP = traction[0] * _normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0] +
+                         traction[1] * _normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1] +
+                         traction[2] * _normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2];
+    tractionNormalComponent[0] = traction_normal_SP * _normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][0];
+    tractionNormalComponent[1] = traction_normal_SP * _normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][1];
+    tractionNormalComponent[2] = traction_normal_SP * _normal[input[0]+_overlap-1][input[1]+_overlap-1][input[2]+_overlap-1][2];
+
+    T WSS[3];
+    WSS[0] = traction[0] - tractionNormalComponent[0];
+    WSS[1] = traction[1] - tractionNormalComponent[1];
+    WSS[2] = traction[2] - tractionNormalComponent[2];
+    // magnitude of the wall shear stress vector
+    output[0] = sqrt(WSS[0]*WSS[0] + WSS[1]*WSS[1] + WSS[2]*WSS[2]);
+
+    return true;
+  }
+  else {
+    return true;
+  }
+}
+
+
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysCorrBoundaryForce3D<T,DESCRIPTOR>::BlockLatticePhysCorrBoundaryForce3D(
+  BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice,
+  BlockIndicatorF3D<T>&                  indicatorF,
+  const UnitConverter<T,DESCRIPTOR>&     converter)
+  : BlockLatticePhysF3D<T,DESCRIPTOR>(blockLattice, converter, 3),
+    _indicatorF(indicatorF),
+    _blockGeometry(indicatorF.getBlockGeometryStructure())
+{
+  this->getName() = "physCorrBoundaryForce";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticePhysCorrBoundaryForce3D<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
+      const Vector<int,3> c = descriptors::c<DESCRIPTOR>(iPop);
+      // Get next cell located in the current direction
+      // Check if the next cell is a fluid node
+      if (_blockGeometry.get(input[0] + c[0], input[1] + c[1], input[2] + c[2]) == 1) {
+        // Get f_q of next fluid cell where l = opposite(q)
+        T f = this->_blockLattice.get(input[0] + c[0], input[1] + c[1], input[2] + c[2])[iPop];
+        // Get f_l of the boundary cell
+        // Add f_q and f_opp
+        f += this->_blockLattice.get(input)[util::opposite<DESCRIPTOR>(iPop)];
+        // Update force
+        for (int i = 0; i < this->getTargetDim(); ++i) {
+          output[i] -= c[i] * (f - 2. * descriptors::t<T,DESCRIPTOR>(iPop));
+        }
+      }
+    }
+    for (int i = 0; i < this->getTargetDim(); ++i) {
+      output[i] = this->_converter.getPhysForce(output[i]);
+    }
+  }
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR, typename FIELD>
+BlockLatticeField3D<T,DESCRIPTOR,FIELD>::BlockLatticeField3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, DESCRIPTOR::template size<FIELD>())
+{
+  this->getName() = "externalField";
+}
+
+template<typename T, typename DESCRIPTOR, typename FIELD>
+bool BlockLatticeField3D<T,DESCRIPTOR,FIELD>::operator()(T output[], const int input[])
+{
+  this->_blockLattice.get(input[0], input[1], input[2]).template computeField<FIELD>(output);
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticePorosity3D<T, DESCRIPTOR>::BlockLatticePorosity3D(BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 1)
+{
+  this->getName() = "porosity";
+}
+
+// under construction
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticePorosity3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+#ifndef excludeDualDynamics
+  this->_blockLattice.get(input[0], input[1], input[2]).template computeField<descriptors::POROSITY>(
+    output);
+#else
+  this->_blockLattice.get(input[0], input[1], input[2]).template computeField<descriptors::POROSITY>(
+    output);
+
+#endif
+  return true;
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeVolumeFractionApproximation3D<T, DESCRIPTOR>::BlockLatticeVolumeFractionApproximation3D(BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice,
+    BlockGeometryStructure3D<T>& blockGeometry,
+    IndicatorF3D<T>& indicator,
+    int refinementLevel,
+    const UnitConverter<T,DESCRIPTOR>& converter,
+    bool insideOut)
+  : BlockLatticeF3D<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 * _refinementLevel))
+{
+  this->getName() = "volumeFractionApproximation";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeVolumeFractionApproximation3D<T, DESCRIPTOR>::operator()(T output[], const int input[])
+{
+  output[0] = 0.;
+  T physR[3];
+  bool inside[1];
+  _blockGeometry.getPhysR(physR, input[0], input[1], input[2]);
+
+  T subGridMinPhysR[3];
+  subGridMinPhysR[0] = physR[0] + _physSubGridMinPhysRshift;
+  subGridMinPhysR[1] = physR[1] + _physSubGridMinPhysRshift;
+  subGridMinPhysR[2] = physR[2] + _physSubGridMinPhysRshift;
+
+  for (int x = 0; x < _refinementLevel; x++) {
+    for (int y = 0; y < _refinementLevel; y++) {
+      for (int z = 0; z < _refinementLevel; z++) {
+        physR[0] = subGridMinPhysR[0] + x * _physSubGridDeltaX;
+