diff --git a/src/field/all.h b/src/field/all.h
index 12a69dfe7ba2ba592cd30ccb6636ffcd93ce422d..05ea39cc7534d32cc00dd84dca8cd72ded61119a 100644
--- a/src/field/all.h
+++ b/src/field/all.h
@@ -45,6 +45,7 @@
 #include "allocation/all.h"
 #include "blockforest/all.h"
 #include "communication/all.h"
+#include "distributors/all.h"
 #include "interpolators/all.h"
 #include "iterators/all.h"
 #include "refinement/all.h"
diff --git a/src/field/communication/UniformPullReductionPackInfo.h b/src/field/communication/UniformPullReductionPackInfo.h
new file mode 100644
index 0000000000000000000000000000000000000000..71ceb7e8e6fba381d5623f80e3ca6f9bb2b8ae47
--- /dev/null
+++ b/src/field/communication/UniformPullReductionPackInfo.h
@@ -0,0 +1,153 @@
+//======================================================================================================================
+//
+//  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 UniformPullReductionPackInfo.h
+//! \ingroup field
+//! \author Tobias Schruff <schruff@iww.rwth-aachen.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "communication/UniformPackInfo.h"
+
+#include "core/debug/Debug.h"
+
+#include "field/GhostLayerField.h"
+
+#include "stencil/Directions.h"
+
+
+namespace walberla {
+namespace field {
+namespace communication {
+
+/**
+ * Data packing/unpacking for ghost layer based communication of a single walberla::field::Field
+ *
+ * \ingroup comm
+ *
+ * template ReduceOperation is e.g. std::plus
+ *
+ * This pack info is used to apply a given reduce operation (e.g. +) to the values in the interior of a ghost layer field
+ * together with the values coming from the sender's ghost layer.
+ */
+template< template<typename> class ReduceOperation, typename GhostLayerField_T >
+class UniformPullReductionPackInfo : public walberla::communication::UniformPackInfo
+{
+public:
+   typedef typename GhostLayerField_T::value_type T;
+
+   UniformPullReductionPackInfo( const BlockDataID & bdID ) : bdID_( bdID ), communicateAllGhostLayers_( true ),
+                                 numberOfGhostLayers_( 0 ) {}
+
+   UniformPullReductionPackInfo( const BlockDataID & bdID, const uint_t numberOfGhostLayers ) : bdID_( bdID ),
+                                 communicateAllGhostLayers_( false ), numberOfGhostLayers_(  numberOfGhostLayers ) {}
+
+   virtual ~UniformPullReductionPackInfo() {}
+
+   bool constantDataExchange() const { return mpi::BufferSizeTrait<T>::constantSize; }
+   bool threadsafeReceiving()  const { return true; }
+
+   void unpackData(IBlock * receiver, stencil::Direction dir, mpi::RecvBuffer & buffer);
+
+   void communicateLocal(const IBlock * sender, IBlock * receiver, stencil::Direction dir);
+
+protected:
+   void packDataImpl(const IBlock * sender, stencil::Direction dir, mpi::SendBuffer & outBuffer) const;
+
+   uint_t numberOfGhostLayersToCommunicate( const GhostLayerField_T * const field ) const;
+
+   const  BlockDataID bdID_;
+   bool   communicateAllGhostLayers_;
+   uint_t numberOfGhostLayers_;
+   ReduceOperation<T> op_;
+};
+
+
+
+template< template<typename> class ReduceOperation, typename GhostLayerField_T >
+void UniformPullReductionPackInfo<ReduceOperation, GhostLayerField_T>::unpackData( IBlock * receiver, stencil::Direction dir, mpi::RecvBuffer & buffer )
+{
+   GhostLayerField_T * f = receiver->getData< GhostLayerField_T >( bdID_ );
+   WALBERLA_ASSERT_NOT_NULLPTR(f);
+
+   cell_idx_t nrOfGhostLayers = cell_idx_c( numberOfGhostLayersToCommunicate( f ) );
+
+   T buf( 0 );
+   for( auto i = f->beginSliceBeforeGhostLayer( dir, nrOfGhostLayers ); i != f->end(); ++i ) {
+      buffer >> buf;
+      *i = op_( *i, buf );
+   }
+}
+
+
+
+template< template<typename> class ReduceOperation, typename GhostLayerField_T>
+void UniformPullReductionPackInfo<ReduceOperation, GhostLayerField_T>::communicateLocal( const IBlock * sender, IBlock * receiver, stencil::Direction dir )
+{
+   const GhostLayerField_T * sf = sender  ->getData< GhostLayerField_T >( bdID_ );
+         GhostLayerField_T * rf = receiver->getData< GhostLayerField_T >( bdID_ );
+
+   WALBERLA_ASSERT_EQUAL(sf->xSize(), rf->xSize());
+   WALBERLA_ASSERT_EQUAL(sf->ySize(), rf->ySize());
+   WALBERLA_ASSERT_EQUAL(sf->zSize(), rf->zSize());
+
+   uint_t nrOfGhostLayers = numberOfGhostLayersToCommunicate( sf );
+   auto srcIter = sf->beginGhostLayerOnly( nrOfGhostLayers, dir );
+   auto dstIter = rf->beginSliceBeforeGhostLayer( stencil::inverseDir[dir], cell_idx_c( nrOfGhostLayers ) );
+
+   while( srcIter != sf->end() ) {
+      *dstIter = op_( *srcIter, *dstIter );
+      ++srcIter;
+      ++dstIter;
+   }
+
+   WALBERLA_ASSERT(srcIter == sf->end() && dstIter == rf->end());
+}
+
+
+
+template< template<typename> class ReduceOperation, typename GhostLayerField_T>
+void UniformPullReductionPackInfo<ReduceOperation, GhostLayerField_T>::packDataImpl( const IBlock * sender, stencil::Direction dir, mpi::SendBuffer & outBuffer ) const
+{
+   const GhostLayerField_T * f = sender->getData< GhostLayerField_T >( bdID_ );
+   WALBERLA_ASSERT_NOT_NULLPTR(f);
+
+   uint_t nrOfGhostLayers = numberOfGhostLayersToCommunicate( f );
+   for( auto i = f->beginGhostLayerOnly( nrOfGhostLayers, dir ); i != f->end(); ++i )
+      outBuffer << *i;
+}
+
+
+
+template< template<typename> class ReduceOperation, typename GhostLayerField_T>
+uint_t UniformPullReductionPackInfo<ReduceOperation, GhostLayerField_T>::numberOfGhostLayersToCommunicate( const GhostLayerField_T * const field ) const
+{
+   if( communicateAllGhostLayers_ )
+   {
+      return field->nrOfGhostLayers();
+   }
+   else
+   {
+      WALBERLA_ASSERT_LESS_EQUAL( numberOfGhostLayers_, field->nrOfGhostLayers() );
+      return numberOfGhostLayers_;
+   }
+}
+
+} // namespace communication
+} // namespace field
+} // namespace walberla
diff --git a/src/field/communication/all.h b/src/field/communication/all.h
index 55782d6e21d8503f473b12a58e10f74ac1d67a8b..a98de888790240fbf77b0454adaf5266ede9ec6b 100644
--- a/src/field/communication/all.h
+++ b/src/field/communication/all.h
@@ -27,3 +27,4 @@
 #include "MPIDatatypes.h"
 #include "ReducePackInfo.h"
 #include "UniformMPIDatatypeInfo.h"
+#include "UniformPullReductionPackInfo.h"
diff --git a/src/field/distributors/DistributorCreators.h b/src/field/distributors/DistributorCreators.h
new file mode 100644
index 0000000000000000000000000000000000000000..e87b907cebac9b1aabf92dba0be34f730e16fffb
--- /dev/null
+++ b/src/field/distributors/DistributorCreators.h
@@ -0,0 +1,108 @@
+//======================================================================================================================
+//
+//  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 DistributorCreators.h
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/Set.h"
+
+#include "blockforest/BlockDataHandling.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+namespace walberla {
+namespace field {
+
+//**********************************************************************************************************************
+/*! DistributorCreators
+ *
+ * \ingroup field
+ *
+ * Distributor_T: A distributor that has a constructor
+ *  ( const weak_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block, const BaseField_T & baseField,
+ *    const FlagField_T & flagField, const flag_t & evaluationMask )
+ * and distribution functions:
+ *  template< typename ForwardIterator_T > inline void distribute( const Vector3<real_t> & position, ForwardIterator_T distributeValueBegin )
+ *  template< typename ForwardIterator_T > inline void distribute( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T distributeValueBegin )
+ *
+ * See NearestNeighborDistributor for an example implementation.
+ *
+ * A distributor is aware of the flag field (FlagField_T) and distributes values only to cells flagged by a given mask.
+ *
+ * Distributors are used to spread a given value to the corresponding destination field.
+ * E.g. if a certain force has to be applied at some specific position onto the fluid, a distributor can be used
+ * to do so by distributing this force value (and conservation fo this force value is ensured) onto the force field.
+ *
+ */
+//**********************************************************************************************************************
+template< typename Distributor_T, typename FlagField_T >
+class DistributorHandling : public blockforest::AlwaysInitializeBlockDataHandling< Distributor_T >
+{
+public:
+
+   DistributorHandling( const weak_ptr<StructuredBlockStorage> & blockStorage,
+                        const BlockDataID & distributionDestinationFieldID,
+                        const ConstBlockDataID & flagFieldID,
+                        const Set< FlagUID > & cellsToEvaluate ) :
+   blockStorage_( blockStorage ), distributionDestinationFieldID_( distributionDestinationFieldID ), flagFieldID_( flagFieldID ), cellsToEvaluate_( cellsToEvaluate )
+   {}
+
+   Distributor_T * initialize( IBlock * const block )
+   {
+      typedef typename Distributor_T::BaseField_T DistributionDestinationField_T;
+      typedef typename FlagField_T::flag_t flag_t;
+      DistributionDestinationField_T * distributionDestinationField = block->getData< DistributionDestinationField_T >( distributionDestinationFieldID_ );
+      const FlagField_T * flagField = block->getData< FlagField_T >( flagFieldID_ );
+
+      WALBERLA_ASSERT_NOT_NULLPTR( distributionDestinationField );
+      WALBERLA_ASSERT_NOT_NULLPTR( flagField );
+
+      const flag_t evaluationMask = flagField->getMask( cellsToEvaluate_ );
+
+      return new Distributor_T( blockStorage_, *block, *distributionDestinationField, *flagField, evaluationMask );
+   }
+
+private:
+
+   weak_ptr<StructuredBlockStorage> blockStorage_;
+   BlockDataID distributionDestinationFieldID_;
+   ConstBlockDataID flagFieldID_;
+   Set< FlagUID > cellsToEvaluate_;
+   
+}; // class DistributorHandling
+
+
+
+template< typename Distributor_T, typename FlagField_T  >
+inline BlockDataID addDistributor( const shared_ptr< StructuredBlockStorage > & blocks,
+                                   const BlockDataID & distributionDestinationFieldID,
+                                   const ConstBlockDataID & flagFieldID,
+                                   const Set< FlagUID > & cellsToEvaluate,
+                                   const std::string & identifier = std::string(),
+                                   const Set<SUID> & requiredSelectors     = Set<SUID>::emptySet(),
+                                   const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() )
+{
+   return blocks->addBlockData( make_shared< DistributorHandling< Distributor_T, FlagField_T > >( blocks, distributionDestinationFieldID, flagFieldID, cellsToEvaluate ), identifier, requiredSelectors, incompatibleSelectors );
+}
+
+
+} // namespace field
+} // namespace walberla
diff --git a/src/field/distributors/KernelDistributor.h b/src/field/distributors/KernelDistributor.h
new file mode 100644
index 0000000000000000000000000000000000000000..1aed2f2ea3c801998114723c8a8ac127d708ac33
--- /dev/null
+++ b/src/field/distributors/KernelDistributor.h
@@ -0,0 +1,159 @@
+//======================================================================================================================
+//
+//  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 KernelDistributor.h
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/debug/Debug.h"
+#include "core/math/Vector3.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/interpolators/KernelFieldInterpolator.h"
+#include "field/GhostLayerField.h"
+
+#include <vector>
+
+namespace walberla {
+namespace field {
+
+/*! Distributor for walberla::field::GhostLayerField with kernel strategy
+ *
+ * \ingroup field
+ *
+ * This distributor uses a smoothed dirac kernel function to distribute values to the field at the given position.
+ * The applied weights are given in the namespace field::kernelweights.
+ * Needs the full neighborhood of the containing cell, i.e. 27 cells.
+ * Never construct this distributor directly, but use the functionality from DistributorCreators.h instead.
+ */
+template< typename Field_T, typename FlagField_T >
+class KernelDistributor
+{
+public:
+
+   static const uint_t F_SIZE = Field_T::F_SIZE;
+
+   typedef Field_T                                BaseField_T;
+   typedef typename FlagField_T::flag_t           flag_t;
+   typedef KernelDistributor<Field_T,FlagField_T> OwnType;
+
+   KernelDistributor( const weak_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block,
+                      BaseField_T & baseField, const FlagField_T & flagField,
+                      const flag_t & evaluationMask )
+   : blockStorage_( blockStorage ), block_( block ), baseField_( baseField ), flagField_( flagField ), evaluationMask_( evaluationMask )
+   {
+      WALBERLA_ASSERT(baseField.nrOfGhostLayers() > uint_t(0), "field for kernel distribution needs at least one ghost layer");
+   }
+
+   inline bool operator==( const OwnType & other ){ return baseField_ == other.baseField_; }
+
+   template< typename ForwardIterator_T >
+   inline void distribute( const Vector3<real_t> & position, ForwardIterator_T distributeValueBegin )
+   {
+      distribute( position[0], position[1], position[2], distributeValueBegin );
+   }
+
+   template< typename ForwardIterator_T >
+   inline void distribute( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T distributeValueBegin )
+   {
+      WALBERLA_ASSERT(block_.getAABB().contains(x,y,z),
+                      "Distribution position <" << x << ", " << y << ", " << z << "> is not contained inside the block of this distributor with AABB " << block_.getAABB() << " !");
+
+      WALBERLA_CHECK( !blockStorage_.expired() );
+      auto blockStorage = blockStorage_.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(blockStorage);
+
+      Cell centerCell = blockStorage->getBlockLocalCell( block_, x, y, z );
+
+      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_ ) );
+      
+      const uint_t neighborhoodSize = cell_idx_t(1);
+
+      CellInterval cellNeighborhood( centerCell[0] - cell_idx_c(neighborhoodSize), centerCell[1] - cell_idx_c(neighborhoodSize), centerCell[2] - cell_idx_c(neighborhoodSize),
+                                     centerCell[0] + cell_idx_c(neighborhoodSize), centerCell[1] + cell_idx_c(neighborhoodSize), centerCell[2] + cell_idx_c(neighborhoodSize) );
+
+      const uint_t kernelSizeOneDirection = uint_t(2) * neighborhoodSize + uint_t(1);
+      std::vector<real_t> weights( kernelSizeOneDirection * kernelSizeOneDirection * kernelSizeOneDirection, real_t(0) );
+
+      uint_t counter = uint_t(0);
+      real_t sumOfWeights = real_t(0);
+      real_t sumOfWeightsUnavailable = real_t(0);
+
+      // get distribution weights and count available cells in surrounding cells
+      for( auto cellIt = cellNeighborhood.begin(); cellIt != cellNeighborhood.end(); ++cellIt )
+      {
+         Vector3<real_t> curCellCenter = blockStorage->getBlockLocalCellCenter( block_, *cellIt );
+         if( flagField_.isPartOfMaskSet( *cellIt, evaluationMask_ ) )
+         {
+            weights[counter] = kernelweights::kernelWeightFunction( x, y, z, curCellCenter[0], curCellCenter[1], curCellCenter[2], dx, dy, dz );
+            sumOfWeights += weights[counter];
+         }
+         else
+         {
+            weights[counter] = real_t(0);
+            sumOfWeightsUnavailable += kernelweights::kernelWeightFunction( x, y, z, curCellCenter[0], curCellCenter[1], curCellCenter[2], dx, dy, dz );
+         }
+         ++counter;
+      }
+
+      // check if at least one cell was available, to prevent division by 0
+      if( sumOfWeights <= real_t(0) )
+         return;
+
+      // scale domain weights if some non-domain cells are in neighborhood
+      const real_t scalingFactor = real_t(1) + sumOfWeightsUnavailable / sumOfWeights;
+
+      // distribute the values to the neighboring domain cells with the corresponding (scaled) weighting
+      counter = uint_t(0);
+      for( auto cellIt = cellNeighborhood.begin(); cellIt != cellNeighborhood.end(); ++cellIt )
+      {
+         if ( weights[counter] > real_t(0) )
+         {
+            addWeightedCellValue( distributeValueBegin, *cellIt, scalingFactor * weights[counter] );
+         }
+         ++counter;
+      }
+   }
+
+private:
+
+   template< typename ForwardIterator_T >
+   void addWeightedCellValue( ForwardIterator_T distributeValueBegin, const Cell & curCell, const real_t & weighting )
+   {
+      for( uint_t f = uint_t(0); f < F_SIZE; ++f )
+      {
+         baseField_( curCell, f) += weighting * (*distributeValueBegin);
+         ++distributeValueBegin;
+      }
+   }
+
+   weak_ptr<StructuredBlockStorage> blockStorage_;
+   const IBlock & block_;
+   BaseField_T & baseField_;
+   const FlagField_T & flagField_;
+   flag_t evaluationMask_;
+
+};
+
+
+} // namespace field
+} // namespace walberla
diff --git a/src/field/distributors/NearestNeighborDistributor.h b/src/field/distributors/NearestNeighborDistributor.h
new file mode 100644
index 0000000000000000000000000000000000000000..6aa4ba2958864e0be6de91c487566b787f5eaaf9
--- /dev/null
+++ b/src/field/distributors/NearestNeighborDistributor.h
@@ -0,0 +1,142 @@
+//======================================================================================================================
+//
+//  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 NearestNeighborDistributor.h
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/math/Vector3.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/GhostLayerField.h"
+
+#include <numeric>
+#include <vector>
+
+namespace walberla {
+namespace field {
+
+/*! Distributor for walberla::field::Field with nearest neighbor strategy
+ *
+ * \ingroup field
+ *
+ * This distributor distributes the given value to the nearest cell, flagged in the evaluation mask, of the given position.
+ * This is usually the cell that contains the distribution position.
+ * If this cell is a cell that is not within the evaluation mask, the direct neighborhood (8 cells) will be searched for an available cell.
+ * Never construct this distributor directly, but use the functionality from DistributorCreators.h instead.
+ */
+template< typename Field_T, typename FlagField_T >
+class NearestNeighborDistributor
+{
+public:
+
+   static const uint_t F_SIZE = Field_T::F_SIZE;
+
+   typedef Field_T                                         BaseField_T;
+   typedef typename FlagField_T::flag_t                    flag_t;
+   typedef NearestNeighborDistributor<Field_T,FlagField_T> OwnType;
+
+   NearestNeighborDistributor( const weak_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block,
+                               BaseField_T & baseField, const FlagField_T & flagField,
+                               const flag_t & evaluationMask )
+   : blockStorage_( blockStorage ), block_( block ), baseField_( baseField ), flagField_( flagField ), evaluationMask_( evaluationMask )
+   {}
+
+   inline bool operator==( const OwnType & other ){ return baseField_ == other.baseField_; }
+
+   template< typename ForwardIterator_T >
+   inline void distribute( const Vector3<real_t> & position, ForwardIterator_T distributeValueBegin )
+   {
+      distribute( position[0], position[1], position[2], distributeValueBegin );
+   }
+
+   template< typename ForwardIterator_T >
+   inline void distribute( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T distributeValueBegin )
+   {
+      WALBERLA_ASSERT(block_.getAABB().contains(x,y,z),
+                      "Distribution position <" << x << ", " << y << ", " << z << "> is not contained inside the block of this distributor with AABB " << block_.getAABB() << " !");
+
+      WALBERLA_CHECK( !blockStorage_.expired() );
+      auto blockStorage = blockStorage_.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(blockStorage);
+
+      Cell nearestCell = blockStorage->getBlockLocalCell( block_, x, y, z );
+
+      if( flagField_.isPartOfMaskSet( nearestCell, evaluationMask_ ) )
+      {
+         for( uint_t f = uint_t(0); f < F_SIZE; ++f )
+         {
+            baseField_( nearestCell, f) += *distributeValueBegin;
+            ++distributeValueBegin;
+         }
+         return;
+      }
+      else
+      {
+         // look for available cell in direct neighborhood
+
+         CellInterval fieldXYZSize = baseField_.xyzSize();
+
+         Vector3<real_t> nearestCellCenter = blockStorage->getBlockLocalCellCenter( block_, nearestCell );
+         const cell_idx_t xNeighbor = cell_idx_c( floor( x - nearestCellCenter[0] ) );
+         const cell_idx_t yNeighbor = cell_idx_c( floor( y - nearestCellCenter[1] ) );
+         const cell_idx_t zNeighbor = cell_idx_c( floor( z - nearestCellCenter[2] ) );
+
+         const cell_idx_t xMin = nearestCell.x() + xNeighbor;
+         const cell_idx_t yMin = nearestCell.y() + yNeighbor;
+         const cell_idx_t zMin = nearestCell.z() + zNeighbor;
+
+         for( cell_idx_t zC = zMin; zC <= zMin + cell_idx_t(1); ++zC)
+         {
+            for( cell_idx_t yC = yMin; yC <= yMin + cell_idx_t(1); ++yC)
+            {
+               for( cell_idx_t xC = xMin; xC <= xMin + cell_idx_t(1); ++xC)
+               {
+                  Cell curCell(xC,yC,zC);
+                  if( fieldXYZSize.contains(curCell) )
+                  {
+                     if (flagField_.isPartOfMaskSet(curCell, evaluationMask_))
+                     {
+                        for (uint_t f = uint_t(0); f < F_SIZE; ++f)
+                        {
+                           baseField_(curCell, f) += *distributeValueBegin;
+                           ++distributeValueBegin;
+                        }
+                        return;
+                     }
+                  }
+               }
+            }
+         }
+      }
+   }
+
+private:
+
+   weak_ptr<StructuredBlockStorage> blockStorage_;
+   const IBlock & block_;
+   BaseField_T & baseField_;
+   const FlagField_T & flagField_;
+   flag_t evaluationMask_;
+
+};
+
+} // namespace field
+} // namespace walberla
diff --git a/src/field/distributors/all.h b/src/field/distributors/all.h
new file mode 100644
index 0000000000000000000000000000000000000000..8a31a22b6e70b8fd6d7647a6a0840220f8d02386
--- /dev/null
+++ b/src/field/distributors/all.h
@@ -0,0 +1,26 @@
+//======================================================================================================================
+//
+//  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 all.h
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "DistributorCreators.h"
+#include "KernelDistributor.h"
+#include "NearestNeighborDistributor.h"
diff --git a/src/field/interpolators/FieldInterpolatorCreators.h b/src/field/interpolators/FieldInterpolatorCreators.h
new file mode 100644
index 0000000000000000000000000000000000000000..560df7879d03c5b4180b39094bf0f991614a7299
--- /dev/null
+++ b/src/field/interpolators/FieldInterpolatorCreators.h
@@ -0,0 +1,107 @@
+//======================================================================================================================
+//
+//  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 FieldInterpolatorCreators.h
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/Set.h"
+
+#include "blockforest/BlockDataHandling.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+namespace walberla {
+namespace field {
+
+//**********************************************************************************************************************
+/*! FieldInterpolatorCreators
+ *
+ * \ingroup field
+ *
+ * Interpolator_T: A field interpolator that has a constructor
+ *  ( const weak_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block, const BaseField_T & baseField,
+ *    const FlagField_T & flagField, const flag_t & evaluationMask )
+ * and getter functions:
+ *  template< typename ForwardIterator_T > inline void get( const Vector3<real_t> & position, ForwardIterator_T interpolationResultBegin )
+ *  template< typename ForwardIterator_T > inline void get( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T interpolationResultBegin )
+ *
+ * See TrilinearFieldInterpolator for an example implementation.
+ *
+ * A field interpolator is aware of the flag field (FlagField_T) and uses only values that contain flags from a given mask.
+ *
+ * Field Interpolators can be used to sample the underlying field at certain positions.
+ * E.g. the fluid velocity can be interpolated to a desired global position from a velocity field.
+ *
+ */
+//**********************************************************************************************************************
+template< typename Interpolator_T, typename FlagField_T >
+class InterpolatorHandling : public blockforest::AlwaysInitializeBlockDataHandling< Interpolator_T >
+{
+public:
+
+   InterpolatorHandling( const weak_ptr<StructuredBlockStorage> & blockStorage,
+                         const ConstBlockDataID & interpolatedFieldID,
+                         const ConstBlockDataID & flagFieldID,
+                         const Set< FlagUID > & cellsToEvaluate ) :
+   blockStorage_( blockStorage ), interpolatedFieldID_( interpolatedFieldID ), flagFieldID_( flagFieldID ), cellsToEvaluate_( cellsToEvaluate )
+   {}
+
+   Interpolator_T * initialize( IBlock * const block )
+   {
+      typedef typename Interpolator_T::BaseField_T InterpolatedField_T;
+      typedef typename FlagField_T::flag_t flag_t;
+      const InterpolatedField_T * interpolatedField = block->getData< InterpolatedField_T >( interpolatedFieldID_ );
+      const FlagField_T * flagField = block->getData< FlagField_T >( flagFieldID_ );
+
+      WALBERLA_ASSERT_NOT_NULLPTR( interpolatedField );
+      WALBERLA_ASSERT_NOT_NULLPTR( flagField );
+
+      const flag_t evaluationMask = flagField->getMask( cellsToEvaluate_ );
+
+      return new Interpolator_T( blockStorage_, *block, *interpolatedField, *flagField, evaluationMask );
+   }
+
+private:
+
+   weak_ptr<StructuredBlockStorage> blockStorage_;
+   ConstBlockDataID interpolatedFieldID_;
+   ConstBlockDataID flagFieldID_;
+   Set< FlagUID > cellsToEvaluate_;
+   
+}; // class InterpolatorHandling
+
+
+
+template< typename Interpolator_T, typename FlagField_T  >
+inline BlockDataID addFieldInterpolator( const shared_ptr< StructuredBlockStorage > & blocks,
+                                         const ConstBlockDataID & interpolatedFieldID,
+                                         const ConstBlockDataID & flagFieldID,
+                                         const Set< FlagUID > & cellsToEvaluate,
+                                         const std::string & identifier = std::string(),
+                                         const Set<SUID> & requiredSelectors     = Set<SUID>::emptySet(),
+                                         const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() )
+{
+   return blocks->addBlockData( make_shared< InterpolatorHandling< Interpolator_T, FlagField_T > >( blocks, interpolatedFieldID, flagFieldID, cellsToEvaluate ), identifier, requiredSelectors, incompatibleSelectors );
+}
+
+
+} // namespace field
+} // namespace walberla
diff --git a/src/field/interpolators/KernelFieldInterpolator.h b/src/field/interpolators/KernelFieldInterpolator.h
new file mode 100644
index 0000000000000000000000000000000000000000..eff9e888a6fa03df0d6b06cb015beb158ee5dd73
--- /dev/null
+++ b/src/field/interpolators/KernelFieldInterpolator.h
@@ -0,0 +1,274 @@
+//======================================================================================================================
+//
+//  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 KernelFieldInterpolator.h
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/debug/Debug.h"
+#include "core/math/Vector3.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/FlagField.h"
+
+#include "stencil/D3Q27.h"
+
+#include <numeric>
+
+namespace walberla {
+namespace field {
+
+namespace kernelweights
+{
+
+// corresponds to the smoothed dirac delta function from Roma et al - An Adaptive Version of the Immersed Boundary Method
+// f(r) != 0 for abs(r) < 1.5 -> requires three neighborhood cells
+real_t smoothedDeltaFunction( const real_t & r )
+{
+   real_t rAbs = std::fabs(r);
+   if( rAbs <= real_t(0.5) )
+   {
+      return (real_t(1) + std::sqrt( - real_t(3) * r * r + real_t(1) ) ) * (real_t(1) / real_t(3));
+   } else if ( rAbs < real_t(1.5) )
+   {
+      return (real_t(5) - real_t(3) * rAbs - std::sqrt( - real_t(3) * ( real_t(1) - rAbs ) * ( real_t(1) - rAbs ) + real_t(1) ) ) * ( real_t(1) / real_t(6) );
+   } else
+   {
+      return real_t(0);
+   }
+
+}
+
+// X: Lagrangian position, x: Eulerian position (usually cell center), global coordinates
+// dx, dy, dz: mesh spacing
+real_t kernelWeightFunction( const real_t & X, const real_t & Y, const real_t & Z,
+                             const real_t & x, const real_t & y, const real_t & z,
+                             const real_t & dx = real_t(1), const real_t & dy = real_t(1), const real_t & dz = real_t(1) )
+{
+   return smoothedDeltaFunction( ( X - x ) / dx ) * smoothedDeltaFunction( ( Y - y ) / dy ) * smoothedDeltaFunction( ( Z - z ) / dz );
+}
+
+// X: Lagrangian position, x: Eulerian position (usually cell center), global coordinates
+// dx, dy, dz: mesh spacing
+real_t kernelWeightFunction( const Vector3<real_t> & X, const Vector3<real_t> & x,
+                             const real_t & dx = real_t(1), const real_t & dy = real_t(1), const real_t & dz = real_t(1) )
+{
+   return smoothedDeltaFunction( ( X[0] - x[0] ) / dx ) * smoothedDeltaFunction( ( X[1] - x[1] ) / dy ) * smoothedDeltaFunction( ( X[2] - x[2] ) / dz );
+}
+
+} // namespace kernelweights
+
+
+/*! Interpolator for walberla::field::GhostLayerField with kernel strategy
+ *
+ * \ingroup field
+ *
+ * This interpolator uses a smoothed dirac kernel function to interpolate values.
+ * The applied weights are given in the namespace field::kernelweights.
+ * Needs the full neighborhood of the containing cell, i.e. 27 cells.
+ * Never construct this interpolator directly, but use the functionality from the FieldInterpolatorCreator.h instead.
+ */
+template< typename Field_T, typename FlagField_T >
+class KernelFieldInterpolator
+{
+public:
+
+   static const uint_t F_SIZE = Field_T::F_SIZE;
+
+   typedef Field_T                                      BaseField_T;
+   typedef typename FlagField_T::flag_t                 flag_t;
+   typedef KernelFieldInterpolator<Field_T,FlagField_T> OwnType;
+
+   KernelFieldInterpolator( const weak_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block,
+                            const BaseField_T & baseField, const FlagField_T & flagField,
+                            const flag_t & evaluationMask )
+   : blockStorage_( blockStorage ), block_( block ), baseField_( baseField ), flagField_( flagField ), evaluationMask_( evaluationMask )
+   {
+      WALBERLA_ASSERT(baseField.nrOfGhostLayers() > uint_t(0), "field for kernel interpolation needs at least one ghost layer");
+   }
+
+
+   inline bool operator==( const OwnType & other ){ return baseField_ == other.baseField_; }
+
+   template< typename ForwardIterator_T >
+   inline void get( const Vector3<real_t> & position, ForwardIterator_T interpolationResultBegin )
+   {
+      get( position[0], position[1], position[2], interpolationResultBegin );
+   }
+
+   template< typename ForwardIterator_T >
+   inline void get( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T interpolationResultBegin )
+   {
+
+      WALBERLA_ASSERT(block_.getAABB().contains(x,y,z),
+                      "Interpolation position <" << x << ", " << y << ", " << z << "> is not contained inside the block of this interpolator with AABB " << block_.getAABB() << " !");
+
+      WALBERLA_CHECK( !blockStorage_.expired() );
+      auto blockStorage = blockStorage_.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(blockStorage);
+
+      Cell centerCell = blockStorage->getBlockLocalCell( block_, x, y, z );
+
+      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_ ) );
+
+      const uint_t neighborhoodSize = uint_t(1);
+
+      CellInterval cellNeighborhood( centerCell[0] - cell_idx_c(neighborhoodSize), centerCell[1] - cell_idx_c(neighborhoodSize), centerCell[2] - cell_idx_c(neighborhoodSize),
+                                     centerCell[0] + cell_idx_c(neighborhoodSize), centerCell[1] + cell_idx_c(neighborhoodSize), centerCell[2] + cell_idx_c(neighborhoodSize) );
+
+      const uint_t kernelSizeOneDirection = uint_t(2) * neighborhoodSize + uint_t(1);
+      // store the calculated weighting factors of the available cells, i.e. cells flagged by the evaluation mask
+      std::vector<real_t> weights( kernelSizeOneDirection * kernelSizeOneDirection * kernelSizeOneDirection, real_t(0) );
+
+      // store the calculated weighting factors of the UNavailable cells, i.e. cells not included in the evaluation mask
+      std::vector<real_t> weightsUnavailable( kernelSizeOneDirection * kernelSizeOneDirection * kernelSizeOneDirection, real_t(0) );
+
+      cell_idx_t cellIdx0x = cellNeighborhood.xMin();
+      cell_idx_t cellIdx0y = cellNeighborhood.yMin();
+      cell_idx_t cellIdx0z = cellNeighborhood.zMin();
+      uint_t nx = kernelSizeOneDirection;
+      uint_t nxy = kernelSizeOneDirection * kernelSizeOneDirection;
+      Vector3<real_t> cellCenter0 = blockStorage->getBlockLocalCellCenter( block_, Cell(cellIdx0x, cellIdx0y, cellIdx0z) ); // = cell in neighborhood with smallest x-, y-, and z-indices
+
+      // calculate kernel weights of all cells in neighborhood
+      for( uint_t k = uint_t(0); k < kernelSizeOneDirection; ++k)
+      {
+         for( uint_t j = uint_t(0); j < kernelSizeOneDirection; ++j)
+         {
+            for( uint_t i = uint_t(0); i < kernelSizeOneDirection; ++i)
+            {
+               Vector3<real_t> curCellCenter( cellCenter0[0] + real_c(i) * dx, cellCenter0[1] + real_c(j) * dy, cellCenter0[2] + real_c(k) * dz);
+               if( flagField_.isPartOfMaskSet( cellIdx0x + cell_idx_c(i), cellIdx0y + cell_idx_c(j), cellIdx0z + cell_idx_c(k), evaluationMask_ ) )
+               {
+                  weights[k*nxy+j*nx+i] = kernelweights::kernelWeightFunction( x, y, z, curCellCenter[0], curCellCenter[1], curCellCenter[2], dx, dy, dz );
+               } 
+               else
+               {
+                  weightsUnavailable[k*nxy+j*nx+i] = kernelweights::kernelWeightFunction( x, y, z, curCellCenter[0], curCellCenter[1], curCellCenter[2], dx, dy, dz );
+               }
+            }
+         }
+      }
+
+/*
+      // zero entries in weights array (and non-zero entries in weightsUnavailable) indicate non-fluid nodes
+      // weight correction to incorporate extrapolation for kernel interpolation, center element can not be redistributed
+      // without this part, the kernel interpolator can not extrapolate values
+      // this is based on experimental findings by the author and general correctness is not ensured, thus it is commented out
+      for( auto stenIt = stencil::D3Q27::beginNoCenter(); stenIt != stencil::D3Q27::end(); ++stenIt )
+      {
+         int cx = stenIt.cx();
+         int cy = stenIt.cy();
+         int cz = stenIt.cz();
+         int i = cx + 1;
+         int j = cy + 1;
+         int k = cz + 1;
+
+         if( weightsUnavailable[k*nxy+j*nx+i] > real_t(0) )
+         {
+            // check if we can distribute the non-fluid weight to nearby fluid weights that lie in one line inside the neighborhood
+            // it should generally not matter in which direction it is distributed
+            // check x direction
+            if( weights[k*nxy+j*nx+i-cx] > real_t(0) && weights[k*nxy+j*nx+i-2*cx] > real_t(0) )
+            {
+               weights[k*nxy+j*nx+i-cx] += real_t(2)*weightsUnavailable[k*nxy+j*nx+i];
+               weights[k*nxy+j*nx+i-2*cx] -= weightsUnavailable[k*nxy+j*nx+i];
+               weightsUnavailable[k*nxy+j*nx+i] = real_t(0);
+               continue;
+            }
+            // check y direction
+            if( weights[k*nxy+(j-cy)*nx+i] > real_t(0) && weights[k*nxy+(j-2*cy)*nx+i] > real_t(0) )
+            {
+               weights[k*nxy+(j-cy)*nx+i] += real_t(2)*weightsUnavailable[k*nxy+j*nx+i];
+               weights[k*nxy+(j-2*cy)*nx+i] -= weightsUnavailable[k*nxy+j*nx+i];
+               weightsUnavailable[k*nxy+j*nx+i] = real_t(0);
+               continue;
+            }
+            // check z direction
+            if( weights[(k-cz)*nxy+j*nx+i] > real_t(0) && weights[(k-2*cz)*nxy+j*nx+i] > real_t(0) )
+            {
+               weights[(k-cz)*nxy+j*nx+i] += real_t(2)*weightsUnavailable[k*nxy+j*nx+i];
+               weights[(k-2*cz)*nxy+j*nx+i] -= weightsUnavailable[k*nxy+j*nx+i];
+               weightsUnavailable[k*nxy+j*nx+i] = real_t(0);
+               continue;
+            }
+         }
+      }
+*/
+      // scale available weights by the total amount of unavailable weights such that afterwards sum over all weights is 1
+      const real_t sumOfWeightsUnavailable = std::accumulate(weightsUnavailable.begin(), weightsUnavailable.end(), real_t(0) );
+      const real_t sumOfWeights = std::accumulate(weights.begin(), weights.end(), real_t(0) );
+
+      // check if at least one cell was available, to prevent division by 0
+      if( sumOfWeights <= real_t(0) )
+         return;
+
+      const real_t scalingFactor = real_t(1) + sumOfWeightsUnavailable / sumOfWeights;
+
+      for( auto weightIt = weights.begin(); weightIt != weights.end(); ++weightIt )
+      {
+         *weightIt *= scalingFactor;
+      }
+
+      // interpolate value to interpolation position using the previously calculated weights
+      // values to not be included have a weight of 0
+      for( uint_t k = uint_t(0); k < kernelSizeOneDirection; ++k)
+      {
+         for( uint_t j = uint_t(0); j < kernelSizeOneDirection; ++j)
+         {
+            for( uint_t i = uint_t(0); i < kernelSizeOneDirection; ++i)
+            {
+               if ( weights[k*nxy+j*nx+i] > real_t(0) )
+               {
+                  Cell curCell( cellIdx0x + cell_idx_c(i), cellIdx0y + cell_idx_c(j), cellIdx0z + cell_idx_c(k) );
+                  addWeightedCellValue( interpolationResultBegin, curCell, weights[k*nxy+j*nx+i] );
+               }
+            }
+         }
+      }
+   }
+
+private:
+
+   template< typename ForwardIterator_T >
+   void addWeightedCellValue( ForwardIterator_T interpolationResultBegin, const Cell & curCell, const real_t & weighting )
+   {
+      for( uint_t f = uint_t(0); f < F_SIZE; ++f )
+      {
+         WALBERLA_ASSERT( !math::isnan(baseField_( curCell, f)), "NaN found in component " << f << " when interpolating from cell " << curCell );
+
+         *interpolationResultBegin += weighting * baseField_( curCell, f);
+         ++interpolationResultBegin;
+      }
+   }
+
+   weak_ptr<StructuredBlockStorage> blockStorage_;
+   const IBlock & block_;
+   const BaseField_T & baseField_;
+   const FlagField_T & flagField_;
+   flag_t evaluationMask_;
+
+};
+
+
+} // namespace field
+} // namespace walberla
diff --git a/src/field/interpolators/NearestNeighborFieldInterpolator.h b/src/field/interpolators/NearestNeighborFieldInterpolator.h
new file mode 100644
index 0000000000000000000000000000000000000000..45e9110de725173af89a4aed20e1d6290b5278a5
--- /dev/null
+++ b/src/field/interpolators/NearestNeighborFieldInterpolator.h
@@ -0,0 +1,143 @@
+//======================================================================================================================
+//
+//  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 NearestNeighborFieldInterpolator.h
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/debug/Debug.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/FlagField.h"
+
+namespace walberla {
+namespace field {
+
+/*! Interpolator for walberla::field::Field with nearest neighbor strategy
+ *
+ * \ingroup field
+ *
+ * This interpolator obtains the value from the nearest cell, flagged in the evaluation mask, of the interpolation position.
+ * This is usually the cell that contains the interpolation position.
+ * If this cell is a cell that is not within the evaluation mask, the direct neighborhood (8 cells) will be searched for an available cell.
+ * Never construct this interpolator directly, but use the functionality from FieldInterpolatorCreators.h instead.
+ */
+template< typename Field_T, typename FlagField_T >
+class NearestNeighborFieldInterpolator
+{
+public:
+
+   static const uint_t F_SIZE = Field_T::F_SIZE;
+
+   typedef Field_T                                               BaseField_T;
+   typedef typename FlagField_T::flag_t                          flag_t;
+   typedef NearestNeighborFieldInterpolator<Field_T,FlagField_T> OwnType;
+
+   NearestNeighborFieldInterpolator( const weak_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block,
+                                     const BaseField_T & baseField, const FlagField_T & flagField,
+                                     const flag_t & evaluationMask )
+   : blockStorage_( blockStorage ), block_( block ), baseField_( baseField ), flagField_( flagField ), evaluationMask_( evaluationMask )
+   {}
+
+
+   inline bool operator==( const OwnType & other ){ return baseField_ == other.baseField_; }
+
+   template< typename ForwardIterator_T >
+   inline void get( const Vector3<real_t> & position, ForwardIterator_T interpolationResultBegin )
+   {
+      get( position[0], position[1], position[2], interpolationResultBegin );
+   }
+
+   template< typename ForwardIterator_T >
+   inline void get( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T interpolationResultBegin )
+   {
+
+      WALBERLA_ASSERT(block_.getAABB().contains(x,y,z),
+                      "Interpolation position <" << x << ", " << y << ", " << z << "> is not contained inside the block of this interpolator with AABB " << block_.getAABB() << " !");
+
+      WALBERLA_CHECK( !blockStorage_.expired() );
+      auto blockStorage = blockStorage_.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(blockStorage);
+
+      Cell nearestCell = blockStorage->getBlockLocalCell( block_, x, y, z );
+
+      if( flagField_.isPartOfMaskSet( nearestCell, evaluationMask_ ) )
+      {
+         for( uint_t f = uint_t(0); f < F_SIZE; ++f )
+         {
+            WALBERLA_ASSERT( !math::isnan( baseField_( nearestCell, f) ), "NaN found in component " << f << " when interpolating from cell " << nearestCell );
+            *interpolationResultBegin = baseField_( nearestCell, f);
+            ++interpolationResultBegin;
+         }
+      }
+      else
+      {
+         // look for available cell in direct neighborhood
+
+         CellInterval fieldXYZSize = baseField_.xyzSize();
+
+         Vector3<real_t> nearestCellCenter = blockStorage->getBlockLocalCellCenter( block_, nearestCell );
+         const cell_idx_t xNeighbor = cell_idx_c( floor( x - nearestCellCenter[0] ) );
+         const cell_idx_t yNeighbor = cell_idx_c( floor( y - nearestCellCenter[1] ) );
+         const cell_idx_t zNeighbor = cell_idx_c( floor( z - nearestCellCenter[2] ) );
+
+         const cell_idx_t xMin = nearestCell.x() + xNeighbor;
+         const cell_idx_t yMin = nearestCell.y() + yNeighbor;
+         const cell_idx_t zMin = nearestCell.z() + zNeighbor;
+
+         for( cell_idx_t zC = zMin; zC <= zMin + cell_idx_t(1); ++zC)
+         {
+            for( cell_idx_t yC = yMin; yC <= yMin + cell_idx_t(1); ++yC)
+            {
+               for( cell_idx_t xC = xMin; xC <= xMin + cell_idx_t(1); ++xC)
+               {
+                  Cell curCell(xC,yC,zC);
+                  if( fieldXYZSize.contains(curCell) )
+                  {
+                     if (flagField_.isPartOfMaskSet(curCell, evaluationMask_))
+                     {
+                        for (uint_t f = uint_t(0); f < F_SIZE; ++f)
+                        {
+                           WALBERLA_ASSERT(!math::isnan(baseField_(curCell, f)),
+                                           "NaN found in component " << f << " when interpolating from cell "
+                                                                     << curCell);
+                           *interpolationResultBegin = baseField_(curCell, f);
+                           ++interpolationResultBegin;
+                        }
+                        break;
+                     }
+                  }
+               }
+            }
+         }
+      }
+   }
+
+private:
+
+   weak_ptr<StructuredBlockStorage> blockStorage_;
+   const IBlock & block_;
+   const BaseField_T & baseField_;
+   const FlagField_T & flagField_;
+   flag_t evaluationMask_;
+};
+
+} // namespace field
+} // namespace walberla
diff --git a/src/field/interpolators/TrilinearFieldInterpolator.h b/src/field/interpolators/TrilinearFieldInterpolator.h
new file mode 100644
index 0000000000000000000000000000000000000000..48d6b2ce9436bb5a87871a558bf0904c196c4e1d
--- /dev/null
+++ b/src/field/interpolators/TrilinearFieldInterpolator.h
@@ -0,0 +1,192 @@
+//======================================================================================================================
+//
+//  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 TrilinearFieldInterpolator.h
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/debug/Debug.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/FlagField.h"
+
+
+namespace walberla {
+namespace field {
+
+/*! Interpolator for walberla::field::GhostLayerField with trilinear strategy
+ *
+ * \ingroup field
+ *
+ * This interpolator uses trilinear interpolation to obtain the interpolated value from the interpolation position.
+ * If at least one of the cells, that are required for trilinear interpolation, i.e. the 8 closest cells,
+ * is not contained in the evaluation mask, the nearest-neighbor interpolation (see NearestNeighborFieldInterpolator.h)
+ * will be used instead.
+ * Never construct this interpolator directly, but use the functionality from FieldInterpolatorCreators.h instead.
+ */
+template< typename Field_T, typename FlagField_T >
+class TrilinearFieldInterpolator
+{
+public:
+
+   static const uint_t F_SIZE = Field_T::F_SIZE;
+
+   typedef Field_T                                         BaseField_T;
+   typedef typename FlagField_T::flag_t                    flag_t;
+   typedef TrilinearFieldInterpolator<Field_T,FlagField_T> OwnType;
+
+   TrilinearFieldInterpolator( const weak_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block,
+                               const BaseField_T & baseField, const FlagField_T & flagField,
+                               const flag_t & evaluationMask )
+   : blockStorage_( blockStorage ), block_( block ), baseField_( baseField ), flagField_( flagField ), evaluationMask_( evaluationMask ),
+     nearestNeighborInterpolator_( blockStorage, block, baseField, flagField, evaluationMask )
+   {
+      WALBERLA_ASSERT(baseField.nrOfGhostLayers() > uint_t(0), "field for trilinear interpolation needs at least one ghost layer");
+   }
+
+
+   inline bool operator==( const OwnType & other ){ return baseField_ == other.baseField_; }
+
+   template< typename ForwardIterator_T >
+   inline void get( const Vector3<real_t> & position, ForwardIterator_T interpolationResultBegin )
+   {
+      get( position[0], position[1], position[2], interpolationResultBegin );
+   }
+
+   template< typename ForwardIterator_T >
+   inline void get( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T interpolationResultBegin )
+   {
+      WALBERLA_ASSERT(block_.getAABB().contains(x,y,z),
+                      "Interpolation position <" << x << ", " << y << ", " << z << "> is not contained inside the block of this interpolator with AABB " << block_.getAABB() << " !");
+
+      WALBERLA_CHECK( !blockStorage_.expired() );
+      auto blockStorage = blockStorage_.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(blockStorage);
+
+      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_ ) );
+
+      Cell containingCell = blockStorage->getBlockLocalCell( block_, x, y, z );
+      Vector3<real_t> containingCellCenter = blockStorage->getBlockLocalCellCenter( block_, containingCell );
+
+      const cell_idx_t xNeighbor1 = cell_idx_c( floor( x - containingCellCenter[0] ) );
+      const cell_idx_t xNeighbor2 = xNeighbor1 + cell_idx_t(1);
+      const cell_idx_t yNeighbor1 = cell_idx_c( floor( y - containingCellCenter[1] ) );
+      const cell_idx_t yNeighbor2 = yNeighbor1 + cell_idx_t(1);
+      const cell_idx_t zNeighbor1 = cell_idx_c( floor( z - containingCellCenter[2] ) );
+      const cell_idx_t zNeighbor2 = zNeighbor1 + cell_idx_t(1);
+
+      // define the 8 nearest cells required for the trilienar interpolation
+      // the cell 'ccc' is the one with the smallest x-, y-, and z-indices
+      Cell ccc( containingCell.x() + xNeighbor1, containingCell.y() + yNeighbor1, containingCell.z() + zNeighbor1 );
+      Cell hcc( containingCell.x() + xNeighbor2, containingCell.y() + yNeighbor1, containingCell.z() + zNeighbor1 );
+      Cell chc( containingCell.x() + xNeighbor1, containingCell.y() + yNeighbor2, containingCell.z() + zNeighbor1 );
+      Cell hhc( containingCell.x() + xNeighbor2, containingCell.y() + yNeighbor2, containingCell.z() + zNeighbor1 );
+      Cell cch( containingCell.x() + xNeighbor1, containingCell.y() + yNeighbor1, containingCell.z() + zNeighbor2 );
+      Cell hch( containingCell.x() + xNeighbor2, containingCell.y() + yNeighbor1, containingCell.z() + zNeighbor2 );
+      Cell chh( containingCell.x() + xNeighbor1, containingCell.y() + yNeighbor2, containingCell.z() + zNeighbor2 );
+      Cell hhh( containingCell.x() + xNeighbor2, containingCell.y() + yNeighbor2, containingCell.z() + zNeighbor2 );
+
+      // check if all neighboring cells are domain cells
+      if( flagField_.isPartOfMaskSet( ccc, evaluationMask_ ) &&
+          flagField_.isPartOfMaskSet( hcc, evaluationMask_ ) &&
+          flagField_.isPartOfMaskSet( chc, evaluationMask_ ) &&
+          flagField_.isPartOfMaskSet( hhc, evaluationMask_ ) &&
+          flagField_.isPartOfMaskSet( cch, evaluationMask_ ) &&
+          flagField_.isPartOfMaskSet( hch, evaluationMask_ ) &&
+          flagField_.isPartOfMaskSet( chh, evaluationMask_ ) &&
+          flagField_.isPartOfMaskSet( hhh, evaluationMask_ ) )
+      {
+         // trilinear interpolation can be applied
+
+         const real_t inv_totalVolume = real_t(1) / ( dx * dy * dz );
+         Vector3<real_t> cccCellCenter = blockStorage->getBlockLocalCellCenter( block_, ccc );
+
+         // weighting = volume of opposing volume element / total volume
+         real_t weighting(0.0);
+         // x: (c)---dxc--(X)--------------dxh-------------(h)
+         const real_t dxc = x - cccCellCenter[0];
+         const real_t dxh = cccCellCenter[0] + dx - x;
+         const real_t dyc = y - cccCellCenter[1];
+         const real_t dyh = cccCellCenter[1] + dy - y;
+         const real_t dzc = z - cccCellCenter[2];
+         const real_t dzh = cccCellCenter[2] + dz - z;
+
+
+         weighting = dxh * dyh * dzh * inv_totalVolume;
+         addWeightedCellValue( interpolationResultBegin, ccc, weighting );
+
+         weighting = dxc * dyh * dzh * inv_totalVolume;
+         addWeightedCellValue( interpolationResultBegin, hcc, weighting );
+
+         weighting = dxh * dyc * dzh * inv_totalVolume;
+         addWeightedCellValue( interpolationResultBegin, chc, weighting );
+
+         weighting = dxc * dyc * dzh * inv_totalVolume;
+         addWeightedCellValue( interpolationResultBegin, hhc, weighting );
+
+         weighting = dxh * dyh * dzc * inv_totalVolume;
+         addWeightedCellValue( interpolationResultBegin, cch, weighting );
+
+         weighting = dxc * dyh * dzc * inv_totalVolume;
+         addWeightedCellValue( interpolationResultBegin, hch, weighting );
+
+         weighting = dxh * dyc * dzc * inv_totalVolume;
+         addWeightedCellValue( interpolationResultBegin, chh, weighting );
+
+         weighting = dxc * dyc * dzc * inv_totalVolume;
+         addWeightedCellValue( interpolationResultBegin, hhh, weighting );
+
+      }
+      else
+      {
+         // revert to nearest neighbor interpolation
+         nearestNeighborInterpolator_.get( x, y, z, interpolationResultBegin );
+      }
+   }
+
+private:
+
+   template< typename ForwardIterator_T >
+   void addWeightedCellValue( ForwardIterator_T interpolationResultBegin, const Cell & curCell, const real_t & weighting )
+   {
+      for( uint_t f = uint_t(0); f < F_SIZE; ++f )
+      {
+         WALBERLA_ASSERT( !math::isnan(baseField_( curCell, f)), "NaN found in component " << f << " when interpolating from cell " << curCell );
+
+         *interpolationResultBegin += weighting * baseField_( curCell, f);
+         ++interpolationResultBegin;
+      }
+   }
+
+   weak_ptr<StructuredBlockStorage> blockStorage_;
+   const IBlock & block_;
+   const BaseField_T & baseField_;
+   const FlagField_T & flagField_;
+   flag_t evaluationMask_;
+
+   NearestNeighborFieldInterpolator<Field_T,FlagField_T> nearestNeighborInterpolator_;
+
+};
+
+
+} // namespace field
+} // namespace walberla
diff --git a/src/field/interpolators/all.h b/src/field/interpolators/all.h
index 0e7cfa80d337842b0320ea205f0dc109e63094a7..92055c83158120b94f671722f30d2e1d46b01ce4 100644
--- a/src/field/interpolators/all.h
+++ b/src/field/interpolators/all.h
@@ -24,3 +24,8 @@
 
 #include "NearestNeighborInterpolator.h"
 #include "TrilinearInterpolator.h"
+
+#include "FieldInterpolatorCreators.h"
+#include "KernelFieldInterpolator.h"
+#include "NearestNeighborFieldInterpolator.h"
+#include "TrilinearFieldInterpolator.h"
diff --git a/tests/field/CMakeLists.txt b/tests/field/CMakeLists.txt
index 4307644a0dcbbdd7d803f6a91d91b8e02e1c869f..d7534a81375d2e85498a6c8fae5b51c0c5b0a68d 100644
--- a/tests/field/CMakeLists.txt
+++ b/tests/field/CMakeLists.txt
@@ -10,6 +10,9 @@ waLBerla_execute_test( NAME AccuracyEvaluationTest4 COMMAND $<TARGET_FILE:Accura
 waLBerla_compile_test( FILES communication/FieldPackInfoTest.cpp DEPENDS blockforest )
 waLBerla_execute_test( NAME  FieldPackInfoTest )
 
+waLBerla_compile_test( FILES distributors/DistributionTest.cpp)
+waLBerla_execute_test( NAME DistributionTest )
+
 waLBerla_compile_test( FILES FieldTest.cpp )
 waLBerla_execute_test( NAME FieldTest ) 
 
@@ -25,6 +28,9 @@ waLBerla_execute_test( NAME FlagFieldTest )
 waLBerla_compile_test( FILES interpolators/InterpolationTest.cpp)
 waLBerla_execute_test( NAME InterpolationTest )
 
+waLBerla_compile_test( FILES interpolators/FieldInterpolationTest.cpp)
+waLBerla_execute_test( NAME FieldInterpolationTest )
+
 waLBerla_compile_test( FILES adaptors/AdaptorTest.cpp DEPENDS blockforest gui lbm )
 waLBerla_execute_test( NAME AdaptorTest )
 
diff --git a/tests/field/communication/FieldPackInfoTest.cpp b/tests/field/communication/FieldPackInfoTest.cpp
index 45f9b00e51e56ada7bf184b095d49c6b88fbe161..d7b94789b8699ffd1b659b9cfdc192ac3fdac4c0 100644
--- a/tests/field/communication/FieldPackInfoTest.cpp
+++ b/tests/field/communication/FieldPackInfoTest.cpp
@@ -16,6 +16,7 @@
 //! \file FieldPackInfoTest.cpp
 //! \ingroup field
 //! \author Martin Bauer <martin.bauer@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
 //! \brief Tests if a Field is correctly packed into buffers
 //
 //======================================================================================================================
@@ -23,6 +24,7 @@
 #include "field/AddToStorage.h"
 #include "field/GhostLayerField.h"
 #include "field/communication/PackInfo.h"
+#include "field/communication/UniformPullReductionPackInfo.h"
 
 #include "blockforest/Initialization.h"
 
@@ -85,9 +87,60 @@ void testScalarField( IBlock * block, BlockDataID fieldId )
    WALBERLA_CHECK_EQUAL ( field(1,1,+2), 4 );
 }
 
+void testScalarFieldPullReduction( IBlock * block, BlockDataID fieldId )
+{
+   GhostLayerField<int,1> & field = *(block->getData<GhostLayerField<int,1> > (fieldId));
+   field.setWithGhostLayer( 0 );
 
+   WALBERLA_CHECK_EQUAL(field.xSize(), 2);
+   WALBERLA_CHECK_EQUAL(field.ySize(), 2);
+   WALBERLA_CHECK_EQUAL(field.zSize(), 2);
 
+   // initialize the bottom ghost layer cells
+   field(0,0,-1) = 1;
+   field(0,1,-1) = 2;
+   field(1,0,-1) = 3;
+   field(1,1,-1) = 4;
+
+   // initialize the top interior cells
+   field(0,0,1) = 1;
+   field(0,1,1) = 1;
+   field(1,0,1) = 1;
+   field(1,1,1) = 1;
+
+   // communicate periodic from bottom to top with uniform pull scheme
+   field::communication::UniformPullReductionPackInfo<std::plus, GhostLayerField<int,1> > pi1 (fieldId);
+   pi1.communicateLocal( block, block, stencil::B );
+
+   // check values in top ghost layer
+   WALBERLA_CHECK_EQUAL ( field(0,0,2), 0 );
+   WALBERLA_CHECK_EQUAL ( field(0,1,2), 0 );
+   WALBERLA_CHECK_EQUAL ( field(1,0,2), 0 );
+   WALBERLA_CHECK_EQUAL ( field(1,1,2), 0 );
+
+   // check values in top interior cells
+   WALBERLA_CHECK_EQUAL ( field(0,0,1), 2 );
+   WALBERLA_CHECK_EQUAL ( field(0,1,1), 3 );
+   WALBERLA_CHECK_EQUAL ( field(1,0,1), 4 );
+   WALBERLA_CHECK_EQUAL ( field(1,1,1), 5 );
+
+   // communicate periodic from top to bottom with standard form to sync ghost layers
+   field::communication::PackInfo< GhostLayerField<int,1> > pi2 (fieldId);
+   pi2.communicateLocal( block, block, stencil::T );
+
+   // check values in bottom ghost layer
+   WALBERLA_CHECK_EQUAL ( field(0,0,-1), 2 );
+   WALBERLA_CHECK_EQUAL ( field(0,1,-1), 3 );
+   WALBERLA_CHECK_EQUAL ( field(1,0,-1), 4 );
+   WALBERLA_CHECK_EQUAL ( field(1,1,-1), 5 );
+
+   // check values in top interior cells
+   WALBERLA_CHECK_EQUAL ( field(0,0,1), 2 );
+   WALBERLA_CHECK_EQUAL ( field(0,1,1), 3 );
+   WALBERLA_CHECK_EQUAL ( field(1,0,1), 4 );
+   WALBERLA_CHECK_EQUAL ( field(1,1,1), 5 );
 
+}
 
 int main(int argc, char **argv)
 {
@@ -127,7 +180,8 @@ int main(int argc, char **argv)
    for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) // block loop
       testScalarField( &(*blockIt), scalarFieldId );
 
-
+   for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) // block loop
+      testScalarFieldPullReduction( &(*blockIt), scalarFieldId );
 
    return 0;
 }
diff --git a/tests/field/distributors/DistributionTest.cpp b/tests/field/distributors/DistributionTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4d0e4f644d153c055c9d61c0648568a13df94eb0
--- /dev/null
+++ b/tests/field/distributors/DistributionTest.cpp
@@ -0,0 +1,536 @@
+//======================================================================================================================
+//
+//  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 DistributionTest.cpp
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/all.h"
+
+#include "core/debug/TestSubsystem.h"
+#include "core/mpi/Environment.h"
+#include "core/math/all.h"
+
+#include "field/AddToStorage.h"
+#include "field/FlagField.h"
+#include "field/GhostLayerField.h"
+#include "field/distributors/all.h"
+
+#include <vector>
+
+using namespace walberla;
+
+namespace distribution_tests {
+
+const uint_t FieldGhostLayers( 1 );
+
+typedef walberla::uint8_t                    flag_t;
+typedef FlagField< flag_t >                  FlagField_T;
+
+typedef GhostLayerField< real_t, 1>          ScalarField_T;
+typedef GhostLayerField< Vector3<real_t>, 1> Vec3Field_T;
+typedef GhostLayerField< real_t, 3>          MultiComponentField_T;
+
+
+const FlagUID Domain_Flag ( "domain" );
+const FlagUID Boundary_Flag ( "boundary" );
+
+void initFlagField( FlagField_T * field, IBlock * const /*block*/ )
+{
+   auto domainFlag = field->getOrRegisterFlag( Domain_Flag );
+   WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( field, field->addFlag( x, y, z, domainFlag ); );
+}
+
+void resetScalarField( const shared_ptr<StructuredBlockStorage> & blocks,
+                       const BlockDataID & scalarFieldID )
+{
+   for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+   {
+      auto sField = blockIt->getData<ScalarField_T>( scalarFieldID );
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(sField,
+                                 sField->get(x,y,z) = real_t(0);
+      );
+   }
+}
+
+void resetVectorField( const shared_ptr<StructuredBlockStorage> & blocks,
+                       const BlockDataID & vectorFieldID )
+{
+   for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+   {
+      auto vField = blockIt->getData<Vec3Field_T>( vectorFieldID );
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(vField,
+                                                       vField->get(x,y,z) = Vector3<real_t>(real_t(0));
+      );
+   }
+}
+
+void resetMultiCompField( const shared_ptr<StructuredBlockStorage> & blocks,
+                          const BlockDataID & multiComponentFieldID )
+{
+   for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+   {
+      auto mField = blockIt->getData<MultiComponentField_T>( multiComponentFieldID );
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(mField,
+                                                       mField->get(x,y,z,0) = real_t(0);
+                                                       mField->get(x,y,z,1) = real_t(0);
+                                                       mField->get(x,y,z,2) = real_t(0);
+      );
+   }
+}
+
+
+void setBoundaryFlags( const shared_ptr<StructuredBlockStorage> & blocks,
+                       const BlockDataID & flagFieldID )
+{
+   for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+   {
+      auto flagField = blockIt->getData<FlagField_T>( flagFieldID );
+      auto domainFlag = flagField->getOrRegisterFlag( Domain_Flag );
+      auto boundaryFlag = flagField->getOrRegisterFlag( Boundary_Flag );
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField,
+                                                       if( x == 2)
+                                                       {
+                                                          flagField->removeFlag(x,y,z,domainFlag);
+                                                          flagField->addFlag(x,y,z,boundaryFlag);
+                                                       }
+      );
+   }
+}
+
+void getScalarFieldQuantities( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID,
+                               const BlockDataID & fieldID,
+                               real_t & summarizedValue, real_t & minValue, real_t & maxValue)
+{
+   real_t sum( real_t(0) );
+   real_t min( real_t(0) );
+   real_t max( real_t(0) );
+   for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+   {
+      auto field = blockIt->getData<ScalarField_T>( fieldID );
+      auto flagField = blockIt->getData<FlagField_T>( flagFieldID );
+      auto domainFlag = flagField->getFlag( Domain_Flag );
+
+      CellInterval xyzSizeWithGhostLayers = field->xyzSizeWithGhostLayer();
+      for( auto cellIt = xyzSizeWithGhostLayers.begin(); cellIt != xyzSizeWithGhostLayers.end(); ++cellIt )
+      {
+         if( flagField->isFlagSet(*cellIt,domainFlag))
+         {
+            real_t value = field->get(*cellIt);
+            sum += value;
+            min = std::min(min, value);
+            max = std::max(max, value);
+         }
+      }
+   }
+   WALBERLA_MPI_SECTION()
+   {
+      mpi::allReduceInplace( sum, mpi::SUM );
+      mpi::allReduceInplace( min, mpi::MIN );
+      mpi::allReduceInplace( max, mpi::MAX );
+   }
+
+   summarizedValue = sum;
+   minValue = min;
+   maxValue = max;
+}
+
+void getVectorFieldQuantities( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID,
+                               const BlockDataID & fieldID,
+                               Vector3<real_t> & summarizedValue, Vector3<real_t> & minValue, Vector3<real_t> & maxValue)
+{
+   Vector3<real_t> sum( real_t(0) );
+   Vector3<real_t> min( real_t(0) );
+   Vector3<real_t> max( real_t(0) );
+   for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+   {
+      auto field = blockIt->getData<Vec3Field_T>( fieldID );
+      auto flagField = blockIt->getData<FlagField_T>( flagFieldID );
+      auto domainFlag = flagField->getFlag( Domain_Flag );
+      CellInterval xyzSizeWithGhostLayers = field->xyzSizeWithGhostLayer();
+      for( auto cellIt = xyzSizeWithGhostLayers.begin(); cellIt != xyzSizeWithGhostLayers.end(); ++cellIt )
+      {
+         if( flagField->isFlagSet(*cellIt,domainFlag))
+         {
+            Vector3<real_t> value = field->get(*cellIt);
+            sum += value;
+            for (size_t i = 0; i < 3; ++i) {
+               min[i] = std::min(min[i], value[i]);
+               max[i] = std::max(max[i], value[i]);
+            }
+         }
+      }
+   }
+   WALBERLA_MPI_SECTION()
+   {
+      for( size_t i = 0; i < 3; ++i) {
+         sum[i] = mpi::allReduce( sum[i], mpi::SUM );
+         min[i] = mpi::allReduce( min[i], mpi::MIN );
+         max[i] = mpi::allReduce( max[i], mpi::MAX );
+      }
+   }
+
+   summarizedValue = sum;
+   minValue = min;
+   maxValue = max;
+}
+
+void getMultiCompFieldQuantities( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID,
+                                  const BlockDataID & fieldID,
+                                  std::vector<real_t> & summarizedValue, std::vector<real_t> & minValue, std::vector<real_t> & maxValue )
+{
+   std::vector<real_t> sum( 3 );
+   std::vector<real_t> min( 3 );
+   std::vector<real_t> max( 3 );
+   for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+   {
+      auto field = blockIt->getData<MultiComponentField_T>( fieldID );
+      auto flagField = blockIt->getData<FlagField_T>( flagFieldID );
+      auto domainFlag = flagField->getFlag( Domain_Flag );
+      CellInterval xyzSizeWithGhostLayers = field->xyzSizeWithGhostLayer();
+      for( auto cellIt = xyzSizeWithGhostLayers.begin(); cellIt != xyzSizeWithGhostLayers.end(); ++cellIt )
+      {
+         if( flagField->isFlagSet(*cellIt,domainFlag))
+         {
+            real_t value0 = field->get(*cellIt, 0);
+            sum[0] += value0;
+            min[0] = std::min(min[0], value0);
+            max[0] = std::max(max[0], value0);
+            real_t value1 = field->get(*cellIt, 1);
+            sum[1] += value1;
+            min[1] = std::min(min[1], value1);
+            max[1] = std::max(max[1], value1);
+            real_t value2 = field->get(*cellIt, 2);
+            sum[2] += value2;
+            min[2] = std::min(min[2], value2);
+            max[2] = std::max(max[2], value2);
+         }
+      }
+   }
+   WALBERLA_MPI_SECTION()
+   {
+      mpi::allReduceInplace( sum, mpi::SUM );
+      mpi::allReduceInplace( min, mpi::MIN );
+      mpi::allReduceInplace( max, mpi::MAX );
+   }
+
+   summarizedValue = sum;
+   minValue = min;
+   maxValue = max;
+}
+
+void testNearestNeighborDistributor( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID,
+                                     const BlockDataID & scalarFieldID, const BlockDataID & vectorFieldID, const BlockDataID & multiComponentFieldID )
+{
+   // distributors
+   typedef field::NearestNeighborDistributor<ScalarField_T, FlagField_T>         ScalarDistributor_T;
+   typedef field::NearestNeighborDistributor<Vec3Field_T, FlagField_T>           Vec3Distributor_T;
+   typedef field::NearestNeighborDistributor<MultiComponentField_T, FlagField_T> MultiComponentDistributor_T;
+   BlockDataID scalarDistributorID         = field::addDistributor< ScalarDistributor_T, FlagField_T >( blocks, scalarFieldID, flagFieldID, Domain_Flag );
+   BlockDataID vectorDistributorID         = field::addDistributor< Vec3Distributor_T, FlagField_T >( blocks, vectorFieldID, flagFieldID, Domain_Flag );
+   BlockDataID multiComponentDistributorID = field::addDistributor< MultiComponentDistributor_T, FlagField_T >( blocks, multiComponentFieldID, flagFieldID, Domain_Flag );
+
+   // check scalar distribution
+   {
+      Vector3<real_t> distributionPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+      real_t distributionValue(real_t(100));
+      auto containingBlockID = blocks->getBlock(distributionPoint);
+      if( containingBlockID != NULL )
+      {
+         auto distPtr = containingBlockID->getData<ScalarDistributor_T>(scalarDistributorID);
+         distPtr->distribute(distributionPoint, &distributionValue);
+      }
+      real_t sum(real_t(0)), min(real_t(0)), max(real_t(0));
+      getScalarFieldQuantities(blocks, flagFieldID, scalarFieldID, sum, min, max);
+      WALBERLA_CHECK_FLOAT_EQUAL(sum, distributionValue, "NearestNeighborDistributor: Sum of scalar distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(min, real_t(0), "NearestNeighborDistributor: Min of scalar distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(max, distributionValue, "NearestNeighborDistributor: Max of scalar distribution failed!" );
+
+      resetScalarField(blocks, scalarFieldID);
+   }
+
+   // check vector distribution
+   {
+      Vector3<real_t> distributionPoint(real_t(5.4),real_t(2.1),real_t(3.2));
+      Vector3<real_t> distributionValue(real_t(100), real_t(-10), real_t(1));
+      auto containingBlockID = blocks->getBlock(distributionPoint);
+      if( containingBlockID != NULL )
+      {
+         auto distPtr = containingBlockID->getData<Vec3Distributor_T>(vectorDistributorID);
+         distPtr->distribute(distributionPoint, &distributionValue);
+      }
+      Vector3<real_t> sum(real_t(0)), min(real_t(0)), max(real_t(0));
+      getVectorFieldQuantities(blocks, flagFieldID, vectorFieldID, sum, min, max);
+      WALBERLA_CHECK_FLOAT_EQUAL(sum[0], distributionValue[0], "NearestNeighborDistributor: Sum of vec[0] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(sum[1], distributionValue[1], "NearestNeighborDistributor: Sum of vec[1] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(sum[2], distributionValue[2], "NearestNeighborDistributor: Sum of vec[2] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(min[0], real_t(0), "NearestNeighborDistributor: Min of vec[0] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(min[1], distributionValue[1], "NearestNeighborDistributor: Min of vec[1] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(min[2], real_t(0), "NearestNeighborDistributor: Min of vec[2] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(max[0], distributionValue[0], "NearestNeighborDistributor: Max of vec[0] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(max[1], real_t(0), "NearestNeighborDistributor: Max of vec[1] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(max[2], distributionValue[2], "NearestNeighborDistributor: Max of vec[2] distribution failed!" );
+
+      resetVectorField(blocks, vectorFieldID);
+   }
+
+   // check multi component distribution
+   {
+      Vector3<real_t> distributionPoint(real_t(4.4),real_t(2.1),real_t(3.2));
+      std::vector<real_t> distributionValue(3);
+      distributionValue[0] = real_t(100);
+      distributionValue[1] = real_t(-10);
+      distributionValue[2] = real_t(1);
+      auto containingBlockID = blocks->getBlock(distributionPoint);
+      if( containingBlockID != NULL )
+      {
+         auto distPtr = containingBlockID->getData<MultiComponentDistributor_T>(multiComponentDistributorID);
+         distPtr->distribute(distributionPoint, distributionValue.begin());
+      }
+      std::vector<real_t> sum(3), min(3), max(3);
+      getMultiCompFieldQuantities(blocks, flagFieldID, multiComponentFieldID, sum, min, max);
+      WALBERLA_CHECK_FLOAT_EQUAL(sum[0], distributionValue[0], "NearestNeighborDistributor: Sum of Multi Component[0] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(sum[1], distributionValue[1], "NearestNeighborDistributor: Sum of Multi Component[1] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(sum[2], distributionValue[2], "NearestNeighborDistributor: Sum of Multi Component[2] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(min[0], real_t(0), "NearestNeighborDistributor: Min of Multi Component[0] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(min[1], distributionValue[1], "NearestNeighborDistributor: Min of Multi Component[1] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(min[2], real_t(0), "NearestNeighborDistributor: Min of Multi Component[2] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(max[0], distributionValue[0], "NearestNeighborDistributor: Max of Multi Component[0] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(max[1], real_t(0), "NearestNeighborDistributor: Max of Multi Component[1] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(max[2], distributionValue[2], "NearestNeighborDistributor: Max of Multi Component[2] distribution failed!" );
+
+      resetMultiCompField(blocks, multiComponentFieldID);
+   }
+}
+
+
+void testKernelDistributor( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID,
+                            const BlockDataID & scalarFieldID, const BlockDataID & vectorFieldID, const BlockDataID & multiComponentFieldID )
+{
+   // distributors
+   typedef field::KernelDistributor<ScalarField_T, FlagField_T>         ScalarDistributor_T;
+   typedef field::KernelDistributor<Vec3Field_T, FlagField_T>           Vec3Distributor_T;
+   typedef field::KernelDistributor<MultiComponentField_T, FlagField_T> MultiComponentDistributor_T;
+   BlockDataID scalarDistributorID         = field::addDistributor< ScalarDistributor_T, FlagField_T >( blocks, scalarFieldID, flagFieldID, Domain_Flag );
+   BlockDataID vectorDistributorID         = field::addDistributor< Vec3Distributor_T, FlagField_T >( blocks, vectorFieldID, flagFieldID, Domain_Flag );
+   BlockDataID multiComponentDistributorID = field::addDistributor< MultiComponentDistributor_T, FlagField_T >( blocks, multiComponentFieldID, flagFieldID, Domain_Flag );
+
+   // check scalar distribution
+   {
+      Vector3<real_t> distributionPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+      real_t distributionValue(real_t(100));
+      auto containingBlockID = blocks->getBlock(distributionPoint);
+      if( containingBlockID != NULL )
+      {
+         auto distPtr = containingBlockID->getData<ScalarDistributor_T>(scalarDistributorID);
+         distPtr->distribute(distributionPoint, &distributionValue);
+      }
+      real_t sum(real_t(0)), min(real_t(0)), max(real_t(0));
+      getScalarFieldQuantities(blocks, flagFieldID, scalarFieldID, sum, min, max);
+      WALBERLA_CHECK_FLOAT_EQUAL(sum, distributionValue, "KernelDistributor: Sum of scalar distribution failed!" );
+      WALBERLA_CHECK(min >= real_t(0), "KernelDistributor: Min of scalar distribution failed!" );
+      WALBERLA_CHECK(max <= distributionValue, "KernelDistributor: Max of scalar distribution failed!" );
+
+      resetScalarField(blocks, scalarFieldID);
+   }
+
+   // check vector distribution
+   {
+      Vector3<real_t> distributionPoint(real_t(5.4),real_t(2.1),real_t(3.2));
+      Vector3<real_t> distributionValue(real_t(100), real_t(-10), real_t(1));
+      auto containingBlockID = blocks->getBlock(distributionPoint);
+      if( containingBlockID != NULL )
+      {
+         auto distPtr = containingBlockID->getData<Vec3Distributor_T>(vectorDistributorID);
+         distPtr->distribute(distributionPoint, &distributionValue);
+      }
+      Vector3<real_t> sum(real_t(0)), min(real_t(0)), max(real_t(0));
+      getVectorFieldQuantities(blocks, flagFieldID, vectorFieldID, sum, min, max);
+      WALBERLA_CHECK_FLOAT_EQUAL(sum[0], distributionValue[0], "KernelDistributor: Sum of vec[0] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(sum[1], distributionValue[1], "KernelDistributor: Sum of vec[1] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(sum[2], distributionValue[2], "KernelDistributor: Sum of vec[2] distribution failed!" );
+      WALBERLA_CHECK(min[0] >= real_t(0), "KernelDistributor: Min of vec[0] distribution failed!" );
+      WALBERLA_CHECK(min[1] >= distributionValue[1], "KernelDistributor: Min of vec[1] distribution failed!" );
+      WALBERLA_CHECK(min[2] >= real_t(0), "KernelDistributor: Min of vec[2] distribution failed!" );
+      WALBERLA_CHECK(max[0] <= distributionValue[0], "KernelDistributor: Max of vec[0] distribution failed!" );
+      WALBERLA_CHECK(max[1] <= real_t(0), "KernelDistributor: Max of vec[1] distribution failed!" );
+      WALBERLA_CHECK(max[2] <= distributionValue[2], "KernelDistributor: Max of vec[2] distribution failed!" );
+
+      resetVectorField(blocks, vectorFieldID);
+   }
+
+   // check multi component distribution
+   {
+      Vector3<real_t> distributionPoint(real_t(4.4),real_t(2.1),real_t(3.2));
+      std::vector<real_t> distributionValue(3);
+      distributionValue[0] = real_t(100);
+      distributionValue[1] = real_t(-10);
+      distributionValue[2] = real_t(1);
+      auto containingBlockID = blocks->getBlock(distributionPoint);
+      if( containingBlockID != NULL )
+      {
+         auto distPtr = containingBlockID->getData<MultiComponentDistributor_T>(multiComponentDistributorID);
+         distPtr->distribute(distributionPoint, distributionValue.begin());
+      }
+      std::vector<real_t> sum(3), min(3), max(3);
+      getMultiCompFieldQuantities(blocks, flagFieldID, multiComponentFieldID, sum, min, max);
+      WALBERLA_CHECK_FLOAT_EQUAL(sum[0], distributionValue[0], "KernelDistributor: Sum of Multi Component[0] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(sum[1], distributionValue[1], "KernelDistributor: Sum of Multi Component[1] distribution failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(sum[2], distributionValue[2], "KernelDistributor: Sum of Multi Component[2] distribution failed!" );
+      WALBERLA_CHECK(min[0] >= real_t(0), "KernelDistributor: Min of Multi Component[0] distribution failed!" );
+      WALBERLA_CHECK(min[1] >= distributionValue[1], "KernelDistributor: Min of Multi Component[1] distribution failed!" );
+      WALBERLA_CHECK(min[2] >= real_t(0), "KernelDistributor: Min of Multi Component[2] distribution failed!" );
+      WALBERLA_CHECK(max[0] <= distributionValue[0], "KernelDistributor: Max of Multi Component[0] distribution failed!" );
+      WALBERLA_CHECK(max[1] <= real_t(0), "KernelDistributor: Max of Multi Component[1] distribution failed!" );
+      WALBERLA_CHECK(max[2] <= distributionValue[2], "KernelDistributor: Max of Multi Component[2] distribution failed!" );
+
+      resetMultiCompField(blocks, multiComponentFieldID);
+   }
+}
+
+
+void testNearestNeighborDistributorAtBoundary( const shared_ptr<StructuredBlockStorage> & blocks,
+                                               const BlockDataID & flagFieldID, const BlockDataID & scalarFieldID )
+{
+   // distributor
+   typedef field::NearestNeighborDistributor<ScalarField_T, FlagField_T> ScalarDistributor_T;
+   BlockDataID scalarDistributorID = field::addDistributor<ScalarDistributor_T, FlagField_T>(blocks, scalarFieldID, flagFieldID, Domain_Flag);
+
+   // check scalar interpolation close to boundary
+   {
+      Vector3<real_t> distributionPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+      real_t distributionValue(real_t(100));
+      auto containingBlockID = blocks->getBlock(distributionPoint);
+      if (containingBlockID != NULL) {
+         auto distPtr = containingBlockID->getData<ScalarDistributor_T>(scalarDistributorID);
+         distPtr->distribute(distributionPoint, &distributionValue);
+      }
+      real_t sum(real_t(0)), min(real_t(0)), max(real_t(0));
+      getScalarFieldQuantities(blocks, flagFieldID, scalarFieldID, sum, min, max);
+      WALBERLA_CHECK_FLOAT_EQUAL(sum, distributionValue, "NearestNeighborDistributor: Sum of scalar distribution near boundary failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(min, real_t(0), "NearestNeighborDistributor: Min of scalar distribution near boundary failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(max, distributionValue, "NearestNeighborDistributor: Max of scalar distribution near boundary failed!" );
+
+      resetScalarField(blocks, scalarFieldID);
+   }
+
+   // check scalar interpolation inside boundary
+   {
+      Vector3<real_t> distributionPoint(real_t(2.7), real_t(2.1), real_t(1.1));
+      real_t distributionValue(real_t(100));
+      auto containingBlockID = blocks->getBlock(distributionPoint);
+      if (containingBlockID != NULL) {
+         auto distPtr = containingBlockID->getData<ScalarDistributor_T>(scalarDistributorID);
+         distPtr->distribute(distributionPoint, &distributionValue);
+      }
+      real_t sum(real_t(0)), min(real_t(0)), max(real_t(0));
+      getScalarFieldQuantities(blocks, flagFieldID, scalarFieldID, sum, min, max);
+      WALBERLA_CHECK_FLOAT_EQUAL(sum, distributionValue, "NearestNeighborDistributor: Sum of scalar distribution inside boundary failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(min, real_t(0), "NearestNeighborDistributor: Min of scalar distribution inside boundary failed!" );
+      WALBERLA_CHECK_FLOAT_EQUAL(max, distributionValue, "NearestNeighborDistributor: Max of scalar distribution inside boundary failed!" );
+
+      resetScalarField(blocks, scalarFieldID);
+   }
+}
+
+void testKernelDistributorAtBoundary( const shared_ptr<StructuredBlockStorage> & blocks,
+                                      const BlockDataID & flagFieldID, const BlockDataID & scalarFieldID )
+{
+   // distributor
+   typedef field::KernelDistributor<ScalarField_T, FlagField_T> ScalarDistributor_T;
+   BlockDataID scalarDistributorID = field::addDistributor<ScalarDistributor_T, FlagField_T>(blocks, scalarFieldID, flagFieldID, Domain_Flag);
+
+   // check scalar interpolation close to boundary
+   {
+      Vector3<real_t> distributionPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+      real_t distributionValue(real_t(100));
+      auto containingBlockID = blocks->getBlock(distributionPoint);
+      if (containingBlockID != NULL) {
+         auto distPtr = containingBlockID->getData<ScalarDistributor_T>(scalarDistributorID);
+         distPtr->distribute(distributionPoint, &distributionValue);
+      }
+      real_t sum(real_t(0)), min(real_t(0)), max(real_t(0));
+      getScalarFieldQuantities(blocks, flagFieldID, scalarFieldID, sum, min, max);
+      WALBERLA_CHECK_FLOAT_EQUAL(sum, distributionValue, "KernelDistributor: Sum of scalar distribution near boundary failed!" );
+      WALBERLA_CHECK(min >= real_t(0), "KernelDistributor: Min of scalar distribution near boundary failed!" );
+      WALBERLA_CHECK(max <= distributionValue, "KernelDistributor: Max of scalar distribution near boundary failed!" );
+
+      resetScalarField(blocks, scalarFieldID);
+   }
+
+   // check scalar interpolation inside boundary
+   {
+      Vector3<real_t> distributionPoint(real_t(2.7), real_t(2.1), real_t(1.1));
+      real_t distributionValue(real_t(100));
+      auto containingBlockID = blocks->getBlock(distributionPoint);
+      if (containingBlockID != NULL) {
+         auto distPtr = containingBlockID->getData<ScalarDistributor_T>(scalarDistributorID);
+         distPtr->distribute(distributionPoint, &distributionValue);
+      }
+      real_t sum(real_t(0)), min(real_t(0)), max(real_t(0));
+      getScalarFieldQuantities(blocks, flagFieldID, scalarFieldID, sum, min, max);
+      WALBERLA_CHECK_FLOAT_EQUAL(sum, distributionValue, "KernelDistributor: Sum of scalar distribution inside boundary failed!" );
+      WALBERLA_CHECK(min >= real_t(0), "KernelDistributor: Min of scalar distribution inside boundary failed!" );
+      WALBERLA_CHECK(max <= distributionValue, "KernelDistributor: Max of scalar distribution inside boundary failed!" );
+
+      resetScalarField(blocks, scalarFieldID);
+   }
+}
+
+
+int main(int argc, char **argv) {
+
+   mpi::Environment mpiEnv(argc, argv);
+   debug::enterTestMode();
+
+   const uint_t numberOfBlocksInDirection = 2;
+   const uint_t numberOfCellsPerBlockInDirection = 4;
+   const real_t dx = real_t(1);
+
+   // block storage
+   auto blocks = blockforest::createUniformBlockGrid( numberOfBlocksInDirection, numberOfBlocksInDirection, numberOfBlocksInDirection,
+                                                      numberOfCellsPerBlockInDirection, numberOfCellsPerBlockInDirection, numberOfCellsPerBlockInDirection,
+                                                      dx, 0, false, false,
+                                                      false, false, false,
+                                                      false );
+   // flag field
+   BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field", FieldGhostLayers, false, initFlagField );
+
+   // data fields
+   BlockDataID scalarFieldID         = field::addToStorage< ScalarField_T >( blocks, "scalar field", real_t(0), field::zyxf, FieldGhostLayers );
+   BlockDataID vectorFieldID         = field::addToStorage< Vec3Field_T >( blocks, "vec3 field", Vector3<real_t>(real_t(0)), field::zyxf, FieldGhostLayers );
+   BlockDataID multiComponentFieldID = field::addToStorage< MultiComponentField_T >( blocks, "multi component field", real_t(0), field::zyxf, FieldGhostLayers );
+
+   // test all distributors with domain flags everywhere, i.e. without special boundary treatment necessary
+   testNearestNeighborDistributor(blocks, flagFieldID, scalarFieldID, vectorFieldID, multiComponentFieldID);
+   testKernelDistributor(blocks, flagFieldID, scalarFieldID, vectorFieldID, multiComponentFieldID);
+
+   // set some boundary flags in flag field and invalidate the corresponding scalar field values
+   setBoundaryFlags(blocks, flagFieldID );
+
+   // test all distributors' behavior close to boundary cells
+   testNearestNeighborDistributorAtBoundary(blocks, flagFieldID, scalarFieldID);
+   testKernelDistributorAtBoundary(blocks, flagFieldID, scalarFieldID);
+
+   return 0;
+}
+
+} // namespace field_distribution_tests
+
+int main( int argc, char **argv ){
+   distribution_tests::main(argc, argv);
+}
\ No newline at end of file
diff --git a/tests/field/interpolators/FieldInterpolationTest.cpp b/tests/field/interpolators/FieldInterpolationTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..32adf2898cf57a2ccbc988d01d0f9de02c8ebfd9
--- /dev/null
+++ b/tests/field/interpolators/FieldInterpolationTest.cpp
@@ -0,0 +1,430 @@
+//======================================================================================================================
+//
+//  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 FieldInterpolationTest.cpp
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/all.h"
+
+#include "core/debug/TestSubsystem.h"
+#include "core/mpi/Environment.h"
+#include "core/math/all.h"
+
+#include "field/AddToStorage.h"
+#include "field/FlagField.h"
+#include "field/GhostLayerField.h"
+#include "field/interpolators/all.h"
+
+#include <vector>
+
+using namespace walberla;
+
+namespace field_interpolation_tests {
+
+const uint_t FieldGhostLayers( 1 );
+
+typedef walberla::uint8_t                    flag_t;
+typedef FlagField< flag_t >                  FlagField_T;
+
+typedef GhostLayerField< real_t, 1>          ScalarField_T;
+typedef GhostLayerField< Vector3<real_t>, 1> Vec3Field_T;
+typedef GhostLayerField< real_t, 3>          MultiComponentField_T;
+
+
+const FlagUID Domain_Flag ( "domain" );
+const FlagUID Boundary_Flag ( "boundary" );
+
+void initFlagField( FlagField_T * field, IBlock * const /*block*/ )
+{
+   auto domainFlag = field->getOrRegisterFlag( Domain_Flag );
+   WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( field, field->addFlag( x, y, z, domainFlag ); );
+}
+
+void initScalarField( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & scalarFieldID )
+{
+   for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+   {
+      auto field = blockIt->getData<ScalarField_T>( scalarFieldID );
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(field,
+                                 const Vector3<real_t> p = blocks->getBlockLocalCellCenter(*blockIt, Cell(x,y,z)) - Vector3<real_t>(real_t(0.5));
+                                 //field->get(x,y,z) = real_t(2);
+                                 field->get(x,y,z) = p[0];
+      );
+   }
+}
+
+void initVectorField( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & vectorFieldID )
+{
+   for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+   {
+      auto field = blockIt->getData<Vec3Field_T>( vectorFieldID );
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(field,
+                                 const Vector3<real_t> p = blocks->getBlockLocalCellCenter(*blockIt, Cell(x,y,z)) - Vector3<real_t>(real_t(0.5));
+                                 field->get(x,y,z) = Vector3<real_t>(p[0], real_t(2), p[1]);
+      );
+   }
+}
+
+void initMultiComponentField( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & multiComponentFieldID )
+{
+   for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+   {
+      auto field = blockIt->getData<MultiComponentField_T>( multiComponentFieldID );
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(field,
+                                 const Vector3<real_t> p = blocks->getBlockLocalCellCenter(*blockIt, Cell(x,y,z)) - Vector3<real_t>(real_t(0.5));
+                                 field->get(x,y,z,0) = p[0];
+                                 field->get(x,y,z,1) = real_t(2);
+                                 field->get(x,y,z,2) = p[1];
+      );
+   }
+}
+
+void setBoundaryFlags( const shared_ptr<StructuredBlockStorage> & blocks,
+                       const BlockDataID & flagFieldID, const BlockDataID & scalarFieldID )
+{
+   for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+   {
+      auto flagField = blockIt->getData<FlagField_T>( flagFieldID );
+      auto valueField = blockIt->getData<ScalarField_T>( scalarFieldID );
+      auto domainFlag = flagField->getOrRegisterFlag( Domain_Flag );
+      auto boundaryFlag = flagField->getOrRegisterFlag( Boundary_Flag );
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField,
+                                                       if( x == 2)
+                                                       {
+                                                          flagField->removeFlag(x,y,z,domainFlag);
+                                                          flagField->addFlag(x,y,z,boundaryFlag);
+                                                          valueField->get(x,y,z) = std::numeric_limits<real_t>::max();
+                                                       }
+      );
+   }
+}
+
+void testNearestNeighborFieldInterpolator( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID,
+                                           const BlockDataID & scalarFieldID, const BlockDataID & vectorFieldID, const BlockDataID & multiComponentFieldID )
+{
+   // field interpolators
+   typedef field::NearestNeighborFieldInterpolator<ScalarField_T, FlagField_T>         ScalarFieldInterpolator_T;
+   typedef field::NearestNeighborFieldInterpolator<Vec3Field_T, FlagField_T>           Vec3FieldInterpolator_T;
+   typedef field::NearestNeighborFieldInterpolator<MultiComponentField_T, FlagField_T> MultiComponentFieldInterpolator_T;
+   BlockDataID scalarFieldInterpolatorID         = field::addFieldInterpolator< ScalarFieldInterpolator_T, FlagField_T >( blocks, scalarFieldID, flagFieldID, Domain_Flag );
+   BlockDataID vectorFieldInterpolatorID         = field::addFieldInterpolator< Vec3FieldInterpolator_T, FlagField_T >( blocks, vectorFieldID, flagFieldID, Domain_Flag );
+   BlockDataID multiComponentFieldInterpolatorID = field::addFieldInterpolator< MultiComponentFieldInterpolator_T, FlagField_T >( blocks, multiComponentFieldID, flagFieldID, Domain_Flag );
+
+   // check scalar interpolation
+   {
+      Vector3<real_t> interpolationPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+      auto containingBlockID = blocks->getBlock(interpolationPoint);
+      if( containingBlockID != NULL )
+      {
+         real_t interpolatedValue(real_t(0));
+         auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+         interPtr->get(interpolationPoint, &interpolatedValue);
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(1), "NearestNeighborFieldInterpolator: Scalar interpolation failed");
+      }
+   }
+
+   // check vector interpolation
+   {
+      Vector3<real_t> interpolationPoint(real_t(5.4),real_t(2.1),real_t(3.2));
+      auto containingBlockID = blocks->getBlock( interpolationPoint );
+      if( containingBlockID != NULL ) {
+         Vector3<real_t> interpolatedValue(real_t(0));
+         auto interPtr = containingBlockID->getData<Vec3FieldInterpolator_T>(vectorFieldInterpolatorID);
+         interPtr->get(interpolationPoint, &interpolatedValue);
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[0], real_t(5), "NearestNeighborFieldInterpolator: Vec3[0] interpolation failed");
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[1], real_t(2), "NearestNeighborFieldInterpolator: Vec3[1] interpolation failed");
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[2], real_t(2), "NearestNeighborFieldInterpolator: Vec3[2] interpolation failed");
+      }
+   }
+
+   // check multi component interpolation
+   {
+      Vector3<real_t> interpolationPoint(real_t(4.4),real_t(2.1),real_t(3.2));
+      auto containingBlockID = blocks->getBlock( interpolationPoint );
+      if( containingBlockID != NULL ) {
+         std::vector<real_t> interpolatedValue(3, real_t(0));
+         auto interPtr = containingBlockID->getData<MultiComponentFieldInterpolator_T>(multiComponentFieldInterpolatorID);
+         interPtr->get(interpolationPoint, interpolatedValue.begin());
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[0], real_t(4), "NearestNeighborFieldInterpolator: Multi Component[0] interpolation failed");
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[1], real_t(2), "NearestNeighborFieldInterpolator: Multi Component[1] interpolation failed");
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[2], real_t(2), "NearestNeighborFieldInterpolator: Multi Component[2] interpolation failed");
+      }
+   }
+}
+
+void testTrilinearFieldInterpolator( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID,
+                                     const BlockDataID & scalarFieldID, const BlockDataID & vectorFieldID, const BlockDataID & multiComponentFieldID )
+{
+   // field interpolators
+   typedef field::TrilinearFieldInterpolator<ScalarField_T, FlagField_T>         ScalarFieldInterpolator_T;
+   typedef field::TrilinearFieldInterpolator<Vec3Field_T, FlagField_T>           Vec3FieldInterpolator_T;
+   typedef field::TrilinearFieldInterpolator<MultiComponentField_T, FlagField_T> MultiComponentFieldInterpolator_T;
+   BlockDataID scalarFieldInterpolatorID         = field::addFieldInterpolator< ScalarFieldInterpolator_T, FlagField_T >( blocks, scalarFieldID, flagFieldID, Domain_Flag );
+   BlockDataID vectorFieldInterpolatorID         = field::addFieldInterpolator< Vec3FieldInterpolator_T, FlagField_T >( blocks, vectorFieldID, flagFieldID, Domain_Flag );
+   BlockDataID multiComponentFieldInterpolatorID = field::addFieldInterpolator< MultiComponentFieldInterpolator_T, FlagField_T >( blocks, multiComponentFieldID, flagFieldID, Domain_Flag );
+
+   // check scalar interpolation
+   {
+      Vector3<real_t> interpolationPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+      auto containingBlockID = blocks->getBlock(interpolationPoint);
+      if( containingBlockID != NULL )
+      {
+         real_t interpolatedValue(real_t(0));
+         auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+         interPtr->get(interpolationPoint, &interpolatedValue);
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(1.4), "TrilinearFieldInterpolator: Scalar interpolation failed");
+      }
+   }
+
+   // check vector interpolation
+   {
+      Vector3<real_t> interpolationPoint(real_t(5.4),real_t(2.1),real_t(3.2));
+      auto containingBlockID = blocks->getBlock( interpolationPoint );
+      if( containingBlockID != NULL ) {
+         Vector3<real_t> interpolatedValue(real_t(0));
+         auto interPtr = containingBlockID->getData<Vec3FieldInterpolator_T>(vectorFieldInterpolatorID);
+         interPtr->get(interpolationPoint, &interpolatedValue);
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[0], real_t(4.9), "TrilinearFieldInterpolator: Vec3[0] interpolation failed");
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[1], real_t(2.0), "TrilinearFieldInterpolator: Vec3[1] interpolation failed");
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[2], real_t(1.6), "TrilinearFieldInterpolator: Vec3[2] interpolation failed");
+      }
+   }
+
+   // check multi component interpolation
+   {
+      Vector3<real_t> interpolationPoint(real_t(4.4),real_t(2.1),real_t(3.2));
+      auto containingBlockID = blocks->getBlock( interpolationPoint );
+      if( containingBlockID != NULL ) {
+         std::vector<real_t> interpolatedValue(3, real_t(0));
+         auto interPtr = containingBlockID->getData<MultiComponentFieldInterpolator_T>(multiComponentFieldInterpolatorID);
+         interPtr->get(interpolationPoint, interpolatedValue.begin());
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[0], real_t(3.9), "TrilinearFieldInterpolator: Multi Component[0] interpolation failed");
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[1], real_t(2.0), "TrilinearFieldInterpolator: Multi Component[1] interpolation failed");
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[2], real_t(1.6), "TrilinearFieldInterpolator: Multi Component[2] interpolation failed");
+      }
+   }
+}
+
+void testKernelFieldInterpolator( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID,
+                                  const BlockDataID & scalarFieldID, const BlockDataID & vectorFieldID, const BlockDataID & multiComponentFieldID )
+{
+   // field interpolators
+   typedef field::KernelFieldInterpolator<ScalarField_T, FlagField_T>         ScalarFieldInterpolator_T;
+   typedef field::KernelFieldInterpolator<Vec3Field_T, FlagField_T>           Vec3FieldInterpolator_T;
+   typedef field::KernelFieldInterpolator<MultiComponentField_T, FlagField_T> MultiComponentFieldInterpolator_T;
+   BlockDataID scalarFieldInterpolatorID         = field::addFieldInterpolator< ScalarFieldInterpolator_T, FlagField_T >( blocks, scalarFieldID, flagFieldID, Domain_Flag );
+   BlockDataID vectorFieldInterpolatorID         = field::addFieldInterpolator< Vec3FieldInterpolator_T, FlagField_T >( blocks, vectorFieldID, flagFieldID, Domain_Flag );
+   BlockDataID multiComponentFieldInterpolatorID = field::addFieldInterpolator< MultiComponentFieldInterpolator_T, FlagField_T >( blocks, multiComponentFieldID, flagFieldID, Domain_Flag );
+
+   // check scalar interpolation
+   {
+      Vector3<real_t> interpolationPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+      auto containingBlockID = blocks->getBlock(interpolationPoint);
+      if( containingBlockID != NULL )
+      {
+         real_t interpolatedValue(real_t(0));
+         auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+         interPtr->get(interpolationPoint, &interpolatedValue);
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(1.4), "KernelFieldInterpolator: Scalar interpolation failed");
+      }
+   }
+
+   // check vector interpolation
+   {
+      Vector3<real_t> interpolationPoint(real_t(5.4),real_t(2.1),real_t(3.2));
+      auto containingBlockID = blocks->getBlock( interpolationPoint );
+      if( containingBlockID != NULL ) {
+         Vector3<real_t> interpolatedValue(real_t(0));
+         auto interPtr = containingBlockID->getData<Vec3FieldInterpolator_T>(vectorFieldInterpolatorID);
+         interPtr->get(interpolationPoint, &interpolatedValue);
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[0], real_t(4.9), "KernelFieldInterpolator: Vec3[0] interpolation failed");
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[1], real_t(2.0), "KernelFieldInterpolator: Vec3[1] interpolation failed");
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[2], real_t(1.6), "KernelFieldInterpolator: Vec3[2] interpolation failed");
+      }
+   }
+
+   // check multi component interpolation
+   {
+      Vector3<real_t> interpolationPoint(real_t(4.4),real_t(2.1),real_t(3.2));
+      auto containingBlockID = blocks->getBlock( interpolationPoint );
+      if( containingBlockID != NULL ) {
+         std::vector<real_t> interpolatedValue(3, real_t(0));
+         auto interPtr = containingBlockID->getData<MultiComponentFieldInterpolator_T>(multiComponentFieldInterpolatorID);
+         interPtr->get(interpolationPoint, interpolatedValue.begin());
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[0], real_t(3.9), "KernelFieldInterpolator: Multi Component[0] interpolation failed");
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[1], real_t(2.0), "KernelFieldInterpolator: Multi Component[1] interpolation failed");
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[2], real_t(1.6), "KernelFieldInterpolator: Multi Component[2] interpolation failed");
+      }
+   }
+}
+
+
+void testNearestNeighborFieldInterpolatorAtBoundary( const shared_ptr<StructuredBlockStorage> & blocks,
+                                                     const BlockDataID & flagFieldID, const BlockDataID & scalarFieldID ) {
+   // field interpolators
+   typedef field::NearestNeighborFieldInterpolator<ScalarField_T, FlagField_T> ScalarFieldInterpolator_T;
+   BlockDataID scalarFieldInterpolatorID = field::addFieldInterpolator<ScalarFieldInterpolator_T, FlagField_T>(blocks, scalarFieldID, flagFieldID, Domain_Flag);
+
+   // check scalar interpolation close to boundary
+   {
+      Vector3<real_t> interpolationPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+      auto containingBlockID = blocks->getBlock(interpolationPoint);
+      if (containingBlockID != NULL) {
+         real_t interpolatedValue(real_t(0));
+         auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+         interPtr->get(interpolationPoint, &interpolatedValue);
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(1),
+                                    "NearestNeighborFieldInterpolator: Scalar interpolation near boundary failed");
+      }
+   }
+
+   // check scalar interpolation inside boundary
+   {
+      Vector3<real_t> interpolationPoint(real_t(2.7), real_t(2.1), real_t(1.1));
+      auto containingBlockID = blocks->getBlock(interpolationPoint);
+      if (containingBlockID != NULL) {
+         real_t interpolatedValue(real_t(0));
+         auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+         interPtr->get(interpolationPoint, &interpolatedValue);
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(3),
+                                    "NearestNeighborFieldInterpolator: Scalar interpolation inside boundary failed");
+      }
+   }
+}
+
+void testTrilinearFieldInterpolatorAtBoundary( const shared_ptr<StructuredBlockStorage> & blocks,
+                                               const BlockDataID & flagFieldID, const BlockDataID & scalarFieldID ) {
+   // field interpolators
+   typedef field::TrilinearFieldInterpolator<ScalarField_T, FlagField_T> ScalarFieldInterpolator_T;
+   BlockDataID scalarFieldInterpolatorID = field::addFieldInterpolator<ScalarFieldInterpolator_T, FlagField_T>(blocks, scalarFieldID, flagFieldID, Domain_Flag);
+
+   // check scalar interpolation close to boundary
+   {
+      Vector3<real_t> interpolationPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+      auto containingBlockID = blocks->getBlock(interpolationPoint);
+      if (containingBlockID != NULL) {
+         real_t interpolatedValue(real_t(0));
+         auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+         interPtr->get(interpolationPoint, &interpolatedValue);
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(1),
+                                    "TrilinearFieldInterpolator: Scalar interpolation near boundary failed");
+      }
+   }
+
+   // check scalar interpolation inside boundary
+   {
+      Vector3<real_t> interpolationPoint(real_t(2.7), real_t(2.1), real_t(1.1));
+      auto containingBlockID = blocks->getBlock(interpolationPoint);
+      if (containingBlockID != NULL) {
+         real_t interpolatedValue(real_t(0));
+         auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+         interPtr->get(interpolationPoint, &interpolatedValue);
+         WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(3),
+                                    "TrilinearFieldInterpolator: Scalar interpolation inside boundary failed");
+      }
+   }
+}
+
+void testKernelFieldInterpolatorAtBoundary( const shared_ptr<StructuredBlockStorage> & blocks,
+                                            const BlockDataID & flagFieldID, const BlockDataID & scalarFieldID ) {
+   // field interpolators
+   typedef field::KernelFieldInterpolator<ScalarField_T, FlagField_T> ScalarFieldInterpolator_T;
+   BlockDataID scalarFieldInterpolatorID = field::addFieldInterpolator<ScalarFieldInterpolator_T, FlagField_T>(blocks, scalarFieldID, flagFieldID, Domain_Flag);
+
+   // check scalar interpolation close to boundary
+   {
+      Vector3<real_t> interpolationPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+      auto containingBlockID = blocks->getBlock(interpolationPoint);
+      if (containingBlockID != NULL) {
+         real_t interpolatedValue(real_t(0));
+         auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+         interPtr->get(interpolationPoint, &interpolatedValue);
+         // kernel interpolation can not extrapolate values from the available ones (see comments in KernelFieldInterpolator.h)
+         // it will thus yield a value between the available ones, which are 0 and 1
+         WALBERLA_CHECK(interpolatedValue < real_t(1),
+                        "KernelFieldInterpolator: Scalar interpolation near boundary failed");
+         WALBERLA_CHECK(interpolatedValue > real_t(0),
+                        "KernelFieldInterpolator: Scalar interpolation near boundary failed");
+      }
+   }
+
+   // check scalar interpolation inside boundary
+   {
+      Vector3<real_t> interpolationPoint(real_t(2.7), real_t(2.1), real_t(1.1));
+      auto containingBlockID = blocks->getBlock(interpolationPoint);
+      if (containingBlockID != NULL) {
+         real_t interpolatedValue(real_t(0));
+         auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+         interPtr->get(interpolationPoint, &interpolatedValue);
+         // values has to be between the available ones, i.e. 1 and 3
+         WALBERLA_CHECK(interpolatedValue > real_t(1),
+                        "KernelFieldInterpolator: Scalar interpolation inside boundary failed");
+         WALBERLA_CHECK(interpolatedValue < real_t(3),
+                        "KernelFieldInterpolator: Scalar interpolation inside boundary failed");
+      }
+   }
+}
+
+int main(int argc, char **argv) {
+
+   mpi::Environment mpiEnv(argc, argv);
+   debug::enterTestMode();
+
+   const uint_t numberOfBlocksInDirection = 2;
+   const uint_t numberOfCellsPerBlockInDirection = 4;
+   const real_t dx = real_t(1);
+
+   // block storage
+   auto blocks = blockforest::createUniformBlockGrid( numberOfBlocksInDirection, numberOfBlocksInDirection, numberOfBlocksInDirection,
+                                                      numberOfCellsPerBlockInDirection, numberOfCellsPerBlockInDirection, numberOfCellsPerBlockInDirection,
+                                                      dx, 0, false, false,
+                                                      false, false, false,
+                                                      false );
+   // flag field
+   BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field", FieldGhostLayers, false, initFlagField );
+
+   // data fields
+   BlockDataID scalarFieldID         = field::addToStorage< ScalarField_T >( blocks, "scalar field", real_t(0), field::zyxf, FieldGhostLayers );
+   BlockDataID vectorFieldID         = field::addToStorage< Vec3Field_T >( blocks, "vec3 field", Vector3<real_t>(real_t(0)), field::zyxf, FieldGhostLayers );
+   BlockDataID multiComponentFieldID = field::addToStorage< MultiComponentField_T >( blocks, "multi component field", real_t(0), field::zyxf, FieldGhostLayers );
+
+   initScalarField(blocks, scalarFieldID);
+   initVectorField(blocks, vectorFieldID );
+   initMultiComponentField(blocks, multiComponentFieldID );
+
+   // test all interpolators with domain flags everywhere, i.e. without special boundary treatment necessary
+   testNearestNeighborFieldInterpolator(blocks, flagFieldID, scalarFieldID, vectorFieldID, multiComponentFieldID);
+   testTrilinearFieldInterpolator(blocks, flagFieldID, scalarFieldID, vectorFieldID, multiComponentFieldID);
+   testKernelFieldInterpolator(blocks, flagFieldID, scalarFieldID, vectorFieldID, multiComponentFieldID);
+
+   // set some boundary flags in flag field and invalidate the corresponding scalar field values
+   setBoundaryFlags(blocks, flagFieldID, scalarFieldID);
+
+   // test all interpolators' behavior close to boundary cells
+   testNearestNeighborFieldInterpolatorAtBoundary(blocks, flagFieldID, scalarFieldID);
+   testTrilinearFieldInterpolatorAtBoundary(blocks, flagFieldID, scalarFieldID);
+   testKernelFieldInterpolatorAtBoundary(blocks, flagFieldID, scalarFieldID);
+
+   return 0;
+}
+
+} // namespace field_interpolation_tests
+
+int main( int argc, char **argv ){
+   field_interpolation_tests::main(argc, argv);
+}
\ No newline at end of file