From e88f391fb0d4da4843b0b8f98b71c382e6812be4 Mon Sep 17 00:00:00 2001
From: Konstantin Kutscher <kutscher@irmb.tu-bs.de>
Date: Mon, 1 Oct 2018 16:43:52 +0200
Subject: [PATCH] add ShadowCopy, add WritePeBlocksCoProcessor, add nozzles,
stable simulation, InteractorsHelper add delete flag
---
source/Applications/Thermoplast/config.txt | 17 +-
.../Applications/Thermoplast/thermoplast.cpp | 216 +++++++++---
source/DemCoupling/DemCoProcessor.cpp | 310 ++++++++++++++----
source/DemCoupling/DemCoProcessor.h | 7 +-
.../DemCoupling/MovableObjectInteractor.cpp | 65 +++-
.../DemCoupling/WritePeBlocksCoProcessor.cpp | 88 +++++
source/DemCoupling/WritePeBlocksCoProcessor.h | 43 +++
.../pe/PeLoadBalancerAdapter.cpp | 11 +-
.../pe/PePhysicsEngineGeometryAdapter.cpp | 42 +--
.../pe/PePhysicsEngineGeometryAdapter.h | 5 +
.../pe/PePhysicsEngineSolverAdapter.cpp | 35 +-
.../pe/PePhysicsEngineSolverAdapter.h | 5 +-
.../CoProcessors/WriteBlocksCoProcessor.cpp | 2 +-
.../CoProcessors/WriteBlocksCoProcessor.h | 6 +-
.../Interactors/InteractorsHelper.cpp | 17 +-
.../Interactors/InteractorsHelper.h | 3 +-
16 files changed, 679 insertions(+), 193 deletions(-)
create mode 100644 source/DemCoupling/WritePeBlocksCoProcessor.cpp
create mode 100644 source/DemCoupling/WritePeBlocksCoProcessor.h
diff --git a/source/Applications/Thermoplast/config.txt b/source/Applications/Thermoplast/config.txt
index c681fec35..654263c4d 100644
--- a/source/Applications/Thermoplast/config.txt
+++ b/source/Applications/Thermoplast/config.txt
@@ -1,20 +1,29 @@
#simulation parameters
#boundingBox = 0 0 0 300 1520 2320
-boundingBox = 60 1370 10 190 1530 320 #test bb
+boundingBox = 60 1370 130 190 1530 320 #test bb
+
+#boundingBox = 60 0 10 190 1530 750 #test bb 2
#boundingBox = 60 0 10 190 1530 2320 #production bb
blocknx = 10 10 10
#blocknx = 300 420 320
endTime = 100000
-outTime = 500
+outTime = 1
availMem = 25e9
uLB = 0.1
#PE parameters
+#test pe offset
peMinOffset = 46 2 2
+#peMaxOffset = -30 -60 -12
peMaxOffset = -8 -25 -2
+
+#production pe offset
+#peMinOffset = 46 18 14
+#peMaxOffset = -8 -25 -23
+
sphereTime = 10
#geometry files
@@ -22,11 +31,11 @@ pathGeo = d:/Projects/ThermoPlast/SimGeo
michel = /michel.stl
plexiglas = /plexiglas.stl
-pathOut = g:/temp/thermoplastCluster6
+pathOut = g:/temp/thermoplastCluster1
#restart
cpStart = 20000000
cpStep = 20000000
restart = false
-nupsStep = 1000 1000 10000000
\ No newline at end of file
+nupsTime = 10 10 1000000
\ No newline at end of file
diff --git a/source/Applications/Thermoplast/thermoplast.cpp b/source/Applications/Thermoplast/thermoplast.cpp
index d8f1e10cc..484795506 100644
--- a/source/Applications/Thermoplast/thermoplast.cpp
+++ b/source/Applications/Thermoplast/thermoplast.cpp
@@ -30,6 +30,7 @@
#include <DummyPhysicsEngineMaterialAdapter.h>
#include <DummyPhysicsEngineGeometryAdapter.h>
#include <WriteDemObjectsCoProcessor.h>
+#include <WritePeBlocksCoProcessor.h>
#include "CreateDemObjectsCoProcessor.h"
#include "RestartDemObjectsCoProcessor.h"
@@ -51,6 +52,61 @@ vector<double> peMaxOffset;
string pathOut;// = "d:/temp/thermoplastCluster";
string pathGeo;// = "d:/Projects/ThermoPlast/Geometrie";
+void addNozzle(SPtr<Grid3D> grid, SPtr<Communicator> comm, SPtr<BCAdapter> noSlipBCAdapter, InteractorsHelper& intHelper)
+{
+ int myid = comm->getProcessID();
+ if (myid==0) UBLOG(logINFO, "Add nozzles:start");
+
+ std::vector< SPtr<Interactor3D> > interactors;
+
+ for (int i = 0; i <= 387; i++)
+ {
+ SPtr<GbTriFaceMesh3D> bbGeo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/Nozzle/bb"+UbSystem::toString(i)+".stl", "bb", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
+ SPtr<Interactor3D> bbInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(bbGeo, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES));
+ //intHelper.addInteractor(bbInt);
+ interactors.push_back(bbInt);
+ }
+
+ for (int i = 0; i <= 51; i++)
+ {
+ SPtr<GbTriFaceMesh3D> bsGeo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/Nozzle/bs"+UbSystem::toString(i)+".stl", "bs", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
+ SPtr<Interactor3D> bsInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(bsGeo, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES));
+ //intHelper.addInteractor(bsInt);
+ interactors.push_back(bsInt);
+ }
+
+ std::array<int,5> n = {0,1,2,3,6};
+
+ for (int i = 0; i < n.size(); i++)
+ {
+ SPtr<GbTriFaceMesh3D> biGeo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/Nozzle/bi"+UbSystem::toString(n[i])+".stl", "bs", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
+ SPtr<Interactor3D> biInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(biGeo, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES));
+ //intHelper.addInteractor(biInt);
+ interactors.push_back(biInt);
+ }
+
+
+ for (SPtr<Interactor3D> interactor : interactors)
+ {
+ std::vector< std::shared_ptr<Block3D> > blockVector;
+ UbTupleInt3 blockNX=grid->getBlockNX();
+ SPtr<GbObject3D> geoObject(interactor->getGbObject3D());
+ double ext = 0.0;
+ std::array<double, 6> AABB ={ geoObject->getX1Minimum(),geoObject->getX2Minimum(),geoObject->getX3Minimum(),geoObject->getX1Maximum(),geoObject->getX2Maximum(),geoObject->getX3Maximum() };
+ grid->getBlocksByCuboid(AABB[0]-(double)val<1>(blockNX)*ext, AABB[1]-(double)val<2>(blockNX)*ext, AABB[2]-(double)val<3>(blockNX)*ext, AABB[3]+(double)val<1>(blockNX)*ext, AABB[4]+(double)val<2>(blockNX)*ext, AABB[5]+(double)val<3>(blockNX)*ext, blockVector);
+ for (std::shared_ptr<Block3D> block : blockVector)
+ {
+ if (block->getKernel())
+ {
+ interactor->setBCBlock(block);
+ }
+ }
+ interactor->initInteractor();
+ }
+
+ if (myid==0) UBLOG(logINFO, "Add nozzles:end");
+}
+
std::shared_ptr<DemCoProcessor> makePeCoProcessor(SPtr<Grid3D> grid, SPtr<Communicator> comm, const SPtr<UbScheduler> peScheduler, const std::shared_ptr<LBMUnitConverter> lbmUnitConverter, int maxpeIterations)
{
double peRelaxtion = 0.7;
@@ -58,17 +114,18 @@ std::shared_ptr<DemCoProcessor> makePeCoProcessor(SPtr<Grid3D> grid, SPtr<Commun
//Beschleunigung g
double g = 9.81 * lbmUnitConverter->getFactorAccWToLb();
//Vector3D globalLinearAcc(0.0, -g, 0.0);
- Vector3D globalLinearAcc(0.0, 0.0, -g);
+ //Vector3D globalLinearAcc(0.0, 0.0, -g);
+ Vector3D globalLinearAcc(0.0, 0.0, 0.0);
std::shared_ptr<PePhysicsEngineMaterialAdapter> planeMaterial = std::make_shared<PePhysicsEngineMaterialAdapter>("granular", 1.0, 0, 0.1 / 2, 0.1 / 2, 0.5, 1, 1, 0, 0);
- const int gridNx = val<1>(grid->getBlockNX()) * grid->getNX1();
- const int gridNy = val<2>(grid->getBlockNX()) * grid->getNX2();
- const int gridNz = val<3>(grid->getBlockNX()) * grid->getNX3();
+ const int gridNX1 = val<1>(grid->getBlockNX()) * grid->getNX1();
+ const int gridNX2 = val<2>(grid->getBlockNX()) * grid->getNX2();
+ const int gridNX3 = val<3>(grid->getBlockNX()) * grid->getNX3();
//UbTupleInt3 simulationDomain(gridNx, gridNy, gridNz);
//std::array<double, 6> simulationDomain = {1, 1, 1, 30, 30, 30};
- std::array<double, 6> simulationDomain ={ g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3 };
+ std::array<double, 6> simulationDomain ={ g_minX1, g_minX2, g_minX3, g_minX1+gridNX1, g_minX2+gridNX2, g_minX3+gridNX3 };
UbTupleInt3 numberOfBlocks(grid->getNX1(), grid->getNX2(), grid->getNX3());
//UbTupleInt3 numberOfBlocks((simulationDomain[3]-simulationDomain[0])/val<1>(grid->getBlockNX()), (simulationDomain[4]-simulationDomain[1])/val<2>(grid->getBlockNX()), (simulationDomain[5]-simulationDomain[2])/val<3>(grid->getBlockNX()));
UbTupleBool3 isPeriodic(grid->isPeriodicX1(), grid->isPeriodicX2(), grid->isPeriodicX3());
@@ -83,12 +140,30 @@ std::shared_ptr<DemCoProcessor> makePeCoProcessor(SPtr<Grid3D> grid, SPtr<Commun
std::shared_ptr<PeLoadBalancerAdapter> loadBalancer(new PeLoadBalancerAdapter(grid, comm->getNumberOfProcesses(), comm->getProcessID()));
std::shared_ptr<PhysicsEngineSolverAdapter> peSolver = std::make_shared<PePhysicsEngineSolverAdapter>(peParamter, loadBalancer);
+ SPtr<CoProcessor> peblocks(new WritePeBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm, std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(peSolver)->getBlockForest()));
+ peblocks->process(0);
+ peblocks.reset();
+
const std::shared_ptr<ForceCalculator> forceCalculator = std::make_shared<ForceCalculator>(comm);
return std::make_shared<DemCoProcessor>(grid, peScheduler, comm, forceCalculator, peSolver);
}
-//pipe flow with forcing
+void createSpheres(double radius, Vector3D origin, double uLB, SPtr<CreateDemObjectsCoProcessor> createSphereCoProcessor)
+{
+ double d = 2.0*radius;
+ int maxX2 = 5;
+ int maxX3 = 6;
+ for (int x3 = 0; x3 < maxX3; x3++)
+ for (int x2 = 0; x2 < maxX2; x2++)
+ for (int x1 = 0; x1 < 1; x1++)
+ {
+ SPtr<GbObject3D> sphere(new GbSphere3D(origin[0]+2.0*d*x1, origin[1]+x2*1.5*d, origin[2]+x3*1.5*d, radius));
+
+ createSphereCoProcessor->addGeoObject(sphere, Vector3D(uLB, -uLB+uLB/2.0*x2, -uLB+uLB/2.5*x3));
+ }
+}
+
void thermoplast(string configname)
{
SPtr<Communicator> comm = MPICommunicator::getInstance();
@@ -119,6 +194,8 @@ void thermoplast(string configname)
pathOut = config.getValue<string>("pathOut");
pathGeo = config.getValue<string>("pathGeo");
+ vector<int> nupsTime = config.getVector<int>("nupsTime");
+
//parameters
//string pathOut = "d:/temp/thermoplast3";
//string pathGeo = "d:/Projects/ThermoPlast/Geometrie";
@@ -131,7 +208,7 @@ void thermoplast(string configname)
double rhoLB = 0.0;
//double uLB = 0.1;
double radius = 5;
- double Re = 300;
+ double Re = 900;
double nuLB = (uLB*2.0*radius)/Re;
//geometry definition
@@ -200,6 +277,17 @@ void thermoplast(string configname)
grid->accept(genBlocks);
+ /////////////////////////////////////////////////////
+ ////PE domain test
+ //std::array<double, 6> simulationDomain ={ g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3 };
+ //Vector3D minOffset(peMinOffset[0], peMinOffset[1], peMinOffset[2]);
+ //Vector3D maxOffset(peMaxOffset[0], peMaxOffset[1], peMaxOffset[2]);
+ //SPtr<GbObject3D> boxPE(new GbCuboid3D(simulationDomain[0]+minOffset[0], simulationDomain[1]+minOffset[1], simulationDomain[2]+minOffset[2], simulationDomain[3]+maxOffset[0], simulationDomain[4]+maxOffset[1], simulationDomain[5]+maxOffset[2]));
+ //GbSystem3D::writeGeoObject(boxPE.get(), pathOut + "/geo/boxPE", WbWriterVtkXmlBinary::getInstance());
+ //return;
+ //////////////////////////////////////////////////////
+
+
if (myid == 0)
{
UBLOG(logINFO, "Parameters:");
@@ -215,8 +303,18 @@ void thermoplast(string configname)
UBLOG(logINFO, "Preprocess - start");
}
- GbCuboid3DPtr geoInflow1(new GbCuboid3D(g_minX1-blockLength, g_maxX2-120.0, g_minX3+190.0, g_minX1+1, g_maxX2+20.0, g_minX3+320.0));
- if (myid == 0) GbSystem3D::writeGeoObject(geoInflow1.get(), pathOut + "/geo/geoInflow1", WbWriterVtkXmlASCII::getInstance());
+ //GbCuboid3DPtr geoInflow1(new GbCuboid3D(g_minX1-blockLength, g_maxX2-120.0, g_minX3+190.0, g_minX1+1, g_maxX2+20.0, g_minX3+130.0));
+ GbCuboid3DPtr geoInjector5(new GbCuboid3D(-12, 1415, 205, 63, 1525, 315));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoInjector5.get(), pathOut + "/geo/geoInjector5", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr geoInjector4(new GbCuboid3D(-12, -5, 205, 63, 105, 315));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoInjector4.get(), pathOut + "/geo/geoInjector4", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr geoInjector7(new GbCuboid3D(28, 705, 542, 103, 815, 652));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoInjector7.get(), pathOut + "/geo/geoInjector7", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr testWallGeo(new GbCuboid3D(g_minX1-blockLength, g_minX2 - blockLength, g_maxX3, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(testWallGeo.get(), pathOut + "/geo/testWallGeo", WbWriterVtkXmlASCII::getInstance());
if (!restart)
{
@@ -236,19 +334,19 @@ void thermoplast(string configname)
if (myid==0) UBLOG(logINFO, "Read plexiglasGeo:end");
if (myid==0) GbSystem3D::writeGeoObject(plexiglasGeo.get(), pathOut+"/geo/plexiglasGeo", WbWriterVtkXmlBinary::getInstance());
- //Duese
- if (myid==0) UBLOG(logINFO, "Read Duese:start");
- SPtr<GbTriFaceMesh3D> s1Geo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/Duese/s1.stl", "plexiglasGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- SPtr<GbTriFaceMesh3D> b1Geo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/Duese/b1.stl", "plexiglasGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- SPtr<GbTriFaceMesh3D> p1Geo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/Duese/p1.stl", "plexiglasGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- SPtr<GbTriFaceMesh3D> p2Geo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/Duese/p2.stl", "plexiglasGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- if (myid==0) UBLOG(logINFO, "Read Duese:end");
+ ////Duese
+ //if (myid==0) UBLOG(logINFO, "Read Duese:start");
+ //SPtr<GbTriFaceMesh3D> s1Geo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/Duese/s1.stl", "plexiglasGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
+ //SPtr<GbTriFaceMesh3D> b1Geo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/Duese/b1.stl", "plexiglasGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
+ //SPtr<GbTriFaceMesh3D> p1Geo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/Duese/p1.stl", "plexiglasGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
+ //SPtr<GbTriFaceMesh3D> p2Geo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/Duese/p2.stl", "plexiglasGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
+ //if (myid==0) UBLOG(logINFO, "Read Duese:end");
//inflow
- GbCuboid3DPtr geoInflow2(new GbCuboid3D(g_minX1-blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_minX1, g_maxX2 + blockLength, g_maxX3 + blockLength));
- if (myid == 0) GbSystem3D::writeGeoObject(geoInflow2.get(), pathOut + "/geo/geoInflow2", WbWriterVtkXmlASCII::getInstance());
+ GbCuboid3DPtr geoOutflowMichel(new GbCuboid3D(g_minX1-blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_minX1, g_maxX2 + blockLength, g_maxX3 + blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoOutflowMichel.get(), pathOut + "/geo/geoOutflowMichel", WbWriterVtkXmlASCII::getInstance());
//GbCuboid3DPtr geoInflow3(new GbCuboid3D(g_minX1-blockLength, g_minX2-radius, g_maxX3-4.0*radius-1, g_minX1+1, g_maxX2+radius, g_maxX3-4.0*radius));
//if (myid == 0) GbSystem3D::writeGeoObject(geoInflow3.get(), pathOut + "/geo/geoInflow3", WbWriterVtkXmlASCII::getInstance());
@@ -257,21 +355,24 @@ void thermoplast(string configname)
//if (myid == 0) GbSystem3D::writeGeoObject(geoInflow4.get(), pathOut + "/geo/geoInflow4", WbWriterVtkXmlASCII::getInstance());
//outflow
- GbCuboid3DPtr geoOutflow(new GbCuboid3D(g_maxX1, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
- if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), pathOut + "/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance());
+ GbCuboid3DPtr geoOutflowPlexiglas(new GbCuboid3D(g_maxX1, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoOutflowPlexiglas.get(), pathOut + "/geo/geoOutflowPlexiglas", WbWriterVtkXmlASCII::getInstance());
//set boundary conditions for blocks and create process decomposition for MPI
SPtr<D3Q27Interactor> boxInt(new D3Q27Interactor(box, grid, noSlipBCAdapter, Interactor3D::INVERSESOLID));
//inflow
- SPtr<D3Q27Interactor> inflowInt1 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow1, grid, inflowAdapter, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> inflowInt2 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow2, grid, outflowAdapter, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> inflowInjector5Int = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInjector5, grid, inflowAdapter, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> inflowInjector4Int = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInjector4, grid, inflowAdapter, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> inflowInjector7Int = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInjector7, grid, inflowAdapter, Interactor3D::SOLID));
+
+ SPtr<D3Q27Interactor> outflowMichelInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoOutflowMichel, grid, outflowAdapter, Interactor3D::SOLID));
//SPtr<D3Q27Interactor> inflowInt3 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow3, grid, noSlipBCAdapter, Interactor3D::SOLID));
//SPtr<D3Q27Interactor> inflowInt4 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow4, grid, noSlipBCAdapter, Interactor3D::SOLID));
//outflow
- SPtr<D3Q27Interactor> outflowInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoOutflow, grid, outflowAdapter, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> outflowPlexiglasInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoOutflowPlexiglas, grid, outflowAdapter, Interactor3D::SOLID));
//michel
SPtr<Interactor3D> michelInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(michelGeo, grid, noSlipBCAdapter, Interactor3D::SOLID));
@@ -279,26 +380,32 @@ void thermoplast(string configname)
//plexiglas
SPtr<Interactor3D> plexiglasInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(plexiglasGeo, grid, noSlipBCAdapter, Interactor3D::SOLID));
- //Duese
- SPtr<Interactor3D> s1Int = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(s1Geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
- SPtr<Interactor3D> b1Int = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(b1Geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
- SPtr<Interactor3D> p1Int = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(p1Geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
- SPtr<Interactor3D> p2Int = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(p2Geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ ////Duese
+ //SPtr<Interactor3D> s1Int = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(s1Geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ //SPtr<Interactor3D> b1Int = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(b1Geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ //SPtr<Interactor3D> p1Int = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(p1Geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ //SPtr<Interactor3D> p2Int = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(p2Geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
+ SPtr<D3Q27Interactor> testWallInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(testWallGeo, grid, inflowAdapter, Interactor3D::SOLID));
//////////////////////////////////////////////////////////////////////////
//SPtr<Grid3DVisitor> peVisitor(new PePartitioningGridVisitor(comm, demCoProcessor));
SPtr<Grid3DVisitor> peVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::KWAY));
- InteractorsHelper intHelper(grid, peVisitor);
+ InteractorsHelper intHelper(grid, peVisitor, true);
intHelper.addInteractor(boxInt);
intHelper.addInteractor(michelInt);
intHelper.addInteractor(plexiglasInt);
- intHelper.addInteractor(s1Int);
- intHelper.addInteractor(b1Int);
- intHelper.addInteractor(p1Int);
- intHelper.addInteractor(p2Int);
- intHelper.addInteractor(inflowInt1);
- intHelper.addInteractor(outflowInt);
- intHelper.addInteractor(inflowInt2);
+ //addNozzle(grid,comm,noSlipBCAdapter,intHelper);
+ //////intHelper.addInteractor(s1Int);
+ //////intHelper.addInteractor(b1Int);
+ //////intHelper.addInteractor(p1Int);
+ //////intHelper.addInteractor(p2Int);
+ intHelper.addInteractor(inflowInjector5Int);
+ intHelper.addInteractor(inflowInjector4Int);
+ intHelper.addInteractor(inflowInjector7Int);
+ intHelper.addInteractor(outflowPlexiglasInt);
+ intHelper.addInteractor(outflowMichelInt);
+ intHelper.addInteractor(testWallInt);
intHelper.selectBlocks();
//write data for visualization of block grid
@@ -335,6 +442,8 @@ void thermoplast(string configname)
SetKernelBlockVisitor kernelVisitor(kernel, nuLB, availMem, needMem);
grid->accept(kernelVisitor);
+ addNozzle(grid,comm,noSlipBCAdapter,intHelper);
+
intHelper.setBC();
@@ -358,7 +467,7 @@ void thermoplast(string configname)
//PE initialization
double refLengthLb = radius*2.0;
- double refLengthWorld = refLengthLb * deltax;
+ double refLengthWorld = 2e-3;
const std::shared_ptr<LBMUnitConverter> lbmUnitConverter = std::make_shared<LBMUnitConverter>(refLengthWorld, LBMUnitConverter::WORLD_MATERIAL::AIR_20C, refLengthLb);
if (myid == 0) std::cout << lbmUnitConverter->toString() << std::endl;
double rhoSphere = 915 * lbmUnitConverter->getFactorDensityWToLb(); // kg/m^3
@@ -408,16 +517,27 @@ void thermoplast(string configname)
//sphere prototypes
//UBLOG(logINFO, "sphere prototypes - start, rank="<<myid);
double d = 2.0*radius;
- Vector3D origin(g_minX1+peMinOffset[0]+radius, geoInflow1->getX2Minimum()+4.0*d, geoInflow1->getX3Minimum()+2.0*d);
- for (int x3 = 0; x3 < 2; x3++)
- for (int x2 = 0; x2 < 2; x2++)
- for (int x1 = 0; x1 < 1; x1++)
- {
- //SPtr<GbObject3D> sphere(new GbSphere3D(origin[0]+x1*d, origin[1]+x2*2.0*d, origin[2]+x3*2.0*d, radius));
- SPtr<GbObject3D> sphere(new GbSphere3D(origin[0]+2.0*d, origin[1]+x2*1.5*d, origin[2]+x3*1.5*d, radius));
- if (myid == 0) GbSystem3D::writeGeoObject(sphere.get(), pathOut + "/geo/sphere"+UbSystem::toString(x1)+UbSystem::toString(x2)+UbSystem::toString(x3), WbWriterVtkXmlASCII::getInstance());
- createSphereCoProcessor->addGeoObject(sphere, Vector3D(uLB, 0.0, 0.0));
- }
+ Vector3D origin1(g_minX1+peMinOffset[0]+radius, geoInjector5->getX2Minimum()+2.0*d, geoInjector5->getX3Minimum()+2.0*d);
+ createSpheres(radius,origin1,uLB,createSphereCoProcessor);
+
+ //Vector3D origin2(g_minX1+peMinOffset[0]+radius, geoInjector4->getX2Minimum()+3.0*d, geoInjector4->getX3Minimum()+2.0*d);
+ //createSpheres(radius, origin2, uLB, createSphereCoProcessor);
+
+ //Vector3D origin3(g_minX1+peMinOffset[0]+radius, geoInjector7->getX2Minimum()+2.0*d, geoInjector7->getX3Minimum()+2.0*d);
+ //createSpheres(radius, origin3, uLB, createSphereCoProcessor);
+
+ //for (int x3 = 0; x3 < 6; x3++)
+ // for (int x2 = 0; x2 < 5; x2++)
+ // for (int x1 = 0; x1 < 1; x1++)
+ // {
+ // //SPtr<GbObject3D> sphere(new GbSphere3D(origin[0]+x1*d, origin[1]+x2*2.0*d, origin[2]+x3*2.0*d, radius));
+ // SPtr<GbObject3D> sphere(new GbSphere3D(origin[0]+2.0*d, origin[1]+x2*1.5*d, origin[2]+x3*1.5*d, radius));
+ // if (myid == 0) GbSystem3D::writeGeoObject(sphere.get(), pathOut + "/geo/sphere"+UbSystem::toString(x1)+UbSystem::toString(x2)+UbSystem::toString(x3), WbWriterVtkXmlASCII::getInstance());
+ // createSphereCoProcessor->addGeoObject(sphere, Vector3D(uLB, uLB, uLB));
+ // }
+
+
+
//UBLOG(logINFO, "sphere prototypes - stop, rank="<<myid);
//Vector3D origin(106+radius, 1372+radius, 12+radius);
@@ -447,7 +567,7 @@ void thermoplast(string configname)
writeDemObjectsCoProcessor->process(0);
////performance control
- SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
+ SPtr<UbScheduler> nupsSch(new UbScheduler(nupsTime[0], nupsTime[1], nupsTime[2]));
SPtr<NUPSCounterCoProcessor> npr(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
//start simulation
diff --git a/source/DemCoupling/DemCoProcessor.cpp b/source/DemCoupling/DemCoProcessor.cpp
index b779534df..dcf59afe5 100644
--- a/source/DemCoupling/DemCoProcessor.cpp
+++ b/source/DemCoupling/DemCoProcessor.cpp
@@ -10,6 +10,7 @@
#include "DistributionArray3D.h"
#include "BCProcessor.h"
#include "DataSet3D.h"
+#include "SetBcBlocksBlockVisitor.h"
#include "PhysicsEngineMaterialAdapter.h"
#include "PhysicsEngineGeometryAdapter.h"
@@ -43,44 +44,60 @@ DemCoProcessor::DemCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, SPtr<Comm
for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
{
walberla::pe::Storage* storage = blockIt->getData< walberla::pe::Storage >(*storageId.get());
- walberla::pe::BodyStorage* bodyStorage = &(*storage)[0];
+ bodyStorage = &(*storage)[0];
+ bodyStorageShadowCopies = &(*storage)[1];
bodyStorage->registerAddCallback("DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::addPeGeo), this));
- bodyStorage->registerRemoveCallback("DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::removePeGeo), this));
+ //bodyStorage->registerRemoveCallback("DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::removePeGeo), this));
+
+ bodyStorageShadowCopies->registerAddCallback("DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::addPeShadowGeo), this));
+ //bodyStorageShadowCopies->registerRemoveCallback("DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::removePeShadowGeo), this));
}
}
DemCoProcessor::~DemCoProcessor()
{
+ bodyStorage->deregisterAddCallback("DemCoProcessor");
+ //bodyStorage->deregisterRemoveCallback("DemCoProcessor");
+ if (&bodyStorage != &bodyStorageShadowCopies)
+ {
+ bodyStorageShadowCopies->deregisterAddCallback("DemCoProcessor");
+ //bodyStorageShadowCopies->deregisterRemoveCallback("DemCoProcessor");
+ }
}
void DemCoProcessor::addInteractor(std::shared_ptr<MovableObjectInteractor> interactor, std::shared_ptr<PhysicsEngineMaterialAdapter> physicsEngineMaterial, Vector3D initalVelocity)
{
interactors.push_back(interactor);
const int id = static_cast<int>(interactors.size()) - 1;
+ //if (comm->getProcessID()==2) UBLOG(logINFO, "DemCoProcessor::addInteractor() id = "<<id);
interactor->setID(id);
const auto peGeometryAdapter = this->createPhysicsEngineGeometryAdapter(interactor, physicsEngineMaterial);
if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(peGeometryAdapter)->isActive())
{
peGeometryAdapter->setLinearVelolocity(initalVelocity);
-
- std::vector< std::shared_ptr<Block3D> > blockVector;
- UbTupleInt3 blockNX=grid->getBlockNX();
- SPtr<GbObject3D> geoObject(interactor->getGbObject3D());
- double ext = 0.0;
- std::array<double, 6> AABB ={ geoObject->getX1Minimum(),geoObject->getX2Minimum(),geoObject->getX3Minimum(),geoObject->getX1Maximum(),geoObject->getX2Maximum(),geoObject->getX3Maximum() };
- grid->getBlocksByCuboid(AABB[0]-(double)val<1>(blockNX)*ext, AABB[1]-(double)val<2>(blockNX)*ext, AABB[2]-(double)val<3>(blockNX)*ext, AABB[3]+(double)val<1>(blockNX)*ext, AABB[4]+(double)val<2>(blockNX)*ext, AABB[5]+(double)val<3>(blockNX)*ext, blockVector);
- for (std::shared_ptr<Block3D> block : blockVector)
- {
- if (block->getKernel())
- {
- interactor->setBCBlock(block);
- //UBLOG(logINFO, "DemCoProcessor::addInteractor() rank = "<<comm->getProcessID());
- }
- }
- interactor->initInteractor();
}
+ SetBcBlocksBlockVisitor setBcVisitor(interactor);
+ grid->accept(setBcVisitor);
+
+ //std::vector< std::shared_ptr<Block3D> > blockVector;
+ //UbTupleInt3 blockNX=grid->getBlockNX();
+ //SPtr<GbObject3D> geoObject(interactor->getGbObject3D());
+ //double ext = 0.0;
+ //std::array<double, 6> AABB ={ geoObject->getX1Minimum(),geoObject->getX2Minimum(),geoObject->getX3Minimum(),geoObject->getX1Maximum(),geoObject->getX2Maximum(),geoObject->getX3Maximum() };
+ //grid->getBlocksByCuboid(AABB[0]-(double)val<1>(blockNX)*ext, AABB[1]-(double)val<2>(blockNX)*ext, AABB[2]-(double)val<3>(blockNX)*ext, AABB[3]+(double)val<1>(blockNX)*ext, AABB[4]+(double)val<2>(blockNX)*ext, AABB[5]+(double)val<3>(blockNX)*ext, blockVector);
+ //for (std::shared_ptr<Block3D> block : blockVector)
+ //{
+ // if (block->getKernel())
+ // {
+ // interactor->setBCBlock(block);
+ // //UBLOG(logINFO, "DemCoProcessor::addInteractor() rank = "<<comm->getProcessID());
+ // }
+ //}
+
+ interactor->initInteractor();
+
physicsEngineGeometrieAdapters.push_back(peGeometryAdapter);
}
@@ -94,8 +111,8 @@ std::shared_ptr<PhysicsEngineGeometryAdapter> DemCoProcessor::createPhysicsEngin
auto peGeometryAdapter = this->physicsEngineSolver->createPhysicsEngineGeometryAdapter(id, position, vfSphere->getRadius(), physicsEngineMaterial);
//if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(peGeometry)->isActive())
//{
- interactor->setPhysicsEngineGeometry(peGeometryAdapter);
- return peGeometryAdapter;
+ interactor->setPhysicsEngineGeometry(peGeometryAdapter);
+ return peGeometryAdapter;
//}
//else
//{
@@ -132,24 +149,24 @@ void DemCoProcessor::process(double actualTimeStep)
if (this->intermediateDemSteps == 1)
this->calculateDemTimeStep(demTimeStepsPerIteration);
-
-//#ifdef TIMING
-// if (comm->isRoot()) UBLOG(logINFO, "DemCoProcessor::calculateDemTimeStep() time = "<<timer.stop()<<" s");
-//#endif
- //if ((int)actualTimeStep % 100 == 0)
- //{
- // if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[0])->isActive())
- // {
- // //UBLOG(logINFO, "v: (x,y,z) " << physicsEngineGeometries[0]->getLinearVelocity() << " actualTimeStep = " << UbSystem::toString(actualTimeStep));
- // }
- //}
-
- // during the intermediate time steps of the collision response, the currently acting forces
- // (interaction forces, gravitational force, ...) have to remain constant.
- // Since they are reset after the call to collision response, they have to be stored explicitly before.
- // Then they are set again after each intermediate step.
- this->moveVfGeoObject();
+ //#ifdef TIMING
+ // if (comm->isRoot()) UBLOG(logINFO, "DemCoProcessor::calculateDemTimeStep() time = "<<timer.stop()<<" s");
+ //#endif
+ //if ((int)actualTimeStep % 100 == 0)
+ //{
+ // if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[0])->isActive())
+ // {
+ // //UBLOG(logINFO, "v: (x,y,z) " << physicsEngineGeometries[0]->getLinearVelocity() << " actualTimeStep = " << UbSystem::toString(actualTimeStep));
+ // }
+ //}
+
+ // during the intermediate time steps of the collision response, the currently acting forces
+ // (interaction forces, gravitational force, ...) have to remain constant.
+ // Since they are reset after the call to collision response, they have to be stored explicitly before.
+ // Then they are set again after each intermediate step.
+
+ this->moveVfGeoObjects();
#ifdef TIMING
if (comm->isRoot()) UBLOG(logINFO, "DemCoProcessor::moveVfGeoObject() time = "<<timer.stop()<<" s");
@@ -265,13 +282,23 @@ void DemCoProcessor::calculateDemTimeStep(double step)
//#endif
}
-void DemCoProcessor::moveVfGeoObject()
+void DemCoProcessor::moveVfGeoObjects()
{
for (int i = 0; i < interactors.size(); i++)
{
if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive())
{
- interactors[i]->moveGbObjectTo(physicsEngineGeometrieAdapters[i]->getPosition());
+ walberla::pe::RigidBody* peGeoObject = getPeGeoObject(std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->getId());
+ if (peGeoObject)
+ {
+ std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->setGeometry(peGeoObject);
+ interactors[i]->moveGbObjectTo(physicsEngineGeometrieAdapters[i]->getPosition());
+ }
+ else
+ {
+ std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->setInactive();
+ }
+
//UBLOG(logINFO, "DemCoProcessor::moveVfGeoObject() id = "<<std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->getId()<<" position="<<physicsEngineGeometrieAdapters[i]->getPosition()<<" rank="<<comm->getProcessID());
}
}
@@ -303,14 +330,14 @@ bool DemCoProcessor::isDemObjectInAABB(std::array<double, 6> AABB)
}
}
}
-
+
std::vector<int> values;
values.push_back((int)result);
std::vector<int> rvalues = comm->gather(values);
if (comm->isRoot())
{
- for (int i = 0; i < (int)rvalues.size(); i ++)
+ for (int i = 0; i < (int)rvalues.size(); i++)
{
result = result || (bool)rvalues[i];
}
@@ -354,52 +381,157 @@ void DemCoProcessor::getObjectsPropertiesVector(std::vector<double>& p)
void DemCoProcessor::addPeGeo(walberla::pe::RigidBody * peGeo)
{
- //UBLOG(logINFO, "DemCoProcessor::addPeGeo()");
- //for (int i = 0; i < physicsEngineGeometries.size(); i++)
+
+
+
+ if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
+
+ auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
+ geometry->counter++;
+ //if (comm->getProcessID()==3)UBLOG(logINFO, "DemCoProcessor::addPeGeo() geoID = "<<peGeo->getID()<<" counter="<<geometry->counter<<" shadowCounter="<<geometry->shadowCounter);
+ if (geometry->getId() == peGeo->getID())
+ {
+ geometry->setActive();
+ geometry->setGeometry(peGeo);
+ return;
+ }
+ else
+ throw UbException(UB_EXARGS, "PeGeo ID is not matching!");
+}
+
+void DemCoProcessor::removePeGeo(walberla::pe::RigidBody * peGeo)
+{
+ //if (comm->getProcessID()==2) UBLOG(logINFO, "DemCoProcessor::removePeGeo() rank = "<<comm->getProcessID());
+
+
+
+ if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
+
+ auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
+ geometry->counter--;
+
+ if (geometry->getId() == peGeo->getID())
+ {
+ if (geometry->shadowCounter == 0)
+ {
+ geometry->setInactive();
+ geometry->setGeometry(NULL);
+ }
+ //if (comm->getProcessID()==3)UBLOG(logINFO, "DemCoProcessor::removePeGeo() geoID = "<<peGeo->getID()<<" counter="<<geometry->counter<<" shadowCounter="<<geometry->shadowCounter);
+ return;
+ }
+ else
+ throw UbException(UB_EXARGS, "PeGeo ID is not matching!");
+
+
+ //if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
+
+ //std::shared_ptr<walberla::blockforest::BlockForest> forest = std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getBlockForest();
+ //std::shared_ptr<walberla::domain_decomposition::BlockDataID> storageId = std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getStorageId();
+
+ //for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
//{
- // auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i]);
- // if (geometry->getId() == peGeo->getID())
+ // for (auto bodyIt = walberla::pe::BodyIterator::begin(*blockIt, *storageId.get()); bodyIt != walberla::pe::BodyIterator::end(); ++bodyIt)
// {
- // geometry->setActive();
- // geometry->setGeometry(peGeo);
- // return;
+ // if (bodyIt->getID() == peGeo->getID())
+ // {
+ // return;
+ // }
// }
//}
+ //auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
+ //if (geometry->getId() == peGeo->getID())
+ //{
+ // geometry->setInactive();
+ // geometry->setGeometry(NULL);
+ // return;
+ //}
+ //else
+ // throw UbException(UB_EXARGS, "PeGeo ID is not matching!");
+}
+
+void DemCoProcessor::addPeShadowGeo(walberla::pe::RigidBody * peGeo)
+{
if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
+ geometry->shadowCounter++;
+ //if (comm->getProcessID()==3) UBLOG(logINFO, "DemCoProcessor::addPeShadowGeo() geoID = "<<peGeo->getID()<<" counter="<<geometry->counter<<" shadowCounter="<<geometry->shadowCounter);
if (geometry->getId() == peGeo->getID())
{
- geometry->setActive();
- geometry->setGeometry(peGeo);
+ if (geometry->shadowCounter == 1)
+ {
+ geometry->setActive();
+ geometry->setGeometry(peGeo);
+ }
+
return;
}
else
throw UbException(UB_EXARGS, "PeGeo ID is not matching!");
}
-void DemCoProcessor::removePeGeo(walberla::pe::RigidBody * peGeo)
+void DemCoProcessor::removePeShadowGeo(walberla::pe::RigidBody * peGeo)
{
- //UBLOG(logINFO, "DemCoProcessor::removePeGeo()");
- //for (int i = 0; i < physicsEngineGeometries.size(); i++)
+ //if (comm->getProcessID()==2) UBLOG(logINFO, "DemCoProcessor::removePeShadowGeo() rank = "<<comm->getProcessID());
+
+ //if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
+
+ //std::shared_ptr<walberla::blockforest::BlockForest> forest = std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getBlockForest();
+ //std::shared_ptr<walberla::domain_decomposition::BlockDataID> storageId = std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getStorageId();
+
+ //for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
+ //{
+ // for (auto bodyIt = walberla::pe::BodyIterator::begin(*blockIt, *storageId.get()); bodyIt != walberla::pe::BodyIterator::end(); ++bodyIt)
+ // {
+ // if (bodyIt->getID() == peGeo->getID())
+ // {
+ // return;
+ // }
+ // }
+ //}
+
+
+
+ //std::vector< std::shared_ptr<Block3D> > blockVector;
+ //UbTupleInt3 blockNX=grid->getBlockNX();
+ //SPtr<GbObject3D> geoObject(interactors[peGeo->getID()]->getGbObject3D());
+ //double ext = 0.0;
+ //std::array<double, 6> AABB ={ geoObject->getX1Minimum(),geoObject->getX2Minimum(),geoObject->getX3Minimum(),geoObject->getX1Maximum(),geoObject->getX2Maximum(),geoObject->getX3Maximum() };
+ //grid->getBlocksByCuboid(AABB[0]-(double)val<1>(blockNX)*ext, AABB[1]-(double)val<2>(blockNX)*ext, AABB[2]-(double)val<3>(blockNX)*ext, AABB[3]+(double)val<1>(blockNX)*ext, AABB[4]+(double)val<2>(blockNX)*ext, AABB[5]+(double)val<3>(blockNX)*ext, blockVector);
+ //if (blockVector.size() == 0)
//{
- // auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i]);
+ // auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
// if (geometry->getId() == peGeo->getID())
// {
// geometry->setInactive();
// geometry->setGeometry(NULL);
// return;
// }
+ // else
+ // throw UbException(UB_EXARGS, "PeGeo ID is not matching!");
//}
+
if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
+ geometry->shadowCounter--;
+
if (geometry->getId() == peGeo->getID())
{
- geometry->setInactive();
- geometry->setGeometry(NULL);
+ if (geometry->shadowCounter == 0)
+ {
+ geometry->setInactive();
+ geometry->setGeometry(NULL);
+ }
+ else
+ {
+ geometry->setActive();
+ geometry->setGeometry(peGeo);
+ }
+ //if (comm->getProcessID()==3)UBLOG(logINFO, "DemCoProcessor::removePeShadowGeo() geoID = "<<peGeo->getID()<<" counter="<<geometry->counter<<" shadowCounter="<<geometry->shadowCounter);
return;
}
else
@@ -414,7 +546,7 @@ bool DemCoProcessor::isSpheresIntersection(double centerX1, double centerX2, dou
if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive())
{
SPtr<GbObject3D> sphere = interactors[i]->getGbObject3D();
- result = result || ( sqrt(pow(sphere->getX1Centroid()-centerX1, 2)+pow(sphere->getX2Centroid()-centerX2, 2)+pow(sphere->getX3Centroid()-centerX3, 2)) < d );
+ result = result || (sqrt(pow(sphere->getX1Centroid()-centerX1, 2)+pow(sphere->getX2Centroid()-centerX2, 2)+pow(sphere->getX3Centroid()-centerX3, 2)) < d);
}
}
std::vector<int> values;
@@ -464,12 +596,12 @@ void DemCoProcessor::distributeIDs()
for (int i = 0; i < interactors.size(); i++)
{
- std::map<int, int>::const_iterator it;
+ std::map<int, int>::const_iterator it;
if ((it=idMap.find(interactors[i]->getID())) == idMap.end())
{
throw UbException(UB_EXARGS, "Interactor ID is invalid! The DEM object may be not in PE domain!");
}
-
+
std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->setId(it->second);
//std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->setId(idMap.find(interactors[i]->getID())->second);
@@ -477,11 +609,11 @@ void DemCoProcessor::distributeIDs()
for (int i = 0; i < physicsEngineGeometrieAdapters.size(); i++)
{
- //physicsEngineSolver->updateGeometry(physicsEngineGeometries[i]);
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive())
- {
- interactors[i]->setPhysicsEngineGeometry(physicsEngineGeometrieAdapters[i]);
- }
+ //physicsEngineSolver->updateGeometry(physicsEngineGeometries[i]);
+ if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive())
+ {
+ interactors[i]->setPhysicsEngineGeometry(physicsEngineGeometrieAdapters[i]);
+ }
}
}
@@ -489,3 +621,53 @@ void DemCoProcessor::setBlockVisitor(std::shared_ptr<BoundaryConditionsBlockVisi
{
this->boundaryConditionsBlockVisitor = boundaryConditionsBlockVisitor;
}
+//////////////////////////////////////////////////////////////////////////
+walberla::pe::RigidBody* DemCoProcessor::getPeGeoObject(walberla::id_t id)
+{
+ std::shared_ptr<walberla::blockforest::BlockForest> forest = std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getBlockForest();
+ std::shared_ptr<walberla::domain_decomposition::BlockDataID> storageId = std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getStorageId();
+
+ for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
+ {
+
+ for (auto bodyIt = walberla::pe::BodyIterator::begin(*blockIt, *storageId.get()); bodyIt != walberla::pe::BodyIterator::end(); ++bodyIt)
+ {
+ if (bodyIt->getID() == id)
+ {
+ walberla::pe::RigidBody* geo = *(bodyIt);
+ return geo;
+ }
+ }
+ }
+ return NULL;
+
+ //walberla::pe::RigidBody* peBody = NULL;
+
+ //std::shared_ptr<walberla::pe::BodyStorage> globalStorage = std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getGlobalBodyStorage();
+
+ //auto bodyIt = globalStorage->find(id);
+ //if (bodyIt != globalStorage->end())
+ //{
+ // return *bodyIt;
+ //}
+
+ //for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
+ //{
+ // walberla::pe::Storage* storage = blockIt->getData< walberla::pe::Storage >(*storageId.get());
+ // walberla::pe::BodyStorage* bodyStorage = &(*storage)[0];
+ // walberla::pe::BodyStorage* bodystorageShadowCopies = &(*storage)[1];
+ //
+ // auto bodyIt = bodyStorage->find(id);
+ // if (bodyIt != bodyStorage->end())
+ // {
+ // return *bodyIt;
+ // }
+
+ // bodyIt = bodystorageShadowCopies->find(id);
+ // if (bodyIt != bodystorageShadowCopies->end())
+ // {
+ // return *bodyIt;
+ // }
+ //}
+ //return NULL;
+}
\ No newline at end of file
diff --git a/source/DemCoupling/DemCoProcessor.h b/source/DemCoupling/DemCoProcessor.h
index 9071084e3..cda842793 100644
--- a/source/DemCoupling/DemCoProcessor.h
+++ b/source/DemCoupling/DemCoProcessor.h
@@ -52,6 +52,8 @@ public:
void getObjectsPropertiesVector(std::vector<double>& p);
void addPeGeo(walberla::pe::RigidBody* peGeo);
void removePeGeo(walberla::pe::RigidBody* peGeo);
+ void addPeShadowGeo(walberla::pe::RigidBody* peGeo);
+ void removePeShadowGeo(walberla::pe::RigidBody* peGeo);
bool isSpheresIntersection(double centerX1, double centerX2, double centerX3, double d);
private:
@@ -60,7 +62,8 @@ private:
void setForcesToObject(SPtr<Grid3D> grid, std::shared_ptr<MovableObjectInteractor> interactor, std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry);
void scaleForcesAndTorques(double scalingFactor);
void calculateDemTimeStep(double step);
- void moveVfGeoObject();
+ void moveVfGeoObjects();
+ walberla::pe::RigidBody* getPeGeoObject(walberla::id_t id);
private:
std::shared_ptr<Communicator> comm;
std::vector<std::shared_ptr<MovableObjectInteractor> > interactors;
@@ -69,6 +72,8 @@ private:
std::vector<std::shared_ptr<PhysicsEngineGeometryAdapter> > physicsEngineGeometrieAdapters;
double intermediateDemSteps;
SPtr<BoundaryConditionsBlockVisitor> boundaryConditionsBlockVisitor;
+ walberla::pe::BodyStorage* bodyStorage; //!< Reference to the central body storage.
+ walberla::pe::BodyStorage* bodyStorageShadowCopies; //!< Reference to the body storage containing body shadow copies.
#ifdef TIMING
UbTimer timer;
diff --git a/source/DemCoupling/MovableObjectInteractor.cpp b/source/DemCoupling/MovableObjectInteractor.cpp
index 562d577cf..49a9ece51 100644
--- a/source/DemCoupling/MovableObjectInteractor.cpp
+++ b/source/DemCoupling/MovableObjectInteractor.cpp
@@ -10,6 +10,7 @@
#include "BCAdapter.h"
#include "BCProcessor.h"
#include "ILBMKernel.h"
+#include "CoordinateTransformation3D.h"
#include "SetBcBlocksBlockVisitor.h"
#include "BoundaryConditionsBlockVisitor.h"
@@ -184,22 +185,54 @@ void MovableObjectInteractor::setBcNodesToFluid()
void MovableObjectInteractor::setBcBlocks()
{
- //SetBcBlocksBlockVisitor v(shared_from_this());
- //this->grid.lock()->accept(v);
- SPtr<GbObject3D> geoObject = this->getGbObject3D();
- std::array<double, 6> AABB ={ geoObject->getX1Minimum(),geoObject->getX2Minimum(),geoObject->getX3Minimum(),geoObject->getX1Maximum(),geoObject->getX2Maximum(),geoObject->getX3Maximum() };
- blockVector.clear();
- UbTupleInt3 blockNX=grid.lock()->getBlockNX();
- double ext = 0.0;
- grid.lock()->getBlocksByCuboid(AABB[0]-(double)val<1>(blockNX)*ext, AABB[1]-(double)val<2>(blockNX)*ext, AABB[2]-(double)val<3>(blockNX)*ext, AABB[3]+(double)val<1>(blockNX)*ext, AABB[4]+(double)val<2>(blockNX)*ext, AABB[5]+(double)val<3>(blockNX)*ext, blockVector);
-
- for(std::shared_ptr<Block3D> block : this->blockVector)
- {
- if (block->getKernel())
- {
- setBCBlock(block);
- }
- }
+ SetBcBlocksBlockVisitor v(shared_from_this());
+ this->grid.lock()->accept(v);
+
+ //////////////////////////////////////////////////////////////////////////
+ //SPtr<GbObject3D> geoObject = this->getGbObject3D();
+ //std::array<double, 6> AABB ={ geoObject->getX1Minimum(),geoObject->getX2Minimum(),geoObject->getX3Minimum(),geoObject->getX1Maximum(),geoObject->getX2Maximum(),geoObject->getX3Maximum() };
+ //blockVector.clear();
+ //UbTupleInt3 blockNX=grid.lock()->getBlockNX();
+ //double ext = 0.0;
+ //grid.lock()->getBlocksByCuboid(AABB[0]-(double)val<1>(blockNX)*ext, AABB[1]-(double)val<2>(blockNX)*ext, AABB[2]-(double)val<3>(blockNX)*ext, AABB[3]+(double)val<1>(blockNX)*ext, AABB[4]+(double)val<2>(blockNX)*ext, AABB[5]+(double)val<3>(blockNX)*ext, blockVector);
+
+ //for(std::shared_ptr<Block3D> block : this->blockVector)
+ //{
+ // if (block->getKernel())
+ // {
+ // setBCBlock(block);
+ // }
+ //}
+ //////////////////////////////////////////////////////////////////////////
+ //SPtr<GbObject3D> geoObject = this->getGbObject3D();
+ //std::array <double, 2> minMax1;
+ //std::array <double, 2> minMax2;
+ //std::array <double, 2> minMax3;
+ //minMax1[0] = geoObject->getX1Minimum();
+ //minMax2[0] = geoObject->getX2Minimum();
+ //minMax3[0] = geoObject->getX3Minimum();
+ //minMax1[1] = geoObject->getX1Maximum();
+ //minMax2[1] = geoObject->getX2Maximum();
+ //minMax3[1] = geoObject->getX3Maximum();
+
+ //SPtr<CoordinateTransformation3D> trafo = grid.lock()->getCoordinateTransformator();
+
+ //for (int x3 = 0; x3 < 2; x3++)
+ // for (int x2 = 0; x2 < 2; x2++)
+ // for (int x1 = 0; x1 < 2; x1++)
+ // {
+ // int ix1 = (int)trafo->transformForwardToX1Coordinate(minMax1[x1], minMax2[x2], minMax3[x3]);
+ // int ix2 = (int)trafo->transformForwardToX2Coordinate(minMax1[x1], minMax2[x2], minMax3[x3]);
+ // int ix3 = (int)trafo->transformForwardToX3Coordinate(minMax1[x1], minMax2[x2], minMax3[x3]);
+ // blockVector.push_back(grid.lock()->getBlock(ix1, ix2, ix3, 0));
+ // }
+ //for(std::shared_ptr<Block3D> block : this->blockVector)
+ //{
+ // if (block->getKernel())
+ // {
+ // setBCBlock(block);
+ // }
+ //}
}
void MovableObjectInteractor::updateVelocityBc()
diff --git a/source/DemCoupling/WritePeBlocksCoProcessor.cpp b/source/DemCoupling/WritePeBlocksCoProcessor.cpp
new file mode 100644
index 000000000..f46208ac9
--- /dev/null
+++ b/source/DemCoupling/WritePeBlocksCoProcessor.cpp
@@ -0,0 +1,88 @@
+#include "WritePeBlocksCoProcessor.h"
+
+#include "basics/writer/WbWriterVtkXmlASCII.h"
+
+#include "D3Q27System.h"
+#include "Block3D.h"
+#include "UbScheduler.h"
+#include "Communicator.h"
+#include "Grid3D.h"
+
+WritePeBlocksCoProcessor::WritePeBlocksCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string& path, WbWriter* const writer, SPtr<Communicator> comm, SPtr<walberla::blockforest::BlockForest> forest) :
+ CoProcessor(grid, s),
+ path(path),
+ writer(writer),
+ comm(comm),
+ forest(forest)
+{
+
+}
+
+WritePeBlocksCoProcessor::~WritePeBlocksCoProcessor()
+{
+
+}
+
+void WritePeBlocksCoProcessor::process(double step)
+{
+ if (scheduler->isDue(step))
+ collectData(step);
+}
+
+void WritePeBlocksCoProcessor::collectData(double step)
+{
+ if (comm->getProcessID() == comm->getRoot())
+ {
+ int istep = int(step);
+ std::vector<std::string> filenames;
+ std::vector< UbTupleFloat3 > nodes;
+ std::vector< UbTupleInt8 > cells;
+ std::vector<std::string> celldatanames;
+
+ celldatanames.push_back("ID");
+ celldatanames.push_back("rank");
+
+ walberla::uint_t rank = 0;
+
+
+ std::vector< std::vector< double > > celldata(celldatanames.size());
+
+ int nr=0;
+
+ for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
+ {
+ walberla::AABB aabb = blockIt->getAABB();
+
+ nodes.push_back(makeUbTuple((float)aabb.xMin(), (float)aabb.yMin(), (float)aabb.zMin()));
+ nodes.push_back(makeUbTuple((float)aabb.xMax(), (float)aabb.yMin(), (float)aabb.zMin()));
+ nodes.push_back(makeUbTuple((float)aabb.xMax(), (float)aabb.yMax(), (float)aabb.zMin()));
+ nodes.push_back(makeUbTuple((float)aabb.xMin(), (float)aabb.yMax(), (float)aabb.zMin()));
+ nodes.push_back(makeUbTuple((float)aabb.xMin(), (float)aabb.yMin(), (float)aabb.zMax()));
+ nodes.push_back(makeUbTuple((float)aabb.xMax(), (float)aabb.yMin(), (float)aabb.zMax()));
+ nodes.push_back(makeUbTuple((float)aabb.xMax(), (float)aabb.yMax(), (float)aabb.zMax()));
+ nodes.push_back(makeUbTuple((float)aabb.xMin(), (float)aabb.yMax(), (float)aabb.zMax()));
+ cells.push_back(makeUbTuple(nr, nr+1, nr+2, nr+3, nr+4, nr+5, nr+6, nr+7));
+ nr += 8;
+
+ //data
+ celldata[0].push_back((double)blockIt->getId().getID());
+ forest->getProcessRank(rank,blockIt->getId());
+ celldata[1].push_back((double)rank);
+ }
+
+ filenames.push_back(writer->writeOctsWithCellData(path+"/peBlocks/peBlocks_" + UbSystem::toString(grid->getRank()) + "_" + UbSystem::toString(istep), nodes, cells, celldatanames, celldata));
+
+ if (istep == CoProcessor::scheduler->getMinBegin())
+ {
+ WbWriterVtkXmlASCII::getInstance()->writeCollection(path+"/peBlocks/peBlocks_collection", filenames, istep, false);
+ }
+ else
+ {
+ WbWriterVtkXmlASCII::getInstance()->addFilesToCollection(path + "/peBlocks/peBlocks_collection", filenames, istep, false);
+ }
+
+ UBLOG(logINFO, "WritePeBlocksCoProcessor step: " << istep);
+
+ }
+
+}
\ No newline at end of file
diff --git a/source/DemCoupling/WritePeBlocksCoProcessor.h b/source/DemCoupling/WritePeBlocksCoProcessor.h
new file mode 100644
index 000000000..6d1022566
--- /dev/null
+++ b/source/DemCoupling/WritePeBlocksCoProcessor.h
@@ -0,0 +1,43 @@
+/*
+* WritePeBlocksCoProcessor.h
+*
+* Created on: 07.09.2018
+* Author: K. Kutscher
+*/
+
+#ifndef WritePeBlocksCoProcessor_H_
+#define WritePeBlocksCoProcessor_H_
+
+#include <PointerDefinitions.h>
+#include <string>
+
+#include "CoProcessor.h"
+
+#include <pe/basic.h>
+
+class Communicator;
+class Grid3D;
+class UbScheduler;
+class WbWriter;
+
+class WritePeBlocksCoProcessor : public CoProcessor
+{
+public:
+ WritePeBlocksCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string& path, WbWriter* const writer, SPtr<Communicator> comm, SPtr<walberla::blockforest::BlockForest> forest);
+ virtual ~WritePeBlocksCoProcessor();
+
+ void process(double step) override;
+
+protected:
+ void collectData(double step);
+
+ std::string path;
+ WbWriter* writer;
+ SPtr<Communicator> comm;
+ SPtr<walberla::blockforest::BlockForest> forest;
+};
+
+
+
+#endif
+
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.cpp b/source/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.cpp
index 9283b3b15..451dc027b 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.cpp
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.cpp
@@ -2,6 +2,7 @@
#include "Grid3D.h"
#include "Block3D.h"
#include "CoordinateTransformation3D.h"
+#include "UbLogger.h"
#include "core/debug/CheckFunctions.h"
@@ -19,11 +20,17 @@ walberla::uint_t PeLoadBalancerAdapter::operator()(walberla::SetupBlockForest &
for (auto peBlock = peBlocks.begin(); peBlock != peBlocks.end(); ++peBlock)
{
walberla::AABB aabb = (*peBlock)->getAABB();
- SPtr<Block3D> block = getBlockByMinUniform((double)aabb.xMin(), (double)aabb.yMin(), (double)aabb.zMin(), grid);
+ SPtr<Block3D> block = getBlockByMinUniform(aabb.xMin()+0.5*(aabb.xMax()-aabb.xMin()), aabb.yMin()+0.5*(aabb.yMax()-aabb.yMin()), aabb.zMin()+0.5*(aabb.zMax()-aabb.zMin()), grid);
if (block)
{
(*peBlock)->assignTargetProcess((walberla::uint_t)block->getRank());
}
+ //else
+ //{
+ ////TODO: the rank of pe blocks is not consistent with VF blocks
+ //(*peBlock)->assignTargetProcess(0);
+ ////UBLOG(logINFO, "PeLoadBalancerAdapter::operator() peBlockId="<<(*peBlock)->getId());
+ //}
}
return numberOfProcesses;
@@ -38,4 +45,4 @@ SPtr<Block3D> PeLoadBalancerAdapter::getBlockByMinUniform(double minX1, double m
int ix3 = (int)trafo->transformForwardToX3Coordinate(minX1, minX2, minX3);
return grid->getBlock(ix1, ix2, ix3, 0);
-}
\ No newline at end of file
+}
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp
index ba87d3bfa..fe23383f0 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp
@@ -15,87 +15,89 @@ PePhysicsEngineGeometryAdapter::PePhysicsEngineGeometryAdapter()
{
this->id = -999;
this->active = false;
+ shadowCounter = 0;
+ counter = 0;
}
void PePhysicsEngineGeometryAdapter::addForce(const Vector3D& force)
{
- peGeoObject->addForce(PeConverter::convert(force));
+ peGeoObject->addForce(PeConverter::convert(force));
}
void PePhysicsEngineGeometryAdapter::addTorque(const Vector3D& torque)
{
- peGeoObject->addTorque(PeConverter::convert(torque));
+ peGeoObject->addTorque(PeConverter::convert(torque));
}
void PePhysicsEngineGeometryAdapter::setForce(const Vector3D& force)
{
- peGeoObject->setForce(PeConverter::convert(force));
+ peGeoObject->setForce(PeConverter::convert(force));
}
void PePhysicsEngineGeometryAdapter::setTorque(const Vector3D& torque)
{
- peGeoObject->setTorque(PeConverter::convert(torque));
+ peGeoObject->setTorque(PeConverter::convert(torque));
}
void PePhysicsEngineGeometryAdapter::addForceAtPosition(const Vector3D& force, const Vector3D& position)
{
- peGeoObject->addForceAtPos(PeConverter::convert(force), PeConverter::convert(position));
+ peGeoObject->addForceAtPos(PeConverter::convert(force), PeConverter::convert(position));
}
void PePhysicsEngineGeometryAdapter::setLinearVelolocity(const Vector3D& velocity)
{
- peGeoObject->setLinearVel(PeConverter::convert(velocity));
+ peGeoObject->setLinearVel(PeConverter::convert(velocity));
}
void PePhysicsEngineGeometryAdapter::setAngularVelocity(const Vector3D& velocity)
{
- peGeoObject->setAngularVel(PeConverter::convert(velocity));
+ peGeoObject->setAngularVel(PeConverter::convert(velocity));
}
void PePhysicsEngineGeometryAdapter::resetForceAndTorque()
{
- peGeoObject->resetForceAndTorque();
+ peGeoObject->resetForceAndTorque();
}
Vector3D PePhysicsEngineGeometryAdapter::getVelocityAtPosition(const Vector3D& position) const
{
- return PeConverter::convert(peGeoObject->velFromWF(PeConverter::convert(position)));
+ return PeConverter::convert(peGeoObject->velFromWF(PeConverter::convert(position)));
}
Vector3D PePhysicsEngineGeometryAdapter::getLinearVelocity() const
{
- return PeConverter::convert(peGeoObject->getLinearVel());
+ return PeConverter::convert(peGeoObject->getLinearVel());
}
Vector3D PePhysicsEngineGeometryAdapter::getAngularVelocity() const
{
- return PeConverter::convert(peGeoObject->getAngularVel());
+ return PeConverter::convert(peGeoObject->getAngularVel());
}
Vector3D PePhysicsEngineGeometryAdapter::getPosition() const
{
- return PeConverter::convert(peGeoObject->getPosition());
+ return PeConverter::convert(peGeoObject->getPosition());
}
Vector3D PePhysicsEngineGeometryAdapter::getForce() const
{
- return PeConverter::convert(peGeoObject->getForce());
+ return PeConverter::convert(peGeoObject->getForce());
}
Vector3D PePhysicsEngineGeometryAdapter::getTorque() const
{
- return PeConverter::convert(peGeoObject->getTorque());
+ return PeConverter::convert(peGeoObject->getTorque());
}
void PePhysicsEngineGeometryAdapter::changeState(State state)
{
- if (state == State::PIN)
- peGeoObject->setMassAndInertiaToInfinity();
+ if (state == State::PIN)
+ peGeoObject->setMassAndInertiaToInfinity();
}
int PePhysicsEngineGeometryAdapter::getId() const
{
- return id;
+ return id;
}
void PePhysicsEngineGeometryAdapter::setId(int id)
@@ -105,7 +107,7 @@ void PePhysicsEngineGeometryAdapter::setId(int id)
void PePhysicsEngineGeometryAdapter::setGeometry(walberla::pe::RigidBody* peGeoObject)
{
- this->peGeoObject = peGeoObject;
+ this->peGeoObject = peGeoObject;
}
//////////////////////////////////////////////////////////////////////////
@@ -121,5 +123,7 @@ void PePhysicsEngineGeometryAdapter::setInactive()
//////////////////////////////////////////////////////////////////////////
bool PePhysicsEngineGeometryAdapter::isActive()
{
+
+
return active;
-}
\ No newline at end of file
+}
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h
index 415b73a60..46a81f4fa 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h
@@ -51,12 +51,17 @@ public:
void setActive();
void setInactive();
bool isActive();
+ //void increaseShadowCounter();
+ //void decreaseShad
+ int shadowCounter;
+ int counter;
private:
walberla::pe::RigidBody* peGeoObject;
//unsigned long long id;
int id;
bool active;
+
};
#endif
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
index cb4ee4bfd..5c667d33f 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
@@ -24,11 +24,11 @@ PePhysicsEngineSolverAdapter::PePhysicsEngineSolverAdapter(std::shared_ptr<PePar
void PePhysicsEngineSolverAdapter::initalizePeEnvironment()
{
this->initialPeBodyStorage();
- this->initalPeBlockForest();
+ this->initialPeBlockForest();
this->initalBlockData();
this->initalPeIntegrator();
this->executePeBodyTypeTuple();
- this->initalPeChannel();
+ this->initialPeChannel();
}
@@ -54,22 +54,15 @@ std::shared_ptr<PhysicsEngineGeometryAdapter> PePhysicsEngineSolverAdapter::crea
return peGeometryAdapter;
}
- //if (peGeometry)
- //{
- // return std::static_pointer_cast<PhysicsEngineGeometryAdapter>(std::make_shared<PePhysicsEngineGeometryAdapter>(peGeometry));
- //}
- //else
- //{
- // return std::shared_ptr<PhysicsEngineGeometryAdapter>(new PePhysicsEngineGeometryAdapter());
- //}
-
walberla::pe::syncNextNeighbors<BodyTypeTuple>(*forest, *storageId);
+ //walberla::pe::syncShadowOwners<BodyTypeTuple>(*forest, *storageId);
}
void PePhysicsEngineSolverAdapter::runTimestep(double step)
{
cr->timestep(walberla::real_c(step));
walberla::pe::syncNextNeighbors<BodyTypeTuple>(*forest, *storageId);
+ //walberla::pe::syncShadowOwners<BodyTypeTuple>(*forest, *storageId);
}
@@ -79,19 +72,8 @@ void PePhysicsEngineSolverAdapter::initialPeBodyStorage()
globalBodyStorage = std::make_shared<walberla::pe::BodyStorage>();
}
-void PePhysicsEngineSolverAdapter::initalPeBlockForest()
+void PePhysicsEngineSolverAdapter::initialPeBlockForest()
{
- //forest = walberla::pe::createBlockForest
- //(
- // //walberla::AABB(0, 0, 0, val<1>(peParameter->simulationDomain), val<2>(peParameter->simulationDomain), val<3>(peParameter->simulationDomain)), // simulationDomain
- // walberla::AABB(peParameter->simulationDomain[0], peParameter->simulationDomain[1], peParameter->simulationDomain[2],
- // peParameter->simulationDomain[3], peParameter->simulationDomain[4], peParameter->simulationDomain[5]), // simulationDomain
- // walberla::Vector3<walberla::uint_t>(val<1>(peParameter->numberOfBlocks), val<2>(peParameter->numberOfBlocks), val<3>(peParameter->numberOfBlocks)), // blocks in each direction
- // walberla::Vector3<bool>(val<1>(peParameter->isPeriodic), val<2>(peParameter->isPeriodic), val<3>(peParameter->isPeriodic)) // periodicity
- //);
-
-
-
walberla::SetupBlockForest sforest;
//sforest.addWorkloadMemorySUIDAssignmentFunction( uniformWorkloadAndMemoryAssignment );
sforest.init(walberla::AABB(peParameter->simulationDomain[0], peParameter->simulationDomain[1], peParameter->simulationDomain[2],
@@ -132,7 +114,7 @@ void PePhysicsEngineSolverAdapter::executePeBodyTypeTuple()
walberla::pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
}
-void PePhysicsEngineSolverAdapter::initalPeChannel() const
+void PePhysicsEngineSolverAdapter::initialPeChannel() const
{
const walberla::pe::MaterialID material = peParameter->planes->getPeMaterial();
@@ -207,3 +189,8 @@ std::shared_ptr<walberla::domain_decomposition::BlockDataID> PePhysicsEngineSolv
{
return storageId;
}
+
+std::shared_ptr<walberla::pe::BodyStorage> PePhysicsEngineSolverAdapter::getGlobalBodyStorage()
+{
+ return globalBodyStorage;
+}
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h
index 1f49f0632..21f674edf 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h
@@ -75,16 +75,17 @@ public:
void loadFromFile(const std::string& path);
std::shared_ptr<walberla::blockforest::BlockForest> getBlockForest();
std::shared_ptr<walberla::domain_decomposition::BlockDataID> getStorageId();
+ std::shared_ptr<walberla::pe::BodyStorage> getGlobalBodyStorage();
private:
void initalizePeEnvironment();
void initialPeBodyStorage();
- void initalPeBlockForest();
+ void initialPeBlockForest();
void initalBlockData();
void initalPeIntegrator();
static void executePeBodyTypeTuple();
- void initalPeChannel() const;
+ void initialPeChannel() const;
private:
std::shared_ptr<PeParameter> peParameter;
diff --git a/source/VirtualFluidsCore/CoProcessors/WriteBlocksCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/WriteBlocksCoProcessor.cpp
index d8af5f547..f27d2d59f 100644
--- a/source/VirtualFluidsCore/CoProcessors/WriteBlocksCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/WriteBlocksCoProcessor.cpp
@@ -148,6 +148,6 @@ void WriteBlocksCoProcessor::collectData(double step)
WbWriterVtkXmlASCII::getInstance()->addFilesToCollection(path + "/blocks/blocks_collection", filenames, istep, false);
}
- UBLOG(logINFO,"BlocksCoProcessor step: " << istep);
+ UBLOG(logINFO,"WriteBlocksCoProcessor step: " << istep);
}
}
diff --git a/source/VirtualFluidsCore/CoProcessors/WriteBlocksCoProcessor.h b/source/VirtualFluidsCore/CoProcessors/WriteBlocksCoProcessor.h
index 6ed84d9a1..041d29586 100644
--- a/source/VirtualFluidsCore/CoProcessors/WriteBlocksCoProcessor.h
+++ b/source/VirtualFluidsCore/CoProcessors/WriteBlocksCoProcessor.h
@@ -1,12 +1,12 @@
/*
-* BlocksCoProcessor.h
+* WriteBlocksCoProcessor.h
*
* Created on: 24.09.2012
* Author: K. Kucher
*/
-#ifndef BlocksCoProcessor_H_
-#define BlocksCoProcessor_H_
+#ifndef WriteBlocksCoProcessor_H_
+#define WriteBlocksCoProcessor_H_
#include <PointerDefinitions.h>
#include <string>
diff --git a/source/VirtualFluidsCore/Interactors/InteractorsHelper.cpp b/source/VirtualFluidsCore/Interactors/InteractorsHelper.cpp
index c8aea0034..2d11b43f9 100644
--- a/source/VirtualFluidsCore/Interactors/InteractorsHelper.cpp
+++ b/source/VirtualFluidsCore/Interactors/InteractorsHelper.cpp
@@ -9,8 +9,8 @@
#include "SetBcBlocksBlockVisitor.h"
-InteractorsHelper::InteractorsHelper(SPtr<Grid3D> grid, SPtr<Grid3DVisitor> visitor) :
- grid(grid), visitor(visitor)
+InteractorsHelper::InteractorsHelper(SPtr<Grid3D> grid, SPtr<Grid3DVisitor> visitor, bool deleteBlocks) :
+ grid(grid), visitor(visitor), deleteBlocks(deleteBlocks)
{
}
@@ -50,16 +50,17 @@ void InteractorsHelper::deleteSolidBlocks()
{
for(SPtr<Interactor3D> interactor : interactors)
{
- //setBlocks(interactor, SetBcBlocksBlockVisitor::SOLID);
SetSolidBlocksBlockVisitor v(interactor);
grid->accept(v);
-
- std::vector<SPtr<Block3D>>& sb = interactor->getSolidBlockSet();
- solidBlocks.insert(solidBlocks.end(), sb.begin(), sb.end());
- interactor->removeSolidBlocks();
+ if (deleteBlocks)
+ {
+ std::vector<SPtr<Block3D>>& sb = interactor->getSolidBlockSet();
+ solidBlocks.insert(solidBlocks.end(), sb.begin(), sb.end());
+ interactor->removeSolidBlocks();
+ }
}
- updateGrid();
+ if (deleteBlocks) updateGrid();
}
//////////////////////////////////////////////////////////////////////////
void InteractorsHelper::setBcBlocks()
diff --git a/source/VirtualFluidsCore/Interactors/InteractorsHelper.h b/source/VirtualFluidsCore/Interactors/InteractorsHelper.h
index 2b10c626d..6380ba132 100644
--- a/source/VirtualFluidsCore/Interactors/InteractorsHelper.h
+++ b/source/VirtualFluidsCore/Interactors/InteractorsHelper.h
@@ -13,7 +13,7 @@ class Grid3DVisitor;
class InteractorsHelper
{
public:
- InteractorsHelper(SPtr<Grid3D> grid, SPtr<Grid3DVisitor> visitor);
+ InteractorsHelper(SPtr<Grid3D> grid, SPtr<Grid3DVisitor> visitor, bool deleteBlocks=true);
~InteractorsHelper();
void addInteractor(SPtr<Interactor3D> interactor);
@@ -32,6 +32,7 @@ private:
SPtr<Grid3D> grid;
std::vector<SPtr<Block3D> > solidBlocks;
SPtr<Grid3DVisitor> visitor;
+ bool deleteBlocks;
};
#endif
--
GitLab