diff --git a/src/pe_coupling/mapping/BodyBBMapping.cpp b/src/pe_coupling/mapping/BodyBBMapping.cpp
index 9defb26ecd2fc688ff0432903e38dda314691a16..7cec1f1fd591d25f6c02762d85039711267a8a9c 100644
--- a/src/pe_coupling/mapping/BodyBBMapping.cpp
+++ b/src/pe_coupling/mapping/BodyBBMapping.cpp
@@ -31,9 +31,12 @@ namespace pe_coupling {
 CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, const shared_ptr<StructuredBlockStorage> & blockStorage,
                         const uint_t numberOfGhostLayersToInclude )
 {
+
+   WALBERLA_ASSERT_NOT_NULLPTR( body );
+
    CellInterval cellBB;
 
-   if (body->isFinite())
+   if( body->isFinite() )
    {
       blockStorage->getCellBBFromAABB( cellBB, body->getAABB(), blockStorage->getLevel(block) );
    } else
@@ -41,8 +44,8 @@ CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, const
       blockStorage->getCellBBFromAABB( cellBB, body->getAABB().getIntersection( blockStorage->getDomain() ), blockStorage->getLevel(block) );
    }
 
-   cellBB.xMin() -= cell_idx_c(1); cellBB.yMin() -= cell_idx_c(1); cellBB.zMin() -= cell_idx_c(1);
-   cellBB.xMax() += cell_idx_c(1); cellBB.yMax() += cell_idx_c(1); cellBB.zMax() += cell_idx_c(1);
+   cellBB.xMin() -= cell_idx_t(1); cellBB.yMin() -= cell_idx_t(1); cellBB.zMin() -= cell_idx_t(1);
+   cellBB.xMax() += cell_idx_t(1); cellBB.yMax() += cell_idx_t(1); cellBB.zMax() += cell_idx_t(1);
 
    CellInterval blockBB = blockStorage->getBlockCellBB( block );
 
diff --git a/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.h b/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.h
index 73cf2303eb5f6a6096e9038ff0a85f9e612b830b..e23d38be85023d5d2976b62d0a995cb09a517de5 100644
--- a/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.h
+++ b/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.h
@@ -45,7 +45,9 @@ namespace pe_coupling {
 *   Different variants are available:
 *
 *   - FlagFieldNormalExtrapolationDirectionFinder:
-*     Uses information from the direct (N,S,E,W,T,B) neighbors' flags to find the extrapolation direction.
+*     Uses information from the direct (N,S,E,W,T,B) neighbors' flags to find the extrapolation direction
+*     by checking for the domain (here: fluid) flag.
+*     This can lead to incorrect directions in close packing scenarios or in the vicinity of other boundaries.
 *
 *   - SphereNormalExtrapolationDirectionFinder:
 *     Calculates the local normal of the body which is in its current form only correct for spherical bodies.
@@ -59,7 +61,7 @@ template < typename Stencil_T >
 void findCorrespondingLatticeDirection( const Vector3<real_t> & direction, Vector3<cell_idx_t> & correspondingLatticeDirection )
 {
    stencil::Direction correspondingDirection = stencil::C;
-   real_t innerProduct = real_c(0);
+   real_t innerProduct = real_t(0);
    for( auto d = Stencil_T::beginNoCenter(); d != Stencil_T::end(); ++d )
    {
       // compute inner product <dir,c_i>
@@ -138,11 +140,15 @@ void SphereNormalExtrapolationDirectionFinder
    BodyField_T * bodyField = block->getData< BodyField_T >( bodyFieldID_ );
 
    WALBERLA_ASSERT_NOT_NULLPTR( bodyField );
+   WALBERLA_ASSERT_NOT_NULLPTR( (*bodyField)(x,y,z) );
 
    real_t cx, cy, cz;
    blockStorage_->getBlockLocalCellCenter( *block, Cell(x,y,z), cx, cy, cz );
-   Vector3<real_t> bodyVelocity = (*bodyField)(x,y,z)->getPosition();
-   Vector3<real_t> bodyNormal( cx - bodyVelocity[0], cy - bodyVelocity[1], cz - bodyVelocity[2] );
+
+   Vector3<real_t> bodyCenterPosition = (*bodyField)(x,y,z)->getPosition();
+   WALBERLA_ASSERT( !math::isnan(bodyCenterPosition) );
+
+   Vector3<real_t> bodyNormal( cx - bodyCenterPosition[0], cy - bodyCenterPosition[1], cz - bodyCenterPosition[2] );
 
    findCorrespondingLatticeDirection< stencil::D3Q27 >( bodyNormal, extrapolationDirection );
 }
diff --git a/src/pe_coupling/momentum_exchange_method/restoration/PDFRestoration.h b/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h
similarity index 99%
rename from src/pe_coupling/momentum_exchange_method/restoration/PDFRestoration.h
rename to src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h
index eb5e7b7613ee8a6f58c133181d23cf98c4769ec6..9d751bf331336010acfab7c74f1f8a602449f841 100644
--- a/src/pe_coupling/momentum_exchange_method/restoration/PDFRestoration.h
+++ b/src/pe_coupling/momentum_exchange_method/restoration/PDFReconstruction.h
@@ -118,7 +118,7 @@ void PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructer_T >
    // reconstruct all missing PDFs (only inside the domain)
    if( optimizeForSmallObstacleFraction_ )
    {
-      const uint_t numberOfGhostLayersToInclude = uint_c(0);
+      const uint_t numberOfGhostLayersToInclude = uint_t(0);
 
       for( auto bodyIt = pe::BodyIterator::begin(*block, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
       {
diff --git a/src/pe_coupling/momentum_exchange_method/restoration/Reconstructor.h b/src/pe_coupling/momentum_exchange_method/restoration/Reconstructor.h
index f6234295821421e268ecf8ab4c44bd889227333e..8acb17a1e97fb993156de1050e72ac25dc89ba85 100644
--- a/src/pe_coupling/momentum_exchange_method/restoration/Reconstructor.h
+++ b/src/pe_coupling/momentum_exchange_method/restoration/Reconstructor.h
@@ -56,6 +56,7 @@ namespace pe_coupling {
 *
 *   - ExtrapolationReconstructor:
 *     Extrapolates the PDFs of three or two neighboring cells that lie in extrapolation direction.
+*     Optionally, a special treatment after the extrapolation can be used that sets certain moments directly (only D3Q19).
 *     Defaults to EquilibriumAndNonEquilibriumReconstructor if not enough cells in this direction are available.
 *
 *   EquilibriumAndNonEquilibriumReconstructor and ExtrapolationReconstructor need an extrapolation direction
@@ -139,7 +140,7 @@ real_t EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T >
    CellInterval localDomain = pdfField->xyzSize();
 
    uint_t nAverage = uint_t(0);
-   real_t averageDensity = real_c(0);
+   real_t averageDensity = real_t(0);
    for( auto neighborDir = stencil::D3Q27::beginNoCenter(); neighborDir != stencil::D3Q27::end(); ++neighborDir )
    {
       Cell neighbor( x + neighborDir.cx(), y + neighborDir.cy(), z + neighborDir.cz() );
@@ -179,6 +180,7 @@ public:
    { }
 
    void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block );
+   void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, const Vector3<cell_idx_t> & extrapolationDirection );
 
 private:
 
@@ -201,11 +203,24 @@ void EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling
 {
    WALBERLA_ASSERT_NOT_NULLPTR( block );
 
-   equilibriumReconstructor_( x, y, z, block );
-
    Vector3<cell_idx_t> extrapolationDirection(0);
    extrapolationDirectionFinder_.getDirection( x, y, z, block, extrapolationDirection );
-   setNonEquilibrium( x, y, z, block, extrapolationDirection);
+
+   (*this)(x, y, z, block, extrapolationDirection);
+}
+
+template< typename LatticeModel_T, typename BoundaryHandling_T, typename ExtrapolationDirectionFinder_T >
+void EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationDirectionFinder_T >
+::operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, const Vector3<cell_idx_t> & extrapolationDirection  )
+{
+   WALBERLA_ASSERT_NOT_NULLPTR( block );
+
+   equilibriumReconstructor_( x, y, z, block );
+
+   if( extrapolationDirection != cell_idx_t(0) )
+   {
+      setNonEquilibrium( x, y, z, block, extrapolationDirection);
+   }
 }
 
 template< typename LatticeModel_T, typename BoundaryHandling_T, typename ExtrapolationDirectionFinder_T >
@@ -219,7 +234,6 @@ void EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling
 
    WALBERLA_ASSERT_NOT_NULLPTR( pdfField );
    WALBERLA_ASSERT_NOT_NULLPTR( boundaryHandling );
-   WALBERLA_ASSERT( extrapolationDirection != cell_idx_t(0) );
 
    CellInterval localDomain = pdfField->xyzSize();
    Cell normalNeighbor( x + extrapolationDirection[0], y + extrapolationDirection[1], z + extrapolationDirection[2] );
@@ -246,20 +260,25 @@ class ExtrapolationReconstructor
 {
 public:
 
-   static_assert( (boost::is_same< typename LatticeModel_T::Stencil, stencil::D3Q19 >::value), "Currently only works with D3Q19 stencil for enforcing the no-slip condition!" );
-
    typedef lbm::PdfField< LatticeModel_T > PdfField_T;
    typedef Field< pe::BodyID, 1 >          BodyField_T;
 
    ExtrapolationReconstructor( const shared_ptr<StructuredBlockStorage> blockStorage, const BlockDataID & boundaryHandlingID,
                                const BlockDataID & pdfFieldID, const BlockDataID & bodyFieldID,
-                               const ExtrapolationDirectionFinder_T & extrapolationDirectionFinder )
+                               const ExtrapolationDirectionFinder_T & extrapolationDirectionFinder,
+                               const bool & enforceNoSlipConstraintAfterExtrapolation = false )
    : blockStorage_( blockStorage ), boundaryHandlingID_( boundaryHandlingID ),
      pdfFieldID_( pdfFieldID ), bodyFieldID_( bodyFieldID ),
+     enforceNoSlipConstraintAfterExtrapolation_( enforceNoSlipConstraintAfterExtrapolation ),
      extrapolationDirectionFinder_( extrapolationDirectionFinder ),
      alternativeReconstructor_( EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationDirectionFinder_T >
       ( blockStorage, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationDirectionFinder ) )
-   { }
+   {
+      if( enforceNoSlipConstraintAfterExtrapolation_ ) {
+         WALBERLA_CHECK((boost::is_same<typename LatticeModel_T::Stencil, stencil::D3Q19>::value),
+                        "Enforcing no-slip constraint after extrapolation currently only works with D3Q19 stencil!");
+      }
+   }
 
    void operator()( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block );
 
@@ -279,6 +298,8 @@ private:
    const BlockDataID pdfFieldID_;
    const BlockDataID bodyFieldID_;
 
+   const bool enforceNoSlipConstraintAfterExtrapolation_;
+
    ExtrapolationDirectionFinder_T extrapolationDirectionFinder_;
    EquilibriumAndNonEquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationDirectionFinder_T > alternativeReconstructor_;
 
@@ -296,12 +317,17 @@ void ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, Extrapolati
 
    uint_t numberOfCellsForExtrapolation = getNumberOfExtrapolationCells( x, y, z, block, extrapolationDirection );
 
-   if( numberOfCellsForExtrapolation < uint_c(2) )
+   if( numberOfCellsForExtrapolation < uint_t(2) )
+   {
+      alternativeReconstructor_( x, y, z, block, extrapolationDirection );
+   }
+   else
    {
-      alternativeReconstructor_( x, y, z, block );
-   } else {
       extrapolatePDFs( x, y, z, block, extrapolationDirection, numberOfCellsForExtrapolation );
-      enforceNoSlipConstraint( x, y, z, block );
+      if( enforceNoSlipConstraintAfterExtrapolation_ )
+      {
+         enforceNoSlipConstraint( x, y, z, block );
+      }
    }
 }
 
@@ -318,9 +344,11 @@ uint_t ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, Extrapola
    WALBERLA_ASSERT_NOT_NULLPTR( pdfField );
    WALBERLA_ASSERT_NOT_NULLPTR( boundaryHandling );
 
+   if( extrapolationDirection == cell_idx_t(0) ) return uint_t(0);
+
    CellInterval localDomain = pdfField->xyzSize();
 
-   uint_t desiredCellsInExtrapolationDirection = uint_c(3);
+   uint_t desiredCellsInExtrapolationDirection = uint_t(3);
 
    for( uint_t numCells = uint_t(1); numCells <= desiredCellsInExtrapolationDirection; ++numCells )
    {
@@ -345,23 +373,22 @@ void ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, Extrapolati
    WALBERLA_ASSERT_NOT_NULLPTR( block );
 
    PdfField_T * pdfField = block->getData< PdfField_T >( pdfFieldID_ );
-
    WALBERLA_ASSERT_NOT_NULLPTR( pdfField );
 
-   if( numberOfCellsForExtrapolation == uint_c(3) )
+   if( numberOfCellsForExtrapolation == uint_t(3) )
    {
       // quadratic normal extrapolation
       for( auto d = LatticeModel_T::Stencil::begin(); d != LatticeModel_T::Stencil::end(); ++d )
       {
-         pdfField->get( x, y, z, d.toIdx() ) = real_c(3) * pdfField->get( x +   extrapolationDirection[0], y +   extrapolationDirection[1], z +   extrapolationDirection[2], d.toIdx() )
-                                             - real_c(3) * pdfField->get( x + 2*extrapolationDirection[0], y + 2*extrapolationDirection[1], z + 2*extrapolationDirection[2], d.toIdx() )
+         pdfField->get( x, y, z, d.toIdx() ) = real_t(3) * pdfField->get( x +   extrapolationDirection[0], y +   extrapolationDirection[1], z +   extrapolationDirection[2], d.toIdx() )
+                                             - real_t(3) * pdfField->get( x + 2*extrapolationDirection[0], y + 2*extrapolationDirection[1], z + 2*extrapolationDirection[2], d.toIdx() )
                                              +             pdfField->get( x + 3*extrapolationDirection[0], y + 3*extrapolationDirection[1], z + 3*extrapolationDirection[2], d.toIdx() );
       }
    } else { // numberOfCellsForExtrapolation == 2
       // linear normal extrapolation
       for( auto d = LatticeModel_T::Stencil::begin(); d != LatticeModel_T::Stencil::end(); ++d )
       {
-         pdfField->get( x, y, z, d.toIdx() ) = real_c(2) * pdfField->get( x +   extrapolationDirection[0], y +   extrapolationDirection[1], z +   extrapolationDirection[2], d.toIdx() )
+         pdfField->get( x, y, z, d.toIdx() ) = real_t(2) * pdfField->get( x +   extrapolationDirection[0], y +   extrapolationDirection[1], z +   extrapolationDirection[2], d.toIdx() )
                                              -             pdfField->get( x + 2*extrapolationDirection[0], y + 2*extrapolationDirection[1], z + 2*extrapolationDirection[2], d.toIdx() );
       }
    }
@@ -371,6 +398,8 @@ template< typename LatticeModel_T, typename BoundaryHandling_T, typename Extrapo
 void ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, ExtrapolationDirectionFinder_T >
 ::enforceNoSlipConstraint( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block )
 {
+   //NOTE: this currently works only for D3Q19 stencils!
+
    WALBERLA_ASSERT_NOT_NULLPTR( block );
 
    PdfField_T *  pdfField  = block->getData< PdfField_T >( pdfFieldID_ );
@@ -385,6 +414,7 @@ void ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, Extrapolati
 
    WALBERLA_ASSERT_NOT_NULLPTR( (*bodyField)(x,y,z) );
    Vector3<real_t> bodyVelocity = (*bodyField)(x,y,z)->velFromWF(cx,cy,cz);
+   WALBERLA_ASSERT( !math::isnan(bodyVelocity) );
 
    // transforms to moment space (see MRT collision model) to set the body's velocity in cell without affecting other moments
    const real_t _1_2  = real_t(1) / real_t(2);
diff --git a/src/pe_coupling/momentum_exchange_method/restoration/all.h b/src/pe_coupling/momentum_exchange_method/restoration/all.h
index 28229f903f5d77edabf9060e755f431767598a6a..f88dacd21d64f60c5068a63b446580b3713b04e9 100644
--- a/src/pe_coupling/momentum_exchange_method/restoration/all.h
+++ b/src/pe_coupling/momentum_exchange_method/restoration/all.h
@@ -23,5 +23,5 @@
 #pragma once
 
 #include "ExtrapolationDirectionFinder.h"
-#include "PDFRestoration.h"
+#include "PDFReconstruction.h"
 #include "Reconstructor.h"
diff --git a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
index 7bd4e73a2e892b2bf4aceb6b0caa3107adf909c4..752d57d302bea9ca299ad22171f38254154cf704 100644
--- a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
@@ -713,7 +713,7 @@ int main( int argc, char **argv )
    {
       pe_coupling::FlagFieldNormalExtrapolationDirectionFinder< BoundaryHandling_T > extrapolationFinder( blocks, boundaryHandlingID );
       typedef pe_coupling::ExtrapolationReconstructor< LatticeModel_T, BoundaryHandling_T, pe_coupling::FlagFieldNormalExtrapolationDirectionFinder<BoundaryHandling_T> > Reconstructor_T;
-      Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID, extrapolationFinder );
+      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" );