From bc1801eca4c5a7278bcef266459a3e261e0c0e48 Mon Sep 17 00:00:00 2001
From: Konstantin Kutscher <kutscher@irmb.tu-bs.de>
Date: Fri, 8 Jun 2018 17:34:37 +0200
Subject: [PATCH] CreateGeoObjectsCoProcessor works now correct

---
 .../Applications/Thermoplast/thermoplast.cpp  |  18 +--
 .../CreateGeoObjectsCoProcessor.cpp           | 107 +++++++++---------
 source/DemCoupling/DemCoProcessor.cpp         |   8 --
 .../WriteMacroscopicQuantitiesCoProcessor.cpp |   1 -
 4 files changed, 62 insertions(+), 72 deletions(-)

diff --git a/source/Applications/Thermoplast/thermoplast.cpp b/source/Applications/Thermoplast/thermoplast.cpp
index b881d3378..bd6fc4a61 100644
--- a/source/Applications/Thermoplast/thermoplast.cpp
+++ b/source/Applications/Thermoplast/thermoplast.cpp
@@ -67,11 +67,11 @@ void pf1()
    int myid = comm->getProcessID();
 
    //parameters
-   string          pathname = "d:/temp/thermoplast";
+   string          pathname = "d:/temp/thermoplast3";
    int             numOfThreads = 1;
    int             blocknx[3] ={ 16,16,16 };
    double          endTime = 1000000;
-   double          outTime = 10;
+   double          outTime = 100;
    double          availMem = 8e9;
    double          deltax = 1;
    double          rhoLB = 0.0;
@@ -109,6 +109,7 @@ void pf1()
       UBLOG(logINFO, "rhoLB = " << rhoLB);
       UBLOG(logINFO, "nuLB = " << nuLB);
       UBLOG(logINFO, "deltaX = " << deltax);
+      UBLOG(logINFO, "Re = " << uLB*(g_maxX3-g_minX3)/nuLB);
       UBLOG(logINFO, "Preprocess - start");
    }
 
@@ -185,7 +186,7 @@ void pf1()
    //const std::shared_ptr<LBMUnitConverter> lbmUnitConverter = std::make_shared<LBMUnitConverter>(refLengthWorld, LBMUnitConverter::WORLD_MATERIAL::AIR_20C, refLengthLb);
    const SPtr<LBMUnitConverter> lbmUnitConverter(new LBMUnitConverter(refLengthWorld, LBMUnitConverter::WORLD_MATERIAL::OIL, refLengthLb));
    if (myid == 0) std::cout << lbmUnitConverter->toString() << std::endl;
-   double rhoSphere = 915 * lbmUnitConverter->getFactorDensityWToLb(); // 2 // kg/m^3
+   double rhoSphere = 0.7625;//915 * lbmUnitConverter->getFactorDensityWToLb(); // 2 // kg/m^3
    if (myid == 0) UBLOG(logINFO, "rhoSphere = "<<rhoSphere);
    SPtr<PhysicsEngineMaterialAdapter> sphereMaterial(new PePhysicsEngineMaterialAdapter("Polyethylen", rhoSphere, 0, 0.15, 0.1, 0.45, 0.5, 1, 0, 0));
    const int timestep = 2;
@@ -272,12 +273,15 @@ void pf1()
       SPtr<UbScheduler> geoSch(new UbScheduler(1));
       WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), comm);
       ppgeo.process(0);
+
+      WriteMacroscopicQuantitiesCoProcessor ppInit(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm);
+      ppInit.process(0);
    }
    
    if (myid == 0) UBLOG(logINFO, "Preprocess - end");
 
    //write data for visualization of macroscopic quantities
-   SPtr<UbScheduler> visSch(new UbScheduler(outTime,1));
+   SPtr<UbScheduler> visSch(new UbScheduler(outTime));
    SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, pathname, 
       WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
 
@@ -291,14 +295,14 @@ void pf1()
 
    //////////////////////////////////////////////////////////////////////////
 
-   SPtr<UbScheduler> sphereScheduler(new UbScheduler(1,1,2));
+   SPtr<UbScheduler> sphereScheduler(new UbScheduler(1000,10,1010));
    SPtr<CreateGeoObjectsCoProcessor> createSphereCoProcessor(new CreateGeoObjectsCoProcessor(grid,sphereScheduler,demCoProcessor,sphere1,sphereMaterial));
 
 
    //start simulation 
    //omp_set_num_threads(numOfThreads);
-   SPtr<UbScheduler> stepGhostLayer(new UbScheduler(outTime));
-   SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
+   //SPtr<UbScheduler> stepGhostLayer(new UbScheduler(outTime));
+   SPtr<Calculator> calculator(new BasicCalculator(grid, peScheduler, endTime));
    //SPtr<Calculator> calculator(new DemCalculator(grid, stepGhostLayer, endTime));
    calculator->addCoProcessor(npr);
    
diff --git a/source/DemCoupling/CreateGeoObjectsCoProcessor.cpp b/source/DemCoupling/CreateGeoObjectsCoProcessor.cpp
index 911debdb1..c1eab16c0 100644
--- a/source/DemCoupling/CreateGeoObjectsCoProcessor.cpp
+++ b/source/DemCoupling/CreateGeoObjectsCoProcessor.cpp
@@ -25,7 +25,31 @@ CreateGeoObjectsCoProcessor::CreateGeoObjectsCoProcessor(SPtr<Grid3D> grid, SPtr
 //////////////////////////////////////////////////////////////////////////
 void CreateGeoObjectsCoProcessor::process(double step)
 {
-   //if (scheduler->isDue(step))
+   if (scheduler->isDue(step))
+   {
+      mu::Parser fct2;
+      fct2.SetExpr("U");
+      fct2.DefineConst("U", 0.0);
+      SPtr<BCAdapter> velocityBcParticleAdapter(new VelocityBCAdapter(true, false, false, fct2, 0, BCFunction::INFCONST));
+      velocityBcParticleAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
+
+      SPtr<MovableObjectInteractor> geoObjectInt;
+      const std::shared_ptr<Reconstructor> velocityReconstructor(new VelocityBcReconstructor());
+      const std::shared_ptr<Reconstructor> extrapolationReconstructor(new ExtrapolationReconstructor(velocityReconstructor));
+      
+      SPtr<GbObject3D> geoObject((GbObject3D*)(geoObjectPrototype->clone()));
+      
+      geoObjectInt = SPtr<MovableObjectInteractor>(new MovableObjectInteractor(geoObject, grid, velocityBcParticleAdapter, Interactor3D::SOLID, extrapolationReconstructor, State::UNPIN));
+      //SetSolidOrBoundaryBlockVisitor setSolidVisitor(geoObjectInt, SetSolidOrBoundaryBlockVisitor::SOLID);
+      //grid->accept(setSolidVisitor);
+      SetSolidOrBoundaryBlockVisitor setBcVisitor(geoObjectInt, SetSolidOrBoundaryBlockVisitor::BC);
+      grid->accept(setBcVisitor);
+      geoObjectInt->initInteractor();
+      demCoProcessor->addInteractor(geoObjectInt, geoObjectMaterial, Vector3D(0.0, 0.0, 0.0));
+      demCoProcessor->distributeIDs();
+   }
+
+   //if (step == 1)
    //{
    //   mu::Parser fct2;
    //   fct2.SetExpr("U");
@@ -36,11 +60,11 @@ void CreateGeoObjectsCoProcessor::process(double step)
    //   SPtr<MovableObjectInteractor> geoObjectInt;
    //   const std::shared_ptr<Reconstructor> velocityReconstructor(new VelocityBcReconstructor());
    //   const std::shared_ptr<Reconstructor> extrapolationReconstructor(new ExtrapolationReconstructor(velocityReconstructor));
-   //   
+
    //   //SPtr<GbObject3D> geoObject((GbObject3D*)(geoObjectPrototype->clone()));
-   //   
+
    //   double radius = 5;
-   //   SPtr<GbObject3D> geoObject(new GbSphere3D(10, 16, 16, radius));
+   //   SPtr<GbObject3D> geoObject(new GbSphere3D(10, 16, 10, radius));
 
    //   geoObjectInt = SPtr<MovableObjectInteractor>(new MovableObjectInteractor(geoObject, grid, velocityBcParticleAdapter, Interactor3D::SOLID, extrapolationReconstructor, State::UNPIN));
    //   //SetSolidOrBoundaryBlockVisitor setSolidVisitor(geoObjectInt, SetSolidOrBoundaryBlockVisitor::SOLID);
@@ -54,62 +78,33 @@ void CreateGeoObjectsCoProcessor::process(double step)
    //   UBLOG(logINFO, "CreateGeoObjectsCoProcessor::process() step = "<<step);
    //}
 
-   if (step == 1)
-   {
-      mu::Parser fct2;
-      fct2.SetExpr("U");
-      fct2.DefineConst("U", 0.0);
-      SPtr<BCAdapter> velocityBcParticleAdapter(new VelocityBCAdapter(true, false, false, fct2, 0, BCFunction::INFCONST));
-      velocityBcParticleAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
-
-      SPtr<MovableObjectInteractor> geoObjectInt;
-      const std::shared_ptr<Reconstructor> velocityReconstructor(new VelocityBcReconstructor());
-      const std::shared_ptr<Reconstructor> extrapolationReconstructor(new ExtrapolationReconstructor(velocityReconstructor));
-
-      //SPtr<GbObject3D> geoObject((GbObject3D*)(geoObjectPrototype->clone()));
-
-      double radius = 5;
-      SPtr<GbObject3D> geoObject(new GbSphere3D(10, 16, 10, radius));
-
-      geoObjectInt = SPtr<MovableObjectInteractor>(new MovableObjectInteractor(geoObject, grid, velocityBcParticleAdapter, Interactor3D::SOLID, extrapolationReconstructor, State::UNPIN));
-      //SetSolidOrBoundaryBlockVisitor setSolidVisitor(geoObjectInt, SetSolidOrBoundaryBlockVisitor::SOLID);
-      //grid->accept(setSolidVisitor);
-      SetSolidOrBoundaryBlockVisitor setBcVisitor(geoObjectInt, SetSolidOrBoundaryBlockVisitor::BC);
-      grid->accept(setBcVisitor);
-      geoObjectInt->initInteractor();
-      demCoProcessor->addInteractor(geoObjectInt, geoObjectMaterial, Vector3D(0.0, 0.0, 0.0));
-      demCoProcessor->distributeIDs();
-
-      UBLOG(logINFO, "CreateGeoObjectsCoProcessor::process() step = "<<step);
-   }
-
-   if (step==1)
-   {
-      mu::Parser fct2;
-      fct2.SetExpr("U");
-      fct2.DefineConst("U", 0.0);
-      SPtr<BCAdapter> velocityBcParticleAdapter(new VelocityBCAdapter(true, false, false, fct2, 0, BCFunction::INFCONST));
-      velocityBcParticleAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
+   //if (step==1)
+   //{
+   //   mu::Parser fct2;
+   //   fct2.SetExpr("U");
+   //   fct2.DefineConst("U", 0.0);
+   //   SPtr<BCAdapter> velocityBcParticleAdapter(new VelocityBCAdapter(true, false, false, fct2, 0, BCFunction::INFCONST));
+   //   velocityBcParticleAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
 
-      SPtr<MovableObjectInteractor> geoObjectInt;
-      const std::shared_ptr<Reconstructor> velocityReconstructor(new VelocityBcReconstructor());
-      const std::shared_ptr<Reconstructor> extrapolationReconstructor(new ExtrapolationReconstructor(velocityReconstructor));
+   //   SPtr<MovableObjectInteractor> geoObjectInt;
+   //   const std::shared_ptr<Reconstructor> velocityReconstructor(new VelocityBcReconstructor());
+   //   const std::shared_ptr<Reconstructor> extrapolationReconstructor(new ExtrapolationReconstructor(velocityReconstructor));
 
-      //SPtr<GbObject3D> geoObject((GbObject3D*)(geoObjectPrototype->clone()));
+   //   //SPtr<GbObject3D> geoObject((GbObject3D*)(geoObjectPrototype->clone()));
 
-      double radius = 5;
-      SPtr<GbObject3D> geoObject(new GbSphere3D(10, 16, 23, radius));
+   //   double radius = 5;
+   //   SPtr<GbObject3D> geoObject(new GbSphere3D(10, 16, 23, radius));
 
-      geoObjectInt = SPtr<MovableObjectInteractor>(new MovableObjectInteractor(geoObject, grid, velocityBcParticleAdapter, Interactor3D::SOLID, extrapolationReconstructor, State::UNPIN));
-      //SetSolidOrBoundaryBlockVisitor setSolidVisitor(geoObjectInt, SetSolidOrBoundaryBlockVisitor::SOLID);
-      //grid->accept(setSolidVisitor);
-      SetSolidOrBoundaryBlockVisitor setBcVisitor(geoObjectInt, SetSolidOrBoundaryBlockVisitor::BC);
-      grid->accept(setBcVisitor);
-      geoObjectInt->initInteractor();
-      demCoProcessor->addInteractor(geoObjectInt, geoObjectMaterial, Vector3D(0.0, 0.0, 0.0));
-      demCoProcessor->distributeIDs();
+   //   geoObjectInt = SPtr<MovableObjectInteractor>(new MovableObjectInteractor(geoObject, grid, velocityBcParticleAdapter, Interactor3D::SOLID, extrapolationReconstructor, State::UNPIN));
+   //   //SetSolidOrBoundaryBlockVisitor setSolidVisitor(geoObjectInt, SetSolidOrBoundaryBlockVisitor::SOLID);
+   //   //grid->accept(setSolidVisitor);
+   //   SetSolidOrBoundaryBlockVisitor setBcVisitor(geoObjectInt, SetSolidOrBoundaryBlockVisitor::BC);
+   //   grid->accept(setBcVisitor);
+   //   geoObjectInt->initInteractor();
+   //   demCoProcessor->addInteractor(geoObjectInt, geoObjectMaterial, Vector3D(0.0, 0.0, 0.0));
+   //   demCoProcessor->distributeIDs();
 
-      UBLOG(logINFO, "CreateGeoObjectsCoProcessor::process() step = "<<step);
-   }
+   //   UBLOG(logINFO, "CreateGeoObjectsCoProcessor::process() step = "<<step);
+   //}
 }
 
diff --git a/source/DemCoupling/DemCoProcessor.cpp b/source/DemCoupling/DemCoProcessor.cpp
index e3ca8b87a..4a6ac2c6c 100644
--- a/source/DemCoupling/DemCoProcessor.cpp
+++ b/source/DemCoupling/DemCoProcessor.cpp
@@ -43,12 +43,10 @@ void DemCoProcessor::addInteractor(std::shared_ptr<MovableObjectInteractor> inte
    {
       peGeometry->setLinearVelolocity(initalVelocity);
       physicsEngineGeometries.push_back(peGeometry);
-      //interactor->setActive();
    }
    else
    {
       physicsEngineGeometries.push_back(peGeometry);
-      //interactor->setInactive();
    }
 }
 
@@ -189,14 +187,11 @@ void DemCoProcessor::calculateDemTimeStep(double step) const
 
    for (int i = 0; i < physicsEngineGeometries.size(); i++)
    {
-      //if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
-      {
       physicsEngineSolver->updateGeometry(physicsEngineGeometries[i]);
       if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
       {
          interactors[i]->setPhysicsEngineGeometry(physicsEngineGeometries[i]);
       }
-      }
    }
 }
 
@@ -245,14 +240,11 @@ void DemCoProcessor::distributeIDs()
 
    for (int i = 0; i < physicsEngineGeometries.size(); i++)
    {
-      //if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
-      {
          physicsEngineSolver->updateGeometry(physicsEngineGeometries[i]);
          if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->isActive())
          {
             interactors[i]->setPhysicsEngineGeometry(physicsEngineGeometries[i]);
          }
-      }
    }
 }
 
diff --git a/source/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
index 13ea06bce..5ee510009 100644
--- a/source/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
@@ -82,7 +82,6 @@ void WriteMacroscopicQuantitiesCoProcessor::collectData(double step)
    piece = subfolder + "/" + piece;
 
    std::vector<std::string> cellDataNames;
-   SPtr<Communicator> comm = Communicator::getInstance();
    std::vector<std::string> pieces = comm->gather(piece);
    if (comm->getProcessID() == comm->getRoot())
    {
-- 
GitLab