diff --git a/source/Applications/Thermoplast/config.txt b/source/Applications/Thermoplast/config.txt
index 81b95d7599743b6f970b5d7982fa93589f1252a7..79fc3464f6477a06c460ab4730d854f21c83a86d 100644
--- a/source/Applications/Thermoplast/config.txt
+++ b/source/Applications/Thermoplast/config.txt
@@ -7,26 +7,26 @@ boundingBox = 60 1370 10 190 1530 320 #test bb
blocknx = 10 10 10
#blocknx = 300 420 320
-endTime = 100
-outTime = 100
+endTime = 100000
+outTime = 1000
availMem = 25e9
uLB = 0.1
#PE parameters
peMinOffset = 46 2 2
peMaxOffset = -8 -25 -2
-sphereTime = 100
+sphereTime = 10
#geometry files
pathGeo = d:/Projects/ThermoPlast/SimGeo
michel = /michel.stl
plexiglas = /plexiglas.stl
-pathOut = g:/temp/thermoplastCluster3
+pathOut = g:/temp/thermoplastCluster5
#restart
-cpStart = 200
-cpStep = 200
+cpStart = 20000000
+cpStep = 20000000
restart = false
nupsStep = 1000 1000 10000000
\ No newline at end of file
diff --git a/source/Applications/Thermoplast/thermoplast.cpp b/source/Applications/Thermoplast/thermoplast.cpp
index 0ffd6c3e445efbbc3d2407d13b7c541ca5eaf016..3ae47c737c53f8e2c0f669ab6931e31ba37ede91 100644
--- a/source/Applications/Thermoplast/thermoplast.cpp
+++ b/source/Applications/Thermoplast/thermoplast.cpp
@@ -56,7 +56,8 @@ std::shared_ptr<DemCoProcessor> makePeCoProcessor(SPtr<Grid3D> grid, SPtr<Commun
//int maxpeIterations = 10000;
//Beschleunigung g
double g = 9.81 * lbmUnitConverter->getFactorAccWToLb();
- Vector3D globalLinearAcc(0.0, -g, 0.0);
+ //Vector3D globalLinearAcc(0.0, -g, 0.0);
+ Vector3D globalLinearAcc(0.0, 0.0, -g);
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);
@@ -162,9 +163,11 @@ void thermoplast(string configname)
fct.DefineConst("U", uLB);
SPtr<BCAdapter> inflowAdapter(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
inflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityBCAlgorithm()));
+ //inflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
SPtr<BCAdapter> outflowAdapter(new DensityBCAdapter(rhoLB));
- outflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonEqDensityBCAlgorithm()));
+ outflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new EqDensityBCAlgorithm()));
+ //outflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonEqDensityBCAlgorithm()));
//outflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonReflectingOutflowBCAlgorithm()));
//sphere BC
@@ -244,15 +247,25 @@ 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");
+
+
+
//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 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());
+ //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());
- GbCuboid3DPtr geoInflow4(new GbCuboid3D(g_minX1-blockLength, g_minX2+4.0*radius, g_maxX3-4.0*radius-1.0, g_minX1+1, g_minX2+4.0*radius+1.0, g_maxX3+radius));
- if (myid == 0) GbSystem3D::writeGeoObject(geoInflow4.get(), pathOut + "/geo/geoInflow4", WbWriterVtkXmlASCII::getInstance());
+ //GbCuboid3DPtr geoInflow4(new GbCuboid3D(g_minX1-blockLength, g_minX2+4.0*radius, g_maxX3-4.0*radius-1.0, g_minX1+1, g_minX2+4.0*radius+1.0, g_maxX3+radius));
+ //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));
@@ -265,8 +278,8 @@ void thermoplast(string configname)
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> inflowInt3 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow3, grid, noSlipBCAdapter, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> inflowInt4 = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow4, grid, noSlipBCAdapter, 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));
@@ -277,12 +290,23 @@ 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));
+
//////////////////////////////////////////////////////////////////////////
- SPtr<Grid3DVisitor> peVisitor(new PePartitioningGridVisitor(comm, demCoProcessor));
+ //SPtr<Grid3DVisitor> peVisitor(new PePartitioningGridVisitor(comm, demCoProcessor));
+ SPtr<Grid3DVisitor> peVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::KWAY));
InteractorsHelper intHelper(grid, peVisitor);
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);
@@ -331,14 +355,14 @@ void thermoplast(string configname)
grid->accept(initVisitor);
//write data for visualization of boundary conditions
- //{
- // SPtr<UbScheduler> geoSch(new UbScheduler(1));
- // WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), comm);
- // ppgeo.process(0);
+ {
+ SPtr<UbScheduler> geoSch(new UbScheduler(1));
+ WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), comm);
+ ppgeo.process(0);
- // WriteMacroscopicQuantitiesCoProcessor ppInit(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm);
- // ppInit.process(0);
- //}
+ WriteMacroscopicQuantitiesCoProcessor ppInit(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm);
+ ppInit.process(0);
+ }
if (myid == 0) UBLOG(logINFO, "Preprocess - end");
}
@@ -381,28 +405,28 @@ 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));
- // }
- //UBLOG(logINFO, "sphere prototypes - stop, rank="<<myid);
-
- Vector3D origin(106+radius, 1372+radius, 12+radius);
- for (int x3 = 0; x3 < 28; x3++)
- for (int x2 = 0; x2 < 12; x2++)
- for (int x1 = 0; x1 < 7; x1++)
+ 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]+x1*1.1*d, origin[1]+x2*1.1*d, origin[2]+x3*1.1*d, radius));
- //if (myid == 0) GbSystem3D::writeGeoObject(sphere.get(), pathOut + "/geo/sphere"+UbSystem::toString(x1)+UbSystem::toString(x2)+UbSystem::toString(x3), WbWriterVtkXmlASCII::getInstance());
+ 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));
}
+ //UBLOG(logINFO, "sphere prototypes - stop, rank="<<myid);
+
+ //Vector3D origin(106+radius, 1372+radius, 12+radius);
+ //for (int x3 = 0; x3 < 28; x3++)
+ // for (int x2 = 0; x2 < 12; x2++)
+ // for (int x1 = 0; x1 < 7; 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]+x1*1.1*d, origin[1]+x2*1.1*d, origin[2]+x3*1.1*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));
+ // }
createSphereCoProcessor->process(0);
@@ -410,6 +434,7 @@ void thermoplast(string configname)
SPtr<UbScheduler> visSch(new UbScheduler(outTime));
SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, pathOut,
WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
+ writeMQCoProcessor->process(0);
SPtr<WriteBoundaryConditionsCoProcessor> writeBCCoProcessor(new WriteBoundaryConditionsCoProcessor(grid, visSch, pathOut,
WbWriterVtkXmlBinary::getInstance(), comm));
diff --git a/source/DemCoupling/CreateDemObjectsCoProcessor.cpp b/source/DemCoupling/CreateDemObjectsCoProcessor.cpp
index 87d7235c03b39237d84030d1ff3d417005623bf7..5c6a43e6c02d4ce6d8cb57f123b60a244f85775c 100644
--- a/source/DemCoupling/CreateDemObjectsCoProcessor.cpp
+++ b/source/DemCoupling/CreateDemObjectsCoProcessor.cpp
@@ -56,11 +56,11 @@ void CreateDemObjectsCoProcessor::process(double step)
if (comm->isRoot()) UBLOG(logINFO, "number of objects = "<<(int)(geoObjectPrototypeVector.size()));
#endif
- demCoProcessor->distributeIDs();
-
-#ifdef TIMING
- if (comm->isRoot()) UBLOG(logINFO, "demCoProcessor->distributeIDs() time = "<<timer.stop()<<" s");
-#endif
+// demCoProcessor->distributeIDs();
+//
+//#ifdef TIMING
+// if (comm->isRoot()) UBLOG(logINFO, "demCoProcessor->distributeIDs() time = "<<timer.stop()<<" s");
+//#endif
if (comm->isRoot())
UBLOG(logINFO, "CreateDemObjectsCoProcessor::process stop step: " << istep);
@@ -83,19 +83,36 @@ void CreateDemObjectsCoProcessor::createGeoObjects()
{
int size = (int)(geoObjectPrototypeVector.size());
+ std::vector< std::shared_ptr<Block3D> > blockVector;
+
for (int i = 0; i < size; i++)
{
//timer.resetAndStart();
std::array<double, 6> AABB ={ geoObjectPrototypeVector[i]->getX1Minimum(),geoObjectPrototypeVector[i]->getX2Minimum(),geoObjectPrototypeVector[i]->getX3Minimum(),geoObjectPrototypeVector[i]->getX1Maximum(),geoObjectPrototypeVector[i]->getX2Maximum(),geoObjectPrototypeVector[i]->getX3Maximum() };
//UBLOG(logINFO, "demCoProcessor->isGeoObjectInAABB(AABB) = " << demCoProcessor->isGeoObjectInAABB(AABB));
- if (demCoProcessor->isDemObjectInAABB(AABB))
+ //if (demCoProcessor->isDemObjectInAABB(AABB))
+ //{
+ // continue;
+ //}
+ SPtr<GbSphere3D> sphere = std::dynamic_pointer_cast<GbSphere3D>(geoObjectPrototypeVector[i]);
+ if (demCoProcessor->isSpheresIntersection(sphere->getX1Centroid(), sphere->getX2Centroid(), sphere->getX3Centroid(), sphere->getRadius()*2.0))
{
continue;
}
+
SPtr<GbObject3D> geoObject((GbObject3D*)(geoObjectPrototypeVector[i]->clone()));
SPtr<MovableObjectInteractor> geoObjectInt = SPtr<MovableObjectInteractor>(new MovableObjectInteractor(geoObject, grid, velocityBcParticleAdapter, Interactor3D::SOLID, extrapolationReconstructor, State::UNPIN));
- SetBcBlocksBlockVisitor setBcVisitor(geoObjectInt);
- grid->accept(setBcVisitor);
+
+ //SetBcBlocksBlockVisitor setBcVisitor(geoObjectInt);
+ //grid->accept(setBcVisitor);
+
+ //std::vector< std::shared_ptr<Block3D> > blockVector;
+ blockVector.clear();
+ grid->getBlocksByCuboid(AABB[0], AABB[1], AABB[2], AABB[3], AABB[4], AABB[5], blockVector);
+ for (std::shared_ptr<Block3D> block : blockVector)
+ geoObjectInt->setBCBlock(block);
+
+
//UBLOG(logINFO, "grid->accept(setBcVisitor) time = "<<timer.stop());
geoObjectInt->initInteractor();
//UBLOG(logINFO, "geoObjectInt->initInteractor() time = "<<timer.stop());
diff --git a/source/DemCoupling/DemCoProcessor.cpp b/source/DemCoupling/DemCoProcessor.cpp
index 88ea7f910ed0a2e1c557b3672dba6730db170b46..8ac7a686f55e570b87fa6870d977969a85f530da 100644
--- a/source/DemCoupling/DemCoProcessor.cpp
+++ b/source/DemCoupling/DemCoProcessor.cpp
@@ -73,15 +73,15 @@ std::shared_ptr<PhysicsEngineGeometryAdapter> DemCoProcessor::createPhysicsEngin
const Vector3D position(vfSphere->getX1Centroid(), vfSphere->getX2Centroid(), vfSphere->getX3Centroid());
auto peGeometry = this->physicsEngineSolver->createPhysicsEngineGeometryAdapter(id, position, vfSphere->getRadius(), physicsEngineMaterial);
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(peGeometry)->isActive())
- {
+ //if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(peGeometry)->isActive())
+ //{
interactor->setPhysicsEngineGeometry(peGeometry);
return peGeometry;
- }
- else
- {
- return peGeometry;
- }
+ //}
+ //else
+ //{
+ // return peGeometry;
+ //}
}
@@ -233,17 +233,17 @@ void DemCoProcessor::calculateDemTimeStep(double step)
if (comm->isRoot()) UBLOG(logINFO, " physicsEngineSolver->runTimestep() time = "<< timer.stop() <<" s");
#endif
- for (int i = 0; i < physicsEngineGeometries.size(); i++)
- {
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
- {
- interactors[i]->setPhysicsEngineGeometry(physicsEngineGeometries[i]);
- }
- }
+ //for (int i = 0; i < physicsEngineGeometries.size(); i++)
+ //{
+ // if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
+ // {
+ // interactors[i]->setPhysicsEngineGeometry(physicsEngineGeometries[i]);
+ // }
+ //}
-#ifdef TIMING
- if (comm->isRoot()) UBLOG(logINFO, " physicsEngineSolver->updateGeometry() time = "<<timer.stop()<<" s");
-#endif
+//#ifdef TIMING
+// if (comm->isRoot()) UBLOG(logINFO, " physicsEngineSolver->updateGeometry() time = "<<timer.stop()<<" s");
+//#endif
}
void DemCoProcessor::moveVfGeoObject()
@@ -386,56 +386,85 @@ void DemCoProcessor::removePeGeo(walberla::pe::RigidBody * peGeo)
throw UbException(UB_EXARGS, "PeGeo ID is not matching!");
}
-void DemCoProcessor::distributeIDs()
+bool DemCoProcessor::isSpheresIntersection(double centerX1, double centerX2, double centerX3, double d)
{
- std::vector<int> peIDsSend;
- std::vector<int> vfIDsSend;
-
+ bool result = false;
for (int i = 0; i < interactors.size(); i++)
{
if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
{
- peIDsSend.push_back(std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->getId());
- vfIDsSend.push_back(interactors[i]->getID());
+ 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 );
}
}
+ std::vector<int> values;
+ values.push_back((int)result);
+ std::vector<int> rvalues = comm->gather(values);
- std::vector<int> peIDsRecv;
- std::vector<int> vfIDsRecv;
-
- comm->allGather(peIDsSend, peIDsRecv);
- comm->allGather(vfIDsSend, vfIDsRecv);
-
- std::map<int, int> idMap;
-
- for (int i = 0; i < peIDsRecv.size(); i++)
- {
- idMap.insert(std::make_pair(vfIDsRecv[i], peIDsRecv[i]));
- }
-
- for (int i = 0; i < interactors.size(); i++)
+ if (comm->isRoot())
{
- std::map<int, int>::const_iterator it;
- if ((it=idMap.find(interactors[i]->getID())) == idMap.end())
+ for (int i = 0; i < (int)rvalues.size(); i++)
{
- throw UbException(UB_EXARGS, "Interactor ID is invalid! The DEM object may be not in PE domain!");
+ result = result || (bool)rvalues[i];
}
-
-
- std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->setId(it->second);
- //std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->setId(idMap.find(interactors[i]->getID())->second);
}
+ int iresult = (int)result;
+ comm->broadcast(iresult);
+ result = (bool)iresult;
- for (int i = 0; i < physicsEngineGeometries.size(); i++)
- {
- //physicsEngineSolver->updateGeometry(physicsEngineGeometries[i]);
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
- {
- interactors[i]->setPhysicsEngineGeometry(physicsEngineGeometries[i]);
- }
- }
+ return result;
}
+//void DemCoProcessor::distributeIDs()
+//{
+// std::vector<int> peIDsSend;
+// std::vector<int> vfIDsSend;
+//
+// for (int i = 0; i < interactors.size(); i++)
+// {
+// if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
+// {
+// peIDsSend.push_back(std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->getId());
+// vfIDsSend.push_back(interactors[i]->getID());
+// }
+// }
+//
+// std::vector<int> peIDsRecv;
+// std::vector<int> vfIDsRecv;
+//
+// comm->allGather(peIDsSend, peIDsRecv);
+// comm->allGather(vfIDsSend, vfIDsRecv);
+//
+// std::map<int, int> idMap;
+//
+// for (int i = 0; i < peIDsRecv.size(); i++)
+// {
+// idMap.insert(std::make_pair(vfIDsRecv[i], peIDsRecv[i]));
+// }
+//
+// for (int i = 0; i < interactors.size(); i++)
+// {
+// 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>(physicsEngineGeometries[i])->setId(it->second);
+// //std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->setId(idMap.find(interactors[i]->getID())->second);
+// }
+//
+// for (int i = 0; i < physicsEngineGeometries.size(); i++)
+// {
+// //physicsEngineSolver->updateGeometry(physicsEngineGeometries[i]);
+// if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
+// {
+// interactors[i]->setPhysicsEngineGeometry(physicsEngineGeometries[i]);
+// }
+// }
+//}
+
void DemCoProcessor::setBlockVisitor(std::shared_ptr<BoundaryConditionsBlockVisitor> boundaryConditionsBlockVisitor)
{
this->boundaryConditionsBlockVisitor = boundaryConditionsBlockVisitor;
diff --git a/source/DemCoupling/DemCoProcessor.h b/source/DemCoupling/DemCoProcessor.h
index f85e0e7d130b67909ba024ba02736782af9ae34f..dcc8074f375e734e23ccecc65a9886f5f4168c86 100644
--- a/source/DemCoupling/DemCoProcessor.h
+++ b/source/DemCoupling/DemCoProcessor.h
@@ -45,13 +45,14 @@ public:
void addInteractor(std::shared_ptr<MovableObjectInteractor> interactor, std::shared_ptr<PhysicsEngineMaterialAdapter> physicsEngineMaterial, Vector3D initalVelocity = Vector3D(0.0, 0.0, 0.0));
void process(double step) override;
std::shared_ptr<PhysicsEngineSolverAdapter> getPhysicsEngineSolver();
- void distributeIDs();
+ //void distributeIDs();
void setBlockVisitor(std::shared_ptr<BoundaryConditionsBlockVisitor> blockVisitor);
bool isDemObjectInAABB(std::array<double,6> AABB);
void addSurfaceTriangleSet(std::vector<UbTupleFloat3>& nodes, std::vector<UbTupleInt3>& triangles);
void getObjectsPropertiesVector(std::vector<double>& p);
void addPeGeo(walberla::pe::RigidBody* peGeo);
void removePeGeo(walberla::pe::RigidBody* peGeo);
+ bool isSpheresIntersection(double centerX1, double centerX2, double centerX3, double d);
private:
std::shared_ptr<PhysicsEngineGeometryAdapter> createPhysicsEngineGeometryAdapter(std::shared_ptr<MovableObjectInteractor> interactor, std::shared_ptr<PhysicsEngineMaterialAdapter> physicsEngineMaterial) const;
diff --git a/source/DemCoupling/MovableObjectInteractor.cpp b/source/DemCoupling/MovableObjectInteractor.cpp
index 91156903a42be156e3ab00b3e4ee8db063151103..bf658073def24203853ed1b8f7c4808fd84c9792 100644
--- a/source/DemCoupling/MovableObjectInteractor.cpp
+++ b/source/DemCoupling/MovableObjectInteractor.cpp
@@ -17,11 +17,19 @@
#include "PhysicsEngineGeometryAdapter.h"
#include "Reconstructor.h"
+#include <array>
+
+//#define TIMING
+
+#ifdef TIMING
+ #include "UbTiming.h"
+#endif
+
MovableObjectInteractor::MovableObjectInteractor(std::shared_ptr<GbObject3D> geoObject3D, std::shared_ptr<Grid3D> grid, std::shared_ptr<BCAdapter> bcAdapter, int type, std::shared_ptr<Reconstructor> reconstructor, State state)
: D3Q27Interactor(geoObject3D, grid, bcAdapter, type), reconstructor(reconstructor), state(state)
{
-
+ //grid->getBlocks(0, grid->getRank(), true, blockVector);
}
MovableObjectInteractor::~MovableObjectInteractor()
@@ -45,19 +53,62 @@ void MovableObjectInteractor::moveGbObjectTo(const Vector3D& position)
void MovableObjectInteractor::rearrangeGrid()
{
+#ifdef TIMING
+ UbTimer timer;
+ timer.resetAndStart();
+#endif
+
+#ifdef TIMING
+ UBLOG(logINFO, "MovableObjectInteractor::rearrangeGrid():start");
+#endif
+
this->reconstructDistributionOnSolidNodes();
-
+
+#ifdef TIMING
+ UBLOG(logINFO, "reconstructDistributionOnSolidNodes() time = "<<timer.stop()<<" s");
+#endif
+
this->setSolidNodesToFluid();
+
+#ifdef TIMING
+ UBLOG(logINFO, "setSolidNodesToFluid() time = "<<timer.stop()<<" s");
+#endif
+
this->setBcNodesToFluid();
+#ifdef TIMING
+ UBLOG(logINFO, "setBcNodesToFluid() time = "<<timer.stop()<<" s");
+#endif
+
this->removeSolidBlocks();
+
+#ifdef TIMING
+ UBLOG(logINFO, "removeSolidBlocks() time = "<<timer.stop()<<" s");
+#endif
+
this->removeBcBlocks();
+#ifdef TIMING
+ UBLOG(logINFO, "removeBcBlocks() time = "<<timer.stop()<<" s");
+#endif
+
this->setBcBlocks();
+#ifdef TIMING
+ UBLOG(logINFO, "setBcBlocks() time = "<<timer.stop()<<" s");
+#endif
+
this->initInteractor();
+#ifdef TIMING
+ UBLOG(logINFO, "initInteractor() time = "<<timer.stop()<<" s");
+#endif
+
this->updateVelocityBc();
+
+#ifdef TIMING
+ UBLOG(logINFO, "updateVelocityBc() time = "<<timer.stop()<<" s");
+#endif
}
void MovableObjectInteractor::reconstructDistributionOnSolidNodes()
@@ -116,8 +167,15 @@ void MovableObjectInteractor::setBcNodesToFluid()
void MovableObjectInteractor::setBcBlocks()
{
- SetBcBlocksBlockVisitor v(shared_from_this());
- this->grid.lock()->accept(v);
+ //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();
+ grid.lock()->getBlocksByCuboid(AABB[0], AABB[1], AABB[2], AABB[3], AABB[4], AABB[5], blockVector);
+
+ for(std::shared_ptr<Block3D> block : this->blockVector)
+ setBCBlock(block);
}
void MovableObjectInteractor::updateVelocityBc()
diff --git a/source/DemCoupling/MovableObjectInteractor.h b/source/DemCoupling/MovableObjectInteractor.h
index 2b43d97dc5fc84e62b949e274a968723b13d36f6..a77fb4502de03f91a85997885de1813232461aa4 100644
--- a/source/DemCoupling/MovableObjectInteractor.h
+++ b/source/DemCoupling/MovableObjectInteractor.h
@@ -47,7 +47,7 @@ private:
std::shared_ptr<Reconstructor> reconstructor;
State state;
-
+ std::vector< std::shared_ptr<Block3D> > blockVector;
};
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp
index b3b687222934c75dd409e3fab22d0a1e8b46d89d..ba87d3bfaf1ee40593ee9a8a1d328ae63073b67d 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp
@@ -5,11 +5,11 @@
#include "PeAdapter.h"
-PePhysicsEngineGeometryAdapter::PePhysicsEngineGeometryAdapter(walberla::pe::RigidBody* peGeoObject) : peGeoObject(peGeoObject)
-{
- this->id = peGeoObject->getID();
- this->active = true;
-}
+//PePhysicsEngineGeometryAdapter::PePhysicsEngineGeometryAdapter(walberla::pe::RigidBody* peGeoObject) : peGeoObject(peGeoObject)
+//{
+// this->id = peGeoObject->getID();
+// this->active = true;
+//}
PePhysicsEngineGeometryAdapter::PePhysicsEngineGeometryAdapter()
{
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h
index 6cb56bc5f3ef999d1fb701b41d84f133421a5ec8..415b73a600265b4e96753c2ea14ddf917e395915 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h
@@ -20,7 +20,7 @@ class PePhysicsEngineGeometryAdapter : public PhysicsEngineGeometryAdapter
{
public:
PePhysicsEngineGeometryAdapter();
- PePhysicsEngineGeometryAdapter(walberla::pe::RigidBody* peGeoObject);
+ //PePhysicsEngineGeometryAdapter(walberla::pe::RigidBody* peGeoObject);
virtual ~PePhysicsEngineGeometryAdapter() {}
void addForce(const Vector3D& force) override;
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
index a58973c1d93620c953effc1bb84b03621fd6e337..82f3ff85c6c2ef88fc172bb82c73b49ed71c9b9b 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
@@ -37,15 +37,31 @@ std::shared_ptr<PhysicsEngineGeometryAdapter> PePhysicsEngineSolverAdapter::crea
walberla::pe::GeomID peGeometry = createSphere(*globalBodyStorage, *forest, *storageId, id, PeConverter::convert(position), radius, peMaterial->getPeMaterial());
+ std::shared_ptr<PePhysicsEngineGeometryAdapter> peGeometryAdapter(new PePhysicsEngineGeometryAdapter());
+
if (peGeometry)
{
- return std::static_pointer_cast<PhysicsEngineGeometryAdapter>(std::make_shared<PePhysicsEngineGeometryAdapter>(peGeometry));
- }
+ peGeometryAdapter->setId(id);
+ peGeometryAdapter->setActive();
+ peGeometryAdapter->setGeometry(peGeometry);
+ return peGeometryAdapter;
+ }
else
{
- return std::shared_ptr<PhysicsEngineGeometryAdapter>(new PePhysicsEngineGeometryAdapter());
+ peGeometryAdapter->setId(id);
+ peGeometryAdapter->setInactive();
+ 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);
}
diff --git a/source/VirtualFluidsCore/BoundaryConditions/NonEqDensityBCAlgorithm.cpp b/source/VirtualFluidsCore/BoundaryConditions/NonEqDensityBCAlgorithm.cpp
index 11cadef1a4029eba5c9143ad7c9c43dbcdffc5bd..865e2f1cf2a522a593890e9d3891ff1f31cc44fa 100644
--- a/source/VirtualFluidsCore/BoundaryConditions/NonEqDensityBCAlgorithm.cpp
+++ b/source/VirtualFluidsCore/BoundaryConditions/NonEqDensityBCAlgorithm.cpp
@@ -44,6 +44,10 @@ void NonEqDensityBCAlgorithm::applyBC()
LBMReal rho, vx1, vx2, vx3;
calcMacrosFct(f, rho, vx1, vx2, vx3);
+ //LBMReal vlimit=0.01;
+ //vx1=(fabs(vx1)>vlimit) ? vx1/fabs(vx1)*vlimit : vx1;
+ //vx2=(fabs(vx2)>vlimit) ? vx2/fabs(vx2)*vlimit : vx2;
+ //vx3=(fabs(vx3)>vlimit) ? vx3/fabs(vx3)*vlimit : vx3;
LBMReal rhoBC = bcPtr->getBoundaryDensity();
for (int fdir = D3Q27System::STARTF; fdir<=D3Q27System::ENDF; fdir++)
{
@@ -52,6 +56,7 @@ void NonEqDensityBCAlgorithm::applyBC()
// Martins NEQ ADDON
////original: 15.2.2013:
LBMReal ftemp = calcFeqsForDirFct(fdir, rho, vx1, vx2, vx3);
+ //rhoBC=(rho>rhoBC)? rhoBC : rho; //Limiter 08.08.2018
ftemp = calcFeqsForDirFct(fdir, rhoBC, vx1, vx2, vx3)+f[fdir]-ftemp;
distributions->setDistributionForDirection(ftemp, nx1, nx2, nx3, fdir);
}
diff --git a/source/VirtualFluidsCore/CoProcessors/MPIIOMigrationCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/MPIIOMigrationCoProcessor.cpp
index 9ce1db8d324aff541a16130e6399c0dbd28661a3..a1e8404ecd6e628577d7b72e7d6df1552441fb80 100644
--- a/source/VirtualFluidsCore/CoProcessors/MPIIOMigrationCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/MPIIOMigrationCoProcessor.cpp
@@ -145,6 +145,8 @@ void MPIIOMigrationCoProcessor::process(double step)
writeDataSet((int)step);
writeBoundaryConds((int)step);
+ writeCpTimeStep(step);
+
if (comm->isRoot()) UBLOG(logINFO, "Save check point - end");
}
}