From 82284392ef09b97d274234d1cd76c02238e4dabb Mon Sep 17 00:00:00 2001
From: kutscher <kutscher@irmb.tu-bs.de>
Date: Mon, 18 Sep 2017 14:49:17 +0200
Subject: [PATCH] - deletes depot folder
---
authors.txt | 2 +
depot/MeanValuesCoProcessor.cpp | 67 --
depot/MeanValuesCoProcessor.h | 26 -
depot/Particles.cpp | 11 -
depot/Particles.h | 29 -
depot/PathLineCoProcessor.cpp | 1505 ---------------------------
depot/PathLineCoProcessor.h | 115 --
depot/PathLineCoProcessorMcpart.cpp | 1301 -----------------------
depot/PathLineCoProcessorMcpart.h | 133 ---
9 files changed, 2 insertions(+), 3187 deletions(-)
create mode 100644 authors.txt
delete mode 100644 depot/MeanValuesCoProcessor.cpp
delete mode 100644 depot/MeanValuesCoProcessor.h
delete mode 100644 depot/Particles.cpp
delete mode 100644 depot/Particles.h
delete mode 100644 depot/PathLineCoProcessor.cpp
delete mode 100644 depot/PathLineCoProcessor.h
delete mode 100644 depot/PathLineCoProcessorMcpart.cpp
delete mode 100644 depot/PathLineCoProcessorMcpart.h
diff --git a/authors.txt b/authors.txt
new file mode 100644
index 000000000..45531bf38
--- /dev/null
+++ b/authors.txt
@@ -0,0 +1,2 @@
+kutscher = Konstantin Kutscher <kutscher@irmb.tu-bs.de>
+kucher = Konstantin Kutscher <kutscher@irmb.tu-bs.de>
\ No newline at end of file
diff --git a/depot/MeanValuesCoProcessor.cpp b/depot/MeanValuesCoProcessor.cpp
deleted file mode 100644
index 73bc3da3a..000000000
--- a/depot/MeanValuesCoProcessor.cpp
+++ /dev/null
@@ -1,67 +0,0 @@
-#include "MeanValuesCoProcessor.h"
-
-MeanValuesCoProcessor::MeanValuesCoProcessor(Grid3DPtr grid, UbSchedulerPtr s, const std::string &path,
- CommunicatorPtr comm, D3Q27IntegrateValuesHelperPtr ih, double factorPressure, double factorVelocity) :
- CoProcessor(grid, s),
- filename(path),
- comm(comm),
- ih(ih),
- factorPressure(factorPressure),
- factorVelocity(factorVelocity)
-{
- if (comm->getProcessID() == comm->getRoot())
- {
- FILE *output;
- output = fopen(filename.c_str(), "w");
- if (output == NULL)
- {
- std::string pathf = UbSystem::getPathFromString(filename);
- if (filename.size()>0) { UbSystem::makeDirectory(pathf); output = fopen(filename.c_str(), "w"); }
- if (output == NULL) throw UbException(UB_EXARGS, "can not open " + filename);
- }
-
- //fprintf(output, "%-18s %-18s %-18s %-18s %-18s %-18s %-18s\n", "time step", "area", "pressure [Pa*m^2]", "velocity [m^3/s]", "velocityX [m^3/s]", "velocityY [m^3/s]", "velocityZ [m^3/s]");
- fprintf(output, "%-18s;%-18s;%-18s;%-18s;%-18s\n", "time step", "area", "velocityX [m^3/s]", "velocityY [m^3/s]", "velocityZ [m^3/s]");
- fclose(output);
- }
-}
-//////////////////////////////////////////////////////////////////////////
-MeanValuesCoProcessor::~MeanValuesCoProcessor()
-{
-
-}
-//////////////////////////////////////////////////////////////////////////
-void MeanValuesCoProcessor::process(double step)
-{
- if (scheduler->isDue(step))
- collectPostprocessData(step);
-}
-//////////////////////////////////////////////////////////////////////////
-void MeanValuesCoProcessor::collectPostprocessData(double step)
-{
- ih->calculateMQ();
-
- if (comm->getProcessID() == comm->getRoot())
- {
- int istep = static_cast<int>(step);
- std::ofstream ostr;
- double nn = ih->getNumberOfFluidsNodes();
- //double p = (ih->getPress())*factorPressure;
- double vx1 = (ih->getVx1())*factorVelocity;
- double vx2 = (ih->getVx2())*factorVelocity;
- double vx3 = (ih->getVx3())*factorVelocity;
- //double vm = (ih->getVm())*factorVelocity;
-
- FILE *output;
- output = fopen(filename.c_str(), "a");
- if (output == NULL)
- throw UbException(UB_EXARGS, "can not open " + filename);
-
- //fprintf(output, "%-18d %-18d %-18.3f %-18.3f %-18.3f %-18.3f %-18.3f\n", istep, (int)nn, p, vm, vx1, vx2, vx3);
- fprintf(output, "%-18d %-18d %-18.3f %-18.3f %-18.3f\n", istep, (int)nn, vx1, vx2, vx3);
-
- fclose(output);
-
- UBLOG(logINFO, "D3Q27MeanValuesPostprocessor step: " << istep);
- }
-}
diff --git a/depot/MeanValuesCoProcessor.h b/depot/MeanValuesCoProcessor.h
deleted file mode 100644
index eb3b38381..000000000
--- a/depot/MeanValuesCoProcessor.h
+++ /dev/null
@@ -1,26 +0,0 @@
-#ifndef D3Q27MEANVALUESPOSTPROCESSOR_H
-#define D3Q27MEANVALUESPOSTPROCESSOR_H
-
-#include <CoProcessor.h>
-#include <D3Q27IntegrateValuesHelper.h>
-#include <LBMUnitConverter.h>
-
-class MeanValuesCoProcessor : public CoProcessor
-{
-public:
- MeanValuesCoProcessor(Grid3DPtr grid, UbSchedulerPtr s, const std::string &path,
- CommunicatorPtr comm, D3Q27IntegrateValuesHelperPtr ih, double factorPressure, double factorVelocity);
-
- virtual ~MeanValuesCoProcessor();
- void process(double step);
-
-private:
- D3Q27IntegrateValuesHelperPtr ih;
- CommunicatorPtr comm;
- std::string filename;
- double factorPressure, factorVelocity;
- void collectPostprocessData(double step);
-};
-
-#endif
-
diff --git a/depot/Particles.cpp b/depot/Particles.cpp
deleted file mode 100644
index 5e6e5141f..000000000
--- a/depot/Particles.cpp
+++ /dev/null
@@ -1,11 +0,0 @@
-#include "Particles.h"
-
-
-Particles::Particles()
-{
-}
-//////////////////////////////////////////////////////////////////////////
-Particles::~Particles()
-{
-}
-//////////////////////////////////////////////////////////////////////////
diff --git a/depot/Particles.h b/depot/Particles.h
deleted file mode 100644
index c133f9de4..000000000
--- a/depot/Particles.h
+++ /dev/null
@@ -1,29 +0,0 @@
-#ifndef PARTICLES_H
-#define PARTICLES_H
-
-#include <boost/shared_ptr.hpp>
-# define numOfvarClass 15
-
-class Particles;
-typedef boost::shared_ptr<Particles> ParticlesPtr;
-
-class Particles
-{
-public:
- double x,y,z,xold,yold,zold;
- double vxold,vyold,vzold;
- double vxoldf,vyoldf,vzoldf;
- int ID;
- int rankOfParticle;
- int level;
- //double dx;
-
-
- Particles();
- ~Particles();
-protected:
-private:
-
-};
-#endif
-
diff --git a/depot/PathLineCoProcessor.cpp b/depot/PathLineCoProcessor.cpp
deleted file mode 100644
index 7ca2741c0..000000000
--- a/depot/PathLineCoProcessor.cpp
+++ /dev/null
@@ -1,1505 +0,0 @@
-#include "PathLineCoProcessor.h"
-#include "LBMKernelETD3Q27.h"
-#include "SimulationParameters.h"
-#include "D3Q27ETBCProcessor.h"
-#include <vector>
-#include <string>
-#include <boost/foreach.hpp>
-#include "WbWriter.h"
-
-using namespace std;
-
-
-PathLineCoProcessor::PathLineCoProcessor(Grid3DPtr grid, const std::string& path,
- WbWriter* const writer, LBMUnitConverterPtr conv,
- UbSchedulerPtr s, CommunicatorPtr comm,
- double x1, double x2, double x3, LBMReal nue, D3Q27InterpolationProcessorPtr iProcessor)
- : CoProcessor(grid, s),
- path(path),
- comm(comm),
- writer(writer),
- conv(conv),
- x1(x1),
- x2(x2),
- x3(x3),
- nue(nue),
- iProcessor(iProcessor),
- vx1old(0.0),
- vx2old(0.0),
- vx3old(0.0),
- vx1oldf(0.1),
- vx2oldf(0.1),
- vx3oldf(0.1),
- isExtrapolation(0),
- istep(0),
- maxtag(-1),
- particleHasMass(true),
- rank(comm->getProcessID())
-{
- iHelper = D3Q27InterpolationHelperPtr(new D3Q27InterpolationHelper(iProcessor));
-}
-//////////////////////////////////////////////////////////////////////////
-PathLineCoProcessor::~PathLineCoProcessor()
-{
-
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::process(double step)
-{
- this->stepcheck=(int)step;
- double aa=scheduler->getMinEnd();
- if (step >= scheduler->getMaxBegin() && step <= scheduler->getMinEnd())
- {
- collectPostprocessData();
- }
-
- UBLOG(logDEBUG3, "PathLinePostprocessor::update:" << step);
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::collectPostprocessData()
-{
- //for (int n=0;n<Positions[][0].size();n++)
- //{
- //}
-
- int minInitLevel = this->grid->getCoarsestInitializedLevel();
- int maxInitLevel = this->grid->getFinestInitializedLevel();
- bool goinside=true;
-
- LBMKernelETD3Q27Ptr kernel;
- DistributionArray3DPtr distributions;
- BCArray3D<D3Q27BoundaryCondition> bcArray;
- Block3DPtr block;
- Block3DPtr block2;
- int isPointSuitable=0;//0 means false
- ///////send tau for printing to proccess 0
- double tau[6];//={tauxx,tauyy,tauzz,tauxy,tauxz,tauyz};
- for(int level = minInitLevel; level<=maxInitLevel; level++)
- {
- UbTupleInt3 blockIndexes = grid->getBlockIndexes(x1, x2, x3,level);
- block = grid->getBlock(val<1>(blockIndexes), val<2>(blockIndexes), val<3>(blockIndexes), level);
- if(!block) continue;
- if(block->isActive())
- {
- this->root=block->getRank();
- if(rank == block->getRank())
- {
- if(!checkNodes(block))
- {
- isPointSuitable=0;//0 meas false
- if(comm->getNumberOfProcesses() > 1) { sendtoOtherIsOk(isPointSuitable);}
- continue;
- }
- initialMovement(block, kernel, distributions, bcArray);
- isPointSuitable=1;//0 meas false
- if(comm->getNumberOfProcesses() > 1) { sendtoOtherIsOk(isPointSuitable);}
- break;
- }
- else
- {
- if(comm->getNumberOfProcesses() > 1)
- {
- MPI_Status status;
- MPI_Recv(&isPointSuitable,1, MPI_INT,this->root,this->rank,MPI_COMM_WORLD,&status);
- if (isPointSuitable)
- {
- break;
- }
- }
- }
- }
- }
- if (isPointSuitable)
- {
- if(comm->getNumberOfProcesses() > 1) updateDistributedValues();
- bool isPointSuitable2;
- checkLevel(block, kernel, distributions, bcArray,isPointSuitable2);
- if (isPointSuitable2)
- {
- if(block->isActive())
- {
- if (rank == block->getRank())
- {
- interpolMovement(block, kernel, distributions, bcArray,isExtrapolation,tau);
- if (isExtrapolation)
- {
- extrapolMovement(block, kernel, distributions, bcArray,tau);
- }
- }
- else
- {
- if(comm->getNumberOfProcesses() > 1)
- {
- //recive isExtrapolation
- if (this->rank!=this->root)
- {
- MPI_Status status;
- MPI_Recv(&isExtrapolation,1, MPI_INT,this->root,this->rank,MPI_COMM_WORLD,&status);
- }
- }
- }
- }
- if(comm->getNumberOfProcesses() > 1)
- {
- if (isExtrapolation)
- {
- if (this->rank!=this->root)
- {
- MPI_Status status;
- MPI_Recv(&maxtag,1, MPI_INT,this->root,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
- }
- }
- if (isExtrapolation)
- {
- if (this->rank!=this->root)
- {
- MPI_Status status;
- double input [7][25];
- double AllCell[7][8*27];
- int sizeinput=25;
- for (int i=0;i<maxtag;i++)
- {
- double vectforTrans[25];
- double *Cell=(double*)malloc(8*27*sizeof(double));
- MPI_Recv(vectforTrans,25, MPI_DOUBLE_PRECISION,this->root,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
- getAllInfoForCompare(vectforTrans,Cell);
-
- for (int j=0;j<sizeinput;j++) { input[i][j]=vectforTrans[j]; }
- for (int k=0;k<8*27;k++) { AllCell[i][k]=Cell[k]; }
- free (Cell);
- }
- for (int i=0;i<maxtag;i++)
- {
- double vectforTrans[25];
- double *Cell=(double*)malloc(8*27*sizeof(double));
- for (int j=0;j<sizeinput;j++) { vectforTrans[j]=input[i][j]; }
- for (int k=0;k<8*27;k++) { Cell[k]=AllCell[i][k]; }
-
- MPI_Send(vectforTrans,25, MPI_DOUBLE_PRECISION,this->root,i,MPI_COMM_WORLD);
- MPI_Send(Cell,8*27, MPI_DOUBLE_PRECISION,this->root,i,MPI_COMM_WORLD);
- free (Cell);
- }
- }
- }
- updateDistributedValues();
- }
-
- if(comm->getNumberOfProcesses() > 1) MPI_Bcast(tau,6, MPI_DOUBLE, this->root, MPI_COMM_WORLD);
- if (rank==0)
- {
- MPI_Status status;
- printPoint(tau);
- }
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::printPoint(double tau[])
-{
- double tauxx=tau[0];
- double tauyy=tau[1];
- double tauzz=tau[2];
- double tauxy=tau[3];
- double tauxz=tau[4];
- double tauyz=tau[5];
- datanames.resize(0);
- datanames.push_back("TimeStep");
- datanames.push_back("Vx");
- datanames.push_back("Vy");
- datanames.push_back("Vz");
- datanames.push_back("tauxx");
- datanames.push_back("tauyy");
- datanames.push_back("tauzz");
- datanames.push_back("tauxy");
- datanames.push_back("tauxz");
- datanames.push_back("tauyz");
- /*datanames.push_back("xoff");
- datanames.push_back("yoff");
- datanames.push_back("zoff");
- */
- data.resize(datanames.size());
-
- int index = 0;
- data[index++].push_back(istep++);
- data[index++].push_back(vx1old);
- data[index++].push_back(vx2old);
- data[index++].push_back(vx3old);
- data[index++].push_back(tauxx);
- data[index++].push_back(tauyy);
- data[index++].push_back(tauzz);
- data[index++].push_back(tauxy);
- data[index++].push_back(tauxz);
- data[index++].push_back(tauyz);
- //data[index++].push_back(xoff);
- //data[index++].push_back(yoff);
- //data[index++].push_back(zoff);
-
- nodes.push_back( makeUbTuple( double(x1), double(x2), double(x3)) );
-
- //if(scheduler->isDue((double)istep) )
- if (istep%1 ==0)
- {
- string test = writer->writeNodesWithNodeDataDouble(path,nodes,datanames,data);
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::clearData()
-{
- nodes.clear();
- datanames.clear();
- data.clear();
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::initialMovement(Block3DPtr block, LBMKernelETD3Q27Ptr& kernel, DistributionArray3DPtr& distributions, BCArray3D<D3Q27BoundaryCondition>& bcArray)
-{
- //UBLOG(logINFO, "addPostprocessData1");
- //root = rank;
-
- double dx = grid->getDeltaX(block);
- kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- distributions = kernel->getDataSet()->getFdistributions();
- bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
-
- x1old=x1;
- x2old=x2;
- x3old=x3;
-
- x1=x1old+vx1old*conv->getFactorTimeLbToW(grid->getDeltaX(block));
- x2=x2old+vx2old*conv->getFactorTimeLbToW(grid->getDeltaX(block));
- x3=x3old+vx3old*conv->getFactorTimeLbToW(grid->getDeltaX(block));
-
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::checkLevel(Block3DPtr& block, LBMKernelETD3Q27Ptr& kernel, DistributionArray3DPtr& distributions, BCArray3D<D3Q27BoundaryCondition>& bcArray,bool &isPointSuitable2)
-{
- int minInitLevel = this->grid->getCoarsestInitializedLevel();
- int maxInitLevel = this->grid->getFinestInitializedLevel();
- bool goinside=true;
- int isPointSuitable=0;//0 means false
- for(int level = minInitLevel; level<=maxInitLevel; level++)
- {
- UbTupleInt3 blockIndexes = grid->getBlockIndexes(x1, x2, x3,level);
- block = grid->getBlock(val<1>(blockIndexes), val<2>(blockIndexes), val<3>(blockIndexes), level);
- if(!block) continue;
- if(block->isActive())
- {
- this->root=block->getRank();
- if(rank == block->getRank())
- {
- if(!checkNodes(block))
- {
- isPointSuitable=0;//0 meas false
- if(comm->getNumberOfProcesses() > 1) { sendtoOtherIsOk(isPointSuitable);}
- continue;
- }
- kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- distributions = kernel->getDataSet()->getFdistributions();
- bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
- isPointSuitable=1;//0 meas false
- if(comm->getNumberOfProcesses() > 1) { sendtoOtherIsOk(isPointSuitable);}
- break;
- }
- else
- {
- if(comm->getNumberOfProcesses() > 1)
- {
- MPI_Status status;
- MPI_Recv(&isPointSuitable,1, MPI_INT,this->root,this->rank,MPI_COMM_WORLD,&status);
- if (isPointSuitable)
- {
- break;
- }
- }
- }
- }
- }
- if (isPointSuitable){isPointSuitable2=true; }
- else isPointSuitable2=false;
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::interpolMovement(Block3DPtr block, LBMKernelETD3Q27Ptr& kernel, DistributionArray3DPtr& distributions, BCArray3D<D3Q27BoundaryCondition>& bcArray, int &isExtrapolation,double tau[])
-{
- LBMReal tauxx, tauyy, tauzz, tauxy, tauxz, tauyz;
- LBMReal vx1,vx2,vx3,rho;
- UbTupleInt3 nodeIndexes = grid->getNodeIndexes(block, x1, x2, x3);
- LBMReal dx = grid->getDeltaX(block);
- D3Q27ICell iCell;
- LBMReal f[D3Q27System::ENDF+1];
- int ix1 = val<1>(nodeIndexes);
- int ix2 = val<2>(nodeIndexes);
- int ix3 = val<3>(nodeIndexes);
- double xoff=0.0;
- double yoff=0.0;
- double zoff=0.0;
- //typedef void (*CalcMacrosFct)(const LBMReal* const& /*feq[27]*/,LBMReal& /*(d)rho*/, LBMReal& /*vx1*/, LBMReal& /*vx2*/, LBMReal& /*vx3*/);
- /*
- CalcMacrosFct calcMacros = NULL;
-
- if(block->getKernel()->getCompressible())
- {
- calcMacros = &D3Q27System::calcCompMacroscopicValues;
- }
- else
- {
- calcMacros = &D3Q27System::calcIncompMacroscopicValues;
- }*/
-
- double Ucor=0,Vcor=0,Wcor=0;
- double tauxxC=00,tauyyC=0,tauzzC=0,tauxyC=0,tauxzC=0,tauyzC=0;
- double X_Wcor,Y_Wcor,Z_Wcor;
- double X_Pcor,Y_Pcor,Z_Pcor;
- double xp000,yp000,zp000;
- double Xw;
- double Yw;
- double Zw;
-
- UbTupleDouble3 orgNodeRW = grid->getNodeCoordinates(block, ix1, ix2, ix3);
- double ix1p= ( val<1>(orgNodeRW));
- double ix2p= ( val<2>(orgNodeRW));
- double ix3p= ( val<3>(orgNodeRW));
- if(!iProcessor->iCellHasSolid(bcArray, ix1, ix2, ix3))
- {
- iProcessor->readICell(distributions, iCell, ix1, ix2, ix3);
- double x1LB, x2LB, x3LB;
- UbTupleDouble3 orgNodeRW = grid->getNodeCoordinates(block, ix1, ix2, ix3);
- double offsetX1RW = abs(x1 - val<1>(orgNodeRW));
- double offsetX2RW = abs(x2 - val<2>(orgNodeRW));
- double offsetX3RW = abs(x3 - val<3>(orgNodeRW));
- x1LB = offsetX1RW / dx;
- x2LB = offsetX2RW / dx;
- x3LB = offsetX3RW / dx;
- x1LB -= 0.5-xoff;
- x2LB -= 0.5-yoff;
- x3LB -= 0.5-zoff;
- LBMReal omega = LBMSystem::calcCollisionFactor(nue, block->getLevel());
- /*iProcessor->interpolate8to1(iCell, f, x1LB, x2LB, x3LB, omega);
- calcMacros(f,rho,vx1,vx2,vx3);*/
- //iProcessor->interpolate8to1WithVelocity(iCell, f, x1LB, x2LB, x3LB, omega,vx1,vx2,vx3);
- iHelper->interpolate8to1WithVelocityWithShearStress(iCell, x1LB, x2LB, x3LB, omega,vx1,vx2,vx3,tauxx,tauyy,tauzz,tauxy,tauxz,tauyz);
- vx1 = vx1 * conv->getFactorVelocityLbToW();
- vx2 = vx2 * conv->getFactorVelocityLbToW();
- vx3 = vx3 * conv->getFactorVelocityLbToW();
- tau[0]=tauxx/grid->getDeltaX(block);
- tau[1]=tauyy/grid->getDeltaX(block);
- tau[2]=tauzz/grid->getDeltaX(block);
- tau[3]=tauxy/grid->getDeltaX(block);
- tau[4]=tauxz/grid->getDeltaX(block);
- tau[5]=tauyz/grid->getDeltaX(block);
-
- //if (particleHasMass) { CalcVelParticle(dx,vx1,vx2,vx3,vx1old,vx2old,vx3old,vx1oldf,vx2oldf,vx3oldf);}
- if (particleHasMass) { CalcVelParticle2(dx,vx1,vx2,vx3,vx1old,vx2old,vx3old,vx1oldf,vx2oldf,vx3oldf);}
- //heuns method
- x1 = (x1old) + 1.0/2.0*(vx1 + vx1old)*conv->getFactorTimeLbToW(dx);
- x2 = (x2old) + 1.0/2.0*(vx2 + vx2old)*conv->getFactorTimeLbToW(dx);
- x3 = (x3old) + 1.0/2.0*(vx3 + vx3old)*conv->getFactorTimeLbToW(dx);
-
- vx1old = vx1;
- vx2old = vx2;
- vx3old = vx3;
-
- //send isExtrapolation false
- isExtrapolation=0;//1 means true
- if(comm->getNumberOfProcesses() > 1)
- {
- int size=comm->getNumberOfProcesses();
- for (int i=0;i<size;i++)
- {
- if (i!=rank)
- {
- MPI_Send(&isExtrapolation,1, MPI_INT,i,i,MPI_COMM_WORLD);
- }
- }
- }
- }
- else
- {
- //send isExtrapolation true
- isExtrapolation=1;//1 means true
- if(comm->getNumberOfProcesses() > 1)
- {
- int size=comm->getNumberOfProcesses();
- for (int i=0;i<size;i++)
- {
- if (i!=rank)
- {
- MPI_Send(&isExtrapolation,1, MPI_INT,i,i,MPI_COMM_WORLD);
- }
- }
- }
-
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::extrapolMovement(Block3DPtr block, LBMKernelETD3Q27Ptr& kernel, DistributionArray3DPtr& distributions, BCArray3D<D3Q27BoundaryCondition>& bcArray,double tau[])
-{
- LBMReal tauxx, tauyy, tauzz, tauxy, tauxz, tauyz;
- LBMReal vx1,vx2,vx3,rho;
- UbTupleInt3 nodeIndexes = grid->getNodeIndexes(block, x1, x2, x3);
- LBMReal dx = grid->getDeltaX(block);
- int level = block->getLevel();
- double deltaT= 1.0/(double)(1<<level);
- LBMReal omega = LBMSystem::calcCollisionFactor(nue,level);
- D3Q27ICell iCell;
- LBMReal f[D3Q27System::ENDF+1];
- int ix1 = val<1>(nodeIndexes);
- int ix2 = val<2>(nodeIndexes);
- int ix3 = val<3>(nodeIndexes);
- double xoff=0.0;
- double yoff=0.0;
- double zoff=0.0;
- double Ucor=0,Vcor=0,Wcor=0;
- double tauxxC=00,tauyyC=0,tauzzC=0,tauxyC=0,tauxzC=0,tauyzC=0;
- double X_Wcor,Y_Wcor,Z_Wcor;
- double X_Pcor,Y_Pcor,Z_Pcor;
- double xp000,yp000,zp000;
- double Xw=999.0;
- double Yw=999.0;
- double Zw=999.0;
-
- UbTupleDouble3 orgNodeRW = grid->getNodeCoordinates(block, ix1, ix2, ix3);
- double ix1p= ( val<1>(orgNodeRW));
- double ix2p= ( val<2>(orgNodeRW));
- double ix3p= ( val<3>(orgNodeRW));
- {
- UbTupleDouble3 orgNodeRW = grid->getNodeCoordinates(block, ix1, ix2, ix3);
- double ix1ph = ( val<1>(orgNodeRW));
- double ix2ph = ( val<2>(orgNodeRW));
- double ix3ph = ( val<3>(orgNodeRW));
- double cofWeightx=999.0,cofWeighty=999.0,cofWeightz=999.0;
- double A,B,C,D,ii=0.0;
- findPlane(ix1,ix2,ix3,grid,block,A,B,C,D,ii);
- double s = A*x1 + B*x2 + C*x3 + D;//The sign of s = Ax + By + Cz + D determines which side the point (x,y,z) lies with respect to the plane. If s > 0 then the point lies on the same side as the normal (A,B,C). If s < 0 then it lies on the opposite side, if s = 0 then the point (x,y,z) lies on the plane.
- if (s>0){s=1;} else if (s<0){s=-1;}else {s=0;}
- double normalDis=((A*x1 + B*x2 + C*x3 + D)/sqrt(A*A+B*B+C*C));///distance point to plane xp-Xw=distance
- double di=A/sqrt(A*A+B*B+C*C); double dj=B/sqrt(A*A+B*B+C*C); double dk=C/sqrt(A*A+B*B+C*C);
- double XXw=x1-di*normalDis; double YYw=x2-dj*normalDis; double ZZw=x3-dk*normalDis;
- findInitialCell(s,di,dj,dk,dx,ix1ph,ix2ph,ix3ph,xp000,yp000,zp000);//find initial cell with xp000,yp000,zp000
- /////////////////////////////////////////////////////////////////////////
- int level=block->getLevel();
- int Rankx,Ranky,Rankz,Rankxy,Rankxz,Rankyz,Rankxyz;
- double disx=9990,disy=9990,disz=9990,disxy=9990,disxz=9990,disyz=9990,disxyz=9990;
- getAllRankWithAllPositions(ix1ph,ix2ph,ix3ph,xp000,yp000,zp000,Rankx,Ranky,Rankz,Rankxy,Rankxz,Rankyz,Rankxyz,level);
- //getAllRankWithAllPositions2(ix1ph,ix2ph,ix3ph,xp000,yp000,zp000,Rankx,Ranky,Rankz,Rankxy,Rankxz,Rankyz,Rankxyz);
-
- int allRank[7]={ Rankx,Ranky,Rankz,Rankxy,Rankxz,Rankyz,Rankxyz};
- int allTag[7];
- getAllTag(allRank,allTag);
- int numberproccessused=0;
- int numberpro=comm->getNumberOfProcesses();
- int *mainproccess=new int[numberpro] ;
- sendMaxtagToAllProccess(allRank,allTag,numberpro,numberproccessused,mainproccess);
-
- //{0=xp000, 1=ix2ph, 2=ix3ph,3=x1,4=x2,5=x3,6=Xw,7=Yw,8=Zw,9=A,10=B,11=C,12=D,
- //13=normalDis,14=di,15=dj,16=dk,17=cofWeightx,18=cofWeighty,19=cofWeightz,,20=dx,21=omega,22=level,23=disx,24=dir };
- const int sizeinput=25;
- double input [7][sizeinput]=
- {{xp000, ix2ph, ix3ph,x1,x2,x3,Xw,Yw,Zw,A,B,C,D,normalDis,di,dj,dk,cofWeightx,cofWeighty,cofWeightz,dx,omega,level,disx ,dirx }
- ,{ix1ph, yp000, ix3ph,x1,x2,x3,Xw,Yw,Zw,A,B,C,D,normalDis,di,dj,dk,cofWeightx,cofWeighty,cofWeightz,dx,omega,level,disy ,diry }
- ,{ix1ph, ix2ph, zp000,x1,x2,x3,Xw,Yw,Zw,A,B,C,D,normalDis,di,dj,dk,cofWeightx,cofWeighty,cofWeightz,dx,omega,level,disz ,dirz }
- ,{xp000, yp000, ix3ph,x1,x2,x3,Xw,Yw,Zw,A,B,C,D,normalDis,di,dj,dk,cofWeightx,cofWeighty,cofWeightz,dx,omega,level,disxy ,dirxy }
- ,{xp000, ix2ph, zp000,x1,x2,x3,Xw,Yw,Zw,A,B,C,D,normalDis,di,dj,dk,cofWeightx,cofWeighty,cofWeightz,dx,omega,level,disxz ,dirxz }
- ,{ix1ph, yp000, zp000,x1,x2,x3,Xw,Yw,Zw,A,B,C,D,normalDis,di,dj,dk,cofWeightx,cofWeighty,cofWeightz,dx,omega,level,disyz ,diryz }
- ,{xp000, yp000, zp000,x1,x2,x3,Xw,Yw,Zw,A,B,C,D,normalDis,di,dj,dk,cofWeightx,cofWeighty,cofWeightz,dx,omega,level,disxyz,dirxyz}};
-
- double AllCell[7][8*27];
- int index;
- SendAndReceiveData(input,AllCell,allRank,allTag);
-
- double dis[7]={input[0][23],input[1][23],input[2][23],input[3][23],input[4][23],input[5][23],input[6][23]};
- rearangedDouble(dis);
- UbTupleInt3 blockIndexes;
-
- getAfterCompare(dis,input,AllCell,xp000, yp000,zp000,Xw,Yw,Zw,cofWeightx,cofWeighty,cofWeightz,dx,omega,iCell);
-
- double x1LB, x2LB, x3LB;
- double offsetX1RW000 = (x1 - xp000);
- double offsetX2RW000 = (x2 - yp000);
- double offsetX3RW000 = (x3 - zp000);
-
- x1LB = offsetX1RW000 / dx;
- x2LB = offsetX2RW000 / dx;
- x3LB = offsetX3RW000 / dx;
- x1LB -= 0.5;
- x2LB -= 0.5;
- x3LB -= 0.5;
- // outICell(iCell);
- //iProcessor->interpolate8to1WithVelocity(iCell, f, x1LB, x2LB, x3LB, omega,vx1,vx2,vx3);
- iHelper->interpolate8to1WithVelocityWithShearStress(iCell, x1LB, x2LB, x3LB, omega,vx1,vx2,vx3,tauxx,tauyy,tauzz,tauxy,tauxz,tauyz);
- addCorrection(vx1,vx2,vx3,tauxx,tauyy,tauzz,tauxy,tauxz,tauyz,dx,iCell,Xw, Yw, Zw, omega,cofWeightx,cofWeighty,cofWeightz,ii, x1LB, x2LB, x3LB,di,dj,dk);
- tau[0]=tauxx;
- tau[1]=tauyy;
- tau[2]=tauzz;
- tau[3]=tauxy;
- tau[4]=tauxz;
- tau[5]=tauyz;
- // if (particleHasMass) { CalcVelParticle(dx,vx1,vx2,vx3,vx1old,vx2old,vx3old,vx1oldf,vx2oldf,vx3oldf);}
- if (particleHasMass) { CalcVelParticle2(dx,vx1,vx2,vx3,vx1old,vx2old,vx3old,vx1oldf,vx2oldf,vx3oldf);}
- //heuns method
- x1 = (x1old) + 1.0/2.0*(vx1 + vx1old)*conv->getFactorTimeLbToW(dx);
- x2 = (x2old) + 1.0/2.0*(vx2 + vx2old)*conv->getFactorTimeLbToW(dx);
- x3 = (x3old) + 1.0/2.0*(vx3 + vx3old)*conv->getFactorTimeLbToW(dx);
- collWall(A,B,C,D,x1,x2,x3,x1old,x2old,x3old,dx,vx1,vx2,vx3,ii);
- vx1old = vx1;
- vx2old = vx2;
- vx3old = vx3;
- }
-}
-//////////////////////////////////////////////////////////////////////////
-bool PathLineCoProcessor::checkNodes( Block3DPtr block)
-{
- bool result = true;
- LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- BCArray3D<D3Q27BoundaryCondition> bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
- D3Q27BoundaryConditionPtr bcPtr;
-
- double x1_ch = x1;// + vx1old*conv->getfactorTimeShouldMultiplebyDx()*grid->getDeltaX(block);
- double x2_ch = x2;// + vx2old*conv->getfactorTimeShouldMultiplebyDx()*grid->getDeltaX(block);
- double x3_ch = x3;// + vx3old*conv->getfactorTimeShouldMultiplebyDx()*grid->getDeltaX(block);
-
- int nx1 = static_cast<int>(bcArray.getNX1());
- int nx2 = static_cast<int>(bcArray.getNX2());
- int nx3 = static_cast<int>(bcArray.getNX3());
-
- UbTupleInt3 nodeIndexes = grid->getNodeIndexes(block, x1_ch, x2_ch, x3_ch);
-
- int ix1 = val<1>(nodeIndexes);
- int ix2 = val<2>(nodeIndexes);
- int ix3 = val<3>(nodeIndexes);
- bool outlet=false;
-
- for (int xx3 = ix3; xx3 <= ix3 + 1 && ix3 + 1 < nx3; xx3++)
- for(int xx2 = ix2; xx2 <= ix2 + 1 && ix2 + 1 < nx2; xx2++)
- for(int xx1 = ix1; xx1 <= ix1 + 1 && ix1 + 1 < nx1; xx1++)
- {
- bcPtr=bcArray.getBC(xx1,xx2,xx3);
- for(int fdir=D3Q27System::STARTF; fdir<=D3Q27System::ENDF; fdir++)
- {
- if ( bcPtr != NULL)
- {
- if (bcPtr->hasDensityBoundaryFlag(fdir))
- {
- if(outlet)break;
- outlet=true;
- break;
- }
- }
- }
- result = result && !bcArray.isUndefined(xx1, xx2, xx3) && !outlet;
- }
-
- return result;
-
-}
-//////////////////////////////////////////////////////////////////////////
-bool PathLineCoProcessor::checkNodes2( Block3DPtr block,double x11,double x22,double x33)
-{
- bool result = true;
- LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- BCArray3D<D3Q27BoundaryCondition> bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
- D3Q27BoundaryConditionPtr bcPtr;
-
- double x1_ch = x11;// + vx1old*conv->getfactorTimeShouldMultiplebyDx()*grid->getDeltaX(block);
- double x2_ch = x22;// + vx2old*conv->getfactorTimeShouldMultiplebyDx()*grid->getDeltaX(block);
- double x3_ch = x33;// + vx3old*conv->getfactorTimeShouldMultiplebyDx()*grid->getDeltaX(block);
-
- int nx1 = static_cast<int>(bcArray.getNX1());
- int nx2 = static_cast<int>(bcArray.getNX2());
- int nx3 = static_cast<int>(bcArray.getNX3());
-
- UbTupleInt3 nodeIndexes = grid->getNodeIndexes(block, x1_ch, x2_ch, x3_ch);
-
- int ix1 = val<1>(nodeIndexes);
- int ix2 = val<2>(nodeIndexes);
- int ix3 = val<3>(nodeIndexes);
- bool outlet=false;
-
- for (int xx3 = ix3; xx3 <= ix3 + 1 && ix3 + 1 < nx3; xx3++)
- for(int xx2 = ix2; xx2 <= ix2 + 1 && ix2 + 1 < nx2; xx2++)
- for(int xx1 = ix1; xx1 <= ix1 + 1 && ix1 + 1 < nx1; xx1++)
- {
- bcPtr=bcArray.getBC(xx1,xx2,xx3);
- for(int fdir=D3Q27System::STARTF; fdir<=D3Q27System::ENDF; fdir++)
- {
- if ( bcPtr != NULL)
- {
- if (bcPtr->hasDensityBoundaryFlag(fdir))
- {
- if(outlet)break;
- outlet=true;
- break;
- }
- }
- }
- result = result && !bcArray.isUndefined(xx1, xx2, xx3) && !outlet;
- }
-
- return result;
-
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::outICell(D3Q27ICell& iCell)
-{
- std::ofstream ostr;
- string fname = path+"_iCell";
- ostr.open(fname.c_str(), std::ios_base::out);
- if(!ostr)
- {
- ostr.clear();
- string path = UbSystem::getPathFromString(fname);
- if(path.size()>0){ UbSystem::makeDirectory(path); ostr.open(fname.c_str(), std::ios_base::out);}
- if(!ostr) throw UbException(UB_EXARGS,"couldn't open file "+fname);
- }
- ostr << "step = "<< stepcheck << endl;
- for (int i = 0; i < 27; i++)
- {
- ostr << "iCell.BNE[" << i <<"] = "<< iCell.BNE[i] << "\t"<<"iCell.BNW[" << i <<"] = "<< iCell.BNW[i] << "\t"<<"iCell.BSE[" << i <<"] = "<< iCell.BSE[i]<<"\t"<<"iCell.BSW[" << i <<"] = "<< iCell.BSW[i]<<"\t"<<"iCell.TNE[" << i <<"] = "<< iCell.TNE[i] << "\t"<<"iCell.TNW[" << i <<"] = "<< iCell.TNW[i] << "\t"<<"iCell.TSE[" << i <<"] = "<< iCell.TSE[i]<<"\t"<<"iCell.TSW[" << i <<"] = "<< iCell.TSW[i]<<"\t"<< endl;
-
- }
- ostr.close();
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::updateDistributedValues()
-{
- int size;
- MPI_Status status;
- MPI_Comm_size(MPI_COMM_WORLD,&size);
- const int row=12;
- double values[row];//
- if (this->root == this->rank)
- {
- values[0]=x1;
- values[1]=x2;
- values[2]=x3;
- values[3]=x1old;
- values[4]=x2old;
- values[5]=x3old;
- values[6]=vx1old;
- values[7]=vx2old;
- values[8]=vx3old;
- values[9]=vx1oldf;
- values[10]=vx2oldf;
- values[11]=vx3oldf;
- for (int i=0;i<size;i++)
- {
- if (i!=this->root)
- {
- MPI_Send(&values,row, MPI_DOUBLE_PRECISION,i,i,MPI_COMM_WORLD);
- }
- }
- }
- else{
- MPI_Recv(&values,row, MPI_DOUBLE_PRECISION,this->root,this->rank,MPI_COMM_WORLD,&status);
- x1 = values[0];
- x2 = values[1];
- x3 = values[2];
- x1old = values[3];
- x2old = values[4];
- x3old = values[5];
- vx1old = values[6];
- vx2old = values[7];
- vx3old = values[8];
- vx1oldf =values[9];
- vx2oldf=values[10];
- vx3oldf=values[11];
- }
-}
-//////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::findPlane(int ix1,int ix2,int ix3,Grid3DPtr grid,Block3DPtr block,double &A,double &B,double &C,double &D,double &ii)
-{
- double x1plane=0.0,y1plane=0.0,z1plane=0.0;
- double x2plane=0.0,y2plane=0.0,z2plane=0.0;
- double x3plane=0.0,y3plane=0.0,z3plane=0.0;
- D3Q27BoundaryConditionPtr bcPtr;
- double dx = grid->getDeltaX(block);
- LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- DistributionArray3DPtr distributions = kernel->getDataSet()->getFdistributions();
- BCArray3D<D3Q27BoundaryCondition> bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
-
- for (int k = ix3; k <= ix3 + 1; k++){
- for(int j = ix2; j <= ix2 + 1; j++){
- for(int i = ix1; i <= ix1 + 1; i++)
- {
- UbTupleDouble3 pointplane1 = grid->getNodeCoordinates(block, i, j, k);
-
- double iph=val<1>(pointplane1);
- double jph=val<2>(pointplane1);
- double kph=val<3>(pointplane1);
-
- if(!bcArray.isSolid(i, j, k))
- {
- bcPtr=bcArray.getBC(i,j,k);
- if(bcPtr)
- {
- for(int fdir=D3Q27System::FSTARTDIR; fdir<=D3Q27System::FENDDIR; fdir++)
- {
- if( ii<=2)
- {
- LBMReal q = bcPtr->getQ(fdir);
- if (q!=999.00000)
- {
-
- if ( fdir==D3Q27System::E )
- {
- if (i+q<=ix1+1)
- {
- if (ii==0) { x1plane=iph+q*dx; y1plane=jph; z1plane=kph; ii++; }
- else if (ii==1) { x2plane=iph+q*dx; y2plane=jph; z2plane=kph; if (x1plane!=x2plane||y1plane!=y2plane||z1plane!=z2plane) ii++; }
- else if(ii==2) { x3plane=iph+q*dx; y3plane=jph; z3plane=kph; if ((x3plane!=x1plane||y3plane!=y1plane||z3plane!=z1plane)&&(x2plane!=x3plane||y2plane!=y3plane||z2plane!=z3plane)) ii++;}
- }
- }
- if ( fdir==D3Q27System::W )
- {
- if (i-q>=ix1)
- {
- if (ii==0) { x1plane=iph-q*dx; y1plane=jph; z1plane=kph; ii++; }
- else if (ii==1) { x2plane=iph-q*dx; y2plane=jph; z2plane=kph; if (x1plane!=x2plane||y1plane!=y2plane||z1plane!=z2plane) ii++; }
- else if(ii==2) { x3plane=iph-q*dx; y3plane=jph; z3plane=kph; if ((x3plane!=x1plane||y3plane!=y1plane||z3plane!=z1plane)&&(x2plane!=x3plane||y2plane!=y3plane||z2plane!=z3plane)) ii++;}
- }
- }
- if ( fdir==D3Q27System::N )
- {
- if(j+q<=ix2+1)
- {
- if (ii==0) { x1plane=iph; y1plane=jph+q*dx; z1plane=kph; ii++; }
- else if (ii==1) { x2plane=iph; y2plane=jph+q*dx; z2plane=kph; if (x1plane!=x2plane||y1plane!=y2plane||z1plane!=z2plane) ii++; }
- else if (ii==2) { x3plane=iph; y3plane=jph+q*dx; z3plane=kph; if ((x3plane!=x1plane||y3plane!=y1plane||z3plane!=z1plane)&&(x2plane!=x3plane||y2plane!=y3plane||z2plane!=z3plane)) ii++;}
- }
- }
- if ( fdir==D3Q27System::S )
- {
- if (j-q>=ix2)
- {
- if (ii==0) { x1plane=iph; y1plane=jph-q*dx; z1plane=kph; ii++; }
- else if (ii==1) { x2plane=iph; y2plane=jph-q*dx; z2plane=kph; if (x1plane!=x2plane||y1plane!=y2plane||z1plane!=z2plane) ii++; }
- else if (ii==2) { x3plane=iph; y3plane=jph-q*dx; z3plane=kph; if ((x3plane!=x1plane||y3plane!=y1plane||z3plane!=z1plane)&&(x2plane!=x3plane||y2plane!=y3plane||z2plane!=z3plane)) ii++;}
- }
- }
-
- if ( fdir==D3Q27System::T )
- {
- if(k+q<=ix3+1)
- {
- if (ii==0) { x1plane=iph; y1plane=jph; z1plane=kph+q*dx; ii++; }
- else if (ii==1) { x2plane=iph; y2plane=jph; z2plane=kph+q*dx; if (x1plane!=x2plane||y1plane!=y2plane||z1plane!=z2plane) ii++; }
- else if (ii==2) { x3plane=iph; y3plane=jph; z3plane=kph+q*dx; if ((x3plane!=x1plane||y3plane!=y1plane||z3plane!=z1plane)&&(x2plane!=x3plane||y2plane!=y3plane||z2plane!=z3plane)) ii++;}
- }
- }
- if ( fdir==D3Q27System::B )
- {
- if (k-q>=ix3)
- {
- if (ii==0) { x1plane=iph; y1plane=jph; z1plane=kph-q*dx; ii++; }
- else if (ii==1) { x2plane=iph; y2plane=jph; z2plane=kph-q*dx; if (x1plane!=x2plane||y1plane!=y2plane||z1plane!=z2plane) ii++; }
- else if (ii==2) { x3plane=iph; y3plane=jph; z3plane=kph-q*dx; if ((x3plane!=x1plane||y3plane!=y1plane||z3plane!=z1plane)&&(x2plane!=x3plane||y2plane!=y3plane||z2plane!=z3plane)) ii++;}
- }
- }
-
- }
- }
- }
-
- }
-
- }
- }
- }
- }
-
- A = y1plane* (z2plane - z3plane) + y2plane*(z3plane - z1plane) + y3plane* (z1plane - z2plane);
- B = z1plane* (x2plane - x3plane) + z2plane*(x3plane - x1plane) + z3plane* (x1plane - x2plane) ;
- C = x1plane* (y2plane - y3plane) + x2plane*(y3plane - y1plane) + x3plane* (y1plane - y2plane) ;
- D =-( x1plane*(y2plane*z3plane - y3plane*z2plane)+x2plane*(y3plane*z1plane - y1plane*z3plane) + x3plane* (y1plane* z2plane - y2plane* z1plane));
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::finddisPointToCentercell(double x1,double x2,double x3,double xp000, double yp000, double zp000,double &dis,double &dx,
- double &deltaT,Block3DPtr &block,int level)
-{
- UbTupleInt3 blockIndexes = grid->getBlockIndexes(xp000, yp000, zp000,level);
- block = grid->getBlock(val<1>(blockIndexes), val<2>(blockIndexes), val<3>(blockIndexes), level);
- LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- if (!kernel) cout<<__FILE__<<" "<<__LINE__<<endl;
- else
- {
- DistributionArray3DPtr distributions = kernel->getDataSet()->getFdistributions();
- BCArray3D<D3Q27BoundaryCondition> bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
- if(block && block->isActive())
- {
- deltaT= 1.0/(double)(1<<level);
- dx = grid->getDeltaX(block);
- UbTupleInt3 nodeIndexes = grid->getNodeIndexes(block, xp000, yp000, zp000);
- int ix1Se = val<1>(nodeIndexes);
- int ix2Se = val<2>(nodeIndexes);
- int ix3Se = val<3>(nodeIndexes);
- if (!iProcessor->iCellHasSolid(bcArray, ix1Se, ix2Se, ix3Se ))
- {
- UbTupleDouble3 pointCell = grid->getNodeCoordinates(block, ix1Se, ix2Se, ix3Se); double iph=val<1>(pointCell)+dx/2.0; double jph=val<2>(pointCell)+dx/2.0; double kph=val<3>(pointCell)+dx/2.0;
- dis=sqrt((iph-x1)*(iph-x1)+(jph-x2)*(jph-x2)+(kph-x3)*(kph-x3));
- }
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-//void D3Q27PathLinePostprocessor::finddisPointToCentercell(double x1,double x2,double x3,double xp000, double yp000, double zp000,double &dis,double &dx,
-// double &deltaT,Block3DPtr &block)
-//{
-
-// int minInitLevel = this->grid->getCoarsestInitializedLevel();
-// int maxInitLevel = this->grid->getFinestInitializedLevel();
-// bool goinside=true;
-// Block3DPtr block2;
-// for(int level = minInitLevel; level<=maxInitLevel; level++)
-// {
-// if (goinside )
-// {
-// UbTupleInt3 blockIndexes = grid->getBlockIndexes(xp000, yp000, zp000,level);
-// block = grid->getBlock(val<1>(blockIndexes), val<2>(blockIndexes), val<3>(blockIndexes), level);
-// LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
-// if (kernel)
-// {
-// DistributionArray3DPtr distributions = kernel->getDataSet()->getFdistributions();
-// BCArray3D<D3Q27BoundaryCondition> bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
-// if(block && block->isActive())
-// {
-// if(!checkNodes(block))continue;
-// goinside=false;
-// deltaT= 1.0/(double)(1<<level);
-// dx = grid->getDeltaX(block);
-// UbTupleInt3 nodeIndexes = grid->getNodeIndexes(block, xp000, yp000, zp000);
-// int ix1Se = val<1>(nodeIndexes);
-// int ix2Se = val<2>(nodeIndexes);
-// int ix3Se = val<3>(nodeIndexes);
-// if (!iProcessor->iCellHasSolid(bcArray, ix1Se, ix2Se, ix3Se ))
-// {
-// UbTupleDouble3 pointCell = grid->getNodeCoordinates(block, ix1Se, ix2Se, ix3Se); double iph=val<1>(pointCell)+dx/2.0; double jph=val<2>(pointCell)+dx/2.0; double kph=val<3>(pointCell)+dx/2.0;
-// dis=sqrt((iph-x1)*(iph-x1)+(jph-x2)*(jph-x2)+(kph-x3)*(kph-x3));
-// block2=block;
-// }
-// }
-// }
-// }
-// }
-
-//}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::getWallPosition(double input[])
-{
-
- //{0=xp000, 1=ix2ph, 2=ix3ph,3=x1,4=x2,5=x3,6=Xw,7=Yw,8=Zw,9=A,10=B,11=C,12=D,13=normalDis,14=di,15=dj,16=dk,
- //17=cofWeightx,18=cofWeighty,19=cofWeightz,,20=dx,21=omega,22=nue,23=disx,24=dir };
- double dx=input[20];double xp000=input[0];double yp000=input[1];double zp000=input[2];double Xw=input[6];double Yw=input[7];double Zw=input[8];
- double x1=input[3] ;double x2=input[4] ;double x3=input[5]; double A=input[9];double B=input[10];double C=input[11];double D=input[12];
- double normalDis=input[13]; double di=input[14];double dj=input[15];double dk=input[16]; double cofWeightx=input[17];double cofWeighty=input[18];
- double cofWeightz=input[19];int dir=(int)input[24];
-
- //char dirx,diry,dirz,dirxy,dirxz,diryz,dirxyz;
- if (dir==dirx ) { Xw=(B*x2+C*x3+D)/(-A); Yw=x2; Zw=x3; cofWeightx=1; cofWeighty=0;cofWeightz=0;}
- else if (dir==diry ) { Xw=x1; Yw=(A*x1+C*x3+D)/(-B); Zw=x3; cofWeightx=0; cofWeighty=1;cofWeightz=0;}
- else if (dir==dirz ) { Xw=x1; Yw=x3; Zw=(B*x2+A*x1+D)/(-C); cofWeightx=0; cofWeighty=0;cofWeightz=1; }
- else if (dir==dirxy ) { Xw=x1-di*normalDis; Yw=x2-dj*normalDis; Zw=x3; cofWeightx=1; cofWeighty=1;cofWeightz=0;}
- else if (dir==dirxz ) { Xw=x1-di*normalDis; Yw=x2; Zw=x3-dk*normalDis; cofWeightx=1; cofWeighty=0;cofWeightz=1; }
- else if (dir==diryz ) { Xw=x1; Yw=x2-dj*normalDis; Zw=x3-dk*normalDis; cofWeightx=0; cofWeighty=1;cofWeightz=1; }
- else if (dir==dirxyz) { Xw=x1-di*normalDis; Yw=x2-dj*normalDis; Zw=x3-dk*normalDis; cofWeightx=1; cofWeighty=1;cofWeightz=1; }
-
- Xw=(Xw-xp000)/dx-0.5;
- Yw=(Yw-yp000)/dx-0.5;
- Zw=(Zw-zp000)/dx-0.5;
-
- input[6]=Xw;
- input[7]=Yw;
- input[8]=Zw;
- input[17]=cofWeightx;
- input[18]=cofWeighty;
- input[19]=cofWeightz;
-
-}
-////////////////////////////////////////////////////////////////////////////
-// void D3Q27PathLinePostprocessor::getWallPosition(double dixmax,double disx,double disy,double disz,double disxy,double disxz,double disyz,double disxyz
-// ,double &Xw,double &Yw,double &Zw,double x1,double x2,double x3,double A,double B,double C,double D,double normalDis,double di,double dj,double dk,
-// double &cofWeightx,double &cofWeighty,double &cofWeightz)
-// {
-// if (dixmax==disx) { Xw=(B*x2+C*x3+D)/(-A); Yw=x2; Zw=x3; cofWeightx=1; cofWeighty=0;cofWeightz=0;}
-// else if (dixmax==disy) { Xw=x1; Yw=(A*x1+C*x3+D)/(-B); Zw=x3; cofWeightx=0; cofWeighty=1;cofWeightz=0;}
-// else if (dixmax==disz) { Xw=x1; Yw=x3; Zw=(B*x2+A*x1+D)/(-C); cofWeightx=0; cofWeighty=0;cofWeightz=1;}
-// else if (dixmax==disxy) { Xw=x1-di*normalDis; Yw=x2-dj*normalDis; Zw=x3; cofWeightx=1; cofWeighty=1;cofWeightz=0;}
-// else if (dixmax==disxz) { Xw=x1-di*normalDis; Yw=x2; Zw=x3-dk*normalDis; cofWeightx=1; cofWeighty=0;cofWeightz=1;}
-// else if (dixmax==disyz) { Xw=x1; Yw=x2-dj*normalDis; Zw=x3-dk*normalDis; cofWeightx=0; cofWeighty=1;cofWeightz=1;}
-// else if (dixmax==disxyz){ Xw=x1-di*normalDis; Yw=x2-dj*normalDis; Zw=x3-dk*normalDis; cofWeightx=1; cofWeighty=1;cofWeightz=1;}
-// }
-////////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::findInitialCell(double s,double di,double dj,double dk,double dx,double ix1ph,double ix2ph,double ix3ph,double &xp000,double &yp000,double &zp000)
-{
- if (s*di>0) { xp000=ix1ph+dx; }
- else if(s*di<0) { xp000=ix1ph-dx; }
- else { xp000=ix1ph; }
- if (s*dj>0) { yp000=ix2ph+dx; }
- else if(s*dj<0) { yp000=ix2ph-dx; }
- else { yp000=ix2ph; }
- if (s*dk>0) { zp000=ix3ph+dx; }
- else if(s*dk<0) { zp000=ix3ph-dx; }
- else { zp000=ix3ph; }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::rearangedDouble( double dis[7])
-{
- double item;
- for (int i=6;i>0;i--){
- for(int j=0;j<i;j++){
- if (dis[j]>dis[j+1])
- {
- item=dis[j];
- dis[j]=dis[j+1];
- dis[j+1]=item;
- }
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::rearangeInt( int dis[7])
-{
- double item;
- for (int i=6;i>0;i--){
- for(int j=0;j<i;j++){
- if (dis[j]>dis[j+1])
- {
- item=dis[j];
- dis[j]=dis[j+1];
- dis[j+1]=(int)item;
- }
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::addCorrection(double &vx1,double &vx2,double &vx3,double &tauxx,double &tauyy,double &tauzz,double &tauxy,double &tauxz,double &tauyz,
- double dx,D3Q27ICell iCell,double Xw, double Yw, double Zw,double omega,double cofWeightx,double cofWeighty,double cofWeightz,double ii,double x1LB,double x2LB,double x3LB,double di,double dj,double dk)
-{
- LBMReal f[D3Q27System::ENDF+1];
- double Udif,Vdif,Wdif;
- iHelper->interpolate8to1WithVelocity(iCell, Xw, Yw, Zw, omega,Udif,Vdif,Wdif);
- double Ucx,Ucy,Ucz;
- double Vcx,Vcy,Vcz;
- double Wcx,Wcy,Wcz;
-
- double dxUcx,dyUcy,dzUcz;
- double dxVcx,dyVcy,dzVcz;
- double dxWcx,dyWcy,dzWcz;
- //double di,dj,dk;
-
- if (cofWeightx==0){di=0;}
- if (cofWeighty==0){dj=0;}
- if (cofWeightz==0){dk=0;}
-
- double weightx=(di*di)/(di*di+dj*dj+dk*dk);
- double weighty=(dj*dj)/(di*di+dj*dj+dk*dk);
- double weightz=(dk*dk)/(di*di+dj*dj+dk*dk);
- // /* if (Xw<0.5||Xw>-0.5){weightx=0;}
- // if (Yw<0.5||Yw>-0.5){weighty=0;}
- // if (Zw<0.5||Zw>-0.5){weightz=0;}
- //*/
- if (weightx==0){Ucx=0;Vcx=0;Wcx=0;dxUcx=0;dxVcx=0;dxWcx=0;}else{
- Ucx=-Udif*weightx/(Xw*Xw*Xw-Xw/4.0)*(x1LB*x1LB*x1LB-x1LB/4.0);Vcx=-Vdif*weightx/(Xw*Xw*Xw-Xw/4.0)*(x1LB*x1LB*x1LB-x1LB/4.0);Wcx=-Wdif*weightx/(Xw*Xw*Xw-Xw/4.0)*(x1LB*x1LB*x1LB-x1LB/4.0);
- dxUcx=-Udif*weightx/(Xw*Xw*Xw-Xw/4.0)*(3.0*x1LB*x1LB-1.0/4.0);dxVcx=-Vdif*weightx/(Xw*Xw*Xw-Xw/4.0)*(3.0*x1LB*x1LB-1.0/4.0);dxWcx=-Wdif*weightx/(Xw*Xw*Xw-Xw/4.0)*(3.0*x1LB*x1LB-1.0/4.0);
- }
- if (weighty==0){Ucy=0;Vcy=0;Wcy=0;dyUcy=0;dyVcy=0;dyWcy=0;}else{
- Ucy=-Udif*weighty/(Yw*Yw*Yw-Yw/4.0)*(x2LB*x2LB*x2LB-x2LB/4.0);Vcy=-Vdif*weighty/(Yw*Yw*Yw-Yw/4.0)*(x2LB*x2LB*x2LB-x2LB/4.0);Wcy=-Wdif*weighty/(Yw*Yw*Yw-Yw/4.0)*(x2LB*x2LB*x2LB-x2LB/4.0);
- dyUcy=-Udif*weighty/(Yw*Yw*Yw-Yw/4.0)*(3.0*x2LB*x2LB-1.0/4.0);dyVcy=-Vdif*weighty/(Yw*Yw*Yw-Yw/4.0)*(3.0*x2LB*x2LB-1.0/4.0);dyWcy=-Wdif*weighty/(Yw*Yw*Yw-Yw/4.0)*(3.0*x2LB*x2LB-1.0/4.0);
- }
- if (weightz==0){Ucz=0;Vcz=0;Wcz=0;dzUcz=0;dzVcz=0;dzWcz=0;}else{
- Ucz=-Udif*weightz/(Zw*Zw*Zw-Zw/4.0)*(x3LB*x3LB*x3LB-x3LB/4.0);Vcz=-Vdif*weightz/(Zw*Zw*Zw-Zw/4.0)*(x3LB*x3LB*x3LB-x3LB/4.0);Wcz=-Wdif*weightz/(Zw*Zw*Zw-Zw/4.0)*(x3LB*x3LB*x3LB-x3LB/4.0);
- dzUcz=-Udif*weightz/(Zw*Zw*Zw-Zw/4.0)*(3.0*x3LB*x3LB-1.0/4.0);dzVcz=-Vdif*weightz/(Zw*Zw*Zw-Zw/4.0)*(3.0*x3LB*x3LB-1.0/4.0);dzWcz=-Wdif*weightz/(Zw*Zw*Zw-Zw/4.0)*(3.0*x3LB*x3LB-1.0/4.0);
- }
-
- double Ucor=Ucx+Ucy+Ucz;
- double Vcor=Vcx+Vcy+Vcz;
- double Wcor=Wcx+Wcy+Wcz;
-
- double tauxxC=dxUcx;
- double tauyyC=dyVcy;
- double tauzzC=dzWcz;
- double tauxyC=0.5*(dyUcy+dxVcx);
- double tauxzC=0.5*(dzUcz+dxWcx);
- double tauyzC=0.5*(dzVcz+dyWcy);
-
- if (ii!=3)
- {
- //std::cout<<"there are not 3point for making plane"<<endl<<"Ucor="<<Ucor<<endl<<"vx1="<<vx1;
- Ucor =0.0;
- Vcor =0.0;
- Wcor =0.0;
- tauxxC=0.0;
- tauyyC=0.0;
- tauzzC=0.0;
- tauxyC=0.0;
- tauxzC=0.0;
- tauyzC=0.0;
- }
-
- vx1+=Ucor;
- vx2+=Vcor;
- vx3+=Wcor;
- tauxx+=tauxxC;
- tauyy+=tauyyC;
- tauzz+=tauzzC;
- tauxy+=tauxyC;
- tauxz+=tauxzC;
- tauyz+=tauyzC;
-
- vx1 = vx1 * conv->getFactorVelocityLbToW();
- vx2 = vx2 * conv->getFactorVelocityLbToW();
- vx3 = vx3 * conv->getFactorVelocityLbToW();
- tauxx=tauxx/dx;
- tauyy=tauyy/dx;
- tauzz=tauzz/dx;
- tauxy=tauxy/dx;
- tauxz=tauxz/dx;
- tauyz=tauyz/dx;
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::collWall(double A,double B,double C,double D,double &x1,double &x2,double &x3,double x1old,double x2old,double x3old,double dx,double &vx1,double &vx2,double &vx3,double ii)
-{
- if (particleHasMass)
- {
- double s = A*x1 + B*x2 + C*x3 + D;//The sign of s = Ax + By + Cz + D determines which side the point (x,y,z) lies with respect to the plane. If s > 0 then the point lies on the same side as the normal (A,B,C). If s < 0 then it lies on the opposite side, if s = 0 then the point (x,y,z) lies on the plane.
- if (s>0){s=1;} else if (s<0){s=-1;}else {s=0;}
- double normalDis=((A*x1 + B*x2 + C*x3 + D)/sqrt(A*A+B*B+C*C));///distance point to plane xp-Xw=distance
- double di=A/sqrt(A*A+B*B+C*C); double dj=B/sqrt(A*A+B*B+C*C); double dk=C/sqrt(A*A+B*B+C*C);
- A*=s;B*=s;C*=s;
- if (abs(normalDis)<0.05*dx)
- {
- if (ii==3)
- {
- ///reflection without e factor
- /* double vn=(A*vx1 + B*vx2 + C*vx3)/(A*A+B*B+C*C);
- double vx1L=vx1-2*vn*A;
- double vx2L=vx2-2*vn*B;
- double vx3L=vx3-2*vn*C;*/
- double vnx=A*(A*vx1 + B*vx2 + C*vx3)/(A*A+B*B+C*C);
- double vny=B*(A*vx1 + B*vx2 + C*vx3)/(A*A+B*B+C*C);
- double vnz=C*(A*vx1 + B*vx2 + C*vx3)/(A*A+B*B+C*C);
- ////do collision
- double CollitionFactor_n=0.01;
- double CollitionFactor_bt=0.01;
- vnx=-1*CollitionFactor_n*vnx;
- vny=-1*CollitionFactor_n*vny;
- vnz=-1*CollitionFactor_n*vnz;
- double vbtx=(B*B*vx1 + C*C*vx1 - A*B*vx2 - A*C*vx3)/(A*A+B*B+C*C);
- double vbty=(-(A*B*vx1) + A*A*vx2 + C*C*vx2 - B*C*vx3)/(A*A+B*B+C*C);
- double vbtz=(-(A*C*vx1) - B*C*vx2 + A*A*vx3 + B*B*vx3)/(A*A+B*B+C*C);
- vbtx=vbtx*CollitionFactor_bt;
- vbty=vbty*CollitionFactor_bt;
- vbtz=vbtz*CollitionFactor_bt;
- vx1=vnx+vbtx;
- vx2=vny+vbty;
- vx3=vnz+vbtz;
-
- x1 = (x1old);
- x2 = (x2old);
- x3 = (x3old);
- }
- else std::cout<<"there is no plane to reflect";
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::CalcVelParticle(double dx,double &vx1,double &vx2,double &vx3,double vx1old,double vx2old,
- double vx3old,double &vx1oldf,double &vx2oldf,double &vx3oldf)
-{
- double signDragx,signDragy,signDragz;
- double dragCoff=1.0;
- if (vx1old>=vx1oldf) { signDragx=-1.0; }
- else if(vx1old<vx1oldf) { signDragx= 1.0; }
- if (vx2old>=vx2oldf) { signDragy=-1.0; }
- else if(vx2old< vx2oldf) { signDragy= 1.0; }
- if (vx3old>=vx3oldf) { signDragz=-1.0; }
- else if(vx3old< vx3oldf) { signDragz= 1.0; }
- double vx1p=vx1old+signDragx*(vx1old-vx1oldf)*(vx1old-vx1oldf)*conv->getFactorTimeLbToW(dx);
- double vx2p=vx2old+signDragy*(vx2old-vx2oldf)*(vx2old-vx2oldf)*conv->getFactorTimeLbToW(dx);
- double vx3p=vx3old+signDragz*(vx3old-vx3oldf)*(vx3old-vx3oldf)*conv->getFactorTimeLbToW(dx);
-
- double signDragNextx,signDragNexty,signDragNextz;
- if (vx1p>=vx1) { signDragNextx=-1.0; }
- else if(vx1p< vx1) { signDragNextx= 1.0; }
- if (vx2p>=vx2) { signDragNexty=-1.0; }
- else if(vx2p< vx2) { signDragNexty= 1.0; }
- if (vx3p>=vx3) { signDragNextz=-1.0; }
- else if(vx3p< vx3) { signDragNextz= 1.0; }
- ////////////////////velocity of particle//////////////////////////////////////////////////////
- double velx1oldf=vx1;
- double velx2oldf=vx2;
- double velx3oldf=vx3;
-
- vx1=vx1old+1.0/2.0*( signDragNextx*(vx1p-vx1)*(vx1p-vx1)+signDragx*(vx1old-vx1oldf)*(vx1old-vx1oldf) )*conv->getFactorTimeLbToW(dx);
- vx2=vx2old+1.0/2.0*( signDragNexty*(vx2p-vx2)*(vx2p-vx2)+signDragy*(vx2old-vx2oldf)*(vx2old-vx2oldf) )*conv->getFactorTimeLbToW(dx);
- vx3=vx3old+1.0/2.0*( signDragNextz*(vx3p-vx3)*(vx3p-vx3)+signDragz*(vx3old-vx3oldf)*(vx3old-vx3oldf) )*conv->getFactorTimeLbToW(dx);
-
- vx1oldf=velx1oldf;
- vx2oldf=velx2oldf;
- vx3oldf=velx3oldf;
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::CalcVelParticle2(double dx,double &vx1,double &vx2,double &vx3,double vx1old,double vx2old,
- double vx3old,double &vx1oldf,double &vx2oldf,double &vx3oldf)
-{
- double signDragx,signDragy,signDragz;
-
- LBMReal muRE = 1.002*1e-3;//m2/s
- double mass=2.3*1e-17;// kg
- double Diameter=600*1e-9;//m
- double deltatime=conv->getFactorTimeLbToW(dx)*0.001;
- double Coff=3*PI*Diameter*muRE*deltatime/mass;
- double exCoff=exp(-Coff);
-
-
- //////////////////////velocity of particle//////////////////////////////////////////////////////
- double velx1oldf=vx1;
- double velx2oldf=vx2;
- double velx3oldf=vx3;
-
- //air
-
- vx1=vx1old*exCoff+(velx1oldf-(velx1oldf)*exCoff);
- vx2=vx2old*exCoff+(velx2oldf-(velx2oldf)*exCoff);
- vx3=vx3old*exCoff+(velx3oldf-(velx3oldf)*exCoff);
-
-
- vx1oldf=velx1oldf;
- vx2oldf=velx2oldf;
- vx3oldf=velx3oldf;
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::getRankWithPosition(Grid3DPtr grid,double xp000, double yp000, double zp000,int &RankNeigber,int level)
-{
- int minInitLevel = this->grid->getCoarsestInitializedLevel();
- int maxInitLevel = this->grid->getFinestInitializedLevel();
- bool goinside=true;
-
- LBMKernelETD3Q27Ptr kernel;
- DistributionArray3DPtr distributions;
- BCArray3D<D3Q27BoundaryCondition> bcArray;
- Block3DPtr block;
-
- //for(int level = minInitLevel; level<=maxInitLevel; level++)
- {
- UbTupleInt3 blockIndexes = grid->getBlockIndexes(xp000,yp000, zp000,level);
- block = grid->getBlock(val<1>(blockIndexes), val<2>(blockIndexes), val<3>(blockIndexes), level);
- //LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- // if(block )
- {
- if(block->isActive())
- {
- // if(!checkNodes2(block,xp000,yp000,zp000))continue;
- RankNeigber= block->getRank();
- //break;
- }
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::getRankWithPosition2(Grid3DPtr grid,double xp000, double yp000, double zp000,int &RankNeigber)
-{
- int minInitLevel = this->grid->getCoarsestInitializedLevel();
- int maxInitLevel = this->grid->getFinestInitializedLevel();
- bool goinside=true;
-
- LBMKernelETD3Q27Ptr kernel;
- DistributionArray3DPtr distributions;
- BCArray3D<D3Q27BoundaryCondition> bcArray;
- Block3DPtr block;
-
- for(int level = minInitLevel; level<=maxInitLevel; level++)
- {
- UbTupleInt3 blockIndexes = grid->getBlockIndexes(xp000,yp000, zp000,level);
- block = grid->getBlock(val<1>(blockIndexes), val<2>(blockIndexes), val<3>(blockIndexes), level);
- LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- // if(kernel )
- {
- if(block->isActive())
- {
- // if(!checkNodes2(block,xp000,yp000,zp000))continue;
- RankNeigber= block->getRank();
- break;
- }
- }
- }
-}
-
-
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::getAllTag(int allRank[7],int allTag[7])
-{
-
- for (int i=0;i<7;i++)
- {
- int tagme=0;
- for (int j=0;j<=i;j++)
- {
- if(allRank[j]==allRank[i])
- {
- tagme++;
- }
- allTag[i]=tagme;
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::getAllInfoForCompare(double input[],double *Vect)
-{
- //{0=xp000, 1=ix2ph, 2=ix3ph,3=x1,4=x2,5=x3,6=Xw,7=Yw,8=Zw,9=A,10=B,11=C,12=D,13=normalDis,14=di,15=dj,16=dk,
- //17=cofWeightx,18=cofWeighty,19=cofWeightz,,20=dx,21=omega,22=level,23=disx,24=dir };
-
- double deltaT;
- Block3DPtr block;
- finddisPointToCentercell(input[3],input[4],input[5],input[0],input[1],input[2],input[23],input[20],deltaT,block,(int)input[22]);//get double &dis,double &dx,double &deltaT
- getWallPosition(input); //get wall position,Xw,Yw,Zw and cofWeightx,cofWeighty,cofWeightz
-
- double xp000=input[0];
- double yp000=input[1];
- double zp000=input[2];
- // double nue=input[22];
- LBMReal omega = LBMSystem::calcCollisionFactor(nue,block->getLevel());
- input[21]=omega;
- UbTupleInt3 nodeIndexes000 = grid->getNodeIndexes(block, xp000, yp000, zp000);
- double xp000id = val<1>(nodeIndexes000);
- double yp000id = val<2>(nodeIndexes000);
- double zp000id = val<3>(nodeIndexes000);
- LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- if (kernel)
- {
- DistributionArray3DPtr distributions = kernel->getDataSet()->getFdistributions();
- D3Q27ICell iCell;
- iProcessor->readICell(distributions, iCell, (int)xp000id, (int)yp000id, (int)zp000id);
- getVectfromiCell(iCell,Vect);
- }
- else
- {
- cout<<__FILE__<<" "<<__LINE__<<endl;
- for (int k=0;k<8*27;k++) { Vect[k]=999.0; }
- }
-
-
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::sendMaxtagToAllProccess(int allRank[],int allTag[],int numberpro,int &numberproccessused,int *mainproccess)
-{
-
- int allRank_[7];
- for (int i=0;i<7;i++)
- {
- allRank_[i]=allRank[i];
- }
- rearangeInt(allRank_);
-
- for (int i=0;i<7;i++)
- {
- int maxtag_=0;
- if (i==6)
- {
- if (allRank_[i]!=this->root)
- {
- for (int k=0;k<=i;k++)
- {
- if (allRank_[i]==allRank_[k])
- {
- maxtag_++;
- }
- }
- MPI_Send(&maxtag_,1, MPI_INT,allRank_[i],i,MPI_COMM_WORLD);
- }
- numberproccessused++;
- mainproccess[numberproccessused-1]=allRank_[i];
- }
- else
- {
- if(allRank_[i]!=allRank_[i+1])
- {
- if (allRank_[i]!=this->root)
- {
- for (int k=0;k<=i;k++)
- {
- if (allRank_[i]==allRank_[k])
- {
- maxtag_++;
- }
- }
- MPI_Send(&maxtag_,1, MPI_INT,allRank_[i],i,MPI_COMM_WORLD);
- }
- numberproccessused++;
- mainproccess[numberproccessused-1]=allRank_[i];
- }
- }
- }
- for (int i=0;i<numberpro;i++)
- {
- bool goinside=true;
- for (int j=0;j<numberproccessused;j++)
- {
- if (i==mainproccess[j]){goinside=false;break;}
- }
- if(goinside)
- {
- int maxtag_=0;
- if (i!=this->root)
- {
- MPI_Send(&maxtag_,1, MPI_INT,i,i,MPI_COMM_WORLD);
- }
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::getAllRankWithAllPositions(double ix1ph,double ix2ph,double ix3ph,double xp000,double yp000,double zp000,
- int &Rankx,int &Ranky,int &Rankz,int &Rankxy,int &Rankxz,int &Rankyz,int &Rankxyz,int level)
-{
-
- getRankWithPosition(grid,xp000, ix2ph, ix3ph,Rankx ,level);
- getRankWithPosition(grid,ix1ph, yp000, ix3ph,Ranky ,level);
- getRankWithPosition(grid,ix1ph, ix2ph, zp000,Rankz ,level);
- getRankWithPosition(grid,xp000, yp000, ix3ph,Rankxy ,level);
- getRankWithPosition(grid,xp000, ix2ph, zp000,Rankxz ,level);
- getRankWithPosition(grid,ix1ph, yp000, zp000,Rankyz ,level);
- getRankWithPosition(grid,xp000, yp000, zp000,Rankxyz,level);
-
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::getAllRankWithAllPositions2(double ix1ph,double ix2ph,double ix3ph,double xp000,double yp000,double zp000,
- int &Rankx,int &Ranky,int &Rankz,int &Rankxy,int &Rankxz,int &Rankyz,int &Rankxyz)
-{
-
- getRankWithPosition2(grid,xp000, ix2ph, ix3ph,Rankx );
- getRankWithPosition2(grid,ix1ph, yp000, ix3ph,Ranky );
- getRankWithPosition2(grid,ix1ph, ix2ph, zp000,Rankz );
- getRankWithPosition2(grid,xp000, yp000, ix3ph,Rankxy );
- getRankWithPosition2(grid,xp000, ix2ph, zp000,Rankxz );
- getRankWithPosition2(grid,ix1ph, yp000, zp000,Rankyz );
- getRankWithPosition2(grid,xp000, yp000, zp000,Rankxyz);
-
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::getVectfromiCell(D3Q27ICell iCell,double *Vect)
-{
-
- for (int i=0*27;i<1*27;i++) { Vect[i]=iCell.TSW[i-0*27]; }
- for (int i=1*27;i<2*27;i++) { Vect[i]=iCell.TNW[i-1*27]; }
- for (int i=2*27;i<3*27;i++) { Vect[i]=iCell.TNE[i-2*27]; }
- for (int i=3*27;i<4*27;i++) { Vect[i]=iCell.TSE[i-3*27]; }
- for (int i=4*27;i<5*27;i++) { Vect[i]=iCell.BSW[i-4*27]; }
- for (int i=5*27;i<6*27;i++) { Vect[i]=iCell.BNW[i-5*27]; }
- for (int i=6*27;i<7*27;i++) { Vect[i]=iCell.BNE[i-6*27]; }
- for (int i=7*27;i<8*27;i++) { Vect[i]=iCell.BSE[i-7*27]; }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::getiCellfromVect(D3Q27ICell &iCell,double *Vect)
-{
-
- for (int i=0*27;i<1*27;i++) { iCell.TSW[i-0*27]=Vect[i]; }
- for (int i=1*27;i<2*27;i++) { iCell.TNW[i-1*27]=Vect[i]; }
- for (int i=2*27;i<3*27;i++) { iCell.TNE[i-2*27]=Vect[i]; }
- for (int i=3*27;i<4*27;i++) { iCell.TSE[i-3*27]=Vect[i]; }
- for (int i=4*27;i<5*27;i++) { iCell.BSW[i-4*27]=Vect[i]; }
- for (int i=5*27;i<6*27;i++) { iCell.BNW[i-5*27]=Vect[i]; }
- for (int i=6*27;i<7*27;i++) { iCell.BNE[i-6*27]=Vect[i]; }
- for (int i=7*27;i<8*27;i++) { iCell.BSE[i-7*27]=Vect[i]; }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::getAfterCompare(double dis[],double input [7][25],double AllCell[7][8*27],double &xp000,double &yp000,double &zp000,
- double &Xw,double &Yw,double &Zw,double &cofWeightx,double &cofWeighty,double &cofWeightz,double &dx,double &omega,D3Q27ICell &iCell)
-{
-
- double *Cell=(double*)malloc(8*27*sizeof(double));
- bool goinside=true;
-
- for (int i=0;i<7;i++)
- {
- double dist=dis[i];
- if(goinside )
- { //17=cofWeightx,18=cofWeighty,19=cofWeightz,,20=dx,21=omega
- if (dist==input[0][23]){ Xw=input[0][6]; Yw=input[0][7]; Zw=input[0][8]; xp000=input[0][0]; yp000=input[0][1]; zp000=input[0][2]; cofWeightx=input[0][17]; cofWeighty=input[0][18]; cofWeightz=input[0][19]; dx=input[0][20]; omega=input[0][21]; for (int i=0;i<8*27;i++){Cell[i]=AllCell[0][i];} } //get suitable block from suitable cell
- else if (dist==input[1][23]){ Xw=input[1][6]; Yw=input[1][7]; Zw=input[1][8]; xp000=input[1][0]; yp000=input[1][1]; zp000=input[1][2]; cofWeightx=input[1][17]; cofWeighty=input[1][18]; cofWeightz=input[1][19]; dx=input[1][20]; omega=input[1][21]; for (int i=0;i<8*27;i++){Cell[i]=AllCell[1][i];} }
- else if (dist==input[2][23]){ Xw=input[2][6]; Yw=input[2][7]; Zw=input[2][8]; xp000=input[2][0]; yp000=input[2][1]; zp000=input[2][2]; cofWeightx=input[2][17]; cofWeighty=input[2][18]; cofWeightz=input[2][19]; dx=input[2][20]; omega=input[2][21]; for (int i=0;i<8*27;i++){Cell[i]=AllCell[2][i];} }
- else if (dist==input[3][23]){ Xw=input[3][6]; Yw=input[3][7]; Zw=input[3][8]; xp000=input[3][0]; yp000=input[3][1]; zp000=input[3][2]; cofWeightx=input[3][17]; cofWeighty=input[3][18]; cofWeightz=input[3][19]; dx=input[3][20]; omega=input[3][21]; for (int i=0;i<8*27;i++){Cell[i]=AllCell[3][i];} }
- else if (dist==input[4][23]){ Xw=input[4][6]; Yw=input[4][7]; Zw=input[4][8]; xp000=input[4][0]; yp000=input[4][1]; zp000=input[4][2]; cofWeightx=input[4][17]; cofWeighty=input[4][18]; cofWeightz=input[4][19]; dx=input[4][20]; omega=input[4][21]; for (int i=0;i<8*27;i++){Cell[i]=AllCell[4][i];} }
- else if (dist==input[5][23]){ Xw=input[5][6]; Yw=input[5][7]; Zw=input[5][8]; xp000=input[5][0]; yp000=input[5][1]; zp000=input[5][2]; cofWeightx=input[5][17]; cofWeighty=input[5][18]; cofWeightz=input[5][19]; dx=input[5][20]; omega=input[5][21]; for (int i=0;i<8*27;i++){Cell[i]=AllCell[5][i];} }
- else if (dist==input[6][23]){ Xw=input[6][6]; Yw=input[6][7]; Zw=input[6][8]; xp000=input[6][0]; yp000=input[6][1]; zp000=input[6][2]; cofWeightx=input[6][17]; cofWeighty=input[6][18]; cofWeightz=input[6][19]; dx=input[6][20]; omega=input[6][21]; for (int i=0;i<8*27;i++){Cell[i]=AllCell[6][i];} }
- if (Xw<2.5&&Xw>-2.5&&Yw<2.50&&Yw>-2.50&&Zw<2.50&&Zw>-2.50){goinside=false;}
- }
- }
-
- getiCellfromVect(iCell,Cell);
- free (Cell);
-}
-
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::SendAndReceiveData(double input[7][25],double AllCell[7][8*27],int allRank[7],int allTag[7])
-{
-
- int sizeinput =25;
- double *Cell=(double*)malloc(8*27*sizeof(double));
-
-
- for (int i=0;i<7;i++)
- {
- if (allRank[i]==this->root)
- {
- double vectforTrans[25];
- for (int j=0;j<sizeinput;j++)
- {
- vectforTrans[j]=input[i][j];
- }
- getAllInfoForCompare(vectforTrans,Cell);
- int index=(int)vectforTrans[24];
- for (int k=0;k<8*27;k++)
- {
- AllCell[index][k]=Cell[k];
- }
- for (int j=0;j<sizeinput;j++)
- {
- input[index][j]=vectforTrans[j];
- }
- }
- else
- {
- double vectforTrans[25];
- for (int j=0;j<sizeinput;j++)
- {
- vectforTrans[j]=input[i][j];
- }
- MPI_Send(vectforTrans,sizeinput, MPI_DOUBLE_PRECISION,allRank[i],allTag[i],MPI_COMM_WORLD);
- }
- }
- for (int i=0;i<7;i++)
- {
- if (allRank[i]!=this->root)
- {
- MPI_Status status;
- double receiveVector[25];
- double *receiceCell=(double*)malloc(8*27*sizeof(double));
- MPI_Recv(receiveVector,25, MPI_DOUBLE_PRECISION,allRank[i],MPI_ANY_TAG,MPI_COMM_WORLD,&status);
- MPI_Recv(receiceCell,8*27, MPI_DOUBLE_PRECISION,allRank[i],MPI_ANY_TAG,MPI_COMM_WORLD,&status);
- int index=(int)receiveVector[24];
- for (int j=0;j<sizeinput;j++)
- {
- input[index][j]=receiveVector[j];
- }
- for (int k=0;k<8*27;k++)
- {
- AllCell[index][k]=receiceCell[k];
- }
- free (receiceCell);
- }
- }
- free (Cell);
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessor::sendtoOtherIsOk(int isPointSuitable)
-{
- if(comm->getNumberOfProcesses() > 1)
- {
- int size=comm->getNumberOfProcesses();
- for (int i=0;i<size;i++)
- {
- if (i!=rank)
- {
- MPI_Send(&isPointSuitable,1, MPI_INT,i,i,MPI_COMM_WORLD);
- }
- }
- }
-}
diff --git a/depot/PathLineCoProcessor.h b/depot/PathLineCoProcessor.h
deleted file mode 100644
index 9677ab3cf..000000000
--- a/depot/PathLineCoProcessor.h
+++ /dev/null
@@ -1,115 +0,0 @@
-/*
-* PathLinePostprocessor.h
-*
-* Created on: 04.24.2012
-* Author: Fard
-*/
-
-
-#ifndef PATHLINEPOSTPROCESSOR_H
-#define PATHLINEPOSTPROCESSOR_H
-
-# define dirx 0
-# define diry 1
-# define dirz 2
-# define dirxy 3
-# define dirxz 4
-# define diryz 5
-# define dirxyz 6
-#include <mpi.h>
-#include "CoProcessor.h"
-#include "Grid3D.h"
-#include "Block3D.h"
-#include "LBMUnitConverter.h"
-#include "Communicator.h"
-#include "D3Q27InterpolationProcessor.h"
-#include "LBMKernelETD3Q27.h"
-#include "D3Q27ETBCProcessor.h"
-#include <boost/shared_ptr.hpp>
-#include "WbWriter.h"
-
-class PathLineCoProcessor;
-typedef boost::shared_ptr<PathLineCoProcessor> PathLineCoProcessorPtr;
-
-class PathLineCoProcessor : public CoProcessor
-{
-public:
- PathLineCoProcessor(Grid3DPtr grid, const std::string& path, WbWriter* const writer,
- LBMUnitConverterPtr conv, UbSchedulerPtr s, CommunicatorPtr comm,
- double x1, double x2, double x3, LBMReal nue, D3Q27InterpolationProcessorPtr iProcessor);
- ~PathLineCoProcessor();
- void process(double step);
-
-protected:
- void collectPostprocessData();
- void initialMovement(Block3DPtr block, LBMKernelETD3Q27Ptr& kernel, DistributionArray3DPtr& distributions, BCArray3D<D3Q27BoundaryCondition>& bcArray);
- void interpolMovement(Block3DPtr block, LBMKernelETD3Q27Ptr& kernel, DistributionArray3DPtr& distributions, BCArray3D<D3Q27BoundaryCondition>& bcArray,int &isExtrapolation,double tau[]);
- void extrapolMovement(Block3DPtr block, LBMKernelETD3Q27Ptr& kernel, DistributionArray3DPtr& distributions, BCArray3D<D3Q27BoundaryCondition>& bcArray,double tau[]);
- void clearData();
- void outICell(D3Q27ICell& iCell);
- void updateDistributedValues();
- void printPoint(double tau[]);
- void findPlane(int ix1,int ix2,int ix3,Grid3DPtr grid,Block3DPtr block,double &A,double &B,double &C,double &D,double &ii);
- void finddisPointToCentercell(double x1,double x2,double x3,double xp000, double yp000, double zp000,double &dis,double &dx,double &deltaT,Block3DPtr &block,int level);
-
- //void getWallPosition(Grid3DPtr grid,Block3DPtr &block,double xp000,double yp000,double zp000,double &Xw,double &Yw,double &Zw,double x1,double x2,double x3,
- // double &xp000T,double &yp000T,double &zp000T,double A,double B,double C,double D,double normalDis,double di,double dj,double dk,double &cofWeightx
- //,double &cofWeighty,double &cofWeightz,double dismax,double dis,int level);
- void getWallPosition(double input[]);
- void findInitialCell(double s,double di,double dj,double dk,double dx,double ix1ph,double ix2ph,double ix3ph,double &xp000,double &yp000,double &zp000);
- void rearangeInt( int dis[7]);
- void rearangedDouble( double dis[7]);
- void addCorrection(double &vx1,double &vx2,double &vx3,double &tauxx,double &tauyy,double &tauzz,double &tauxy,double &tauxz,double &tauyz,
- double dx,D3Q27ICell iCell,double Xw, double Yw, double Zw,double omega,double cofWeightx,double cofWeighty,double cofWeightz,double ii,double x1LB,double x2LB,double x3LB,double di,double dj,double dk);
- void collWall(double A,double B,double C,double D,double &x1,double &x2,double &x3,double x1old,double x2old,double x3old,double dx,double &vx1,double &vx2,double &vx3,double ii);
- void CalcVelParticle(double dx,double &vx1,double &vx2,double &vx3,double vx1old,double vx2old,double vx3old,double &vx1oldf,double &vx2oldf,double &vx3oldf);
- void CalcVelParticle2(double dx,double &vx1,double &vx2,double &vx3,double vx1old,double vx2old,double vx3old,double &vx1oldf,double &vx2oldf,double &vx3oldf);
-
- // void getWallPosition(double dixmax,double disx,double disy,double disz,double disxy,double disxz,double disyz,double disxyz,double &Xw,double &Yw,double &Zw,double x1,double x2,double x3,double A,double B,double C,double D,double normalDis,double di,double dj,double dk,
- // double &cofWeightx,double &cofWeighty,double &cofWeightz);
- void getRankWithPosition(Grid3DPtr grid,double xp000, double yp000, double zp000,int &RankNeigber,int level);
- void getRankWithPosition(Grid3DPtr grid,double xp000, double yp000, double zp000,int &RankNeigber);
- void getRankWithPosition2(Grid3DPtr grid,double xp000, double yp000, double zp000,int &RankNeigber);
- void getAllTag(int allRank[7],int allTag[7]);
- void getAllInfoForCompare(double input[],double *Vect);
- void sendMaxtagToAllProccess(int allRank[],int allTag[],int numberpro,int &numberproccessused,int *mainproccess);
- void getAllRankWithAllPositions(double ix1ph,double ix2ph,double ix3ph,double xp000,double yp000,double zp000,int &Rankx,int &Ranky,int &Rankz,int &Rankxy,int &Rankxz,int &Rankyz,int &Rankxyz,int level);
- void getAllRankWithAllPositions2(double ix1ph,double ix2ph,double ix3ph,double xp000,double yp000,double zp000,int &Rankx,int &Ranky,int &Rankz,int &Rankxy,int &Rankxz,int &Rankyz,int &Rankxyz);
- void getVectfromiCell(D3Q27ICell iCell,double *Vect);
- void getiCellfromVect(D3Q27ICell &iCell,double *Vect);
-
- void getAfterCompare(double dis[],double input [7][25],double AllCell[7][8*27],double &xp000,double &yp000,double &zp000,
- double &Xw,double &Yw,double &Zw,double &cofWeightx,double &cofWeighty,double &cofWeightz,double &dx,double &omega,D3Q27ICell &iCell);
- void SendAndReceiveData(double input [7][25],double AllCell[7][8*27],int allRank[7],int allTag[7]);
- void sendtoOtherIsOk(int isPointSuitable);
-private:
- bool checkNodes(Block3DPtr block);
- bool checkNodes2( Block3DPtr block,double x11,double x22,double x33);
- void checkLevel(Block3DPtr& block, LBMKernelETD3Q27Ptr& kernel, DistributionArray3DPtr& distributions, BCArray3D<D3Q27BoundaryCondition>& bcArray,bool &isPointSuitable2);
-
- std::vector<UbTupleDouble3> nodes;
- std::vector<std::string> datanames;
- std::vector<std::vector<double> > data;
- std::string path;
- WbWriter* writer;
- LBMUnitConverterPtr conv;
- CommunicatorPtr comm;
- D3Q27InterpolationProcessorPtr interpolation;
- double x1, x2, x3;
- double x1old, x2old, x3old;
- LBMReal vx1old,vx2old,vx3old;
- LBMReal vx1oldf,vx2oldf,vx3oldf;
- bool particleHasMass;
- int isExtrapolation;
- int istep;
- int stepcheck;
- LBMReal nue;
- D3Q27InterpolationProcessorPtr iProcessor;
- D3Q27InterpolationHelperPtr iHelper;
- int rank;
- int root;
- int maxtag;
- std::vector<std::vector<double> > Positions;
-};
-#endif
-
diff --git a/depot/PathLineCoProcessorMcpart.cpp b/depot/PathLineCoProcessorMcpart.cpp
deleted file mode 100644
index e3acc7b3e..000000000
--- a/depot/PathLineCoProcessorMcpart.cpp
+++ /dev/null
@@ -1,1301 +0,0 @@
-#include "PathLineCoProcessorMcpart.h"
-#include "LBMKernelETD3Q27.h"
-#include "SimulationParameters.h"
-#include "D3Q27ETBCProcessor.h"
-#include <vector>
-#include <string>
-#include <boost/foreach.hpp>
-#include "basics/writer/WbWriterVtkXmlASCII.h"
-//#include "D3Q27ETFCoffVectorConnector.h"
-using namespace std;
-
-
-PathLineCoProcessorMcpart::PathLineCoProcessorMcpart(Grid3DPtr grid, const std::string& path, WbWriter* const writer,
- LBMUnitConverterPtr conv, UbSchedulerPtr s, CommunicatorPtr comm,
- std::vector<UbTupleDouble3 > _Positions, LBMReal nue, D3Q27InterpolationProcessorPtr iProcessor)
- : CoProcessor(grid, s),
- path(path),
- comm(comm),
- writer(writer),
- conv(conv),
- nue(nue),
- iProcessor(iProcessor),
- istep(0),
- particleHasMass(true),
- rank(comm->getProcessID())
-{
- iHelper = D3Q27InterpolationHelperPtr(new D3Q27InterpolationHelper(iProcessor));
- getParticlePerProccess(_Positions);
- getNeighborsRank();
- initializeForPrinting((int)_Positions.size());
-}
-//////////////////////////////////////////////////////////////////////////
-PathLineCoProcessorMcpart::~PathLineCoProcessorMcpart()
-{
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::process(double step)
-{
- this->stepcheck=(int)step;
- double aa=scheduler->getMinEnd();
- if (step >= scheduler->getMaxBegin() && step <= scheduler->getMinEnd())
- {
- collectPostprocessData();
- }
-
- UBLOG(logDEBUG3, "PathLinePostprocessor::update:" << step);
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::collectPostprocessData()
-{
-
- int ran=rank;
- updateParticles();
- checkParticles();
- initialMovement();
- updateParticles();
- checkParticles();
- BOOST_FOREACH(ParticlesPtr particle, particles)
- {
- double tau[6];//={tauxx,tauyy,tauzz,tauxy,tauxz,tauyz};
- finalMovement(particle,tau);
- //printParticle(particle);
- gatherData(particle);
- }
-
- if(scheduler->isDue((double)istep) )
- {
- string partName = writer->writeNodesWithNodeDataDouble(path+ UbSystem::toString(rank),nodes,datanames,data);
- size_t found=partName.find_last_of("//");
- string piece = partName.substr(found+1);
-
- vector<string> cellDataNames;
-
- vector<string> pieces = comm->gather(piece);
- if (comm->getProcessID() == comm->getRoot())
- {
- string pname = WbWriterVtkXmlASCII::getInstance()->writeParallelFile(path,pieces,datanames,cellDataNames);
-
- //vector<string> filenames;
- //filenames.push_back(pname);
- //if (istep == Postprocessor::scheduler->getMinBegin())
- //{
- // WbWriterVtkXmlASCII::getInstance()->writeCollection(path+"_collection",filenames,0,false);
- //}
- //else
- //{
- // WbWriterVtkXmlASCII::getInstance()->addFilesToCollection(path+"_collection",filenames,0,false);
- //}
- UBLOG(logINFO,"D3Q27MacroscopicQuantitiesPostprocessor step: " << istep);
- }
- }
-
- istep++;
-}
-////////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::initializeForPrinting(int _number)
-{
- for (int i=0;i<_number;i++)
- {
- string fileName=path + UbSystem::toString(i)+ ".txt";
- files.push_back(boost::shared_ptr<std::ofstream>(new std::ofstream(fileName.c_str())));
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::getParticlePerProccess(std::vector<UbTupleDouble3 >particlesVec)
-{
- for (int i=0;i<particlesVec.size();i++)
- {
- int minInitLevel = this->grid->getCoarsestInitializedLevel();
- int maxInitLevel = this->grid->getFinestInitializedLevel();
- LBMKernelETD3Q27Ptr kernel;
- DistributionArray3DPtr distributions;
- BCArray3D<D3Q27BoundaryCondition> bcArray;
- Block3DPtr block;
-
- for(int level = minInitLevel; level<=maxInitLevel; level++)
- {
- UbTupleInt3 blockIndexes = grid->getBlockIndexes(val<1>(particlesVec[i]), val<2>(particlesVec[i]), val<3>(particlesVec[i]),level);
- block = grid->getBlock(val<1>(blockIndexes), val<2>(blockIndexes), val<3>(blockIndexes), level);
- if(!block) continue;
- if(block->isActive())
- {
- if (block->hasInterpolationFlagCF() || block->hasInterpolationFlagFC())
- {
- if(rank == block->getRank())
- {
- if(!checkNodes(block,val<1>(particlesVec[i]), val<2>(particlesVec[i]), val<3>(particlesVec[i])))
- {
- sendStatusOfPoint(false,val<1>(particlesVec[i]), val<2>(particlesVec[i]), val<3>(particlesVec[i]),level,i);
- continue;
- }
- initializeParticle ( val<1>(particlesVec[i]), val<2>(particlesVec[i]), val<3>(particlesVec[i]),i ,level);
- //getParticle(particle,grid->getDeltaX(block),particle.x,particle.y,particle.z,level,i);
- sendStatusOfPoint(true,val<1>(particlesVec[i]), val<2>(particlesVec[i]), val<3>(particlesVec[i]),level,i);
- break;
- }
- else
- {
- bool isPointSuitable=true;
- receiveStatusOfPoint(isPointSuitable,block->getRank(),val<1>(particlesVec[i]), val<2>(particlesVec[i]), val<3>(particlesVec[i]),level) ;
- if (isPointSuitable){break;}
- }
- }
- else
- {
- if(rank == block->getRank())
- {
- if(!checkNodes(block,val<1>(particlesVec[i]), val<2>(particlesVec[i]), val<3>(particlesVec[i])))
- {
- continue;
- }
- initializeParticle ( val<1>(particlesVec[i]), val<2>(particlesVec[i]), val<3>(particlesVec[i]),i,level );
- Particles particle;
- //getParticle(particle,grid->getDeltaX(block),particle.x,particle.y,particle.z,level,i);
- //particlesVec.push_back(particle);
- //particles.erase(t.first);
- break;
- }
-
- }
- }
- }
- }
-}
-////////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::initializeParticle(double x, double y,double z,int i,int level)
-{
- ParticlesPtr particle = ParticlesPtr(new Particles() );
-
- particle->xold=x; particle->yold=y; particle->zold=z;
- particle->x =x; particle->y =y; particle->z =z;
- particle->vxold =0.1; particle->vyold =0.1; particle->vzold=0.1;
- particle->vxoldf=0.0; particle->vyoldf=0.0; particle->vzoldf=0.0;
- particle->rankOfParticle=rank; particle->ID=i; particle->level=level;
- particles.push_back(particle);
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::getNeighborsRank()
-{
- int gridRank = grid->getRank();
- int minInitLevel = this->grid->getCoarsestInitializedLevel();
- int maxInitLevel = this->grid->getFinestInitializedLevel();
-
- for(int level = minInitLevel; level<=maxInitLevel;level++)
- {
- vector<Block3DPtr> blockVector;
- grid->getBlocks(level, gridRank, blockVector);
- BOOST_FOREACH(Block3DPtr block, blockVector)
- {
- int blockRank = block->getRank();
- if (gridRank == blockRank && block->isActive())
- {
- //std::vector<Block3DPtr> neighbors;
- int ix1 = block->getX1();
- int ix2 = block->getX2();
- int ix3 = block->getX3();
- int level = block->getLevel();
- //grid->getAllNeighbors(ix1, ix2, ix3, level, level, neighbors);
- int dirs=26;//number of directions
- for( int dir = 0; dir < dirs; dir++)
- {
- Block3DPtr neighBlock = grid->getNeighborBlock(dir, ix1, ix2, ix3, level);
- if(neighBlock)
- {
- int neighBlockRank = neighBlock->getRank();
- if(blockRank != neighBlockRank && neighBlock->isActive())
- {
- if (neighbors.size()==0)
- {
- neighbors.insert((std::make_pair(neighBlockRank,0)));
- }
- else
- {
- for ( std::map< int, int >::const_iterator iter = neighbors.begin();iter != neighbors.end(); ++iter )
- {
- if (iter->first!=neighBlockRank)
- {
- neighbors.insert((std::make_pair(neighBlockRank,0)));
- }
- }
- }
- }
- }
- }
- }
- }
- }
-}
-
-//////////////////////////////////////////////////////////////////////////
-//void D3Q27PathLinePostprocessorMcpart::getParticle(Particles &particle, double dx,double x, double y,double z,int level,int i)
-//{
-// particle.particleNumber=i;
-// particle.xold=val<1>(xold[i]); particle.yold=val<2>(xold[i]); particle.zold=val<3>(xold[i]);
-// particle.x =x; particle.y =y; particle.z =z;
-// particle.level=level;
-// particle.dx=dx;
-// particle.vxold=val<1>(vxold[i]); particle.vyold=val<2>(vxold[i]); particle.vzold=val<3>(vxold[i]);
-// particle.vxoldf=val<1>(vxoldf[i]); particle.vyoldf=val<2>(vxoldf[i]); particle.vzoldf=val<3>(vxoldf[i]);
-//
-//}
-//////////////////////////////////////////////////////////////////////////
-bool PathLineCoProcessorMcpart::checkNodes( Block3DPtr block,double _x1,double _x2,double _x3)
-{
- bool result = true;
- LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- BCArray3D<D3Q27BoundaryCondition> bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
- D3Q27BoundaryConditionPtr bcPtr;
-
- double x1_ch = _x1;// + vx1old*conv->getfactorTimeShouldMultiplebyDx()*grid->getDeltaX(block);
- double x2_ch = _x2;// + vx2old*conv->getfactorTimeShouldMultiplebyDx()*grid->getDeltaX(block);
- double x3_ch = _x3;// + vx3old*conv->getfactorTimeShouldMultiplebyDx()*grid->getDeltaX(block);
-
- int nx1 = static_cast<int>(bcArray.getNX1());
- int nx2 = static_cast<int>(bcArray.getNX2());
- int nx3 = static_cast<int>(bcArray.getNX3());
-
- UbTupleInt3 nodeIndexes = grid->getNodeIndexes(block, x1_ch, x2_ch, x3_ch);
-
- int ix1 = val<1>(nodeIndexes);
- int ix2 = val<2>(nodeIndexes);
- int ix3 = val<3>(nodeIndexes);
- bool outlet=false;
-
- for (int xx3 = ix3; xx3 <= ix3 + 1 && ix3 + 1 < nx3; xx3++)
- for(int xx2 = ix2; xx2 <= ix2 + 1 && ix2 + 1 < nx2; xx2++)
- for(int xx1 = ix1; xx1 <= ix1 + 1 && ix1 + 1 < nx1; xx1++)
- {
- bcPtr=bcArray.getBC(xx1,xx2,xx3);
- for(int fdir=D3Q27System::STARTF; fdir<=D3Q27System::ENDF; fdir++)
- {
- if ( bcPtr != NULL)
- {
- if (bcPtr->hasDensityBoundaryFlag(fdir))
- {
- if(outlet)break;
- outlet=true;
- break;
- }
- }
- }
- result = result && !bcArray.isUndefined(xx1, xx2, xx3) && !outlet;
- }
-
- return result;
-
-}
-// //////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::initialMovement()
-{
- BOOST_FOREACH(ParticlesPtr particle, particles)
- {
- double dx=grid->getDeltaX(particle->level);
- particle->xold=particle->x;
- particle->yold=particle->y;
- particle->zold=particle->z;
- particle->x =particle->xold+particle->vxold*dx*conv->getFactorTimeLbToW(dx);
- particle->y =particle->yold+particle->vyold*dx*conv->getFactorTimeLbToW(dx);
- particle->z =particle->zold+particle->vzold*dx*conv->getFactorTimeLbToW(dx);
- }
-}
-///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-//void D3Q27PathLinePostprocessorMcpart::updateDistributedValues(std::vector<Particles> *particlesVec,std::vector<UbTupleDouble3> &x)
-//{
-//
-// int numberOFVariable=18;
-// std::vector<double> sendParticlesInfomation;
-// std::vector<double> receiveParticlesInfomation;
-// getVectorFromParticles(sendParticlesInfomation,particlesVec,numberOFVariable);
-// allGatherDoubles(sendParticlesInfomation, receiveParticlesInfomation);
-// fillOutPositions(receiveParticlesInfomation,x);
-// //getParticlesFromVector(receiveParticlesInfomation,particlesVec,numberOFVariable);
-// //std::vector<Particles> *particlesVec;
-//}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::sendBool(bool _bool,int distinationRank ,int i)
-{
- int variable;
- if (_bool==false)
- {
- variable=0;
- }
- else variable=1;
- if(comm->getNumberOfProcesses() > 1)
- {
- MPI_Send(&variable,1, MPI_INT,distinationRank,i,MPI_COMM_WORLD);
- //std::cout<<"\n"<<variable<<"send from rank "<<rank<<" to rank "<< distinationRank;
- }
-}
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::sendInt(int svalue,int distinationRank ,int tag)
-{
- if(comm->getNumberOfProcesses() > 1)
- {
- MPI_Send(&svalue,1, MPI_INT,distinationRank,tag,MPI_COMM_WORLD);
- }
-}
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::receiveBool(bool &variable,int rankBlock)
-{
- if(comm->getNumberOfProcesses() > 1)
- {
- int _bool;
- MPI_Status status;
- MPI_Recv(&_bool,1, MPI_INT,rankBlock,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
- //std::cout<<"\n"<<_bool<<"received in rank "<<rank<<" from rank "<< rankBlock;
- if (_bool==false)
- {
- variable=false;
- }
- else variable=true;
- }
-}
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::receiveInt(int &rvalue,int root)
-{
- if(comm->getNumberOfProcesses() > 1)
- {
- MPI_Status status;
- MPI_Recv(&rvalue,1, MPI_INT,root,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
- }
-}
-////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::sendStatusOfPoint(bool status,double x1,double x2,double x3,double level,int i)
-{
- for(int lev = (int)level-1; lev<=(int)level+1; lev++)
- {
- UbTupleInt3 blockIndexesDistination = grid->getBlockIndexes(x1,x2,x3,lev);
- Block3DPtr blockDistination= grid->getBlock(val<1>(blockIndexesDistination), val<2>(blockIndexesDistination), val<3>(blockIndexesDistination), lev);
- if(!blockDistination) continue;
- if(blockDistination->isActive())
- {
- if(rank!=blockDistination->getRank())
- {
- sendBool(status,blockDistination->getRank(),i );
- }
-
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::receiveStatusOfPoint(bool &status,int rankRoot,double x1,double x2,double x3,double level)
-{
- for(int lev = (int)level-1; lev<=(int)level+1; lev++)
- {
- UbTupleInt3 blockIndexes = grid->getBlockIndexes(x1,x2,x3,lev);
- Block3DPtr block= grid->getBlock(val<1>(blockIndexes), val<2>(blockIndexes), val<3>(blockIndexes), lev);
- if (block)
- {
- if(rankRoot!=block->getRank())
- {
- LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- if(kernel)
- {
- if(checkNodes(block,x1,x2,x3))
- {
- receiveBool(status,rankRoot);
- }
- }
-
- }
-
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::getVectorFromParticles(std::vector<double> &particlesInfomation,ParticlesPtr particle)
-{
- particlesInfomation.push_back(particle->x );
- particlesInfomation.push_back(particle->y );
- particlesInfomation.push_back(particle->z );
- particlesInfomation.push_back(particle->xold );
- particlesInfomation.push_back(particle->yold );
- particlesInfomation.push_back(particle->zold );
- particlesInfomation.push_back(particle->vxold );
- particlesInfomation.push_back(particle->vyold );
- particlesInfomation.push_back(particle->vzold );
- particlesInfomation.push_back(particle->vxoldf);
- particlesInfomation.push_back(particle->vyoldf);
- particlesInfomation.push_back(particle->vzoldf);
- particlesInfomation.push_back(particle->rankOfParticle);
- particlesInfomation.push_back(particle->ID);
- particlesInfomation.push_back(particle->level);
-
-}
-////////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::getParticlesFromVector(std::vector<double> particlesInfomation,int numberOFVariable)
-{
- int numParticles= (int)particlesInfomation.size()/numberOFVariable;
- for (int i=0;i<numParticles;i++)
- {
- ParticlesPtr particle = ParticlesPtr(new Particles() );
-
- particle->x =particlesInfomation[i*numberOFVariable+0];
- particle->y =particlesInfomation[i*numberOFVariable+1];
- particle->z =particlesInfomation[i*numberOFVariable+2];
- particle->xold =particlesInfomation[i*numberOFVariable+3];
- particle->yold =particlesInfomation[i*numberOFVariable+4];
- particle->zold =particlesInfomation[i*numberOFVariable+5];
- particle->vxold =particlesInfomation[i*numberOFVariable+6];
- particle->vyold =particlesInfomation[i*numberOFVariable+7];
- particle->vzold =particlesInfomation[i*numberOFVariable+8];
- particle->vxoldf =particlesInfomation[i*numberOFVariable+9];
- particle->vyoldf =particlesInfomation[i*numberOFVariable+10];
- particle->vzoldf =particlesInfomation[i*numberOFVariable+11];
- particle->rankOfParticle =(int)particlesInfomation[i*numberOFVariable+12];
- particle->ID =(int)particlesInfomation[i*numberOFVariable+13];
- particle->level =(int)particlesInfomation[i*numberOFVariable+14];
- particles.push_back(particle);
- }
-}
-////////////////////////////////////////////////////////////////////////////
-//void D3Q27PathLinePostprocessorMcpart::updateinfo(std::vector<Particles> *particlesVec,std::vector<UbTupleDouble3> &x)
-//{
-// std::vector<Particles>&particles=*particlesVec;
-// for (int i=0;i<particlesVec->size();i++)
-// {
-// Particles &particle=particles[i];
-// int index=particle.particleNumber;
-// val<1>(x[index])=particle.xold;
-// val<2>(x[index])=particle.yold;
-// val<3>(x[index])=particle.zold;
-// }
-//}
-////////////////////////////////////////////////////////////////////////////
-//void D3Q27PathLinePostprocessorMcpart::allGatherDoubles(std::vector<double>& svalues, std::vector<double>& rvalues)
-//{
-// int scount;
-// vector<int> displs, rcounts;
-//
-// scount = (int)(svalues.size());
-// rcounts.resize(comm->getNumberOfProcesses());
-// MPI_Allgather(&scount, 1, MPI_INT, &rcounts[0], 1, MPI_INT, MPI_COMM_WORLD);
-// displs.resize(comm->getNumberOfProcesses());
-//
-// displs[0] = 0;
-// for (int i=1; i<comm->getNumberOfProcesses(); ++i) {
-// displs[i] = displs[i-1]+rcounts[i-1];
-// }
-//
-// rvalues.resize(displs[comm->getNumberOfProcesses()-1]+rcounts[comm->getNumberOfProcesses()-1]);
-//
-// if(rvalues.size() == 0)
-// {
-// rvalues.resize(1);
-// rvalues[0] = -999;
-// }
-// if(scount == 0)
-// {
-// svalues.resize(1);
-// svalues[0] = -999;
-// }
-//
-// MPI_Allgatherv(&svalues[0], scount, MPI_DOUBLE, &rvalues[0], &rcounts[0], &displs[0], MPI_DOUBLE, MPI_COMM_WORLD);
-//}
-////////////////////////////////////////////////////////////////////////////
-//void D3Q27PathLinePostprocessorMcpart::fillOutPositions(std::vector<double> particlesInfomation,std::vector<UbTupleDouble3> &x,double numberOFVariable)
-//{
-// if (particlesInfomation.size()!=x.size()*numberOFVariable)
-// {
-// std::cout<<"number of particle is wrong";
-// }
-// for (int i=0;i<x.size();i++)
-// {
-// Particles particle;
-// particle.x =particlesInfomation[i*numberOFVariable+0];
-// particle.y =particlesInfomation[i*numberOFVariable+1];
-// particle.z =particlesInfomation[i*numberOFVariable+2];
-// particle.xold =particlesInfomation[i*numberOFVariable+3];
-// particle.yold =particlesInfomation[i*numberOFVariable+4];
-// particle.zold =particlesInfomation[i*numberOFVariable+5];
-// particle.vx =particlesInfomation[i*numberOFVariable+6];
-// particle.vy =particlesInfomation[i*numberOFVariable+7];
-// particle.vz =particlesInfomation[i*numberOFVariable+8];
-// particle.vxold =particlesInfomation[i*numberOFVariable+9];
-// particle.vyold =particlesInfomation[i*numberOFVariable+10];
-// particle.vzold =particlesInfomation[i*numberOFVariable+11];
-// particle.vxoldf =particlesInfomation[i*numberOFVariable+12];
-// particle.vyoldf =particlesInfomation[i*numberOFVariable+13];
-// particle.vzoldf =particlesInfomation[i*numberOFVariable+14];
-// particle.dx =particlesInfomation[i*numberOFVariable+15];
-// particle.level =particlesInfomation[i*numberOFVariable+16];
-// particle.particleNumber =particlesInfomation[i*numberOFVariable+17];
-//
-// int index=particle.particleNumber;
-// val<1>(x[index])=particle.x; val<2>(x[index])=particle.y; val<3>(x[index])=particle.z;
-// val<1>(xold[index])=particle.xold; val<2>(xold[index])=particle.yold; val<3>(xold[index])=particle.zold;
-// val<1>(vxold[index])=particle.vxold; val<2>(vxold[index])=particle.vyold; val<3>(vxold[index])=particle.vzold;
-// val<1>(vxoldf[index])=particle.vxoldf; val<2>(vxoldf[index])=particle.vyoldf; val<3>(vxoldf[index])=particle.vzoldf;
-// }
-//
-//}
-//
-//
-//
-
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::updateParticles()
-{
- //find particles which want to go to another process
- BOOST_FOREACH(ParticlesPtr particle, particles)
- {
- int minInitLevel = this->grid->getCoarsestInitializedLevel();
- int maxInitLevel = this->grid->getFinestInitializedLevel();
- LBMKernelETD3Q27Ptr kernel;
- DistributionArray3DPtr distributions;
- BCArray3D<D3Q27BoundaryCondition> bcArray;
- Block3DPtr block;
-
- for(int level = minInitLevel; level<=maxInitLevel; level++)
- {
- UbTupleInt3 blockIndexes = grid->getBlockIndexes(particle->x,particle->y,particle->z,level);
- block = grid->getBlock(val<1>(blockIndexes), val<2>(blockIndexes), val<3>(blockIndexes), level);
- if(!block) continue;
- if(block->isActive())
- {
- if(rank != block->getRank())
- {
- std::map<int,int>::iterator p;
- p=neighbors.find(block->getRank());
- p->second++;
- particle->rankOfParticle=block->getRank();
- particle->level=level;
- break;//break should be here
- }
- //break;//break should be here
- }
- }
- }
- //send number of particle moving to another process
- for ( std::map< int, int >::const_iterator iter = neighbors.begin();iter != neighbors.end(); ++iter )
- {
- int svalue=iter->second;
- if(comm->getNumberOfProcesses() > 1)
- {
- MPI_Send(&svalue,1, MPI_INT,iter->first,iter->first,MPI_COMM_WORLD);
- //std::cout<<"\n"<<svalue<<" send from rank "<<rank<<" to rank "<< iter->first;
- }
- }
- std::map<int,int> receiveNeighbors;
- //receive number of particle coming form another process
- for ( std::map< int, int >::const_iterator iter = neighbors.begin();iter != neighbors.end(); ++iter )
- {
- int rvalue;
- if(comm->getNumberOfProcesses() > 1)
- {
- MPI_Status status;
- MPI_Recv(&rvalue,1, MPI_INT,iter->first,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
- //std::cout<<"\n"<<"rank "<< rank<<" receive "<<rvalue<<" from rank "<<iter->first;
- receiveNeighbors.insert((std::make_pair(iter->first,rvalue)));
- }
- }
- //int numOfvarClass=18;
-
- ///send data
- for ( std::map< int, int >::const_iterator iter = neighbors.begin();iter != neighbors.end(); ++iter )
- {
- if (iter->second!=0)
- {
- std::vector<double> sendParticlesInfomation;
- int lengthSending=numOfvarClass*iter->second;
-
- BOOST_FOREACH(ParticlesPtr particle, particles)
- {
- if (particle->rankOfParticle==iter->first)
- {
- getVectorFromParticles(sendParticlesInfomation,particle);
- }
- }
- MPI_Send(&sendParticlesInfomation[0],lengthSending, MPI_DOUBLE,iter->first,iter->first,MPI_COMM_WORLD);
- std::cout<<"\n send from rank "<<rank<<" to rank "<< iter->first;
- }
- }
- //delete particle
- std::list<ParticlesPtr> newparticles;
- BOOST_FOREACH(ParticlesPtr particle, particles)
- {
- if (particle->rankOfParticle==rank)
- {
- newparticles.push_back(particle);
- }
- else
- {
- // int idNumber =particle->ID;
- // printParticle(idNumber);
- //delete data
- }
- }
- particles=newparticles;
- ///receive data
- for ( std::map< int, int >::const_iterator iter = receiveNeighbors.begin();iter != receiveNeighbors.end(); ++iter )
- {
- if (iter->second!=0)
- {
- int lengthReceiving=numOfvarClass*iter->second;
- std::vector<double> rvalue(lengthReceiving);
- MPI_Status status;
- MPI_Recv(&rvalue[0],lengthReceiving, MPI_DOUBLE,iter->first,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
- std::cout<<"\n"<<"rank "<< rank<<" receive from rank "<<iter->first;
- getParticlesFromVector(rvalue,numOfvarClass);
- }
- }
- //reset neighbors
- std::map<int,int> Neighbors2;
- Neighbors2=neighbors;
- neighbors.clear();
- for ( std::map< int, int >::const_iterator iter = Neighbors2.begin();iter != Neighbors2.end(); ++iter )
- {
- neighbors.insert((std::make_pair(iter->first,0)));
- }
-
-
-
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::checkParticles()
-{
- std::list<ParticlesPtr> newparticles;
- BOOST_FOREACH(ParticlesPtr particle, particles)
- {
- int minInitLevel = this->grid->getCoarsestInitializedLevel();
- int maxInitLevel = this->grid->getFinestInitializedLevel();
- Block3DPtr block;
-
- for(int level = minInitLevel; level<=maxInitLevel; level++)
- {
- UbTupleInt3 blockIndexes = grid->getBlockIndexes(particle->x,particle->y,particle->z,level);
- block = grid->getBlock(val<1>(blockIndexes), val<2>(blockIndexes), val<3>(blockIndexes), level);
- if(!block) continue;
- if(block->isActive())
- {
- LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- if(!kernel) continue;
- if(rank == block->getRank())
- {
- if(!checkNodes(block,particle->x,particle->y,particle->z))
- {
- std::cout<<"particle number "<<particle->ID <<"is gone";
- continue;
- }
- newparticles.push_back(particle);
- break;
- }
- }
- }
- }
- particles=newparticles;
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::finalMovement(ParticlesPtr particle,double tau[])
-{
- Block3DPtr block;
- UbTupleInt3 blockIndexes = grid->getBlockIndexes(particle->x, particle->y, particle->z,particle->level);
- block = grid->getBlock(val<1>(blockIndexes), val<2>(blockIndexes), val<3>(blockIndexes), particle->level);
- LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- DistributionArray3DPtr distributions = kernel->getDataSet()->getFdistributions();
- BCArray3D<D3Q27BoundaryCondition> bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
-
- UbTupleInt3 nodeIndexes = grid->getNodeIndexes(block, particle->x, particle->y, particle->z);
- int ix1 = val<1>(nodeIndexes);
- int ix2 = val<2>(nodeIndexes);
- int ix3 = val<3>(nodeIndexes);
-
- if(!iProcessor->iCellHasSolid(bcArray, ix1, ix2, ix3))
- {
- interpolMovement(block,distributions,particle,tau,ix1, ix2, ix3);
- }
- else
- {
- extrapolMovement( block, particle, tau, ix1, ix2, ix3);
- }
-
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::interpolMovement(Block3DPtr block,DistributionArray3DPtr& distributions,ParticlesPtr particle,double tau[],int ix1, int ix2, int ix3)
-{
- D3Q27ICell iCell;
- //LBMReal f[D3Q27System::ENDF+1];
- LBMReal dx =grid->getDeltaX(particle->level);
- LBMReal vx1,vx2,vx3;
- iProcessor->readICell(distributions, iCell, ix1, ix2, ix3);
- double x1LB, x2LB, x3LB;
- UbTupleDouble3 orgNodeRW = grid->getNodeCoordinates(block, ix1, ix2, ix3);
- double offsetX1RW = abs(particle->x - val<1>(orgNodeRW));
- double offsetX2RW = abs(particle->y - val<2>(orgNodeRW));
- double offsetX3RW = abs(particle->z - val<3>(orgNodeRW));
- x1LB = offsetX1RW / dx;
- x2LB = offsetX2RW / dx;
- x3LB = offsetX3RW / dx;
-
- x1LB -= 0.5;
- x2LB -= 0.5;
- x3LB -= 0.5;
- LBMReal omega = LBMSystem::calcCollisionFactor(nue, block->getLevel());
- LBMReal tauxx, tauyy, tauzz, tauxy, tauxz, tauyz;
- iHelper->interpolate8to1WithVelocityWithShearStress(iCell, x1LB, x2LB, x3LB, omega,vx1,vx2,vx3,tauxx,tauyy,tauzz,tauxy,tauxz,tauyz);
- vx1 = vx1 * conv->getFactorVelocityLbToW();
- vx2 = vx2 * conv->getFactorVelocityLbToW();
- vx3 = vx3 * conv->getFactorVelocityLbToW();
- tau[0]=tauxx/dx;
- tau[1]=tauyy/dx;
- tau[2]=tauzz/dx;
- tau[3]=tauxy/dx;
- tau[4]=tauxz/dx;
- tau[5]=tauyz/dx;
-
- if (particleHasMass) { CalcVelParticle(dx,vx1,vx2,vx3,particle->vxold,particle->vyold,particle->vzold,particle->vxoldf,particle->vyoldf,particle->vzoldf);}
- //if (particleHasMass) { CalcVelParticle2(dx,vx1,vx2,vx3,particle->vxold,particle->vyold,particle->vzold,particle->vxoldf,particle->vyoldf,particle->vzoldf);}
- //heuns method
- particle->x= (particle->xold) + 1.0/2.0*(vx1 + particle->vxold)*conv->getFactorTimeLbToW(dx);
- particle->y= (particle->yold) + 1.0/2.0*(vx2 + particle->vyold)*conv->getFactorTimeLbToW(dx);
- particle->z= (particle->zold) + 1.0/2.0*(vx3 + particle->vzold)*conv->getFactorTimeLbToW(dx);
-
- particle->vxold = vx1;
- particle->vyold = vx2;
- particle->vzold = vx3;
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::extrapolMovement(Block3DPtr block,ParticlesPtr particle,double tau[],int ix1, int ix2, int ix3)
-{
-
- LBMReal vx1,vx2,vx3;
- LBMReal dx =grid->getDeltaX(particle->level);
- int level =particle->level;
- LBMReal omega = LBMSystem::calcCollisionFactor(nue,level);
- D3Q27ICell iCell;
- //LBMReal f[D3Q27System::ENDF+1];
- double xp000,yp000,zp000;
- double Xw=999.0,Yw=999.0,Zw=999.0;
- double cofWeightx=999.0,cofWeighty=999.0,cofWeightz=999.0;
- double A,B,C,D,ii=0.0;
- double disx=9990,disy=9990,disz=9990,disxy=9990,disxz=9990,disyz=9990,disxyz=9990;
-
- UbTupleDouble3 orgNodeRW = grid->getNodeCoordinates(block, ix1, ix2, ix3);
- double ix1ph = ( val<1>(orgNodeRW));
- double ix2ph = ( val<2>(orgNodeRW));
- double ix3ph = ( val<3>(orgNodeRW));
- findPlane(ix1,ix2,ix3,grid,block,A,B,C,D,ii);
- double s = A*particle->x + B*particle->y + C*particle->z + D;//The sign of s = Ax + By + Cz + D determines which side the point (x,y,z) lies with respect to the plane. If s > 0 then the point lies on the same side as the normal (A,B,C). If s < 0 then it lies on the opposite side, if s = 0 then the point (x,y,z) lies on the plane.
- if (s>0){s=1;} else if (s<0){s=-1;}else {s=0;}
- double normalDis=((A*particle->x + B*particle->y + C*particle->z + D)/sqrt(A*A+B*B+C*C));///distance point to plane xp-Xw=distance
- double di=A/sqrt(A*A+B*B+C*C); double dj=B/sqrt(A*A+B*B+C*C); double dk=C/sqrt(A*A+B*B+C*C);
- double XXw=particle->x-di*normalDis; double YYw=particle->y-dj*normalDis; double ZZw=particle->z-dk*normalDis;
- findInitialCell(s,di,dj,dk,dx,ix1ph,ix2ph,ix3ph,xp000,yp000,zp000);//find initial cell with xp000,yp000,zp000
-
- double dis[7]={disx,disy,disz,disxy,disxz,disyz,disxyz};
- getAllDis(dis,particle->x,particle->y,particle->z,ix1ph,ix2ph,ix3ph,xp000,yp000,zp000,level);
- rearangedDouble(dis);
- getAfterCompare(dis,dx,ix1ph,ix2ph,ix3ph,xp000,yp000,zp000,Xw,Yw,Zw,particle->x,particle->y,particle->z,A,B,C, D, normalDis,di,dj,dk,cofWeightx,cofWeighty,cofWeightz,iCell,level);
-
- double x1LB, x2LB, x3LB;
- double offsetX1RW000 = (particle->x - xp000);
- double offsetX2RW000 = (particle->y - yp000);
- double offsetX3RW000 = (particle->z - zp000);
-
- x1LB = offsetX1RW000 / dx;
- x2LB = offsetX2RW000 / dx;
- x3LB = offsetX3RW000 / dx;
- x1LB -= 0.5;
- x2LB -= 0.5;
- x3LB -= 0.5;
- // outICell(iCell);
- //iProcessor->interpolate8to1WithVelocity(iCell, f, x1LB, x2LB, x3LB, omega,vx1,vx2,vx3);
- LBMReal tauxx, tauyy, tauzz, tauxy, tauxz, tauyz;
- iHelper->interpolate8to1WithVelocityWithShearStress(iCell, x1LB, x2LB, x3LB, omega,vx1,vx2,vx3,tauxx,tauyy,tauzz,tauxy,tauxz,tauyz);
- addCorrection(vx1,vx2,vx3,tauxx,tauyy,tauzz,tauxy,tauxz,tauyz,dx,iCell,Xw, Yw, Zw, omega,cofWeightx,cofWeighty,cofWeightz,ii, x1LB, x2LB, x3LB,di,dj,dk);
- tau[0]=tauxx;
- tau[1]=tauyy;
- tau[2]=tauzz;
- tau[3]=tauxy;
- tau[4]=tauxz;
- tau[5]=tauyz;
- // if (particleHasMass) { CalcVelParticle(dx,vx1,vx2,vx3,vx1old,vx2old,vx3old,vx1oldf,vx2oldf,vx3oldf);}
- if (particleHasMass) { CalcVelParticle2(dx,vx1,vx2,vx3,particle->vxold,particle->vyold,particle->vzold,particle->vxoldf,particle->vyoldf,particle->vzoldf);}
- //heuns method
- particle->x= (particle->xold) + 1.0/2.0*(vx1 + particle->vxold)*conv->getFactorTimeLbToW(dx);
- particle->y= (particle->yold) + 1.0/2.0*(vx2 + particle->vyold)*conv->getFactorTimeLbToW(dx);
- particle->z= (particle->zold) + 1.0/2.0*(vx3 + particle->vzold)*conv->getFactorTimeLbToW(dx);
- collWall(A,B,C,D,particle->x,particle->y,particle->z,particle->xold,particle->yold,particle->zold,dx,vx1,vx2,vx3,ii);
- particle->vxold = vx1;
- particle->vyold = vx2;
- particle->vzold = vx3;
-
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::CalcVelParticle(double dx,double &vx1,double &vx2,double &vx3,double vx1old,double vx2old,
- double vx3old,double &vx1oldf,double &vx2oldf,double &vx3oldf)
-{
- double signDragx,signDragy,signDragz;
- double dragCoff=1.0;
- if (vx1old>=vx1oldf) { signDragx=-1.0; }
- else if(vx1old<vx1oldf) { signDragx= 1.0; }
- if (vx2old>=vx2oldf) { signDragy=-1.0; }
- else if(vx2old< vx2oldf) { signDragy= 1.0; }
- if (vx3old>=vx3oldf) { signDragz=-1.0; }
- else if(vx3old< vx3oldf) { signDragz= 1.0; }
- double vx1p=vx1old+signDragx*(vx1old-vx1oldf)*(vx1old-vx1oldf)*conv->getFactorTimeLbToW(dx);
- double vx2p=vx2old+signDragy*(vx2old-vx2oldf)*(vx2old-vx2oldf)*conv->getFactorTimeLbToW(dx);
- double vx3p=vx3old+signDragz*(vx3old-vx3oldf)*(vx3old-vx3oldf)*conv->getFactorTimeLbToW(dx);
-
- double signDragNextx,signDragNexty,signDragNextz;
- if (vx1p>=vx1) { signDragNextx=-1.0; }
- else if(vx1p< vx1) { signDragNextx= 1.0; }
- if (vx2p>=vx2) { signDragNexty=-1.0; }
- else if(vx2p< vx2) { signDragNexty= 1.0; }
- if (vx3p>=vx3) { signDragNextz=-1.0; }
- else if(vx3p< vx3) { signDragNextz= 1.0; }
- ////////////////////velocity of particle//////////////////////////////////////////////////////
- double velx1oldf=vx1;
- double velx2oldf=vx2;
- double velx3oldf=vx3;
-
- vx1=vx1old+1.0/2.0*( signDragNextx*(vx1p-vx1)*(vx1p-vx1)+signDragx*(vx1old-vx1oldf)*(vx1old-vx1oldf) )*conv->getFactorTimeLbToW(dx);
- vx2=vx2old+1.0/2.0*( signDragNexty*(vx2p-vx2)*(vx2p-vx2)+signDragy*(vx2old-vx2oldf)*(vx2old-vx2oldf) )*conv->getFactorTimeLbToW(dx);
- vx3=vx3old+1.0/2.0*( signDragNextz*(vx3p-vx3)*(vx3p-vx3)+signDragz*(vx3old-vx3oldf)*(vx3old-vx3oldf) )*conv->getFactorTimeLbToW(dx);
-
- vx1oldf=velx1oldf;
- vx2oldf=velx2oldf;
- vx3oldf=velx3oldf;
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::CalcVelParticle2(double dx,double &vx1,double &vx2,double &vx3,double vx1old,double vx2old,
- double vx3old,double &vx1oldf,double &vx2oldf,double &vx3oldf)
-{
- LBMReal muRE = 1.002*1e-3;//m2/s
- double mass=2.3*1e-17;// kg
- double Diameter=600*1e-9;//m
- double deltatime=conv->getFactorTimeLbToW(dx)*0.001;
- double Coff=3*PI*Diameter*muRE*deltatime/mass;
- double exCoff=exp(-Coff);
- //velocity of particle///
- double velx1oldf=vx1;
- double velx2oldf=vx2;
- double velx3oldf=vx3;
- //air
- vx1=vx1old*exCoff+(velx1oldf-(velx1oldf)*exCoff);
- vx2=vx2old*exCoff+(velx2oldf-(velx2oldf)*exCoff);
- vx3=vx3old*exCoff+(velx3oldf-(velx3oldf)*exCoff);
-
- vx1oldf=velx1oldf;
- vx2oldf=velx2oldf;
- vx3oldf=velx3oldf;
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::findPlane(int ix1,int ix2,int ix3,Grid3DPtr grid,Block3DPtr block,double &A,double &B,double &C,double &D,double &ii)
-{
- double x1plane=0.0,y1plane=0.0,z1plane=0.0;
- double x2plane=0.0,y2plane=0.0,z2plane=0.0;
- double x3plane=0.0,y3plane=0.0,z3plane=0.0;
- D3Q27BoundaryConditionPtr bcPtr;
- double dx = grid->getDeltaX(block);
- LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- DistributionArray3DPtr distributions = kernel->getDataSet()->getFdistributions();
- BCArray3D<D3Q27BoundaryCondition> bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
-
- for (int k = ix3; k <= ix3 + 1; k++){
- for(int j = ix2; j <= ix2 + 1; j++){
- for(int i = ix1; i <= ix1 + 1; i++)
- {
- UbTupleDouble3 pointplane1 = grid->getNodeCoordinates(block, i, j, k);
-
- double iph=val<1>(pointplane1);
- double jph=val<2>(pointplane1);
- double kph=val<3>(pointplane1);
-
- if(!bcArray.isSolid(i, j, k))
- {
- bcPtr=bcArray.getBC(i,j,k);
- if(bcPtr)
- {
- for(int fdir=D3Q27System::FSTARTDIR; fdir<=D3Q27System::FENDDIR; fdir++)
- {
- if( ii<=2)
- {
- LBMReal q = bcPtr->getQ(fdir);
- if (q!=999.00000)
- {
-
- if ( fdir==D3Q27System::E )
- {
- if (i+q<=ix1+1)
- {
- if (ii==0) { x1plane=iph+q*dx; y1plane=jph; z1plane=kph; ii++; }
- else if (ii==1) { x2plane=iph+q*dx; y2plane=jph; z2plane=kph; if (x1plane!=x2plane||y1plane!=y2plane||z1plane!=z2plane) ii++; }
- else if(ii==2) { x3plane=iph+q*dx; y3plane=jph; z3plane=kph; if ((x3plane!=x1plane||y3plane!=y1plane||z3plane!=z1plane)&&(x2plane!=x3plane||y2plane!=y3plane||z2plane!=z3plane)) ii++;}
- }
- }
- if ( fdir==D3Q27System::W )
- {
- if (i-q>=ix1)
- {
- if (ii==0) { x1plane=iph-q*dx; y1plane=jph; z1plane=kph; ii++; }
- else if (ii==1) { x2plane=iph-q*dx; y2plane=jph; z2plane=kph; if (x1plane!=x2plane||y1plane!=y2plane||z1plane!=z2plane) ii++; }
- else if(ii==2) { x3plane=iph-q*dx; y3plane=jph; z3plane=kph; if ((x3plane!=x1plane||y3plane!=y1plane||z3plane!=z1plane)&&(x2plane!=x3plane||y2plane!=y3plane||z2plane!=z3plane)) ii++;}
- }
- }
- if ( fdir==D3Q27System::N )
- {
- if(j+q<=ix2+1)
- {
- if (ii==0) { x1plane=iph; y1plane=jph+q*dx; z1plane=kph; ii++; }
- else if (ii==1) { x2plane=iph; y2plane=jph+q*dx; z2plane=kph; if (x1plane!=x2plane||y1plane!=y2plane||z1plane!=z2plane) ii++; }
- else if (ii==2) { x3plane=iph; y3plane=jph+q*dx; z3plane=kph; if ((x3plane!=x1plane||y3plane!=y1plane||z3plane!=z1plane)&&(x2plane!=x3plane||y2plane!=y3plane||z2plane!=z3plane)) ii++;}
- }
- }
- if ( fdir==D3Q27System::S )
- {
- if (j-q>=ix2)
- {
- if (ii==0) { x1plane=iph; y1plane=jph-q*dx; z1plane=kph; ii++; }
- else if (ii==1) { x2plane=iph; y2plane=jph-q*dx; z2plane=kph; if (x1plane!=x2plane||y1plane!=y2plane||z1plane!=z2plane) ii++; }
- else if (ii==2) { x3plane=iph; y3plane=jph-q*dx; z3plane=kph; if ((x3plane!=x1plane||y3plane!=y1plane||z3plane!=z1plane)&&(x2plane!=x3plane||y2plane!=y3plane||z2plane!=z3plane)) ii++;}
- }
- }
-
- if ( fdir==D3Q27System::T )
- {
- if(k+q<=ix3+1)
- {
- if (ii==0) { x1plane=iph; y1plane=jph; z1plane=kph+q*dx; ii++; }
- else if (ii==1) { x2plane=iph; y2plane=jph; z2plane=kph+q*dx; if (x1plane!=x2plane||y1plane!=y2plane||z1plane!=z2plane) ii++; }
- else if (ii==2) { x3plane=iph; y3plane=jph; z3plane=kph+q*dx; if ((x3plane!=x1plane||y3plane!=y1plane||z3plane!=z1plane)&&(x2plane!=x3plane||y2plane!=y3plane||z2plane!=z3plane)) ii++;}
- }
- }
- if ( fdir==D3Q27System::B )
- {
- if (k-q>=ix3)
- {
- if (ii==0) { x1plane=iph; y1plane=jph; z1plane=kph-q*dx; ii++; }
- else if (ii==1) { x2plane=iph; y2plane=jph; z2plane=kph-q*dx; if (x1plane!=x2plane||y1plane!=y2plane||z1plane!=z2plane) ii++; }
- else if (ii==2) { x3plane=iph; y3plane=jph; z3plane=kph-q*dx; if ((x3plane!=x1plane||y3plane!=y1plane||z3plane!=z1plane)&&(x2plane!=x3plane||y2plane!=y3plane||z2plane!=z3plane)) ii++;}
- }
- }
-
- }
- }
- }
-
- }
-
- }
- }
- }
- }
-
- A = y1plane* (z2plane - z3plane) + y2plane*(z3plane - z1plane) + y3plane* (z1plane - z2plane);
- B = z1plane* (x2plane - x3plane) + z2plane*(x3plane - x1plane) + z3plane* (x1plane - x2plane) ;
- C = x1plane* (y2plane - y3plane) + x2plane*(y3plane - y1plane) + x3plane* (y1plane - y2plane) ;
- D =-( x1plane*(y2plane*z3plane - y3plane*z2plane)+x2plane*(y3plane*z1plane - y1plane*z3plane) + x3plane* (y1plane* z2plane - y2plane* z1plane));
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::findInitialCell(double s,double di,double dj,double dk,double dx,double ix1ph,double ix2ph,double ix3ph,double &xp000,double &yp000,double &zp000)
-{
- if (s*di>0) { xp000=ix1ph+dx; }
- else if(s*di<0) { xp000=ix1ph-dx; }
- else { xp000=ix1ph; }
- if (s*dj>0) { yp000=ix2ph+dx; }
- else if(s*dj<0) { yp000=ix2ph-dx; }
- else { yp000=ix2ph; }
- if (s*dk>0) { zp000=ix3ph+dx; }
- else if(s*dk<0) { zp000=ix3ph-dx; }
- else { zp000=ix3ph; }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::finddisPointToCentercell(double x1,double x2,double x3,double xp000, double yp000, double zp000,double &dis,int level)
-{
- UbTupleInt3 blockIndexes = grid->getBlockIndexes(xp000, yp000, zp000,level);
- Block3DPtr block = grid->getBlock(val<1>(blockIndexes), val<2>(blockIndexes), val<3>(blockIndexes), level);
- if(block&& block->isActive())
- {
- LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- if (kernel)
- {
- {
- DistributionArray3DPtr distributions = kernel->getDataSet()->getFdistributions();
- BCArray3D<D3Q27BoundaryCondition> bcArray = boost::dynamic_pointer_cast<D3Q27ETBCProcessor>(kernel->getBCProcessor())->getBCArray();
- if(block && block->isActive())
- {
- double dx = grid->getDeltaX(block);
- UbTupleInt3 nodeIndexes = grid->getNodeIndexes(block, xp000, yp000, zp000);
- int ix1Se = val<1>(nodeIndexes);
- int ix2Se = val<2>(nodeIndexes);
- int ix3Se = val<3>(nodeIndexes);
- if (!iProcessor->iCellHasSolid(bcArray, ix1Se, ix2Se, ix3Se ))
- {
- UbTupleDouble3 pointCell = grid->getNodeCoordinates(block, ix1Se, ix2Se, ix3Se); double iph=val<1>(pointCell)+dx/2.0; double jph=val<2>(pointCell)+dx/2.0; double kph=val<3>(pointCell)+dx/2.0;
- dis=sqrt((iph-x1)*(iph-x1)+(jph-x2)*(jph-x2)+(kph-x3)*(kph-x3));
- }
- }
- }
- }
-
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::collWall(double A,double B,double C,double D,double &x1,double &x2,double &x3,double x1old,double x2old,double x3old,double dx,double &vx1,double &vx2,double &vx3,double ii)
-{
- if (particleHasMass)
- {
- double s = A*x1 + B*x2 + C*x3 + D;//The sign of s = Ax + By + Cz + D determines which side the point (x,y,z) lies with respect to the plane. If s > 0 then the point lies on the same side as the normal (A,B,C). If s < 0 then it lies on the opposite side, if s = 0 then the point (x,y,z) lies on the plane.
- if (s>0){s=1;} else if (s<0){s=-1;}else {s=0;}
- double normalDis=((A*x1 + B*x2 + C*x3 + D)/sqrt(A*A+B*B+C*C));///distance point to plane xp-Xw=distance
- double di=A/sqrt(A*A+B*B+C*C); double dj=B/sqrt(A*A+B*B+C*C); double dk=C/sqrt(A*A+B*B+C*C);
- A*=s;B*=s;C*=s;
- if (abs(normalDis)<0.05*dx)
- {
- if (ii==3)
- {
- ///reflection without e factor
- /* double vn=(A*vx1 + B*vx2 + C*vx3)/(A*A+B*B+C*C);
- double vx1L=vx1-2*vn*A;
- double vx2L=vx2-2*vn*B;
- double vx3L=vx3-2*vn*C;*/
- double vnx=A*(A*vx1 + B*vx2 + C*vx3)/(A*A+B*B+C*C);
- double vny=B*(A*vx1 + B*vx2 + C*vx3)/(A*A+B*B+C*C);
- double vnz=C*(A*vx1 + B*vx2 + C*vx3)/(A*A+B*B+C*C);
- ////do collision
- double CollitionFactor_n=0.01;
- double CollitionFactor_bt=0.01;
- vnx=-1*CollitionFactor_n*vnx;
- vny=-1*CollitionFactor_n*vny;
- vnz=-1*CollitionFactor_n*vnz;
- double vbtx=(B*B*vx1 + C*C*vx1 - A*B*vx2 - A*C*vx3)/(A*A+B*B+C*C);
- double vbty=(-(A*B*vx1) + A*A*vx2 + C*C*vx2 - B*C*vx3)/(A*A+B*B+C*C);
- double vbtz=(-(A*C*vx1) - B*C*vx2 + A*A*vx3 + B*B*vx3)/(A*A+B*B+C*C);
- vbtx=vbtx*CollitionFactor_bt;
- vbty=vbty*CollitionFactor_bt;
- vbtz=vbtz*CollitionFactor_bt;
- vx1=vnx+vbtx;
- vx2=vny+vbty;
- vx3=vnz+vbtz;
-
- x1 = (x1old);
- x2 = (x2old);
- x3 = (x3old);
- }
- else std::cout<<"there is no plane to reflect";
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::addCorrection(double &vx1,double &vx2,double &vx3,double &tauxx,double &tauyy,double &tauzz,double &tauxy,double &tauxz,double &tauyz,
- double dx,D3Q27ICell iCell,double Xw, double Yw, double Zw,double omega,double cofWeightx,double cofWeighty,double cofWeightz,double ii,double x1LB,double x2LB,double x3LB,double di,double dj,double dk)
-{
- //LBMReal f[D3Q27System::ENDF+1];
- double Udif,Vdif,Wdif;
- iHelper->interpolate8to1WithVelocity(iCell, Xw, Yw, Zw, omega,Udif,Vdif,Wdif);
-
- double Ucx,Ucy,Ucz;
- double Vcx,Vcy,Vcz;
- double Wcx,Wcy,Wcz;
-
- double dxUcx,dyUcy,dzUcz;
- double dxVcx,dyVcy,dzVcz;
- double dxWcx,dyWcy,dzWcz;
- //double di,dj,dk;
-
- if (cofWeightx==0){di=0;}
- if (cofWeighty==0){dj=0;}
- if (cofWeightz==0){dk=0;}
-
- double weightx=(di*di)/(di*di+dj*dj+dk*dk);
- double weighty=(dj*dj)/(di*di+dj*dj+dk*dk);
- double weightz=(dk*dk)/(di*di+dj*dj+dk*dk);
- // /* if (Xw<0.5||Xw>-0.5){weightx=0;}
- // if (Yw<0.5||Yw>-0.5){weighty=0;}
- // if (Zw<0.5||Zw>-0.5){weightz=0;}
- //*/
- if (weightx==0){Ucx=0;Vcx=0;Wcx=0;dxUcx=0;dxVcx=0;dxWcx=0;}else{
- Ucx=-Udif*weightx/(Xw*Xw*Xw-Xw/4.0)*(x1LB*x1LB*x1LB-x1LB/4.0);Vcx=-Vdif*weightx/(Xw*Xw*Xw-Xw/4.0)*(x1LB*x1LB*x1LB-x1LB/4.0);Wcx=-Wdif*weightx/(Xw*Xw*Xw-Xw/4.0)*(x1LB*x1LB*x1LB-x1LB/4.0);
- dxUcx=-Udif*weightx/(Xw*Xw*Xw-Xw/4.0)*(3.0*x1LB*x1LB-1.0/4.0);dxVcx=-Vdif*weightx/(Xw*Xw*Xw-Xw/4.0)*(3.0*x1LB*x1LB-1.0/4.0);dxWcx=-Wdif*weightx/(Xw*Xw*Xw-Xw/4.0)*(3.0*x1LB*x1LB-1.0/4.0);
- }
- if (weighty==0){Ucy=0;Vcy=0;Wcy=0;dyUcy=0;dyVcy=0;dyWcy=0;}else{
- Ucy=-Udif*weighty/(Yw*Yw*Yw-Yw/4.0)*(x2LB*x2LB*x2LB-x2LB/4.0);Vcy=-Vdif*weighty/(Yw*Yw*Yw-Yw/4.0)*(x2LB*x2LB*x2LB-x2LB/4.0);Wcy=-Wdif*weighty/(Yw*Yw*Yw-Yw/4.0)*(x2LB*x2LB*x2LB-x2LB/4.0);
- dyUcy=-Udif*weighty/(Yw*Yw*Yw-Yw/4.0)*(3.0*x2LB*x2LB-1.0/4.0);dyVcy=-Vdif*weighty/(Yw*Yw*Yw-Yw/4.0)*(3.0*x2LB*x2LB-1.0/4.0);dyWcy=-Wdif*weighty/(Yw*Yw*Yw-Yw/4.0)*(3.0*x2LB*x2LB-1.0/4.0);
- }
- if (weightz==0){Ucz=0;Vcz=0;Wcz=0;dzUcz=0;dzVcz=0;dzWcz=0;}else{
- Ucz=-Udif*weightz/(Zw*Zw*Zw-Zw/4.0)*(x3LB*x3LB*x3LB-x3LB/4.0);Vcz=-Vdif*weightz/(Zw*Zw*Zw-Zw/4.0)*(x3LB*x3LB*x3LB-x3LB/4.0);Wcz=-Wdif*weightz/(Zw*Zw*Zw-Zw/4.0)*(x3LB*x3LB*x3LB-x3LB/4.0);
- dzUcz=-Udif*weightz/(Zw*Zw*Zw-Zw/4.0)*(3.0*x3LB*x3LB-1.0/4.0);dzVcz=-Vdif*weightz/(Zw*Zw*Zw-Zw/4.0)*(3.0*x3LB*x3LB-1.0/4.0);dzWcz=-Wdif*weightz/(Zw*Zw*Zw-Zw/4.0)*(3.0*x3LB*x3LB-1.0/4.0);
- }
-
- double Ucor=Ucx+Ucy+Ucz;
- double Vcor=Vcx+Vcy+Vcz;
- double Wcor=Wcx+Wcy+Wcz;
-
- double tauxxC=dxUcx;
- double tauyyC=dyVcy;
- double tauzzC=dzWcz;
- double tauxyC=0.5*(dyUcy+dxVcx);
- double tauxzC=0.5*(dzUcz+dxWcx);
- double tauyzC=0.5*(dzVcz+dyWcy);
-
- if (ii!=3)
- {
- //std::cout<<"there are not 3point for making plane"<<endl<<"Ucor="<<Ucor<<endl<<"vx1="<<vx1;
- Ucor =0.0;
- Vcor =0.0;
- Wcor =0.0;
- tauxxC=0.0;
- tauyyC=0.0;
- tauzzC=0.0;
- tauxyC=0.0;
- tauxzC=0.0;
- tauyzC=0.0;
- }
-
- vx1+=Ucor;
- vx2+=Vcor;
- vx3+=Wcor;
- tauxx+=tauxxC;
- tauyy+=tauyyC;
- tauzz+=tauzzC;
- tauxy+=tauxyC;
- tauxz+=tauxzC;
- tauyz+=tauyzC;
-
- vx1 = vx1 * conv->getFactorVelocityLbToW();
- vx2 = vx2 * conv->getFactorVelocityLbToW();
- vx3 = vx3 * conv->getFactorVelocityLbToW();
- tauxx=tauxx/dx;
- tauyy=tauyy/dx;
- tauzz=tauzz/dx;
- tauxy=tauxy/dx;
- tauxz=tauxz/dx;
- tauyz=tauyz/dx;
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::getAllDis(double dis[],double x,double y,double z,double ix1ph,double ix2ph,double ix3ph,double xp000,double yp000,double zp000,int level)
-{
- finddisPointToCentercell(x,y,z,xp000, ix2ph, ix3ph,dis[0],level);//disx
- finddisPointToCentercell(x,y,z,ix1ph, yp000, ix3ph,dis[1],level);//disy
- finddisPointToCentercell(x,y,z,ix1ph, ix2ph, zp000,dis[2],level);//disz
- finddisPointToCentercell(x,y,z,xp000, yp000, ix3ph,dis[3],level);//disxy
- finddisPointToCentercell(x,y,z,xp000, ix2ph, zp000,dis[4],level);//disxz
- finddisPointToCentercell(x,y,z,ix1ph, yp000, zp000,dis[5],level);//disyz
- finddisPointToCentercell(x,y,z,xp000, yp000, zp000,dis[6],level);//disxyz
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::getAfterCompare(double dis[],double dx,double ix1ph,double ix2ph,double ix3ph,double &xp000,double &yp000,double &zp000,double &Xw,
- double &Yw,double &Zw,double x1,double x2,double x3,double A,double B,double C,double D,double normalDis, double di,double dj,double dk, double &cofWeightx
- ,double &cofWeighty,double &cofWeightz,D3Q27ICell &iCell,int level)
-{
- for (int i=0;i<7;i++)
- {
- double dist=dis[i];
- if (dist==dirx ){ Xw=(B*x2+C*x3+D)/(-A); Yw=x2; Zw=x3; cofWeightx=1; cofWeighty=0; cofWeightz=0; xp000=xp000; yp000=ix2ph; zp000=ix3ph; getIcell(iCell,xp000, ix2ph, ix3ph,level); } //get suitable block from suitable cell
- else if (dist==diry ){ Xw=x1; Yw=(A*x1+C*x3+D)/(-B); Zw=x3; cofWeightx=0; cofWeighty=1; cofWeightz=0; xp000=ix1ph; yp000=yp000; zp000=ix3ph; getIcell(iCell,ix1ph, yp000, ix3ph,level); }
- else if (dist==dirz ){ Xw=x1; Yw=x3; Zw=(B*x2+A*x1+D)/(-C);cofWeightx=0; cofWeighty=0; cofWeightz=1; xp000=ix1ph; yp000=ix2ph; zp000=zp000; getIcell(iCell,ix1ph, ix2ph, zp000,level); }
- else if (dist==dirxy ){ Xw=x1-di*normalDis; Yw=x2-dj*normalDis; Zw=x3; cofWeightx=1; cofWeighty=1; cofWeightz=0; xp000=xp000; yp000=yp000; zp000=ix3ph; getIcell(iCell,xp000, yp000, ix3ph,level); }
- else if (dist==dirxz ){ Xw=x1-di*normalDis; Yw=x2; Zw=x3-dk*normalDis; cofWeightx=1; cofWeighty=0; cofWeightz=1; xp000=xp000; yp000=ix2ph; zp000=zp000; getIcell(iCell,xp000, ix2ph, zp000,level); }
- else if (dist==diryz ){ Xw=x1; Yw=x2-dj*normalDis; Zw=x3-dk*normalDis; cofWeightx=0; cofWeighty=1; cofWeightz=1; xp000=ix1ph; yp000=yp000; zp000=zp000; getIcell(iCell,ix1ph, yp000, zp000,level); }
- else if (dist==dirxyz){ Xw=x1-di*normalDis; Yw=x2-dj*normalDis; Zw=x3-dk*normalDis; cofWeightx=1; cofWeighty=1; cofWeightz=1; xp000=xp000; yp000=yp000; zp000=zp000; getIcell(iCell,xp000, yp000, zp000,level); }
-
- Xw=(Xw-xp000)/dx-0.5;
- Yw=(Yw-yp000)/dx-0.5;
- Zw=(Zw-zp000)/dx-0.5;
- if (Xw<2.5&&Xw>-2.5&&Yw<2.50&&Yw>-2.50&&Zw<2.50&&Zw>-2.50){break;}
- }
-}
-
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::rearangedDouble( double dis[7])
-{
- double item;
- for (int i=6;i>0;i--){
- for(int j=0;j<i;j++){
- if (dis[j]>dis[j+1])
- {
- item=dis[j];
- dis[j]=dis[j+1];
- dis[j+1]=item;
- }
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::getIcell(D3Q27ICell &iCell,double xp000, double yp000, double zp000,int level)
-{
- UbTupleInt3 blockIndexes = grid->getBlockIndexes(xp000, yp000, zp000,level);
- Block3DPtr block= grid->getBlock(val<1>(blockIndexes), val<2>(blockIndexes), val<3>(blockIndexes), level);
- LBMKernelETD3Q27Ptr kernel = boost::dynamic_pointer_cast<LBMKernelETD3Q27>(block->getKernel());
- DistributionArray3DPtr distributions = kernel->getDataSet()->getFdistributions();
- iProcessor->readICell(distributions, iCell, (int)xp000, (int)yp000, (int)zp000);
-}
-////////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::printParticle(ParticlesPtr particle)
-
-{
- /*BOOST_FOREACH(ParticlesPtr particle, particles)
- {*/
- int index=particle->ID;
- (*files[index]).precision (std::numeric_limits<double>::digits10 + 1);
- (*files[index])<<index<<" "<<istep<<" "<<particle->x<<" "<< particle->y<<" "<<particle->z<<" "<<particle->vxold<<" "<<particle->vyold<<" "<<particle->vzold<< std::endl;
- //}
-}
-
-///////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::printParticle(int index)
-
-{
-
-}
-//////////////////////////////////////////////////////////////////////////
-void PathLineCoProcessorMcpart::gatherData(ParticlesPtr particle)
-{
- //int index=particle->ID;
- //std::vector<std::vector<double> > tempData;
- //tempData.resize(7);
- //if (!data.count(index ))
- //{
- // data.insert((std::make_pair(index,tempData)));
- //}
- //tempData=data.at(index);
- //int t=0;
- //tempData[t++].push_back(istep);
- //tempData[t++].push_back(particle->x);
- //tempData[t++].push_back(particle->y);
- //tempData[t++].push_back(particle->z);
- //tempData[t++].push_back(particle->vxold);
- //tempData[t++].push_back(particle->vyold);
- //tempData[t++].push_back(particle->vzold);
-
- // data.at(index)=tempData;
-
- datanames.resize(0);
- datanames.push_back("ID");
- datanames.push_back("TimeStep");
- datanames.push_back("Vx");
- datanames.push_back("Vy");
- datanames.push_back("Vz");
- //datanames.push_back("tauxx");
- //datanames.push_back("tauyy");
- //datanames.push_back("tauzz");
- //datanames.push_back("tauxy");
- //datanames.push_back("tauxz");
- //datanames.push_back("tauyz");
- //datanames.push_back("xoff");
- //datanames.push_back("yoff");
- //datanames.push_back("zoff");
-
- data.resize(datanames.size());
-
- int index = 0;
- data[index++].push_back(particle->ID);
- data[index++].push_back(istep);
- data[index++].push_back(particle->vxold);
- data[index++].push_back(particle->vyold);
- data[index++].push_back(particle->vzold);
- //data[index++].push_back(tauxx);
- //data[index++].push_back(tauyy);
- //data[index++].push_back(tauzz);
- //data[index++].push_back(tauxy);
- //data[index++].push_back(tauxz);
- //data[index++].push_back(tauyz);
- //data[index++].push_back(xoff);
- //data[index++].push_back(yoff);
- //data[index++].push_back(zoff);
-
- nodes.push_back( makeUbTuple( double(particle->x), double(particle->y), double(particle->z)) );
-
-}
diff --git a/depot/PathLineCoProcessorMcpart.h b/depot/PathLineCoProcessorMcpart.h
deleted file mode 100644
index 0661640ee..000000000
--- a/depot/PathLineCoProcessorMcpart.h
+++ /dev/null
@@ -1,133 +0,0 @@
-/*
-* PathLinePostprocessor.h
-*
-* Created on: 3.6.2012
-* Author: Fard
-*/
-
-
-#ifndef PATHLINEPOSTPROCESSORMCPART_H
-#define PATHLINEPOSTPROCESSORMCPART_H
-
-# define dirx 0
-# define diry 1
-# define dirz 2
-# define dirxy 3
-# define dirxz 4
-# define diryz 5
-# define dirxyz 6
-#include <mpi.h>
-#include "CoProcessor.h"
-#include "Grid3D.h"
-#include "Block3D.h"
-#include "LBMUnitConverter.h"
-#include "Communicator.h"
-#include "D3Q27InterpolationProcessor.h"
-#include "LBMKernelETD3Q27.h"
-#include "D3Q27ETBCProcessor.h"
-#include <boost/shared_ptr.hpp>
-//#include <boost/tr1/memory.hpp>
-#include "Particles.h"
-#include "basics/writer/WbWriterVtkXmlASCII.h"
-
-class PathLineCoProcessorMcpart;
-typedef boost::shared_ptr<PathLineCoProcessorMcpart> PathLineCoProcessorMcpartPtr;
-
-class PathLineCoProcessorMcpart : public CoProcessor
-{
-public:
- //typedef std::map<int,ParticlesPtr> ParticlesMap;
- std::map<int,int> neighbors;
- //std::map<ofstream,std::vector<int> > pathOfstream;
- PathLineCoProcessorMcpart(Grid3DPtr grid, const std::string& path, WbWriter* const writer,
- LBMUnitConverterPtr conv, UbSchedulerPtr s, CommunicatorPtr comm,
- std::vector<UbTupleDouble3 > Positions, LBMReal nue, D3Q27InterpolationProcessorPtr iProcessor);
- ~PathLineCoProcessorMcpart();
- void process(double step);
- void getNeighborsRank();
-
-protected:
- void collectPostprocessData();
- //void getParticles(std::vector<Particles> *particles,std::vector<UbTupleDouble3> x);
- void initializeParticle(double x, double y,double z,int i,int level);
- void getParticlePerProccess(std::vector<UbTupleDouble3 >);
- void getParticle(Particles &particle, double dx,double x, double y,double z,int level,int i);
- void initialMovement();
- void updateDistributedValues(std::vector<Particles> *particles,std::vector<UbTupleDouble3> &x);
- void receiveStatusOfPoint(bool &status,int rankRoot,double x1,double x2,double x3,double level);
- void sendStatusOfPoint(bool status,double x1,double x2,double x3,double level,int i);
- void receiveBool(bool &variable,int rankBlock);
- void receiveInt(int &rvalue,int root);
- void sendBool(bool _bool,int distinationRank,int i );
- void sendInt(int svalue,int distinationRank ,int tag);
- void getVectorFromParticles(std::vector<double> &particlesInfomation,ParticlesPtr particles);
- void getParticlesFromVector(std::vector<double> particlesInfomation,int numberOFVariable);
- void updateinfo(std::vector<Particles> *particlesVec,std::vector<UbTupleDouble3> &x);
- void allGatherDoubles(std::vector<double>& svalues, std::vector<double>& rvalues);
- void fillOutPositions(std::vector<double> particlesInfomation,std::vector<UbTupleDouble3> &x,double numberOFVariable);
- void updateParticles();
- void checkParticles();
- void finalMovement(ParticlesPtr particle,double tau[]);
- void interpolMovement(Block3DPtr block,DistributionArray3DPtr& distributions,ParticlesPtr particle,double tau[],int ix1, int ix2, int ix3);
- void extrapolMovement(Block3DPtr block, ParticlesPtr particle,double tau[],int ix1, int ix2, int ix3);
- void findPlane(int ix1,int ix2,int ix3,Grid3DPtr grid,Block3DPtr block,double &A,double &B,double &C,double &D,double &ii);
- void CalcVelParticle2(double dx,double &vx1,double &vx2,double &vx3,double vx1old,double vx2old,double vx3old,double &vx1oldf,double &vx2oldf,double &vx3oldf);
- void CalcVelParticle(double dx,double &vx1,double &vx2,double &vx3,double vx1old,double vx2old,double vx3old,double &vx1oldf,double &vx2oldf,double &vx3oldf);
- void findInitialCell(double s,double di,double dj,double dk,double dx,double ix1ph,double ix2ph,double ix3ph,double &xp000,double &yp000,double &zp000);
- void finddisPointToCentercell(double x1,double x2,double x3,double xp000, double yp000, double zp000,double &dis,int level);
- void collWall(double A,double B,double C,double D,double &x1,double &x2,double &x3,double x1old,double x2old,double x3old,double dx,double &vx1,double &vx2,double &vx3,double ii);
- void addCorrection(double &vx1,double &vx2,double &vx3,double &tauxx,double &tauyy,double &tauzz,double &tauxy,double &tauxz,double &tauyz,
- double dx,D3Q27ICell iCell,double Xw, double Yw, double Zw,double omega,double cofWeightx,double cofWeighty,double cofWeightz,double ii,double x1LB,double x2LB,double x3LB,double di,double dj,double dk);
- void getAllDis(double dis[],double x,double y,double z,double ix1ph,double ix2ph,double ix3ph,double xp000,double yp000,double zp000,int level);
- void rearangedDouble( double dis[7]);
- void getWallPosition(double dx,double xp000,double yp000,double zp000,double &Xw,double &Yw,double &Zw,double x1,double x2,double x3,double A
- ,double B,double C,double D,double normalDis, double di,double dj,double dk, double &cofWeightx,double &cofWeighty,double &cofWeightz,int dir);
- void getIcell(D3Q27ICell &iCell,double xp000, double yp000, double zp000,int level);
- void getAfterCompare(double dis[],double dx,double ix1ph,double ix2ph,double ix3ph,double &xp000,double &yp000,double &zp000,double &Xw,
- double &Yw,double &Zw,double x1,double x2,double x3,double A,double B,double C,double D,double normalDis, double di,double dj,double dk, double &cofWeightx
- ,double &cofWeighty,double &cofWeightz,D3Q27ICell &iCell,int level);
- void printParticle(ParticlesPtr particle);
- void printParticle(int index);
- void initializeForPrinting(int number);
- void gatherData(ParticlesPtr particle);
-private:
- bool checkNodes(Block3DPtr block,double _x1,double _x2,double _x3);
- CommunicatorPtr comm;
- //std::map<int ,std::vector<UbTupleDouble3> > nodes;
- //std::map<int ,std::vector<std::string> >datanames;
- //std::map<int ,std::vector<std::vector<double> > >data;
-
- std::vector<UbTupleDouble3> nodes;
- std::vector<std::string> datanames;
- std::vector<std::vector<double> > data;
-
- std::vector<boost::shared_ptr<std::ofstream> > files;
- std::string path;
- WbWriter* writer;
- LBMUnitConverterPtr conv;
-
- D3Q27InterpolationProcessorPtr interpolation;
- //std::vector<UbTupleDouble3> x;
- //std::vector<UbTupleDouble3> xold;
- //std::vector<UbTupleDouble3 >vxold;
- //std::vector<UbTupleDouble3 > vxoldf;
- bool particleHasMass;
- int isExtrapolation;
- int istep;
- int stepcheck;
- LBMReal nue;
- D3Q27InterpolationProcessorPtr iProcessor;
- D3Q27InterpolationHelperPtr iHelper;
- //Particles particles;
- int rank;
- int root;
- int maxtag;
- std::list<ParticlesPtr> particles;
- //std::vector<std::vector<UbTupleDouble3> > Positions;
- //Particles particles;
-
- //ParticlesMap particles;
-
-};
-#endif
-
--
GitLab