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/blockLatticeLocalF2D.hh | 1229 ++++++++++++++++++++++++++
1 file changed, 1229 insertions(+)
create mode 100644 src/functors/lattice/blockLatticeLocalF2D.hh
(limited to 'src/functors/lattice/blockLatticeLocalF2D.hh')
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
+ *
+ *
+ * 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
+#include
+
+#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
+BlockLatticeDissipation2D::BlockLatticeDissipation2D
+(BlockLatticeStructure2D& blockLattice, const UnitConverter& converter)
+ : BlockLatticeF2D(blockLattice,1), _converter(converter)
+{
+ this->getName() = "dissipation";
+}
+
+// todo: get functor working
+template
+bool BlockLatticeDissipation2D::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::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::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(),2)/rho/2.;
+
+ // return std::vector(1,finalResult);
+ // } else {
+ // return std::vector(); // empty vector
+ // }
+
+ return false;
+}
+
+template
+BlockLatticePhysDissipation2D::BlockLatticePhysDissipation2D
+(BlockLatticeStructure2D& blockLattice, const UnitConverter& converter)
+ : BlockLatticePhysF2D(blockLattice,converter,1)
+{
+ this->getName() = "physDissipation";
+}
+
+template
+bool BlockLatticePhysDissipation2D::operator() (T output[], const int input[])
+{
+ T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal::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::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()/this->_converter.getLatticeViscosity()/dt/dt;
+
+ return true;
+}
+
+template
+BlockLatticeDensity2D::BlockLatticeDensity2D
+(BlockLatticeStructure2D& blockLattice)
+ : BlockLatticeF2D(blockLattice,1)
+{
+ this->getName() = "density";
+}
+
+
+template
+bool BlockLatticeDensity2D::operator() (T output[], const int input[])
+{
+ output[0] = this->_blockLattice.get( input[0], input[1] ).computeRho();
+ return true;
+}
+
+
+template
+BlockLatticeVelocity2D::BlockLatticeVelocity2D
+(BlockLatticeStructure2D& blockLattice)
+ : BlockLatticeF2D(blockLattice,2)
+{
+ this->getName() = "velocity";
+}
+
+template
+bool BlockLatticeVelocity2D::operator() (T output[], const int input[])
+{
+ T rho;
+ this->_blockLattice.get(input[0], input[1]).computeRhoU(rho,output);
+ return true;
+}
+
+/*template
+BlockLatticeStrainRate2D::BlockLatticeStrainRate2D
+ (BlockLatticeStructure2D& blockLattice, const UnitConverter& converter)
+ : BlockLatticePhysF2D(blockLattice,converter,4)
+{ this->getName() = "strainRate"; }
+
+template
+std::vector BlockLatticeStrainRate2D::operator() (std::vector input) {
+
+ T strainRate[4];
+ T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal::n];
+ this->_blockLattice.get( input[0], input[1] ).computeAllMomenta(rho, uTemp, pi);
+
+ T omega = this->_converter.getOmega();
+
+ strainRate[0] = -pi[0]*omega*descriptors::invCs2()/rho/2.;
+ strainRate[1] = -pi[1]*omega*descriptors::invCs2()/rho/2.;
+ strainRate[2] = -pi[1]*omega*descriptors::invCs2()/rho/2.;
+ strainRate[3] = -pi[2]*omega*descriptors::invCs2()/rho/2.;
+
+ //cout << pi[0] << " " << pi[1] << " " << pi[2] << " " << descriptors::invCs2() << endl;
+
+ std::vector output(strainRate, strainRate+4); // first adress, last adress
+ return output;
+}*/
+
+template
+BlockLatticePhysStrainRate2D::BlockLatticePhysStrainRate2D
+(BlockLatticeStructure2D& blockLattice, const UnitConverter& converter)
+ : BlockLatticePhysF2D(blockLattice,converter,4)
+{
+ this->getName() = "strainRate";
+}
+
+template
+bool BlockLatticePhysStrainRate2D::operator() (T output[], const int input[])
+{
+ T rho, uTemp[DESCRIPTOR::d], pi[util::TensorVal::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()/rho/2./dt;
+ output[1] = -pi[1]*omega*descriptors::invCs2()/rho/2./dt;
+ output[2] = -pi[1]*omega*descriptors::invCs2()/rho/2./dt;
+ output[3] = -pi[2]*omega*descriptors::invCs2()/rho/2./dt;
+
+ return true;
+}
+
+template
+BlockLatticeGeometry2D::BlockLatticeGeometry2D
+(BlockLatticeStructure2D& blockLattice, BlockGeometryStructure2D& blockGeometry, int material)
+ : BlockLatticeF2D(blockLattice,1), _blockGeometry(blockGeometry), _material(material)
+{
+ this->getName() = "geometry";
+}
+
+template
+bool BlockLatticeGeometry2D::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
+BlockLatticeRank2D::BlockLatticeRank2D
+(BlockLatticeStructure2D& blockLattice)
+ : BlockLatticeF2D(blockLattice,1)
+{
+ this->getName() = "rank";
+}
+
+template
+bool BlockLatticeRank2D::operator() (T output[], const int input[])
+{
+ output[0] = singleton::mpi().getRank() + 1;
+ return false;
+}
+
+
+template
+BlockLatticeCuboid2D::BlockLatticeCuboid2D
+(BlockLatticeStructure2D& blockLattice, const int iC)
+ : BlockLatticeF2D(blockLattice,1), _iC(iC)
+{
+ this->getName() = "cuboid";
+}
+
+template
+bool BlockLatticeCuboid2D::operator() (T output[], const int input[])
+{
+ output[0] = _iC + 1;
+ return false;
+}
+
+
+template
+BlockLatticePhysPressure2D::BlockLatticePhysPressure2D
+(BlockLatticeStructure2D& blockLattice, const UnitConverter& converter)
+ : BlockLatticePhysF2D(blockLattice,converter,1)
+{
+ this->getName() = "physPressure";
+}
+
+
+template
+bool BlockLatticePhysPressure2D::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();
+ output[0] = this->_converter.getPhysPressure(latticePressure);
+
+ return true;
+}
+
+
+template
+BlockLatticePhysVelocity2D::BlockLatticePhysVelocity2D
+(BlockLatticeStructure2D& blockLattice, const UnitConverter& converter)
+ : BlockLatticePhysF2D(blockLattice,converter,2)
+{
+ this->getName() = "physVelocity";
+}
+
+template
+bool BlockLatticePhysVelocity2D::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
+BlockLatticeField2D::BlockLatticeField2D(
+ BlockLatticeStructure2D& blockLattice)
+ : BlockLatticeF2D(blockLattice, DESCRIPTOR::template size())
+{
+ this->getName() = "extField";
+}
+
+template
+bool BlockLatticeField2D::operator()(
+ T output[], const int input[])
+{
+ this->_blockLattice.get(input[0], input[1]).template computeField(output);
+ return true;
+}
+
+
+template
+BlockLatticePhysExternalPorosity2D::BlockLatticePhysExternalPorosity2D(
+ BlockLatticeStructure2D& blockLattice,
+ int overlap,
+ const UnitConverter& converter)
+ : BlockLatticePhysF2D(blockLattice, converter, 2),
+ _overlap(overlap)
+{
+ this->getName() = "ExtPorosityField";
+}
+
+template
+bool BlockLatticePhysExternalPorosity2D::operator() (T output[], const int input[])
+{
+ this->_blockLattice.get(
+ input[0]+_overlap, input[1]+_overlap
+ ).template computeField(output);
+ return true;
+}
+
+template
+BlockLatticePhysExternalVelocity2D::BlockLatticePhysExternalVelocity2D(
+ BlockLatticeStructure2D& blockLattice,
+ int overlap,
+ const UnitConverter& converter)
+ : BlockLatticePhysF2D(blockLattice, converter, 2),
+ _overlap(overlap)
+{
+ this->getName() = "ExtVelocityField";
+}
+
+template
+bool BlockLatticePhysExternalVelocity2D::operator() (T output[], const int input[])
+{
+ this->_blockLattice.get(
+ input[0]+_overlap, input[1]+_overlap
+ ).template computeField(output);
+ output[0] = this->_converter.getPhysVelocity(output[0]);
+ output[1] = this->_converter.getPhysVelocity(output[1]);
+ return true;
+}
+
+template
+BlockLatticePhysExternalParticleVelocity2D::BlockLatticePhysExternalParticleVelocity2D
+(BlockLatticeStructure2D& blockLattice, const UnitConverter& converter)
+ : BlockLatticePhysF2D(blockLattice,converter,2)
+{
+ this->getName() = "ExtParticleVelocityField";
+}
+
+template
+bool BlockLatticePhysExternalParticleVelocity2D::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]);
+ return true;
+ }
+ output[0]=this->_converter.getPhysVelocity(velocity_numerator[0]);
+ output[1]=this->_converter.getPhysVelocity(velocity_numerator[1]);
+ return true;
+}
+
+template
+BlockLatticePhysWallShearStress2D::BlockLatticePhysWallShearStress2D(
+ BlockLatticeStructure2D& blockLattice,
+ BlockGeometryStructure2D& blockGeometry,
+ int overlap,
+ int material,
+ const UnitConverter& converter,
+ IndicatorF2D& indicator)
+ : BlockLatticePhysF2D(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(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 origin(physR[0],physR[1]);
+ Vector direction(-_discreteNormal[iX-1][iY-1][0] * scaling,
+ -_discreteNormal[iX-1][iY-1][1] * scaling);
+ Vector 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
+bool BlockLatticePhysWallShearStress2D::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
+BlockLatticePhysBoundaryForce2D::BlockLatticePhysBoundaryForce2D(
+ BlockLatticeStructure2D& blockLattice,
+ BlockIndicatorF2D& indicatorF,
+ const UnitConverter& converter)
+ : BlockLatticePhysF2D(blockLattice, converter, 2),
+ _indicatorF(indicatorF),
+ _blockGeometry(indicatorF.getBlockGeometryStructure())
+{
+ this->getName() = "physBoundaryForce";
+}
+
+template
+bool BlockLatticePhysBoundaryForce2D::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(iPop,0), input[1] + descriptors::c(iPop,1)) == 1) {
+ // Get f_q of next fluid cell where l = opposite(q)
+ T f = this->_blockLattice.get(input[0] + descriptors::c(iPop,0), input[1] + descriptors::c(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(iPop)];
+ // Update force
+ for (int i = 0; i < this->getTargetDim(); ++i) {
+ output[i] -= descriptors::c(iPop,i) * f;
+ }
+ }
+ }
+ for (int i = 0; i < this->getTargetDim(); ++i) {
+ output[i] = this->_converter.getPhysForce(output[i]);
+ }
+ }
+ return true;
+}
+
+template
+BlockLatticePSMPhysForce2D::BlockLatticePSMPhysForce2D(
+ BlockLatticeStructure2D& blockLattice,
+ const UnitConverter& converter,
+ int mode_)
+ : BlockLatticePhysF2D(blockLattice, converter, 2)
+{
+ this->getName() = "physPSMForce";
+ mode = (Mode) mode_;
+}
+
+template
+bool BlockLatticePSMPhysForce2D::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());
+
+ //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())[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(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])[iPop])
+ + (1 - omega)
+ * (this->_blockLattice.get(input[0], input[1])[iPop]
+ - lbDynamicsHelpers::equilibrium(iPop, rho, u, uSqr));
+ break;
+ case M3:
+ omega_s =
+ (this->_blockLattice.get(input[0], input[1])[descriptors::opposite(iPop)]
+ - lbDynamicsHelpers::equilibrium(
+ descriptors::opposite(iPop), rho, u_s, uSqr_s))
+ - (this->_blockLattice.get(input[0], input[1])[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
+BlockLatticePhysCorrBoundaryForce2D::BlockLatticePhysCorrBoundaryForce2D
+(BlockLatticeStructure2D& blockLattice,BlockGeometry2D& blockGeometry,
+ int material, const UnitConverter& converter)
+ : BlockLatticePhysF2D(blockLattice,converter,2),
+ _blockGeometry(blockGeometry), _material(material)
+{
+ this->getName() = "physCorrBoundaryForce";
+}
+
+template
+bool BlockLatticePhysCorrBoundaryForce2D::operator() (T output[], const int input[])
+{
+ // int globIC = input[0];
+ // int locix= input[1];
+ // int lociy= input[2];
+ // int lociz= input[3];
+
+ // std::vector 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(iPop)];
+ // // Update force
+ // force[0] -= c[0]*(f-2.*descriptors::t(iPop));
+ // force[1] -= c[1]*(f-2.*descriptors::t(iPop));
+ // force[2] -= c[2]*(f-2.*descriptors::t(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
+BlockLatticePorosity2D::BlockLatticePorosity2D
+(BlockLatticeStructure2D& blockLattice, BlockGeometryStructure2D& blockGeometry,
+ int material, const UnitConverter& converter)
+ : BlockLatticePhysF2D(blockLattice,converter,1), _blockGeometry(blockGeometry), _material(material)
+{
+ this->getName() = "porosity";
+}
+
+// under construction
+template
+bool BlockLatticePorosity2D::operator()(T output[], const int input[])
+{
+ return false;
+}
+
+template
+BlockLatticeVolumeFractionApproximation2D::BlockLatticeVolumeFractionApproximation2D(BlockLatticeStructure2D& blockLattice,
+ BlockGeometryStructure2D& blockGeometry,
+ IndicatorF2D& indicator,
+ int refinementLevel,
+ const UnitConverter& converter,
+ bool insideOut)
+ : BlockLatticeF2D(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
+bool BlockLatticeVolumeFractionApproximation2D::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
+BlockLatticeVolumeFractionPolygonApproximation2D::BlockLatticeVolumeFractionPolygonApproximation2D(BlockLatticeStructure2D& blockLattice,
+ BlockGeometryStructure2D& blockGeometry,
+ IndicatorF2D& indicator,
+ const UnitConverter& converter,
+ bool insideOut)
+ : BlockLatticeF2D(blockLattice, 1),
+ _blockGeometry(blockGeometry), _indicator(indicator), _converter(converter), _insideOut(insideOut)
+{
+ this->getName() = "volumeFractionPolygonApproximation";
+}
+
+template
+bool BlockLatticeVolumeFractionPolygonApproximation2D::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 cornerXMYM_vec(physR[0] - 0.5*_converter.getPhysDeltaX(), physR[1] - 0.5*_converter.getPhysDeltaX());
+ Vector cornerXPYP_vec(physR[0] + 0.5*_converter.getPhysDeltaX(), physR[1] + 0.5*_converter.getPhysDeltaX());
+
+ Vector directionXP(_converter.getPhysDeltaX(), 0);
+ Vector directionXM(-_converter.getPhysDeltaX(), 0);
+ Vector directionYP(0, _converter.getPhysDeltaX());
+ Vector 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
+BlockLatticePhysPermeability2D::BlockLatticePhysPermeability2D
+(BlockLatticeStructure2D& blockLattice, BlockGeometryStructure2D& blockGeometry,
+ int material, const UnitConverter& converter)
+ : BlockLatticePhysF2D(blockLattice,converter,1), _blockGeometry(blockGeometry), _material(material)
+{
+ this->getName() = "permeability";
+}
+
+template
+bool BlockLatticePhysPermeability2D::operator()(T output[], const int input[])
+{
+ return false;
+}
+
+
+template
+BlockLatticePhysDarcyForce2D::BlockLatticePhysDarcyForce2D
+(BlockLatticeStructure2D& blockLattice, BlockGeometry2D& blockGeometry,
+ int material, const UnitConverter& converter)
+ : BlockLatticePhysF2D(blockLattice,converter,2),
+ _blockGeometry(blockGeometry), _material(material)
+{
+ this->getName() = "alphaU";
+}
+
+template
+bool BlockLatticePhysDarcyForce2D::operator()(T output[], const int input[])
+{
+ BlockLatticePhysPermeability2D permeability(this->_blockLattice,this->_blockGeometry,this->_material,this->_converter);
+ BlockLatticeVelocity2D velocity(this->_blockLattice);
+
+ T nu = this->_converter.getPhysViscosity();
+ permeability(output,input);
+ T K = output[0];
+ velocity(output,input);
+
+ output[0] *= -nu/K;
+ output[1] *= -nu/K;
+
+ return true;
+}
+
+
+template
+BlockLatticeAverage2D::BlockLatticeAverage2D
+(BlockLatticeF2D& f, BlockGeometry2D& blockGeometry,
+ int material, T radius)
+ : BlockLatticeF2D(f.getBlockLattice(), f.getTargetDim()),
+ _f(f), _blockGeometry(blockGeometry), _material(material), _radius(radius)
+{
+ this->getName() = "Average("+f.getName()+")";
+}
+
+
+template
+bool BlockLatticeAverage2D::operator() (T output[], const int input[])
+{
+ // CuboidGeometry2D& cGeometry = f.getBlockLattice2D().get_cGeometry();
+ // loadBalancer& load = f.getBlockLattice2D().get_load();
+
+ // //create boolean indicator functor isInSphere
+ // std::vector center(3,0);
+ // center[0] = (int)cGeometry.get(load.glob(input[0])).get_globPosX() + input[1];
+ // center[1] = (int)cGeometry.get(load.glob(input[0])).get_globPosY() + input[2];
+ // center[2] = (int)cGeometry.get(load.glob(input[0])).get_globPosZ() + input[3];
+ // SphereAnalyticalF2D isInSphere(center,radius);
+
+ // iterate over all cuboids & points and test for material && isInSphere
+ // std::vector tmp( this->_n, T() );
+ // int numVoxels(0);
+ // if (this->_blockGeometry.getMaterial(center[0],center[1],center[2]) == material) {
+ // for (int iC=0; iC glob(3,0);
+ // glob[0] = (int)cGeometry.get(load.glob(iC)).get_globPosX() + iX;
+ // glob[1] = (int)cGeometry.get(load.glob(iC)).get_globPosY() + iY;
+ // glob[2] = (int)cGeometry.get(load.glob(iC)).get_globPosZ() + iZ;
+ // if (this->_blockGeometry.getMaterial(glob[0],glob[1],glob[2]) == material
+ // && isInSphere(glob)[0]==true) {
+ // for (unsigned iD=0; iD0) {
+ // tmp[iD] /= numVoxels;
+ // }
+ // }
+ // }
+ // return tmp;
+
+ return false;
+}
+
+
+template