fixed update for peGeoObjects pointers, exchange the itaration of forest with callback function

source/Applications/DLR-F16-Solid/f16-solid.cfg

-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
 

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);
 
 

source/Applications/Thermoplast/config.txt

 #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

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);
 

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]);

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;

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;
+}

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();