diff --git a/source/Applications/Thermoplast/config.txt b/source/Applications/Thermoplast/config.txt
index bce2bb29f67465bb4492d7ca4023ca12eab3d009..44596f0d0c83c76123bd51631e3c210519db10ca 100644
--- a/source/Applications/Thermoplast/config.txt
+++ b/source/Applications/Thermoplast/config.txt
@@ -9,8 +9,8 @@ boundingBox = 60 1370 130 190 1530 320 #test bb
  
 blocknx = 10 10 10 
 #blocknx = 300 420 320
-endTime = 100000
-outTime = 100
+endTime = 10000
+outTime = 1000
 availMem = 25e9
 uLB = 0.1
 
diff --git a/source/Applications/Thermoplast/thermoplast.cpp b/source/Applications/Thermoplast/thermoplast.cpp
index 353dc56051d8ad1678cf69d4d8fac37d52140fca..481664758b6885dad2e6729e8479dae7545bbd8d 100644
--- a/source/Applications/Thermoplast/thermoplast.cpp
+++ b/source/Applications/Thermoplast/thermoplast.cpp
@@ -482,7 +482,7 @@ void thermoplast(string configname)
    ////////////////////////////////////////////////////////////////////////////
    ////generating spheres 
    //UBLOG(logINFO, "generating spheres - start, rank="<<myid);
-   SPtr<UbScheduler> sphereScheduler(new UbScheduler(sphereTime));
+   SPtr<UbScheduler> sphereScheduler(new UbScheduler(sphereTime/*10,10,10*/));
    SPtr<CreateDemObjectsCoProcessor> createSphereCoProcessor(new CreateDemObjectsCoProcessor(grid, sphereScheduler, comm, demCoProcessor, sphereMaterial));
    //UBLOG(logINFO, "generating spheres - stop, rank="<<myid);
 
diff --git a/source/DemCoupling/CreateDemObjectsCoProcessor.cpp b/source/DemCoupling/CreateDemObjectsCoProcessor.cpp
index b7401d07f04b1a17a2959b7475cd796d30a3cac7..bd24c8ecd7a19caffafb1699a4089e666771e90f 100644
--- a/source/DemCoupling/CreateDemObjectsCoProcessor.cpp
+++ b/source/DemCoupling/CreateDemObjectsCoProcessor.cpp
@@ -89,7 +89,7 @@ void CreateDemObjectsCoProcessor::createGeoObjects()
    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() };
+      //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))
       //{
@@ -119,6 +119,7 @@ void CreateDemObjectsCoProcessor::createGeoObjects()
       //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/DemCoProcessor.cpp b/source/DemCoupling/DemCoProcessor.cpp
index be3d7207075e4601c612f232ba7ca070b0f22b0d..f6a7f7fe1ee78e786e13fc87a9e20d942b166ede 100644
--- a/source/DemCoupling/DemCoProcessor.cpp
+++ b/source/DemCoupling/DemCoProcessor.cpp
@@ -31,7 +31,7 @@
 #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)
+   CoProcessor(grid, s), comm(comm), forceCalculator(forceCalculator), physicsEngineSolver(std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)), intermediateDemSteps(intermediatePeSteps)
 {
 #ifdef TIMING
    timer.resetAndStart();
@@ -48,7 +48,7 @@ DemCoProcessor::DemCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, SPtr<Comm
       walberla::pe::BodyStorage* bodyStorageShadowCopies = &(*storage)[1];
 
       bodyStorage->registerAddCallback("DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::addPeGeo), this));
-      //bodyStorage->registerRemoveCallback("DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::removePeGeo), this));
+      bodyStorage->registerRemoveCallback("DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::removePeGeo), this));
 
       bodyStorageShadowCopies->registerAddCallback("DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::addPeShadowGeo), this));
       bodyStorageShadowCopies->registerRemoveCallback("DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::removePeShadowGeo), this));
@@ -94,6 +94,7 @@ void DemCoProcessor::addInteractor(std::shared_ptr<MovableObjectInteractor> inte
    if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(peGeometryAdapter)->isActive())
    {
       peGeometryAdapter->setLinearVelolocity(initalVelocity);
+      geoIdMap.insert(std::make_pair(std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(peGeometryAdapter)->getSystemID(), std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(peGeometryAdapter)));
    }
    SetBcBlocksBlockVisitor setBcVisitor(interactor);
    grid->accept(setBcVisitor);
@@ -305,18 +306,36 @@ void DemCoProcessor::moveVfGeoObjects()
    {
       if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive())
       {
-         walberla::pe::RigidBody* peGeoObject = getPeGeoObject(std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->getSystemId());
-         if (peGeoObject != nullptr)
+         ////walberla::pe::RigidBody* peGeoObject = getPeGeoObject(std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->getSystemId());
+         ////if (peGeoObject != nullptr)
+         ////{
+         ////   std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->setGeometry(peGeoObject);
+         //   interactors[i]->moveGbObjectTo(physicsEngineGeometrieAdapters[i]->getPosition());
+         ////}
+         ////else
+         ////{
+         ////   std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->setInactive();
+         ////}
+
+         ////UBLOG(logINFO, "DemCoProcessor::moveVfGeoObject() id = "<<std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->getId()<<"  position="<<physicsEngineGeometrieAdapters[i]->getPosition()<<" rank="<<comm->getProcessID());
+         if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->getSemiactive())
          {
-            std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->setGeometry(peGeoObject);
-            interactors[i]->moveGbObjectTo(physicsEngineGeometrieAdapters[i]->getPosition());
+            walberla::pe::RigidBody* peGeoObject = getPeGeoObject(std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->getSystemID());
+            if (peGeoObject != nullptr)
+            {
+               std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->setGeometry(peGeoObject);
+               interactors[i]->moveGbObjectTo(physicsEngineGeometrieAdapters[i]->getPosition());
+               std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->setSemiactive(false);
+            }
+            else
+            {
+               std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->setInactive();
+            }
          }
          else
          {
-            std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->setInactive();
+            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());
       }
    }
 }
@@ -399,21 +418,23 @@ void DemCoProcessor::getObjectsPropertiesVector(std::vector<double>& p)
 void DemCoProcessor::addPeGeo(walberla::pe::RigidBody * peGeo)
 {
 
+   //UBLOG(logINFO, "DemCoProcessor::addPeGeo() SystemID="<<peGeo->getSystemID());
 
+   //if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
 
-   if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
-
-   auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
-   geometry->counter++;
+   //auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
+   auto geometry = getPeGeoAdapter(peGeo->getSystemID());
+   //geometry->counter++;
    //if (comm->getProcessID()==3)UBLOG(logINFO, "DemCoProcessor::addPeGeo() geoID = "<<peGeo->getID()<<" counter="<<geometry->counter<<" shadowCounter="<<geometry->shadowCounter);
-   if (geometry->getSystemId() == peGeo->getSystemID())
+   if (geometry != nullptr) //(geometry->getSystemId() == peGeo->getSystemID())
    {
       geometry->setActive();
       geometry->setGeometry(peGeo);
+      //geometry->counter++;
       return;
    }
    else
-      throw UbException(UB_EXARGS, "PeGeo ID="+UbSystem::toString(peGeo->getID())+" is not matching geometry ID="+UbSystem::toString(geometry->getId()));
+      return;//throw UbException(UB_EXARGS, "PeGeo SystemId="+UbSystem::toString(peGeo->getSystemID())+" is not matching geometry ID");
 }
 
 void DemCoProcessor::removePeGeo(walberla::pe::RigidBody * peGeo)
@@ -422,70 +443,64 @@ void DemCoProcessor::removePeGeo(walberla::pe::RigidBody * peGeo)
 
    //if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
 
+   //UBLOG(logINFO, "DemCoProcessor::removePeGeo() SystemID="<<peGeo->getSystemID());
 
-   auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
-   geometry->counter--;
+   //auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
+   auto geometry = getPeGeoAdapter(peGeo->getSystemID());
+   //geometry->counter--;
 
-   if (geometry->getSystemId() == peGeo->getSystemID())
+   if (geometry != nullptr)//(geometry->getSystemID() == peGeo->getSystemID())
    {
-      if (geometry->shadowCounter == 0)
-      {
-         geometry->setInactive();
-         geometry->setGeometry(NULL);
-      }
-      //if (comm->getProcessID()==3)UBLOG(logINFO, "DemCoProcessor::removePeGeo() geoID = "<<peGeo->getID()<<" counter="<<geometry->counter<<" shadowCounter="<<geometry->shadowCounter);
-      return;
+      geometry->setSemiactive(true);
+      //walberla::pe::RigidBody* peGeoObject = getPeGeoObject(peGeo->getSystemID());
+      //if (peGeoObject != nullptr)
+      //{
+      //   geometry->setGeometry(peGeoObject);
+      //}
+      //else
+      //{
+      //   geometry->setInactive();
+      //   geometry->setGeometry(nullptr);
+      //}
+      //////if (geometry->shadowCounter == 0)
+      //////{
+      ////   geometry->setInactive();
+      ////   geometry->setGeometry(nullptr);
+      //////   geometry->counter--;
+      //////}
+      ////if (comm->getProcessID()==3)UBLOG(logINFO, "DemCoProcessor::removePeGeo() geoID = "<<peGeo->getID()<<" counter="<<geometry->counter<<" shadowCounter="<<geometry->shadowCounter);
+      //return;
    }
    else
-      throw UbException(UB_EXARGS, "PeGeo ID="+UbSystem::toString(peGeo->getID())+" is not matching geometry ID="+UbSystem::toString(geometry->getId()));
-
+      throw UbException(UB_EXARGS, "PeGeo SystemId="+UbSystem::toString(peGeo->getSystemID())+" is not matching geometry ID");
+}
 
+void DemCoProcessor::addPeShadowGeo(walberla::pe::RigidBody * peGeo)
+{
    //if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
 
-   //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();
+   //UBLOG(logINFO, "DemCoProcessor::addPeShadowGeo()");
+   //UBLOG(logINFO, "DemCoProcessor::addPeShadowGeo() SystemID="<<peGeo->getSystemID());
 
-   //for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
-   //{
-   //   for (auto bodyIt = walberla::pe::BodyIterator::begin(*blockIt, *storageId.get()); bodyIt != walberla::pe::BodyIterator::end(); ++bodyIt)
-   //   {
-   //      if (bodyIt->getID() == peGeo->getID())
-   //      {
-   //         return;
-   //      }
-   //   }
-   //}
 
    //auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
-   //if (geometry->getId() == peGeo->getID())
-   //{
-   //   geometry->setInactive();
-   //   geometry->setGeometry(NULL);
-   //   return;
-   //}
-   //else
-   //   throw UbException(UB_EXARGS, "PeGeo ID is not matching!");
-}
-
-void DemCoProcessor::addPeShadowGeo(walberla::pe::RigidBody * peGeo)
-{
-   if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
+   auto geometry = getPeGeoAdapter(peGeo->getSystemID());
 
-   auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
-   geometry->shadowCounter++;
    //if (comm->getProcessID()==3) UBLOG(logINFO, "DemCoProcessor::addPeShadowGeo() geoID = "<<peGeo->getID()<<" counter="<<geometry->counter<<" shadowCounter="<<geometry->shadowCounter);
-   if (geometry->getSystemId() == peGeo->getSystemID())
+   if (geometry != nullptr)//(geometry->getSystemId() == peGeo->getSystemID())
    {
-      if (geometry->shadowCounter == 1)
-      {
-         geometry->setActive();
-         geometry->setGeometry(peGeo);
-      }
+      //geometry->shadowCounter++;
+      //if (geometry->shadowCounter == 1)
+      //{
+      geometry->setActive();
+      geometry->setGeometry(peGeo);
+      //geometry->shadowCounter++;
+   //}
 
       return;
    }
    else
-      throw UbException(UB_EXARGS, "PeGeo ID="+UbSystem::toString(peGeo->getID())+" is not matching geometry ID="+UbSystem::toString(geometry->getId()));
+      throw UbException(UB_EXARGS, "PeGeo ID="+UbSystem::toString(peGeo->getSystemID())+" is not matching geometry ID");
 }
 
 void DemCoProcessor::removePeShadowGeo(walberla::pe::RigidBody * peGeo)
@@ -494,64 +509,50 @@ void DemCoProcessor::removePeShadowGeo(walberla::pe::RigidBody * peGeo)
 
    //if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
 
-   //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)
-   //{
-   //   for (auto bodyIt = walberla::pe::BodyIterator::begin(*blockIt, *storageId.get()); bodyIt != walberla::pe::BodyIterator::end(); ++bodyIt)
-   //   {
-   //      if (bodyIt->getID() == peGeo->getID())
-   //      {
-   //         return;
-   //      }
-   //   }
-   //}
-
-
-
-   //std::vector< std::shared_ptr<Block3D> > blockVector;
-   //UbTupleInt3 blockNX=grid->getBlockNX();
-   //SPtr<GbObject3D> geoObject(interactors[peGeo->getID()]->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);
-   //if (blockVector.size() == 0)
-   //{
-   //   auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
-   //   if (geometry->getId() == peGeo->getID())
-   //   {
-   //      geometry->setInactive();
-   //      geometry->setGeometry(NULL);
-   //      return;
-   //   }
-   //   else
-   //      throw UbException(UB_EXARGS, "PeGeo ID is not matching!");
-   //}
-
+   //auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
 
-   //if (peGeo->getID() < 0 || peGeo->getID() >= physicsEngineGeometrieAdapters.size()) return;
+   //UBLOG(logINFO, "DemCoProcessor::removePeShadowGeo() SystemID="<<peGeo->getSystemID());
 
-   auto geometry = std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[peGeo->getID()]);
-   geometry->shadowCounter--;
+   auto geometry = getPeGeoAdapter(peGeo->getSystemID());
 
-   if (geometry->getSystemId() == peGeo->getSystemID())
+   if (geometry != nullptr) //(geometry->getSystemId() == peGeo->getSystemID())
    {
-      if (geometry->shadowCounter == 0)
-      {
-         geometry->setInactive();
-         geometry->setGeometry(NULL);
-      }
-      else
-      {
-         geometry->setActive();
-         geometry->setGeometry(peGeo);
-      }
-      //if (comm->getProcessID()==3)UBLOG(logINFO, "DemCoProcessor::removePeShadowGeo() geoID = "<<peGeo->getID()<<" counter="<<geometry->counter<<" shadowCounter="<<geometry->shadowCounter);
-      return;
+      geometry->setSemiactive(true);
+      //walberla::pe::RigidBody* peGeoObject = getPeGeoObject(peGeo->getSystemID());
+      //if (peGeoObject != nullptr)
+      //{
+      //   geometry->setGeometry(peGeoObject);
+      //}
+      //else
+      //{
+      //   geometry->setInactive();
+      //   geometry->setGeometry(nullptr);
+      //}
+      //////if (geometry->shadowCounter == 0)
+      //////{
+      ////   geometry->setInactive();
+      ////   geometry->setGeometry(nullptr);
+      //////   geometry->counter--;
+      //////}
+      ////if (comm->getProcessID()==3)UBLOG(logINFO, "DemCoProcessor::removePeGeo() geoID = "<<peGeo->getID()<<" counter="<<geometry->counter<<" shadowCounter="<<geometry->shadowCounter);
+      //return;
+      //////geometry->shadowCounter--;
+      //////if (geometry->shadowCounter == 0)
+      //////{
+      ////   geometry->setInactive();
+      ////   geometry->setGeometry(nullptr);
+      //////}
+      //////else
+      //////{
+      //////   geometry->setActive();
+      //////   geometry->setGeometry(peGeo);
+      ////   //geometry->shadowCounter--;
+      //////}
+      //////if (comm->getProcessID()==3)UBLOG(logINFO, "DemCoProcessor::removePeShadowGeo() geoID = "<<peGeo->getID()<<" counter="<<geometry->counter<<" shadowCounter="<<geometry->shadowCounter);
+      ////return;
    }
    else
-      throw UbException(UB_EXARGS, "PeGeo ID="+UbSystem::toString(peGeo->getID())+" is not matching geometry ID="+UbSystem::toString(geometry->getId()));
+      throw UbException(UB_EXARGS, "PeGeo ID="+UbSystem::toString(peGeo->getSystemID())+" is not matching geometry ID");
 }
 
 bool DemCoProcessor::isSpheresIntersection(double centerX1, double centerX2, double centerX3, double d)
@@ -592,7 +593,7 @@ void DemCoProcessor::distributeIDs()
    {
       if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive())
       {
-         peIDsSend.push_back(std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->getSystemId());
+         peIDsSend.push_back(std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->getSystemID());
          vfIDsSend.push_back(interactors[i]->getID());
       }
    }
@@ -615,11 +616,13 @@ void DemCoProcessor::distributeIDs()
       std::map<int, unsigned long long>::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!");
+         throw UbException(UB_EXARGS, "Interactor ID = "+UbSystem::toString(interactors[i]->getID())+" is invalid! The DEM object may be not in PE domain!");
       }
 
+      std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->setSystemID(it->second);
+
+      geoIdMap.insert(std::make_pair(it->second, std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])));
 
-      std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->setSystemId(it->second);
       //std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[i])->setId(idMap.find(interactors[i]->getID())->second);
    }
 
@@ -632,7 +635,7 @@ void DemCoProcessor::distributeIDs()
    //   }
    //}
 }
-
+//////////////////////////////////////////////////////////////////////////
 void DemCoProcessor::setBlockVisitor(std::shared_ptr<BoundaryConditionsBlockVisitor> boundaryConditionsBlockVisitor)
 {
    this->boundaryConditionsBlockVisitor = boundaryConditionsBlockVisitor;
@@ -642,60 +645,19 @@ walberla::pe::RigidBody* DemCoProcessor::getPeGeoObject(walberla::id_t id)
 {
    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)
-   //{
-
-   //   for (auto bodyIt = walberla::pe::BodyIterator::begin(*blockIt, *storageId.get()); bodyIt != walberla::pe::BodyIterator::end(); ++bodyIt)
-   //   {
-   //      if (bodyIt->getID() == id)
-   //      {
-   //         walberla::pe::RigidBody* geo = *(bodyIt);
-   //         return geo;
-   //      }
-   //   }
-   //}
-
    std::shared_ptr<walberla::pe::BodyStorage> globalBodyStorage = std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getGlobalBodyStorage();
-   //walberla::pe::BodyID peGeo =  walberla::pe::getBody(*globalBodyStorage, *forest, *storageId, id);
-   //if (peGeo != nullptr)
-   //{
-   //   return peGeo;
-   //}
-   //return NULL;
-
-
-   return walberla::pe::getBody(*globalBodyStorage, *forest, *storageId, id);
-
-
-
-   //walberla::pe::RigidBody* peBody = NULL;
-
-   //std::shared_ptr<walberla::pe::BodyStorage> globalStorage = std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getGlobalBodyStorage();
-
-   //auto bodyIt = globalStorage->find(id);
-   //if (bodyIt != globalStorage->end())
-   //{
-   //   return *bodyIt;
-   //}
-
-   //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];
-   //   walberla::pe::BodyStorage* bodystorageShadowCopies = &(*storage)[1];
-   //   
-   //   auto bodyIt = bodyStorage->find(id);
-   //   if (bodyIt != bodyStorage->end())
-   //   {
-   //      return *bodyIt;
-   //   }
 
-   //   bodyIt = bodystorageShadowCopies->find(id);
-   //   if (bodyIt != bodystorageShadowCopies->end())
-   //   {
-   //      return *bodyIt;
-   //   }
-   //}
-   //return NULL;
-}
\ No newline at end of file
+   return walberla::pe::getBody(*globalBodyStorage, *forest, *storageId, id, walberla::pe::StorageSelect::LOCAL | walberla::pe::StorageSelect::SHADOW);
+}
+////////////////////////////////////////////////////////////////////////////
+std::shared_ptr<PePhysicsEngineGeometryAdapter> DemCoProcessor::getPeGeoAdapter(unsigned long long systemId)
+{
+   std::map< unsigned long long, std::shared_ptr<PePhysicsEngineGeometryAdapter> >::const_iterator it;
+   if ((it=geoIdMap.find(systemId)) == geoIdMap.end())
+   {
+      //throw UbException(UB_EXARGS, "pe SystemId can not be found!");
+      return nullptr;
+   }
+   else
+      return it->second;
+}
diff --git a/source/DemCoupling/DemCoProcessor.h b/source/DemCoupling/DemCoProcessor.h
index 932bade259c20df9e1df494dee04979422afc101..74693a5de13252100c9828cffe69b71b44656272 100644
--- a/source/DemCoupling/DemCoProcessor.h
+++ b/source/DemCoupling/DemCoProcessor.h
@@ -25,7 +25,9 @@
 
 class PhysicsEngineGeometryAdapter;
 class PhysicsEngineSolverAdapter;
+class PePhysicsEngineSolverAdapter;
 class PhysicsEngineMaterialAdapter;
+class PePhysicsEngineGeometryAdapter;
 
 class UbScheduler;
 class Grid3D;
@@ -64,17 +66,20 @@ private:
     void calculateDemTimeStep(double step);
     void moveVfGeoObjects();
     walberla::pe::RigidBody* getPeGeoObject(walberla::id_t id);
+    std::shared_ptr<PePhysicsEngineGeometryAdapter> getPeGeoAdapter(unsigned long long systemId);
 private:
     std::shared_ptr<Communicator> comm;
     std::vector<std::shared_ptr<MovableObjectInteractor> > interactors;
     std::shared_ptr<ForceCalculator> forceCalculator;
-    std::shared_ptr<PhysicsEngineSolverAdapter> physicsEngineSolver;
+    std::shared_ptr<PePhysicsEngineSolverAdapter> physicsEngineSolver;
     std::vector<std::shared_ptr<PhysicsEngineGeometryAdapter> > physicsEngineGeometrieAdapters;
     double intermediateDemSteps;
     SPtr<BoundaryConditionsBlockVisitor> boundaryConditionsBlockVisitor;
     //walberla::pe::BodyStorage* bodyStorage;    //!< Reference to the central body storage.
     //walberla::pe::BodyStorage* bodyStorageShadowCopies;    //!< Reference to the body storage containing body shadow copies.
 
+    std::map< unsigned long long, std::shared_ptr<PePhysicsEngineGeometryAdapter> > geoIdMap;
+
 #ifdef TIMING
     UbTimer timer;
 #endif
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp
index cbe0a325874a54d6ec0b1fb66a38b8227270eebc..a044a6ea2dde5c3c6d352d20edb19e8e6e378360 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp
@@ -14,8 +14,9 @@
 PePhysicsEngineGeometryAdapter::PePhysicsEngineGeometryAdapter()
 {
    this->id = -999;
-   this->systemId = -999;
+   this->systemID = -999;
    this->active = false;
+   this->semiactive = false;
    shadowCounter = 0;
    counter = 0;
 }
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h
index feea998f902724217bd88fffcea928820bf3bb3d..c5427ab823c23bd559d0b3bf6cee6c5d9c0c0da2 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h
@@ -57,16 +57,18 @@ public:
     int shadowCounter;
     int counter;
 
-    unsigned long long getSystemId() const { return systemId; }
-    void setSystemId(unsigned long long val) { systemId = val; }
+    unsigned long long getSystemID() const { return systemID; }
+    void setSystemID(unsigned long long val) { systemID = val; }
+    bool getSemiactive() const { return semiactive; }
+    void setSemiactive(bool val) { semiactive = val; }
 private:
     walberla::pe::RigidBody* peGeoObject;
     //unsigned long long id;
     int id;
     //walberla::id_t systemId;
-    unsigned long long systemId;
+    unsigned long long systemID;
     bool active;
-
+    bool semiactive;
 };
 
 #endif
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
index 7e5f5fccef7f7b4c9da9498d2c30ab284479715f..02ce4a783b66c0c77b40084c6c7201ccc8cfc2c1 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
@@ -36,15 +36,19 @@ void PePhysicsEngineSolverAdapter::initalizePeEnvironment()
 std::shared_ptr<PhysicsEngineGeometryAdapter> PePhysicsEngineSolverAdapter::createPhysicsEngineGeometryAdapter(int id, const Vector3D& position, double radius, std::shared_ptr<PhysicsEngineMaterialAdapter> material) const
 {
     const std::shared_ptr<PePhysicsEngineMaterialAdapter> peMaterial = std::dynamic_pointer_cast<PePhysicsEngineMaterialAdapter>(material);
+    std::shared_ptr<PePhysicsEngineGeometryAdapter> peGeometryAdapter(new PePhysicsEngineGeometryAdapter());
 
+    //UBLOG(logINFO, "PePhysicsEngineSolverAdapter::createSphere():start");
     walberla::pe::GeomID peGeometry = createSphere(*globalBodyStorage, *forest, *storageId, id, PeConverter::convert(position), radius, peMaterial->getPeMaterial());
-
-    std::shared_ptr<PePhysicsEngineGeometryAdapter> peGeometryAdapter(new PePhysicsEngineGeometryAdapter());
+    //UBLOG(logINFO, "PePhysicsEngineSolverAdapter::createSphere():end");
 
     if (peGeometry)
     {
        peGeometryAdapter->setId(id);
-       peGeometryAdapter->setSystemId(peGeometry->getSystemID());
+       peGeometryAdapter->setSystemID(peGeometry->getSystemID());
+       //unsigned long long sid = peGeometryAdapter->getSystemID();
+       //geoIdMap.insert(std::make_pair(sid, peGeometryAdapter));
+       //peGeometryAdapter->counter++;
        
        //if(peGeometry->getSystemID() != id+1)
        //   UB_THROW(UbException(UB_EXARGS, "id="+UbSystem::toString(id)+" does not match pe::SystemId="+UbSystem::toString(peGeometry->getSystemID())));
@@ -80,6 +84,10 @@ void PePhysicsEngineSolverAdapter::initialPeBodyStorage()
 
 void PePhysicsEngineSolverAdapter::initialPeBlockForest()
 {
+
+   //walberla::SetupBlockForest sforest = walberla::blockforest::createUniformBlockGrid(walberla::AABB(peParameter->simulationDomain[0], peParameter->simulationDomain[1], peParameter->simulationDomain[2],
+   //   peParameter->simulationDomain[3], peParameter->simulationDomain[4], peParameter->simulationDomain[5]), // simulationDomain
+   //   walberla::uint_t(val<1>(peParameter->numberOfBlocks)), walberla::uint_t(val<2>(peParameter->numberOfBlocks)), walberla::uint_t(val<3>(peParameter->numberOfBlocks)),walberla::uint_t(10),walberla::uint_t(10),walberla::uint_t(10), 5.0,false);
     walberla::SetupBlockForest sforest;
     //sforest.addWorkloadMemorySUIDAssignmentFunction( uniformWorkloadAndMemoryAssignment );
     sforest.init(walberla::AABB(peParameter->simulationDomain[0], peParameter->simulationDomain[1], peParameter->simulationDomain[2],
@@ -200,3 +208,15 @@ std::shared_ptr<walberla::pe::BodyStorage> PePhysicsEngineSolverAdapter::getGlob
 {
    return globalBodyStorage;
 }
+//////////////////////////////////////////////////////////////////////////
+//std::shared_ptr<PePhysicsEngineGeometryAdapter> PePhysicsEngineSolverAdapter::getPeGeoAdapter(unsigned long long systemId)
+//{
+//   std::map< unsigned long long, std::shared_ptr<PePhysicsEngineGeometryAdapter> >::const_iterator it;
+//   if ((it=geoIdMap.find(systemId)) == geoIdMap.end())
+//   {
+//      //throw UbException(UB_EXARGS, "pe SystemId can not be found!");
+//      return nullptr;
+//   }
+//   else
+//      return it->second;
+//}
\ No newline at end of file
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h
index 21f674edfd703c2f5a888ac9c586b4cbcca9410a..7ddeed3fb948dcbed20592114fd94b8bcc6e42b4 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h
@@ -17,6 +17,7 @@
 
 class PePhysicsEngineMaterialAdapter;
 class PhysicsEngineGeometryAdapter;
+class PePhysicsEngineGeometryAdapter;
 class PeLoadBalancerAdapter;
 
 namespace walberla
@@ -76,6 +77,7 @@ public:
     std::shared_ptr<walberla::blockforest::BlockForest> getBlockForest();
     std::shared_ptr<walberla::domain_decomposition::BlockDataID> getStorageId();
     std::shared_ptr<walberla::pe::BodyStorage> getGlobalBodyStorage();
+    //std::shared_ptr<PePhysicsEngineGeometryAdapter> getPeGeoAdapter(unsigned long long systemId);
 
 private:
     void initalizePeEnvironment();
@@ -86,6 +88,7 @@ private:
     void initalPeIntegrator();
     static void executePeBodyTypeTuple();
     void initialPeChannel() const;
+    
 
 private:
     std::shared_ptr<PeParameter> peParameter;
@@ -95,6 +98,7 @@ private:
     std::shared_ptr< walberla::blockforest::BlockForest > forest;
     std::shared_ptr<walberla::domain_decomposition::BlockDataID> storageId;
     std::shared_ptr<walberla::pe::cr::HardContactSemiImplicitTimesteppingSolvers> cr;
+    // std::map< unsigned long long, std::shared_ptr<PePhysicsEngineGeometryAdapter> > geoIdMap;
 };
 
 #endif