From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Mon, 24 Jun 2019 14:43:36 +0200
Subject: Initialize at openlb-1-3
---
src/functors/lattice/blockLatticeLocalF3D.hh | 1556 ++++++++++++++++++++++++++
1 file changed, 1556 insertions(+)
create mode 100644 src/functors/lattice/blockLatticeLocalF3D.hh
(limited to 'src/functors/lattice/blockLatticeLocalF3D.hh')
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
+ *
+ *
+ * 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
+#include
+
+#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
+BlockLatticeFpop3D::BlockLatticeFpop3D(BlockLatticeStructure3D& blockLattice)
+ : BlockLatticeF3D(blockLattice, DESCRIPTOR::q)
+{
+ this->getName() = "fPop";
+}
+
+template
+bool BlockLatticeFpop3D::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(iPop);
+ }
+ return true;
+}
+
+template
+BlockLatticeDissipation3D::BlockLatticeDissipation3D(
+ BlockLatticeStructure3D& blockLattice, const UnitConverter& converter)
+ : BlockLatticeF3D(blockLattice, 1), _converter(converter)
+{
+ this->getName() = "dissipation";
+}
+
+template
+bool BlockLatticeDissipation3D::operator()(T output[], const int input[])
+{
+ T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal::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::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(), 2) / rho / 2.;
+
+ return true;
+}
+
+template
+BlockLatticePhysDissipation3D::BlockLatticePhysDissipation3D(
+ BlockLatticeStructure3D& blockLattice,
+ int overlap,
+ const UnitConverter& converter)
+ : BlockLatticeF3D(blockLattice, 1),
+ _overlap(overlap),
+ _converter(converter)
+{
+ this->getName() = "physDissipation";
+}
+
+template
+bool BlockLatticePhysDissipation3D::operator()(T output[], const int input[])
+{
+ T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal::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::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() / rho, 2) / 2.
+ * this->_converter.getPhysViscosity() / nuLattice / dt / dt;
+
+ return true;
+}
+
+template
+BlockLatticeEffevtiveDissipation3D::BlockLatticeEffevtiveDissipation3D(
+ BlockLatticeStructure3D& blockLattice, const UnitConverter& converter, T smagoConst,
+ LESDynamics& LESdynamics)
+ : BlockLatticeF3D(blockLattice, 1),
+ _converter(converter), _smagoConst(smagoConst), _LESdynamics(LESdynamics)
+{
+ this->getName() = "EffevtiveDissipation";
+}
+
+template
+bool BlockLatticeEffevtiveDissipation3D::operator()(T output[], const int input[])
+{
+ T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal::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::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();
+
+ output[0] = PiNeqNormSqr * (nuEff)
+ * pow(omegaEff * descriptors::invCs2() / rho, 2) / 2.;
+
+ return true;
+}
+
+template
+BlockLatticePhysEffevtiveDissipation3D::BlockLatticePhysEffevtiveDissipation3D(
+ BlockLatticeStructure3D& blockLattice, const UnitConverter& converter, T smagoConst,
+ LESDynamics& LESdynamics)
+ : BlockLatticeF3D(blockLattice, 1),
+ _converter(converter), _smagoConst(smagoConst), _LESdynamics(LESdynamics)
+{
+ this->getName() = "physEffevtiveDissipation";
+}
+
+template
+bool BlockLatticePhysEffevtiveDissipation3D::operator()(T output[], const int input[])
+{
+ T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal::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::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(); // BGK shear viscosity definition
+
+ output[0] = PiNeqNormSqr * nuEff
+ * pow(omegaEff * descriptors::invCs2() / rho, 2) / 2.
+ * _converter.getPhysViscosity() / _converter.getLatticeViscosity() / dt / dt;
+
+ return true;
+}
+
+template
+BlockLatticeDensity3D::BlockLatticeDensity3D(
+ BlockLatticeStructure3D& blockLattice)
+ : BlockLatticeF3D(blockLattice, 1)
+{
+ this->getName() = "density";
+}
+
+template
+bool BlockLatticeDensity3D::operator()(T output[], const int input[])
+{
+ output[0] = this->_blockLattice.get(input[0], input[1], input[2]).computeRho();
+ return true;
+}
+
+template
+BlockLatticeVelocity3D::BlockLatticeVelocity3D(
+ BlockLatticeStructure3D& blockLattice)
+ : BlockLatticeF3D(blockLattice, 3)
+{
+ this->getName() = "velocity";
+}
+
+template
+bool BlockLatticeVelocity3D::operator()(T output[], const int input[])
+{
+ T rho;
+ this->_blockLattice.get(input[0], input[1], input[2]).computeRhoU(rho, output);
+ return true;
+}
+
+template
+BlockLatticeExternalVelocity3D::BlockLatticeExternalVelocity3D(
+ BlockLatticeStructure3D& blockLattice)
+ : BlockLatticeF3D(blockLattice, 3)
+{
+ this->getName() = "externalVelocity";
+}
+
+template
+bool BlockLatticeExternalVelocity3D::operator()(T output[], const int input[])
+{
+ T* ExtVel = this->_blockLattice.get(input[0], input[1], input[2]).template getFieldPointer();
+ for (int iVel=0; iVel
+BlockLatticeFlux3D::BlockLatticeFlux3D(
+ BlockLatticeStructure3D& blockLattice)
+ : BlockLatticeF3D(blockLattice, 3)
+{
+ this->getName() = "flux";
+}
+
+template
+bool BlockLatticeFlux3D::operator()(T output[], const int input[])
+{
+ this->_blockLattice.get(input[0], input[1], input[2]).computeJ(output);
+ return true;
+}
+
+template
+BlockLatticeStrainRate3D::BlockLatticeStrainRate3D(
+ BlockLatticeStructure3D& blockLattice, const UnitConverter& converter)
+ : BlockLatticePhysF3D(blockLattice, converter, 9)
+{
+ this->getName() = "strainRate";
+}
+
+template
+bool BlockLatticeStrainRate3D::operator()(T output[], const int input[])
+{
+ T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal::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() / rho / 2.;
+ output[1] = -pi[1] * omega * descriptors::invCs2() / rho / 2.;
+ output[2] = -pi[2] * omega * descriptors::invCs2() / rho / 2.;
+ output[3] = -pi[1] * omega * descriptors::invCs2() / rho / 2.;
+ output[4] = -pi[3] * omega * descriptors::invCs2() / rho / 2.;
+ output[5] = -pi[4] * omega * descriptors::invCs2() / rho / 2.;
+ output[6] = -pi[2] * omega * descriptors::invCs2() / rho / 2.;
+ output[7] = -pi[4] * omega * descriptors::invCs2() / rho / 2.;
+ output[8] = -pi[5] * omega * descriptors::invCs2() / rho / 2.;
+
+ return true;
+}
+
+template
+BlockLatticePhysStrainRate3D::BlockLatticePhysStrainRate3D(
+ BlockLatticeStructure3D& blockLattice,
+ int overlap,
+ const UnitConverter& converter)
+ : BlockLatticePhysF3D(blockLattice, converter, 9),
+ _overlap(overlap)
+{
+ this->getName() = "strainRate";
+}
+
+template
+bool BlockLatticePhysStrainRate3D::operator()(T output[], const int input[])
+{
+ T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal::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() / rho / 2. / dt;
+ output[1] = -pi[1] * omega * descriptors::invCs2() / rho / 2. / dt;
+ output[2] = -pi[2] * omega * descriptors::invCs2() / rho / 2. / dt;
+ output[3] = -pi[1] * omega * descriptors::invCs2() / rho / 2. / dt;
+ output[4] = -pi[3] * omega * descriptors::invCs2() / rho / 2. / dt;
+ output[5] = -pi[4] * omega * descriptors::invCs2() / rho / 2. / dt;
+ output[6] = -pi[2] * omega * descriptors::invCs2() / rho / 2. / dt;
+ output[7] = -pi[4] * omega * descriptors::invCs2() / rho / 2. / dt;
+ output[8] = -pi[5] * omega * descriptors::invCs2() / rho / 2. / dt;
+
+ return true;
+}
+
+template
+BlockLatticeGeometry3D::BlockLatticeGeometry3D(BlockLatticeStructure3D& blockLattice,
+ BlockGeometryStructure3D& blockGeometry, int material)
+ : BlockLatticeF3D(blockLattice, 1),
+ _blockGeometry(blockGeometry),
+ _material(material)
+{
+ this->getName() = "geometry";
+}
+
+template
+bool BlockLatticeGeometry3D::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
+BlockLatticeRank3D::BlockLatticeRank3D(
+ BlockLatticeStructure3D& blockLattice)
+ : BlockLatticeF3D(blockLattice, 1)
+{
+ this->getName() = "rank";
+}
+
+template
+bool BlockLatticeRank3D::operator()(T output[], const int input[])
+{
+ output[0] = singleton::mpi().getRank() + 1;
+ return true;
+}
+
+
+template
+BlockLatticeCuboid3D::BlockLatticeCuboid3D(
+ BlockLatticeStructure3D& blockLattice, int iC)
+ : BlockLatticeF3D(blockLattice, 1), _iC(iC)
+{
+ this->getName() = "cuboid";
+}
+
+template
+bool BlockLatticeCuboid3D::operator()(T output[], const int input[])
+{
+ output[0] = _iC + 1;
+ return true;
+}
+
+
+
+template
+BlockLatticePhysPressure3D::BlockLatticePhysPressure3D(
+ BlockLatticeStructure3D& blockLattice,
+ int overlap,
+ const UnitConverter& converter)
+ : BlockLatticePhysF3D(blockLattice, converter, 1),
+ _overlap(overlap)
+{
+ this->getName() = "physPressure";
+}
+
+template
+bool BlockLatticePhysPressure3D::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();
+ output[0] = this->_converter.getPhysPressure(latticePressure);
+
+ return true;
+}
+
+template
+BlockLatticePhysVelocity3D::BlockLatticePhysVelocity3D(
+ BlockLatticeStructure3D& blockLattice,
+ int overlap,
+ const UnitConverter& converter,
+ bool print)
+ : BlockLatticePhysF3D(blockLattice, converter, 3),
+ _overlap(overlap),
+ _print(print)
+{
+ this->getName() = "physVelocity";
+}
+
+template
+bool BlockLatticePhysVelocity3D::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
+BlockLatticePhysExternalPorosity3D::BlockLatticePhysExternalPorosity3D(
+ BlockLatticeStructure3D& blockLattice,
+ int overlap,
+ const UnitConverter& converter)
+ : BlockLatticePhysF3D(blockLattice,converter,2),
+ _overlap(overlap)
+{
+ this->getName() = "ExtPorosityField";
+}
+
+template
+bool BlockLatticePhysExternalPorosity3D::operator() (T output[], const int input[])
+{
+ this->_blockLattice.get(
+ input[0]+_overlap, input[1]+_overlap, input[2]+_overlap
+ ).template computeField(output);
+ return true;
+}
+
+template
+BlockLatticePhysExternalVelocity3D::BlockLatticePhysExternalVelocity3D(
+ BlockLatticeStructure3D& blockLattice, const UnitConverter& converter)
+ : BlockLatticePhysF3D(blockLattice, converter, 3)
+{
+ this->getName() = "physVelExtField";
+}
+
+template
+bool BlockLatticePhysExternalVelocity3D::operator()(
+ T output[], const int input[])
+{
+ this->_blockLattice.get(input[0], input[1], input[2]).template computeField(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
+BlockLatticePhysExternalParticleVelocity3D::BlockLatticePhysExternalParticleVelocity3D
+(BlockLatticeStructure3D& blockLattice, const UnitConverter& converter)
+ : BlockLatticePhysF3D(blockLattice,converter,2)
+{
+ this->getName() = "ExtParticleVelocityField";
+}
+
+template
+bool BlockLatticePhysExternalParticleVelocity3D::operator() (T output[], const int input[])
+{
+ const T* velocity_numerator = this->blockLattice.get(input).template getFieldPointer();
+ const T* velocity_denominator = this->blockLattice.get(input).template getFieldPointer();
+
+ if (velocity_denominator[0] > std::numeric_limits::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
+BlockLatticePhysExternal3D::BlockLatticePhysExternal3D(
+ BlockLatticeStructure3D& blockLattice,
+ T convFactorToPhysUnits, int offset, int size)
+ : BlockLatticeF3D(blockLattice, 3),
+ _convFactorToPhysUnits(convFactorToPhysUnits),
+ _offset(offset), _size(size)
+{
+ this->getName() = "physExtField";
+}
+
+template
+bool BlockLatticePhysExternal3D::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
+BlockLatticePhysBoundaryForce3D::BlockLatticePhysBoundaryForce3D(
+ BlockLatticeStructure3D& blockLattice,
+ BlockIndicatorF3D& indicatorF,
+ const UnitConverter& converter)
+ : BlockLatticePhysF3D(blockLattice, converter, 3),
+ _indicatorF(indicatorF),
+ _blockGeometry(indicatorF.getBlockGeometryStructure())
+{
+ this->getName() = "physBoundaryForce";
+}
+
+template
+bool BlockLatticePhysBoundaryForce3D::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 c = descriptors::c(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(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
+BlockLatticePSMPhysForce3D::BlockLatticePSMPhysForce3D(
+ BlockLatticeStructure3D& blockLattice,
+ const UnitConverter& converter,
+ int mode_)
+ : BlockLatticePhysF3D(blockLattice, converter, 3)
+{
+ this->getName() = "physPSMForce";
+ mode = (Mode) mode_;
+}
+
+template
+bool BlockLatticePSMPhysForce3D::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());
+
+ //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())[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(u_s);
+ T uSqr = util::normSqr(u);
+ for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
+ switch (mode) {
+ case M2:
+ omega_s = (lbDynamicsHelpers::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::equilibrium(iPop, rho, u, uSqr));
+ break;
+ case M3:
+ omega_s =
+ (this->_blockLattice.get(input[0], input[1], input[2])[descriptors::opposite(iPop)]
+ - lbDynamicsHelpers::equilibrium(
+ descriptors::opposite(iPop), rho, u_s, uSqr_s))
+ - (this->_blockLattice.get(input[0], input[1], input[2])[iPop]
+ - lbDynamicsHelpers::equilibrium(iPop, rho, u_s, uSqr_s));
+ }
+
+ for (int i = 0; i < this->getTargetDim(); ++i) {
+ output[i] -= descriptors::c(iPop,i) * omega_s;
+ }
+ }
+
+ for (int i = 0; i < this->getTargetDim(); ++i) {
+ output[i] = this->_converter.getPhysForce(output[i] * paramB);
+ }
+ }
+ return true;
+}
+
+template
+BlockLatticePhysWallShearStress3D::BlockLatticePhysWallShearStress3D(
+ BlockLatticeStructure3D& blockLattice,
+ BlockGeometryStructure3D& blockGeometry,
+ int overlap,
+ int material,
+ const UnitConverter& converter,
+ IndicatorF3D& indicator)
+ : BlockLatticePhysF3D(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() / dt * this->_converter.getPhysDensity() * this->_converter.getPhysViscosity();
+ std::vector 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 origin(physR[0],physR[1],physR[2]);
+ Vector 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 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
+bool BlockLatticePhysWallShearStress3D::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
+BlockLatticePhysCorrBoundaryForce3D::BlockLatticePhysCorrBoundaryForce3D(
+ BlockLatticeStructure3D& blockLattice,
+ BlockIndicatorF3D& indicatorF,
+ const UnitConverter& converter)
+ : BlockLatticePhysF3D(blockLattice, converter, 3),
+ _indicatorF(indicatorF),
+ _blockGeometry(indicatorF.getBlockGeometryStructure())
+{
+ this->getName() = "physCorrBoundaryForce";
+}
+
+template
+bool BlockLatticePhysCorrBoundaryForce3D::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 c = descriptors::c(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(iPop)];
+ // Update force
+ for (int i = 0; i < this->getTargetDim(); ++i) {
+ output[i] -= c[i] * (f - 2. * descriptors::t(iPop));
+ }
+ }
+ }
+ for (int i = 0; i < this->getTargetDim(); ++i) {
+ output[i] = this->_converter.getPhysForce(output[i]);
+ }
+ }
+ return true;
+}
+
+template
+BlockLatticeField3D::BlockLatticeField3D(
+ BlockLatticeStructure3D& blockLattice)
+ : BlockLatticeF3D(blockLattice, DESCRIPTOR::template size())
+{
+ this->getName() = "externalField";
+}
+
+template
+bool BlockLatticeField3D::operator()(T output[], const int input[])
+{
+ this->_blockLattice.get(input[0], input[1], input[2]).template computeField(output);
+ return true;
+}
+
+template
+BlockLatticePorosity3D::BlockLatticePorosity3D(BlockLatticeStructure3D& blockLattice)
+ : BlockLatticeF3D(blockLattice, 1)
+{
+ this->getName() = "porosity";
+}
+
+// under construction
+template
+bool BlockLatticePorosity3D::operator()(T output[], const int input[])
+{
+#ifndef excludeDualDynamics
+ this->_blockLattice.get(input[0], input[1], input[2]).template computeField(
+ output);
+#else
+ this->_blockLattice.get(input[0], input[1], input[2]).template computeField(
+ output);
+
+#endif
+ return true;
+}
+
+template
+BlockLatticeVolumeFractionApproximation3D::BlockLatticeVolumeFractionApproximation3D(BlockLatticeStructure3D& blockLattice,
+ BlockGeometryStructure3D& blockGeometry,
+ IndicatorF3D& indicator,
+ int refinementLevel,
+ const UnitConverter