From bd3476c0b29752abdfbac67526efcec368f3a734 Mon Sep 17 00:00:00 2001
From: Konstantin Kutscher <kutscher@irmb.tu-bs.de>
Date: Mon, 27 Aug 2018 14:27:35 +0200
Subject: [PATCH] Optimisation of discretisation of VF-blocks
---
source/Applications/DLR-F16-Porous/f16.cpp | 31 +++-
.../Applications/DLR-F16-Solid/f16-solid.cfg | 8 +-
source/Applications/DLR-F16-Solid/f16.cpp | 17 +-
source/Applications/Thermoplast/config.txt | 2 +-
.../Applications/Thermoplast/thermoplast.cpp | 7 +-
.../CreateDemObjectsCoProcessor.cpp | 22 ++-
.../DemCoupling/CreateDemObjectsCoProcessor.h | 2 +-
source/DemCoupling/DemCoProcessor.cpp | 172 ++++++++++--------
source/DemCoupling/DemCoProcessor.h | 4 +-
.../DemCoupling/MovableObjectInteractor.cpp | 11 +-
source/VirtualFluidsCore/Grid/Grid3D.cpp | 55 ++++++
source/VirtualFluidsCore/Grid/Grid3D.h | 1 +
12 files changed, 218 insertions(+), 114 deletions(-)
diff --git a/source/Applications/DLR-F16-Porous/f16.cpp b/source/Applications/DLR-F16-Porous/f16.cpp
index 39b920152..5f08df3e1 100644
--- a/source/Applications/DLR-F16-Porous/f16.cpp
+++ b/source/Applications/DLR-F16-Porous/f16.cpp
@@ -174,9 +174,9 @@ void initPteBC(SPtr<Grid3D> grid, vector<SPtr<Block3D>>& vectorTE, string fngFil
//////////////////////////////////////////////////////////////////////////
void initPte(SPtr<Grid3D> grid, SPtr<Interactor3D> fngIntrTE, SPtr<Interactor3D> fngIntrTEmesh, string pathOut, SPtr<Communicator> comm)
{
- SetSolidOrBoundaryBlockVisitor v1(fngIntrTE, SetSolidOrBoundaryBlockVisitor::SOLID);
+ SetSolidBlocksBlockVisitor v1(fngIntrTE);
grid->accept(v1);
- SetSolidOrBoundaryBlockVisitor v2(fngIntrTE, SetSolidOrBoundaryBlockVisitor::BC);
+ SetBcBlocksBlockVisitor v2(fngIntrTE);
grid->accept(v2);
std::vector<SPtr<Block3D>>& sb = fngIntrTE->getSolidBlockSet();
std::vector<SPtr<Block3D>>& bb = fngIntrTE->getBcBlocks();
@@ -649,7 +649,7 @@ void run(string configname)
SPtr<Interactor3D> fngIntrNoTapeBody = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshNoTapeBody, grid, noSlipBCAdapter, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy));//, Interactor3D::POINTS));
- SetSolidOrBoundaryBlockVisitor v(fngIntrNoTapeBody, SetSolidOrBoundaryBlockVisitor::SOLID);
+ SetSolidBlocksBlockVisitor v(fngIntrNoTapeBody);
grid->accept(v);
std::vector<SPtr<Block3D>>& sb = fngIntrNoTapeBody->getSolidBlockSet();
for (SPtr<Block3D> block : sb)
@@ -892,9 +892,9 @@ void run(string configname)
if (myid==0) UBLOG(logINFO, "PID = "<<myid<<" Physical Memory currently used by current process: "<<Utilities::getPhysMemUsedByMe()/1073741824.0<<" GB");
if (myid==0) UBLOG(logINFO, "initPteBC:start");
- SetSolidOrBoundaryBlockVisitor v1(fngIntrTE, SetSolidOrBoundaryBlockVisitor::SOLID);
+ SetSolidBlocksBlockVisitor v1(fngIntrTE);
grid->accept(v1);
- SetSolidOrBoundaryBlockVisitor v2(fngIntrTE, SetSolidOrBoundaryBlockVisitor::BC);
+ SetBcBlocksBlockVisitor v2(fngIntrTE);
grid->accept(v2);
std::vector<SPtr<Block3D>>& vectorTE = fngIntrTE->getSolidBlockSet();
std::vector<SPtr<Block3D>>& bb = fngIntrTE->getBcBlocks();
@@ -959,7 +959,7 @@ void run(string configname)
//Post process
{
SPtr<UbScheduler> geoSch(new UbScheduler(1));
- WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), conv, comm);
+ WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), comm);
ppgeo.process(0);
}
@@ -1008,7 +1008,22 @@ void run(string configname)
SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
grid->accept(setConnsVisitor);
- if (reinit)
{
//InitDistributionsBlockVisitor initVisitor1;
//grid->accept(initVisitor1);
SPtr<Grid3D> oldGrid(new Grid3D(comm));
SPtr<UbScheduler> iSch(new UbScheduler());
SPtr<MPIIORestartCoProcessor> rcp(new MPIIORestartCoProcessor(oldGrid, iSch, pathReInit, comm));
rcp->setLBMKernel(kernel);
rcp->setBCProcessor(bcProc);
rcp->restart(stepReInit);
InitDistributionsWithInterpolationGridVisitor initVisitor2(oldGrid, iProcessor, nuLB);
grid->accept(initVisitor2);
//if (myid==0) UBLOG(logINFO, "reinitGrid:start");
//reinitGrid(grid);
//if (myid==0) UBLOG(logINFO, "reinitGrid:end");
}
+ if (reinit)
+ {
+ //InitDistributionsBlockVisitor initVisitor1;
+ //grid->accept(initVisitor1);
+ SPtr<Grid3D> oldGrid(new Grid3D(comm));
+ SPtr<UbScheduler> iSch(new UbScheduler());
+ SPtr<MPIIORestartCoProcessor> rcp(new MPIIORestartCoProcessor(oldGrid, iSch, pathReInit, comm));
+ rcp->setLBMKernel(kernel);
+ rcp->setBCProcessor(bcProc);
+ rcp->restart(stepReInit);
+ InitDistributionsWithInterpolationGridVisitor initVisitor2(oldGrid, iProcessor, nuLB);
+ grid->accept(initVisitor2);
+ //if (myid==0) UBLOG(logINFO, "reinitGrid:start");
+ //reinitGrid(grid);
+ //if (myid==0) UBLOG(logINFO, "reinitGrid:end");
+ }
//if (myid==0) UBLOG(logINFO, "setPointsTE:start");
//SPtr<GbTriFaceMesh3D> fngMeshTE;
@@ -1040,7 +1055,7 @@ void run(string configname)
fngMeshBody->translate(0.0, 0.0, -0.00011);
if (myid==0) GbSystem3D::writeGeoObject(fngMeshBody.get(), pathOut+"/geo/fngMeshBody2", WbWriterVtkXmlBinary::getInstance());
SPtr<Interactor3D> fngIntrBody = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshBody, grid, noSlipBCAdapter, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy));
- SetSolidOrBoundaryBlockVisitor v(fngIntrBody, SetSolidOrBoundaryBlockVisitor::BC);
+ SetBcBlocksBlockVisitor v(fngIntrBody);
grid->accept(v);
fngIntrBody->initInteractor();
//////////////////////////////////////////////////////////////////////////
diff --git a/source/Applications/DLR-F16-Solid/f16-solid.cfg b/source/Applications/DLR-F16-Solid/f16-solid.cfg
index fb55f3e85..304613a23 100644
--- a/source/Applications/DLR-F16-Solid/f16-solid.cfg
+++ b/source/Applications/DLR-F16-Solid/f16-solid.cfg
@@ -1,4 +1,4 @@
-pathOut = d:/temp/DLR-F16-Solid-mic7
+pathOut = d:/temp/DLR-F16-Solid-comp
pathGeo = d:/Projects/SFB880/DLR-F16/Geometry
fngFileWhole1 = F16_broad_Quad_noTape_full.stl
@@ -11,7 +11,7 @@ pathReInit = /work/koskuche/DLR-F16_L7
stepReInit = 10000
numOfThreads = 4
-availMem = 3.5e9
+availMem = 10e9
logToFile = false
@@ -23,9 +23,9 @@ blockNx = 10 10 10
refineLevel = 0
-#deltaXfine = 0.003
+deltaXfine = 0.003
#deltaXfine = 0.0015
-deltaXfine = 0.00075 #level 0
+#deltaXfine = 0.00075 #level 0
#deltaXfine = 0.000375 #level 1
#deltaXfine = 0.0001875 #level 2
#deltaXfine = 0.00009375 #level 3
diff --git a/source/Applications/DLR-F16-Solid/f16.cpp b/source/Applications/DLR-F16-Solid/f16.cpp
index 96f6143a6..27b7bc623 100644
--- a/source/Applications/DLR-F16-Solid/f16.cpp
+++ b/source/Applications/DLR-F16-Solid/f16.cpp
@@ -169,9 +169,12 @@ void run(string configname)
SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CompressibleCumulant4thOrderViscosityLBMKernel());
//dynamicPointerCast<CompressibleCumulant4thOrderViscosityLBMKernel>(kernel)->setBulkViscosity(nuLB*2.0e3);
+ //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new IncompressibleCumulantLBMKernel());
+
kernel->setBCProcessor(bcProc);
SPtr<LBMKernel> spKernel = SPtr<LBMKernel>(new CompressibleCumulantLBMKernel());
+ //SPtr<LBMKernel> spKernel = SPtr<LBMKernel>(new IncompressibleCumulantLBMKernel());
spKernel->setBCProcessor(bcProc);
//////////////////////////////////////////////////////////////////////////
//restart
@@ -745,15 +748,15 @@ void run(string configname)
SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
calculator->addCoProcessor(nupsCoProcessor);
calculator->addCoProcessor(restartCoProcessor);
- calculator->addCoProcessor(writeMQSelectCoProcessor);
+ //calculator->addCoProcessor(writeMQSelectCoProcessor);
calculator->addCoProcessor(writeMQCoProcessor);
- calculator->addCoProcessor(tsp1);
- calculator->addCoProcessor(tsp2);
- calculator->addCoProcessor(tsp3);
- calculator->addCoProcessor(tsp4);
- calculator->addCoProcessor(tsp5);
+ //calculator->addCoProcessor(tsp1);
+ //calculator->addCoProcessor(tsp2);
+ //calculator->addCoProcessor(tsp3);
+ //calculator->addCoProcessor(tsp4);
+ //calculator->addCoProcessor(tsp5);
//calculator->addCoProcessor(tsp6);
- calculator->addCoProcessor(tav);
+ //calculator->addCoProcessor(tav);
if (myid==0) UBLOG(logINFO, "Simulation-start");
diff --git a/source/Applications/Thermoplast/config.txt b/source/Applications/Thermoplast/config.txt
index 79fc3464f..954ed6821 100644
--- a/source/Applications/Thermoplast/config.txt
+++ b/source/Applications/Thermoplast/config.txt
@@ -8,7 +8,7 @@ boundingBox = 60 1370 10 190 1530 320 #test bb
blocknx = 10 10 10
#blocknx = 300 420 320
endTime = 100000
-outTime = 1000
+outTime = 500
availMem = 25e9
uLB = 0.1
diff --git a/source/Applications/Thermoplast/thermoplast.cpp b/source/Applications/Thermoplast/thermoplast.cpp
index 3ae47c737..fac556f3c 100644
--- a/source/Applications/Thermoplast/thermoplast.cpp
+++ b/source/Applications/Thermoplast/thermoplast.cpp
@@ -297,8 +297,8 @@ void thermoplast(string configname)
SPtr<Interactor3D> p2Int = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(p2Geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
//////////////////////////////////////////////////////////////////////////
- //SPtr<Grid3DVisitor> peVisitor(new PePartitioningGridVisitor(comm, demCoProcessor));
- SPtr<Grid3DVisitor> peVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::KWAY));
+ 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);
@@ -438,6 +438,7 @@ void thermoplast(string configname)
SPtr<WriteBoundaryConditionsCoProcessor> writeBCCoProcessor(new WriteBoundaryConditionsCoProcessor(grid, visSch, pathOut,
WbWriterVtkXmlBinary::getInstance(), comm));
+ writeBCCoProcessor->process(0);
SPtr<WriteDemObjectsCoProcessor> writeDemObjectsCoProcessor(new WriteDemObjectsCoProcessor(grid, visSch, pathOut, WbWriterVtkXmlBinary::getInstance(), demCoProcessor, comm));
writeDemObjectsCoProcessor->process(0);
@@ -454,7 +455,7 @@ void thermoplast(string configname)
calculator->addCoProcessor(createSphereCoProcessor);
calculator->addCoProcessor(demCoProcessor);
- //calculator->addCoProcessor(writeBCCoProcessor);
+ calculator->addCoProcessor(writeBCCoProcessor);
calculator->addCoProcessor(writeDemObjectsCoProcessor);
calculator->addCoProcessor(writeMQCoProcessor);
//calculator->addCoProcessor(restartDemObjectsCoProcessor);
diff --git a/source/DemCoupling/CreateDemObjectsCoProcessor.cpp b/source/DemCoupling/CreateDemObjectsCoProcessor.cpp
index 5c6a43e6c..146e61160 100644
--- a/source/DemCoupling/CreateDemObjectsCoProcessor.cpp
+++ b/source/DemCoupling/CreateDemObjectsCoProcessor.cpp
@@ -43,9 +43,10 @@ void CreateDemObjectsCoProcessor::process(double step)
if (scheduler->isDue(step))
{
int istep = static_cast<int>(step);
- if (comm->isRoot()) UBLOG(logINFO, "CreateDemObjectsCoProcessor::process start step: " << istep);
+
#ifdef TIMING
+ if (comm->isRoot()) UBLOG(logINFO, "CreateDemObjectsCoProcessor::process start step: " << istep);
timer.resetAndStart();
#endif
@@ -54,16 +55,16 @@ void CreateDemObjectsCoProcessor::process(double step)
#ifdef TIMING
if (comm->isRoot()) UBLOG(logINFO, "createGeoObjects() time = "<<timer.stop()<<" s");
if (comm->isRoot()) UBLOG(logINFO, "number of objects = "<<(int)(geoObjectPrototypeVector.size()));
+ if (comm->isRoot()) UBLOG(logINFO, "CreateDemObjectsCoProcessor::process stop step: " << istep);
#endif
-// demCoProcessor->distributeIDs();
+ //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);
+
}
}
//////////////////////////////////////////////////////////////////////////
@@ -107,14 +108,15 @@ void CreateDemObjectsCoProcessor::createGeoObjects()
//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);
+ //blockVector.clear();
+ //UbTupleInt3 blockNX=grid->getBlockNX();
+ //grid->getAllBlocksByCuboid(AABB[0]-(double)val<1>(blockNX)*2.0, AABB[1]-(double)val<2>(blockNX)*2.0, AABB[2]-(double)val<3>(blockNX)*2.0, AABB[3]+(double)val<1>(blockNX)*2.0, AABB[4]+(double)val<2>(blockNX)*2.0, AABB[5]+(double)val<3>(blockNX)*2.0, blockVector);
+ //for (std::shared_ptr<Block3D> block : blockVector)
+ // geoObjectInt->setBCBlock(block);
- //UBLOG(logINFO, "grid->accept(setBcVisitor) time = "<<timer.stop());
- geoObjectInt->initInteractor();
+ ////UBLOG(logINFO, "grid->accept(setBcVisitor) time = "<<timer.stop());
+ //geoObjectInt->initInteractor();
//UBLOG(logINFO, "geoObjectInt->initInteractor() time = "<<timer.stop());
demCoProcessor->addInteractor(geoObjectInt, demObjectMaterial, initalVelocity[i]);
//UBLOG(logINFO, "demCoProcessor->addInteractor() time = "<<timer.stop());
diff --git a/source/DemCoupling/CreateDemObjectsCoProcessor.h b/source/DemCoupling/CreateDemObjectsCoProcessor.h
index bfd0d0c64..ffb26a35d 100644
--- a/source/DemCoupling/CreateDemObjectsCoProcessor.h
+++ b/source/DemCoupling/CreateDemObjectsCoProcessor.h
@@ -7,7 +7,7 @@
#include <array>
-#define TIMING
+//#define TIMING
#ifdef TIMING
#include "UbTiming.h"
diff --git a/source/DemCoupling/DemCoProcessor.cpp b/source/DemCoupling/DemCoProcessor.cpp
index 8ac7a686f..b779534df 100644
--- a/source/DemCoupling/DemCoProcessor.cpp
+++ b/source/DemCoupling/DemCoProcessor.cpp
@@ -60,9 +60,28 @@ void DemCoProcessor::addInteractor(std::shared_ptr<MovableObjectInteractor> inte
interactors.push_back(interactor);
const int id = static_cast<int>(interactors.size()) - 1;
interactor->setID(id);
- const auto peGeometry = this->createPhysicsEngineGeometryAdapter(interactor, physicsEngineMaterial);
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(peGeometry)->isActive()) peGeometry->setLinearVelolocity(initalVelocity);
- physicsEngineGeometries.push_back(peGeometry);
+ 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();
+ }
+ physicsEngineGeometrieAdapters.push_back(peGeometryAdapter);
}
@@ -72,11 +91,11 @@ std::shared_ptr<PhysicsEngineGeometryAdapter> DemCoProcessor::createPhysicsEngin
SPtr<GbSphere3D> vfSphere = std::static_pointer_cast<GbSphere3D>(interactor->getGbObject3D());
const Vector3D position(vfSphere->getX1Centroid(), vfSphere->getX2Centroid(), vfSphere->getX3Centroid());
- auto peGeometry = this->physicsEngineSolver->createPhysicsEngineGeometryAdapter(id, position, vfSphere->getRadius(), physicsEngineMaterial);
+ auto peGeometryAdapter = this->physicsEngineSolver->createPhysicsEngineGeometryAdapter(id, position, vfSphere->getRadius(), physicsEngineMaterial);
//if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(peGeometry)->isActive())
//{
- interactor->setPhysicsEngineGeometry(peGeometry);
- return peGeometry;
+ interactor->setPhysicsEngineGeometry(peGeometryAdapter);
+ return peGeometryAdapter;
//}
//else
//{
@@ -157,11 +176,11 @@ std::shared_ptr<PhysicsEngineSolverAdapter> DemCoProcessor::getPhysicsEngineSolv
void DemCoProcessor::applyForcesOnGeometries()
{
- for (int i = 0; i < physicsEngineGeometries.size(); i++)
+ for (int i = 0; i < physicsEngineGeometrieAdapters.size(); i++)
{
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
+ if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive())
{
- this->setForcesToObject(grid, interactors[i], physicsEngineGeometries[i]);
+ this->setForcesToObject(grid, interactors[i], physicsEngineGeometrieAdapters[i]);
//physicsEngineGeometries[i]->setLinearVelolocity(Vector3D(-0.001, 0.0, 0.0));
//physicsEngineGeometries[i]->setAngularVelocity(Vector3D(0.01, 0.01, 0.01));
@@ -206,17 +225,17 @@ void DemCoProcessor::setForcesToObject(SPtr<Grid3D> grid, SPtr<MovableObjectInte
void DemCoProcessor::scaleForcesAndTorques(double scalingFactor)
{
- for (int i = 0; i < physicsEngineGeometries.size(); i++)
+ for (int i = 0; i < physicsEngineGeometrieAdapters.size(); i++)
{
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
+ if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive())
{
- const Vector3D force = physicsEngineGeometries[i]->getForce() * scalingFactor;
- const Vector3D torque = physicsEngineGeometries[i]->getTorque() * scalingFactor;
+ const Vector3D force = physicsEngineGeometrieAdapters[i]->getForce() * scalingFactor;
+ const Vector3D torque = physicsEngineGeometrieAdapters[i]->getTorque() * scalingFactor;
- physicsEngineGeometries[i]->resetForceAndTorque();
+ physicsEngineGeometrieAdapters[i]->resetForceAndTorque();
- physicsEngineGeometries[i]->setForce(force);
- physicsEngineGeometries[i]->setTorque(torque);
+ physicsEngineGeometrieAdapters[i]->setForce(force);
+ physicsEngineGeometrieAdapters[i]->setTorque(torque);
//UBLOG(logINFO, "F: (x,y,z) " << force);
//UBLOG(logINFO, "T: (x,y,z) " << torque);
@@ -250,9 +269,10 @@ void DemCoProcessor::moveVfGeoObject()
{
for (int i = 0; i < interactors.size(); i++)
{
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
+ if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive())
{
- interactors[i]->moveGbObjectTo(physicsEngineGeometries[i]->getPosition());
+ interactors[i]->moveGbObjectTo(physicsEngineGeometrieAdapters[i]->getPosition());
+ //UBLOG(logINFO, "DemCoProcessor::moveVfGeoObject() id = "<<std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->getId()<<" position="<<physicsEngineGeometrieAdapters[i]->getPosition()<<" rank="<<comm->getProcessID());
}
}
}
@@ -262,7 +282,7 @@ bool DemCoProcessor::isDemObjectInAABB(std::array<double, 6> AABB)
bool result = false;
for (int i = 0; i < interactors.size(); i++)
{
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
+ if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive())
{
SPtr<GbObject3D> geoObject = interactors[i]->getGbObject3D();
std::array <double, 2> minMax1;
@@ -307,7 +327,7 @@ void DemCoProcessor::addSurfaceTriangleSet(std::vector<UbTupleFloat3>& nodes, st
for (int i = 0; i < interactors.size(); i++)
{
//UBLOG(logINFO, "DemCoProcessor::addSurfaceTriangleSet()1");
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
+ if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive())
{
//UBLOG(logINFO, "DemCoProcessor::addSurfaceTriangleSet()2");
interactors[i]->getGbObject3D()->addSurfaceTriangleSet(nodes, triangles);
@@ -319,12 +339,12 @@ void DemCoProcessor::getObjectsPropertiesVector(std::vector<double>& p)
{
for (int i = 0; i < interactors.size(); i++)
{
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
+ if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive())
{
p.push_back(interactors[i]->getGbObject3D()->getX1Centroid());
p.push_back(interactors[i]->getGbObject3D()->getX2Centroid());
p.push_back(interactors[i]->getGbObject3D()->getX3Centroid());
- Vector3D v = physicsEngineGeometries[i]->getLinearVelocity();
+ Vector3D v = physicsEngineGeometrieAdapters[i]->getLinearVelocity();
p.push_back(v[0]);
p.push_back(v[1]);
p.push_back(v[2]);
@@ -346,9 +366,9 @@ void DemCoProcessor::addPeGeo(walberla::pe::RigidBody * peGeo)
// }
//}
- if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometries.size()) return;
+ if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
- auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[peGeo->getID()]);
+ auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
if (geometry->getId() == peGeo->getID())
{
geometry->setActive();
@@ -373,9 +393,9 @@ void DemCoProcessor::removePeGeo(walberla::pe::RigidBody * peGeo)
// }
//}
- if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometries.size()) return;
+ if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
- auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[peGeo->getID()]);
+ auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
if (geometry->getId() == peGeo->getID())
{
geometry->setInactive();
@@ -391,7 +411,7 @@ bool DemCoProcessor::isSpheresIntersection(double centerX1, double centerX2, dou
bool result = false;
for (int i = 0; i < interactors.size(); i++)
{
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
+ 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 );
@@ -415,55 +435,55 @@ bool DemCoProcessor::isSpheresIntersection(double centerX1, double centerX2, dou
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::distributeIDs()
+{
+ std::vector<int> peIDsSend;
+ std::vector<int> vfIDsSend;
+
+ for (int i = 0; i < interactors.size(); i++)
+ {
+ if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive())
+ {
+ peIDsSend.push_back(std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[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>(physicsEngineGeometrieAdapters[i])->setId(it->second);
+ //std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->setId(idMap.find(interactors[i]->getID())->second);
+ }
+
+ 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]);
+ }
+ }
+}
void DemCoProcessor::setBlockVisitor(std::shared_ptr<BoundaryConditionsBlockVisitor> boundaryConditionsBlockVisitor)
{
diff --git a/source/DemCoupling/DemCoProcessor.h b/source/DemCoupling/DemCoProcessor.h
index dcc8074f3..9071084e3 100644
--- a/source/DemCoupling/DemCoProcessor.h
+++ b/source/DemCoupling/DemCoProcessor.h
@@ -45,7 +45,7 @@ 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);
@@ -66,7 +66,7 @@ private:
std::vector<std::shared_ptr<MovableObjectInteractor> > interactors;
std::shared_ptr<ForceCalculator> forceCalculator;
std::shared_ptr<PhysicsEngineSolverAdapter> physicsEngineSolver;
- std::vector<std::shared_ptr<PhysicsEngineGeometryAdapter> > physicsEngineGeometries;
+ std::vector<std::shared_ptr<PhysicsEngineGeometryAdapter> > physicsEngineGeometrieAdapters;
double intermediateDemSteps;
SPtr<BoundaryConditionsBlockVisitor> boundaryConditionsBlockVisitor;
diff --git a/source/DemCoupling/MovableObjectInteractor.cpp b/source/DemCoupling/MovableObjectInteractor.cpp
index bf658073d..c6a2b6933 100644
--- a/source/DemCoupling/MovableObjectInteractor.cpp
+++ b/source/DemCoupling/MovableObjectInteractor.cpp
@@ -172,10 +172,17 @@ void MovableObjectInteractor::setBcBlocks()
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);
+ 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)
- setBCBlock(block);
+ {
+ if (block->getKernel())
+ {
+ setBCBlock(block);
+ }
+ }
}
void MovableObjectInteractor::updateVelocityBc()
diff --git a/source/VirtualFluidsCore/Grid/Grid3D.cpp b/source/VirtualFluidsCore/Grid/Grid3D.cpp
index 88e6f9d8a..fa991fd84 100644
--- a/source/VirtualFluidsCore/Grid/Grid3D.cpp
+++ b/source/VirtualFluidsCore/Grid/Grid3D.cpp
@@ -1831,6 +1831,61 @@ void Grid3D::getBlocksByCuboid(int level, double minX1, double minX2, double min
std::copy(blockset.begin(), blockset.end(), blocks.begin());
}
//////////////////////////////////////////////////////////////////////////
+void Grid3D::getAllBlocksByCuboid(double minX1, double minX2, double minX3, double maxX1, double maxX2, double maxX3, std::vector<SPtr<Block3D>>& blocks)
+{
+ int coarsestLevel = this->getCoarsestInitializedLevel();
+ int finestLevel = this->getFinestInitializedLevel();
+
+ //////////////////////////////////////////////////////////////////////////
+ //MINIMALE BLOCK-INDIZES BESTIMMEN
+ //
+ //min:
+ double dMinX1 = trafo->transformForwardToX1Coordinate(minX1, minX2, minX3)*(1<<finestLevel);
+ double dMinX2 = trafo->transformForwardToX2Coordinate(minX1, minX2, minX3)*(1<<finestLevel);
+ double dMinX3 = trafo->transformForwardToX3Coordinate(minX1, minX2, minX3)*(1<<finestLevel);
+
+ //Achtung, wenn minX1 genau auf grenze zwischen zwei bloecken -> der "kleinere" muss genommen werden,
+ //da beim Transformieren der "groessere" Index rauskommt
+ int iMinX1 = (int)dMinX1; if (UbMath::zero(dMinX1-iMinX1)) iMinX1-=1;
+ int iMinX2 = (int)dMinX2; if (UbMath::zero(dMinX2-iMinX2)) iMinX2-=1;
+ int iMinX3 = (int)dMinX3; if (UbMath::zero(dMinX3-iMinX3)) iMinX3-=1;
+
+ //max (hier kann die Zusatzabfrage vernachlaessigt werden):
+ int iMaxX1 = (int)(trafo->transformForwardToX1Coordinate(maxX1, maxX2, maxX3)*(1<<finestLevel));
+ int iMaxX2 = (int)(trafo->transformForwardToX2Coordinate(maxX1, maxX2, maxX3)*(1<<finestLevel));
+ int iMaxX3 = (int)(trafo->transformForwardToX3Coordinate(maxX1, maxX2, maxX3)*(1<<finestLevel));
+
+ SPtr<Block3D> block;
+
+ //set, um doppelte bloecke zu vermeiden, die u.U. bei periodic auftreten koennen
+ std::set<SPtr<Block3D>> blockset;
+ for (int level=coarsestLevel; level<=finestLevel; level++)
+ {
+ //damit bei negativen werten auch der "kleinere" genommen wird -> floor!
+ int minx1 = (int)std::floor((double)iMinX1/(1<<(finestLevel-level)));
+ int minx2 = (int)std::floor((double)iMinX2/(1<<(finestLevel-level)));
+ int minx3 = (int)std::floor((double)iMinX3/(1<<(finestLevel-level)));
+
+ int maxx1 = iMaxX1/(1<<(finestLevel-level));
+ int maxx2 = iMaxX2/(1<<(finestLevel-level));
+ int maxx3 = iMaxX3/(1<<(finestLevel-level));
+
+ for (int ix1=minx1; ix1<=maxx1; ix1++)
+ for (int ix2=minx2; ix2<=maxx2; ix2++)
+ for (int ix3=minx3; ix3<=maxx3; ix3++)
+ if ((block=this->getBlock(ix1, ix2, ix3, level)))
+ {
+ if (block)
+ {
+ blockset.insert(block);
+ }
+ }
+ }
+
+ blocks.resize(blockset.size());
+ std::copy(blockset.begin(), blockset.end(), blocks.begin());
+}
+//////////////////////////////////////////////////////////////////////////
void Grid3D::calcStartCoordinatesAndDelta(SPtr<Block3D> block, double& worldX1, double& worldX2, double& worldX3, double& deltaX)
{
int blocklevel = block->getLevel();
diff --git a/source/VirtualFluidsCore/Grid/Grid3D.h b/source/VirtualFluidsCore/Grid/Grid3D.h
index bf58fcabc..35188bc0a 100644
--- a/source/VirtualFluidsCore/Grid/Grid3D.h
+++ b/source/VirtualFluidsCore/Grid/Grid3D.h
@@ -52,6 +52,7 @@ public:
void getBlocksByCuboid(int level, double minX1, double minX2, double minX3,
double maxX1, double maxX2, double maxX3,
std::vector<SPtr<Block3D>>& blocks);
+ void getAllBlocksByCuboid(double minX1, double minX2, double minX3, double maxX1, double maxX2, double maxX3, std::vector<SPtr<Block3D>>& blocks);
//!get blocks for level
void getBlocks(int level, std::vector<SPtr<Block3D>>& blockVector);
//!get blocks for level with current rank
--
GitLab