From cbe233f3d64feb7c1ea616033ab4e6dd78ced9ce Mon Sep 17 00:00:00 2001 From: Konstantin Kutscher <kutscher@irmb.tu-bs.de> Date: Fri, 3 Aug 2018 13:03:19 +0200 Subject: [PATCH] fixed update for peGeoObjects pointers, exchange the itaration of forest with callback function --- .../Applications/DLR-F16-Solid/f16-solid.cfg | 8 +-- source/Applications/DLR-F16-Solid/f16.cpp | 23 +++++---- source/Applications/Thermoplast/config.txt | 8 +-- .../Applications/Thermoplast/thermoplast.cpp | 38 +++++++------- source/DemCoupling/DemCoProcessor.cpp | 51 ++++++++++++++++++- source/DemCoupling/DemCoProcessor.h | 6 ++- .../pe/PePhysicsEngineSolverAdapter.cpp | 13 +++++ .../pe/PePhysicsEngineSolverAdapter.h | 2 + 8 files changed, 110 insertions(+), 39 deletions(-) diff --git a/source/Applications/DLR-F16-Solid/f16-solid.cfg b/source/Applications/DLR-F16-Solid/f16-solid.cfg index fab078fe8..fb55f3e85 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-mic +pathOut = d:/temp/DLR-F16-Solid-mic7 pathGeo = d:/Projects/SFB880/DLR-F16/Geometry fngFileWhole1 = F16_broad_Quad_noTape_full.stl @@ -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 @@ -44,7 +44,7 @@ cpStep = 1000 cpStart = 3000 outTimeStep = 100 -outTimeStart = 3000 +outTimeStart = 100 endTime = 10000 diff --git a/source/Applications/DLR-F16-Solid/f16.cpp b/source/Applications/DLR-F16-Solid/f16.cpp index 763ff4f70..96f6143a6 100644 --- a/source/Applications/DLR-F16-Solid/f16.cpp +++ b/source/Applications/DLR-F16-Solid/f16.cpp @@ -361,7 +361,7 @@ void run(string configname) /////delete solid blocks if (myid==0) UBLOG(logINFO, "deleteSolidBlocks - start"); - SetSolidOrBoundaryBlockVisitor v(fngIntrWhole1, SetSolidOrBoundaryBlockVisitor::SOLID); + SetSolidBlocksBlockVisitor v(fngIntrWhole1); grid->accept(v); std::vector<SPtr<Block3D>>& sb = fngIntrWhole1->getSolidBlockSet(); for (SPtr<Block3D> block : sb) @@ -487,6 +487,11 @@ void run(string configname) ppblocks.process(2); } + GbCuboid3DPtr mic6(new GbCuboid3D( 0.3, 0.015, -0.46+4.25*deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.46+5.25*deltaXcoarse)); + if (myid==0) GbSystem3D::writeGeoObject(mic6.get(), pathOut+"/geo/mic6", WbWriterVtkXmlBinary::getInstance()); + GbCuboid3DPtr mic7(new GbCuboid3D(0.3, 0.015, -0.3+4.25*deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.3+5.25*deltaXcoarse)); + if (myid==0) GbSystem3D::writeGeoObject(mic7.get(), pathOut+"/geo/mic7", WbWriterVtkXmlBinary::getInstance()); + unsigned long long numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks(); int ghostLayer = 3; unsigned long long numberOfNodesPerBlock = (unsigned long long)(blockNx[0])* (unsigned long long)(blockNx[1])* (unsigned long long)(blockNx[2]); @@ -653,7 +658,7 @@ void run(string configname) fngMeshWhole2->rotate(0.0, 0.5, 0.0); if (myid==0) GbSystem3D::writeGeoObject(fngMeshWhole2.get(), pathOut+"/geo/fngMeshWhole3", WbWriterVtkXmlBinary::getInstance()); SPtr<Interactor3D> fngIntrWhole2 = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(fngMeshWhole2, grid, noSlipBCAdapter, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy)); - SetSolidOrBoundaryBlockVisitor v(fngIntrWhole2, SetSolidOrBoundaryBlockVisitor::BC); + SetBcBlocksBlockVisitor v(fngIntrWhole2); grid->accept(v); fngIntrWhole2->initInteractor(); /////////////////////////////////////////////////////////////////////////////////////////////// @@ -727,13 +732,13 @@ void run(string configname) if (myid==0) GbSystem3D::writeGeoObject(mic5->getBoundingBox().get(), pathOut+"/geo/mic5", WbWriterVtkXmlBinary::getInstance()); SPtr<TimeseriesCoProcessor> tsp5(new TimeseriesCoProcessor(grid, stepMV, mic5, pathOut+"/mic/mic5", comm)); - SPtr<IntegrateValuesHelper> mic6(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.46+4.25*deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.46+5.25*deltaXcoarse)); - if (myid==0) GbSystem3D::writeGeoObject(mic6->getBoundingBox().get(), pathOut+"/geo/mic6", WbWriterVtkXmlBinary::getInstance()); - SPtr<TimeseriesCoProcessor> tsp6(new TimeseriesCoProcessor(grid, stepMV, mic6, pathOut+"/mic/mic6", comm)); + //SPtr<IntegrateValuesHelper> mic6(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.46+4.25*deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.46+5.25*deltaXcoarse)); + //if (myid==0) GbSystem3D::writeGeoObject(mic6->getBoundingBox().get(), pathOut+"/geo/mic6", WbWriterVtkXmlBinary::getInstance()); + //SPtr<TimeseriesCoProcessor> tsp6(new TimeseriesCoProcessor(grid, stepMV, mic6, pathOut+"/mic/mic6", comm)); - SPtr<IntegrateValuesHelper> mic7(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.3+4.25*deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.3+5.25*deltaXcoarse)); - if (myid==0) GbSystem3D::writeGeoObject(mic7->getBoundingBox().get(), pathOut+"/geo/mic7", WbWriterVtkXmlBinary::getInstance()); - SPtr<TimeseriesCoProcessor> tsp7(new TimeseriesCoProcessor(grid, stepMV, mic7, pathOut+"/mic/mic7", comm)); + //SPtr<IntegrateValuesHelper> mic7(new IntegrateValuesHelper(grid, comm, 0.3, 0.015, -0.3+4.25*deltaXcoarse, 0.3+deltaXcoarse, 0.015+deltaXcoarse, -0.3+5.25*deltaXcoarse)); + //if (myid==0) GbSystem3D::writeGeoObject(mic7->getBoundingBox().get(), pathOut+"/geo/mic7", WbWriterVtkXmlBinary::getInstance()); + //SPtr<TimeseriesCoProcessor> tsp7(new TimeseriesCoProcessor(grid, stepMV, mic7, pathOut+"/mic/mic7", comm)); omp_set_num_threads(numOfThreads); SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1)); @@ -747,7 +752,7 @@ void run(string configname) calculator->addCoProcessor(tsp3); calculator->addCoProcessor(tsp4); calculator->addCoProcessor(tsp5); - calculator->addCoProcessor(tsp6); + //calculator->addCoProcessor(tsp6); calculator->addCoProcessor(tav); diff --git a/source/Applications/Thermoplast/config.txt b/source/Applications/Thermoplast/config.txt index 7c53a6839..b25b5012c 100644 --- a/source/Applications/Thermoplast/config.txt +++ b/source/Applications/Thermoplast/config.txt @@ -1,21 +1,21 @@ #simulation parameters #boundingBox = 0 0 0 300 1520 2320 -boundingBox = 60 1370 10 190 1530 320 #test bb +boundingBox = 60 1370 150 190 1530 320 #test bb #boundingBox = 60 0 10 190 1530 2320 #production bb blocknx = 10 10 10 #blocknx = 300 420 320 -endTime = 100 -outTime = 100 +endTime = 1000 +outTime = 10000 availMem = 25e9 uLB = 0.1 #PE parameters peMinOffset = 46 2 2 peMaxOffset = -8 -25 -2 -sphereTime = 10 +sphereTime = 100000 #geometry files pathGeo = d:/Projects/ThermoPlast/SimGeo diff --git a/source/Applications/Thermoplast/thermoplast.cpp b/source/Applications/Thermoplast/thermoplast.cpp index 5009794af..9ca15b763 100644 --- a/source/Applications/Thermoplast/thermoplast.cpp +++ b/source/Applications/Thermoplast/thermoplast.cpp @@ -223,7 +223,7 @@ 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)); + GbCuboid3DPtr geoInflow1(new GbCuboid3D(g_minX1-blockLength, g_maxX2-120.0, g_minX3+190.0-150.0, g_minX1+1, g_maxX2+20.0, g_minX3+320.0)); if (myid == 0) GbSystem3D::writeGeoObject(geoInflow1.get(), pathOut + "/geo/geoInflow1", WbWriterVtkXmlASCII::getInstance()); if (!restart) @@ -381,28 +381,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 < 6; x3++) - // for (int x2 = 0; x2 < 4; 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*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); diff --git a/source/DemCoupling/DemCoProcessor.cpp b/source/DemCoupling/DemCoProcessor.cpp index b7bd63185..5634027a9 100644 --- a/source/DemCoupling/DemCoProcessor.cpp +++ b/source/DemCoupling/DemCoProcessor.cpp @@ -14,6 +14,7 @@ #include "PhysicsEngineMaterialAdapter.h" #include "PhysicsEngineGeometryAdapter.h" #include "PhysicsEngineSolverAdapter.h" +#include "PePhysicsEngineSolverAdapter.h" #include "PePhysicsEngineGeometryAdapter.h" #include "BoundaryConditions.h" @@ -26,6 +27,7 @@ #include <array> +#include <functional> DemCoProcessor::DemCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, SPtr<Communicator> comm, std::shared_ptr<ForceCalculator> forceCalculator, std::shared_ptr<PhysicsEngineSolverAdapter> physicsEngineSolver, double intermediatePeSteps) : CoProcessor(grid, s), comm(comm), forceCalculator(forceCalculator), physicsEngineSolver(physicsEngineSolver), intermediateDemSteps(intermediatePeSteps) @@ -33,6 +35,19 @@ DemCoProcessor::DemCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, SPtr<Comm #ifdef TIMING timer.resetAndStart(); #endif + + 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) + { + walberla::pe::Storage* storage = blockIt->getData< walberla::pe::Storage >(*storageId.get()); + walberla::pe::BodyStorage* bodyStorage = &(*storage)[0]; + + bodyStorage->registerAddCallback("DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::addPeGeo), this)); + bodyStorage->registerRemoveCallback("DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::removePeGeo), this)); + } } DemCoProcessor::~DemCoProcessor() @@ -228,7 +243,7 @@ void DemCoProcessor::calculateDemTimeStep(double step) for (int i = 0; i < physicsEngineGeometries.size(); i++) { - physicsEngineSolver->updateGeometry(physicsEngineGeometries[i]); + //physicsEngineSolver->updateGeometry(physicsEngineGeometries[i]); if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive()) { interactors[i]->setPhysicsEngineGeometry(physicsEngineGeometries[i]); @@ -300,8 +315,10 @@ 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()) { + //UBLOG(logINFO, "DemCoProcessor::addSurfaceTriangleSet()2"); interactors[i]->getGbObject3D()->addSurfaceTriangleSet(nodes, triangles); } } @@ -324,6 +341,36 @@ 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++) + { + auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i]); + if (geometry->getId() == peGeo->getID()) + { + geometry->setActive(); + geometry->setGeometry(peGeo); + return; + } + } +} + +void DemCoProcessor::removePeGeo(walberla::pe::RigidBody * peGeo) +{ + //UBLOG(logINFO, "DemCoProcessor::removePeGeo()"); + for (int i = 0; i < physicsEngineGeometries.size(); i++) + { + auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i]); + if (geometry->getId() == peGeo->getID()) + { + geometry->setInactive(); + geometry->setGeometry(NULL); + return; + } + } +} + void DemCoProcessor::distributeIDs() { std::vector<int> peIDsSend; @@ -366,7 +413,7 @@ void DemCoProcessor::distributeIDs() for (int i = 0; i < physicsEngineGeometries.size(); i++) { - physicsEngineSolver->updateGeometry(physicsEngineGeometries[i]); + //physicsEngineSolver->updateGeometry(physicsEngineGeometries[i]); if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive()) { interactors[i]->setPhysicsEngineGeometry(physicsEngineGeometries[i]); diff --git a/source/DemCoupling/DemCoProcessor.h b/source/DemCoupling/DemCoProcessor.h index 898699eb6..f85e0e7d1 100644 --- a/source/DemCoupling/DemCoProcessor.h +++ b/source/DemCoupling/DemCoProcessor.h @@ -14,7 +14,9 @@ #include "CoProcessor.h" #include "UbTuple.h" -#define TIMING +#include <pe/basic.h> + +//#define TIMING #ifdef TIMING #include "UbTiming.h" @@ -48,6 +50,8 @@ public: 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); private: std::shared_ptr<PhysicsEngineGeometryAdapter> createPhysicsEngineGeometryAdapter(std::shared_ptr<MovableObjectInteractor> interactor, std::shared_ptr<PhysicsEngineMaterialAdapter> physicsEngineMaterial) const; diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp index c85219e05..7fdc9cae0 100644 --- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp +++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp @@ -150,6 +150,9 @@ walberla::pe::RigidBody* PePhysicsEngineSolverAdapter::getPeGeoObject(walberla:: } } return NULL; + + //return getBody(*globalBodyStorage.get(), *forest.get(), *storageId.get(), id, walberla::pe::StorageSelect::LOCAL); + } void PePhysicsEngineSolverAdapter::updateGeometry(std::shared_ptr<PhysicsEngineGeometryAdapter> geometryAdapter) @@ -205,3 +208,13 @@ void PePhysicsEngineSolverAdapter::loadFromFile(const std::string & path) // ccd->reloadBodies(); //} } + +std::shared_ptr<walberla::blockforest::BlockForest> PePhysicsEngineSolverAdapter::getBlockForest() +{ + return forest; +} + +std::shared_ptr<walberla::domain_decomposition::BlockDataID> PePhysicsEngineSolverAdapter::getStorageId() +{ + return storageId; +} diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h index f79facea7..6dbdda185 100644 --- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h +++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h @@ -74,6 +74,8 @@ public: std::shared_ptr< walberla::blockforest::BlockForest > getForest(); void saveToFile(const std::string& path); void loadFromFile(const std::string& path); + std::shared_ptr<walberla::blockforest::BlockForest> getBlockForest(); + std::shared_ptr<walberla::domain_decomposition::BlockDataID> getStorageId(); private: void initalizePeEnvironment(); -- GitLab