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