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 +BlockEuklidNorm2D::BlockEuklidNorm2D(BlockF2D& f) + : BlockF2D(f.getBlockStructure(),1), _f(f) +{ + this->getName() = "l2("+f.getName()+")"; +} + +template +bool BlockEuklidNorm2D::operator() (T output[], const int input[]) +{ + output[0] = T(); // flash output, since this methods adds values. + T data[_f.getTargetDim()]; + _f(data,input); + for ( int i = 0; i < _f.getTargetDim(); ++i) { + output[0] += data[i]*data[i]; + } + output[0] = sqrt(output[0]); + return true; +} + +template +BlockLatticePorousMomentumLossForce2D::BlockLatticePorousMomentumLossForce2D( + BlockLatticeStructure2D& blockLattice, BlockGeometryStructure2D& blockGeometry, + std::vector* >& indicator, const UnitConverter& converter) + : BlockLatticePhysF2D(blockLattice, converter, 4*indicator.size()), _blockGeometry(blockGeometry), _vectorOfIndicator(indicator) +{ + this->getName() = "physPorousMomentumLossForce"; +} + +template +bool BlockLatticePorousMomentumLossForce2D::operator()(T output[], const int input[]) +{ + // iterate over all particles in _indicator + for (size_t iInd=0; iInd!=_vectorOfIndicator.size(); iInd++) { + + int numVoxels = 0; + int start[2] = {0}; + int end[2] = {0}; + // check for intersection of cuboid and indicator + Cuboid2D tmpCuboid(_blockGeometry.getOrigin()[0], _blockGeometry.getOrigin()[1], this->_converter.getPhysDeltaX(), _blockGeometry.getNx(), _blockGeometry.getNy()); + T posXmin = _vectorOfIndicator[iInd]->getPos()[0] - _vectorOfIndicator[iInd]->getCircumRadius(); + T posXmax = _vectorOfIndicator[iInd]->getPos()[0] + _vectorOfIndicator[iInd]->getCircumRadius(); + T posYmin = _vectorOfIndicator[iInd]->getPos()[1] - _vectorOfIndicator[iInd]->getCircumRadius(); + T posYmax = _vectorOfIndicator[iInd]->getPos()[1] + _vectorOfIndicator[iInd]->getCircumRadius(); + if (tmpCuboid.checkInters(posXmin, posXmax, posYmin, posYmax, start[0], end[0], start[1], end[1])) { + + T invDeltaX = 1./this->_converter.getPhysDeltaX(); + for(int k=0; k<2; k++) { + start[k] -= 1; + if (start[k] < 0) start[k] = 0; + end[k] += 2; + if (end[k] > _blockGeometry.getExtend()[k]) end[k] = _blockGeometry.getExtend()[k]; + } + // iterate over cells in the constructed intersection box + for (int iX = start[0]; iX < end[0]; iX++) { + for (int iY = start[1]; iY < end[1]; iY++) { + // check if cell belongs to particle + T inside[1] = {0.}; + T posIn[2] = {0.}; + _blockGeometry.getPhysR(posIn, iX, iY); + (*(_vectorOfIndicator[iInd]))( inside, posIn); + if ( !util::nearZero(inside[0]) && this->_blockGeometry.get(iX,iY)==1) { + // compute momentum exchange force on particle + T tmpForce[2] = {0.,0.}; + tmpForce[0] += this->_blockLattice.get(iX, iY).template getFieldPointer()[0]; + tmpForce[1] += this->_blockLattice.get(iX, iY).template getFieldPointer()[1]; + // reset external field for next timestep + T reset_to_zero[2] = {0.,0.}; + this->_blockLattice.get(iX, iY).template setField(reset_to_zero); + // convert force to SI units and compute torque + numVoxels++; + // division bei length of lattice cell necessary due to converter handling of force + tmpForce[0] = this->_converter.getPhysForce(tmpForce[0])*invDeltaX; + tmpForce[1] = this->_converter.getPhysForce(tmpForce[1])*invDeltaX; + output[0+iInd*4] += tmpForce[0]; + output[1+iInd*4] += tmpForce[1]; + output[2+iInd*4] += (posIn[0]-_vectorOfIndicator[iInd]->getPos()[0])*tmpForce[1] - (posIn[1]-_vectorOfIndicator[iInd]->getPos()[1])*tmpForce[0]; + } + } + } + + } + output[3+iInd*4] = numVoxels; + + } + return true; +} + +template +BlockLatticePhysTemperature2D::BlockLatticePhysTemperature2D +(BlockLatticeStructure2D& blockLattice, ThermalUnitConverter const& converter) + : BlockLatticeThermalPhysF2D(blockLattice,converter,1) +{ + this->getName() = "physTemperature"; +} + + +template +bool BlockLatticePhysTemperature2D::operator() (T output[], const int input[]) +{ + T latticeTemperature = this->_blockLattice.get( input[0], input[1] ).computeRho(); + output[0] = this->_converter.getPhysTemperature(latticeTemperature); + + return true; +} + +template +BlockLatticePhysHeatFlux2D::BlockLatticePhysHeatFlux2D +(BlockLatticeStructure2D& blockLattice, ThermalUnitConverter const& converter) + : BlockLatticeThermalPhysF2D(blockLattice,converter,2), + _temp(converter.getLatticeSpecificHeatCapacity(converter.getPhysSpecificHeatCapacity())*(converter.getLatticeThermalRelaxationTime() - 0.5) / converter.getLatticeThermalRelaxationTime()) +{ + this->getName() = "physHeatFlux"; + + // std::cout << _temp << " " << converter.getConversionFactorHeatFlux() << " " << 0.426 / _temp / converter.getConversionFactorHeatFlux() << std::endl; +} + +template +bool BlockLatticePhysHeatFlux2D::operator() (T output[], const int input[]) +{ + T temperature, extVel[2], j[2]; + this->_blockLattice.get( input[0], input[1] ).computeRhoU(temperature,extVel); + this->_blockLattice.get( input[0], input[1] ).computeJ(j); + + output[0] = this->_converter.getPhysHeatFlux((j[0] - temperature * extVel[0])*_temp); + output[1] = this->_converter.getPhysHeatFlux((j[1] - temperature * extVel[1])*_temp); + + return true; +} + +template +BlockLatticeIndicatorSmoothIndicatorIntersection2D::BlockLatticeIndicatorSmoothIndicatorIntersection2D ( + BlockLatticeStructure2D& blockLattice, + BlockGeometryStructure2D& blockGeometry, + IndicatorF2D& normalInd, SmoothIndicatorF2D& smoothInd ) + : BlockLatticeF2D(blockLattice, 1), + _blockGeometry(blockGeometry), _normalInd(normalInd), _smoothInd(smoothInd) +{ + this->getName() = "Indicator-SmoothIndicator Intersection"; +} + +template +bool BlockLatticeIndicatorSmoothIndicatorIntersection2D::operator()(T output[], const int input[]) +{ + output[0] = 0.; + int start[2] = {0}; + int end[2] = {0}; + // check for intersection of cuboid and smoothIndicator + Cuboid2D tmpCuboid(_blockGeometry.getOrigin()[0], _blockGeometry.getOrigin()[1], _blockGeometry.getDeltaR(), _blockGeometry.getNx(), _blockGeometry.getNy()); + T posXmin = _smoothInd.getPos()[0] - _smoothInd.getCircumRadius(); + T posXmax = _smoothInd.getPos()[0] + _smoothInd.getCircumRadius(); + T posYmin = _smoothInd.getPos()[1] - _smoothInd.getCircumRadius(); + T posYmax = _smoothInd.getPos()[1] + _smoothInd.getCircumRadius(); + if (tmpCuboid.checkInters(posXmin, posXmax, posYmin, posYmax, start[0], end[0], start[1], end[1])) { + + for(int k=0; k<2; k++) { + start[k] -= 1; + if (start[k] < 0) start[k] = 0; + end[k] += 2; + if (end[k] > _blockGeometry.getExtend()[k]) end[k] = _blockGeometry.getExtend()[k]; + } + + // iterate over cells in the constructed intersection box + for (int iX = start[0]; iX < end[0]; iX++) { + for (int iY = start[1]; iY < end[1]; iY++) { + + // check if cell belongs to particle + T insideT[1] = {0.}; + T posIn[2] = {0.}; + _blockGeometry.getPhysR(posIn, iX, iY); + _smoothInd( insideT, posIn); + if ( !util::nearZero(insideT[0]) && this->_blockGeometry.get(iX,iY)==1) { + // Return 1 if at least one cell is found to be inside both A and B + bool insideBool[1] = {false}; + _normalInd(insideBool, posIn); + if (insideBool[0]) { + output[0] = 1.; + return true; + } + } + } + } + } + + return true; +} + +/////////// BlockLatticeGuoZhaoEpsilon2D ///////////////////////////////////////////// + +template +BlockLatticeGuoZhaoEpsilon2D::BlockLatticeGuoZhaoEpsilon2D +(BlockLatticeStructure2D& blockLattice) + : BlockLatticeF2D(blockLattice, 1) +{ + this->getName() = "epsilon"; +} + +template +bool BlockLatticeGuoZhaoEpsilon2D::operator() (T output[], const int input[]) +{ + output[0] = *this->_blockLattice.get( input[0], input[1] ).template getFieldPointer(); + return true; +} + +/////////// BlockLatticeGuoZhaoPhysK2D ///////////////////////////////////////////// + +template +BlockLatticeGuoZhaoPhysK2D::BlockLatticeGuoZhaoPhysK2D +(BlockLatticeStructure2D& blockLattice, const UnitConverter& converter) + : BlockLatticePhysF2D(blockLattice, converter, 1) +{ + this->getName() = "physK"; +} + +template +bool BlockLatticeGuoZhaoPhysK2D::operator() (T output[], const int input[]) +{ + output[0] = *this->_blockLattice.get( input[0], input[1] ) + .template getFieldPointer() + * this->_converter.getConversionFactorLength() + * this->_converter.getConversionFactorLength(); + return true; +} + +/////////// BlockLatticeGuoZhaoPhysBodyForce2D ///////////////////////////////////////////// + +template +BlockLatticeGuoZhaoPhysBodyForce2D::BlockLatticeGuoZhaoPhysBodyForce2D +(BlockLatticeStructure2D& blockLattice, const UnitConverter& converter) + : BlockLatticePhysF2D(blockLattice, converter, 2) +{ + this->getName() = "physBodyForce"; +} + +template +bool BlockLatticeGuoZhaoPhysBodyForce2D::operator() (T output[], const int input[]) +{ + output[0] = this->_converter.getPhysForce( *this->_blockLattice.get( input[0], input[1] ) + .template getFieldPointer() ); + output[1] = this->_converter.getPhysForce( *this->_blockLattice.get( input[0], input[1] ) + .template getFieldPointer(1) ); + return true; +} + +} // end namespace olb + +#endif + -- cgit v1.2.3