From 72f8ad800e93bbb976f2546e1c331c08a86adcdb Mon Sep 17 00:00:00 2001
From: Christoph Rettinger <christoph.rettinger@fau.de>
Date: Wed, 17 Jan 2018 20:27:43 +0100
Subject: [PATCH] [API]: changed body mapping functions: removed most of them,
 added mapping-decider function, accomodated changes to test cases, added new
 mapping test case

---
 .../ForcesOnSphereNearPlaneInShearFlow.cpp    |    6 +-
 .../MotionSingleHeavySphere.cpp               |   16 +-
 src/pe_coupling/mapping/BodyMapping.h         |  120 +-
 .../momentum_exchange_method/BodyMapping.h    |  204 +--
 .../restoration/PDFReconstruction.h           |  139 +-
 .../utility/BodySelectorFunctions.h           |   49 +
 src/pe_coupling/utility/all.h                 |    1 +
 tests/pe_coupling/CMakeLists.txt              |    3 +
 .../BodyAtBlockBoarderCheck.cpp               |   53 +-
 .../BodyMappingTest.cpp                       | 1134 +++++++++++++++++
 .../DragForceSphereMEM.cpp                    |   10 +-
 .../DragForceSphereMEMRefinement.cpp          |    8 +-
 ...lobalBodyAsBoundaryMEMStaticRefinement.cpp |    3 +-
 .../LubricationCorrectionMEM.cpp              |   13 +-
 .../PeriodicParticleChannelMEM.cpp            |   29 +-
 .../SegreSilberbergMEM.cpp                    |   20 +-
 .../SettlingSphereMEM.cpp                     |    8 +-
 .../SettlingSphereMEMDynamicRefinement.cpp    |   15 +-
 .../SettlingSphereMEMStaticRefinement.cpp     |   10 +-
 .../momentum_exchange_method/SquirmerTest.cpp |    9 +-
 .../TaylorCouetteFlowMEM.cpp                  |    4 +-
 .../TorqueSphereMEM.cpp                       |   10 +-
 22 files changed, 1491 insertions(+), 373 deletions(-)
 create mode 100644 src/pe_coupling/utility/BodySelectorFunctions.h
 create mode 100644 tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp

diff --git a/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp b/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
index c14c38922..e10ef4cb9 100644
--- a/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
+++ b/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
@@ -617,14 +617,14 @@ int main( int argc, char **argv )
    // map planes into the LBM simulation -> act as no-slip boundaries
    // since we have refinement boundaries along the bottom plane, we use SBB here
    // (see test/pe_coupling/momentum_exchange_method/GlobalBodyAsBoundaryMEMStaticRefinement test for more infos)
-   pe_coupling::mapMovingGlobalBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, *globalBodyStorage, bodyFieldID, MO_SBB_Flag );
+   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_SBB_Flag, pe_coupling::selectGlobalBodies );
 
    // map movable sphere into the LBM simulation
    // the whole sphere resides on the same refinement level, so we can also use CLI instead of SBB
    if( useSBB )
-      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_SBB_Flag );
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_SBB_Flag, pe_coupling::selectRegularBodies );
    else
-      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag );
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_CLI_Flag, pe_coupling::selectRegularBodies );
 
 
    // initialize Couette velocity profile in whole domain
diff --git a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
index 0774932f3..90e2b3479 100644
--- a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
+++ b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
@@ -797,7 +797,7 @@ int main( int argc, char **argv )
    /////////////////
 
    // set up pe functionality
-   shared_ptr<pe::BodyStorage>  globalBodyStorage = make_shared<pe::BodyStorage>();
+   shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
    pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
    auto bodyStorageID  = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "pe Body Storage");
    auto ccdID          = blocks->addBlockData(pe::ccd::createHashGridsDataHandling( globalBodyStorage, bodyStorageID ), "CCD");
@@ -912,11 +912,11 @@ int main( int argc, char **argv )
 
       if( memVariant == MEMVariant::CLI )
       {
-         pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_CLI_Flag );
+         pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MEM_CLI_Flag );
       }else if ( memVariant == MEMVariant::MR ){
-         pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_MR_Flag );
+         pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MEM_MR_Flag );
       }else{
-         pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_BB_Flag );
+         pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MEM_BB_Flag );
       }
    }
 
@@ -1195,11 +1195,11 @@ int main( int argc, char **argv )
 
       // sweep for updating the pe body mapping into the LBM simulation
       if( memVariant == MEMVariant::CLI )
-         timeloop.add() << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MEM_CLI_Flag, FormerMEM_Flag ), "Body Mapping" );
+         timeloop.add() << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MEM_CLI_Flag, FormerMEM_Flag ), "Body Mapping" );
       else if ( memVariant == MEMVariant::MR )
-         timeloop.add() << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MEM_MR_Flag, FormerMEM_Flag ), "Body Mapping" );
+         timeloop.add() << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MEM_MR_Flag, FormerMEM_Flag ), "Body Mapping" );
       else
-         timeloop.add() << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MEM_BB_Flag, FormerMEM_Flag ), "Body Mapping" );
+         timeloop.add() << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MEM_BB_Flag, FormerMEM_Flag ), "Body Mapping" );
 
 
       // reconstruct missing PDFs
@@ -1208,7 +1208,7 @@ int main( int argc, char **argv )
       typedef pe_coupling::ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationFinder_T > Reconstructor_T;
       Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationFinder, true );
       timeloop.add() << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >
-                               ( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, reconstructor, FormerMEM_Flag, Fluid_Flag  ), "PDF Restore" );
+                               ( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMEM_Flag, Fluid_Flag ), "PDF Restore" );
 
       for( uint_t lbmcycle = 0; lbmcycle < numLBMSubCycles; ++lbmcycle ){
 
diff --git a/src/pe_coupling/mapping/BodyMapping.h b/src/pe_coupling/mapping/BodyMapping.h
index 2e70c1439..042c8dd23 100644
--- a/src/pe_coupling/mapping/BodyMapping.h
+++ b/src/pe_coupling/mapping/BodyMapping.h
@@ -28,25 +28,29 @@
 #include "core/debug/Debug.h"
 #include "domain_decomposition/StructuredBlockStorage.h"
 #include "field/FlagField.h"
-#include "field/iterators/IteratorMacros.h"
 
 #include "pe/rigidbody/RigidBody.h"
 #include "pe/rigidbody/BodyStorage.h"
 #include "pe/rigidbody/BodyIterators.h"
 
+#include "pe_coupling/utility/BodySelectorFunctions.h"
+
+#include <boost/function.hpp>
+
 namespace walberla {
 namespace pe_coupling {
 
 // general mapping functions for a given single body on a given single block
 template< typename BoundaryHandling_T >
-void mapBody( const pe::BodyID & body, IBlock & block, StructuredBlockStorage & blockStorage,
-              BoundaryHandling_T * boundaryHandlingPtr, const FlagUID & obstacle )
+void mapBody( pe::BodyID body, IBlock & block, StructuredBlockStorage & blockStorage,
+              const BlockDataID & boundaryHandlingID, const FlagUID & obstacle )
 {
    WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage.getBlockStorage()) );
-   WALBERLA_ASSERT_NOT_NULLPTR( boundaryHandlingPtr );
 
-   auto * flagField = boundaryHandlingPtr->getFlagField();
+   BoundaryHandling_T * boundaryHandling = block.getData< BoundaryHandling_T >( boundaryHandlingID );
+   auto * flagField = boundaryHandling->getFlagField();
 
+   WALBERLA_ASSERT_NOT_NULLPTR( boundaryHandling );
    WALBERLA_ASSERT_NOT_NULLPTR( flagField );
    WALBERLA_ASSERT( flagField->flagExists( obstacle ) );
 
@@ -71,7 +75,9 @@ void mapBody( const pe::BodyID & body, IBlock & block, StructuredBlockStorage &
          for( cell_idx_t x = cellBB.xMin(); x <= cellBB.xMax(); ++x )
          {
             if( body->containsPoint(cx,cy,cz) )
-               boundaryHandlingPtr->forceBoundary( obstacleFlag, x, y, z );
+            {
+               boundaryHandling->forceBoundary( obstacleFlag, x, y, z );
+            }
             cx += dx;
          }
          cy += dy;
@@ -80,110 +86,38 @@ void mapBody( const pe::BodyID & body, IBlock & block, StructuredBlockStorage &
    }
 }
 
-template< typename BoundaryHandling_T >
-void mapBody( const pe::BodyID & body, IBlock & block, StructuredBlockStorage & blockStorage,
-              const BlockDataID & boundaryHandlingID, const FlagUID & obstacle )
-{
-   WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage.getBlockStorage()) );
-
-   BoundaryHandling_T * boundaryHandling = block.getData< BoundaryHandling_T >( boundaryHandlingID );
-
-   mapBody(body, block, blockStorage, boundaryHandling, obstacle );
-}
-
 
 
-// mapping function to map all bodies from the body storage - with certain properties - onto all blocks
+/*!\brief Mapping function to map all bodies - with certain properties - onto all blocks with the 'obstacle' flag
+ *
+ * All bodies (from bodyStorage and globalBodyStorage) are iterated and mapped to all blocks.
+ * Cells that are inside the bodies are set to 'obstacle'.
+ *
+ * Whether or not a body is mapped depends on the return value of the 'mappingBodySelectorFct'.
+ */
 template< typename BoundaryHandling_T >
 void mapBodies( StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
-                const BlockDataID & bodyStorageID, const FlagUID & obstacle,
-                const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
+                const BlockDataID & bodyStorageID, pe::BodyStorage & globalBodyStorage,
+                const FlagUID & obstacle,
+                const boost::function<bool(pe::BodyID)> & mappingBodySelectorFct = selectAllBodies )
 {
    for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
    {
+
       for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
       {
-         WALBERLA_UNUSED(moBodiesOnly); // undo when other coupling algorithms are available
-
-         if( fixedBodiesOnly && bodyIt->isFixed() )
-            continue;
-
-         mapBody<BoundaryHandling_T>( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle );
+         if( mappingBodySelectorFct(*bodyIt))
+            mapBody<BoundaryHandling_T>( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle );
       }
-   }
-}
-
 
-// mapping function to map all global bodies - with certain properties - onto all blocks
-template< typename BoundaryHandling_T >
-void mapGlobalBodies( StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
-                      pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle,
-                      const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
-{
-   for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
-   {
       for( auto bodyIt = globalBodyStorage.begin(); bodyIt != globalBodyStorage.end(); ++bodyIt )
       {
-         WALBERLA_UNUSED(moBodiesOnly); // undo when other coupling algorithms are available
-
-         if( fixedBodiesOnly && bodyIt->isFixed() )
-            continue;
-
-         mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle );
+         if( mappingBodySelectorFct(*bodyIt))
+            mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle );
       }
    }
 }
 
-// mapping function to map a given single global body onto all blocks
-template< typename BoundaryHandling_T >
-void mapGlobalBody( const id_t globalBodySystemID,
-                    StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
-                    pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle )
-{
-   for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
-   {
-      auto bodyIt = globalBodyStorage.find( globalBodySystemID );
-      if( bodyIt != globalBodyStorage.end() )
-      {
-         mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle );
-      }
-   }
-}
-
-
-// mapping function to map all global bodies - with certain properties - onto a given single block
-template< typename BoundaryHandling_T >
-void mapGlobalBodiesOnBlock( IBlock & block,
-                             StructuredBlockStorage & blockStorage, BoundaryHandling_T * boundaryHandlingPtr,
-                             pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle,
-                             const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
-{
-   for( auto bodyIt = globalBodyStorage.begin(); bodyIt != globalBodyStorage.end(); ++bodyIt)
-   {
-      WALBERLA_UNUSED(moBodiesOnly); // undo when other coupling algorithms are available
-
-      if( fixedBodiesOnly && bodyIt->isFixed() )
-         continue;
-
-      mapBody< BoundaryHandling_T >( *bodyIt, block, blockStorage, boundaryHandlingPtr, obstacle );
-   }
-}
-
-
-// mapping function to map a given single global body onto a given single block
-template< typename BoundaryHandling_T >
-void mapGlobalBodyOnBlock( const id_t globalBodySystemID, IBlock & block,
-                           StructuredBlockStorage & blockStorage, BoundaryHandling_T * boundaryHandlingPtr,
-                           pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle )
-{
-   auto bodyIt = globalBodyStorage.find( globalBodySystemID );
-   if( bodyIt != globalBodyStorage.end() )
-   {
-      mapBody< BoundaryHandling_T >( *bodyIt, block, blockStorage, boundaryHandlingPtr, obstacle );
-   }
-
-}
-
 
 } // namespace pe_coupling
 } // namespace walberla
diff --git a/src/pe_coupling/momentum_exchange_method/BodyMapping.h b/src/pe_coupling/momentum_exchange_method/BodyMapping.h
index 77d5cf401..3b4b7ff57 100644
--- a/src/pe_coupling/momentum_exchange_method/BodyMapping.h
+++ b/src/pe_coupling/momentum_exchange_method/BodyMapping.h
@@ -23,6 +23,7 @@
 
 #pragma once
 
+#include "core/ConcatIterator.h"
 #include "core/debug/Debug.h"
 #include "core/cell/Cell.h"
 #include "domain_decomposition/StructuredBlockStorage.h"
@@ -32,88 +33,97 @@
 #include "pe/rigidbody/BodyIterators.h"
 
 #include "pe_coupling/mapping/BodyBBMapping.h"
+#include "pe_coupling/utility/BodySelectorFunctions.h"
+
+#include <boost/function.hpp>
 
 namespace walberla {
 namespace pe_coupling {
 
+/*!\brief Maps the moving bodies into the simulation domain and updates the mapping
+ *
+ * Cells that are inside a body, will be marked with the 'obstacle' flag.
+ * 'Inside' means that the cell center is contained in the body.
+ * Thus, the body has to provide a valid containsPoint(x,y,z) function.
+ *
+ * Cells that in the last time step were inside the body but are now outside of it, i.e. the body has moved,
+ * will be marked with the 'formerObstacle' flag.
+ * The 'formerObstacle' flag is used in a second step by the PDFReconstruction class (see PDFReconstruction.h) to
+ * re-initialize the missing PDFs. Afterwards, the 'formerObstacle' flag is removed and the 'fluid' flag is set.
+ *
+ * It is not required that the mapping has been initialized with one of the free functions from below.
+ *
+ * The 'mappingBodySelectorFct' can be used to decide which bodies should be mapped or not.
+ */
 template< typename BoundaryHandling_T >
 class BodyMapping
 {
 public:
 
    typedef typename BoundaryHandling_T::FlagField FlagField_T;
-   typedef typename FlagField_T::iterator         FlagFieldIterator;
    typedef typename BoundaryHandling_T::flag_t    flag_t;
    typedef Field< pe::BodyID, 1 >                 BodyField_T;
 
-   inline BodyMapping( const shared_ptr<StructuredBlockStorage> & blockStorage,
-                       const BlockDataID & boundaryHandlingID,
-                       const BlockDataID & bodyStorageID, const BlockDataID & bodyFieldID,
-                       const FlagUID & obstacle, const FlagUID & formerObstacle );
-
-   void operator()( IBlock * const block ) const;
-
-protected:
-
-   shared_ptr<StructuredBlockStorage> blockStorage_;
-
-   const BlockDataID boundaryHandlingID_;
-   const BlockDataID bodyStorageID_;
-   const BlockDataID bodyFieldID_;
-
-   const FlagUID obstacle_;
-   const FlagUID formerObstacle_;
-
-}; // class BodyMapping
-
-
-
-template< typename BoundaryHandling_T >
-inline BodyMapping< BoundaryHandling_T >::BodyMapping( const shared_ptr<StructuredBlockStorage> & blockStorage,
-                                                       const BlockDataID & boundaryHandlingID,
-                                                       const BlockDataID & bodyStorageID, const BlockDataID & bodyFieldID,
-                                                       const FlagUID & obstacle, const FlagUID & formerObstacle ) :
-
-   blockStorage_( blockStorage ), boundaryHandlingID_( boundaryHandlingID ), bodyStorageID_(bodyStorageID), bodyFieldID_( bodyFieldID ),
-   obstacle_( obstacle ), formerObstacle_( formerObstacle )
-{}
+   BodyMapping( const shared_ptr<StructuredBlockStorage> & blockStorage,
+                const BlockDataID & boundaryHandlingID,
+                const BlockDataID & bodyStorageID,
+                const shared_ptr<pe::BodyStorage> & globalBodyStorage,
+                const BlockDataID & bodyFieldID,
+                const FlagUID & obstacle, const FlagUID & formerObstacle,
+                const boost::function<bool(pe::BodyID)> & mappingBodySelectorFct = selectRegularBodies )
+   : blockStorage_( blockStorage ), boundaryHandlingID_( boundaryHandlingID ),
+     bodyStorageID_(bodyStorageID), globalBodyStorage_( globalBodyStorage ), bodyFieldID_( bodyFieldID ),
+     obstacle_( obstacle ), formerObstacle_( formerObstacle ), mappingBodySelectorFct_( mappingBodySelectorFct )
+   {}
+
+   void operator()( IBlock * const block )
+   {
+      WALBERLA_ASSERT_NOT_NULLPTR( block );
 
+      BoundaryHandling_T * boundaryHandling = block->getData< BoundaryHandling_T >( boundaryHandlingID_ );
+      FlagField_T *        flagField        = boundaryHandling->getFlagField();
+      BodyField_T *        bodyField        = block->getData< BodyField_T >( bodyFieldID_ );
 
+      WALBERLA_ASSERT_NOT_NULLPTR( boundaryHandling );
+      WALBERLA_ASSERT_NOT_NULLPTR( flagField );
+      WALBERLA_ASSERT_NOT_NULLPTR( bodyField );
 
-template< typename BoundaryHandling_T >
-void BodyMapping< BoundaryHandling_T >::operator()( IBlock * const block ) const
-{
-   WALBERLA_ASSERT_NOT_NULLPTR( block );
+      WALBERLA_ASSERT_EQUAL( flagField->xyzSize(), bodyField->xyzSize() );
 
-   BoundaryHandling_T * boundaryHandling = block->getData< BoundaryHandling_T >( boundaryHandlingID_ );
-   FlagField_T *        flagField        = boundaryHandling->getFlagField();
-   BodyField_T *        bodyField        = block->getData< BodyField_T >( bodyFieldID_ );
+      WALBERLA_ASSERT( flagField->flagExists( obstacle_ ) );
 
-   WALBERLA_ASSERT_NOT_NULLPTR( boundaryHandling );
-   WALBERLA_ASSERT_NOT_NULLPTR( flagField );
-   WALBERLA_ASSERT_NOT_NULLPTR( bodyField );
+      const flag_t       obstacle = flagField->getFlag( obstacle_ );
+      const flag_t formerObstacle = flagField->flagExists( formerObstacle_ ) ? flagField->getFlag( formerObstacle_ ) :
+                                    flagField->registerFlag( formerObstacle_ );
 
-   WALBERLA_ASSERT_EQUAL( flagField->xyzSize(), bodyField->xyzSize() );
+      const real_t dx = blockStorage_->dx( blockStorage_->getLevel(*block) );
+      const real_t dy = blockStorage_->dy( blockStorage_->getLevel(*block) );
+      const real_t dz = blockStorage_->dz( blockStorage_->getLevel(*block) );
 
-   WALBERLA_ASSERT( flagField->flagExists( obstacle_ ) );
+      for( auto bodyIt = pe::BodyIterator::begin(*block, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+      {
+         if( mappingBodySelectorFct_(*bodyIt) )
+            mapBodyAndUpdateMapping(*bodyIt, block, boundaryHandling, flagField , bodyField, obstacle, formerObstacle, dx, dy, dz);
+      }
+      for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt)
+      {
+         if( mappingBodySelectorFct_(*bodyIt))
+            mapBodyAndUpdateMapping(*bodyIt, block, boundaryHandling, flagField , bodyField, obstacle, formerObstacle, dx, dy, dz);
+      }
+   }
 
-   const flag_t       obstacle = flagField->getFlag( obstacle_ );
-   const flag_t formerObstacle = flagField->flagExists( formerObstacle_ ) ? flagField->getFlag( formerObstacle_ ) :
-                                                                            flagField->registerFlag( formerObstacle_ );
+private:
 
-   for( auto bodyIt = pe::BodyIterator::begin(*block, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+   void mapBodyAndUpdateMapping(pe::BodyID body, IBlock * const block,
+                                BoundaryHandling_T * boundaryHandling, FlagField_T * flagField, BodyField_T * bodyField,
+                                const flag_t & obstacle, const flag_t & formerObstacle,
+                                real_t dx, real_t dy, real_t dz)
    {
-      if( bodyIt->isFixed() /*|| !isMOBody( *bodyIt )*/ )
-         continue;
-
       // policy: every body manages only its own flags
 
-      CellInterval cellBB = getCellBB( *bodyIt, *block, *blockStorage_, flagField->nrOfGhostLayers() );
+      CellInterval cellBB = getCellBB( body, *block, *blockStorage_, flagField->nrOfGhostLayers() );
 
       Vector3<real_t> startCellCenter = blockStorage_->getBlockLocalCellCenter( *block, cellBB.min() );
-      const real_t dx = blockStorage_->dx( blockStorage_->getLevel(*block) );
-      const real_t dy = blockStorage_->dy( blockStorage_->getLevel(*block) );
-      const real_t dz = blockStorage_->dz( blockStorage_->getLevel(*block) );
 
       real_t cz = startCellCenter[2];
       for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
@@ -127,34 +137,40 @@ void BodyMapping< BoundaryHandling_T >::operator()( IBlock * const block ) const
 
                flag_t & cellFlagPtr = flagField->get(x,y,z);
 
-               if( bodyIt->containsPoint(cx,cy,cz) )
+               if( body->containsPoint(cx,cy,cz) )
                {
+                  // cell is inside body
                   if( !isFlagSet( cellFlagPtr, obstacle ) )
                   {
-
+                     // cell is not yet an obstacle cell
                      if( isFlagSet( cellFlagPtr, formerObstacle ) )
                      {
+                        // cell was marked as former obstacle, e.g. by another body that has moved away
                         boundaryHandling->setBoundary( obstacle, x, y, z );
                         removeFlag( cellFlagPtr, formerObstacle );
                      }
                      else
                      {
+                        // set obstacle flag
                         boundaryHandling->forceBoundary( obstacle, x, y, z );
                      }
                   }
-
-                  (*bodyField)(x,y,z) = *bodyIt;
+                  // let pointer from body field point to this body
+                  (*bodyField)(x,y,z) = body;
                }
                else
                {
-                  if( isFlagSet( cellFlagPtr, obstacle ) && ((*bodyField)(x, y, z) == *bodyIt) )
+                  // cell is outside body
+                  if( isFlagSet( cellFlagPtr, obstacle ) && ((*bodyField)(x, y, z) == body) )
                   {
+                     // cell was previously occupied by this body
                      boundaryHandling->removeBoundary( obstacle, x, y, z );
                      addFlag( cellFlagPtr, formerObstacle );
                      // entry at (*bodyField)(x,y,z) should still point to the previous body.
                      // If during initialization the overlap between neighboring blocks
                      // was chosen correctly/large enough, the body should still be on this block.
                      // The body information is needed in the PDF restoration step.
+                     // There, the flag will be removed and replaced by a domain flag after restoration.
                   }
                }
 
@@ -165,25 +181,36 @@ void BodyMapping< BoundaryHandling_T >::operator()( IBlock * const block ) const
          cz += dz;
       }
    }
-}
 
+   shared_ptr<StructuredBlockStorage> blockStorage_;
+
+   const BlockDataID boundaryHandlingID_;
+   const BlockDataID bodyStorageID_;
+   shared_ptr<pe::BodyStorage> globalBodyStorage_;
+   const BlockDataID bodyFieldID_;
+
+   const FlagUID obstacle_;
+   const FlagUID formerObstacle_;
+
+   boost::function<bool(pe::BodyID)> mappingBodySelectorFct_;
+
+}; // class BodyMapping
 
 
 ////////////////////////////
 // FREE MAPPING FUNCTIONS //
 ////////////////////////////
 
+
+// general mapping functions for a given single (moving) body on a given single block
 template< typename BoundaryHandling_T >
-void mapMovingBody( const pe::BodyID body, IBlock & block, StructuredBlockStorage & blockStorage,
+void mapMovingBody( pe::BodyID body, IBlock & block, StructuredBlockStorage & blockStorage,
                     const BlockDataID & boundaryHandlingID, const BlockDataID & bodyFieldID, const FlagUID & obstacle )
 {
    typedef Field< pe::BodyID, 1 > BodyField_T;
 
    WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage.getBlockStorage()) );
 
-   if( body->isFixed() /*|| !body->isFinite()*/ )
-      return;
-
    BoundaryHandling_T * boundaryHandling = block.getData< BoundaryHandling_T >( boundaryHandlingID );
    auto *               flagField        = boundaryHandling->getFlagField();
    BodyField_T *        bodyField        = block.getData< BodyField_T >( bodyFieldID );
@@ -229,43 +256,30 @@ void mapMovingBody( const pe::BodyID body, IBlock & block, StructuredBlockStorag
    }
 }
 
-
-
-template< typename BoundaryHandling_T >
-void mapMovingBodies( StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID, const BlockDataID & bodyStorageID,
-                      const BlockDataID & bodyFieldID, const FlagUID & obstacle )
-{
-   for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
-   {
-       for (auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt)
-           mapMovingBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, bodyFieldID, obstacle );
-   }
-}
-
+/*!\brief Mapping function to map all bodies - with certain properties - onto all blocks to the "moving obstacle" boundary condition, that requires the additional bodyField
+ *
+ * All bodies (from bodyStorage and globalBodyStorage) are iterated and mapped to all blocks.
+ * Cells that are inside the bodies are set to 'obstacle' and a pointer to the body is stored in the bodyField.
+ *
+ * Whether or not a body is mapped depends on the return value of the 'mappingBodySelectorFct'.
+ */
 template< typename BoundaryHandling_T >
-void mapMovingGlobalBodies( StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
-                            pe::BodyStorage & globalBodyStorage, const BlockDataID & bodyFieldID, const FlagUID & obstacle )
+void mapMovingBodies( StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
+                      const BlockDataID & bodyStorageID, pe::BodyStorage & globalBodyStorage,
+                      const BlockDataID & bodyFieldID, const FlagUID & obstacle,
+                      const boost::function<bool(pe::BodyID)> & mappingBodySelectorFct = selectAllBodies )
 {
    for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
    {
-      for( auto bodyIt = globalBodyStorage.begin(); bodyIt != globalBodyStorage.end(); ++bodyIt)
+      for (auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt)
       {
-         mapMovingBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, bodyFieldID, obstacle );
+         if( mappingBodySelectorFct(*bodyIt))
+            mapMovingBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, bodyFieldID, obstacle );
       }
-   }
-}
-
-template< typename BoundaryHandling_T >
-void mapMovingGlobalBody( const id_t globalBodySystemID,
-                          StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
-                          pe::BodyStorage & globalBodyStorage, const BlockDataID & bodyFieldID, const FlagUID & obstacle )
-{
-   for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
-   {
-      auto bodyIt = globalBodyStorage.find( globalBodySystemID );
-      if( bodyIt != globalBodyStorage.end() )
+      for( auto bodyIt = globalBodyStorage.begin(); bodyIt != globalBodyStorage.end(); ++bodyIt)
       {
-         mapMovingBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, bodyFieldID, obstacle );
+         if( mappingBodySelectorFct(*bodyIt))
+            mapMovingBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, bodyFieldID, obstacle );
       }
    }
 }
diff --git a/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h b/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h
index f7ca0d80a..ef8289da0 100644
--- a/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h
+++ b/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h
@@ -24,6 +24,7 @@
 
 #include "pe/rigidbody/BodyIterators.h"
 #include "pe_coupling/mapping/BodyBBMapping.h"
+#include "pe_coupling/utility/BodySelectorFunctions.h"
 
 #include "core/debug/Debug.h"
 #include "domain_decomposition/StructuredBlockStorage.h"
@@ -32,23 +33,25 @@
 #include "field/iterators/IteratorMacros.h"
 #include "Reconstructor.h"
 
+#include <boost/function.hpp>
 
 namespace walberla {
 namespace pe_coupling {
 
 
 //**************************************************************************************************************************************
-/*!
-*   \brief Class to manage the reconstruction of PDFs that is needed when cells are becoming uncovered by moving obstacles.
+/*!\brief Class to manage the reconstruction of PDFs that is needed when cells are becoming uncovered by moving obstacles.
 *
-*   Due to the explicit mapping of bodies into the fluid domain via flags, the PDFs of cells that turned from obstacle to fluid
-*   are missing and must be reconstructed in order to continue with the simulation.
-*   This class is to be used as a sweep in a LBM simulation with moving obstacles and calls for each cell that is tagged as
-*   'formerObstacle' the specified reconstructor (see Reconstructor.h for the available variants).
-*   After the successful reconstruction of all PDFs, the flags are updated to 'fluid'.
-*   For small obstacle fractions, an optimized variant is available that only looks for 'formerObstacle' cells in the vicinity
-*   of available bodies. It is activated via the 'optimizeForSmallObstacleFraction' argument in the constructor.
+* Due to the explicit mapping of bodies into the fluid domain via flags, the PDFs of cells that turned from obstacle to fluid
+* are missing and must be reconstructed in order to continue with the simulation.
+* This class is to be used as a sweep in a LBM simulation with moving obstacles and calls for each cell that is tagged as
+* 'formerObstacle' the specified reconstructor (see Reconstructor.h for the available variants).
+* After the successful reconstruction of all PDFs, the flags are updated to 'fluid'.
+* For small obstacle fractions, an optimized variant is available that only looks for 'formerObstacle' cells in the vicinity
+* of available bodies. It is activated via the 'optimizeForSmallObstacleFraction' argument in the constructor.
 *
+* The 'movingBodySelectorFct' can be used to decide which bodies should be check for reconstruction
+* (only used when 'optimizeForSmallObstacleFraction' is chosen).
 */
 //**************************************************************************************************************************************
 
@@ -58,29 +61,65 @@ class PDFReconstruction
 public:
 
    typedef typename BoundaryHandling_T::FlagField FlagField_T;
-   typedef typename FlagField_T::iterator         FlagFieldIterator;
    typedef typename BoundaryHandling_T::flag_t    flag_t;
    typedef Field< pe::BodyID, 1 >                 BodyField_T;
 
    inline PDFReconstruction( const shared_ptr<StructuredBlockStorage> & blockStorage,
                              const BlockDataID & boundaryHandlingID,
-                             const BlockDataID & bodyStorageID, const BlockDataID & bodyFieldID,
+                             const BlockDataID & bodyStorageID,
+                             const shared_ptr<pe::BodyStorage> & globalBodyStorage,
+                             const BlockDataID & bodyFieldID,
                              const Reconstructer_T & reconstructor,
                              const FlagUID & formerObstacle, const FlagUID & fluid,
+                             const boost::function<bool(pe::BodyID)> & movingBodySelectorFct = selectRegularBodies,
                              const bool optimizeForSmallObstacleFraction = false ) :
-      blockStorage_( blockStorage ), boundaryHandlingID_( boundaryHandlingID ), bodyStorageID_(bodyStorageID), bodyFieldID_( bodyFieldID ),
+      blockStorage_( blockStorage ), boundaryHandlingID_( boundaryHandlingID ), bodyStorageID_(bodyStorageID),
+      globalBodyStorage_( globalBodyStorage ), bodyFieldID_( bodyFieldID ),
       reconstructor_ ( reconstructor ), formerObstacle_( formerObstacle ), fluid_( fluid ),
+      movingBodySelectorFct_( movingBodySelectorFct ),
       optimizeForSmallObstacleFraction_( optimizeForSmallObstacleFraction )
    {}
 
    void operator()( IBlock * const block );
 
-protected:
+private:
+
+   void reconstructPDFsInCells( const CellInterval & cells, IBlock * const block,
+                                FlagField_T * flagField, const flag_t & formerObstacle )
+   {
+      for( cell_idx_t z = cells.zMin(); z <= cells.zMax(); ++z ) {
+         for (cell_idx_t y = cells.yMin(); y <= cells.yMax(); ++y) {
+            for (cell_idx_t x = cells.xMin(); x <= cells.xMax(); ++x) {
+               if (isFlagSet(flagField->get(x,y,z), formerObstacle)) {
+                  reconstructor_(x,y,z,block);
+               }
+            }
+         }
+      }
+   }
+
+   void updateFlags( const CellInterval & cells,
+                     BoundaryHandling_T * boundaryHandling, FlagField_T * flagField, BodyField_T * bodyField,
+                     const flag_t & formerObstacle, const flag_t & fluid)
+   {
+      for( cell_idx_t z = cells.zMin(); z <= cells.zMax(); ++z ) {
+         for (cell_idx_t y = cells.yMin(); y <= cells.yMax(); ++y) {
+            for (cell_idx_t x = cells.xMin(); x <= cells.xMax(); ++x) {
+               if (isFlagSet(flagField->get(x,y,z), formerObstacle)) {
+                  boundaryHandling->setDomain( fluid, x, y, z );
+                  removeFlag( flagField->get(x,y,z), formerObstacle );
+                  (*bodyField)(x,y,z) = NULL;
+               }
+            }
+         }
+      }
+   }
 
    shared_ptr<StructuredBlockStorage> blockStorage_;
 
    const BlockDataID boundaryHandlingID_;
    const BlockDataID bodyStorageID_;
+   shared_ptr<pe::BodyStorage> globalBodyStorage_;
    const BlockDataID bodyFieldID_;
 
    Reconstructer_T reconstructor_;
@@ -88,6 +127,8 @@ protected:
    const FlagUID formerObstacle_;
    const FlagUID fluid_;
 
+   boost::function<bool(pe::BodyID)> movingBodySelectorFct_;
+
    const bool optimizeForSmallObstacleFraction_;
 
 };
@@ -115,39 +156,32 @@ void PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructer_T >
    const flag_t formerObstacle = flagField->getFlag( formerObstacle_ );
    const flag_t fluid          = flagField->getFlag( fluid_ );
 
-   // reconstruct all missing PDFs (only inside the domain)
+   // reconstruct all missing PDFs (only inside the domain, ghost layer values get communicated)
    if( optimizeForSmallObstacleFraction_ )
    {
       const uint_t numberOfGhostLayersToInclude = uint_t(0);
 
       for( auto bodyIt = pe::BodyIterator::begin(*block, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
       {
-         if( bodyIt->isFixed() )
+         if( !movingBodySelectorFct_(*bodyIt) )
             continue;
 
          CellInterval cellBB = getCellBB( *bodyIt, *block, *blockStorage_, numberOfGhostLayersToInclude );
+         reconstructPDFsInCells(cellBB, block, flagField, formerObstacle );
+      }
+      for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
+      {
+         if( !movingBodySelectorFct_(*bodyIt) )
+            continue;
 
-         for( auto cellIt = cellBB.begin(); cellIt != cellBB.end(); ++cellIt )
-         {
-            Cell cell( *cellIt );
-            if( isFlagSet( flagField->get(cell), formerObstacle ) )
-            {
-               reconstructor_( cell.x(), cell.y(), cell.z(), block );
-            }
-         }
+         CellInterval cellBB = getCellBB( *bodyIt, *block, *blockStorage_, numberOfGhostLayersToInclude );
+         reconstructPDFsInCells(cellBB, block, flagField, formerObstacle );
       }
    }
    else
    {
-      auto xyzFieldSize = flagField->xyzSize();
-      for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
-      {
-         Cell cell( *fieldIt );
-         if( isFlagSet( flagField->get(cell), formerObstacle ) )
-         {
-            reconstructor_( cell.x(), cell.y(), cell.z(), block );
-         }
-      }
+      CellInterval cells = flagField->xyzSize();
+      reconstructPDFsInCells(cells, block, flagField, formerObstacle );
    }
 
    // update the flags from formerObstacle to fluid (inside domain & in ghost layers)
@@ -157,46 +191,25 @@ void PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructer_T >
 
       for( auto bodyIt = pe::BodyIterator::begin(*block, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
       {
-         if( bodyIt->isFixed() )
+         if( !movingBodySelectorFct_(*bodyIt) )
             continue;
 
          CellInterval cellBB = getCellBB( *bodyIt, *block, *blockStorage_, numberOfGhostLayersToInclude );
+         updateFlags(cellBB, boundaryHandling, flagField, bodyField, formerObstacle, fluid);
+      }
+      for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
+      {
+         if( !movingBodySelectorFct_(*bodyIt) )
+            continue;
 
-         for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
-         {
-            for( cell_idx_t y = cellBB.yMin(); y <= cellBB.yMax(); ++y )
-            {
-               for( cell_idx_t x = cellBB.xMin(); x <= cellBB.xMax(); ++x )
-               {
-                  if( isFlagSet( flagField->get(x,y,z), formerObstacle ) )
-                  {
-                     boundaryHandling->setDomain( fluid, x, y, z );
-                     removeFlag( flagField->get(x,y,z), formerObstacle );
-                     (*bodyField)(x,y,z) = NULL;
-                  }
-               }
-            }
-         }
+         CellInterval cellBB = getCellBB( *bodyIt, *block, *blockStorage_, numberOfGhostLayersToInclude );
+         updateFlags(cellBB, boundaryHandling, flagField, bodyField, formerObstacle, fluid);
       }
    }
    else
    {
       CellInterval cells = flagField->xyzSizeWithGhostLayer();
-      for( cell_idx_t z = cells.zMin(); z <= cells.zMax(); ++z )
-      {
-         for( cell_idx_t y = cells.yMin(); y <= cells.yMax(); ++y )
-         {
-            for( cell_idx_t x = cells.xMin(); x <= cells.xMax(); ++x )
-            {
-               if( isFlagSet( flagField->get(x,y,z), formerObstacle ) )
-               {
-                  boundaryHandling->setDomain( fluid, x, y, z );
-                  removeFlag( flagField->get(x,y,z), formerObstacle );
-                  (*bodyField)(x,y,z) = NULL;
-               }
-            }
-         }
-      }
+      updateFlags(cells, boundaryHandling, flagField, bodyField, formerObstacle, fluid);
    }
 }
 
diff --git a/src/pe_coupling/utility/BodySelectorFunctions.h b/src/pe_coupling/utility/BodySelectorFunctions.h
new file mode 100644
index 000000000..75bc5c93b
--- /dev/null
+++ b/src/pe_coupling/utility/BodySelectorFunctions.h
@@ -0,0 +1,49 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file BodySelectorFunctions.h
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "pe/Types.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+bool selectAllBodies(pe::BodyID /*bodyID*/)
+{
+   return true;
+}
+
+bool selectRegularBodies(pe::BodyID bodyID)
+{
+   return !bodyID->hasInfiniteMass() && !bodyID->isGlobal();
+}
+
+bool selectFixedBodies(pe::BodyID bodyID)
+{
+   return bodyID->hasInfiniteMass() && !bodyID->isGlobal();
+}
+
+bool selectGlobalBodies(pe::BodyID bodyID)
+{
+   return bodyID->isGlobal();
+}
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/utility/all.h b/src/pe_coupling/utility/all.h
index 9bb040405..d72c034b9 100644
--- a/src/pe_coupling/utility/all.h
+++ b/src/pe_coupling/utility/all.h
@@ -23,6 +23,7 @@
 #pragma once
 
 #include "BodiesForceTorqueContainer.h"
+#include "BodySelectorFunctions.h"
 #include "ForceOnBodiesAdder.h"
 #include "ForceTorqueOnBodiesResetter.h"
 #include "ForceTorqueOnBodiesScaler.h"
diff --git a/tests/pe_coupling/CMakeLists.txt b/tests/pe_coupling/CMakeLists.txt
index b675cc52d..1803b3090 100644
--- a/tests/pe_coupling/CMakeLists.txt
+++ b/tests/pe_coupling/CMakeLists.txt
@@ -72,6 +72,9 @@ waLBerla_execute_test( NAME BodyAtBlockBoarderCheckTest1        COMMAND $<TARGET
 waLBerla_execute_test( NAME BodyAtBlockBoarderCheckTest2        COMMAND $<TARGET_FILE:BodyAtBlockBoarderCheck> 2         PROCESSES 8 LABELS longrun )
 waLBerla_execute_test( NAME BodyAtBlockBoarderCheckTest3        COMMAND $<TARGET_FILE:BodyAtBlockBoarderCheck> 3         PROCESSES 8 LABELS verylongrun )
 
+waLBerla_compile_test( FILES momentum_exchange_method/BodyMappingTest.cpp DEPENDS blockforest pe timeloop )
+waLBerla_execute_test( NAME BodyMappingTest COMMAND $<TARGET_FILE:BodyMappingTest> PROCESSES 1 )
+
 waLBerla_compile_test( FILES momentum_exchange_method/DragForceSphereMEM.cpp DEPENDS blockforest pe timeloop )
 waLBerla_execute_test( NAME DragForceSphereMEMFuncTest        COMMAND $<TARGET_FILE:DragForceSphereMEM> --funcTest          PROCESSES 1 )
 waLBerla_execute_test( NAME DragForceSphereMEMSingleTest      COMMAND $<TARGET_FILE:DragForceSphereMEM>                     PROCESSES 1 LABELS longrun     CONFIGURATIONS Release RelWithDbgInfo )
diff --git a/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp b/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
index 16e0afc65..e47a1d4cb 100644
--- a/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
@@ -223,39 +223,6 @@ private:
       const pe::Vec3 referenceVelocity_;
 };
 
-class RefinementFunctorWrapper
-{
-public:
-   RefinementFunctorWrapper( boost::function< void () > fct )
-   : fct_( fct )
-   {
-   }
-
-   void operator()(const uint_t /*level*/, const uint_t /*executionCounter*/ )
-   {
-      fct_();
-   }
-
-private:
-   boost::function< void () > fct_;
-};
-
-class RefinementSweepWrapper
-{
-public:
-   RefinementSweepWrapper( boost::function< void ( IBlock * ) > fct )
-   : fct_( fct )
-   {
-   }
-
-   void operator()( IBlock * block, const uint_t /*level*/, const uint_t /*executionCounter*/ )
-   {
-      fct_( block );
-   }
-
-private:
-   boost::function< void ( IBlock * ) > fct_;
-};
 
 /*
  * This test case evaluates the correctness of the pe and the pe-coupling at the boundary between two blocks.
@@ -354,7 +321,7 @@ int main( int argc, char **argv )
                                     MyBoundaryHandling( flagFieldID, pdfFieldID, bodyFieldID ), "boundary handling" );
 
    // initially map pe bodies into the LBM simulation
-   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MO_Flag );
+   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID,  MO_Flag );
 
    ///////////////
    // TIME LOOP //
@@ -372,25 +339,25 @@ int main( int argc, char **argv )
    BodyEvaluator evalBodies( blocks, bodyStorageID, referenceVelocity );
 
    // add body evaluator
-   refinementTimestep->addPostCollideVoidFunction( RefinementFunctorWrapper( evalBodies ), "Body Evaluation", finestLevel );
+   refinementTimestep->addPostStreamVoidFunction( lbm::refinement::FunctorWrapper( evalBodies ), "Body Evaluation", finestLevel );
 
    // add pe timestep
    const real_t dtPE = real_c(1) / real_c( uint_c(1) << finestLevel );
-   //refinementTimestep->addPostCollideVoidFunction( RefinementFunctorWrapper( pe_coupling::TimeStep( blocks->getBlockStorage(), bodyStorageID, HCSITSWrapper( cr ), syncCall, dtPE, 1 ) ),
-   //                                                "pe Time Step", finestLevel );
 
-   refinementTimestep->addPostCollideVoidFunction( RefinementFunctorWrapper( pe_coupling::TimeStep( blocks, bodyStorageID, cr, syncCall, dtPE, 1 ) ),
-                                                   "pe Time Step", finestLevel );
+   refinementTimestep->addPostStreamVoidFunction( lbm::refinement::FunctorWrapper( pe_coupling::TimeStep( blocks, bodyStorageID, cr, syncCall, dtPE, 1 ) ),
+                                                  "pe Time Step", finestLevel );
 
    // add sweep for updating the pe body mapping into the LBM simulation
-   refinementTimestep->addPostCollideBlockFunction( RefinementSweepWrapper( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MO_Flag, FormerMO_Flag ) ),
-                                                    "Body Mapping", finestLevel );
+   refinementTimestep->addPostStreamVoidFunction( lbm::refinement::SweepAsFunctorWrapper( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID,
+                                                  bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag ), blocks ),
+                                                  "Body Mapping", finestLevel );
 
    // add sweep for restoring PDFs in cells previously occupied by pe bodies
    typedef pe_coupling::EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T > Reconstructor_T;
    Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID );
-   refinementTimestep->addPostCollideBlockFunction( RefinementSweepWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag  ) ),
-                                                    "PDF Restore", finestLevel );
+   refinementTimestep->addPostStreamVoidFunction( lbm::refinement::SweepAsFunctorWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks,
+                                                  boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag  ), blocks ),
+                                                  "PDF Restore", finestLevel );
 
    timeloop.addFuncBeforeTimeStep( makeSharedFunctor( refinementTimestep ), "LBM refinement time step" );
 
diff --git a/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp b/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp
new file mode 100644
index 000000000..0c527317c
--- /dev/null
+++ b/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp
@@ -0,0 +1,1134 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file BodyMappingTest.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "boundary/all.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/debug/Debug.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/math/all.h"
+#include "core/mpi/MPIManager.h"
+
+#include "field/AddToStorage.h"
+
+#include "lbm/boundary/all.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/lattice_model/D3Q19.h"
+
+#include "pe/basic.h"
+#include "pe/utility/DestroyBody.h"
+#include "pe/utility/GetBody.h"
+
+#include "pe_coupling/mapping/all.h"
+#include "pe_coupling/momentum_exchange_method/all.h"
+#include <pe_coupling/utility/all.h>
+
+#include "vtk/all.h"
+#include "field/vtk/all.h"
+
+namespace body_mapping_test
+{
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+using walberla::uint_t;
+
+// PDF field, flag field & body field
+typedef lbm::D3Q19< lbm::collision_model::SRT >  LatticeModel_T;
+
+typedef LatticeModel_T::Stencil                         Stencil_T;
+typedef lbm::PdfField< LatticeModel_T >                 PdfField_T;
+
+const uint_t FieldGhostLayers = 1;
+
+typedef walberla::uint8_t                 flag_t;
+typedef FlagField< flag_t >               FlagField_T;
+typedef GhostLayerField< pe::BodyID, 1 >  BodyField_T;
+
+// boundary handling
+typedef lbm::NoSlip< LatticeModel_T, flag_t > NoSlip_T;
+typedef pe_coupling::SimpleBB< LatticeModel_T, FlagField_T >  MO_T;
+typedef boost::tuples::tuple< NoSlip_T, MO_T > BoundaryConditions_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+
+typedef boost::tuple<pe::Sphere> BodyTypeTuple ;
+
+///////////
+// FLAGS //
+///////////
+
+const FlagUID Fluid_Flag ( "fluid" );
+const FlagUID MO_Flag ( "moving obstacle" );
+const FlagUID FormerMO_Flag ( "former moving obstacle" );
+const FlagUID NoSlip_Flag  ( "no slip" );
+
+
+/////////////////////////////////////
+// BOUNDARY HANDLING CUSTOMIZATION //
+/////////////////////////////////////
+
+class MyBoundaryHandling
+{
+public:
+
+   MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID, const BlockDataID & bodyFieldID ) :
+      flagFieldID_( flagFieldID ), pdfFieldID_( pdfFieldID ), bodyFieldID_ ( bodyFieldID ) {}
+
+   BoundaryHandling_T * operator()( IBlock* const block, const StructuredBlockStorage* const storage ) const;
+
+private:
+
+   const BlockDataID flagFieldID_;
+   const BlockDataID pdfFieldID_;
+   const BlockDataID bodyFieldID_;
+
+}; // class MyBoundaryHandling
+
+BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const StructuredBlockStorage * const storage ) const
+{
+   WALBERLA_ASSERT_NOT_NULLPTR( block );
+   WALBERLA_ASSERT_NOT_NULLPTR( storage );
+
+   FlagField_T * flagField       = block->getData< FlagField_T >( flagFieldID_ );
+   PdfField_T *  pdfField        = block->getData< PdfField_T > ( pdfFieldID_ );
+   BodyField_T * bodyField       = block->getData< BodyField_T >( bodyFieldID_ );
+
+   const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
+
+   BoundaryHandling_T * handling = new BoundaryHandling_T( "fixed obstacle boundary handling", flagField, fluid,
+                                                           boost::tuples::make_tuple( NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+                                                                                      MO_T (  "MO_BB",  MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+
+   handling->fillWithDomain( FieldGhostLayers );
+
+   return handling;
+}
+
+
+class MappingChecker
+{
+public:
+   MappingChecker(const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & bodyStorageID, const shared_ptr<pe::BodyStorage> & globalBodyStorage,
+                  const BlockDataID & boundaryHandlingID, const BlockDataID & bodyFieldID, real_t sphereRadius) :
+         blocks_( blocks ), bodyStorageID_( bodyStorageID ), globalBodyStorage_( globalBodyStorage ),
+         boundaryHandlingID_( boundaryHandlingID ), bodyFieldID_( bodyFieldID ),
+         sphereVolume_( math::M_PI * real_t(4) / real_t(3) * sphereRadius * sphereRadius * sphereRadius )
+   { }
+
+   // check the mapping in the inner domain of the block and check mapped volume against real sphere volume
+   void operator()(std::string & testIdentifier)
+   {
+      uint_t cellCounter( uint_t(0) );
+      for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
+      {
+         BoundaryHandling_T * boundaryHandling = blockIt->getData< BoundaryHandling_T >( boundaryHandlingID_ );
+         FlagField_T *        flagField        = boundaryHandling->getFlagField();
+
+         auto xyzSize = flagField->xyzSize();
+
+         // look for bodies in body storage
+         for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt ) {
+            for (auto cellIt = xyzSize.begin(); cellIt != xyzSize.end(); ++cellIt) {
+               Vector3<real_t> cellCenter = blocks_->getBlockLocalCellCenter(*blockIt, *cellIt);
+               Vector3<real_t> distance = cellCenter - bodyIt->getPosition();
+               if (bodyIt->containsPoint(cellCenter)) {
+                  WALBERLA_CHECK(boundaryHandling->isBoundary(*cellIt),
+                                 testIdentifier << "Invalid mapping in cell " << *cellIt
+                                                << " with center at " << cellCenter
+                                                << " - expected boundary cell. Distance cell center - body center = "
+                                                << distance.length() << ".");
+                  ++cellCounter;
+               } else {
+                  WALBERLA_CHECK(boundaryHandling->isDomain(*cellIt),
+                                 testIdentifier << "Invalid mapping in cell " << *cellIt
+                                                << " with center at " << cellCenter
+                                                << " - expected domain cell. Distance cell center - body center = "
+                                                << distance.length() << ".");
+               }
+            }
+         }
+         // look for bodies in global body storage
+         for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
+         {
+            for (auto cellIt = xyzSize.begin(); cellIt != xyzSize.end(); ++cellIt) {
+               Vector3<real_t> cellCenter = blocks_->getBlockLocalCellCenter(*blockIt, *cellIt);
+               Vector3<real_t> distance = cellCenter - bodyIt->getPosition();
+               if (bodyIt->containsPoint(cellCenter)) {
+                  WALBERLA_CHECK(boundaryHandling->isBoundary(*cellIt),
+                                 testIdentifier << "Invalid mapping in cell " << *cellIt
+                                                << " with center at " << cellCenter
+                                                << " - expected boundary cell. Distance cell center - body center = "
+                                                << distance.length() << ".");
+                  ++cellCounter;
+               } else {
+                  if( globalBodyStorage_->size() <= uint_t(1) ) {
+                     WALBERLA_CHECK(boundaryHandling->isDomain(*cellIt),
+                                    testIdentifier << "Invalid mapping in cell " << *cellIt
+                                                   << " with center at " << cellCenter
+                                                   << " - expected domain cell. Distance cell center - body center = "
+                                                   << distance.length() << ".");
+                  }
+               }
+            }
+         }
+      }
+      // mapped volume has to be - approximately - the same as the real volume
+      real_t mappedVolume = real_c(cellCounter); // dx=1
+      WALBERLA_CHECK(std::fabs( mappedVolume - sphereVolume_ ) / sphereVolume_ <= real_t(0.1),
+                     "Mapped volume " << mappedVolume << " does not fit to real sphere volume " << sphereVolume_ << ".");
+   }
+
+   // checks only the mapping in the ghost layers
+   void checkGhostLayer(std::string & testIdentifier)
+   {
+      for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
+      {
+         BoundaryHandling_T * boundaryHandling = blockIt->getData< BoundaryHandling_T >( boundaryHandlingID_ );
+         FlagField_T *        flagField        = boundaryHandling->getFlagField();
+
+         auto xyzSizeWGl = flagField->xyzSizeWithGhostLayer();
+
+         // look for bodies in body storage
+         for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            for( auto cellIt = xyzSizeWGl.begin(); cellIt != xyzSizeWGl.end(); ++cellIt )
+            {
+               if( flagField->isInInnerPart(*cellIt) ) continue;
+
+               Vector3<real_t> cellCenter = blocks_->getBlockLocalCellCenter( *blockIt, *cellIt );
+               Vector3<real_t> distance = cellCenter - bodyIt->getPosition();
+               if( bodyIt->containsPoint(cellCenter))
+               {
+                  WALBERLA_CHECK( boundaryHandling->isBoundary(*cellIt),
+                                  testIdentifier << ": Invalid mapping in ghost layer cell " << *cellIt
+                                                 << " with center at " << cellCenter
+                                                 << " - expected boundary cell. Distance cell center - body center = "
+                                                 << distance.length() << ".");
+               }
+               else
+               {
+                  WALBERLA_CHECK( boundaryHandling->isDomain(*cellIt),
+                                  testIdentifier << ": Invalid mapping in ghost layer cell " << *cellIt
+                                                 << " with center at " << cellCenter
+                                                 << " - expected domain cell. Distance cell center - body center = "
+                                                 << distance.length() << "." );
+               }
+            }
+         }
+
+         // look for bodies in global body storage
+         for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
+         {
+            for( auto cellIt = xyzSizeWGl.begin(); cellIt != xyzSizeWGl.end(); ++cellIt )
+            {
+
+               if( flagField->isInInnerPart(*cellIt) ) continue;
+
+               Vector3<real_t> cellCenter = blocks_->getBlockLocalCellCenter( *blockIt, *cellIt );
+               Vector3<real_t> distance = cellCenter - bodyIt->getPosition();
+               if( bodyIt->containsPoint(cellCenter))
+               {
+                  WALBERLA_CHECK( boundaryHandling->isBoundary(*cellIt),
+                                  testIdentifier << ": Invalid mapping in ghost layer cell " << *cellIt
+                                                 << " with center at " << cellCenter
+                                                 << " - expected boundary cell. Distance cell center - body center = "
+                                                 << distance.length() << ".");
+               }
+               else
+               {
+                  if( globalBodyStorage_->size() <= uint_t(1) ) {
+                     WALBERLA_CHECK(boundaryHandling->isDomain(*cellIt),
+                                    testIdentifier << ": Invalid mapping in ghost layer cell " << *cellIt
+                                                   << " with center at " << cellCenter
+                                                   << " - expected domain cell. Distance cell center - body center = "
+                                                   << distance.length() << ".");
+                  }
+               }
+            }
+         }
+      }
+   }
+
+   // check if the pointer in the body field is the correct one
+   void checkBodyField(std::string & testIdentifier)
+   {
+      for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
+      {
+
+         BodyField_T * bodyField = blockIt->getData< BodyField_T >( bodyFieldID_ );
+
+         auto xyzSizeWGl = bodyField->xyzSizeWithGhostLayer();
+
+         // look for bodies in body storage
+         for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            for( auto cellIt = xyzSizeWGl.begin(); cellIt != xyzSizeWGl.end(); ++cellIt )
+            {
+               Vector3<real_t> cellCenter = blocks_->getBlockLocalCellCenter( *blockIt, *cellIt );
+               Vector3<real_t> distance = cellCenter - bodyIt->getPosition();
+               if( bodyIt->containsPoint(cellCenter) )
+               {
+                  WALBERLA_CHECK_EQUAL( bodyIt->getSystemID(), (*bodyField)(*cellIt)->getSystemID(),
+                                        testIdentifier << "Invalid mapping in cell " << *cellIt
+                                                       << " with center at " << cellCenter
+                                                       << " - expected body. Distance cell center - body center = "
+                                                       << distance.length() << ".");
+               }
+               else
+               {
+                  WALBERLA_CHECK_NULLPTR( bodyField->get(*cellIt),
+                                          testIdentifier << "Invalid mapping in cell " << *cellIt
+                                                         << " with center at " << cellCenter
+                                                         << " - expected NULL pointer. Distance cell center - body center = "
+                                                         << distance.length() << "." );
+               }
+            }
+         }
+
+         // look for bodies in global body storage
+         for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
+         {
+            for( auto cellIt = xyzSizeWGl.begin(); cellIt != xyzSizeWGl.end(); ++cellIt )
+            {
+               Vector3<real_t> cellCenter = blocks_->getBlockLocalCellCenter( *blockIt, *cellIt );
+               Vector3<real_t> distance = cellCenter - bodyIt->getPosition();
+               if( bodyIt->containsPoint(cellCenter))
+               {
+                  WALBERLA_CHECK_EQUAL( bodyIt->getSystemID(), (*bodyField)(*cellIt)->getSystemID(),
+                                  testIdentifier << "Invalid mapping in cell " << *cellIt
+                                                 << " with center at " << cellCenter
+                                                 << " - expected body. Distance cell center - body center = "
+                                                 << distance.length() << ".");
+               }
+               else
+               {
+                  if( globalBodyStorage_->size() <= uint_t(1) )
+                  {
+                     WALBERLA_CHECK_NULLPTR( bodyField->get(*cellIt),
+                                     testIdentifier << "Invalid mapping in cell " << *cellIt
+                                                    << " with center at " << cellCenter
+                                                    << " - expected NULL pointer. Distance cell center - body center = "
+                                                    << distance.length() << ".");
+                  }
+               }
+            }
+         }
+      }
+   }
+
+
+private:
+   shared_ptr< StructuredBlockStorage > blocks_;
+   const BlockDataID bodyStorageID_;
+   shared_ptr<pe::BodyStorage> globalBodyStorage_;
+   const BlockDataID boundaryHandlingID_;
+   const BlockDataID bodyFieldID_;
+   real_t sphereVolume_;
+
+};
+
+class MappingResetter
+{
+public:
+   MappingResetter(const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & boundaryHandlingID, const BlockDataID & bodyFieldID) :
+         blocks_( blocks ), boundaryHandlingID_( boundaryHandlingID ), bodyFieldID_ ( bodyFieldID )
+   { }
+
+   void operator()()
+   {
+      for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
+      {
+         BoundaryHandling_T * boundaryHandling = blockIt->getData< BoundaryHandling_T >( boundaryHandlingID_ );
+         FlagField_T *        flagField        = boundaryHandling->getFlagField();
+         BodyField_T *        bodyField        = blockIt->getData< BodyField_T >( bodyFieldID_ );
+
+         auto xyzSizeWGl = flagField->xyzSizeWithGhostLayer();
+
+         // reset to domain (fluid)
+         boundaryHandling->forceDomain(xyzSizeWGl);
+
+         for( auto cellIt = xyzSizeWGl.begin(); cellIt != xyzSizeWGl.end(); ++cellIt )
+         {
+            // reset body field
+            bodyField->get(*cellIt) = NULL;
+         }
+      }
+   }
+
+private:
+   shared_ptr< StructuredBlockStorage > blocks_;
+   const BlockDataID boundaryHandlingID_;
+   const BlockDataID bodyFieldID_;
+};
+
+
+/*!\brief Test case for the available body mapping functions
+ *
+ * There are two types of body mapping functions:
+ *  - src/pe_coupling/mapping/BodyMapping.h:
+ *      - to map a standard LBM boundary condition to the body, e.g. NoSlip
+ *  - src/pe_coupling/momentum_exchange_method/BodyMapping.h:
+ *      - to map the coupling boundaries (rc/pe_coupling/momentum_exchange_method/boundary) to the body
+ *      - the surface velocity of the mapped body is seen by the fluid
+ *      - forces acting on the mapped bodies are calculated and set onto the body
+ *
+ * The following mapping functions are tested:
+ *  - mapAllBodies
+ *  - selectRegularBodies
+ *  - selectFixedBodies
+ *  - selectGlobalBodies
+ * The mapping inside a block, at a regular block boarder and at a periodic block boarder is tested.
+ * The No-Slip mapping as well as the moving boundary mapping is tested.
+ * Regular, global and infinite-mass spheres are tested.
+ * After each test, the sphere is destroyed and the mapping is erased from the datastructures.
+ *
+ * Note that global bodies do not see the periodicity and thus the periodic copy has to be created explicitly.
+ *
+ */
+//////////
+// MAIN //
+//////////
+int main( int argc, char **argv )
+{
+   debug::enterTestMode();
+
+   mpi::Environment env( argc, argv );
+
+   ///////////////////////////
+   // SIMULATION PROPERTIES //
+   ///////////////////////////
+
+   bool writeVTK = false;
+   const real_t omega  = real_t(1);
+   const real_t dx     = real_t(1);
+   const real_t radius = real_t(5);
+
+   ///////////////////////////
+   // DATA STRUCTURES SETUP //
+   ///////////////////////////
+
+   Vector3<uint_t> blocksPerDirection(uint_t(3), uint_t(1), uint_t(1));
+   Vector3<uint_t> cellsPerBlock(uint_t(20), uint_t(20), uint_t(20));
+   Vector3<bool> periodicity(true, false, false);
+
+   // create fully periodic domain with refined blocks
+   auto blocks = blockforest::createUniformBlockGrid( blocksPerDirection[0], blocksPerDirection[1], blocksPerDirection[2],
+                                                      cellsPerBlock[0], cellsPerBlock[1], cellsPerBlock[2],
+                                                      dx,
+                                                      0, false, false,
+                                                      periodicity[0], periodicity[1], periodicity[2],
+                                                      false );
+
+   // create the lattice model
+   LatticeModel_T latticeModel = LatticeModel_T( omega );
+
+   // add PDF field ( uInit = <0.1,0,0>, rhoInit = 1 )
+   BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel,
+                                                                         Vector3<real_t>(real_t(0)), real_t(1),
+                                                                         FieldGhostLayers, field::zyxf );
+
+   // add flag field
+   BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field", FieldGhostLayers );
+
+   // add body field
+   BlockDataID bodyFieldID = field::addToStorage<BodyField_T>( blocks, "body field", NULL, field::zyxf, FieldGhostLayers );
+
+   // add boundary handling
+   BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >(MyBoundaryHandling( flagFieldID, pdfFieldID, bodyFieldID ), "boundary handling" );
+
+   // pe body storage
+   pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
+   shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
+   auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "Storage");
+   auto sphereMaterialID = pe::createMaterial( "sphereMat", real_t(1) , real_t(0.3), real_t(0.2), real_t(0.2), real_t(0.24), real_t(200), real_t(200), real_t(0), real_t(0) );
+
+   // pe coupling
+   const real_t overlap = real_t( 1.5 ) * dx;
+
+   // vtk
+   auto flagFieldVTK = vtk::createVTKOutput_BlockData( blocks, "flag_field" );
+   flagFieldVTK->addCellDataWriter( make_shared< field::VTKWriter< FlagField_T > >( flagFieldID, "FlagField" ) );
+
+   // testing utilities
+   MappingChecker mappingChecker(blocks, bodyStorageID, globalBodyStorage, boundaryHandlingID, bodyFieldID, radius);
+   MappingResetter mappingResetter(blocks, boundaryHandlingID, bodyFieldID);
+   pe_coupling::BodyMapping<BoundaryHandling_T> regularBodyMapper(blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies );
+   pe_coupling::BodyMapping<BoundaryHandling_T> globalBodyMapper(blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectGlobalBodies );
+
+
+   // sphere positions for test scenarios
+   Vector3<real_t> positionInsideBlock(real_t(10), real_t(10), real_t(10));
+   Vector3<real_t> positionAtBlockBoarder(real_t(19), real_t(10), real_t(10));
+   Vector3<real_t> positionAtPeriodicBoarder(real_t(1), real_t(10), real_t(10));
+
+   /////////////////////
+   // NO SLIP MAPPING //
+   /////////////////////
+
+   //////////////////////////////////
+   // TEST SPHERE INSIDE ONE BLOCK //
+   //////////////////////////////////
+
+   ////////////////////
+   // REGULAR SPHERE //
+   ////////////////////
+   {
+      std::string testIdentifier("Test: regular sphere inside block with no slip mapping ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionInsideBlock, radius, sphereMaterialID, false, true, false);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      pe_coupling::mapBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectRegularBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   /////////////////////
+   // GLOBAL SPHERE 1 //
+   /////////////////////
+   {
+      std::string testIdentifier("Test: global sphere inside block with no slip mapping, case 1 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionInsideBlock, radius, sphereMaterialID, true, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      pe_coupling::mapBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectGlobalBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   //////////////////
+   // FIXED SPHERE //
+   //////////////////
+   {
+      std::string testIdentifier("Test: fixed sphere inside block with no slip mapping ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionInsideBlock, radius, sphereMaterialID, false, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, true);
+
+      pe_coupling::mapBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectFixedBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, true);
+   }
+
+   //////////////////////////////////
+   // TEST SPHERE AT BLOCK BOARDER //
+   //////////////////////////////////
+
+   ////////////////////
+   // REGULAR SPHERE //
+   ////////////////////
+   {
+      std::string testIdentifier("Test: regular sphere at block boarder with no slip mapping ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtBlockBoarder, radius, sphereMaterialID, false, true, false);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      pe_coupling::mapBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectRegularBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   /////////////////////
+   // GLOBAL SPHERE 1 //
+   /////////////////////
+   {
+      std::string testIdentifier("Test: global sphere at block boarder with no slip mapping, case 1 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtBlockBoarder, radius, sphereMaterialID, true, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      pe_coupling::mapBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectGlobalBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   //////////////////
+   // FIXED SPHERE //
+   //////////////////
+   {
+      std::string testIdentifier("Test: fixed sphere at block boarder with no slip mapping ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtBlockBoarder, radius, sphereMaterialID, false, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, true);
+
+      pe_coupling::mapBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectFixedBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, true);
+   }
+
+   /////////////////////////////////////
+   // TEST SPHERE AT PERIODIC BOARDER //
+   /////////////////////////////////////
+
+   ////////////////////
+   // REGULAR SPHERE //
+   ////////////////////
+   {
+      std::string testIdentifier("Test: regular sphere at periodic boarder with no slip mapping ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtPeriodicBoarder, radius, sphereMaterialID, false, true, false);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      pe_coupling::mapBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectRegularBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   /////////////////////
+   // GLOBAL SPHERE 1 //
+   /////////////////////
+   {
+      std::string testIdentifier("Test: global sphere at periodic boarder with no slip mapping, case 1 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtPeriodicBoarder, radius, sphereMaterialID, true, false, true);
+
+      //NOTE: global bodies are not communicated, thus they do not follow periodicity!!!
+      //workaround: create the periodic copy explicitly
+      Vector3<real_t> positionAtPeriodicBoarderCopy(real_t(1) + real_c(blocksPerDirection[0]) * real_c(cellsPerBlock[0]), real_t(10), real_t(10));
+      pe::SphereID sp2 = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                          positionAtPeriodicBoarderCopy, radius, sphereMaterialID, true, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      pe_coupling::mapBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectGlobalBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp2->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   //////////////////
+   // FIXED SPHERE //
+   //////////////////
+   {
+      std::string testIdentifier("Test: fixed sphere at periodic boarder with no slip mapping ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtPeriodicBoarder, radius, sphereMaterialID, false, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, true);
+
+      pe_coupling::mapBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectFixedBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, true);
+   }
+
+   //////////////////////////
+   // MOVING BODY MAPPING //
+   /////////////////////////
+
+   //////////////////////////////////
+   // TEST SPHERE INSIDE ONE BLOCK //
+   //////////////////////////////////
+
+   //////////////////////
+   // REGULAR SPHERE 1 //
+   //////////////////////
+   {
+      std::string testIdentifier("Test: regular sphere inside block with moving body mapping, case 1 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionInsideBlock, radius, sphereMaterialID, false, true, false);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectRegularBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   //////////////////////
+   // REGULAR SPHERE 2 //
+   //////////////////////
+   {
+      std::string testIdentifier("Test: regular sphere inside block with moving body mapping, case 2 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionInsideBlock, radius, sphereMaterialID, false, true, false);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) regularBodyMapper(&(*blockIt));
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   /////////////////////
+   // GLOBAL SPHERE 1 //
+   /////////////////////
+   {
+      std::string testIdentifier("Test: global sphere inside block with moving body mapping, case 1 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionInsideBlock, radius, sphereMaterialID, true, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectGlobalBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   //////////////////////
+   // GLOBAL SPHERE  2 //
+   //////////////////////
+   {
+      std::string testIdentifier("Test: global sphere inside block with moving body mapping, case 2 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionInsideBlock, radius, sphereMaterialID, true, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) globalBodyMapper(&(*blockIt));
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   //////////////////
+   // FIXED SPHERE //
+   //////////////////
+   {
+      std::string testIdentifier("Test: fixed sphere inside block with moving body mapping ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionInsideBlock, radius, sphereMaterialID, false, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, true);
+
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectFixedBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, true);
+   }
+
+   //////////////////////////////////
+   // TEST SPHERE AT BLOCK BOARDER //
+   //////////////////////////////////
+
+   //////////////////////
+   // REGULAR SPHERE 1 //
+   //////////////////////
+   {
+      std::string testIdentifier("Test: regular sphere at block boarder with moving body mapping, case 1 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtBlockBoarder, radius, sphereMaterialID, false, true, false);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectRegularBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   //////////////////////
+   // REGULAR SPHERE 2 //
+   //////////////////////
+   {
+      std::string testIdentifier("Test: regular sphere at block boarder with moving body mapping, case 2 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtBlockBoarder, radius, sphereMaterialID, false, true, false);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) regularBodyMapper(&(*blockIt));
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   /////////////////////
+   // GLOBAL SPHERE 1 //
+   /////////////////////
+   {
+      std::string testIdentifier("Test: global sphere at block boarder with moving body mapping, case 1 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtBlockBoarder, radius, sphereMaterialID, true, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectGlobalBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   //////////////////////
+   // GLOBAL SPHERE  2 //
+   //////////////////////
+   {
+      std::string testIdentifier("Test: global sphere at block boarder with moving body mapping, case 2 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtBlockBoarder, radius, sphereMaterialID, true, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) globalBodyMapper(&(*blockIt));
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   //////////////////
+   // FIXED SPHERE //
+   //////////////////
+   {
+      std::string testIdentifier("Test: fixed sphere at block boarder with moving body mapping ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtBlockBoarder, radius, sphereMaterialID, false, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, true);
+
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectFixedBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, true);
+   }
+
+   /////////////////////////////////////
+   // TEST SPHERE AT PERIODIC BOARDER //
+   /////////////////////////////////////
+
+   ///////////////////////
+   // REGULAR SPHERE  1 //
+   //////////////////////
+   {
+      std::string testIdentifier("Test: regular sphere at periodic boarder with moving body mapping, case 1 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtPeriodicBoarder, radius, sphereMaterialID, false, true, false);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectRegularBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   //////////////////////
+   // REGULAR SPHERE 2 //
+   //////////////////////
+   {
+      std::string testIdentifier("Test: regular sphere at periodic boarder with moving body mapping, case 2 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtPeriodicBoarder, radius, sphereMaterialID, false, true, false);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) regularBodyMapper(&(*blockIt));
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   /////////////////////
+   // GLOBAL SPHERE 1 //
+   /////////////////////
+   {
+      std::string testIdentifier("Test: global sphere at periodic boarder with moving body mapping, case 1 ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtPeriodicBoarder, radius, sphereMaterialID, true, false, true);
+
+      //NOTE: global bodies are not communicated, thus they do not follow periodicity!!!
+      //workaround: create the periodic copy explicitly
+      Vector3<real_t> positionAtPeriodicBoarderCopy(real_t(1) + real_c(blocksPerDirection[0]) * real_c(cellsPerBlock[0]), real_t(10), real_t(10));
+      pe::SphereID sp2 = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                          positionAtPeriodicBoarderCopy, radius, sphereMaterialID, true, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectGlobalBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp2->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, false);
+   }
+
+   //////////////////
+   // FIXED SPHERE //
+   //////////////////
+   {
+      std::string testIdentifier("Test: fixed sphere at periodic boarder with moving body mapping ");
+      WALBERLA_LOG_DEVEL(testIdentifier << " - started");
+      pe::SphereID sp = pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                                         positionAtPeriodicBoarder, radius, sphereMaterialID, false, false, true);
+
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, true);
+
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectFixedBodies );
+
+      if( writeVTK ) flagFieldVTK->write();
+
+      mappingChecker(testIdentifier);
+      mappingChecker.checkGhostLayer(testIdentifier);
+      mappingChecker.checkBodyField(testIdentifier);
+      mappingResetter();
+
+      WALBERLA_LOG_DEVEL(testIdentifier << " - ended");
+
+      pe::destroyBodyBySID(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, sp->getSystemID());
+      pe::syncNextNeighbors<BodyTypeTuple>(blocks->getBlockForest(), bodyStorageID, static_cast<WcTimingTree *>(NULL), overlap, true);
+   }
+
+
+   return 0;
+
+}
+
+} //namespace body_mapping_test
+
+int main( int argc, char **argv ){
+   body_mapping_test::main(argc, argv);
+}
diff --git a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
index 581f51deb..1e4bf71b5 100644
--- a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
@@ -457,7 +457,7 @@ int main( int argc, char **argv )
    // PE //
    ////////
 
-   pe::BodyStorage globalBodyStorage;
+   shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
    pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
    auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "pe Body Storage");
 
@@ -476,7 +476,7 @@ int main( int argc, char **argv )
 
    // create the sphere in the middle of the domain
    Vector3<real_t> position (real_c(setup.length) * real_c(0.5));
-   pe::createSphere(globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, position, setup.radius );
+   pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, position, setup.radius );
 
    // synchronize often enough for large bodies
    for( uint_t i = 0; i < XBlocks / 2; ++i)
@@ -527,13 +527,13 @@ int main( int argc, char **argv )
    if( method == MEMVariant::CLI )
    {
       // uses a higher order boundary condition (CLI)
-      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag );
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_CLI_Flag );
    }else if ( method == MEMVariant::MR ){
       // uses a higher order boundary condition (MR)
-      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MO_MR_Flag );
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID,  MO_MR_Flag );
    }else{
       // uses standard bounce back boundary conditions
-      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MO_BB_Flag );
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID,  MO_BB_Flag );
    }
 
    // since external forcing is applied, the evaluation of the velocity has to be carried out directly after the streaming step
diff --git a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
index 49ac5f453..a084604c7 100644
--- a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
@@ -504,7 +504,7 @@ int main( int argc, char **argv )
    // PE //
    ////////
 
-   pe::BodyStorage globalBodyStorage;
+   shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
    pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
    auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "pe Body Storage");
 
@@ -517,7 +517,7 @@ int main( int argc, char **argv )
 
    // create the sphere in the middle of the domain
    Vector3<real_t> position ( setup.length * real_c(0.5));
-   pe::createSphere(globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, position, setup.radius );
+   pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, position, setup.radius );
 
    // synchronize often enough for large bodies
    for( uint_t i = 0; i < std::max( uint_c(1), XBlocks / uint_c(2) ) * ( uint_t(1) << ( setup.levels - uint_t(1) ) ); ++i)
@@ -549,10 +549,10 @@ int main( int argc, char **argv )
    if( MO_CLI )
    {
       // uses a higher order boundary condition (CLI)
-      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag );
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_CLI_Flag );
    }else{
       // uses standard bounce back boundary conditions
-      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MO_BB_Flag );
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID,  MO_BB_Flag );
    }
 
 
diff --git a/tests/pe_coupling/momentum_exchange_method/GlobalBodyAsBoundaryMEMStaticRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/GlobalBodyAsBoundaryMEMStaticRefinement.cpp
index f3d40bc3a..e833f9a71 100644
--- a/tests/pe_coupling/momentum_exchange_method/GlobalBodyAsBoundaryMEMStaticRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/GlobalBodyAsBoundaryMEMStaticRefinement.cpp
@@ -341,6 +341,7 @@ int main( int argc, char **argv )
    // set up pe functionality
    shared_ptr<pe::BodyStorage>  globalBodyStorage = make_shared<pe::BodyStorage>();
    pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
+   auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "pe Body Storage");
 
    // create pe bodies
 
@@ -382,7 +383,7 @@ int main( int argc, char **argv )
    BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >( MyBoundaryHandling( flagFieldID, pdfFieldID, bodyFieldID ), "boundary handling" );
 
    // map planes into the LBM simulation -> act as no-slip boundaries
-   pe_coupling::mapMovingGlobalBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, *globalBodyStorage, bodyFieldID, MO_SBB_Flag );
+   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_SBB_Flag );
 
 
    ///////////////
diff --git a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
index 62a3cfd2c..2569114b6 100644
--- a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
@@ -681,8 +681,7 @@ int main( int argc, char **argv )
 
    if( ! ( sphSphTest || sphWallTest || funcTest ) )
    {
-      std::cerr << "Specify either --sphSphTest, --sphWallTest, or --funcTest !" << std::endl;
-      return EXIT_FAILURE;
+      WALBERLA_ABORT("Specify either --sphSphTest, --sphWallTest, or --funcTest !");
    }
 
    ///////////////////////////
@@ -778,7 +777,7 @@ int main( int argc, char **argv )
    // PE COUPLING //
    /////////////////
 
-   shared_ptr<pe::BodyStorage>  globalBodyStorage = make_shared<pe::BodyStorage>();
+   shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
    pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
 
    auto bodyStorageID  = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "pe Body Storage");
@@ -922,7 +921,7 @@ int main( int argc, char **argv )
    }
 
    // map pe bodies into the LBM simulation
-   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag );
+   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag );
 
    ///////////////
    // TIME LOOP //
@@ -932,14 +931,14 @@ int main( int argc, char **argv )
 
    // sweep for updating the pe body mapping into the LBM simulation
    timeloop.add()
-      << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag, FormerMO_Flag ), "Body Mapping" );
+      << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
 
    // sweep for restoring PDFs in cells previously occupied by pe bodies
    typedef pe_coupling::EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T > Reconstructor_T;
    Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID );
    timeloop.add()
-      << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,
-                                                                                                        reconstructor, FormerMO_Flag, Fluid_Flag  ), "PDF Restore" );
+      << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID,
+                                                                                                       reconstructor, FormerMO_Flag, Fluid_Flag ), "PDF Restore" );
    // setup of the LBM communication for synchronizing the pdf field between neighboring blocks
    boost::function< void () > commFunction;
    if( fullPDFSync )
diff --git a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
index 2e351fa79..15d44380b 100644
--- a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
@@ -193,9 +193,6 @@ public:
       {
          for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
          {
-            if( bodyIt->isFixed() || !bodyIt->isFinite() )
-                continue;
-
             WALBERLA_CHECK( aabb_.contains( bodyIt->getPosition()[0], bodyIt->getPosition()[1], bodyIt->getPosition()[2] ) );
          }
       }
@@ -225,9 +222,6 @@ public:
       {
          for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
          {
-            if( bodyIt->isFixed() || !bodyIt->isFinite() )
-               continue;
-
             const auto & u = bodyIt->getLinearVel();
             WALBERLA_CHECK_LESS_EQUAL( (u[0]*u[0] + u[1]*u[1] + u[2]*u[2]), uMax_ );
 
@@ -432,13 +426,13 @@ int main( int argc, char **argv )
    pe::createPlane( *globalBodyStorage, 0, Vector3<real_t>(0,-1,0), Vector3<real_t>(0,real_c(width),0), material );
 
    // spheres as obstacles
-   std::vector<walberla::id_t> globalBodiesToBeMapped;
+   std::vector<pe::BodyID> globalBodiesToBeMapped;
    auto globalSphere1 = pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, Vector3<real_t>(real_c(length) / real_t(2), real_t(50), real_t(110)), real_t(60), material, true, false, true );
-   globalBodiesToBeMapped.push_back(globalSphere1->getSystemID() );
+   globalBodiesToBeMapped.push_back(globalSphere1);
    auto globalSphere2 = pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, Vector3<real_t>(                 real_t(0), real_t(50), -real_t(60)), real_t(80), material, true, false, true );
-   globalBodiesToBeMapped.push_back(globalSphere2->getSystemID() );
+   globalBodiesToBeMapped.push_back(globalSphere2);
    auto globalSphere3 = pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, Vector3<real_t>(            real_c(length), real_t(50), -real_t(60)), real_t(80), material, true, false, true );
-   globalBodiesToBeMapped.push_back(globalSphere3->getSystemID() );
+   globalBodiesToBeMapped.push_back(globalSphere3);
 
    // local bodies: moving spheres
    const real_t radius = real_t(10);
@@ -500,13 +494,16 @@ int main( int argc, char **argv )
    // map pe bodies into the LBM simulation
    // global bodies act as no-slip obstacles and are not treated by the fluid-solid coupling
    // special care has to be taken here that only the global spheres (not the planes) are mapped since the planes would overwrite the already set boundary flags
-   for( auto globalBodyIt = globalBodiesToBeMapped.begin(); globalBodyIt != globalBodiesToBeMapped.end(); ++globalBodyIt )
+   for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
    {
-      pe_coupling::mapGlobalBody< BoundaryHandling_T >( *globalBodyIt, *blocks, boundaryHandlingID, *globalBodyStorage, NoSlip_Flag );
+      for( auto globalBodyIt = globalBodiesToBeMapped.begin(); globalBodyIt != globalBodiesToBeMapped.end(); ++globalBodyIt )
+      {
+         pe_coupling::mapBody< BoundaryHandling_T >( *globalBodyIt, *blockIt, *blocks, boundaryHandlingID, NoSlip_Flag );
+      }
    }
 
    // moving bodies are handled by the momentum exchange method
-   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag );
+   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectRegularBodies );
 
    ///////////////
    // TIME LOOP //
@@ -516,15 +513,15 @@ int main( int argc, char **argv )
 
    // sweep for updating the pe body mapping into the LBM simulation
    timeloop.add()
-      << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag, FormerMO_Flag ),
+      << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ),
                 "Body Mapping" );
 
    // sweep for restoring PDFs in cells previously occupied by pe bodies
    typedef pe_coupling::EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T > Reconstructor_T;
    Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID );
    timeloop.add()
-      << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,
-                                                                                                        reconstructor, FormerMO_Flag, Fluid_Flag  ), "PDF Restore" );
+      << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID,
+                                                                                                       reconstructor, FormerMO_Flag, Fluid_Flag  ), "PDF Restore" );
 
    // setup of the LBM communication for synchronizing the pdf field between neighboring blocks
    boost::function< void () > commFunction;
diff --git a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
index 055c68b0d..6cfeb3dea 100644
--- a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
@@ -534,7 +534,7 @@ int main( int argc, char **argv )
    /////////////////
 
    // set up pe functionality
-   shared_ptr<pe::BodyStorage>  globalBodyStorage = make_shared<pe::BodyStorage>();
+   shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
    pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
 
    auto bodyStorageID  = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "pe Body Storage");
@@ -609,13 +609,13 @@ int main( int argc, char **argv )
    if( MO_CLI )
    {
       // uses a higher order boundary condition (CLI)
-      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag );
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_CLI_Flag, pe_coupling::selectRegularBodies );
    }else if ( MO_MR ){
       // uses a higher order boundary condition (MR)
-      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MO_MR_Flag );
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID,  MO_MR_Flag, pe_coupling::selectRegularBodies );
    }else{
       // uses standard bounce back boundary conditions
-      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MO_BB_Flag );
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID,  MO_BB_Flag, pe_coupling::selectRegularBodies );
    }
 
    ///////////////
@@ -635,13 +635,13 @@ int main( int argc, char **argv )
    if( MO_CLI )
    {
       timeloop.add()
-         << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag, FormerMO_Flag ), "Body Mapping" );
+         << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_CLI_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
    }else if ( MO_MR ){
       timeloop.add()
-         << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MO_MR_Flag, FormerMO_Flag ), "Body Mapping" );
+         << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID,  MO_MR_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
    }else{
       timeloop.add()
-         << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MO_BB_Flag, FormerMO_Flag ), "Body Mapping" );
+         << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID,  MO_BB_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
    }
 
    // sweep for restoring PDFs in cells previously occupied by pe bodies
@@ -652,7 +652,7 @@ int main( int argc, char **argv )
       Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationFinder, true );
       timeloop.add()
          << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >
-                      ( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag  ), "PDF Restore" );
+                      ( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag  ), "PDF Restore" );
    } else if ( eanReconstructor )
    {
       pe_coupling::SphereNormalExtrapolationDirectionFinder extrapolationFinder( blocks, bodyFieldID );
@@ -660,13 +660,13 @@ int main( int argc, char **argv )
       Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationFinder );
       timeloop.add()
          << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >
-                      ( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag  ), "PDF Restore" );
+                      ( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag  ), "PDF Restore" );
    } else {
       typedef pe_coupling::EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T > Reconstructor_T;
       Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID );
       timeloop.add()
          << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >
-                      ( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag  ), "PDF Restore" );
+                      ( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag  ), "PDF Restore" );
    }
 
    shared_ptr<pe_coupling::BodiesForceTorqueContainer> bodiesFTContainer1 = make_shared<pe_coupling::BodiesForceTorqueContainer>(blocks, bodyStorageID);
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
index 8e57d6fa6..0dfd102a5 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
@@ -499,10 +499,10 @@ int main( int argc, char **argv )
    BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >( MyBoundaryHandling( flagFieldID, pdfFieldID, bodyFieldID ), "boundary handling" );
 
    // map planes into the LBM simulation -> act as no-slip boundaries
-   pe_coupling::mapGlobalBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, *globalBodyStorage, NoSlip_Flag, false, true );
+   pe_coupling::mapBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectGlobalBodies );
 
    // map pe bodies into the LBM simulation
-   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag );
+   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectRegularBodies );
 
    ///////////////
    // TIME LOOP //
@@ -519,7 +519,7 @@ int main( int argc, char **argv )
 
    // sweep for updating the pe body mapping into the LBM simulation
    timeloop.add()
-         << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag, FormerMO_Flag ), "Body Mapping" );
+         << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), "Body Mapping" );
 
    // sweep for restoring PDFs in cells previously occupied by pe bodies
    pe_coupling::SphereNormalExtrapolationDirectionFinder extrapolationFinder( blocks, bodyFieldID );
@@ -527,7 +527,7 @@ int main( int argc, char **argv )
    Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationFinder );
    timeloop.add()
       << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >
-                  ( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag  ), "PDF Restore" );
+                  ( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), "PDF Restore" );
 
    shared_ptr<pe_coupling::BodiesForceTorqueContainer> bodiesFTContainer1 = make_shared<pe_coupling::BodiesForceTorqueContainer>(blocks, bodyStorageID);
    boost::function<void(void)> storeForceTorqueInCont1 = boost::bind(&pe_coupling::BodiesForceTorqueContainer::store, bodiesFTContainer1);
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
index 38184ac3c..75b9ff72e 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
@@ -720,10 +720,11 @@ int main( int argc, char **argv )
                                                           "boundary handling" );
 
    // map planes into the LBM simulation -> act as no-slip boundaries
-   pe_coupling::mapGlobalBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, *globalBodyStorage, NoSlip_Flag, false, true );
+   pe_coupling::mapBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectGlobalBodies );
 
    // map pe bodies into the LBM simulation
-   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag );
+   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectRegularBodies );
+
 
    // force averaging functionality
    shared_ptr<pe_coupling::BodiesForceTorqueContainer> bodiesFTContainer1 = make_shared<pe_coupling::BodiesForceTorqueContainer>(blocks, bodyStorageID);
@@ -783,14 +784,15 @@ int main( int argc, char **argv )
                                                   "pe Time Step", finestLevel );
 
    // add sweep for updating the pe body mapping into the LBM simulation
-   refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MO_Flag, FormerMO_Flag ), blocks ),
+   refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID,  MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), blocks ),
                                                  "Body Mapping", finestLevel );
 
    // add sweep for restoring PDFs in cells previously occupied by pe bodies
    pe_coupling::SphereNormalExtrapolationDirectionFinder extrapolationFinder( blocks, bodyFieldID );
    typedef pe_coupling::EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T, pe_coupling::SphereNormalExtrapolationDirectionFinder > Reconstructor_T;
    Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationFinder );
-   refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag  ), blocks ),
+   refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks,
+                                                 boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag ), blocks ),
                                                  "PDF Restore", finestLevel );
 
 
@@ -901,10 +903,11 @@ int main( int argc, char **argv )
          recreateBoundaryHandling(forest, boundaryHandlingID);
 
          // re-set the no-slip flags along the walls
-         pe_coupling::mapGlobalBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, *globalBodyStorage, NoSlip_Flag, false, true );
+         pe_coupling::mapBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectGlobalBodies );
 
          // re-map the body into the domain (initializing the bodyField as well)
-         pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag );
+         pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectRegularBodies );
+
       }
    }
 
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
index 988aff9e9..d93c8731e 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
@@ -588,10 +588,10 @@ int main( int argc, char **argv )
    BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >( MyBoundaryHandling( flagFieldID, pdfFieldID, bodyFieldID ), "boundary handling" );
 
    // map planes into the LBM simulation -> act as no-slip boundaries
-   pe_coupling::mapGlobalBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, *globalBodyStorage, NoSlip_Flag, false, true );
+   pe_coupling::mapBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, NoSlip_Flag, pe_coupling::selectGlobalBodies );
 
    // map pe bodies into the LBM simulation
-   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag );
+   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag, pe_coupling::selectRegularBodies );
 
    // force averaging functionality
    shared_ptr<pe_coupling::BodiesForceTorqueContainer> bodiesFTContainer1 = make_shared<pe_coupling::BodiesForceTorqueContainer>(blocks, bodyStorageID);
@@ -648,14 +648,16 @@ int main( int argc, char **argv )
                                                  "pe Time Step", finestLevel );
 
    // add sweep for updating the pe body mapping into the LBM simulation
-   refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MO_Flag, FormerMO_Flag ), blocks ),
+   refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID,
+                                                 globalBodyStorage, bodyFieldID,  MO_Flag, FormerMO_Flag, pe_coupling::selectRegularBodies ), blocks ),
                                                  "Body Mapping", finestLevel );
 
    // add sweep for restoring PDFs in cells previously occupied by pe bodies
    pe_coupling::SphereNormalExtrapolationDirectionFinder extrapolationFinder( blocks, bodyFieldID );
    typedef pe_coupling::EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T, pe_coupling::SphereNormalExtrapolationDirectionFinder > Reconstructor_T;
    Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationFinder );
-   refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag  ), blocks ),
+   refinementTimestep->addPostStreamVoidFunction(lbm::refinement::SweepAsFunctorWrapper( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T > ( blocks,
+                                                 boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID, reconstructor, FormerMO_Flag, Fluid_Flag  ), blocks ),
                                                  "PDF Restore", finestLevel );
 
 
diff --git a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
index 5680cb5ab..00ce16ed0 100644
--- a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
@@ -339,8 +339,8 @@ int main(int argc, char **argv) {
 
    // sweep for updating the pe body mapping into the LBM simulation
    timeloop.add()
-         << Sweep(pe_coupling::BodyMapping<BoundaryHandling_T>(blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,
-                                                               MO_BB_Flag, FormerMO_BB_Flag),
+         << Sweep(pe_coupling::BodyMapping<BoundaryHandling_T>(blocks, boundaryHandlingID, bodyStorageID, globalBodyStorage, bodyFieldID,
+                                                               MO_BB_Flag, FormerMO_BB_Flag, pe_coupling::selectRegularBodies),
                   "Body Mapping");
 
    // sweep for restoring PDFs in cells previously occupied by pe bodies
@@ -350,6 +350,7 @@ int main(int argc, char **argv) {
          << Sweep(pe_coupling::PDFReconstruction<LatticeModel_T, BoundaryHandling_T, Reconstructor_T>(blocks,
                                                                                                       boundaryHandlingID,
                                                                                                       bodyStorageID,
+                                                                                                      globalBodyStorage,
                                                                                                       bodyFieldID,
                                                                                                       reconstructor,
                                                                                                       FormerMO_BB_Flag,
@@ -364,8 +365,8 @@ int main(int argc, char **argv) {
    commFunction = scheme;
 
    // uses standard bounce back boundary conditions
-   pe_coupling::mapMovingBodies<BoundaryHandling_T>(*blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,
-                                                    MO_BB_Flag);
+   pe_coupling::mapMovingBodies<BoundaryHandling_T>(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID,
+                                                    MO_BB_Flag, pe_coupling::selectRegularBodies);
 
    auto sweep = lbm::makeCellwiseSweep<LatticeModel_T, FlagField_T>(pdfFieldID, flagFieldID, Fluid_Flag);
 
diff --git a/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp b/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
index f4c2df525..1b8d4677e 100644
--- a/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
@@ -297,7 +297,7 @@ int main( int argc, char **argv )
    /////////////////
 
    // set up pe functionality
-   shared_ptr<pe::BodyStorage>  globalBodyStorage = make_shared<pe::BodyStorage>();
+   shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
    pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
    auto bodyStorageID  = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "pe Body Storage");
 
@@ -339,7 +339,7 @@ int main( int argc, char **argv )
 
    // map pe bodies into the LBM simulation
    // moving bodies are handled by the momentum exchange method, thus act here as velocity boundary condition
-   pe_coupling::mapMovingGlobalBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, *globalBodyStorage, bodyFieldID, MO_Flag );
+   pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_Flag );
 
    ///////////////
    // TIME LOOP //
diff --git a/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
index 96bb5a31e..85cc5083e 100644
--- a/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
@@ -401,7 +401,7 @@ int main( int argc, char **argv )
    // PE //
    ////////
 
-   pe::BodyStorage globalBodyStorage;
+   shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
    pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
    auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "pe Body Storage");
 
@@ -420,7 +420,7 @@ int main( int argc, char **argv )
 
    // create the sphere in the middle of the domain
    Vector3<real_t> position (real_c(setup.length) * real_c(0.5));
-   auto sphere = pe::createSphere( globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, position, setup.radius );
+   auto sphere = pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, position, setup.radius );
    if ( sphere != NULL )
    {
       sphere->setAngularVel( real_c(0), setup.angularVel, real_c(0) );
@@ -476,13 +476,13 @@ int main( int argc, char **argv )
    if( MO_CLI )
    {
       // uses a higher order boundary condition (CLI)
-      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_CLI_Flag );
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, MO_CLI_Flag );
    }else if ( MO_MR ){
       // uses a higher order boundary condition (MR)
-      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MO_MR_Flag );
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID,  MO_MR_Flag );
    }else{
       // uses standard bounce back boundary conditions
-      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,  MO_BB_Flag );
+      pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID,  MO_BB_Flag );
    }
 
    shared_ptr< TorqueEval > torqueEval = make_shared< TorqueEval >( &timeloop, &setup, blocks, bodyStorageID, fileIO );
-- 
GitLab