diff --git a/source/Applications/Thermoplast/config.txt b/source/Applications/Thermoplast/config.txt
index 954ed682104a5a01602f668369a953a0d4f7315d..c681fec35a9f03f83e1394028341494f6ee6578f 100644
--- a/source/Applications/Thermoplast/config.txt
+++ b/source/Applications/Thermoplast/config.txt
@@ -22,7 +22,7 @@ pathGeo = d:/Projects/ThermoPlast/SimGeo
 michel = /michel.stl
 plexiglas = /plexiglas.stl
 
-pathOut = g:/temp/thermoplastCluster5
+pathOut = g:/temp/thermoplastCluster6
 
 #restart
 cpStart = 20000000
diff --git a/source/Applications/Thermoplast/thermoplast.cpp b/source/Applications/Thermoplast/thermoplast.cpp
index fac556f3c0fefc2f1ca7fbe34793ec89b6dd9655..d8f1e10cce7d287c25e1fef2bdd07623f4a3e22c 100644
--- a/source/Applications/Thermoplast/thermoplast.cpp
+++ b/source/Applications/Thermoplast/thermoplast.cpp
@@ -20,6 +20,7 @@
 #include <PePhysicsEngineMaterialAdapter.h>
 #include <PePhysicsEngineGeometryAdapter.h>
 #include <PePhysicsEngineSolverAdapter.h>
+#include "PeLoadBalancerAdapter.h"
 
 #include <VelocityBcReconstructor.h>
 #include <EquilibriumReconstructor.h>
@@ -79,7 +80,8 @@ std::shared_ptr<DemCoProcessor> makePeCoProcessor(SPtr<Grid3D> grid, SPtr<Commun
 
    std::shared_ptr<PeParameter> peParamter = std::make_shared<PeParameter>(peRelaxtion, maxpeIterations, globalLinearAcc,
       planeMaterial, simulationDomain, numberOfBlocks, isPeriodic, minOffset, maxOffset);
-   std::shared_ptr<PhysicsEngineSolverAdapter> peSolver = std::make_shared<PePhysicsEngineSolverAdapter>(peParamter);
+   std::shared_ptr<PeLoadBalancerAdapter> loadBalancer(new PeLoadBalancerAdapter(grid, comm->getNumberOfProcesses(), comm->getProcessID()));
+   std::shared_ptr<PhysicsEngineSolverAdapter> peSolver = std::make_shared<PePhysicsEngineSolverAdapter>(peParamter, loadBalancer);
 
    const std::shared_ptr<ForceCalculator> forceCalculator = std::make_shared<ForceCalculator>(comm);
 
@@ -197,19 +199,6 @@ void thermoplast(string configname)
    GenBlocksGridVisitor genBlocks(gridCube);
    grid->accept(genBlocks);
 
-   //PE initialization
-   double refLengthLb = radius*2.0;
-   double refLengthWorld = refLengthLb * deltax;
-   const std::shared_ptr<LBMUnitConverter> lbmUnitConverter = std::make_shared<LBMUnitConverter>(refLengthWorld, LBMUnitConverter::WORLD_MATERIAL::AIR_20C, refLengthLb);
-   if (myid == 0) std::cout << lbmUnitConverter->toString() << std::endl;
-   double rhoSphere = 915 * lbmUnitConverter->getFactorDensityWToLb();  // kg/m^3
-   if (myid == 0) UBLOG(logINFO, "rhoSphere = "<<rhoSphere);
-   SPtr<PhysicsEngineMaterialAdapter> sphereMaterial(new PePhysicsEngineMaterialAdapter("Polypropylen", rhoSphere, 0, 0.15, 0.1, 0.45, 0.5, 1, 0, 0));
-   const int timestep = 2;
-   const SPtr<UbScheduler> peScheduler(new UbScheduler(timestep));
-   int maxpeIterations = 10;//endTime/2;
-   SPtr<DemCoProcessor> demCoProcessor = makePeCoProcessor(grid, comm, peScheduler, lbmUnitConverter, maxpeIterations);
-   demCoProcessor->setBlockVisitor(bcVisitor);
 
    if (myid == 0)
    {
@@ -297,8 +286,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);
@@ -367,6 +356,20 @@ void thermoplast(string configname)
       if (myid == 0) UBLOG(logINFO, "Preprocess - end");
    }
 
+   //PE initialization
+   double refLengthLb = radius*2.0;
+   double refLengthWorld = refLengthLb * deltax;
+   const std::shared_ptr<LBMUnitConverter> lbmUnitConverter = std::make_shared<LBMUnitConverter>(refLengthWorld, LBMUnitConverter::WORLD_MATERIAL::AIR_20C, refLengthLb);
+   if (myid == 0) std::cout << lbmUnitConverter->toString() << std::endl;
+   double rhoSphere = 915 * lbmUnitConverter->getFactorDensityWToLb();  // kg/m^3
+   if (myid == 0) UBLOG(logINFO, "rhoSphere = "<<rhoSphere);
+   SPtr<PhysicsEngineMaterialAdapter> sphereMaterial(new PePhysicsEngineMaterialAdapter("Polypropylen", rhoSphere, 0, 0.15, 0.1, 0.45, 0.5, 1, 0, 0));
+   const int timestep = 2;
+   const SPtr<UbScheduler> peScheduler(new UbScheduler(timestep));
+   int maxpeIterations = 10;//endTime/2;
+   SPtr<DemCoProcessor> demCoProcessor = makePeCoProcessor(grid, comm, peScheduler, lbmUnitConverter, maxpeIterations);
+   demCoProcessor->setBlockVisitor(bcVisitor);
+
    ////////////////////////////////////////////////////////////////////////////
    ////generating spheres 
    //UBLOG(logINFO, "generating spheres - start, rank="<<myid);
diff --git a/source/DemCoupling/MovableObjectInteractor.cpp b/source/DemCoupling/MovableObjectInteractor.cpp
index c6a2b6933c108c0f906e2e8190d049e1047c8b27..562d577cfacf7ab0f8cf5708e4707f80d5eac4e9 100644
--- a/source/DemCoupling/MovableObjectInteractor.cpp
+++ b/source/DemCoupling/MovableObjectInteractor.cpp
@@ -111,6 +111,23 @@ void MovableObjectInteractor::rearrangeGrid()
 #endif
 }
 
+void MovableObjectInteractor::updateNodeLists()
+{
+   //for (BcNodeIndicesMap::value_type t : bcNodeIndicesMap)
+   //{
+   //   SPtr<Block3D> block = t.first;
+   //   std::set< UbTupleInt3 >& bcNodeIndices = t.second;
+
+
+   //   SPtr<ILBMKernel> kernel = block->getKernel();
+
+   //   for (UbTupleInt3 node : bcNodeIndices)
+   //   {
+
+   //   }
+   //}
+}
+
 void MovableObjectInteractor::reconstructDistributionOnSolidNodes()
 {
     for(SolidNodeIndicesMap::value_type t : solidNodeIndicesMap)
diff --git a/source/DemCoupling/MovableObjectInteractor.h b/source/DemCoupling/MovableObjectInteractor.h
index a77fb4502de03f91a85997885de1813232461aa4..0aeb22bbdbd8003efb9d2a245aa013160114247e 100644
--- a/source/DemCoupling/MovableObjectInteractor.h
+++ b/source/DemCoupling/MovableObjectInteractor.h
@@ -25,6 +25,9 @@ class Reconstructor;
 
 class MovableObjectInteractor : public D3Q27Interactor
 {
+public:
+   typedef std::map<SPtr<Block3D>, std::set< std::array<int,3> > > InBcNodeIndicesMap;
+   typedef std::map<SPtr<Block3D>, std::set< std::array<int,3> > > OutBcNodeIndicesMap;
 public:
     MovableObjectInteractor(std::shared_ptr<GbObject3D> geoObject3D, std::shared_ptr<Grid3D> grid, std::shared_ptr<BCAdapter> bcAdapter, int type, std::shared_ptr<Reconstructor> reconstructor, State isPinned);
     virtual ~MovableObjectInteractor();
@@ -35,6 +38,7 @@ public:
 
 private:
     void rearrangeGrid();
+    void updateNodeLists();
     void setSolidNodesToFluid();
     void setBcNodesToFluid();
     void reconstructDistributionOnSolidNodes();
diff --git a/source/DemCoupling/PePartitioningGridVisitor.cpp b/source/DemCoupling/PePartitioningGridVisitor.cpp
index df05d5ef2aa9241ddeec6ebee5d506d7304c2663..b52d3ba575c361763b27eedf07903d68e82d08c6 100644
--- a/source/DemCoupling/PePartitioningGridVisitor.cpp
+++ b/source/DemCoupling/PePartitioningGridVisitor.cpp
@@ -38,8 +38,8 @@ void PePartitioningGridVisitor::visit(SPtr<Grid3D> grid)
 //////////////////////////////////////////////////////////////////////////
 void PePartitioningGridVisitor::collectData(SPtr<Grid3D> grid)
 {
-   int minInitLevel = grid->getCoarsestInitializedLevel();
-   int maxInitLevel = grid->getFinestInitializedLevel();
+   //int minInitLevel = grid->getCoarsestInitializedLevel();
+   //int maxInitLevel = grid->getFinestInitializedLevel();
 
    walberla::uint_t peRank;
 
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.cpp b/source/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9283b3b15ff8676356a3ee75b006fd22c29fdfd2
--- /dev/null
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.cpp
@@ -0,0 +1,41 @@
+#include "PeLoadBalancerAdapter.h"
+#include "Grid3D.h"
+#include "Block3D.h"
+#include "CoordinateTransformation3D.h"
+
+#include "core/debug/CheckFunctions.h"
+
+
+PeLoadBalancerAdapter::PeLoadBalancerAdapter(SPtr<Grid3D> grid, unsigned numberOfProcesses, int rank) : grid(grid), numberOfProcesses(numberOfProcesses), rank(rank)
+{
+
+}
+
+walberla::uint_t PeLoadBalancerAdapter::operator()(walberla::SetupBlockForest & forest, const walberla::uint_t numberOfProcesses, const walberla::memory_t perProcessMemoryLimit)
+{
+   std::vector< walberla::SetupBlock * > peBlocks;
+   forest.getBlocks(peBlocks);
+
+   for (auto peBlock = peBlocks.begin(); peBlock != peBlocks.end(); ++peBlock)
+   {
+      walberla::AABB aabb = (*peBlock)->getAABB();
+      SPtr<Block3D> block = getBlockByMinUniform((double)aabb.xMin(), (double)aabb.yMin(), (double)aabb.zMin(), grid);
+      if (block)
+      {
+         (*peBlock)->assignTargetProcess((walberla::uint_t)block->getRank());
+      }
+   }
+
+   return numberOfProcesses;
+}
+
+SPtr<Block3D> PeLoadBalancerAdapter::getBlockByMinUniform(double minX1, double minX2, double minX3, SPtr<Grid3D> grid)
+{
+   SPtr<CoordinateTransformation3D> trafo = grid->getCoordinateTransformator();
+
+   int ix1 = (int)trafo->transformForwardToX1Coordinate(minX1, minX2, minX3);
+   int ix2 = (int)trafo->transformForwardToX2Coordinate(minX1, minX2, minX3);
+   int ix3 = (int)trafo->transformForwardToX3Coordinate(minX1, minX2, minX3);
+
+   return grid->getBlock(ix1, ix2, ix3, 0);
+}
\ No newline at end of file
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.h b/source/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.h
new file mode 100644
index 0000000000000000000000000000000000000000..6b904d994de933052435190cb4436ada882e023f
--- /dev/null
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.h
@@ -0,0 +1,25 @@
+#ifndef PeLoadBalancerAdapter_h__
+#define PeLoadBalancerAdapter_h__
+
+#include "blockforest/SetupBlockForest.h"
+#include "PointerDefinitions.h"
+
+class Grid3D;
+class Block3D;
+
+class PeLoadBalancerAdapter
+{
+public:
+   PeLoadBalancerAdapter(SPtr<Grid3D> grid, unsigned numberOfProcesses, int rank);
+   walberla::uint_t operator()( walberla::SetupBlockForest & forest, const walberla::uint_t numberOfProcesses, const walberla::memory_t perProcessMemoryLimit );
+   unsigned getNumberOfProcesses() const { return numberOfProcesses; }
+   int getRank() const { return rank; }
+protected:
+   SPtr<Block3D> getBlockByMinUniform(double minX1, double minX2, double minX3, SPtr<Grid3D> grid);
+private:
+   SPtr<Grid3D> grid;
+   unsigned numberOfProcesses;
+   int rank;
+};
+
+#endif // PeLoadBalancerAdapter_h__
\ No newline at end of file
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
index 82f3ff85c6c2ef88fc172bb82c73b49ed71c9b9b..cb4ee4bfdba3ab3d601037ac078a9925716de028 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
@@ -6,6 +6,7 @@
 #include "PeAdapter.h"
 #include "PePhysicsEngineGeometryAdapter.h"
 #include "PePhysicsEngineMaterialAdapter.h"
+#include "PeLoadBalancerAdapter.h"
 #include "Communicator.h"
 #include "UbLogger.h"
 #include <boost/tuple/tuple.hpp>
@@ -15,7 +16,7 @@
 typedef boost::tuple<walberla::pe::Sphere, walberla::pe::Plane> BodyTypeTuple;
 
 
-PePhysicsEngineSolverAdapter::PePhysicsEngineSolverAdapter(std::shared_ptr<PeParameter> peParameter) : peParameter(peParameter)
+PePhysicsEngineSolverAdapter::PePhysicsEngineSolverAdapter(std::shared_ptr<PeParameter> peParameter, std::shared_ptr<PeLoadBalancerAdapter> loadBalancer) : peParameter(peParameter), loadBalancer(loadBalancer)
 {
     this->initalizePeEnvironment();
 }
@@ -80,17 +81,30 @@ void PePhysicsEngineSolverAdapter::initialPeBodyStorage()
 
 void PePhysicsEngineSolverAdapter::initalPeBlockForest()
 {
-    forest = walberla::pe::createBlockForest
-    (
-        //walberla::AABB(0, 0, 0, val<1>(peParameter->simulationDomain), val<2>(peParameter->simulationDomain), val<3>(peParameter->simulationDomain)), // simulationDomain
-       walberla::AABB(peParameter->simulationDomain[0], peParameter->simulationDomain[1], peParameter->simulationDomain[2], 
-          peParameter->simulationDomain[3], peParameter->simulationDomain[4], peParameter->simulationDomain[5]), // simulationDomain
-        walberla::Vector3<walberla::uint_t>(val<1>(peParameter->numberOfBlocks), val<2>(peParameter->numberOfBlocks), val<3>(peParameter->numberOfBlocks)), // blocks in each direction
-        walberla::Vector3<bool>(val<1>(peParameter->isPeriodic), val<2>(peParameter->isPeriodic), val<3>(peParameter->isPeriodic)) // periodicity
-    ); 
-
+    //forest = walberla::pe::createBlockForest
+    //(
+    //    //walberla::AABB(0, 0, 0, val<1>(peParameter->simulationDomain), val<2>(peParameter->simulationDomain), val<3>(peParameter->simulationDomain)), // simulationDomain
+    //   walberla::AABB(peParameter->simulationDomain[0], peParameter->simulationDomain[1], peParameter->simulationDomain[2], 
+    //      peParameter->simulationDomain[3], peParameter->simulationDomain[4], peParameter->simulationDomain[5]), // simulationDomain
+    //    walberla::Vector3<walberla::uint_t>(val<1>(peParameter->numberOfBlocks), val<2>(peParameter->numberOfBlocks), val<3>(peParameter->numberOfBlocks)), // blocks in each direction
+    //    walberla::Vector3<bool>(val<1>(peParameter->isPeriodic), val<2>(peParameter->isPeriodic), val<3>(peParameter->isPeriodic)) // periodicity
+    //); 
+
+
+
+    walberla::SetupBlockForest sforest;
+    //sforest.addWorkloadMemorySUIDAssignmentFunction( uniformWorkloadAndMemoryAssignment );
+    sforest.init(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)), // blocks in each direction
+       val<1>(peParameter->isPeriodic), val<2>(peParameter->isPeriodic), val<3>(peParameter->isPeriodic));
+    sforest.balanceLoad(*loadBalancer.get(), loadBalancer->getNumberOfProcesses());
+    forest = std::shared_ptr< walberla::blockforest::BlockForest >( new walberla::blockforest::BlockForest( walberla::uint_c( loadBalancer->getRank() ), sforest) );
+
+     auto mpiManager = walberla::MPIManager::instance();
+     mpiManager->useWorldComm();
     if (!forest)
-        throw std::runtime_error("No PE BlockForest created ... ");
+       throw std::runtime_error("No PE BlockForest created ... ");
 }
 
 void PePhysicsEngineSolverAdapter::initalBlockData()
diff --git a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h
index a7b699252d397080dc9e1bfeceac5edbeb02337e..1f49f063285440ec6148ff38920d9af2824dee9f 100644
--- a/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h
+++ b/source/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h
@@ -17,6 +17,7 @@
 
 class PePhysicsEngineMaterialAdapter;
 class PhysicsEngineGeometryAdapter;
+class PeLoadBalancerAdapter;
 
 namespace walberla
 {
@@ -64,7 +65,7 @@ struct PeParameter
 class PePhysicsEngineSolverAdapter : public PhysicsEngineSolverAdapter
 {
 public:
-    PePhysicsEngineSolverAdapter(std::shared_ptr<PeParameter> peParameter);
+    PePhysicsEngineSolverAdapter(std::shared_ptr<PeParameter> peParameter, std::shared_ptr<PeLoadBalancerAdapter> loadBalancer);
     virtual ~PePhysicsEngineSolverAdapter() {}
 
     std::shared_ptr<PhysicsEngineGeometryAdapter> createPhysicsEngineGeometryAdapter(int id, const Vector3D& position, double radius, std::shared_ptr<PhysicsEngineMaterialAdapter> material) const override;
@@ -87,6 +88,7 @@ private:
 
 private:
     std::shared_ptr<PeParameter> peParameter;
+    std::shared_ptr<PeLoadBalancerAdapter> loadBalancer;
 
     std::shared_ptr<walberla::pe::BodyStorage> globalBodyStorage;
     std::shared_ptr< walberla::blockforest::BlockForest > forest;