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/boundary/offBoundaryInstantiator2D.h | 715 +++++++++++++++++++++++++++++++ 1 file changed, 715 insertions(+) create mode 100644 src/boundary/offBoundaryInstantiator2D.h (limited to 'src/boundary/offBoundaryInstantiator2D.h') diff --git a/src/boundary/offBoundaryInstantiator2D.h b/src/boundary/offBoundaryInstantiator2D.h new file mode 100644 index 0000000..f72eae8 --- /dev/null +++ b/src/boundary/offBoundaryInstantiator2D.h @@ -0,0 +1,715 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2012 Jonas Kratzke, Mathias J. Krause + * 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. + */ + +/** \file + * A helper for initialising 2D boundaries -- header file. + */ + +#ifndef OFF_BOUNDARY_INSTANTIATOR_2D_H +#define OFF_BOUNDARY_INSTANTIATOR_2D_H + +#include "offBoundaryCondition2D.h" +#include "offBoundaryPostProcessors2D.h" +#include "geometry/blockGeometry2D.h" +#include "geometry/blockGeometryStatistics2D.h" +#include "dynamics/dynamics.h" +#include "core/cell.h" +#include "io/ostreamManager.h" +#include "core/util.h" +#include "io/stlReader.h" +#include "functors/lattice/indicator/blockIndicatorF2D.h" + +namespace olb { + +/** +* This class gets the needed processors from BoundaryManager and adds them +* to the Processor Vector of the DESCRIPTOR +*/ + +template +class OffBoundaryConditionInstantiator2D: public OffLatticeBoundaryCondition2D { +public: + OffBoundaryConditionInstantiator2D(BlockLatticeStructure2D& block_, T epsFraction_ = 0.0001); + ~OffBoundaryConditionInstantiator2D() override; + + void addOnePointZeroVelocityBoundary(int x, int y, int iPop, T dist) override; + void addTwoPointZeroVelocityBoundary(int x, int y, int iPop, T dist) override; + + void addOnePointVelocityBoundary(int x, int y, int iPop, T dist) override; + void addTwoPointVelocityBoundary(int x, int y, int iPop, T dist) override; + + void addOffDynamics(int x, int y, T location[DESCRIPTOR::d]) override; + void addOffDynamics(int x, int y, T location[DESCRIPTOR::d], T distances[DESCRIPTOR::q]) override; + void addOffDynamics(BlockIndicatorF2D& indicator) override; + + void addZeroVelocityBoundary(BlockGeometryStructure2D& blockGeometryStructure, int x, int y, int iPop, T dist) override; + void addZeroVelocityBoundary(BlockGeometryStructure2D& blockGeometryStructure, int x, int y, T distances[DESCRIPTOR::q]); + void addZeroVelocityBoundary(BlockGeometryStructure2D& blockGeometryStructure, int iX, int iY, BlockIndicatorF2D& bulkIndicator, IndicatorF2D& geometryIndicator) override; + void addZeroVelocityBoundary(BlockIndicatorF2D& boundaryIndicator, BlockIndicatorF2D& bulkIndicator, IndicatorF2D& geometryIndicator) override; + void addZeroVelocityBoundary(BlockIndicatorF2D& boundaryIndicator, BlockIndicatorF2D& bulkIndicator) override; + + void addVelocityBoundary(BlockGeometryStructure2D& blockGeometryStructure, int x, int y, int iPop, T dist) override; + void addVelocityBoundary(BlockGeometryStructure2D& blockGeometryStructure, int x, int y, T distances[DESCRIPTOR::q]); + void addVelocityBoundary(BlockGeometryStructure2D& blockGeometryStructure, int iX, int iY, BlockIndicatorF2D& bulkIndicator, IndicatorF2D& geometryIndicator) override; + void addVelocityBoundary(BlockIndicatorF2D& boundaryIndicator, BlockIndicatorF2D& bulkIndicator, IndicatorF2D& geometryIndicator) override; + void addVelocityBoundary(BlockIndicatorF2D& boundaryIndicator, BlockIndicatorF2D& bulkIndicator) override; + + void addPressureBoundary(BlockIndicatorF2D& boundaryIndicator, BlockIndicatorF2D& bulkIndicator, IndicatorF2D& geometryIndicator) override; + void addPressureBoundary(BlockIndicatorF2D& boundaryIndicator, BlockIndicatorF2D& bulkIndicator) override; + + void setBoundaryIntersection(int iX, int iY, int iPop, T distance); + bool getBoundaryIntersection(int iX, int iY, int iPop, T point[DESCRIPTOR::d]); + + void defineU(int iX, int iY, int iPop, const T u[DESCRIPTOR::d]) override; + void defineU(BlockIndicatorF2D& indicator, BlockIndicatorF2D& bulkIndicator, AnalyticalF2D& u) override; + + void defineRho(int iX, int iY, int iPop, const T rho) override; + void defineRho(BlockIndicatorF2D& indicator, BlockIndicatorF2D& bulkIndicator, AnalyticalF2D& rho) override; + + void outputOn() override; + void outputOff() override; + + BlockLatticeStructure2D& getBlock() override; + BlockLatticeStructure2D const& getBlock() const override; +private: + BlockLatticeStructure2D& block; + //std::vector*> momentaVector; + std::vector*> dynamicsVector; + bool _output; + mutable OstreamManager clout; +}; + +///////// class OffBoundaryConditionInstantiator2D //////////////////////// + +template +BlockLatticeStructure2D& OffBoundaryConditionInstantiator2D::getBlock() +{ + return block; +} + +template +BlockLatticeStructure2D const& OffBoundaryConditionInstantiator2D::getBlock() const +{ + return block; +} + +template +OffBoundaryConditionInstantiator2D::OffBoundaryConditionInstantiator2D( + BlockLatticeStructure2D& block_, T epsFraction_) : + block(block_), _output(false), clout(std::cout,"BoundaryConditionInstantiator2D") +{ + this->_epsFraction = epsFraction_; +} + +template +OffBoundaryConditionInstantiator2D::~OffBoundaryConditionInstantiator2D() +{ + for (unsigned iDynamics = 0; iDynamics < dynamicsVector.size(); ++iDynamics) { + delete dynamicsVector[iDynamics]; + } + /* + for (unsigned iMomenta = 0; iMomenta < dynamicsVector.size(); ++iMomenta) { + delete momentaVector[iMomenta]; + }*/ +} + +template +void OffBoundaryConditionInstantiator2D::addOnePointZeroVelocityBoundary( + int x, int y, int iPop, T dist) +{ + PostProcessorGenerator2D* postProcessor = + BoundaryManager::getOnePointZeroVelocityBoundaryProcessor + (x, y, iPop, dist); + if (postProcessor) { + this->getBlock().addPostProcessor(*postProcessor); + } +} + +template +void OffBoundaryConditionInstantiator2D::addTwoPointZeroVelocityBoundary( + int x, int y, int iPop, T dist) +{ + PostProcessorGenerator2D* postProcessor = + BoundaryManager::getTwoPointZeroVelocityBoundaryProcessor + (x, y, iPop, dist); + if (postProcessor) { + this->getBlock().addPostProcessor(*postProcessor); + } +} + +template +void OffBoundaryConditionInstantiator2D::addOnePointVelocityBoundary( + int x, int y, int iPop, T dist) +{ + PostProcessorGenerator2D* postProcessor = + BoundaryManager::getOnePointVelocityBoundaryProcessor + (x, y, iPop, dist); + if (postProcessor) { + this->getBlock().addPostProcessor(*postProcessor); + } +} + +template +void OffBoundaryConditionInstantiator2D::addTwoPointVelocityBoundary( + int x, int y, int iPop, T dist) +{ + PostProcessorGenerator2D* postProcessor = + BoundaryManager::getTwoPointVelocityBoundaryProcessor + (x, y, iPop, dist); + if (postProcessor) { + this->getBlock().addPostProcessor(*postProcessor); + } +} + +template +void OffBoundaryConditionInstantiator2D::addOffDynamics( + int x, int y, T location[DESCRIPTOR::d]) +{ + Dynamics* dynamics + = BoundaryManager::getOffDynamics(location); + this->getBlock().defineDynamics(x,x,y,y, dynamics); + dynamicsVector.push_back(dynamics); +} + +template +void OffBoundaryConditionInstantiator2D::addOffDynamics( + int x, int y, T location[DESCRIPTOR::d], T distances[DESCRIPTOR::q]) +{ + Dynamics* dynamics + = BoundaryManager::getOffDynamics(location, distances); + this->getBlock().defineDynamics(x,x,y,y, dynamics); + dynamicsVector.push_back(dynamics); +} + +template +void OffBoundaryConditionInstantiator2D::addOffDynamics( + BlockIndicatorF2D& indicator) +{ + if ( !indicator.isEmpty() ) { + const Vector min = indicator.getMin(); + const Vector max = indicator.getMax(); + + for (int iX = min[0]; iX <= max[0]; ++iX) { + for (int iY = min[1]; iY <= max[1]; ++iY) { + if (indicator(iX,iY)) { + T location[DESCRIPTOR::d]; + indicator.getBlockGeometryStructure().getPhysR(location, iX,iY); + addOffDynamics(iX, iY, location); + } + } + } + } +} + +template +void OffBoundaryConditionInstantiator2D::addZeroVelocityBoundary( + BlockGeometryStructure2D& blockGeometryStructure, int x, int y, int iPop, T dist) +{ + if (blockGeometryStructure.getMaterial(x-descriptors::c(iPop,0), y-descriptors::c(iPop,1)) != 1) { + addOnePointZeroVelocityBoundary(x, y, iPop, dist); + } + else { + addTwoPointZeroVelocityBoundary(x, y, iPop, dist); + } +} + +template +void OffBoundaryConditionInstantiator2D:: +addZeroVelocityBoundary(BlockGeometryStructure2D& blockGeometryStructure, int x, int y, T distances[DESCRIPTOR::q]) +{ + typedef DESCRIPTOR L; + //T location[DESCRIPTOR::d]; + //location[0] = blockGeometryStructure.physCoordX(x); + //location[1] = blockGeometryStructure.physCoordY(y); + //location[2] = blockGeometryStructure.physCoordZ(z); + //T distancesCopy[L::q]; + //T spacing = blockGeometryStructure.getDeltaR(); + //for (int iPop = 1; iPop < L::q ; ++iPop) { + // distancesCopy[iPop] = spacing*(1.-distances[iPop]); + // if (distances[iPop] == -1) + // distancesCopy[iPop] = -1; + //} + //addOffDynamics(x, y, z, location, distancesCopy); + + for (int iPop = 1; iPop < L::q ; ++iPop) { + if ( !util::nearZero(distances[iPop]+1) ) { + addZeroVelocityBoundary(blockGeometryStructure, x-descriptors::c(iPop,0), y-descriptors::c(iPop,1), iPop, distances[iPop]); + } + } +} + +template +void OffBoundaryConditionInstantiator2D::addZeroVelocityBoundary( + BlockGeometryStructure2D& blockGeometryStructure, int iX, int iY, + BlockIndicatorF2D& bulkIndicator, IndicatorF2D& geometryIndicator) +{ + T distances[DESCRIPTOR::q]; + for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) { + distances[iPop] = -1; + } + + for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) { + const int iXn = iX + descriptors::c(iPop,0); + const int iYn = iY + descriptors::c(iPop,1); + if (bulkIndicator(iXn,iYn)) { + T dist = -1; + + T physR[2]; + blockGeometryStructure.getPhysR(physR,iXn,iYn); + T voxelSize=blockGeometryStructure.getDeltaR(); + + Vector physC(physR); + + Vector direction(-voxelSize*descriptors::c(iPop,0),-voxelSize*descriptors::c(iPop,1)) ; + T cPhysNorm = voxelSize*sqrt(descriptors::c(iPop,0)*descriptors::c(iPop,0)+descriptors::c(iPop,1)*descriptors::c(iPop,1)); + + if (!geometryIndicator.distance(dist,physC,direction,blockGeometryStructure.getIcGlob() ) ) { + T epsX = voxelSize*descriptors::c(iPop,0)*this->_epsFraction; + T epsY = voxelSize*descriptors::c(iPop,1)*this->_epsFraction; + + Vector physC2(physC); + physC2[0] += epsX; + physC2[1] += epsY; + Vector direction2(direction); + direction2[0] -= 2.*epsX; + direction2[1] -= 2.*epsY; + + if ( !geometryIndicator.distance(dist,physC2,direction2,blockGeometryStructure.getIcGlob())) { + clout << "ERROR: no boundary found at (" << iXn << "," << iYn <<") ~ (" + << physR[0] << "," << physR[1] << "), " + << "in direction " << util::opposite(iPop) + << std::endl; + } + T distNew = (dist - sqrt(epsX*epsX+epsY*epsY) )/cPhysNorm; + if (distNew < 0.5) { + dist = 0; + } + else { + dist = 0.5 * cPhysNorm; + clout << "WARNING: distance at (" << iXn << "," << iYn << ") ~ (" + << physR[0] << "," << physR[1] << "), " + << "in direction " << util::opposite(iPop) << ": " + << distNew + << " rounded to " + << dist/cPhysNorm + << std::endl; + } + } + distances[util::opposite(iPop)] = dist/cPhysNorm; + } // bulk indicator if + } // iPop loop + addZeroVelocityBoundary(blockGeometryStructure, iX, iY, distances); +} + +template +void OffBoundaryConditionInstantiator2D::addZeroVelocityBoundary( + BlockIndicatorF2D& boundaryIndicator, + BlockIndicatorF2D& bulkIndicator, + IndicatorF2D& geometryIndicator) +{ + if ( !boundaryIndicator.isEmpty() ) { + const Vector min = boundaryIndicator.getMin(); + const Vector max = boundaryIndicator.getMax(); + + for (int iX = min[0]; iX <= max[0]; ++iX) { + for (int iY = min[1]; iY <= max[1]; ++iY) { + if (boundaryIndicator(iX,iY)) { + addZeroVelocityBoundary(bulkIndicator.getBlockGeometryStructure(), iX, iY, + bulkIndicator, geometryIndicator); + } + } + } + } +} + +template +void OffBoundaryConditionInstantiator2D::addZeroVelocityBoundary( + BlockIndicatorF2D& boundaryIndicator, BlockIndicatorF2D& bulkIndicator) +{ + if ( boundaryIndicator.isEmpty() ) { + return; + } + + auto& blockGeometryStructure = boundaryIndicator.getBlockGeometryStructure(); + const Vector min = boundaryIndicator.getMin(); + const Vector max = boundaryIndicator.getMax(); + + for (int iX = min[0]; iX <= max[0]; ++iX) { + for (int iY = min[1]; iY <= max[1]; ++iY) { + if (boundaryIndicator(iX, iY)) { + bool streamDirections[DESCRIPTOR::q]; + // check if (ix,iY) is really a boundary node and compute the stream direction + if (blockGeometryStructure.template findStreamDirections(iX, iY, boundaryIndicator, bulkIndicator, streamDirections) ) { + T distances[DESCRIPTOR::q]; + for (int iPop=0; iPop(iPop)] = .5; + } + else { + distances[descriptors::opposite(iPop)] = -1; + } + } + addZeroVelocityBoundary(blockGeometryStructure, iX, iY, distances); + } + } + } + } +} + +template +void OffBoundaryConditionInstantiator2D::addVelocityBoundary( + BlockGeometryStructure2D& blockGeometryStructure, int x, int y, int iPop, T dist) +{ + if (blockGeometryStructure.getMaterial(x-descriptors::c(iPop,0), y-descriptors::c(iPop,1)) != 1) { + addOnePointVelocityBoundary(x, y, iPop, dist); + } + else { + addTwoPointVelocityBoundary(x, y, iPop, dist); + } +} + +template +void OffBoundaryConditionInstantiator2D:: +addVelocityBoundary(BlockGeometryStructure2D& blockGeometryStructure, int x, int y, T distances[DESCRIPTOR::q]) +{ + typedef DESCRIPTOR L; + T location[DESCRIPTOR::d]; + blockGeometryStructure.getPhysR(location, x,y); + + T distancesCopy[L::q]; + T spacing = blockGeometryStructure.getDeltaR(); + for (int iPop = 1; iPop < L::q ; ++iPop) { + distancesCopy[iPop] = spacing*(1.-distances[iPop]); + if ( !util::nearZero(distances[iPop]+1) ) { + distancesCopy[iPop] = -1; + } + } + addOffDynamics(x, y, location, distancesCopy); + + for (int iPop = 1; iPop < L::q ; ++iPop) { + if ( !util::nearZero(distances[iPop]+1) ) { + addVelocityBoundary(blockGeometryStructure, x-descriptors::c(iPop,0), y-descriptors::c(iPop,1), iPop, distances[iPop]); + } + } +} + +template +void OffBoundaryConditionInstantiator2D::addVelocityBoundary( + BlockGeometryStructure2D& blockGeometryStructure, int iX, int iY, + BlockIndicatorF2D& bulkIndicator, IndicatorF2D& geometryIndicator) +{ + T distances[DESCRIPTOR::q]; + for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) { + distances[iPop] = -1; + } + + for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) { + const int iXn = iX + descriptors::c(iPop,0); + const int iYn = iY + descriptors::c(iPop,1); + if (bulkIndicator(iXn,iYn)) { + T dist = -1; + T physR[2]; + blockGeometryStructure.getPhysR(physR,iXn,iYn); + T voxelSize=blockGeometryStructure.getDeltaR(); + Vector physC(physR); + + Vector direction(-voxelSize*descriptors::c(iPop,0),-voxelSize*descriptors::c(iPop,1)); + T cPhysNorm = voxelSize*sqrt(descriptors::c(iPop,0)*descriptors::c(iPop,0)+descriptors::c(iPop,1)*descriptors::c(iPop,1)); + + if (!geometryIndicator.distance(dist,physC,direction,blockGeometryStructure.getIcGlob() ) ) { + T epsX = voxelSize*descriptors::c(iPop,0)*this->_epsFraction; + T epsY = voxelSize*descriptors::c(iPop,1)*this->_epsFraction; + + Vector physC2(physC); + physC2[0] += epsX; + physC2[1] += epsY; + Vector direction2(direction); + direction2[0] -= 2.*epsX; + direction2[1] -= 2.*epsY; + + if ( !geometryIndicator.distance(dist,physC2,direction2,blockGeometryStructure.getIcGlob())) { + clout << "ERROR: no boundary found at (" << iXn << "," << iYn <<") ~ (" + << physR[0] << "," << physR[1] << "), " + << "in direction " << util::opposite(iPop) + << std::endl; + } + T distNew = (dist - sqrt(epsX*epsX+epsY*epsY))/cPhysNorm; + if (distNew < 0.5) { + dist = 0; + } + else { + dist = 0.5 * cPhysNorm; + clout << "WARNING: distance at (" << iXn << "," << iYn <<") ~ (" + << physR[0] << "," << physR[1] <<"), " + << "in direction " << util::opposite(iPop) << ": " + << distNew + << " rounded to " + << dist/cPhysNorm + << std::endl; + } + } + distances[util::opposite(iPop)] = dist/cPhysNorm; + } // bulk indicator + } // iPop loop + addVelocityBoundary(blockGeometryStructure, iX, iY, distances); +} + +template +void OffBoundaryConditionInstantiator2D::addVelocityBoundary( + BlockIndicatorF2D& boundaryIndicator, BlockIndicatorF2D& bulkIndicator, IndicatorF2D& geometryIndicator) +{ + if ( !boundaryIndicator.isEmpty() ) { + const Vector min = boundaryIndicator.getMin(); + const Vector max = boundaryIndicator.getMax(); + + for (int iX = min[0]; iX <= max[0]; ++iX) { + for (int iY = min[1]; iY <= max[1]; ++iY) { + if (boundaryIndicator(iX, iY)) { + addVelocityBoundary(bulkIndicator.getBlockGeometryStructure(), iX, iY, + bulkIndicator, geometryIndicator); + } + } + } + } +} + +template +void OffBoundaryConditionInstantiator2D::addVelocityBoundary( + BlockIndicatorF2D& boundaryIndicator, BlockIndicatorF2D& bulkIndicator) +{ + if ( boundaryIndicator.isEmpty() ) { + return; + } + + auto& blockGeometryStructure = boundaryIndicator.getBlockGeometryStructure(); + const Vector min = boundaryIndicator.getMin(); + const Vector max = boundaryIndicator.getMax(); + + for (int iX = min[0]; iX <= max[0]; ++iX) { + for (int iY = min[1]; iY <= max[1]; ++iY) { + if (boundaryIndicator(iX,iY)) { + bool streamDirections[DESCRIPTOR::q]; + // check if (ix,iY) is really a boundary node and compute the stream direction + if (blockGeometryStructure.template findStreamDirections(iX, iY, boundaryIndicator, bulkIndicator, streamDirections) ) { + T distances[DESCRIPTOR::q]; + for (int iPop=0; iPop(iPop)] = .5; + } + else { + distances[descriptors::opposite(iPop)] = -1; + } + } + addVelocityBoundary(blockGeometryStructure, iX, iY, distances); + } + } + } + } +} + +template +void OffBoundaryConditionInstantiator2D::addPressureBoundary( + BlockIndicatorF2D& boundaryIndicator, BlockIndicatorF2D& bulkIndicator, IndicatorF2D& geometryIndicator) +{ + auto& blockGeometryStructure = boundaryIndicator.getBlockGeometryStructure(); + if ( boundaryIndicator.isEmpty() ) { + const Vector min = boundaryIndicator.getMin(); + const Vector max = boundaryIndicator.getMax(); + + for (int iX = min[0]; iX <= max[0]; ++iX) { + for (int iY = min[1]; iY <= max[1]; ++iY) { + bool streamDirections[DESCRIPTOR::q]; + // check if (ix,iY) is really a boundary node and compute the stream direction + if (blockGeometryStructure.template findStreamDirections(iX, iY, boundaryIndicator, bulkIndicator, streamDirections) ) { + T location[DESCRIPTOR::d]; + blockGeometryStructure.getPhysR(location, iX, iX); + block.defineDynamics(iX,iX,iY,iY,&instances::getBounceBackAnti(1.)); + } + } + } + } + clout << "ERROR: OffBoundaryConditionInstantiator2D::addPressureBoundary NOT YET IMPLEMENTED" << std::endl; + exit(-1); +} + +template +void OffBoundaryConditionInstantiator2D::addPressureBoundary( + BlockIndicatorF2D& boundaryIndicator, BlockIndicatorF2D& bulkIndicator) +{ + if ( boundaryIndicator.isEmpty() ) { + return; + } + + auto& blockGeometryStructure = boundaryIndicator.getBlockGeometryStructure(); + const Vector min = boundaryIndicator.getMin(); + const Vector max = boundaryIndicator.getMax(); + + for (int iX = min[0]; iX <= max[0]; ++iX) { + for (int iY = min[1]; iY <= max[1]; ++iY) { + bool streamDirections[DESCRIPTOR::q]; + // check if (ix,iY) is really a boundary node and compute the stream direction + if (blockGeometryStructure.template findStreamDirections(iX, iY, boundaryIndicator, bulkIndicator, streamDirections) ) { + T location[DESCRIPTOR::d]; + blockGeometryStructure.getPhysR(location, iX, iX); + block.defineDynamics(iX,iX,iY,iY,&instances::getBounceBackAnti(1.)); + for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) { + if (streamDirections[iPop]) { + PostProcessorGenerator2D* postProcessor = new AntiBounceBackPostProcessorGenerator2D(iX+descriptors::c(iPop,0), iY+descriptors::c(iPop,1), descriptors::opposite(iPop)); + this->getBlock().addPostProcessor(*postProcessor); + } + } + } + } + } +} + +template +void OffBoundaryConditionInstantiator2D:: +setBoundaryIntersection(int iX, int iY, int iPop, T distance) +{ + + this->getBlock().getDynamics(iX, iY)->setBoundaryIntersection(iPop, distance); + if (_output) { + clout << "setBoundaryIntersection(" << iX << ", " << iY << " )" << std::endl; + } +} + +template +bool OffBoundaryConditionInstantiator2D:: +getBoundaryIntersection(int iX, int iY, int iPop, T point[DESCRIPTOR::d]) +{ + + return this->getBlock().getDynamics(iX, iY)->getBoundaryIntersection(iPop, point); +} + + +template +void OffBoundaryConditionInstantiator2D:: +defineU(int iX, int iY, int iPop, const T u[DESCRIPTOR::d]) +{ + + this->getBlock().getDynamics(iX, iY)->defineU(iPop, u); + if (_output) { + clout << "defineU(" << iX << ", " << iY << " )" << std::endl; + } +} + +template +void OffBoundaryConditionInstantiator2D:: +defineU(BlockIndicatorF2D& indicator, + BlockIndicatorF2D& bulkIndicator, + AnalyticalF2D& u) +{ + if ( indicator.isEmpty() ) { + return; + } + + const Vector min = indicator.getMin(); + const Vector max = indicator.getMax(); + + for (int iX = min[0]; iX <= max[0]; ++iX) { + for (int iY = min[1]; iY <= max[1]; ++iY) { + if (indicator(iX, iY)) { + for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) { + // Get direction + int iXn = iX + descriptors::c(iPop,0); + int iYn = iY + descriptors::c(iPop,1); + if (bulkIndicator(iXn, iYn)) { + T intersection[] = { T(), T() }; // coord. of intersection + int opp = util::opposite(iPop); + if (this->getBoundaryIntersection(iX, iY, opp, intersection) ) { + T vel[]= {T(),T()}; + u(vel,intersection); + this->defineU(iX, iY, opp, vel); + } + } + } + } + } + } +} + +template +void OffBoundaryConditionInstantiator2D:: +defineRho(int iX, int iY, int iPop, const T rho) +{ + this->getBlock().getDynamics(iX, iY)->defineRho(iPop, rho); + if (_output) { + clout << "defineRho(" << iX << ", " << iY << " )" << std::endl; + } +} + +template +void OffBoundaryConditionInstantiator2D:: +defineRho(BlockIndicatorF2D& indicator, + BlockIndicatorF2D& bulkIndicator, + AnalyticalF2D& rho) +{ + if ( indicator.isEmpty() ) { + return; + } + + const Vector min = indicator.getMin(); + const Vector max = indicator.getMax(); + + for (int iX = min[0]; iX <= max[0]; ++iX) { + for (int iY = min[1]; iY <= max[1]; ++iY) { + if (indicator(iX, iY)) { + for (int iPop = 1; iPop < DESCRIPTOR::q ; ++iPop) { + // Get direction + int iXn = iX + descriptors::c(iPop,0); + int iYn = iY + descriptors::c(iPop,1); + if (bulkIndicator(iXn, iYn)) { + T intersection[] = { T(), T() }; // coord. of intersection + int opp = util::opposite(iPop); + if (this->getBoundaryIntersection(iX, iY, opp, intersection) ) { + T rhoLocal[]= {T(1)}; + rho(rhoLocal,intersection); + this->defineRho(iX, iY, opp, rhoLocal[0]); + } + } + } + } + } + } +} + + +template +void OffBoundaryConditionInstantiator2D::outputOn() +{ + _output = true; +} + +template +void OffBoundaryConditionInstantiator2D::outputOff() +{ + _output = false; +} + +} + +#endif -- cgit v1.2.3