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/boundaryInstantiator3D.h | 1535 +++++++++++++++++++++++++++++++++ 1 file changed, 1535 insertions(+) create mode 100644 src/boundary/boundaryInstantiator3D.h (limited to 'src/boundary/boundaryInstantiator3D.h') diff --git a/src/boundary/boundaryInstantiator3D.h b/src/boundary/boundaryInstantiator3D.h new file mode 100644 index 0000000..4e40803 --- /dev/null +++ b/src/boundary/boundaryInstantiator3D.h @@ -0,0 +1,1535 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2007 Jonas Latt + * 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 3D boundaries -- header file. + */ + +#ifndef BOUNDARY_INSTANTIATOR_3D_H +#define BOUNDARY_INSTANTIATOR_3D_H + +#include "boundaryCondition3D.h" +#include "boundaryPostProcessors3D.h" +#include "wallFunctionBoundaryPostProcessors3D.h" +#include "io/ostreamManager.h" +#include "dynamics/wallFunctionLatticeDescriptors.h" +#include "functors/lattice/indicator/blockIndicatorF3D.h" +#include "dynamics/freeEnergyDynamics.h" + +namespace olb { + + +template +class BoundaryConditionInstantiator3D : public OnLatticeBoundaryCondition3D { +public: + BoundaryConditionInstantiator3D( BlockLatticeStructure3D& block_ ); + ~BoundaryConditionInstantiator3D() override; + + void addVelocityBoundary0N(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addVelocityBoundary0P(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addVelocityBoundary1N(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addVelocityBoundary1P(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addVelocityBoundary2N(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addVelocityBoundary2P(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + + void addSlipBoundary(int x0, int x1, int y0, int y1, int z0, int z1, int discreteNormalX, int discreteNormalY, int discreteNormalZ); + void addPartialSlipBoundary(T tuner, int x0, int x1, int y0, int y1, int z0, int z1, int discreteNormalX, int discreteNormalY, int discreteNormalZ); + + void addWallFunctionBoundary(int x0, int x1, int y0, int y1, int z0, int z1, BlockGeometryStructure3D& blockGeometryStructure, + std::vector discreteNormal, std::vector missingIndices, + UnitConverter const& converter, wallFunctionParam const& wallFunctionParam, + IndicatorF3D* geoIndicator); + + void addPressureBoundary0N(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addPressureBoundary0P(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addPressureBoundary1N(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addPressureBoundary1P(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addPressureBoundary2N(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addPressureBoundary2P(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + + void addConvectionBoundary0N(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv=NULL) override; + void addConvectionBoundary0P(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv=NULL) override; + void addConvectionBoundary1N(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv=NULL) override; + void addConvectionBoundary1P(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv=NULL) override; + void addConvectionBoundary2N(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv=NULL) override; + void addConvectionBoundary2P(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv=NULL) override; + + void addExternalVelocityEdge0NN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addExternalVelocityEdge0NP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addExternalVelocityEdge0PN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addExternalVelocityEdge0PP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addExternalVelocityEdge1NN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addExternalVelocityEdge1NP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addExternalVelocityEdge1PN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addExternalVelocityEdge1PP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addExternalVelocityEdge2NN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addExternalVelocityEdge2NP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addExternalVelocityEdge2PN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addExternalVelocityEdge2PP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + + void addInternalVelocityEdge0NN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addInternalVelocityEdge0NP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addInternalVelocityEdge0PN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addInternalVelocityEdge0PP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addInternalVelocityEdge1NN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addInternalVelocityEdge1NP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addInternalVelocityEdge1PN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addInternalVelocityEdge1PP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addInternalVelocityEdge2NN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addInternalVelocityEdge2NP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addInternalVelocityEdge2PN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + void addInternalVelocityEdge2PP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) override; + + void addExternalVelocityCornerNNN(int x, int y, int z, T omega) override; + void addExternalVelocityCornerNNP(int x, int y, int z, T omega) override; + void addExternalVelocityCornerNPN(int x, int y, int z, T omega) override; + void addExternalVelocityCornerNPP(int x, int y, int z, T omega) override; + void addExternalVelocityCornerPNN(int x, int y, int z, T omega) override; + void addExternalVelocityCornerPNP(int x, int y, int z, T omega) override; + void addExternalVelocityCornerPPN(int x, int y, int z, T omega) override; + void addExternalVelocityCornerPPP(int x, int y, int z, T omega) override; + + void addInternalVelocityCornerNNN(int x, int y, int z, T omega) override; + void addInternalVelocityCornerNNP(int x, int y, int z, T omega) override; + void addInternalVelocityCornerNPN(int x, int y, int z, T omega) override; + void addInternalVelocityCornerNPP(int x, int y, int z, T omega) override; + void addInternalVelocityCornerPNN(int x, int y, int z, T omega) override; + void addInternalVelocityCornerPNP(int x, int y, int z, T omega) override; + void addInternalVelocityCornerPPN(int x, int y, int z, T omega) override; + void addInternalVelocityCornerPPP(int x, int y, int z, T omega) override; + + void addVelocityBoundary(BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1, + T omega) override; + void addSlipBoundary(BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1) override; + void addPartialSlipBoundary(T tuner, BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1) override; + void addPressureBoundary(BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1, + T omega) override; + void addConvectionBoundary(BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1, + T omega, T* uAv=NULL) override; + void addWallFunctionBoundary(BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1, + UnitConverter const& converter, + wallFunctionParam const& wallFunctionParam, + IndicatorF3D* geoIndicator=NULL) override; + void addFreeEnergyWallBoundary(BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1, + T addend, int latticeNumber) override; + void addFreeEnergyInletBoundary(BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1, + T omega, std::string type, int latticeNumber) override; + void addFreeEnergyOutletBoundary(BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1, + T omega, std::string type, int latticeNumber) override; + + void outputOn() override; + void outputOff() override; + +private: + template + void addVelocityBoundary(int x0, int x1, int y0, int y1, int z0, int z1, T omega); + template + void addPressureBoundary(int x0, int x1, int y0, int y1, int z0, int z1, T omega); + template + void addConvectionBoundary(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv=NULL); + template + void addExternalVelocityEdge(int x0, int x1, int y0, int y1, int z0, int z1, T omega); + template + void addInternalVelocityEdge(int x0, int x1, int y0, int y1, int z0, int z1, T omega); + template + void addExternalVelocityCorner(int x, int y, int z, T omega); + template + void addInternalVelocityCorner(int x, int y, int z, T omega); +private: + BlockLatticeStructure3D& _block; + std::vector*> momentaVector; + std::vector*> dynamicsVector; + bool _output; + mutable OstreamManager clout; +}; + + +///////// class BoundaryConditionInstantiator3D //////////////////////// + +template +BoundaryConditionInstantiator3D::BoundaryConditionInstantiator3D ( + BlockLatticeStructure3D& block) + : _block(block), _output(false), clout(std::cout,"BoundaryConditionInstantiator3D") +{ } + +template +BoundaryConditionInstantiator3D::~BoundaryConditionInstantiator3D() +{ + for (auto &iDynamics : dynamicsVector) { + delete iDynamics; + } + for (auto &iMomenta : momentaVector) { + delete iMomenta; + } +} + +// Velocity BC + +template +template +void BoundaryConditionInstantiator3D:: +addVelocityBoundary(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + OLB_PRECONDITION( x0==x1 || y0==y1 || z0==z1 ); + + for (int iX=x0; iX<=x1; ++iX) { + for (int iY=y0; iY<=y1; ++iY) { + for (int iZ=z0; iZ<=z1; ++iZ) { + Momenta* momenta + = BoundaryManager::template getVelocityBoundaryMomenta(); + Dynamics* dynamics + = BoundaryManager::template getVelocityBoundaryDynamics(omega, *momenta); + _block.defineDynamics(iX,iX,iY,iY,iZ,iZ, dynamics); + momentaVector.push_back(momenta); + dynamicsVector.push_back(dynamics); + if (_output) { + clout << "addVelocityBoundary<" << direction << ", " << orientation << ">(" << iX << ", " << iX << ", "<< iY << ", " << iY << ", " << iZ << ", " << iZ << ", "<< omega << " )" << std::endl; + } + } + } + } + + PostProcessorGenerator3D* postProcessor + = BoundaryManager::template getVelocityBoundaryProcessor(x0,x1, y0,y1, z0,z1); + if (postProcessor) { + _block.addPostProcessor(*postProcessor); + } +} + +// Slip BC + +template +void BoundaryConditionInstantiator3D::addSlipBoundary( + int x0, int x1, int y0, int y1, int z0, int z1, int discreteNormalX, int discreteNormalY, int discreteNormalZ) +{ + OLB_PRECONDITION(x0==x1 || y0==y1 || z0==z1); + + for (int iX = x0; iX <= x1; ++iX) { + for (int iY = y0; iY <= y1; ++iY) { + for (int iZ = z0; iZ <= z1; ++iZ) { + if (_output) { + clout << "addSlipBoundary<" << discreteNormalX << ","<< discreteNormalY << ","<< discreteNormalZ << ">(" << x0 << ", "<< x1 << ", " << y0 << ", " << y1 << ", " << z0 << ", " << z1 << " )" << std::endl; + } + } + } + } + + PostProcessorGenerator3D* postProcessor = new SlipBoundaryProcessorGenerator3D(x0, x1, y0, y1, z0, z1, discreteNormalX, discreteNormalY, discreteNormalZ); + if (postProcessor) { + _block.addPostProcessor(*postProcessor); + } +} + +// Partial slip BC + +template +void BoundaryConditionInstantiator3D::addPartialSlipBoundary( + T tuner, int x0, int x1, int y0, int y1, int z0, int z1, int discreteNormalX, int discreteNormalY, int discreteNormalZ) +{ + OLB_PRECONDITION(x0==x1 || y0==y1 || z0==z1); + + for (int iX = x0; iX <= x1; ++iX) { + for (int iY = y0; iY <= y1; ++iY) { + for (int iZ = z0; iZ <= z1; ++iZ) { + if (_output) { + clout << "addPartialSlipBoundary<" << discreteNormalX << ","<< discreteNormalY << ","<< discreteNormalZ << ">(" << x0 << ", "<< x1 << ", " << y0 << ", " << y1 << ", " << z0 << ", " << z1 << " )" << std::endl; + } + } + } + } + + PostProcessorGenerator3D* postProcessor = new PartialSlipBoundaryProcessorGenerator3D(tuner, x0, x1, y0, y1, z0, z1, discreteNormalX, discreteNormalY, discreteNormalZ); + if (postProcessor) { + _block.addPostProcessor(*postProcessor); + } +} + +// Wall Function BC + +namespace { + + +template +struct WallFunctionBoundaryProcessorGenerator3DInstatiator { + static PostProcessorGenerator3D* getWallFunctionBoundaryProcessorGenerator(int x0, int x1, int y0, int y1, int z0, int z1, BlockGeometryStructure3D&, std::vector& discreteNormal, const std::vector& missingIndices, UnitConverter const&, wallFunctionParam const& wallFunctionParam, IndicatorF3D* geoIndicator) + { + std::abort(); + } +}; + +template +struct WallFunctionBoundaryProcessorGenerator3DInstatiator< + T, DESCRIPTOR, BoundaryManager, + typename std::enable_if< + DESCRIPTOR::template provides() && + DESCRIPTOR::template provides() && + DESCRIPTOR::template provides() && + DESCRIPTOR::template provides() && + DESCRIPTOR::template provides() + >::type +> { + static PostProcessorGenerator3D* getWallFunctionBoundaryProcessorGenerator(int x0, int x1, int y0, int y1, int z0, int z1, BlockGeometryStructure3D& blockGeometryStructure, std::vector& discreteNormal, const std::vector& missingIndices, UnitConverter const& converter, wallFunctionParam const& wallFunctionParam, IndicatorF3D* geoIndicator) + { + PostProcessorGenerator3D* postProcessor = new WallFunctionBoundaryProcessorGenerator3D(x0, x1, y0, y1, z0, z1, blockGeometryStructure, discreteNormal, missingIndices, + converter, wallFunctionParam, geoIndicator); + + return postProcessor; + } +}; + +} + +template +void BoundaryConditionInstantiator3D::addWallFunctionBoundary( + int x0, int x1, int y0, int y1, int z0, int z1, BlockGeometryStructure3D& blockGeometryStructure, + std::vector discreteNormal, std::vector missingIndices, + UnitConverter const& converter, wallFunctionParam const& wallFunctionParam, IndicatorF3D* geoIndicator) +{ + OLB_PRECONDITION(x0==x1 || y0==y1 || z0==z1); + + for (int iX = x0; iX <= x1; ++iX) { + for (int iY = y0; iY <= y1; ++iY) { + for (int iZ = z0; iZ <= z1; ++iZ) { + if (_output) { + clout << "addWallFunctionBoundary<" << discreteNormal[0] << ","<< discreteNormal[1] << ","<< discreteNormal[2] << ">(" << x0 << ", "<< x1 << ", " << y0 << ", " << y1 << ", " << z0 << ", " << z1 << " )" << std::endl; + } + } + } + } + + PostProcessorGenerator3D* postProcessor = WallFunctionBoundaryProcessorGenerator3DInstatiator::getWallFunctionBoundaryProcessorGenerator(x0, x1, y0, y1, z0, z1, blockGeometryStructure, discreteNormal, missingIndices, converter, wallFunctionParam, geoIndicator); + if (postProcessor) { + _block.addPostProcessor(*postProcessor); + } +} + +// Pressure BC + +template +template +void BoundaryConditionInstantiator3D:: +addPressureBoundary(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + OLB_PRECONDITION( x0==x1 || y0==y1 || z0==z1 ); + + for (int iX=x0; iX<=x1; ++iX) { + for (int iY=y0; iY<=y1; ++iY) { + for (int iZ=z0; iZ<=z1; ++iZ) { + Momenta* momenta + = BoundaryManager::template getPressureBoundaryMomenta(); + Dynamics* dynamics + = BoundaryManager::template getPressureBoundaryDynamics(omega, *momenta); + _block.defineDynamics(iX,iX,iY,iY,iZ,iZ, dynamics); + momentaVector.push_back(momenta); + dynamicsVector.push_back(dynamics); + if (_output) { + clout << "addPressureBoundary<" << direction << ", " << orientation << ">(" << iX << ", " << iX << ", "<< iY << ", " << iY << ", " << iZ << ", " << iZ << ", "<< omega << " )" << std::endl; + } + } + } + } + + PostProcessorGenerator3D* postProcessor + = BoundaryManager::template getPressureBoundaryProcessor(x0,x1, y0,y1, z0,z1); + if (postProcessor) { + _block.addPostProcessor(*postProcessor); + } +} + +// Convection BC + +template +template +void BoundaryConditionInstantiator3D:: +addConvectionBoundary(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv) +{ + OLB_PRECONDITION( x0==x1 || y0==y1 || z0==z1 ); + + for (int iX=x0; iX<=x1; ++iX) { + for (int iY=y0; iY<=y1; ++iY) { + for (int iZ=z0; iZ<=z1; ++iZ) { + if (_output) { + clout << "addConvectionBoundary<" << direction << ", " << orientation << ">(" << iX << ", " << iX << ", "<< iY << ", " << iY << ", " << iZ << ", " << iZ << ", "<< omega << " )" << std::endl; + } + } + } + } + + PostProcessorGenerator3D* postProcessor + = BoundaryManager::template getConvectionBoundaryProcessor(x0,x1, y0,y1, z0,z1, uAv); + if (postProcessor) { + _block.addPostProcessor(*postProcessor); + } +} + +template +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityEdge(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + OLB_PRECONDITION( + ( x0==x1 && y0==y1 ) || + ( x0==x1 && z0==z1 ) || + ( y0==y1 && z0==z1 ) ); + + for (int iX=x0; iX<=x1; ++iX) { + for (int iY=y0; iY<=y1; ++iY) { + for (int iZ=z0; iZ<=z1; ++iZ) { + Momenta* momenta + = BoundaryManager::template getExternalVelocityEdgeMomenta(); + Dynamics* dynamics + = BoundaryManager::template getExternalVelocityEdgeDynamics(omega, *momenta); + _block.defineDynamics(iX,iX,iY,iY,iZ,iZ, dynamics); + momentaVector.push_back(momenta); + dynamicsVector.push_back(dynamics); + if (_output) { + clout << "addExternalVelocityEdge<" << plane << ", " << normal1 << ", " << normal2 << ">(" << iX << ", " << iX << ", "<< iY << ", " << iY << ", " << iZ << ", " << iZ << ", "<< omega << " )" << std::endl; + } + } + } + } + + PostProcessorGenerator3D* postProcessor + = BoundaryManager::template getExternalVelocityEdgeProcessor(x0,x1, y0,y1, z0,z1); + if (postProcessor) { + _block.addPostProcessor(*postProcessor); + } +} + +template +template +void BoundaryConditionInstantiator3D:: +addInternalVelocityEdge(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + if (!(( x0==x1 && y0==y1 ) || + ( x0==x1 && z0==z1 ) || + ( y0==y1 && z0==z1 ) )) { + clout << x0 <<" "<< x1 <<" "<< y0 <<" "<< y1 <<" "<< z0 <<" "<< z1 << std::endl; + } + + OLB_PRECONDITION( + ( x0==x1 && y0==y1 ) || + ( x0==x1 && z0==z1 ) || + ( y0==y1 && z0==z1 ) ); + + for (int iX=x0; iX<=x1; ++iX) { + for (int iY=y0; iY<=y1; ++iY) { + for (int iZ=z0; iZ<=z1; ++iZ) { + Momenta* momenta + = BoundaryManager::template getInternalVelocityEdgeMomenta(); + Dynamics* dynamics + = BoundaryManager::template getInternalVelocityEdgeDynamics(omega, *momenta); + _block.defineDynamics(iX,iX,iY,iY,iZ,iZ, dynamics); + momentaVector.push_back(momenta); + dynamicsVector.push_back(dynamics); + if (_output) { + clout << "addInternalVelocityEdge<" << plane << ", " << normal1 << ", " << normal2 << ">(" << iX << ", " << iX << ", "<< iY << ", " << iY << ", " << iZ << ", " << iZ << ", "<< omega << " )" << std::endl; + } + } + } + } + + PostProcessorGenerator3D* postProcessor + = BoundaryManager::template getInternalVelocityEdgeProcessor(x0,x1, y0,y1, z0,z1); + if (postProcessor) { + _block.addPostProcessor(*postProcessor); + } +} + + +template +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityCorner(int x, int y, int z, T omega) +{ + Momenta* momenta + = BoundaryManager::template getExternalVelocityCornerMomenta(); + Dynamics* dynamics + = BoundaryManager::template getExternalVelocityCornerDynamics(omega, *momenta); + + _block.defineDynamics(x,x,y,y,z,z, dynamics); + + momentaVector.push_back(momenta); + dynamicsVector.push_back(dynamics); + + PostProcessorGenerator3D* postProcessor + = BoundaryManager::template getExternalVelocityCornerProcessor(x, y, z); + if (postProcessor) { + _block.addPostProcessor(*postProcessor); + } + if (_output) { + clout << "addExternalVelocityCorner<" << xNormal << ", " << yNormal << ", " << zNormal << ">(" << x << ", " << y << ", "<< z << omega << " )" << std::endl; + } +} + +template +template +void BoundaryConditionInstantiator3D:: +addInternalVelocityCorner(int x, int y, int z, T omega) +{ + Momenta* momenta + = BoundaryManager::template getInternalVelocityCornerMomenta(); + Dynamics* dynamics + = BoundaryManager::template getInternalVelocityCornerDynamics(omega, *momenta); + + _block.defineDynamics(x,x,y,y,z,z, dynamics); + + momentaVector.push_back(momenta); + dynamicsVector.push_back(dynamics); + + PostProcessorGenerator3D* postProcessor + = BoundaryManager::template getInternalVelocityCornerProcessor(x, y, z); + if (postProcessor) { + _block.addPostProcessor(*postProcessor); + } + if (_output) { + clout << "addInternalVelocityCorner<" << xNormal << ", " << yNormal << ", " << zNormal << ">(" << x << ", " << y << ", "<< z << omega << " )" << std::endl; + } +} + +// Velocity BC + +template +void BoundaryConditionInstantiator3D:: +addVelocityBoundary(BlockIndicatorF3D& indicator, int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + std::vector discreteNormal(4,0); + for (int iX = x0; iX <= x1; iX++) { + for (int iY = y0; iY <= y1; iY++) { + for (int iZ = z0; iZ <= z1; iZ++) { + if (indicator(iX, iY, iZ)) { + discreteNormal = indicator.getBlockGeometryStructure().getStatistics().getType(iX, iY, iZ); + if (discreteNormal[0] == 0) { + if (discreteNormal[1] != 0 && discreteNormal[1] == -1) { + addVelocityBoundary<0,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] != 0 && discreteNormal[1] == 1) { + addVelocityBoundary<0,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[2] != 0 && discreteNormal[2] == -1) { + addVelocityBoundary<1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[2] != 0 && discreteNormal[2] == 1) { + addVelocityBoundary<1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[3] != 0 && discreteNormal[3] == -1) { + addVelocityBoundary<2,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[3] != 0 && discreteNormal[3] == 1) { + addVelocityBoundary<2,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + } + + else if (discreteNormal[0] == 1) { + if (discreteNormal[1] == 1 && discreteNormal[2] == 1 && discreteNormal[3] == 1) { + addExternalVelocityCorner<1,1,1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == -1 && discreteNormal[3] == 1) { + addExternalVelocityCorner<1,-1,1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == 1 && discreteNormal[3] == -1) { + addExternalVelocityCorner<1,1,-1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == -1 && discreteNormal[3] == -1) { + addExternalVelocityCorner<1,-1,-1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == 1 && discreteNormal[3] == 1) { + addExternalVelocityCorner<-1,1,1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == -1 && discreteNormal[3] == 1) { + addExternalVelocityCorner<-1,-1,1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == 1 && discreteNormal[3] == -1) { + addExternalVelocityCorner<-1,1,-1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == -1 && discreteNormal[3] == -1) { + addExternalVelocityCorner<-1,-1,-1>(iX,iY,iZ, omega); + } + /// addExternalVelocityCorner(iX,iY,iZ, omega); + } + + else if (discreteNormal[0] == 2) { + if (discreteNormal[1] == 1 && discreteNormal[2] == 1 && discreteNormal[3] == 1) { + addInternalVelocityCorner<1,1,1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == -1 && discreteNormal[3] == 1) { + addExternalVelocityCorner<1,-1,1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == 1 && discreteNormal[3] == -1) { + addInternalVelocityCorner<1,1,-1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == -1 && discreteNormal[3] == -1) { + addInternalVelocityCorner<1,-1,-1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == 1 && discreteNormal[3] == 1) { + addInternalVelocityCorner<-1,1,1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == -1 && discreteNormal[3] == 1) { + addInternalVelocityCorner<-1,-1,1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == 1 && discreteNormal[3] == -1) { + addInternalVelocityCorner<-1,1,-1>(iX,iY,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == -1 && discreteNormal[3] == -1) { + addInternalVelocityCorner<-1,-1,-1>(iX,iY,iZ, omega); + } + /// addInternalVelocityCorner(iX,iY,iZ, omega); + } + + else if (discreteNormal[0] == 3) { + if (discreteNormal[1] == 0 && discreteNormal[2] == 1 && discreteNormal[3] == 1) { + addExternalVelocityEdge<0,1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 0 && discreteNormal[2] == -1 && discreteNormal[3] == 1) { + addExternalVelocityEdge<0,-1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 0 && discreteNormal[2] == 1 && discreteNormal[3] == -1) { + addExternalVelocityEdge<0,1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 0 && discreteNormal[2] == -1 && discreteNormal[3] == -1) { + addExternalVelocityEdge<0,-1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == 0 && discreteNormal[3] == 1) { + addExternalVelocityEdge<1,1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == 0 && discreteNormal[3] == 1) { + addExternalVelocityEdge<1,1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == 0 && discreteNormal[3] == -1) { + addExternalVelocityEdge<1,-1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == 0 && discreteNormal[3] == -1) { + addExternalVelocityEdge<1,-1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == 1 && discreteNormal[3] == 0) { + addExternalVelocityEdge<2,1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == 1 && discreteNormal[3] == 0) { + addExternalVelocityEdge<2,-1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == -1 && discreteNormal[3] == 0) { + addExternalVelocityEdge<2,1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == -1 && discreteNormal[3] == 0) { + addExternalVelocityEdge<2,-1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + } + + else if (discreteNormal[0] == 4) { + if (discreteNormal[1] == 0 && discreteNormal[2] == 1 && discreteNormal[3] == 1) { + addInternalVelocityEdge<0,1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 0 && discreteNormal[2] == -1 && discreteNormal[3] == 1) { + addInternalVelocityEdge<0,-1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 0 && discreteNormal[2] == 1 && discreteNormal[3] == -1) { + addInternalVelocityEdge<0,1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 0 && discreteNormal[2] == -1 && discreteNormal[3] == -1) { + addInternalVelocityEdge<0,-1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == 0 && discreteNormal[3] == 1) { + addInternalVelocityEdge<1,1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == 0 && discreteNormal[3] == 1) { + addInternalVelocityEdge<1,1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == 0 && discreteNormal[3] == -1) { + addInternalVelocityEdge<1,-1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == 0 && discreteNormal[3] == -1) { + addInternalVelocityEdge<1,-1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == 1 && discreteNormal[3] == 0) { + addInternalVelocityEdge<2,1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == 1 && discreteNormal[3] == 0) { + addInternalVelocityEdge<2,-1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == 1 && discreteNormal[2] == -1 && discreteNormal[3] == 0) { + addInternalVelocityEdge<2,1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] == -1 && discreteNormal[2] == -1 && discreteNormal[3] == 0) { + addInternalVelocityEdge<2,-1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + } + } + } + } + } +} + +// Slip BC + +template +void BoundaryConditionInstantiator3D::addSlipBoundary( + BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1) +{ + auto& blockGeometryStructure = indicator.getBlockGeometryStructure(); + std::vector discreteNormal(4, 0); + for (int iX = x0; iX <= x1; iX++) { + for (int iY = y0; iY <= y1; iY++) { + for (int iZ = z0; iZ <= z1; iZ++) { + if (indicator(iX, iY, iZ)) { + discreteNormal = blockGeometryStructure.getStatistics().getType(iX, iY, iZ); + if (discreteNormal[1]!=0 || discreteNormal[2]!=0 || discreteNormal[3]!=0) { + addSlipBoundary(iX, iX, iY, iY, iZ, iZ, + discreteNormal[1], discreteNormal[2], discreteNormal[3]); + } + else { + clout << "Warning: Could not addSlipBoundary (" << iX << ", " << iY << ", " << iZ << "), discreteNormal=(" << discreteNormal[0] <<","<< discreteNormal[1] <<","<< discreteNormal[2] <<","<< discreteNormal[3] <<"), set to bounceBack" << std::endl; + _block.defineDynamics(iX, iY, iZ, &instances::getBounceBack()); + } + } + } + } + } +} + +// Partial slip BC + +template +void BoundaryConditionInstantiator3D::addPartialSlipBoundary( + T tuner, BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1) +{ + auto& blockGeometryStructure = indicator.getBlockGeometryStructure(); + std::vector discreteNormal(4, 0); + for (int iX = x0; iX <= x1; iX++) { + for (int iY = y0; iY <= y1; iY++) { + for (int iZ = z0; iZ <= z1; iZ++) { + if (indicator(iX, iY, iZ)) { + if (tuner < 0. || tuner > 1.) { + clout << "Warning: Could not addPartialSlipBoundary (" << iX << ", " << iY << ", " << iZ << "), tuner must be between 0.1 and instead is=" << tuner <<", set to bounceBack" << std::endl; + _block.defineDynamics(iX, iY, iZ, &instances::getBounceBack()); + } else { + discreteNormal = blockGeometryStructure.getStatistics().getType(iX, iY, iZ); + if (discreteNormal[1]!=0 || discreteNormal[2]!=0 || discreteNormal[3]!=0) { + addPartialSlipBoundary(tuner, iX, iX, iY, iY, iZ, iZ, + discreteNormal[1], discreteNormal[2], discreteNormal[3]); + } + else { + clout << "Warning: Could not addPartialSlipBoundary (" << iX << ", " << iY << ", " << iZ << "), discreteNormal=(" << discreteNormal[0] <<","<< discreteNormal[1] <<","<< discreteNormal[2] <<","<< discreteNormal[3] <<"), set to bounceBack" << std::endl; + _block.defineDynamics(iX, iY, iZ, &instances::getBounceBack()); + } + } + } + } + } + } +} + +// Pressure BC + +template +void BoundaryConditionInstantiator3D:: +addPressureBoundary(BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + auto& blockGeometryStructure = indicator.getBlockGeometryStructure(); + std::vector discreteNormal(4,0); + for (int iX = x0; iX <= x1; ++iX) { + for (int iY = y0; iY <= y1; ++iY) { + for (int iZ = z0; iZ <= z1; ++iZ) { + if (indicator(iX, iY, iZ)) { + discreteNormal = blockGeometryStructure.getStatistics().getType(iX, iY, iZ); + + if (discreteNormal[0] == 0) { + if (discreteNormal[1] != 0 && discreteNormal[1] == -1) { + addPressureBoundary<0,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[1] != 0 && discreteNormal[1] == 1) { + addPressureBoundary<0,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[2] != 0 && discreteNormal[2] == -1) { + addPressureBoundary<1,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[2] != 0 && discreteNormal[2] == 1) { + addPressureBoundary<1,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[3] != 0 && discreteNormal[3] == -1) { + addPressureBoundary<2,-1>(iX,iX,iY,iY,iZ,iZ, omega); + } + else if (discreteNormal[3] != 0 && discreteNormal[3] == 1) { + addPressureBoundary<2,1>(iX,iX,iY,iY,iZ,iZ, omega); + } + } + } + } + } + } +} + +// Convection BC + +template +void BoundaryConditionInstantiator3D:: +addConvectionBoundary(BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1, + T omega, T* uAv) +{ + auto& blockGeometryStructure = indicator.getBlockGeometryStructure(); + std::vector discreteNormal(4,0); + for (int iX = x0; iX <= x1; ++iX) { + for (int iY = y0; iY <= y1; ++iY) { + for (int iZ = z0; iZ <= z1; ++iZ) { + if (indicator(iX, iY, iZ)) { + discreteNormal = blockGeometryStructure.getStatistics().getType(iX, iY, iZ); + + if (discreteNormal[0] == 0) { + if (discreteNormal[1] != 0 && discreteNormal[1] == -1) { + addConvectionBoundary<0,-1>(iX,iX,iY,iY,iZ,iZ, omega, uAv); + } + else if (discreteNormal[1] != 0 && discreteNormal[1] == 1) { + addConvectionBoundary<0,1>(iX,iX,iY,iY,iZ,iZ, omega, uAv); + } + else if (discreteNormal[2] != 0 && discreteNormal[2] == -1) { + addConvectionBoundary<1,-1>(iX,iX,iY,iY,iZ,iZ, omega, uAv); + } + else if (discreteNormal[2] != 0 && discreteNormal[2] == 1) { + addConvectionBoundary<1,1>(iX,iX,iY,iY,iZ,iZ, omega, uAv); + } + + else if (discreteNormal[3] != 0 && discreteNormal[3] == -1) { + addConvectionBoundary<2,-1>(iX,iX,iY,iY,iZ,iZ, omega, uAv); + } + else if (discreteNormal[3] != 0 && discreteNormal[3] == 1) { + addConvectionBoundary<2,1>(iX,iX,iY,iY,iZ,iZ, omega, uAv); + } + } + } + } + } + } +} + +// Wall Function BC + +template +void BoundaryConditionInstantiator3D::addWallFunctionBoundary( + BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1, + UnitConverter const& converter, + wallFunctionParam const& wallFunctionParam, + IndicatorF3D* geoIndicator) +{ + const auto& blockGeometryStructure = indicator.getBlockGeometryStructure(); + std::vector discreteNormal(4, 0); + for (int iX = x0; iX <= x1; ++iX) { + for (int iY = y0; iY <= y1; ++iY) { + for (int iZ = z0; iZ <= z1; ++iZ) { + if (indicator(iX, iY, iZ)) { + discreteNormal = indicator.getBlockGeometryStructure().getStatistics().getType(iX, iY, iZ); + std::vector missingIndices; + for (int x = -1 ; x < 2; ++x) { + for (int y = -1 ; y < 2; ++y) { + for (int z = -1 ; z < 2; ++z) { + if (blockGeometryStructure.getMaterial(iX + x, iY + y, iZ + z) == 0) { + for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) { + if (descriptors::c(iPop,0) == x && + descriptors::c(iPop,1) == y && + descriptors::c(iPop,2) == z) { + missingIndices.push_back(descriptors::opposite(iPop)); + } + } + } + } + } + } + if (discreteNormal[1]!=0 || discreteNormal[2]!=0 || discreteNormal[3]!=0) { + discreteNormal.erase(discreteNormal.begin()); + addWallFunctionBoundary(iX, iX, iY, iY, iZ, iZ, + indicator.getBlockGeometryStructure(), + discreteNormal, missingIndices, + converter, wallFunctionParam, geoIndicator); + } + else { + //clout << "Warning: Could not add WallFunction (" << iX << ", " << iY << ", " << iZ << "), discreteNormal=(" << discreteNormal[0] <<","<< discreteNormal[1] <<","<< discreteNormal[2] <<","<< discreteNormal[3] <<"), set to bounceBack" << std::endl; + _block.defineDynamics(iX, iY, iZ, &instances::getBounceBack()); + } + } + } + } + } +} + +// Free Energy BC + +template +void BoundaryConditionInstantiator3D::addFreeEnergyWallBoundary( + BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1, T addend, int latticeNumber) +{ + auto& blockGeometryStructure = indicator.getBlockGeometryStructure(); + std::vector discreteNormal(4, 0); + for (int iX = x0; iX <= x1; ++iX) { + for (int iY = y0; iY <= y1; ++iY) { + for (int iZ = z0; iZ <= z1; ++iZ) { + if(indicator(iX,iY,iZ)) { + discreteNormal = blockGeometryStructure.getStatistics().getType(iX, iY, iZ, true); + if (discreteNormal[1]!=0 || discreteNormal[2]!=0 || discreteNormal[3]!=0) { + + Dynamics* dynamics = NULL; + if (latticeNumber == 1) { + dynamics = &instances::getBounceBack(); + } else { + dynamics = new FreeEnergyWallDynamics; + dynamicsVector.push_back(dynamics); + } + _block.get(iX,iY,iZ).defineDynamics(dynamics); + + PostProcessorGenerator3D* wettingPostProcessor = + new FreeEnergyWallProcessorGenerator3D ( iX, iX, iY, iY, iZ, iZ, + discreteNormal[1], discreteNormal[2], discreteNormal[3], addend ); + PostProcessorGenerator3D* chemPotPostProcessor = + new FreeEnergyChemPotBoundaryProcessorGenerator3D ( iX, iX, iY, iY, iZ, iZ, + discreteNormal[1], discreteNormal[2], discreteNormal[3], latticeNumber ); + if (wettingPostProcessor) { + _block.addPostProcessor(*wettingPostProcessor); + } + if (chemPotPostProcessor) { + _block.addPostProcessor(*chemPotPostProcessor); + } + } + if (_output) { + clout << "addFreeEnergyWallBoundary<" << "," << ">(" << x0 << ", "<< x1 << ", " << y0 << ", " << y1 << ", " << z0 << ", " << z1 << ")" << std::endl; + } + } + } + } + } + +} + +template +void BoundaryConditionInstantiator3D::addFreeEnergyInletBoundary( + BlockIndicatorF3D& indicator, + int x0, int x1, int y0, int y1, int z0, int z1, + T omega, std::string type, int latticeNumber) +{ + auto& blockGeometryStructure = indicator.getBlockGeometryStructure(); + std::vector discreteNormal(4, 0); + for (int iX = x0; iX <= x1; ++iX) { + for (int iY = y0; iY <= y1; ++iY) { + for (int iZ = z0; iZ <= z1; ++iZ) { + if(indicator(iX,iY,iZ)) { + discreteNormal = blockGeometryStructure.getStatistics().getType(iX, iY, iZ, true); + + if (discreteNormal[0] == 0) { + Momenta* momenta = NULL; + Dynamics* dynamics = NULL; + + if (discreteNormal[1] != 0 && discreteNormal[1] == -1) { + if (latticeNumber == 1) { + if (type == "density") { + addPressureBoundary<0,-1>(iX, iX, iY, iY, iZ, iZ, omega); + } else { + addVelocityBoundary<0,-1>(iX, iX, iY, iY, iZ, iZ, omega); + } + } else { + momenta = BoundaryManager::template + getPressureBoundaryMomenta<0,-1>(); + dynamics = new FreeEnergyInletOutletDynamics(omega,*momenta); + } + } + + else if (discreteNormal[1] != 0 && discreteNormal[1] == 1) { + if (latticeNumber == 1) { + if (type == "density") { + addPressureBoundary<0,1>(iX, iX, iY, iY, iZ, iZ, omega); + } else { + addVelocityBoundary<0,1>(iX, iX, iY, iY, iZ, iZ, omega); + } + } else { + momenta = BoundaryManager::template + getPressureBoundaryMomenta<0,1>(); + dynamics = new FreeEnergyInletOutletDynamics(omega,*momenta); + } + } + + else if (discreteNormal[2] != 0 && discreteNormal[2] == -1) { + if (latticeNumber == 1) { + if (type == "density") { + addPressureBoundary<1,-1>(iX, iX, iY, iY, iZ, iZ, omega); + } else { + addVelocityBoundary<1,-1>(iX, iX, iY, iY, iZ, iZ, omega); + } + } else { + momenta = BoundaryManager::template + getPressureBoundaryMomenta<1,-1>(); + dynamics = new FreeEnergyInletOutletDynamics(omega,*momenta); + } + } + + else if (discreteNormal[2] != 0 && discreteNormal[2] == 1) { + if (latticeNumber == 1) { + if (type == "density") { + addPressureBoundary<1,1>(iX, iX, iY, iY, iZ, iZ, omega); + } else { + addVelocityBoundary<1,1>(iX, iX, iY, iY, iZ, iZ, omega); + } + } else { + momenta = BoundaryManager::template + getPressureBoundaryMomenta<1,1>(); + dynamics = new FreeEnergyInletOutletDynamics(omega,*momenta); + } + momenta = BoundaryManager::template + getPressureBoundaryMomenta<1,1>(); + if (latticeNumber == 1) { + dynamics = BoundaryManager::template + getPressureBoundaryDynamics<1,1>(omega, *momenta); + } else { + dynamics = new FreeEnergyInletOutletDynamics(omega,*momenta); + } + } + + else if (discreteNormal[3] != 0 && discreteNormal[3] == -1) { + if (latticeNumber == 1) { + if (type == "density") { + addPressureBoundary<2,-1>(iX, iX, iY, iY, iZ, iZ, omega); + } else { + addVelocityBoundary<2,-1>(iX, iX, iY, iY, iZ, iZ, omega); + } + } else { + momenta = BoundaryManager::template + getPressureBoundaryMomenta<2,-1>(); + dynamics = new FreeEnergyInletOutletDynamics(omega,*momenta); + } + } + + else if (discreteNormal[3] != 0 && discreteNormal[3] == 1) { + if (latticeNumber == 1) { + if (type == "density") { + addPressureBoundary<2,1>(iX, iX, iY, iY, iZ, iZ, omega); + } else { + addVelocityBoundary<2,1>(iX, iX, iY, iY, iZ, iZ, omega); + } + } else { + momenta = BoundaryManager::template + getPressureBoundaryMomenta<2,1>(); + dynamics = new FreeEnergyInletOutletDynamics(omega,*momenta); + } + } + + if (latticeNumber != 1) { + _block.defineDynamics(iX,iX,iY,iY,iZ,iZ, dynamics); + momentaVector.push_back(momenta); + dynamicsVector.push_back(dynamics); + } + } + + PostProcessorGenerator3D* chemPotPostProcessor = + new FreeEnergyChemPotBoundaryProcessorGenerator3D ( iX, iX, iY, iY, iZ, iZ, + discreteNormal[1], discreteNormal[2], discreteNormal[3], latticeNumber ); + if (chemPotPostProcessor) { + _block.addPostProcessor(*chemPotPostProcessor); + } + + if (_output) { + clout << "addFreeEnergyInletBoundary<" << "," << ">(" << x0 << ", "<< x1 << ", " << y0 << ", " << y1 << ", " << z0 << ", " << z1 << ")" << std::endl; + } + } + } + } + } + +} + +template +void BoundaryConditionInstantiator3D::addFreeEnergyOutletBoundary( + BlockIndicatorF3D& indicator, int x0, int x1, int y0, int y1, int z0, int z1, + T omega, std::string type, int latticeNumber) +{ + addFreeEnergyInletBoundary(indicator, x0, x1, y0, y1, z0, z1, omega, type, latticeNumber); + + auto& blockGeometryStructure = indicator.getBlockGeometryStructure(); + std::vector discreteNormal(4, 0); + for (int iX = x0; iX <= x1; ++iX) { + for (int iY = y0; iY <= y1; ++iY) { + for (int iZ = z0; iZ <= z1; ++iZ) { + if(indicator(iX,iY,iZ)) { + discreteNormal = blockGeometryStructure.getStatistics().getType(iX, iY, iZ, true); + if (discreteNormal[0] == 0) { + + PostProcessorGenerator3D* convectivePostProcessor = + new FreeEnergyConvectiveProcessorGenerator3D ( + iX, iX, iY, iY, iZ, iZ, discreteNormal[1], discreteNormal[2], discreteNormal[3] ); + if (convectivePostProcessor) { + _block.addPostProcessor(*convectivePostProcessor); + } + + if (_output) { + clout << "addFreeEnergyOutletBoundary<" << "," << ">(" << x0 << ", "<< x1 << ", " << y0 << ", " << y1 << ", " << z0 << ", " << z1 << ")" << std::endl; + } + } + } + } + } + } + +} + + +template +void BoundaryConditionInstantiator3D:: +addVelocityBoundary0N(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addVelocityBoundary<0,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addVelocityBoundary0P(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addVelocityBoundary<0,1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addVelocityBoundary1N(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addVelocityBoundary<1,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addVelocityBoundary1P(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addVelocityBoundary<1, 1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addVelocityBoundary2N(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addVelocityBoundary<2,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addVelocityBoundary2P(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addVelocityBoundary<2, 1>(x0,x1,y0,y1,z0,z1, omega); +} + +// Pressure BC + +template +void BoundaryConditionInstantiator3D:: +addPressureBoundary0N(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addPressureBoundary<0,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addPressureBoundary0P(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addPressureBoundary<0,1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addPressureBoundary1N(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addPressureBoundary<1,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addPressureBoundary1P(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addPressureBoundary<1, 1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addPressureBoundary2N(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addPressureBoundary<2,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addPressureBoundary2P(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addPressureBoundary<2, 1>(x0,x1,y0,y1,z0,z1, omega); +} + +// Convection BC + +template +void BoundaryConditionInstantiator3D:: +addConvectionBoundary0N(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv) +{ + addConvectionBoundary<0,-1>(x0,x1,y0,y1,z0,z1, omega, uAv); +} + +template +void BoundaryConditionInstantiator3D:: +addConvectionBoundary0P(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv) +{ + addConvectionBoundary<0,1>(x0,x1,y0,y1,z0,z1, omega, uAv); +} + +template +void BoundaryConditionInstantiator3D:: +addConvectionBoundary1N(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv) +{ + addConvectionBoundary<1,-1>(x0,x1,y0,y1,z0,z1, omega, uAv); +} + +template +void BoundaryConditionInstantiator3D:: +addConvectionBoundary1P(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv) +{ + addConvectionBoundary<1, 1>(x0,x1,y0,y1,z0,z1, omega, uAv); +} + +template +void BoundaryConditionInstantiator3D:: +addConvectionBoundary2N(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv) +{ + addConvectionBoundary<2,-1>(x0,x1,y0,y1,z0,z1, omega, uAv); +} + +template +void BoundaryConditionInstantiator3D:: +addConvectionBoundary2P(int x0, int x1, int y0, int y1, int z0, int z1, T omega, T* uAv) +{ + addConvectionBoundary<2, 1>(x0,x1,y0,y1,z0,z1, omega, uAv); +} + + +// Velocity BC + +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityEdge0NN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addExternalVelocityEdge<0,-1,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityEdge0NP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addExternalVelocityEdge<0,-1, 1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityEdge0PN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addExternalVelocityEdge<0, 1,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityEdge0PP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addExternalVelocityEdge<0, 1, 1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityEdge1NN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addExternalVelocityEdge<1,-1,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityEdge1NP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addExternalVelocityEdge<1,-1, 1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityEdge1PN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addExternalVelocityEdge<1, 1,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityEdge1PP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addExternalVelocityEdge<1, 1, 1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityEdge2NN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addExternalVelocityEdge<2,-1,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityEdge2NP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addExternalVelocityEdge<2,-1, 1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityEdge2PN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addExternalVelocityEdge<2, 1,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addExternalVelocityEdge2PP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addExternalVelocityEdge<2, 1, 1>(x0,x1,y0,y1,z0,z1, omega); +} + + + +template +void BoundaryConditionInstantiator3D:: +addInternalVelocityEdge0NN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addInternalVelocityEdge<0,-1,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addInternalVelocityEdge0NP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addInternalVelocityEdge<0,-1, 1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addInternalVelocityEdge0PN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addInternalVelocityEdge<0, 1,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addInternalVelocityEdge0PP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addInternalVelocityEdge<0, 1, 1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addInternalVelocityEdge1NN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addInternalVelocityEdge<1,-1,-1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addInternalVelocityEdge1NP(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addInternalVelocityEdge<1,-1, 1>(x0,x1,y0,y1,z0,z1, omega); +} + +template +void BoundaryConditionInstantiator3D:: +addInternalVelocityEdge1PN(int x0, int x1, int y0, int y1, int z0, int z1, T omega) +{ + addInternalVelocityEdge<1, 1