Skip to content
Snippets Groups Projects
BoundaryHandling.h 139 KiB
Newer Older
//======================================================================================================================
//
//  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 BoundaryHandling.h
//! \ingroup boundary
//! \author Florian Schornbaum <florian.schornbaum@fau.de>
//
//======================================================================================================================

#pragma once

#include "Boundary.h"

#include "core/Abort.h"
#include "core/DataTypes.h"
#include "core/cell/CellInterval.h"
#include "core/cell/CellSet.h"
#include "core/cell/CellVector.h"
#include "core/config/Config.h"
#include "core/debug/Debug.h"
#include "core/uid/UID.h"
#include "core/uid/UIDGenerators.h"

#include "domain_decomposition/IBlock.h"

#include "field/FlagField.h"

#include "stencil/Directions.h"

#include <tuple>
#include <vector>


namespace walberla {
namespace boundary {



class BHUIDGenerator : public uid::IndexGenerator< BHUIDGenerator, uint_t >{};
typedef UID< BHUIDGenerator > BoundaryHandlingUID;



template< typename FlagField_T, typename Stencil, typename... Boundaries > // Boundaries: all the boundary classes that are considered by this boundary handler
   template< typename F, typename... T >
   friend class BoundaryHandlingCollection;

   typedef FlagField_T                               FlagField;
   typedef typename FlagField_T::flag_t              flag_t;
   typedef typename FlagField_T::const_base_iterator ConstFlagFieldBaseIterator;

   enum Mode { OPTIMIZED_SPARSE_TRAVERSAL, ENTIRE_FIELD_TRAVERSAL };



   class BlockSweep {
   public:
      BlockSweep( const BlockDataID & handling, const uint_t numberOfGhostLayersToInclude = 0 ) :
         handling_( handling ), numberOfGhostLayersToInclude_( numberOfGhostLayersToInclude ) {}
      void operator()( IBlock * block ) const {
         BoundaryHandling * handling = block->getData< BoundaryHandling >( handling_ );
         (*handling)( numberOfGhostLayersToInclude_ );
      }
   protected:
      const BlockDataID handling_;
      const uint_t numberOfGhostLayersToInclude_;
   };



   BoundaryHandling( const std::string & identifier, FlagField_T * const flagField, const flag_t domain, const Boundaries & ... boundaryConditions,
                     const Mode mode = OPTIMIZED_SPARSE_TRAVERSAL );

   bool operator==( const BoundaryHandling & rhs ) const { WALBERLA_CHECK( false, "You are trying to compare boundary handling " << uid_ <<                        // For testing purposes, block data items must be comparable with operator "==".
                                                                                  " with boundary handling " << rhs.getUID() <<                                    // Since instances of type "BoundaryHandling" are registered as block data items,
                                                                                  ".\nHowever, boundary handling instances are not comparable!" ); return false; } // "BoundaryHandling" must implement operator "==". As of right now, comparing
                                                                                                                                                                   // two boundary handling instances will always fail... :-) TODO: fixit?
   bool operator!=( const BoundaryHandling & rhs ) const { return !operator==( rhs ); }

   const BoundaryHandlingUID & getUID() const { return uid_; }

   const FlagField_T * getFlagField() const { return flagField_; }
         FlagField_T * getFlagField()       { return flagField_; }

   flag_t getNearBoundaryFlag() const { return nearBoundary_; } ///< Never (ever) set near boundary flags manually outside the boundary handler!
   flag_t getBoundaryMask() const { return boundary_; }
   flag_t getDomainMask() const { return domain_; }

   inline bool isEmpty       ( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const;
   inline bool isNearBoundary( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const;
   inline bool isBoundary    ( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const;
   inline bool isDomain      ( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const;

   inline bool isEmpty       ( const Cell & cell ) const;
   inline bool isNearBoundary( const Cell & cell ) const;
   inline bool isBoundary    ( const Cell & cell ) const;
   inline bool isDomain      ( const Cell & cell ) const;

   inline bool isEmpty       ( const ConstFlagFieldBaseIterator & it ) const;
   inline bool isNearBoundary( const ConstFlagFieldBaseIterator & it ) const;
   inline bool isBoundary    ( const ConstFlagFieldBaseIterator & it ) const;
   inline bool isDomain      ( const ConstFlagFieldBaseIterator & it ) const;

   inline bool containsBoundaryCondition( const BoundaryUID & uid ) const;
   inline bool containsBoundaryCondition( const FlagUID & flag ) const;
   inline bool containsBoundaryCondition( const flag_t    flag ) const { return ( boundary_ & flag ) == flag; }

   inline flag_t getBoundaryMask( const BoundaryUID & uid ) const { return getBoundaryMask( boundaryConditions_, uid ); }

   template< typename Boundary_T >
   inline const Boundary_T & getBoundaryCondition( const BoundaryUID & uid ) const; ///< You most likely have to call this function via "handling.template getBoundaryCondition< Boundary_T >( uid )".
   template< typename Boundary_T >
   inline       Boundary_T & getBoundaryCondition( const BoundaryUID & uid );       ///< You most likely have to call this function via "handling.template getBoundaryCondition< Boundary_T >( uid )".

   inline BoundaryUID getBoundaryUID( const FlagUID & flag ) const;
   inline BoundaryUID getBoundaryUID( const flag_t    flag ) const;

   inline uint_t numberOfMatchingBoundaryConditions( const flag_t mask ) const;

   inline bool checkConsistency( const uint_t numberOfGhostLayersToInclude = 0 ) const;
          bool checkConsistency( const CellInterval & cells ) const;

   inline void refresh( const uint_t numberOfGhostLayersToInclude = 0 ); // reset near ...
          void refresh( const CellInterval & cells );                    // ... boundary flags

   void refreshOutermostLayer( cell_idx_t thickness = 1 ); // reset near boundary flags in the outermost "inner" layers

   //** Set Domain Cells ***********************************************************************************************
   /*! \name Set Domain Cells */
   //@{
   inline void setDomain(                             const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
   inline void setDomain( const flag_t domainSubFlag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );

   inline void setDomain(                             const CellInterval & cells );
   inline void setDomain( const flag_t domainSubFlag, const CellInterval & cells );

   template< typename CellIterator >
   inline void setDomain(                             const CellIterator & begin, const CellIterator & end );
   template< typename CellIterator >
   inline void setDomain( const flag_t domainSubFlag, const CellIterator & begin, const CellIterator & end );

   inline void forceDomain(                             const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
   inline void forceDomain( const flag_t domainSubFlag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );

   inline void forceDomain(                             const CellInterval & cells );
   inline void forceDomain( const flag_t domainSubFlag, const CellInterval & cells );

   template< typename CellIterator >
   inline void forceDomain(                             const CellIterator & begin, const CellIterator & end );
   template< typename CellIterator >
   inline void forceDomain( const flag_t domainSubFlag, const CellIterator & begin, const CellIterator & end );

   inline void fillWithDomain(                             const uint_t numberOfGhostLayersToInclude = 0 );
   inline void fillWithDomain( const flag_t domainSubFlag, const uint_t numberOfGhostLayersToInclude     );

   inline void fillWithDomain(                             const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
   inline void fillWithDomain( const flag_t domainSubFlag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );

   inline void fillWithDomain(                             const CellInterval & cells );
   inline void fillWithDomain( const flag_t domainSubFlag, const CellInterval & cells );

   template< typename CellIterator >
   inline void fillWithDomain(                             const CellIterator & begin, const CellIterator & end );
   template< typename CellIterator >
   inline void fillWithDomain( const flag_t domainSubFlag, const CellIterator & begin, const CellIterator & end );
   //@}
   //*******************************************************************************************************************

   //** Set Boundary Cells *********************************************************************************************
   /*! \name Set Boundary Cells */
   //@{
   inline shared_ptr<BoundaryConfiguration> createBoundaryConfiguration( const BoundaryUID & uid, const Config::BlockHandle & config ) const;

   inline void setBoundary( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                            const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );
   inline void setBoundary( const flag_t   flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                            const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );

   inline void setBoundary( const FlagUID & flag, const CellInterval & cells, const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );
   inline void setBoundary( const flag_t    flag, const CellInterval & cells, const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );

   template< typename CellIterator >
   inline void setBoundary( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
                            const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );
   template< typename CellIterator >
   inline void setBoundary( const flag_t    flag, const CellIterator & begin, const CellIterator & end,
                            const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );

   inline void forceBoundary( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                              const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );
   inline void forceBoundary( const flag_t    flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                              const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );

   inline void forceBoundary( const FlagUID & flag, const CellInterval & cells, const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );
   inline void forceBoundary( const flag_t    flag, const CellInterval & cells, const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );

   template< typename CellIterator >
   inline void forceBoundary( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
                              const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );
   template< typename CellIterator >
   inline void forceBoundary( const flag_t    flag, const CellIterator & begin, const CellIterator & end,
                              const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );
   //@}
   //*******************************************************************************************************************

   //** Remove Domain Cells ********************************************************************************************
   /*! \name Remove Domain Cells */
   //@{
   inline void removeDomain(                    const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
   inline void removeDomain( const flag_t mask, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );

   inline void removeDomain( const uint_t numberOfGhostLayersToInclude = 0 );
   inline void removeDomain(                    const CellInterval & cells );
   inline void removeDomain( const flag_t mask, const CellInterval & cells );

   template< typename CellIterator >
   inline void removeDomain(                    const CellIterator & begin, const CellIterator & end );
   template< typename CellIterator >
   inline void removeDomain( const flag_t mask, const CellIterator & begin, const CellIterator & end );
   //@}
   //*******************************************************************************************************************

   //** Remove Boundary Cells ******************************************************************************************
   /*! \name Remove Boundary Cells */
   //@{
   inline void removeBoundary(                      const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
   inline void removeBoundary( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
   inline void removeBoundary( const flag_t    mask, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );

   inline void removeBoundary( const uint_t numberOfGhostLayersToInclude = 0 );
   inline void removeBoundary(                       const CellInterval & cells );
   inline void removeBoundary( const FlagUID & flag, const CellInterval & cells );
   inline void removeBoundary( const flag_t    mask, const CellInterval & cells );

   template< typename CellIterator >
   inline void removeBoundary(                       const CellIterator & begin, const CellIterator & end );
   template< typename CellIterator >
   inline void removeBoundary( const FlagUID & flag, const CellIterator & begin, const CellIterator & end );
   template< typename CellIterator >
   inline void removeBoundary( const flag_t    mask, const CellIterator & begin, const CellIterator & end );
   //@}
   //*******************************************************************************************************************

   //** General Flag Handling ******************************************************************************************
   /*! \name General Flag Handling */
   //@{
   inline void setFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                        const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );
   inline void setFlag( const flag_t    flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                        const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );

   inline void setFlag( const FlagUID & flag, const CellInterval & cells,
                        const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );
   inline void setFlag( const flag_t    flag, const CellInterval & cells,
                        const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );

   template< typename CellIterator >
   inline void setFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
                        const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );
   template< typename CellIterator >
   inline void setFlag( const flag_t    flag, const CellIterator & begin, const CellIterator & end,
                        const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );

   inline void forceFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                          const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );
   inline void forceFlag( const flag_t    flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                          const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );

   inline void forceFlag( const FlagUID & flag, const CellInterval & cells,
                          const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );
   inline void forceFlag( const flag_t    flag, const CellInterval & cells,
                          const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );

   template< typename CellIterator >
   inline void forceFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
                          const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );
   template< typename CellIterator >
   inline void forceFlag( const flag_t    flag, const CellIterator & begin, const CellIterator & end,
                          const BoundaryConfiguration & parameter = BoundaryConfiguration::null() );

   inline void removeFlag( const FlagUID & flag, const uint_t numberOfGhostLayersToInclude = 0 );
   inline void removeFlag( const flag_t    flag, const uint_t numberOfGhostLayersToInclude = 0 );

   inline void removeFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
   inline void removeFlag( const flag_t    flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );

   inline void removeFlag( const FlagUID & flag, const CellInterval & cells );
   inline void removeFlag( const flag_t    flag, const CellInterval & cells );

   template< typename CellIterator >
   inline void removeFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end );
   template< typename CellIterator >
   inline void removeFlag( const flag_t    flag, const CellIterator & begin, const CellIterator & end );
   //@}
   //*******************************************************************************************************************

   //** Clear Cells ****************************************************************************************************
   /*! \name Clear Cells */
   //@{
   inline void clear( const uint_t numberOfGhostLayersToInclude = 0 );
   inline void clear( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
   inline void clear( const CellInterval & cells );
   template< typename CellIterator >
   inline void clear( const CellIterator & begin, const CellIterator & end );
   //@}
   //*******************************************************************************************************************

   //** Boundary Treatment *********************************************************************************************
   /*! \name Boundary Treatment */
   //@{
   static BlockSweep getBlockSweep( const BlockDataID handling, const uint_t numberOfGhostLayersToInclude = 0 )
      { return BlockSweep( handling, numberOfGhostLayersToInclude ); }

   inline void operator()( const uint_t numberOfGhostLayersToInclude = 0 );
   inline void operator()( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
   inline void operator()( const CellInterval & cells );
   template< typename CellIterator >
   inline void operator()( const CellIterator & begin, const CellIterator & end );

   inline void beforeBoundaryTreatment() { beforeBoundaryTreatment( boundaryConditions_ ); }
   inline void  afterBoundaryTreatment() {  afterBoundaryTreatment( boundaryConditions_ ); }
   //@}
   //*******************************************************************************************************************

   //** Pack / Unpack boundary handling ********************************************************************************
   /*! \name Pack / Unpack boundary handling */
   //@{
   template< typename Buffer_T >
   void   pack( Buffer_T & buffer, const CellInterval & interval, const bool assumeIdenticalFlagMapping = true ) const;
   template< typename Buffer_T >
   void unpack( Buffer_T & buffer, const CellInterval & interval, const bool assumeIdenticalFlagMapping = true );

   template< typename Buffer_T >
   void   pack( Buffer_T & buffer, stencil::Direction direction, const uint_t numberOfLayers = 1, const bool assumeIdenticalFlagMapping = true ) const;
   template< typename Buffer_T >
   void unpack( Buffer_T & buffer, stencil::Direction direction, const uint_t numberOfLayers = 1, const bool assumeIdenticalFlagMapping = true );
   //@}
   //*******************************************************************************************************************

   inline void        toStream( std::ostream & os ) const;
   inline std::string toString() const;

private:

   CellInterval getGhostLayerCellInterval( const uint_t numberOfGhostLayersToInclude ) const;

   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   typename std::enable_if<(N!=-1), void>::type setupBoundaryConditions(       BoundariesTuple & boundaryConditions );
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1  >
   typename std::enable_if<(N==-1), void>::type setupBoundaryConditions( const BoundariesTuple & ) const {}

   inline std::vector< BoundaryUID > getBoundaryUIDs() const;
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), void>::type getBoundaryUIDs( const BoundariesTuple & boundaryConditions, std::vector< BoundaryUID > & uids ) const;
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), void>::type getBoundaryUIDs( const BoundariesTuple &, std::vector< BoundaryUID > & ) const {}

   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), BoundaryUID>::type getBoundaryUID( const BoundariesTuple & boundaryConditions, const flag_t flag ) const;
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), BoundaryUID>::type getBoundaryUID( const BoundariesTuple &, const flag_t flagUID ) const;

   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), bool>::type containsBoundaryCondition( const BoundariesTuple & boundaryConditions, const BoundaryUID & uid ) const;
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), bool>::type containsBoundaryCondition( const BoundariesTuple &, const BoundaryUID & ) const { return false; }

   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), flag_t>::type getBoundaryMask( const BoundariesTuple & boundaryConditions, const BoundaryUID & uid ) const;
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), flag_t>::type getBoundaryMask( const BoundariesTuple &, const BoundaryUID & ) const { return numeric_cast<flag_t>(0); }

   //** Get Boundary Class (private helper functions) ******************************************************************
   /*! \name Get Boundary Class (private helper functions) */
   //@{

   // matching type (-> Boundary_T) not yet found ...

   template< typename Boundary_T, typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline const typename std::enable_if<(N!=0), Boundary_T>::type & getBoundaryCondition( const BoundaryUID & uid, const BoundariesTuple & boundaryConditions,
                                                   typename std::enable_if< std::is_same< Boundary_T, typename std::tuple_element<N, BoundariesTuple>::type >::value >::type* /*dummy*/ = 0 ) const
      if( uid == std::get<N>( boundaryConditions ).getUID() )
         return std::get<N>( boundaryConditions );
         return getBoundaryCondition_TypeExists< Boundary_T, BoundariesTuple, N-1 >( uid, boundaryConditions );
   template< typename Boundary_T, typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline const typename std::enable_if<(N==0), Boundary_T>::type & getBoundaryCondition( const BoundaryUID & uid, const BoundariesTuple & boundaryConditions,
                                                   typename std::enable_if< std::is_same< Boundary_T, typename std::tuple_element<N, BoundariesTuple>::type >::value >::type* /*dummy*/ = 0 ) const
      if( uid == std::get<N>( boundaryConditions ).getUID() )
         return std::get<N>( boundaryConditions );
      else
         WALBERLA_ABORT( "The requested boundary condition " << uid.getIdentifier() << " is not part of this boundary handling." );

#ifdef __IBMCPP__
      return *(reinterpret_cast< Boundary_T * >( NULL )); // silencing incorrect IBM compiler warning
#endif
   }

   template< typename Boundary_T, typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline const typename std::enable_if<(N!=0), Boundary_T>::type & getBoundaryCondition( const BoundaryUID & uid, const BoundariesTuple & boundaryConditions,
                                                   typename std::enable_if< std::is_same< typename std::is_same< Boundary_T, typename std::tuple_element<N, BoundariesTuple>::type >::type,
                                                                                              std::false_type >::value >::type* /*dummy*/ = 0,
                                                   typename std::enable_if< (N>0) >::type* /*dummy*/ = 0 ) const
      return getBoundaryCondition< Boundary_T, BoundariesTuple, N-1 >( uid, boundaryConditions );
   template< typename Boundary_T, typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline const typename std::enable_if<(N==0), Boundary_T>::type & getBoundaryCondition( const BoundaryUID & /*uid*/, const BoundariesTuple & /*boundaryConditions*/,
                                                   typename std::enable_if< std::is_same< typename std::is_same< Boundary_T, typename std::tuple_element<0, BoundariesTuple>::type >::type,
                                                                                              std::false_type >::value >::type* /*dummy*/ = 0 ) const
   {
      static_assert( sizeof(Boundary_T) == 0, "The requested boundary class is not part of this boundary handling." );
   }

   // matching type (-> Boundary_T) exists!

   template< typename Boundary_T, typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline const typename std::enable_if<(N!=0), Boundary_T>::type & getBoundaryCondition_TypeExists( const BoundaryUID & uid, const BoundariesTuple & boundaryConditions,
                                                              typename std::enable_if< std::is_same< Boundary_T, typename std::tuple_element<N, BoundariesTuple>::type >::value >::type* /*dummy*/ = 0 ) const
      if( uid == std::get<N>( boundaryConditions ).getUID() )
         return std::get<N>( boundaryConditions );
         return getBoundaryCondition_TypeExists< Boundary_T, BoundariesTuple, N-1 >( uid, boundaryConditions );
   template< typename Boundary_T, typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline const typename std::enable_if<(N==0), Boundary_T>::type & getBoundaryCondition_TypeExists( const BoundaryUID & uid, const BoundariesTuple & boundaryConditions,
                                                              typename std::enable_if< std::is_same< Boundary_T, typename std::tuple_element<0, BoundariesTuple>::type >::value >::type* /*dummy*/ = 0 ) const
      if( uid == std::get<0>( boundaryConditions ).getUID() )
         return std::get<0>( boundaryConditions );
      else
         WALBERLA_ABORT( "The requested boundary condition " << uid.getIdentifier() << " is not part of this boundary handling." );

#ifdef __IBMCPP__
      return *(reinterpret_cast< Boundary_T * >( NULL )); // silencing incorrect IBM compiler warning
#endif
   }

   template< typename Boundary_T, typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline const typename std::enable_if<(N!=0), Boundary_T>::type & getBoundaryCondition_TypeExists( const BoundaryUID & uid, const BoundariesTuple & boundaryConditions,
                                                              typename std::enable_if< std::is_same< typename std::is_same< Boundary_T, typename std::tuple_element<N, BoundariesTuple>::type >::type,
                                                                                                         std::false_type >::value >::type* /*dummy*/ = 0 ) const
      return getBoundaryCondition_TypeExists< Boundary_T, BoundariesTuple, N-1 >( uid, boundaryConditions );
   template< typename Boundary_T, typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline const typename std::enable_if<(N==0), Boundary_T>::type & getBoundaryCondition_TypeExists( const BoundaryUID & uid, const BoundariesTuple & /*boundaryConditions*/,
                                                              typename std::enable_if< std::is_same< typename std::is_same< Boundary_T, typename std::tuple_element<0, BoundariesTuple>::type >::type,
                                                                                                         std::false_type >::value >::type* /*dummy*/ = 0 ) const
   {
      WALBERLA_ABORT( "The requested boundary condition " << uid.getIdentifier() << " is not part of this boundary handling." );

#ifdef __IBMCPP__
      return *(reinterpret_cast< Boundary_T * >( NULL )); // silencing incorrect IBM compiler warning
#endif
   }

   //*******************************************************************************************************************

   inline uint_t numberOfMatchingBoundaryConditions( const BoundaryUID & uid ) const;

   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), uint_t>::type numberOfMatchingBoundaryConditions( const BoundariesTuple & boundaryConditions, const BoundaryUID & uid ) const;
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), uint_t>::type numberOfMatchingBoundaryConditions( const BoundariesTuple &, const BoundaryUID & ) const { return uint_c(0); }
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), uint_t>::type numberOfMatchingBoundaryConditions( const BoundariesTuple & boundaryConditions, const flag_t mask ) const;
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), uint_t>::type numberOfMatchingBoundaryConditions( const BoundariesTuple &, const flag_t ) const { return uint_c(0); }

   inline bool checkFlagField( const uint_t numberOfGhostLayersToInclude = 0 ) const;

   inline void addDomain( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const flag_t domain );

   //** Set Boundary Cells (private helper functions) ******************************************************************
   /*! \name Set Boundary Cells (private helper functions) */
   //@{
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), shared_ptr<BoundaryConfiguration>>::type createBoundaryConfiguration( const BoundariesTuple & boundaryConditions,
                                                                         const BoundaryUID & uid, const Config::BlockHandle & config ) const;
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), shared_ptr<BoundaryConfiguration>>::type createBoundaryConfiguration( const BoundariesTuple &, const BoundaryUID & uid,
                                                                         const Config::BlockHandle & ) const;

   inline void addNearBoundary( const CellInterval & cells );
   inline void addBoundary( const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );

   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), void>::type setBoundary(       BoundariesTuple & boundaryConditions, const flag_t flag,
                                                                                        const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                                        const BoundaryConfiguration & parameter );
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), void>::type setBoundary( const BoundariesTuple &, const flag_t, const cell_idx_t, const cell_idx_t, const cell_idx_t,
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), void>::type setBoundary(       BoundariesTuple & boundaryConditions, const flag_t flag, const CellInterval & cells,
                                                                                        const BoundaryConfiguration & parameter );
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), void>::type setBoundary( const BoundariesTuple &, const flag_t, const CellInterval &, const BoundaryConfiguration & ) const;
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), void>::type setBoundary(       BoundariesTuple & boundaryConditions, const flag_t flag, const CellVector & cells,
                                                                                        const BoundaryConfiguration & parameter );
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), void>::type setBoundary( const BoundariesTuple &, const flag_t, const CellVector &, const BoundaryConfiguration & ) const;
   //@}
   //*******************************************************************************************************************

   //** Remove Boundary Cells (private helper functions) ***************************************************************
   /*! \name Remove Boundary Cells (private helper functions) */
   //@{
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), void>::type removeBoundary(       BoundariesTuple & boundaryConditions, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                                           const bool checkNearBoundaryFlags = true );
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), void>::type removeBoundary( const BoundariesTuple &, const cell_idx_t, const cell_idx_t, const cell_idx_t, const bool ) const { WALBERLA_CHECK( false ); }
   //@}
   //*******************************************************************************************************************

   //** Boundary Treatment (private helper functions) ******************************************************************
   /*! \name Boundary Treatment (private helper functions) */
   //@{
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), void>::type treatDirection( BoundariesTuple & boundaryConditions, const uint_t index,
                               const std::vector< std::vector< std::pair< Cell, stencil::Direction > > > & cellDirectionPairs );
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), void>::type treatDirection( const BoundariesTuple &, const uint_t,
                               const std::vector< std::vector< std::pair< Cell, stencil::Direction > > > & ) const {}

   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), void>::type treatDirection(       BoundariesTuple & boundaryConditions, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                                           const stencil::Direction dir,
                                                                                           const cell_idx_t nx, const cell_idx_t ny, const cell_idx_t nz );
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), void>::type treatDirection( const BoundariesTuple & , const cell_idx_t, const cell_idx_t, const cell_idx_t, const stencil::Direction,
                                                                  const cell_idx_t, const cell_idx_t, const cell_idx_t ) const { WALBERLA_CHECK( false ); }

   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), void>::type beforeBoundaryTreatment(       BoundariesTuple & boundaryConditions );
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), void>::type beforeBoundaryTreatment( const BoundariesTuple & ) const {}
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), void>::type afterBoundaryTreatment(       BoundariesTuple & boundaryConditions );
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), void>::type afterBoundaryTreatment( const BoundariesTuple & ) const {}
   //@}
   //*******************************************************************************************************************

   //** Pack / Unpack boundary handling (private helper functions) *****************************************************
   /*! \name Pack / Unpack boundary handling (private helper functions) */
   //@{
   inline std::map< std::string, flag_t > getFlagMapping() const;

   template< typename Buffer_T >
   std::vector< flag_t > getNeighborFlagMapping( Buffer_T & buffer, const bool assumeIdenticalFlagMapping, bool & identicalFlagMapping ) const;

   void translateMask( flag_t & mask, const std::vector< flag_t > & flagMapping ) const;

   CellInterval   getPackingInterval( stencil::Direction direction, const uint_t numberOfLayers ) const;
   CellInterval getUnpackingInterval( stencil::Direction direction, const uint_t numberOfLayers ) const;

   template< typename Buffer_T >
   inline void pack( Buffer_T & buffer, const flag_t mask, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const;

   template< typename Buffer_T, typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), void>::type pack( const BoundariesTuple & boundaryConditions, Buffer_T & buffer,
                     const flag_t mask, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const;
   template< typename Buffer_T, typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   typename std::enable_if<(N==-1), void>::type pack( const BoundariesTuple &, Buffer_T &, const flag_t, const cell_idx_t, const cell_idx_t, const cell_idx_t )
              const { WALBERLA_CHECK( false ); }

   template< typename Buffer_T >
   inline void unpack( Buffer_T & buffer, const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );

   template< typename Buffer_T >
   inline void unpackBoundary( Buffer_T & buffer, const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );

   template< typename Buffer_T, typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N!=-1), void>::type unpackBoundary( BoundariesTuple & boundaryConditions, Buffer_T & buffer, const flag_t flag,
                               const cell_idx_t x, const cell_idx_t y, const cell_idx_t z );
   template< typename Buffer_T, typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   inline typename std::enable_if<(N==-1), void>::type unpackBoundary( const BoundariesTuple &, Buffer_T &, const flag_t, const cell_idx_t, const cell_idx_t, const cell_idx_t)
                               const { WALBERLA_CHECK( false ); }
   //@}
   //*******************************************************************************************************************

   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   typename std::enable_if<(N!=-1), void>::type getBoundaryConditions( const BoundariesTuple & boundaryConditions, std::vector< std::string > & bcs ) const;
   template< typename BoundariesTuple, int N = std::tuple_size<BoundariesTuple>::value - 1 >
   typename std::enable_if<(N==-1), void>::type getBoundaryConditions( const BoundariesTuple &, std::vector< std::string > & ) const {}

   template< typename T > static void valueToStream( std::ostream & os, const T       value ) { os << value; }
                          static void valueToStream( std::ostream & os, const  int8_t value ) { os <<  int_c( value ); }
                          static void valueToStream( std::ostream & os, const uint8_t value ) { os << uint_c( value ); }



   const BoundaryHandlingUID uid_;

   FlagField_T * const flagField_;

   const CellInterval innerBB_;
   const CellInterval outerBB_;

   const flag_t nearBoundary_; // "near boundary" flags are only set in cells that are marked as "domain" cells and that are ...
         flag_t boundary_;     // ... located next to "boundary" cells. "boundary_" and "domain_" must be disjoint!
   const flag_t domain_;       // Cells that are marked as neither "boundary" cells nor "domain" cells are ignored.

   const Mode mode_;

   bool dirty_; //!< Every time "near boundary" flags are set or removed, dirty_ is set to true.
                /*!< Triggers a rebuild of vector cellDirectionPairs_ in "OPTIMIZED_SPARSE_TRAVERSAL" mode. */

   std::vector< flag_t > bcMaskMapping_; // boundary condition index (position in tuple) to associated mask mapping

   std::vector< bool > rebuildCellDirectionPairs_;
   std::vector< std::vector< std::vector< std::pair< Cell, stencil::Direction > > > > cellDirectionPairs_; // 1st vector: numberOfGhostLayersToInclude
                                                                                                           // 2nd vector: boundary condition index
                                                                                                           // 3rd vector: vector of cell<->direction pairs
   typedef std::tuple<Boundaries...> Tuple;
   Tuple boundaryConditions_;
   bool  threadSafeBCs_;

}; // class BoundaryHandling



template< typename FlagField_T, typename Stencil, typename... Boundaries >
BoundaryHandling< FlagField_T, Stencil, Boundaries... >::BoundaryHandling( const std::string & identifier, FlagField_T * const flagField,
                                                                   const flag_t domain, const Boundaries & ... boundaryConditions, const Mode mode ) :

   uid_( identifier ),
   flagField_( flagField ),
   innerBB_( 1 - cell_idx_c( flagField_->nrOfGhostLayers() ), 1 - cell_idx_c( flagField_->nrOfGhostLayers() ), 1 - cell_idx_c( flagField_->nrOfGhostLayers() ),
             cell_idx_c( flagField_->xSize() + flagField_->nrOfGhostLayers() ) - 2, cell_idx_c( flagField_->ySize() + flagField_->nrOfGhostLayers() ) - 2,
             cell_idx_c( flagField_->zSize() + flagField_->nrOfGhostLayers() ) - 2 ),
   outerBB_( -cell_idx_c( flagField_->nrOfGhostLayers() ), -cell_idx_c( flagField_->nrOfGhostLayers() ), -cell_idx_c( flagField_->nrOfGhostLayers() ),
             cell_idx_c( flagField_->xSize() + flagField_->nrOfGhostLayers() ) - 1, cell_idx_c( flagField_->ySize() + flagField_->nrOfGhostLayers() ) - 1,
             cell_idx_c( flagField_->zSize() + flagField_->nrOfGhostLayers() ) - 1 ),
   nearBoundary_( flagField_->registerFlag( std::string( "near boundary (" ) + identifier + std::string( ")" ) ) ),
   boundary_( 0 ),
   domain_( domain ),
   mode_( mode ),
   dirty_( false ),
   boundaryConditions_( std::make_tuple(boundaryConditions...) ),
   threadSafeBCs_( true )
{
   setupBoundaryConditions( boundaryConditions_ );

   if( flagField_->nrOfGhostLayers() < 1 )
      WALBERLA_ABORT( "The flag field passed to the boundary handling \"" << identifier << "\" must contain at least one ghost layer!" );

   rebuildCellDirectionPairs_.resize( flagField_->nrOfGhostLayers(), false );
   cellDirectionPairs_.resize( flagField_->nrOfGhostLayers() );

   WALBERLA_CHECK_EQUAL( nearBoundary_ & boundary_, flag_t(0), "The near boundary flag must not be identical to a flag used for marking a boundary cell.\n"
                                                               "This check failed for boundary handling " << uid_ << "!" );
   WALBERLA_CHECK_EQUAL( nearBoundary_ & domain_,   flag_t(0), "The near boundary flag must not be identical to a flag used for marking a domain cell.\n"
                                                               "This check failed for boundary handling " << uid_ << "!" );
   WALBERLA_CHECK_EQUAL( boundary_ & domain_,       flag_t(0), "Flags used for marking domain cells must be different to flags used for marking boundary cells.\n"
                                                               "This check failed for boundary handling " << uid_ << "!" );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::isEmpty( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
{
   return !flagField_->isPartOfMaskSet( x, y, z, boundary_ ) && !flagField_->isPartOfMaskSet( x, y, z, domain_ );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::isNearBoundary( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::isBoundary( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
{
   return flagField_->isPartOfMaskSet( x, y, z, boundary_ );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::isDomain( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::isEmpty( const Cell & cell ) const
{
   return !flagField_->isPartOfMaskSet( cell.x(), cell.y(), cell.z(), boundary_ ) && 
          !flagField_->isPartOfMaskSet( cell.x(), cell.y(), cell.z(), domain_ );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::isNearBoundary( const Cell & cell ) const
{
   return flagField_->isFlagSet( cell.x(), cell.y(), cell.z(), nearBoundary_ );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::isBoundary( const Cell & cell ) const
{
   return flagField_->isPartOfMaskSet( cell.x(), cell.y(), cell.z(), boundary_ );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::isDomain( const Cell & cell ) const
{
   return flagField_->isPartOfMaskSet( cell.x(), cell.y(), cell.z(), domain_ );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::isEmpty( const ConstFlagFieldBaseIterator & it ) const
{
   WALBERLA_ASSERT_EQUAL( it.getField(), flagField_ );

   return !isPartOfMaskSet( it, boundary_ ) && !isPartOfMaskSet( it, domain_ );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::isNearBoundary( const ConstFlagFieldBaseIterator & it ) const
{
   WALBERLA_ASSERT_EQUAL( it.getField(), flagField_ );

   return isFlagSet( it, nearBoundary_ );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::isBoundary( const ConstFlagFieldBaseIterator & it ) const
{
   WALBERLA_ASSERT_EQUAL( it.getField(), flagField_ );

   return isPartOfMaskSet( it, boundary_ );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::isDomain( const ConstFlagFieldBaseIterator & it ) const
{
   WALBERLA_ASSERT_EQUAL( it.getField(), flagField_ );

   return isPartOfMaskSet( it, domain_ );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::containsBoundaryCondition( const BoundaryUID & uid ) const
{
   return containsBoundaryCondition( boundaryConditions_, uid );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::containsBoundaryCondition( const FlagUID & flag ) const
{
   if( !flagField_->flagExists( flag ) )
      return false;

   return containsBoundaryCondition( flagField_->getFlag( flag ) );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline const Boundary_T & BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getBoundaryCondition( const BoundaryUID & uid ) const
   return getBoundaryCondition< Boundary_T, std::tuple< Boundaries... > >( uid, boundaryConditions_ );
template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline Boundary_T & BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getBoundaryCondition( const BoundaryUID & uid )
{
   return const_cast< Boundary_T & >( static_cast< const BoundaryHandling * >( this )->template getBoundaryCondition< Boundary_T >( uid ) );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline BoundaryUID BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getBoundaryUID( const FlagUID & flag ) const
{
   WALBERLA_ASSERT( flagField_->flagExists( flag ) );

   return getBoundaryUID( flagField_->getFlag( flag ) );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline BoundaryUID BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getBoundaryUID( const flag_t flag ) const
{
   WALBERLA_ASSERT( field::isFlag( flag ) );
   WALBERLA_ASSERT( flagField_->isRegistered( flag ) );

   return getBoundaryUID( boundaryConditions_, flag );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline uint_t BoundaryHandling< FlagField_T, Stencil, Boundaries... >::numberOfMatchingBoundaryConditions( const flag_t mask ) const
{
   return numberOfMatchingBoundaryConditions( boundaryConditions_, mask );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::checkConsistency( const uint_t numberOfGhostLayersToInclude ) const
{
   CellInterval cells = getGhostLayerCellInterval( numberOfGhostLayersToInclude );
   return checkConsistency( cells );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::checkConsistency( const CellInterval & cells ) const
{
   CellInterval localCells( innerBB_ );
   localCells.intersect( cells );

   for( auto cell = flagField_->beginSliceXYZ( localCells ); cell != flagField_->end(); ++cell )
   {
      if( isPartOfMaskSet( cell, boundary_ ) ) // if the cell is marked as a boundary cell ...
      {
         // ... near boundary & domain must not be set
         WALBERLA_CHECK( !isPartOfMaskSet( cell, domain_ ) && !isPartOfMaskSet( cell, nearBoundary_ ) );
         if( isPartOfMaskSet( cell, domain_ ) || isPartOfMaskSet( cell, nearBoundary_ ) )
            return false;

         // ... exactly one boundary condition must match
         WALBERLA_CHECK_EQUAL( numberOfMatchingBoundaryConditions( *cell ), uint_t(1) );
         if( numberOfMatchingBoundaryConditions( *cell ) != 1 )
            return false;

         // ... only one bit of the mask may be set
         WALBERLA_CHECK( field::isFlag( *cell & boundary_ ) );
         if( !field::isFlag( *cell & boundary_ ) )
            return false;
      }
      else if( isPartOfMaskSet( cell, domain_ ) ) // if the cell is marked as a domain cell ...
      {
         // ... it must not be marked as a boundary cell
         WALBERLA_CHECK( !isPartOfMaskSet( cell, boundary_ ) );
         if( isPartOfMaskSet( cell, boundary_ ) )
            return false;

         bool boundaryNeighbor( false );
         for( auto d = Stencil::beginNoCenter(); d != Stencil::end(); ++d )
         {
            const flag_t & nMask = cell.neighbor(*d);
            if( isPartOfMaskSet( nMask, boundary_ ) )
            {
               // if a neighbor is marked as a boundary cell, exactly one boundary condition must match
               WALBERLA_CHECK_EQUAL( numberOfMatchingBoundaryConditions( nMask ), uint_t(1) );
               if( numberOfMatchingBoundaryConditions( nMask ) != 1 )
                  return false;
               boundaryNeighbor = true;
            }
         }

         // if at least one neighbor is marked as a boundary cell, near boundary must be set - otherwise it must not be set
         WALBERLA_CHECK( ( boundaryNeighbor && isFlagSet( cell, nearBoundary_ ) ) || ( !boundaryNeighbor && !isFlagSet( cell, nearBoundary_ ) ) );
         if( ( boundaryNeighbor && !isFlagSet( cell, nearBoundary_ ) ) || ( !boundaryNeighbor && isFlagSet( cell, nearBoundary_ ) ) )
            return false;

         // ... only one bit of the mask may be set
         WALBERLA_CHECK( field::isFlag( *cell & domain_ ) );
         if( !field::isFlag( *cell & domain_ ) )
            return false;
      }
   }

   return true;
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::refresh( const uint_t numberOfGhostLayersToInclude )
{
   CellInterval cells = getGhostLayerCellInterval( numberOfGhostLayersToInclude );
   refresh( cells );
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::refresh( const CellInterval & cells )
{
   CellInterval localCells( innerBB_ );
   localCells.intersect( cells );

   for( auto cell = flagField_->beginSliceXYZ( localCells ); cell != flagField_->end(); ++cell )
   {
      if( isPartOfMaskSet( cell, domain_ ) )
      {
         field::removeFlag( cell, nearBoundary_ );

         for( auto d = Stencil::beginNoCenter(); d != Stencil::end(); ++d )
         {
            if( isPartOfMaskSet( cell.neighbor(*d), boundary_ ) )
            {
               addFlag( cell, nearBoundary_ );
               break;
            }
         }
      }
   }

   dirty_ = true;
}



template< typename FlagField_T, typename Stencil, typename... Boundaries >
void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::refreshOutermostLayer( cell_idx_t thickness )
{
   uint_t extent = std::min( std::min( innerBB_.xSize(), innerBB_.ySize() ), innerBB_.zSize() );
   WALBERLA_ASSERT_GREATER( extent, uint_t(0) );

   if( extent == 1 )
   {
      refresh();
      return;
   }

   WALBERLA_ASSERT_GREATER_EQUAL( thickness, cell_idx_t(1) );
   WALBERLA_ASSERT_LESS_EQUAL( thickness, cell_idx_c( extent ) / cell_idx_t(2) );

   const cell_idx_t one = numeric_cast<cell_idx_t>(1);
   thickness -= one;

   CellInterval xlow ( innerBB_.xMin(), innerBB_.yMin(), innerBB_.zMin(), innerBB_.xMin() + thickness, innerBB_.yMax(), innerBB_.zMax() );
   CellInterval xhigh( innerBB_.xMax() - thickness, innerBB_.yMin(), innerBB_.zMin(), innerBB_.xMax(), innerBB_.yMax(), innerBB_.zMax() );

   CellInterval ylow ( innerBB_.xMin() + thickness + one, innerBB_.yMin(), innerBB_.zMin(),
                       innerBB_.xMax() - thickness - one, innerBB_.yMin() + thickness, innerBB_.zMax() );
   CellInterval yhigh( innerBB_.xMin() + thickness + one, innerBB_.yMax() - thickness, innerBB_.zMin(),
                       innerBB_.xMax() - thickness - one, innerBB_.yMax(), innerBB_.zMax() );