diff --git a/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp b/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp
index 42cbcd96183881ed58ce5e718fedd0e9d371214b..750a8e96459db775bf75fa9e4ec3c79659429f70 100644
--- a/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp
+++ b/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/AMRSedimentSettling.cpp
@@ -106,8 +106,7 @@ typedef lbm::NoSlip< LatticeModel_T, flag_t > NoSlip_T;
 
 typedef pe_coupling::CurvedLinear< LatticeModel_T, FlagField_T > MO_T;
 
-typedef boost::tuples::tuple< NoSlip_T, MO_T > BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, MO_T > BoundaryHandling_T;
 
 typedef std::tuple<pe::Sphere, pe::Ellipsoid, pe::Plane> BodyTypeTuple;
 
@@ -590,8 +589,8 @@ BoundaryHandling_T * MyBoundaryHandling::initialize( IBlock * const block )
    WALBERLA_CHECK_NOT_NULLPTR( blocksPtr );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
-                                                           boost::tuples::make_tuple( NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
-                                                                                      MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *blocksPtr, *block ) ),
+                                                           NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+                                                           MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *blocksPtr, *block ),
                                                            BoundaryHandling_T::Mode::ENTIRE_FIELD_TRAVERSAL);
 
    handling->fillWithDomain( FieldGhostLayers );
diff --git a/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/WorkloadEvaluation.cpp b/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/WorkloadEvaluation.cpp
index 239c024bb39b958721d44d013ecaec6b43bfcf2f..588007edfb0baa0d18b8d843ec3cc86b3f58477a 100644
--- a/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/WorkloadEvaluation.cpp
+++ b/apps/benchmarks/AdaptiveMeshRefinementFluidParticleCoupling/WorkloadEvaluation.cpp
@@ -71,8 +71,6 @@
 #include "field/vtk/all.h"
 #include "lbm/vtk/all.h"
 
-#include <boost/tuple/tuple.hpp>
-
 #include <vector>
 #include <iomanip>
 #include <iostream>
@@ -102,8 +100,7 @@ const uint_t FieldGhostLayers = 1;
 // boundary handling
 typedef pe_coupling::CurvedLinear< LatticeModel_T, FlagField_T > MO_CLI_T;
 
-typedef boost::tuples::tuple<MO_CLI_T >               BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, MO_CLI_T > BoundaryHandling_T;
 
 typedef std::tuple<pe::Sphere, pe::Ellipsoid, pe::Plane> BodyTypeTuple;
 
@@ -152,7 +149,7 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    BoundaryHandling_T::Mode mode = (useEntireFieldTraversal_) ? BoundaryHandling_T::Mode::ENTIRE_FIELD_TRAVERSAL : BoundaryHandling_T::Mode::OPTIMIZED_SPARSE_TRAVERSAL ;
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "fixed obstacle boundary handling", flagField, fluid,
-                                                           boost::tuples::make_tuple( MO_CLI_T ( "MO_CLI", MO_CLI_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ),
+                                                           MO_CLI_T ( "MO_CLI", MO_CLI_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
                                                            mode);
 
    // boundaries are set by mapping the planes into the domain
diff --git a/apps/benchmarks/CouetteFlow/CouetteFlow.cpp b/apps/benchmarks/CouetteFlow/CouetteFlow.cpp
index 05126f0202fecd9662a03202b8572fa1a526c98b..9ce1da4a306fdfc381026928ce534a14efd17ec3 100644
--- a/apps/benchmarks/CouetteFlow/CouetteFlow.cpp
+++ b/apps/benchmarks/CouetteFlow/CouetteFlow.cpp
@@ -373,9 +373,7 @@ public:
    typedef lbm::NoSlip< LatticeModel_T, flag_t >    NoSlip_T;
    typedef lbm::SimpleUBB< LatticeModel_T, flag_t > UBB_T;
 
-   typedef boost::tuples::tuple< NoSlip_T, UBB_T > BoundaryConditions_T;
-
-   typedef BoundaryHandling< FlagField_T, typename Types<LatticeModel_T>::Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+   typedef BoundaryHandling< FlagField_T, typename Types<LatticeModel_T>::Stencil_T, NoSlip_T, UBB_T > BoundaryHandling_T;
 
 
 
@@ -405,8 +403,8 @@ MyBoundaryHandling<LatticeModel_T>::operator()( IBlock * const block ) const
    const flag_t fluid = flagField->registerFlag( Fluid_Flag );
 
    return new BoundaryHandling_T( "boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
-                                       UBB_T( "velocity bounce back", UBB_Flag, pdfField, topVelocity_, real_t(0), real_t(0) ) ) );
+                                  NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
+                                  UBB_T( "velocity bounce back", UBB_Flag, pdfField, topVelocity_, real_t(0), real_t(0) ) );
 }
 
 
diff --git a/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp b/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
index 52e9bf90a338dcdb7ee0528169b3182da8e61815..18b847692335722c3560a584f3b3a56669dfd6d3 100644
--- a/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
+++ b/apps/benchmarks/ForcesOnSphereNearPlaneInShearFlow/ForcesOnSphereNearPlaneInShearFlow.cpp
@@ -95,8 +95,7 @@ const uint_t FieldGhostLayers = 4;
 typedef pe_coupling::SimpleBB< LatticeModel_T, FlagField_T > MO_SBB_T;
 typedef pe_coupling::CurvedLinear< LatticeModel_T, FlagField_T > MO_CLI_T;
 
-typedef boost::tuples::tuple< MO_SBB_T, MO_CLI_T > BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, MO_SBB_T, MO_CLI_T > BoundaryHandling_T;
 
 typedef std::tuple< pe::Sphere, pe::Plane > BodyTypeTuple;
 
@@ -238,8 +237,8 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( MO_SBB_T( "MO_SBB", MO_SBB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
-                                    MO_CLI_T( "MO_CLI", MO_CLI_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+                                                           MO_SBB_T( "MO_SBB", MO_SBB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
+                                                           MO_CLI_T( "MO_CLI", MO_CLI_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) );
 
    // boundary conditions are set by mapping the (moving) planes into the domain
 
diff --git a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
index e7743f7959ecde538ce55734bbdc56da60371ac6..de0647508960abe3cb40b4e51e8ec8e7e0ee493d 100644
--- a/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
+++ b/apps/benchmarks/MotionSingleHeavySphere/MotionSingleHeavySphere.cpp
@@ -98,8 +98,7 @@ typedef pe_coupling::SimpleBB< LatticeModel_T, FlagField_T >        MEM_BB_T;
 typedef pe_coupling::CurvedLinear< LatticeModel_T, FlagField_T >    MEM_CLI_T;
 typedef pe_coupling::CurvedQuadratic< LatticeModel_T, FlagField_T > MEM_MR_T;
 
-typedef boost::tuples::tuple< UBB_T, Outlet_T, MEM_BB_T, MEM_CLI_T, MEM_MR_T > BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, UBB_T, Outlet_T, MEM_BB_T, MEM_CLI_T, MEM_MR_T > BoundaryHandling_T;
 
 using BodyTypeTuple = std::tuple<pe::Sphere>;
 
@@ -199,11 +198,11 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( UBB_T( "UBB", UBB_Flag, pdfField, velocity_),
+                                    UBB_T( "UBB", UBB_Flag, pdfField, velocity_),
                                     Outlet_T( "Outlet", Outlet_Flag, pdfField, real_t(1) ),
                                     MEM_BB_T (  "MEM_BB",  MEM_BB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
                                     MEM_CLI_T( "MEM_CLI", MEM_CLI_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
-                                    MEM_MR_T (  "MEM_MR",  MEM_MR_Flag, pdfField, flagField, bodyField, fluid, *storage, *block, pdfFieldPreCol ) )  );
+                                    MEM_MR_T (  "MEM_MR",  MEM_MR_Flag, pdfField, flagField, bodyField, fluid, *storage, *block, pdfFieldPreCol ) );
 
    const auto ubb = flagField->getFlag( UBB_Flag );
    const auto outlet = flagField->getFlag( Outlet_Flag );
diff --git a/apps/benchmarks/NonUniformGrid/NonUniformGrid.cpp b/apps/benchmarks/NonUniformGrid/NonUniformGrid.cpp
index 3a50ad87572dae545de03b87c1fded47f6c0d4b6..f4528cbf65db1514641c07921a1e71d3fcbacc91 100644
--- a/apps/benchmarks/NonUniformGrid/NonUniformGrid.cpp
+++ b/apps/benchmarks/NonUniformGrid/NonUniformGrid.cpp
@@ -434,9 +434,7 @@ struct MyBoundaryTypes
    typedef lbm::NoSlip< LatticeModel_T, flag_t >    NoSlip_T;
    typedef lbm::SimpleUBB< LatticeModel_T, flag_t > UBB_T;
 
-   typedef boost::tuples::tuple< NoSlip_T, UBB_T > BoundaryConditions_T;
-
-   typedef BoundaryHandling< FlagField_T, typename Types<LatticeModel_T>::Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+   typedef BoundaryHandling< FlagField_T, typename Types<LatticeModel_T>::Stencil_T, NoSlip_T, UBB_T > BoundaryHandling_T;
 };  
 
 template< typename LatticeModel_T >
@@ -447,8 +445,6 @@ public:
    using NoSlip_T = typename MyBoundaryTypes< LatticeModel_T >::NoSlip_T;
    using UBB_T = typename MyBoundaryTypes< LatticeModel_T >::UBB_T;
 
-   using BoundaryConditions_T = typename MyBoundaryTypes< LatticeModel_T >::BoundaryConditions_T;
-
    using BoundaryHandling_T = typename MyBoundaryTypes< LatticeModel_T >::BoundaryHandling_T;
 
 
@@ -484,8 +480,8 @@ MyBoundaryHandling<LatticeModel_T>::initialize( IBlock * const block )
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
-                                       UBB_T( "velocity bounce back", UBB_Flag, pdfField, velocity_, real_c(0), real_c(0) ) ) );
+                                                           NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
+                                                           UBB_T( "velocity bounce back", UBB_Flag, pdfField, velocity_, real_c(0), real_c(0) ) );
 
    auto forest = forest_.lock();
    WALBERLA_CHECK_NOT_NULLPTR( forest );
diff --git a/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp b/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp
index 15dbd17058095c4aba2b7e680c51ea46dc2d1d31..ff9d0f46beaaeb6c567244fa0e738ca01d4c72ad 100644
--- a/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp
+++ b/apps/benchmarks/PoiseuilleChannel/PoiseuilleChannel.cpp
@@ -458,9 +458,7 @@ public:
    typedef lbm::NoSlip< LatticeModel_T, flag_t >      NoSlip_T;
    typedef lbm::Curved< LatticeModel_T, FlagField_T > Curved_T;
 
-   typedef boost::tuples::tuple< NoSlip_T, Curved_T > BoundaryConditions_T;
-
-   typedef BoundaryHandling< FlagField_T, typename Types<LatticeModel_T>::Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+   typedef BoundaryHandling< FlagField_T, typename Types<LatticeModel_T>::Stencil_T, NoSlip_T, Curved_T > BoundaryHandling_T;
 
 
 
@@ -488,8 +486,8 @@ MyBoundaryHandling<LatticeModel_T>::operator()( IBlock * const block ) const
    const flag_t fluid = flagField->registerFlag( Fluid_Flag );
 
    return new BoundaryHandling_T( "boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
-                                    Curved_T( "curved", Curved_Flag, pdfField, flagField, fluid ) ) );
+                                   NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
+                                   Curved_T( "curved", Curved_Flag, pdfField, flagField, fluid ) );
 }
 
 
diff --git a/apps/benchmarks/SchaeferTurek/SchaeferTurek.cpp b/apps/benchmarks/SchaeferTurek/SchaeferTurek.cpp
index 46bd0c3663a74e75114bbd94635272fca517c877..2f7a5f6b5858cf51664abc1d0938abf590aab18d 100644
--- a/apps/benchmarks/SchaeferTurek/SchaeferTurek.cpp
+++ b/apps/benchmarks/SchaeferTurek/SchaeferTurek.cpp
@@ -784,11 +784,8 @@ struct MyBoundaryTypes
    typedef lbm::Outlet< LatticeModel_T, FlagField_T, 2, 1 >                       Outlet21_T;
    typedef lbm::Outlet< LatticeModel_T, FlagField_T, 4, 3 >                       Outlet43_T;
    typedef lbm::SimplePressure< LatticeModel_T, flag_t >                          PressureOutlet_T;
-   
-
-   typedef boost::tuples::tuple< NoSlip_T, Obstacle_T, Curved_T, DynamicUBB_T, Outlet21_T, Outlet43_T, PressureOutlet_T > BoundaryConditions_T;
 
-   typedef BoundaryHandling< FlagField_T, typename Types<LatticeModel_T>::Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+   typedef BoundaryHandling< FlagField_T, typename Types<LatticeModel_T>::Stencil_T, NoSlip_T, Obstacle_T, Curved_T, DynamicUBB_T, Outlet21_T, Outlet43_T, PressureOutlet_T > BoundaryHandling_T;
 };
 
 
@@ -805,8 +802,6 @@ public:
    using Outlet43_T = typename MyBoundaryTypes<LatticeModel_T>::Outlet43_T;
    using PressureOutlet_T = typename MyBoundaryTypes<LatticeModel_T>::PressureOutlet_T;
 
-   using BoundaryConditions_T = typename MyBoundaryTypes<LatticeModel_T>::BoundaryConditions_T;
-
    using BoundaryHandling_T = typename MyBoundaryTypes<LatticeModel_T>::BoundaryHandling_T;
 
 
@@ -848,13 +843,13 @@ MyBoundaryHandling<LatticeModel_T>::initialize( IBlock * const block )
    SinusInflowVelocity<Is2D< LatticeModel_T >::value> velocity( setup_.inflowVelocity_L, setup_.raisingTime_L, setup_.sinPeriod_L, setup_.H );
 
    return new BoundaryHandling_T( "boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
+                                    NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
                                   Obstacle_T( "obstacle (staircase)", Obstacle_Flag, pdfField ),
                                     Curved_T( "obstacle (curved)", Curved_Flag, pdfField, flagField, fluid ),
                                 DynamicUBB_T( "velocity bounce back", UBB_Flag, pdfField, timeTracker_, blocks->getLevel(*block), velocity, block->getAABB() ),
                                   Outlet21_T( "outlet (2/1)", Outlet21_Flag, pdfField, flagField, fluid ),
                                   Outlet43_T( "outlet (4/3)", Outlet43_Flag, pdfField, flagField, fluid ),
-                            PressureOutlet_T( "pressure outlet", PressureOutlet_Flag, pdfField, real_t(1) ) ) );
+                            PressureOutlet_T( "pressure outlet", PressureOutlet_Flag, pdfField, real_t(1) ) );
 }
 
 
diff --git a/apps/benchmarks/UniformGrid/UniformGrid.cpp b/apps/benchmarks/UniformGrid/UniformGrid.cpp
index f237a4b8f51e342ef205e873393c109b9067b6db..96e50eb129069daa7b5006ca9fe144eb005bdf7a 100644
--- a/apps/benchmarks/UniformGrid/UniformGrid.cpp
+++ b/apps/benchmarks/UniformGrid/UniformGrid.cpp
@@ -353,9 +353,7 @@ public:
    typedef lbm::NoSlip< LatticeModel_T, flag_t >    NoSlip_T;
    typedef lbm::SimpleUBB< LatticeModel_T, flag_t > UBB_T;
 
-   typedef boost::tuples::tuple< NoSlip_T, UBB_T > BoundaryConditions_T;
-
-   typedef BoundaryHandling< FlagField_T, typename Types<LatticeModel_T>::Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+   typedef BoundaryHandling< FlagField_T, typename Types<LatticeModel_T>::Stencil_T, NoSlip_T, UBB_T > BoundaryHandling_T;
 
 
 
@@ -388,8 +386,8 @@ MyBoundaryHandling<LatticeModel_T>::operator()( IBlock * const block, const Stru
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
-                                       UBB_T( "velocity bounce back", UBB_Flag, pdfField, velocity_, real_c(0), real_c(0) ) ) );
+                                                           NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
+                                                           UBB_T( "velocity bounce back", UBB_Flag, pdfField, velocity_, real_c(0), real_c(0) ) );
 
    CellInterval domainBB = storage->getDomainCellBB();
    storage->transformGlobalToBlockLocalCellInterval( domainBB, *block );
diff --git a/apps/showcases/BidisperseFluidizedBed/BidisperseFluidizedBedDPM.cpp b/apps/showcases/BidisperseFluidizedBed/BidisperseFluidizedBedDPM.cpp
index 34cb30ae1d1530649f93ab9fd81a5c62cfb60ba5..b91418447e3a8a5b0e7bcf503aa1bfbe98c4f580 100644
--- a/apps/showcases/BidisperseFluidizedBed/BidisperseFluidizedBedDPM.cpp
+++ b/apps/showcases/BidisperseFluidizedBed/BidisperseFluidizedBedDPM.cpp
@@ -98,8 +98,7 @@ typedef lbm::Outlet< LatticeModel_T, FlagField_T >                     Outflow_T
 typedef lbm::SimplePressure< LatticeModel_T, flag_t >                  Outflow_T;
 #endif
 
-typedef boost::tuples::tuple< NoSlip_T, Inflow_T, Outflow_T >          BoundaryConditions_T;
-typedef BoundaryHandling<FlagField_T, Stencil_T, BoundaryConditions_T> BoundaryHandling_T;
+typedef BoundaryHandling<FlagField_T, Stencil_T, NoSlip_T, Inflow_T, Outflow_T> BoundaryHandling_T;
 
 typedef std::tuple<pe::Plane, pe::Sphere> BodyTypeTuple ;
 
@@ -165,14 +164,14 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
 
 #ifdef OutletBC
    BoundaryHandling_T * handling = new BoundaryHandling_T( "Boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+                                    NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
                                     Inflow_T( "Inflow", Inflow_Flag, pdfField, Vector3<real_t>(real_t(0),real_t(0),uInflow_) ),
-                                    Outflow_T( "Outflow", Outflow_Flag, pdfField, flagField, fluid ) ) );
+                                    Outflow_T( "Outflow", Outflow_Flag, pdfField, flagField, fluid ) );
 #else
    BoundaryHandling_T * handling = new BoundaryHandling_T( "Boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+                                    NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
                                     Inflow_T( "Inflow", Inflow_Flag, pdfField, Vector3<real_t>(real_t(0),real_t(0),uInflow_) ),
-                                    Outflow_T( "Outflow", Outflow_Flag, pdfField, real_t(1) ) ) );
+                                    Outflow_T( "Outflow", Outflow_Flag, pdfField, real_t(1) ) );
 #endif
 
    const auto noslip  = flagField->getFlag( NoSlip_Flag );
diff --git a/apps/tutorials/lbm/03_LBLidDrivenCavity.cpp b/apps/tutorials/lbm/03_LBLidDrivenCavity.cpp
index e28259f8d46d35eb081635ec455fa71be16b8caf..7aa2e29bc21be81d2d2877dcea2abb3af21b5662 100644
--- a/apps/tutorials/lbm/03_LBLidDrivenCavity.cpp
+++ b/apps/tutorials/lbm/03_LBLidDrivenCavity.cpp
@@ -79,10 +79,7 @@ typedef lbm::NoSlip< LatticeModel_T, flag_t >     NoSlip_T; // no slip boundary
 typedef lbm::SimpleUBB< LatticeModel_T, flag_t >  UBB_T;    // velocity bounce back boundary condition that internally works with one ...
                                                             // ... constant velocity that must be set during the setup phase
 
-typedef boost::tuples::tuple< NoSlip_T, UBB_T >  BoundaryConditions_T; // a collection of all boundary conditions - required for
-                                                                       // defining the type of the boundary handling (see below)
-
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T; // the boundary handling
+typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, UBB_T > BoundaryHandling_T; // the boundary handling, includes a collection of all boundary conditions
 
 ///////////
 // FLAGS //
@@ -139,11 +136,11 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block ) cons
    // The fluid flag is used in order to decide whether a certain cell is part of the domain (and therefore has to be
    // treated by the boundary handling if it is located next to a boundary). Also, the no slip and velocity bounce back
    // boundary conditions are initialized (Please note that the order in which these boundary conditions are initialized
-   // must be identical to the order in 'BoundaryConditions_T'!).
+   // must be identical to the order of the template arguments of 'BoundaryHandling_T'!).
 
    return new BoundaryHandling_T( "boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
-                                       UBB_T( "velocity bounce back", UBB_Flag, pdfField, topVelocity_, real_c(0), real_c(0) ) ) );
+                                  NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
+                                     UBB_T( "velocity bounce back", UBB_Flag, pdfField, topVelocity_, real_c(0), real_c(0) ) );
 }
 
 
@@ -420,4 +417,4 @@ int main( int argc, char ** argv )
 int main( int argc, char ** argv )
 {
    return walberla::main(argc, argv);
-}
\ No newline at end of file
+}
diff --git a/src/boundary/BoundaryHandling.h b/src/boundary/BoundaryHandling.h
index 981cc3e4fc34f960951631c0202e428d6a5d3269..e61e7f90b98c51524ba8ac9ce1058dee51d45672 100644
--- a/src/boundary/BoundaryHandling.h
+++ b/src/boundary/BoundaryHandling.h
@@ -56,12 +56,12 @@ typedef UID< BHUIDGenerator > BoundaryHandlingUID;
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple > // Tuple: all the boundary classes that are considered by this boundary handler
+template< typename FlagField_T, typename Stencil, typename... Boundaries > // Boundaries: all the boundary classes that are considered by this boundary handler
 class BoundaryHandling
 {
 public:
 
-   template< typename F, typename T >
+   template< typename F, typename... T >
    friend class BoundaryHandlingCollection;
 
    typedef FlagField_T                               FlagField;
@@ -87,7 +87,7 @@ public:
 
 
 
-   BoundaryHandling( const std::string & identifier, FlagField_T * const flagField, const flag_t domain, const Tuple & boundaryConditions,
+   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 "==".
@@ -649,6 +649,7 @@ private:
    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 boost::tuples::tuple<Boundaries...> Tuple;
    Tuple boundaryConditions_;
    bool  threadSafeBCs_;
 
@@ -656,9 +657,9 @@ private:
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-BoundaryHandling< FlagField_T, Stencil, Tuple >::BoundaryHandling( const std::string & identifier, FlagField_T * const flagField,
-                                                                   const flag_t domain, const Tuple & boundaryConditions, const Mode mode ) :
+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 ),
@@ -673,7 +674,7 @@ BoundaryHandling< FlagField_T, Stencil, Tuple >::BoundaryHandling( const std::st
    domain_( domain ),
    mode_( mode ),
    dirty_( false ),
-   boundaryConditions_( boundaryConditions ),
+   boundaryConditions_( boost::make_tuple(boundaryConditions...) ),
    threadSafeBCs_( true )
 {
    setupBoundaryConditions( boundaryConditions_ );
@@ -694,40 +695,40 @@ BoundaryHandling< FlagField_T, Stencil, Tuple >::BoundaryHandling( const std::st
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isEmpty( 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_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 Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::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... >::isNearBoundary( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
 {
    return flagField_->isFlagSet( x, y, z, nearBoundary_ );
 }
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isBoundary( 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 Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::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... >::isDomain( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
 {
    return flagField_->isPartOfMaskSet( x, y, z, domain_ );
 }
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isEmpty( const Cell & cell ) 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_ );
@@ -735,32 +736,32 @@ inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isEmpty( const Cell
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isNearBoundary( const Cell & cell ) const
+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 Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isBoundary( const Cell & cell ) const
+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 Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isDomain( const Cell & cell ) const
+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 Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isEmpty( const ConstFlagFieldBaseIterator & it ) const
+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_ );
 
@@ -769,8 +770,8 @@ inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isEmpty( const Cons
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isNearBoundary( const ConstFlagFieldBaseIterator & it ) const
+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_ );
 
@@ -779,8 +780,8 @@ inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isNearBoundary( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isBoundary( const ConstFlagFieldBaseIterator & it ) const
+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_ );
 
@@ -789,8 +790,8 @@ inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isBoundary( const C
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isDomain( const ConstFlagFieldBaseIterator & it ) const
+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_ );
 
@@ -799,16 +800,16 @@ inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::isDomain( const Con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::containsBoundaryCondition( const BoundaryUID & uid ) const
+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 Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::containsBoundaryCondition( const FlagUID & flag ) const
+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;
@@ -818,26 +819,26 @@ inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::containsBoundaryCon
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Boundary_T >
-inline const Boundary_T & BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryCondition( const BoundaryUID & uid ) const
+inline const Boundary_T & BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getBoundaryCondition( const BoundaryUID & uid ) const
 {
    return getBoundaryCondition< Boundary_T, typename Tuple::head_type, typename Tuple::tail_type >( uid, boundaryConditions_ );
 }
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Boundary_T >
-inline Boundary_T & BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryCondition( const BoundaryUID & uid )
+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 Tuple >
-inline BoundaryUID BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryUID( const FlagUID & flag ) const
+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 ) );
 
@@ -846,8 +847,8 @@ inline BoundaryUID BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryU
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline BoundaryUID BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryUID( const flag_t flag ) const
+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 ) );
@@ -857,16 +858,16 @@ inline BoundaryUID BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryU
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline uint_t BoundaryHandling< FlagField_T, Stencil, Tuple >::numberOfMatchingBoundaryConditions( const flag_t mask ) const
+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 Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::checkConsistency( const uint_t numberOfGhostLayersToInclude ) const
+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 );
@@ -874,8 +875,8 @@ inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::checkConsistency( c
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-bool BoundaryHandling< FlagField_T, Stencil, Tuple >::checkConsistency( const CellInterval & cells ) const
+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 );
@@ -937,8 +938,8 @@ bool BoundaryHandling< FlagField_T, Stencil, Tuple >::checkConsistency( const Ce
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::refresh( const uint_t numberOfGhostLayersToInclude )
+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 );
@@ -946,8 +947,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::refresh( const uint
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-void BoundaryHandling< FlagField_T, Stencil, Tuple >::refresh( const CellInterval & 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 );
@@ -974,8 +975,8 @@ void BoundaryHandling< FlagField_T, Stencil, Tuple >::refresh( const CellInterva
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-void BoundaryHandling< FlagField_T, Stencil, Tuple >::refreshOutermostLayer( cell_idx_t thickness )
+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) );
@@ -1015,8 +1016,8 @@ void BoundaryHandling< FlagField_T, Stencil, Tuple >::refreshOutermostLayer( cel
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setDomain( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setDomain( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT( field::isFlag( domain_ ) );
    setDomain( domain_, x, y, z );
@@ -1024,8 +1025,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setDomain( const ce
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setDomain( const flag_t domainSubFlag,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setDomain( const flag_t domainSubFlag,
                                                                         const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT_EQUAL( domain_ & domainSubFlag, domainSubFlag );
@@ -1037,8 +1038,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setDomain( const fl
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setDomain( const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setDomain( const CellInterval & cells )
 {
    WALBERLA_ASSERT( field::isFlag( domain_ ) );
    setDomain( domain_, cells );
@@ -1046,8 +1047,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setDomain( const Ce
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setDomain( const flag_t domainSubFlag, const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setDomain( const flag_t domainSubFlag, const CellInterval & cells )
 {
    WALBERLA_ASSERT_EQUAL( domain_ & domainSubFlag, domainSubFlag );
    WALBERLA_ASSERT( field::isFlag( domainSubFlag ) );
@@ -1090,9 +1091,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setDomain( const fl
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setDomain( const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setDomain( const CellIterator & begin, const CellIterator & end )
 {
    WALBERLA_ASSERT( field::isFlag( domain_ ) );
    setDomain( domain_, begin, end );
@@ -1100,9 +1101,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setDomain( const Ce
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setDomain( const flag_t domainSubFlag, const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setDomain( const flag_t domainSubFlag, const CellIterator & begin, const CellIterator & end )
 {
    WALBERLA_ASSERT_EQUAL( domain_ & domainSubFlag, domainSubFlag );
    WALBERLA_ASSERT( field::isFlag( domainSubFlag ) );
@@ -1120,16 +1121,16 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setDomain( const fl
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceDomain( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceDomain( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    forceDomain( domain_, x, y, z );
 }
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceDomain( const flag_t domainSubFlag,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceDomain( const flag_t domainSubFlag,
                                                                           const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    clear( x, y, z );
@@ -1138,16 +1139,16 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceDomain( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceDomain( const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceDomain( const CellInterval & cells )
 {
    forceDomain( domain_, cells );
 }
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceDomain( const flag_t domainSubFlag, const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceDomain( const flag_t domainSubFlag, const CellInterval & cells )
 {
    clear( cells );
    setDomain( domainSubFlag, cells );
@@ -1155,18 +1156,18 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceDomain( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceDomain( const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceDomain( const CellIterator & begin, const CellIterator & end )
 {
    forceDomain( domain_, begin, end );
 }
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceDomain( const flag_t domainSubFlag, const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceDomain( const flag_t domainSubFlag, const CellIterator & begin, const CellIterator & end )
 {
    clear( begin, end );
    setDomain( domainSubFlag, begin, end );
@@ -1174,8 +1175,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceDomain( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( const uint_t numberOfGhostLayersToInclude )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::fillWithDomain( const uint_t numberOfGhostLayersToInclude )
 {
    WALBERLA_ASSERT_GREATER_EQUAL( flagField_->nrOfGhostLayers(), numberOfGhostLayersToInclude );
    WALBERLA_ASSERT( field::isFlag( domain_ ) );
@@ -1184,8 +1185,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( const flag_t domainSubFlag, const uint_t numberOfGhostLayersToInclude )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::fillWithDomain( const flag_t domainSubFlag, const uint_t numberOfGhostLayersToInclude )
 {
    WALBERLA_ASSERT_EQUAL( domain_ & domainSubFlag, domainSubFlag );
    WALBERLA_ASSERT( field::isFlag( domainSubFlag ) );
@@ -1197,8 +1198,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::fillWithDomain( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT( field::isFlag( domain_ ) );
    fillWithDomain( domain_, x, y, z );
@@ -1206,8 +1207,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( const flag_t domainSubFlag,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::fillWithDomain( const flag_t domainSubFlag,
                                                                              const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT_EQUAL( domain_ & domainSubFlag, domainSubFlag );
@@ -1219,8 +1220,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::fillWithDomain( const CellInterval & cells )
 {
    WALBERLA_ASSERT( field::isFlag( domain_ ) );
    fillWithDomain( domain_, cells );
@@ -1228,8 +1229,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( const flag_t domainSubFlag, const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::fillWithDomain( const flag_t domainSubFlag, const CellInterval & cells )
 {
    WALBERLA_ASSERT_EQUAL( domain_ & domainSubFlag, domainSubFlag );
    WALBERLA_ASSERT( field::isFlag( domainSubFlag ) );
@@ -1253,18 +1254,18 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::fillWithDomain( const CellIterator & begin, const CellIterator & end )
 {
    fillWithDomain( domain_, begin, end );
 }
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( const flag_t domainSubFlag,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::fillWithDomain( const flag_t domainSubFlag,
                                                                              const CellIterator & begin, const CellIterator & end )
 {
    WALBERLA_ASSERT_EQUAL( domain_ & domainSubFlag, domainSubFlag );
@@ -1283,8 +1284,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::fillWithDomain( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline shared_ptr<BoundaryConfiguration> BoundaryHandling< FlagField_T, Stencil, Tuple >::createBoundaryConfiguration( const BoundaryUID & uid,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline shared_ptr<BoundaryConfiguration> BoundaryHandling< FlagField_T, Stencil, Boundaries... >::createBoundaryConfiguration( const BoundaryUID & uid,
                                                                                                                        const Config::BlockHandle & config ) const
 {
    return createBoundaryConfiguration( boundaryConditions_, uid, config );
@@ -1292,8 +1293,8 @@ inline shared_ptr<BoundaryConfiguration> BoundaryHandling< FlagField_T, Stencil,
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const FlagUID & flag,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setBoundary( const FlagUID & flag,
                                                                           const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                           const BoundaryConfiguration & parameter )
 {
@@ -1304,8 +1305,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const flag_t flag,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setBoundary( const flag_t flag,
                                                                           const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                           const BoundaryConfiguration & parameter )
 {
@@ -1319,8 +1320,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const FlagUID & flag, const CellInterval & cells,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setBoundary( const FlagUID & flag, const CellInterval & cells,
                                                                           const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -1330,8 +1331,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const flag_t flag, const CellInterval & cells,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setBoundary( const flag_t flag, const CellInterval & cells,
                                                                           const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT_EQUAL( flag & boundary_, flag );
@@ -1349,9 +1350,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setBoundary( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
                                                                           const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -1361,9 +1362,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const flag_t flag, const CellIterator & begin, const CellIterator & end,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setBoundary( const flag_t flag, const CellIterator & begin, const CellIterator & end,
                                                                           const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT_EQUAL( flag & boundary_, flag );
@@ -1379,8 +1380,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceBoundary( const FlagUID & flag,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceBoundary( const FlagUID & flag,
                                                                             const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                             const BoundaryConfiguration & parameter )
 {
@@ -1391,8 +1392,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceBoundary( cons
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceBoundary( const flag_t flag,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceBoundary( const flag_t flag,
                                                                             const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                             const BoundaryConfiguration & parameter )
 {
@@ -1403,8 +1404,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceBoundary( cons
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceBoundary( const FlagUID & flag, const CellInterval & cells,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceBoundary( const FlagUID & flag, const CellInterval & cells,
                                                                             const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -1414,8 +1415,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceBoundary( cons
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceBoundary( const flag_t flag, const CellInterval & cells,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceBoundary( const flag_t flag, const CellInterval & cells,
                                                                             const BoundaryConfiguration & parameter )
 {
    clear( cells );
@@ -1424,9 +1425,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceBoundary( cons
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceBoundary( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceBoundary( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
                                                                             const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -1436,9 +1437,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceBoundary( cons
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceBoundary( const flag_t flag, const CellIterator & begin, const CellIterator & end,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceBoundary( const flag_t flag, const CellIterator & begin, const CellIterator & end,
                                                                             const BoundaryConfiguration & parameter )
 {
    for( auto cell = begin; cell != end; ++cell )
@@ -1450,8 +1451,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceBoundary( cons
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeDomain( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    if( outerBB_.contains(x,y,z) )
    {
@@ -1463,8 +1464,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const flag_t mask,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeDomain( const flag_t mask,
                                                                            const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    if( outerBB_.contains(x,y,z) )
@@ -1480,8 +1481,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const uint_t numberOfGhostLayersToInclude )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeDomain( const uint_t numberOfGhostLayersToInclude )
 {
    WALBERLA_ASSERT_GREATER_EQUAL( flagField_->nrOfGhostLayers(), numberOfGhostLayersToInclude );
 
@@ -1491,8 +1492,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeDomain( const CellInterval & cells )
 {
    CellInterval localCells( outerBB_ );
    localCells.intersect( cells );
@@ -1511,8 +1512,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const flag_t mask, const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeDomain( const flag_t mask, const CellInterval & cells )
 {
    const flag_t dMask = mask & domain_;
    if( dMask == 0 )
@@ -1537,9 +1538,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeDomain( const CellIterator & begin, const CellIterator & end )
 {
    for( auto cell = begin; cell != end; ++cell )
       removeDomain( cell->x(), cell->y(), cell->z() );
@@ -1547,9 +1548,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const flag_t mask, const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeDomain( const flag_t mask, const CellIterator & begin, const CellIterator & end )
 {
    const flag_t dMask = mask & domain_;
    if( dMask == 0 )
@@ -1575,8 +1576,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeDomain( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeBoundary( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    if( outerBB_.contains(x,y,z) && flagField_->isPartOfMaskSet( x, y, z, boundary_ ) )
       removeBoundary( boundaryConditions_, x, y, z );
@@ -1584,8 +1585,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( const FlagUID & flag,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeBoundary( const FlagUID & flag,
                                                                              const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    if( !flagField_->flagExists( flag ) )
@@ -1596,8 +1597,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( const flag_t mask,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeBoundary( const flag_t mask,
                                                                              const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    const flag_t bMask = mask & boundary_;
@@ -1610,8 +1611,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( const uint_t numberOfGhostLayersToInclude )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeBoundary( const uint_t numberOfGhostLayersToInclude )
 {
    WALBERLA_ASSERT_GREATER_EQUAL( flagField_->nrOfGhostLayers(), numberOfGhostLayersToInclude );
 
@@ -1621,8 +1622,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeBoundary( const CellInterval & cells )
 {
    CellInterval localCells( outerBB_ );
    localCells.intersect( cells );
@@ -1675,8 +1676,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( const FlagUID & flag, const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeBoundary( const FlagUID & flag, const CellInterval & cells )
 {
    if( !flagField_->flagExists( flag ) )
       return;
@@ -1686,8 +1687,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( const flag_t mask, const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeBoundary( const flag_t mask, const CellInterval & cells )
 {
    const flag_t bMask = mask & boundary_;
    if( bMask == 0 )
@@ -1703,18 +1704,18 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeBoundary( const CellIterator & begin, const CellIterator & end )
 {
    removeBoundary( boundary_, begin, end );
 }
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( const FlagUID & flag, const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeBoundary( const FlagUID & flag, const CellIterator & begin, const CellIterator & end )
 {
    if( !flagField_->flagExists( flag ) )
       return;
@@ -1724,9 +1725,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( const flag_t mask, const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeBoundary( const flag_t mask, const CellIterator & begin, const CellIterator & end )
 {
    const flag_t bMask = mask & boundary_;
    if( bMask == 0 )
@@ -1745,8 +1746,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                       const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -1756,8 +1757,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setFlag( const Flag
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setFlag( const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setFlag( const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                       const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( field::isFlag(flag) );
@@ -1775,8 +1776,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setFlag( const flag
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setFlag( const FlagUID & flag, const CellInterval & cells,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setFlag( const FlagUID & flag, const CellInterval & cells,
                                                                       const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -1786,8 +1787,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setFlag( const Flag
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setFlag( const flag_t flag, const CellInterval & cells,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setFlag( const flag_t flag, const CellInterval & cells,
                                                                       const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( field::isFlag(flag) );
@@ -1812,9 +1813,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setFlag( const flag
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
                                                                       const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -1824,9 +1825,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setFlag( const Flag
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setFlag( const flag_t flag, const CellIterator & begin, const CellIterator & end,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setFlag( const flag_t flag, const CellIterator & begin, const CellIterator & end,
                                                                       const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( field::isFlag(flag) );
@@ -1848,8 +1849,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setFlag( const flag
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                         const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -1859,8 +1860,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceFlag( const Fl
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceFlag( const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceFlag( const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                         const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( field::isFlag(flag) );
@@ -1875,8 +1876,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceFlag( const fl
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceFlag( const FlagUID & flag, const CellInterval & cells,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceFlag( const FlagUID & flag, const CellInterval & cells,
                                                                         const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -1886,8 +1887,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceFlag( const Fl
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceFlag( const flag_t flag, const CellInterval & cells,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceFlag( const flag_t flag, const CellInterval & cells,
                                                                         const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( field::isFlag(flag) );
@@ -1909,9 +1910,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceFlag( const fl
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
                                                                         const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -1921,9 +1922,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceFlag( const Fl
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceFlag( const flag_t flag, const CellIterator & begin, const CellIterator & end,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::forceFlag( const flag_t flag, const CellIterator & begin, const CellIterator & end,
                                                                         const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( field::isFlag(flag) );
@@ -1942,8 +1943,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::forceFlag( const fl
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const FlagUID & flag, const uint_t numberOfGhostLayersToInclude )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeFlag( const FlagUID & flag, const uint_t numberOfGhostLayersToInclude )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
 
@@ -1952,8 +1953,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const F
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const flag_t flag, const uint_t numberOfGhostLayersToInclude )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeFlag( const flag_t flag, const uint_t numberOfGhostLayersToInclude )
 {
    WALBERLA_ASSERT_GREATER_EQUAL( flagField_->nrOfGhostLayers(), numberOfGhostLayersToInclude );
 
@@ -1963,8 +1964,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const f
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
 
@@ -1973,8 +1974,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const F
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeFlag( const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    if( ( flag & boundary_ ) == flag )
       removeBoundary( flag, x, y, z );
@@ -1986,8 +1987,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const f
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const FlagUID & flag, const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeFlag( const FlagUID & flag, const CellInterval & cells )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
 
@@ -1996,8 +1997,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const F
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const flag_t flag, const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeFlag( const flag_t flag, const CellInterval & cells )
 {
    if( ( flag & boundary_ ) == flag )
       removeBoundary( flag, cells );
@@ -2016,9 +2017,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const f
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
 
@@ -2027,9 +2028,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const F
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const flag_t flag, const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeFlag( const flag_t flag, const CellIterator & begin, const CellIterator & end )
 {
    if( ( flag & boundary_ ) == flag )
       removeBoundary( flag, begin, end );
@@ -2045,8 +2046,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeFlag( const f
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::clear( const uint_t numberOfGhostLayersToInclude )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::clear( const uint_t numberOfGhostLayersToInclude )
 {
    WALBERLA_ASSERT_GREATER_EQUAL( flagField_->nrOfGhostLayers(), numberOfGhostLayersToInclude );
 
@@ -2056,8 +2057,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::clear( const uint_t
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::clear( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::clear( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    if( outerBB_.contains(x,y,z) )
    {
@@ -2073,8 +2074,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::clear( const cell_i
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::clear( const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::clear( const CellInterval & cells )
 {
    CellInterval localCells( outerBB_ );
    localCells.intersect( cells );
@@ -2112,9 +2113,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::clear( const CellIn
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::clear( const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::clear( const CellIterator & begin, const CellIterator & end )
 {
    for( auto cell = begin; cell != end; ++cell )
       clear( cell->x(), cell->y(), cell->z() );
@@ -2122,8 +2123,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::clear( const CellIt
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::operator()( const uint_t numberOfGhostLayersToInclude )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::operator()( const uint_t numberOfGhostLayersToInclude )
 {
    WALBERLA_ASSERT_LESS( numberOfGhostLayersToInclude, flagField_->nrOfGhostLayers() );
    WALBERLA_ASSERT( checkConsistency( numberOfGhostLayersToInclude ) );
@@ -2182,8 +2183,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::operator()( const u
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::operator()( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::operator()( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT( innerBB_.contains(x,y,z) );
 
@@ -2203,8 +2204,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::operator()( const c
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::operator()( const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::operator()( const CellInterval & cells )
 {
    CellInterval localCells( innerBB_ );
    localCells.intersect( cells );
@@ -2234,9 +2235,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::operator()( const C
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename CellIterator >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::operator()( const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::operator()( const CellIterator & begin, const CellIterator & end )
 {
    for( auto cell = begin; cell != end; ++cell )
    {
@@ -2251,9 +2252,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::operator()( const C
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Buffer_T >
-void BoundaryHandling< FlagField_T, Stencil, Tuple >::pack( Buffer_T & buffer, const CellInterval & interval, const bool assumeIdenticalFlagMapping ) const
+void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::pack( Buffer_T & buffer, const CellInterval & interval, const bool assumeIdenticalFlagMapping ) const
 {
    WALBERLA_UNUSED( assumeIdenticalFlagMapping );
 
@@ -2283,9 +2284,9 @@ void BoundaryHandling< FlagField_T, Stencil, Tuple >::pack( Buffer_T & buffer, c
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Buffer_T >
-void BoundaryHandling< FlagField_T, Stencil, Tuple >::unpack( Buffer_T & buffer, const CellInterval & interval, const bool assumeIdenticalFlagMapping )
+void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::unpack( Buffer_T & buffer, const CellInterval & interval, const bool assumeIdenticalFlagMapping )
 {
    bool identicalFlagMapping = false;
    std::vector< flag_t > flagMapping = getNeighborFlagMapping( buffer, assumeIdenticalFlagMapping, identicalFlagMapping ); // neighbor-flag_t -> flag_t
@@ -2322,9 +2323,9 @@ void BoundaryHandling< FlagField_T, Stencil, Tuple >::unpack( Buffer_T & buffer,
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Buffer_T >
-void BoundaryHandling< FlagField_T, Stencil, Tuple >::pack( Buffer_T & buffer, stencil::Direction direction, const uint_t numberOfLayers,
+void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::pack( Buffer_T & buffer, stencil::Direction direction, const uint_t numberOfLayers,
                                                             const bool assumeIdenticalFlagMapping ) const
 {
    CellInterval interval = getPackingInterval( direction, numberOfLayers );
@@ -2333,9 +2334,9 @@ void BoundaryHandling< FlagField_T, Stencil, Tuple >::pack( Buffer_T & buffer, s
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Buffer_T >
-void BoundaryHandling< FlagField_T, Stencil, Tuple >::unpack( Buffer_T & buffer, stencil::Direction direction, const uint_t numberOfLayers,
+void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::unpack( Buffer_T & buffer, stencil::Direction direction, const uint_t numberOfLayers,
                                                               const bool assumeIdenticalFlagMapping )
 {
    CellInterval interval = getUnpackingInterval( direction, numberOfLayers );
@@ -2344,8 +2345,8 @@ void BoundaryHandling< FlagField_T, Stencil, Tuple >::unpack( Buffer_T & buffer,
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::toStream( std::ostream & os ) const {
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::toStream( std::ostream & os ) const {
 
    os << "==================== BoundaryHandling ====================\n\n"
       << "Identifier: " << uid_.getIdentifier() << "\n\n"
@@ -2370,8 +2371,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::toStream( std::ostr
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline std::string BoundaryHandling< FlagField_T, Stencil, Tuple >::toString() const {
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline std::string BoundaryHandling< FlagField_T, Stencil, Boundaries... >::toString() const {
 
    std::ostringstream oss;
    toStream( oss );
@@ -2380,8 +2381,8 @@ inline std::string BoundaryHandling< FlagField_T, Stencil, Tuple >::toString() c
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-CellInterval BoundaryHandling< FlagField_T, Stencil, Tuple >::getGhostLayerCellInterval( const uint_t numberOfGhostLayersToInclude ) const
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+CellInterval BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getGhostLayerCellInterval( const uint_t numberOfGhostLayersToInclude ) const
 {
    CellInterval cells( -cell_idx_c( numberOfGhostLayersToInclude ),
                        -cell_idx_c( numberOfGhostLayersToInclude ),
@@ -2394,9 +2395,9 @@ CellInterval BoundaryHandling< FlagField_T, Stencil, Tuple >::getGhostLayerCellI
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-void BoundaryHandling< FlagField_T, Stencil, Tuple >::setupBoundaryConditions( boost::tuples::cons<Head, Tail> & boundaryConditions )
+void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setupBoundaryConditions( boost::tuples::cons<Head, Tail> & boundaryConditions )
 {
    Head & boundaryCondition = boundaryConditions.get_head();
 
@@ -2430,8 +2431,8 @@ void BoundaryHandling< FlagField_T, Stencil, Tuple >::setupBoundaryConditions( b
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline std::vector< BoundaryUID > BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryUIDs() const
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline std::vector< BoundaryUID > BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getBoundaryUIDs() const
 {
    std::vector< BoundaryUID > uids;
    getBoundaryUIDs( boundaryConditions_, uids );
@@ -2440,9 +2441,9 @@ inline std::vector< BoundaryUID > BoundaryHandling< FlagField_T, Stencil, Tuple
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryUIDs( const boost::tuples::cons<Head, Tail> & boundaryConditions,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getBoundaryUIDs( const boost::tuples::cons<Head, Tail> & boundaryConditions,
                                                                               std::vector< BoundaryUID > & uids ) const
 {
    uids.push_back( boundaryConditions.get_head().getUID() );
@@ -2451,9 +2452,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryUIDs( co
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline BoundaryUID BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryUID( const boost::tuples::cons<Head, Tail> & boundaryConditions,
+inline BoundaryUID BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getBoundaryUID( const boost::tuples::cons<Head, Tail> & boundaryConditions,
                                                                                     const flag_t flag ) const
 {
    const Head & boundaryCondition = boundaryConditions.get_head();
@@ -2470,8 +2471,8 @@ inline BoundaryUID BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryU
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline BoundaryUID BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryUID( const boost::tuples::null_type &,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline BoundaryUID BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getBoundaryUID( const boost::tuples::null_type &,
                                                                                     const flag_t flag ) const
 {
    if( !flagField_->isRegistered( flag ) )
@@ -2489,9 +2490,9 @@ inline BoundaryUID BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryU
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::containsBoundaryCondition( const boost::tuples::cons<Head, Tail> & boundaryConditions,
+inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::containsBoundaryCondition( const boost::tuples::cons<Head, Tail> & boundaryConditions,
                                                                                         const BoundaryUID & uid ) const
 {
    if( boundaryConditions.get_head().getUID() == uid )
@@ -2501,10 +2502,10 @@ inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::containsBoundaryCon
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline typename BoundaryHandling< FlagField_T, Stencil, Tuple >::flag_t
-   BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryMask( const boost::tuples::cons<Head, Tail> & boundaryConditions,
+inline typename BoundaryHandling< FlagField_T, Stencil, Boundaries... >::flag_t
+   BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getBoundaryMask( const boost::tuples::cons<Head, Tail> & boundaryConditions,
                                                                      const BoundaryUID & uid ) const
 {
    const Head & boundaryCondition = boundaryConditions.get_head();
@@ -2516,17 +2517,17 @@ inline typename BoundaryHandling< FlagField_T, Stencil, Tuple >::flag_t
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline uint_t BoundaryHandling< FlagField_T, Stencil, Tuple >::numberOfMatchingBoundaryConditions( const BoundaryUID & uid ) const
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline uint_t BoundaryHandling< FlagField_T, Stencil, Boundaries... >::numberOfMatchingBoundaryConditions( const BoundaryUID & uid ) const
 {
    return numberOfMatchingBoundaryConditions( boundaryConditions_, uid );
 }
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline uint_t BoundaryHandling< FlagField_T, Stencil, Tuple >::numberOfMatchingBoundaryConditions( const boost::tuples::cons<Head, Tail> & boundaryConditions,
+inline uint_t BoundaryHandling< FlagField_T, Stencil, Boundaries... >::numberOfMatchingBoundaryConditions( const boost::tuples::cons<Head, Tail> & boundaryConditions,
                                                                                                  const BoundaryUID & uid ) const
 {
    return ( ( boundaryConditions.get_head().getUID() == uid ) ? uint_c(1) : uint_c(0) ) +
@@ -2535,9 +2536,9 @@ inline uint_t BoundaryHandling< FlagField_T, Stencil, Tuple >::numberOfMatchingB
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline uint_t BoundaryHandling< FlagField_T, Stencil, Tuple >::numberOfMatchingBoundaryConditions( const boost::tuples::cons<Head, Tail> & boundaryConditions,
+inline uint_t BoundaryHandling< FlagField_T, Stencil, Boundaries... >::numberOfMatchingBoundaryConditions( const boost::tuples::cons<Head, Tail> & boundaryConditions,
                                                                                                  const flag_t mask ) const
 {
    return ( ( ( boundaryConditions.get_head().getMask() & mask ) != 0 ) ? uint_c(1) : uint_c(0) ) +
@@ -2546,8 +2547,8 @@ inline uint_t BoundaryHandling< FlagField_T, Stencil, Tuple >::numberOfMatchingB
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::checkFlagField( const uint_t numberOfGhostLayersToInclude ) const
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline bool BoundaryHandling< FlagField_T, Stencil, Boundaries... >::checkFlagField( const uint_t numberOfGhostLayersToInclude ) const
 {
    if( !flagField_->isRegistered( nearBoundary_ ) )
       return false;
@@ -2589,8 +2590,8 @@ inline bool BoundaryHandling< FlagField_T, Stencil, Tuple >::checkFlagField( con
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::addDomain( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const flag_t domain )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::addDomain( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const flag_t domain )
 {
    WALBERLA_ASSERT( outerBB_.contains(x,y,z) );
    WALBERLA_ASSERT_EQUAL( domain_ & domain, domain );
@@ -2617,9 +2618,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::addDomain( const ce
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline shared_ptr<BoundaryConfiguration> BoundaryHandling< FlagField_T, Stencil, Tuple >::createBoundaryConfiguration( const boost::tuples::cons<Head, Tail> & boundaryConditions,
+inline shared_ptr<BoundaryConfiguration> BoundaryHandling< FlagField_T, Stencil, Boundaries... >::createBoundaryConfiguration( const boost::tuples::cons<Head, Tail> & boundaryConditions,
                                                                                                                        const BoundaryUID & uid, const Config::BlockHandle & config ) const
 {
    if( boundaryConditions.get_head().getUID() == uid )
@@ -2628,8 +2629,8 @@ inline shared_ptr<BoundaryConfiguration> BoundaryHandling< FlagField_T, Stencil,
    return createBoundaryConfiguration( boundaryConditions.get_tail(), uid, config );
 }
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline shared_ptr<BoundaryConfiguration> BoundaryHandling< FlagField_T, Stencil, Tuple >::createBoundaryConfiguration( const boost::tuples::null_type &,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline shared_ptr<BoundaryConfiguration> BoundaryHandling< FlagField_T, Stencil, Boundaries... >::createBoundaryConfiguration( const boost::tuples::null_type &,
                                                                                                                        const BoundaryUID & uid, const Config::BlockHandle & ) const
 {
    WALBERLA_CHECK( false, "There is no boundary condition registered at boundary handling " << uid_ << " for a boundary with UID" << uid << "." );
@@ -2639,8 +2640,8 @@ inline shared_ptr<BoundaryConfiguration> BoundaryHandling< FlagField_T, Stencil,
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::addNearBoundary( const CellInterval & cells )
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::addNearBoundary( const CellInterval & cells )
 {
    WALBERLA_ASSERT( innerBB_.contains( cells ) );
 
@@ -2660,8 +2661,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::addNearBoundary( co
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::addBoundary( const flag_t flag,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::addBoundary( const flag_t flag,
                                                                           const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT( outerBB_.contains(x,y,z) );
@@ -2683,9 +2684,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::addBoundary( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( boost::tuples::cons<Head, Tail> & boundaryConditions, const flag_t flag,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setBoundary( boost::tuples::cons<Head, Tail> & boundaryConditions, const flag_t flag,
                                                                           const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                           const BoundaryConfiguration & parameter )
 {
@@ -2702,8 +2703,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( boost:
       setBoundary( boundaryConditions.get_tail(), flag, x, y, z, parameter );
 }
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const boost::tuples::null_type &, const flag_t flag,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setBoundary( const boost::tuples::null_type &, const flag_t flag,
                                                                           const cell_idx_t, const cell_idx_t, const cell_idx_t,
                                                                           const BoundaryConfiguration & ) const
 {
@@ -2719,9 +2720,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( boost::tuples::cons<Head, Tail> & boundaryConditions, const flag_t flag,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setBoundary( boost::tuples::cons<Head, Tail> & boundaryConditions, const flag_t flag,
                                                                           const CellInterval & cells, const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( outerBB_.contains( cells ) );
@@ -2754,8 +2755,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( boost:
       setBoundary( boundaryConditions.get_tail(), flag, cells, parameter );
 }
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const boost::tuples::null_type &, const flag_t flag,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setBoundary( const boost::tuples::null_type &, const flag_t flag,
                                                                           const CellInterval &, const BoundaryConfiguration & ) const
 {
    if( flagField_->isRegistered( flag ) )
@@ -2770,9 +2771,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( boost::tuples::cons<Head, Tail> & boundaryConditions, const flag_t flag,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setBoundary( boost::tuples::cons<Head, Tail> & boundaryConditions, const flag_t flag,
                                                                           const CellVector & cells, const BoundaryConfiguration & parameter )
 {
    if( cells.empty() )
@@ -2790,8 +2791,8 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( boost:
       setBoundary( boundaryConditions.get_tail(), flag, cells, parameter );
 }
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const boost::tuples::null_type &, const flag_t flag,
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::setBoundary( const boost::tuples::null_type &, const flag_t flag,
                                                                           const CellVector &, const BoundaryConfiguration & ) const
 {
    if( flagField_->isRegistered( flag ) )
@@ -2806,9 +2807,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::setBoundary( const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( boost::tuples::cons<Head, Tail> & boundaryConditions,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::removeBoundary( boost::tuples::cons<Head, Tail> & boundaryConditions,
                                                                              const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                              const bool checkNearBoundaryFlags )
 {
@@ -2864,9 +2865,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::removeBoundary( boo
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::treatDirection( boost::tuples::cons<Head, Tail> & boundaryConditions, const uint_t index,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::treatDirection( boost::tuples::cons<Head, Tail> & boundaryConditions, const uint_t index,
                                                                              const std::vector< std::vector< std::pair< Cell, stencil::Direction > > > & cellDirectionPairs )
 {
    Head & boundaryCondition = boundaryConditions.get_head();
@@ -2897,9 +2898,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::treatDirection( boo
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::treatDirection( boost::tuples::cons<Head, Tail> & boundaryConditions,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::treatDirection( boost::tuples::cons<Head, Tail> & 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 )
@@ -2918,9 +2919,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::treatDirection( boo
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::beforeBoundaryTreatment( boost::tuples::cons<Head, Tail> & boundaryConditions )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::beforeBoundaryTreatment( boost::tuples::cons<Head, Tail> & boundaryConditions )
 {
    boundaryConditions.get_head().beforeBoundaryTreatment();
 
@@ -2929,9 +2930,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::beforeBoundaryTreat
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::afterBoundaryTreatment( boost::tuples::cons<Head, Tail> & boundaryConditions )
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::afterBoundaryTreatment( boost::tuples::cons<Head, Tail> & boundaryConditions )
 {
    boundaryConditions.get_head().afterBoundaryTreatment();
 
@@ -2940,9 +2941,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::afterBoundaryTreatm
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline std::map< std::string, typename BoundaryHandling< FlagField_T, Stencil, Tuple >::flag_t >
-BoundaryHandling< FlagField_T, Stencil, Tuple >::getFlagMapping() const
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline std::map< std::string, typename BoundaryHandling< FlagField_T, Stencil, Boundaries... >::flag_t >
+BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getFlagMapping() const
 {
    std::map< std::string, flag_t > mapping;
    const auto uidMapping = flagField_->getMapping();
@@ -2953,10 +2954,10 @@ BoundaryHandling< FlagField_T, Stencil, Tuple >::getFlagMapping() const
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Buffer_T >
-std::vector< typename BoundaryHandling< FlagField_T, Stencil, Tuple >::flag_t >
-BoundaryHandling< FlagField_T, Stencil, Tuple >::getNeighborFlagMapping( Buffer_T & buffer, const bool assumeIdenticalFlagMapping,
+std::vector< typename BoundaryHandling< FlagField_T, Stencil, Boundaries... >::flag_t >
+BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getNeighborFlagMapping( Buffer_T & buffer, const bool assumeIdenticalFlagMapping,
                                                                          bool & identicalFlagMapping ) const
 {
    identicalFlagMapping = assumeIdenticalFlagMapping;
@@ -3000,8 +3001,8 @@ BoundaryHandling< FlagField_T, Stencil, Tuple >::getNeighborFlagMapping( Buffer_
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-void BoundaryHandling< FlagField_T, Stencil, Tuple >::translateMask( flag_t & mask, const std::vector< flag_t > & flagMapping ) const
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::translateMask( flag_t & mask, const std::vector< flag_t > & flagMapping ) const
 {
    flag_t neighbor( mask );
    mask = static_cast< flag_t >(0);
@@ -3012,8 +3013,8 @@ void BoundaryHandling< FlagField_T, Stencil, Tuple >::translateMask( flag_t & ma
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-CellInterval BoundaryHandling< FlagField_T, Stencil, Tuple >::getPackingInterval( stencil::Direction direction, const uint_t numberOfLayers ) const
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+CellInterval BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getPackingInterval( stencil::Direction direction, const uint_t numberOfLayers ) const
 {
    CellInterval interval = getUnpackingInterval( direction, numberOfLayers );
 
@@ -3029,8 +3030,8 @@ CellInterval BoundaryHandling< FlagField_T, Stencil, Tuple >::getPackingInterval
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-CellInterval BoundaryHandling< FlagField_T, Stencil, Tuple >::getUnpackingInterval( stencil::Direction direction, const uint_t numberOfLayers ) const
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+CellInterval BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getUnpackingInterval( stencil::Direction direction, const uint_t numberOfLayers ) const
 {
    WALBERLA_ASSERT_GREATER_EQUAL( numberOfLayers, uint_t(1) );
    WALBERLA_ASSERT( stencil::cx[direction] == 0 || outerBB_.xSize() >= ( uint_c(4) * numberOfLayers ) );
@@ -3058,9 +3059,9 @@ CellInterval BoundaryHandling< FlagField_T, Stencil, Tuple >::getUnpackingInterv
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Buffer_T >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::pack( Buffer_T & buffer, const flag_t mask,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::pack( Buffer_T & buffer, const flag_t mask,
                                                                    const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
 {
    if( field::isPartOfMaskSet( mask, boundary_ ) )
@@ -3069,9 +3070,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::pack( Buffer_T & bu
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Buffer_T, typename Head, typename Tail >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::pack( const boost::tuples::cons<Head, Tail> & boundaryConditions, Buffer_T & buffer,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::pack( const boost::tuples::cons<Head, Tail> & boundaryConditions, Buffer_T & buffer,
                                                                    const flag_t mask, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
 {
    const Head & boundaryCondition = boundaryConditions.get_head();
@@ -3085,9 +3086,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::pack( const boost::
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Buffer_T >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::unpack( Buffer_T & buffer, const flag_t flag,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::unpack( Buffer_T & buffer, const flag_t flag,
                                                                      const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT( field::isFlag(flag) );
@@ -3100,9 +3101,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::unpack( Buffer_T &
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Buffer_T >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::unpackBoundary( Buffer_T & buffer, const flag_t flag,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::unpackBoundary( Buffer_T & buffer, const flag_t flag,
                                                                              const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT_EQUAL( flag & boundary_, flag );
@@ -3115,9 +3116,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::unpackBoundary( Buf
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Buffer_T, typename Head, typename Tail >
-inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::unpackBoundary( boost::tuples::cons<Head, Tail> & boundaryConditions, Buffer_T & buffer,
+inline void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::unpackBoundary( boost::tuples::cons<Head, Tail> & boundaryConditions, Buffer_T & buffer,
                                                                              const flag_t flag,
                                                                              const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
@@ -3136,9 +3137,9 @@ inline void BoundaryHandling< FlagField_T, Stencil, Tuple >::unpackBoundary( boo
 
 
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
 template< typename Head, typename Tail >
-void BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryConditions( const boost::tuples::cons<Head, Tail> & boundaryConditions,
+void BoundaryHandling< FlagField_T, Stencil, Boundaries... >::getBoundaryConditions( const boost::tuples::cons<Head, Tail> & boundaryConditions,
                                                                              std::vector< std::string > & bcs ) const
 {
    const Head & boundaryCondition = boundaryConditions.get_head();
@@ -3172,8 +3173,8 @@ void BoundaryHandling< FlagField_T, Stencil, Tuple >::getBoundaryConditions( con
 
 using boundary::BoundaryHandling;
 
-template< typename FlagField_T, typename Stencil, typename Tuple >
-inline std::ostream & operator<<( std::ostream & os, const BoundaryHandling< FlagField_T, Stencil, Tuple > & boundaryHandling ) {
+template< typename FlagField_T, typename Stencil, typename... Boundaries >
+inline std::ostream & operator<<( std::ostream & os, const BoundaryHandling< FlagField_T, Stencil, Boundaries... > & boundaryHandling ) {
 
    boundaryHandling.toStream( os );
    return os;
diff --git a/src/boundary/BoundaryHandlingCollection.h b/src/boundary/BoundaryHandlingCollection.h
index b30ba815e269e241b76cb730f0c50c82e2503701..b9128661e40b32f559d2cfe7de8373b2f4ee22a6 100644
--- a/src/boundary/BoundaryHandlingCollection.h
+++ b/src/boundary/BoundaryHandlingCollection.h
@@ -53,7 +53,7 @@ typedef UID< BHCUIDGenerator > BoundaryHandlingCollectionUID;
 
 
 
-template< typename FlagField_T, typename Tuple > // Tuple: all the boundary handlers that are considered by this boundary handling collection
+template< typename FlagField_T, typename... Handlers > // Handlers: all the boundary handlers that are considered by this boundary handling collection
 class BoundaryHandlingCollection
 {
 public:
@@ -79,7 +79,7 @@ public:
 
 
 
-   BoundaryHandlingCollection( const std::string & identifier, FlagField_T * const flagField, const Tuple & boundaryHandlers );
+   BoundaryHandlingCollection( const std::string & identifier, FlagField_T * const flagField, const Handlers & ... boundaryHandlers );
 
    bool operator==( const BoundaryHandlingCollection &     ) const { WALBERLA_ASSERT( false ); return false; } // For testing purposes, block data items must be comparable with operator "==".
                                                                                                       // Since instances of type "BoundaryHandlingCollection" are registered as block data items,
@@ -524,22 +524,23 @@ private:
 
    const CellInterval outerBB_;
 
+   typedef boost::tuples::tuple<Handlers...> Tuple;
    Tuple boundaryHandlers_;
 
 }; // class BoundaryHandlingCollection
 
 
 
-template< typename FlagField_T, typename Tuple >
-BoundaryHandlingCollection< FlagField_T, Tuple >::BoundaryHandlingCollection( const std::string & identifier, FlagField_T * const flagField,
-                                                                              const Tuple & boundaryHandlers ) :
+template< typename FlagField_T, typename... Handlers >
+BoundaryHandlingCollection< FlagField_T, Handlers... >::BoundaryHandlingCollection( const std::string & identifier, FlagField_T * const flagField,
+                                                                              const Handlers & ... boundaryHandlers ) :
 
    uid_( identifier ),
    flagField_( flagField ),
    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 ),
-   boundaryHandlers_( boundaryHandlers )
+   boundaryHandlers_( boost::make_tuple(boundaryHandlers...) )
 {
    if( flagField_->nrOfGhostLayers() < 1 )
       WALBERLA_ABORT( "The flag field passed to the boundary handling collection\"" << identifier << "\" must contain at least one ghost layer!" );
@@ -560,8 +561,8 @@ BoundaryHandlingCollection< FlagField_T, Tuple >::BoundaryHandlingCollection( co
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::isEmpty( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
+template< typename FlagField_T, typename... Handlers >
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::isEmpty( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
 {
    WALBERLA_ASSERT( outerBB_.contains(x,y,z) );
 
@@ -570,8 +571,8 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::isEmpty( const cel
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::isEmpty( const ConstFlagFieldBaseIterator & it ) const
+template< typename FlagField_T, typename... Handlers >
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::isEmpty( const ConstFlagFieldBaseIterator & it ) const
 {
    WALBERLA_ASSERT_EQUAL( it.getField(), flagField_ );
    WALBERLA_ASSERT( outerBB_.contains(it.x(),it.y(),it.z()) );
@@ -581,8 +582,8 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::isEmpty( const Con
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHandlers( const uint_t numberOfGhostLayersToInclude ) const
+template< typename FlagField_T, typename... Handlers >
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::consideredByAllHandlers( const uint_t numberOfGhostLayersToInclude ) const
 {
    CellInterval cells = getGhostLayerCellInterval( numberOfGhostLayersToInclude );
    return consideredByAllHandlers( cells );
@@ -590,8 +591,8 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHan
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHandlers( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
+template< typename FlagField_T, typename... Handlers >
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::consideredByAllHandlers( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
 {
    WALBERLA_ASSERT( outerBB_.contains(x,y,z) );
 
@@ -600,8 +601,8 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHan
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHandlers( const ConstFlagFieldBaseIterator & it ) const
+template< typename FlagField_T, typename... Handlers >
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::consideredByAllHandlers( const ConstFlagFieldBaseIterator & it ) const
 {
    WALBERLA_ASSERT_EQUAL( it.getField(), flagField_ );
    WALBERLA_ASSERT( outerBB_.contains(it.x(),it.y(),it.z()) );
@@ -611,8 +612,8 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHan
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHandlers( const CellInterval & cells ) const
+template< typename FlagField_T, typename... Handlers >
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::consideredByAllHandlers( const CellInterval & cells ) const
 {
    WALBERLA_ASSERT( outerBB_.contains( cells ) );
 
@@ -625,9 +626,9 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHan
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename CellIterator >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHandlers( const CellIterator & begin, const CellIterator & end ) const
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::consideredByAllHandlers( const CellIterator & begin, const CellIterator & end ) const
 {
    for( auto cell = begin; cell != end; ++cell )
       if( !consideredByAllHandlers( cell->x(), cell->y(), cell->z() ) )
@@ -637,26 +638,26 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHan
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename BoundaryHandling_T >
-inline const BoundaryHandling_T & BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundaryHandling( const BoundaryHandlingUID & uid ) const
+inline const BoundaryHandling_T & BoundaryHandlingCollection< FlagField_T, Handlers... >::getBoundaryHandling( const BoundaryHandlingUID & uid ) const
 {
    return getBoundaryHandling< BoundaryHandling_T, typename Tuple::head_type, typename Tuple::tail_type >( uid, boundaryHandlers_ );
 }
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename BoundaryHandling_T >
-inline BoundaryHandling_T & BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundaryHandling( const BoundaryHandlingUID & uid )
+inline BoundaryHandling_T & BoundaryHandlingCollection< FlagField_T, Handlers... >::getBoundaryHandling( const BoundaryHandlingUID & uid )
 {
    return const_cast< BoundaryHandling_T & >( static_cast< const BoundaryHandlingCollection * >( this )->template getBoundaryHandling< BoundaryHandling_T >( uid ) );
 }
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatchingHandlers( const flag_t flag ) const
+template< typename FlagField_T, typename... Handlers >
+inline uint_t BoundaryHandlingCollection< FlagField_T, Handlers... >::numberOfMatchingHandlers( const flag_t flag ) const
 {
    WALBERLA_ASSERT( field::isFlag(flag) );
    return numberOfMatchingHandlers( boundaryHandlers_, flag );
@@ -664,8 +665,8 @@ inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatching
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatchingHandlersForDomain( const flag_t flag ) const
+template< typename FlagField_T, typename... Handlers >
+inline uint_t BoundaryHandlingCollection< FlagField_T, Handlers... >::numberOfMatchingHandlersForDomain( const flag_t flag ) const
 {
    WALBERLA_ASSERT( field::isFlag(flag) );
    return numberOfMatchingHandlersForDomain( boundaryHandlers_, flag );
@@ -673,8 +674,8 @@ inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatching
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatchingHandlersForBoundary( const flag_t flag ) const
+template< typename FlagField_T, typename... Handlers >
+inline uint_t BoundaryHandlingCollection< FlagField_T, Handlers... >::numberOfMatchingHandlersForBoundary( const flag_t flag ) const
 {
    WALBERLA_ASSERT( field::isFlag(flag) );
    return numberOfMatchingHandlersForBoundary( boundaryHandlers_, flag );
@@ -682,16 +683,16 @@ inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatching
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::containsBoundaryCondition( const BoundaryUID & uid ) const
+template< typename FlagField_T, typename... Handlers >
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::containsBoundaryCondition( const BoundaryUID & uid ) const
 {
    return containsBoundaryCondition( boundaryHandlers_, uid );
 }
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::containsBoundaryCondition( const FlagUID & flag ) const
+template< typename FlagField_T, typename... Handlers >
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::containsBoundaryCondition( const FlagUID & flag ) const
 {
    if( flagField_->flagExists( flag ) )
       return containsBoundaryCondition( flagField_->getFlag( flag ) );
@@ -700,16 +701,16 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::containsBoundaryCo
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::containsBoundaryCondition( const flag_t flag ) const
+template< typename FlagField_T, typename... Handlers >
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::containsBoundaryCondition( const flag_t flag ) const
 {
    return containsBoundaryCondition( boundaryHandlers_, flag );
 }
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline BoundaryUID BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundaryUID( const FlagUID & flag ) const
+template< typename FlagField_T, typename... Handlers >
+inline BoundaryUID BoundaryHandlingCollection< FlagField_T, Handlers... >::getBoundaryUID( const FlagUID & flag ) const
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
 
@@ -718,8 +719,8 @@ inline BoundaryUID BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundary
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline BoundaryUID BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundaryUID( const flag_t flag ) const
+template< typename FlagField_T, typename... Handlers >
+inline BoundaryUID BoundaryHandlingCollection< FlagField_T, Handlers... >::getBoundaryUID( const flag_t flag ) const
 {
    WALBERLA_ASSERT( field::isFlag( flag ) );
    WALBERLA_ASSERT( flagField_->isRegistered( flag ) );
@@ -729,8 +730,8 @@ inline BoundaryUID BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundary
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::checkConsistency( const uint_t numberOfGhostLayersToInclude ) const
+template< typename FlagField_T, typename... Handlers >
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::checkConsistency( const uint_t numberOfGhostLayersToInclude ) const
 {
    CellInterval cells = getGhostLayerCellInterval( numberOfGhostLayersToInclude );
    return checkConsistency( cells );
@@ -738,16 +739,16 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::checkConsistency(
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::checkConsistency( const CellInterval & cells ) const
+template< typename FlagField_T, typename... Handlers >
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::checkConsistency( const CellInterval & cells ) const
 {
    return checkConsistency( boundaryHandlers_, cells );
 }
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::refresh( const uint_t numberOfGhostLayersToInclude )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::refresh( const uint_t numberOfGhostLayersToInclude )
 {
    CellInterval cells = getGhostLayerCellInterval( numberOfGhostLayersToInclude );
    refresh( cells );
@@ -755,24 +756,24 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::refresh( const uin
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::refresh( const CellInterval & cells )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::refresh( const CellInterval & cells )
 {
    refresh( boundaryHandlers_, cells );
 }
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::refreshOutermostLayer( cell_idx_t thickness  )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::refreshOutermostLayer( cell_idx_t thickness  )
 {
    refreshOutermostLayer( boundaryHandlers_, thickness );
 }
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::setFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                        const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -782,8 +783,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const Fla
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::setFlag( const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                        const BoundaryConfiguration & parameter )
 {
    if( !outerBB_.contains(x,y,z) )
@@ -796,8 +797,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const fla
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const FlagUID & flag, const CellInterval & cells,
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::setFlag( const FlagUID & flag, const CellInterval & cells,
                                                                        const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -807,8 +808,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const Fla
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const flag_t flag, const CellInterval & cells,
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::setFlag( const flag_t flag, const CellInterval & cells,
                                                                        const BoundaryConfiguration & parameter )
 {
    CellInterval localCells( outerBB_ );
@@ -827,9 +828,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const fla
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename CellIterator >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::setFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
                                                                        const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -839,9 +840,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const Fla
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename CellIterator >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const flag_t flag, const CellIterator & begin, const CellIterator & end,
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::setFlag( const flag_t flag, const CellIterator & begin, const CellIterator & end,
                                                                        const BoundaryConfiguration & parameter )
 {
    for( auto cell = begin; cell != end; ++cell )
@@ -850,8 +851,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const fla
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::forceFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                          const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -861,8 +862,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlag( const F
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlag( const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::forceFlag( const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                          const BoundaryConfiguration & parameter )
 {
    if( !outerBB_.contains(x,y,z) )
@@ -873,8 +874,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlag( const f
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlag( const FlagUID & flag, const CellInterval & cells,
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::forceFlag( const FlagUID & flag, const CellInterval & cells,
                                                                          const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -884,8 +885,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlag( const F
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlag( const flag_t flag, const CellInterval & cells,
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::forceFlag( const flag_t flag, const CellInterval & cells,
                                                                          const BoundaryConfiguration & parameter )
 {
    CellInterval localCells( outerBB_ );
@@ -902,9 +903,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlag( const f
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename CellIterator >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::forceFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end,
                                                                          const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
@@ -914,9 +915,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlag( const F
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename CellIterator >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlag( const flag_t flag, const CellIterator & begin, const CellIterator & end,
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::forceFlag( const flag_t flag, const CellIterator & begin, const CellIterator & end,
                                                                          const BoundaryConfiguration & parameter )
 {
    for( auto cell = begin; cell != end; ++cell )
@@ -925,8 +926,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlag( const f
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const FlagUID & flag, const uint_t numberOfGhostLayersToInclude )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::removeFlag( const FlagUID & flag, const uint_t numberOfGhostLayersToInclude )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
 
@@ -935,8 +936,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const flag_t flag, const uint_t numberOfGhostLayersToInclude )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::removeFlag( const flag_t flag, const uint_t numberOfGhostLayersToInclude )
 {
    CellInterval cells = getGhostLayerCellInterval( numberOfGhostLayersToInclude );
    removeFlag( flag, cells );
@@ -944,8 +945,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::removeFlag( const FlagUID & flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
 
@@ -954,8 +955,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::removeFlag( const flag_t flag, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    if( !outerBB_.contains(x,y,z) || !flagField_->isFlagSet( x, y, z, flag ) )
       return;
@@ -965,8 +966,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const FlagUID & flag, const CellInterval & cells )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::removeFlag( const FlagUID & flag, const CellInterval & cells )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
 
@@ -975,8 +976,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const flag_t flag, const CellInterval & cells )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::removeFlag( const flag_t flag, const CellInterval & cells )
 {
    CellInterval localCells( outerBB_ );
    localCells.intersect( cells );
@@ -993,9 +994,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename CellIterator >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::removeFlag( const FlagUID & flag, const CellIterator & begin, const CellIterator & end )
 {
    WALBERLA_ASSERT( flagField_->flagExists( flag ) );
 
@@ -1004,9 +1005,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename CellIterator >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const flag_t flag, const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::removeFlag( const flag_t flag, const CellIterator & begin, const CellIterator & end )
 {
    for( auto cell = begin; cell != end; ++cell )
       removeFlag( flag, cell->x(), cell->y(), cell->z() );
@@ -1014,8 +1015,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::clear( const uint_t numberOfGhostLayersToInclude )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::clear( const uint_t numberOfGhostLayersToInclude )
 {
    CellInterval cells = getGhostLayerCellInterval( numberOfGhostLayersToInclude );
    clear( cells );
@@ -1023,8 +1024,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::clear( const uint_
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::clear( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::clear( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    if( !outerBB_.contains(x,y,z) )
       return;
@@ -1034,8 +1035,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::clear( const cell_
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::clear( const CellInterval & cells )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::clear( const CellInterval & cells )
 {
    CellInterval localCells( outerBB_ );
    localCells.intersect( cells );
@@ -1051,9 +1052,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::clear( const CellI
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename CellIterator >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::clear( const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::clear( const CellIterator & begin, const CellIterator & end )
 {
    for( auto cell = begin; cell != end; ++cell )
       clear( cell->x(), cell->y(), cell->z() );
@@ -1061,42 +1062,42 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::clear( const CellI
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::operator()( const uint_t numberOfGhostLayersToInclude )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::operator()( const uint_t numberOfGhostLayersToInclude )
 {
    execute( boundaryHandlers_, numberOfGhostLayersToInclude );
 }
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::operator()( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::operator()( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    execute( boundaryHandlers_, x, y, z );
 }
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::operator()( const CellInterval & cells )
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::operator()( const CellInterval & cells )
 {
    execute( boundaryHandlers_, cells );
 }
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename CellIterator >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::operator()( const CellIterator & begin, const CellIterator & end )
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::operator()( const CellIterator & begin, const CellIterator & end )
 {
    execute( boundaryHandlers_, begin, end );
 }
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Buffer_T >
-void BoundaryHandlingCollection< FlagField_T, Tuple >::pack( Buffer_T & buffer, stencil::Direction direction, const uint_t numberOfLayers,
+void BoundaryHandlingCollection< FlagField_T, Handlers... >::pack( Buffer_T & buffer, stencil::Direction direction, const uint_t numberOfLayers,
                                                              const bool assumeIdenticalFlagMapping ) const
 {
 #ifdef NDEBUG
@@ -1120,9 +1121,9 @@ void BoundaryHandlingCollection< FlagField_T, Tuple >::pack( Buffer_T & buffer,
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Buffer_T >
-void BoundaryHandlingCollection< FlagField_T, Tuple >::unpack( Buffer_T & buffer, stencil::Direction direction, const uint_t numberOfLayers,
+void BoundaryHandlingCollection< FlagField_T, Handlers... >::unpack( Buffer_T & buffer, stencil::Direction direction, const uint_t numberOfLayers,
                                                                const bool assumeIdenticalFlagMapping )
 {
    bool identicalFlagMapping = false;
@@ -1150,8 +1151,8 @@ void BoundaryHandlingCollection< FlagField_T, Tuple >::unpack( Buffer_T & buffer
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::toStream( std::ostream & os ) const {
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::toStream( std::ostream & os ) const {
 
    os << "========================= BoundaryHandlingCollection =========================\n\n"
       << "Identifier: " << uid_.getIdentifier() << "\n\n"
@@ -1164,8 +1165,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::toStream( std::ost
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline std::string BoundaryHandlingCollection< FlagField_T, Tuple >::toString() const {
+template< typename FlagField_T, typename... Handlers >
+inline std::string BoundaryHandlingCollection< FlagField_T, Handlers... >::toString() const {
 
    std::ostringstream oss;
    toStream( oss );
@@ -1174,8 +1175,8 @@ inline std::string BoundaryHandlingCollection< FlagField_T, Tuple >::toString()
 
 
 
-template< typename FlagField_T, typename Tuple >
-CellInterval BoundaryHandlingCollection< FlagField_T, Tuple >::getGhostLayerCellInterval( const uint_t numberOfGhostLayersToInclude ) const
+template< typename FlagField_T, typename... Handlers >
+CellInterval BoundaryHandlingCollection< FlagField_T, Handlers... >::getGhostLayerCellInterval( const uint_t numberOfGhostLayersToInclude ) const
 {
    CellInterval cells( -cell_idx_c( numberOfGhostLayersToInclude ),
                        -cell_idx_c( numberOfGhostLayersToInclude ),
@@ -1188,9 +1189,9 @@ CellInterval BoundaryHandlingCollection< FlagField_T, Tuple >::getGhostLayerCell
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::isEmpty( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::isEmpty( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                        const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
 {
    if( !(boundaryHandlers.get_head().isEmpty(x,y,z)) )
@@ -1200,9 +1201,9 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::isEmpty( const boo
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::isEmpty( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::isEmpty( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                        const ConstFlagFieldBaseIterator & it ) const
 {
    if( !(boundaryHandlers.get_head().isEmpty(it)) )
@@ -1212,9 +1213,9 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::isEmpty( const boo
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHandlers( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::consideredByAllHandlers( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                                          const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
 {
    if( boundaryHandlers.get_head().isEmpty(x,y,z) )
@@ -1224,9 +1225,9 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHan
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHandlers( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::consideredByAllHandlers( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                                          const ConstFlagFieldBaseIterator & it ) const
 {
    if( boundaryHandlers.get_head().isEmpty(it) )
@@ -1236,9 +1237,9 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::consideredByAllHan
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::checkForUniqueBoundaryHandlingUIDs( const boost::tuples::cons<Head, Tail> & boundaryHandlers ) const
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::checkForUniqueBoundaryHandlingUIDs( const boost::tuples::cons<Head, Tail> & boundaryHandlers ) const
 {
    if( numberOfMatchingBoundaryHandlers( boundaryHandlers.get_head().getUID() ) != uint_c(1) )
       WALBERLA_ABORT( "Every boundary handler registered at the same boundary handling collection must have a unique boundary handling UID!\n"
@@ -1249,8 +1250,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::checkForUniqueBoun
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline std::vector< BoundaryUID > BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundaryUIDs() const
+template< typename FlagField_T, typename... Handlers >
+inline std::vector< BoundaryUID > BoundaryHandlingCollection< FlagField_T, Handlers... >::getBoundaryUIDs() const
 {
    std::vector< BoundaryUID > uids;
    getBoundaryUIDs( boundaryHandlers_, uids );
@@ -1259,9 +1260,9 @@ inline std::vector< BoundaryUID > BoundaryHandlingCollection< FlagField_T, Tuple
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundaryUIDs( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::getBoundaryUIDs( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                                std::vector< BoundaryUID > & uids ) const
 {
    std::vector< BoundaryUID > handlerUIDs = boundaryHandlers.get_head().getBoundaryUIDs();
@@ -1271,9 +1272,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundaryUIDs( c
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::checkForIdenticalFlagFields( const boost::tuples::cons<Head, Tail> & boundaryHandlers ) const
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::checkForIdenticalFlagFields( const boost::tuples::cons<Head, Tail> & boundaryHandlers ) const
 {
    return checkForIdenticalFlagFields( boundaryHandlers.get_tail() ) &&
           boundaryHandlers.get_head().getFlagField() == flagField_ && boundaryHandlers.get_head().outerBB_ == outerBB_;
@@ -1281,17 +1282,17 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::checkForIdenticalF
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatchingBoundaryHandlers( const BoundaryHandlingUID & uid ) const
+template< typename FlagField_T, typename... Handlers >
+inline uint_t BoundaryHandlingCollection< FlagField_T, Handlers... >::numberOfMatchingBoundaryHandlers( const BoundaryHandlingUID & uid ) const
 {
    return numberOfMatchingBoundaryHandlers( boundaryHandlers_, uid );
 }
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatchingBoundaryHandlers( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline uint_t BoundaryHandlingCollection< FlagField_T, Handlers... >::numberOfMatchingBoundaryHandlers( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                                                   const BoundaryHandlingUID & uid ) const
 {
    return ( ( boundaryHandlers.get_head().getUID() == uid ) ? uint_c(1) : uint_c(0) ) +
@@ -1300,9 +1301,9 @@ inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatching
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatchingHandlers( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline uint_t BoundaryHandlingCollection< FlagField_T, Handlers... >::numberOfMatchingHandlers( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                                           const flag_t flag ) const
 {
    return ( ( (boundaryHandlers.get_head().getBoundaryMask() | boundaryHandlers.get_head().getDomainMask()) & flag ) == flag ? 1 : 0 ) +
@@ -1311,9 +1312,9 @@ inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatching
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatchingHandlersForDomain( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline uint_t BoundaryHandlingCollection< FlagField_T, Handlers... >::numberOfMatchingHandlersForDomain( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                                                    const flag_t flag ) const
 {
    return ( ( boundaryHandlers.get_head().getDomainMask() & flag ) == flag ? 1 : 0 ) +
@@ -1322,9 +1323,9 @@ inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatching
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatchingHandlersForBoundary( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline uint_t BoundaryHandlingCollection< FlagField_T, Handlers... >::numberOfMatchingHandlersForBoundary( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                                                      const flag_t flag ) const
 {
    return ( ( boundaryHandlers.get_head().getBoundaryMask() & flag ) == flag ? 1 : 0 ) +
@@ -1333,9 +1334,9 @@ inline uint_t BoundaryHandlingCollection< FlagField_T, Tuple >::numberOfMatching
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::containsBoundaryCondition( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::containsBoundaryCondition( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                                          const BoundaryUID & uid ) const
 {
    if( boundaryHandlers.get_head().containsBoundaryCondition( uid ) )
@@ -1345,9 +1346,9 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::containsBoundaryCo
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::containsBoundaryCondition( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::containsBoundaryCondition( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                                          const flag_t flag ) const
 {
    if( boundaryHandlers.get_head().containsBoundaryCondition( flag ) )
@@ -1357,10 +1358,10 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::containsBoundaryCo
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline typename BoundaryHandlingCollection< FlagField_T, Tuple >::flag_t
-   BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundaryMask( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline typename BoundaryHandlingCollection< FlagField_T, Handlers... >::flag_t
+   BoundaryHandlingCollection< FlagField_T, Handlers... >::getBoundaryMask( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                       const BoundaryUID & uid ) const
 {
    if( boundaryHandlers.get_head().containsBoundaryCondition(uid) )
@@ -1370,9 +1371,9 @@ inline typename BoundaryHandlingCollection< FlagField_T, Tuple >::flag_t
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline BoundaryUID BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundaryUID( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline BoundaryUID BoundaryHandlingCollection< FlagField_T, Handlers... >::getBoundaryUID( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                                      const flag_t flag ) const
 {
    const Head & boundaryHandler = boundaryHandlers.get_head();
@@ -1389,8 +1390,8 @@ inline BoundaryUID BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundary
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline BoundaryUID BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundaryUID( const boost::tuples::null_type &, const flag_t flag ) const
+template< typename FlagField_T, typename... Handlers >
+inline BoundaryUID BoundaryHandlingCollection< FlagField_T, Handlers... >::getBoundaryUID( const boost::tuples::null_type &, const flag_t flag ) const
 {
    if( !flagField_->isRegistered( flag ) )
       WALBERLA_ABORT( "The requested flag with value " << flag << " is not registered at the flag field and is not handled "\
@@ -1403,10 +1404,10 @@ inline BoundaryUID BoundaryHandlingCollection< FlagField_T, Tuple >::getBoundary
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
 inline shared_ptr<BoundaryConfiguration>
-   BoundaryHandlingCollection< FlagField_T, Tuple >::createBoundaryConfiguration( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+   BoundaryHandlingCollection< FlagField_T, Handlers... >::createBoundaryConfiguration( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                                   const BoundaryUID & uid, const Config::BlockHandle & config ) const
 {
    if( boundaryHandlers.get_head().containsBoundaryCondition(uid) )
@@ -1416,9 +1417,9 @@ inline shared_ptr<BoundaryConfiguration>
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::checkConsistency( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline bool BoundaryHandlingCollection< FlagField_T, Handlers... >::checkConsistency( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                                 const CellInterval & cells ) const
 {
    return checkConsistency( boundaryHandlers.get_tail(), cells ) && boundaryHandlers.get_head().checkConsistency(cells);
@@ -1426,9 +1427,9 @@ inline bool BoundaryHandlingCollection< FlagField_T, Tuple >::checkConsistency(
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::refresh( boost::tuples::cons<Head, Tail> & boundaryHandlers, const CellInterval & cells )
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::refresh( boost::tuples::cons<Head, Tail> & boundaryHandlers, const CellInterval & cells )
 {
    boundaryHandlers.get_head().refresh( cells );
    refresh( boundaryHandlers.get_tail(), cells );
@@ -1436,9 +1437,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::refresh( boost::tu
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::refreshOutermostLayer( boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::refreshOutermostLayer( boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                                      cell_idx_t thickness )
 {
    boundaryHandlers.get_head().refreshOutermostLayer( thickness );
@@ -1447,9 +1448,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::refreshOutermostLa
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( boost::tuples::cons<Head, Tail> & boundaryHandlers, const flag_t flag,
+void BoundaryHandlingCollection< FlagField_T, Handlers... >::setFlag( boost::tuples::cons<Head, Tail> & boundaryHandlers, const flag_t flag,
                                                                 const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                 const BoundaryConfiguration & parameter )
 {
@@ -1468,8 +1469,8 @@ void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( boost::tuples::c
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const boost::tuples::null_type & /*boundaryHandlers*/, const flag_t flag,
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::setFlag( const boost::tuples::null_type & /*boundaryHandlers*/, const flag_t flag,
                                                                        const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
                                                                        const BoundaryConfiguration & /*parameter*/ )
 {
@@ -1480,9 +1481,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const boo
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( boost::tuples::cons<Head, Tail> & boundaryHandlers, const flag_t flag,
+void BoundaryHandlingCollection< FlagField_T, Handlers... >::setFlag( boost::tuples::cons<Head, Tail> & boundaryHandlers, const flag_t flag,
                                                                 const CellInterval & cells, const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( outerBB_.contains(cells) );
@@ -1505,8 +1506,8 @@ void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( boost::tuples::c
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const boost::tuples::null_type & /*boundaryHandlers*/, const flag_t flag,
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::setFlag( const boost::tuples::null_type & /*boundaryHandlers*/, const flag_t flag,
                                                                        const CellInterval & cells, const BoundaryConfiguration & /*parameter*/ )
 {
    WALBERLA_ASSERT( outerBB_.contains(cells) );
@@ -1517,8 +1518,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::setFlag( const boo
 
 
 
-template< typename FlagField_T, typename Tuple >
-void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlagHelper( const flag_t flag, const cell_idx_t x, const cell_idx_t y,
+template< typename FlagField_T, typename... Handlers >
+void BoundaryHandlingCollection< FlagField_T, Handlers... >::forceFlagHelper( const flag_t flag, const cell_idx_t x, const cell_idx_t y,
                                                                         const cell_idx_t z, const BoundaryConfiguration & parameter )
 {
    WALBERLA_ASSERT( outerBB_.contains(x,y,z) );
@@ -1540,10 +1541,10 @@ void BoundaryHandlingCollection< FlagField_T, Tuple >::forceFlagHelper( const fl
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-typename BoundaryHandlingCollection< FlagField_T, Tuple >::flag_t
-BoundaryHandlingCollection< FlagField_T, Tuple >::flagsToRemove( boost::tuples::cons<Head, Tail> & boundaryHandlers, const flag_t flag,
+typename BoundaryHandlingCollection< FlagField_T, Handlers... >::flag_t
+BoundaryHandlingCollection< FlagField_T, Handlers... >::flagsToRemove( boost::tuples::cons<Head, Tail> & boundaryHandlers, const flag_t flag,
                                                                  const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT( outerBB_.contains(x,y,z) );
@@ -1561,9 +1562,9 @@ BoundaryHandlingCollection< FlagField_T, Tuple >::flagsToRemove( boost::tuples::
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( boost::tuples::cons<Head, Tail> & boundaryHandlers, const flag_t flag,
+void BoundaryHandlingCollection< FlagField_T, Handlers... >::removeFlag( boost::tuples::cons<Head, Tail> & boundaryHandlers, const flag_t flag,
                                                                    const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT( outerBB_.contains(x,y,z) );
@@ -1581,8 +1582,8 @@ void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( boost::tuples
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const boost::tuples::null_type & /*boundaryHandlers*/, const flag_t flag,
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::removeFlag( const boost::tuples::null_type & /*boundaryHandlers*/, const flag_t flag,
                                                                           const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT( outerBB_.contains(x,y,z) );
@@ -1592,10 +1593,10 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::removeFlag( const
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-typename BoundaryHandlingCollection< FlagField_T, Tuple >::flag_t
-BoundaryHandlingCollection< FlagField_T, Tuple >::clear( boost::tuples::cons<Head, Tail> & boundaryHandlers,
+typename BoundaryHandlingCollection< FlagField_T, Handlers... >::flag_t
+BoundaryHandlingCollection< FlagField_T, Handlers... >::clear( boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                          const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT( outerBB_.contains(x,y,z) );
@@ -1610,9 +1611,9 @@ BoundaryHandlingCollection< FlagField_T, Tuple >::clear( boost::tuples::cons<Hea
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::execute( boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::execute( boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                        const uint_t numberOfGhostLayersToInclude )
 {
    boundaryHandlers.get_head()( numberOfGhostLayersToInclude );
@@ -1621,9 +1622,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::execute( boost::tu
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::execute( boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::execute( boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                        const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    boundaryHandlers.get_head()(x,y,z);
@@ -1632,9 +1633,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::execute( boost::tu
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::execute( boost::tuples::cons<Head, Tail> & boundaryHandlers, const CellInterval & cells )
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::execute( boost::tuples::cons<Head, Tail> & boundaryHandlers, const CellInterval & cells )
 {
    boundaryHandlers.get_head()( cells );
    execute( boundaryHandlers.get_tail(), cells );
@@ -1642,9 +1643,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::execute( boost::tu
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename CellIterator, typename Head, typename Tail >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::execute( boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::execute( boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                        const CellIterator & begin, const CellIterator & end )
 {
    boundaryHandlers.get_head()( begin, end );
@@ -1653,9 +1654,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::execute( boost::tu
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::beforeBoundaryTreatment( boost::tuples::cons<Head, Tail> & boundaryHandlers )
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::beforeBoundaryTreatment( boost::tuples::cons<Head, Tail> & boundaryHandlers )
 {
    boundaryHandlers.get_head().beforeBoundaryTreatment();
    beforeBoundaryTreatment( boundaryHandlers.get_tail() );
@@ -1663,9 +1664,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::beforeBoundaryTrea
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::afterBoundaryTreatment( boost::tuples::cons<Head, Tail> & boundaryHandlers )
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::afterBoundaryTreatment( boost::tuples::cons<Head, Tail> & boundaryHandlers )
 {
    boundaryHandlers.get_head().afterBoundaryTreatment();
    afterBoundaryTreatment( boundaryHandlers.get_tail() );
@@ -1673,16 +1674,16 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::afterBoundaryTreat
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::translateMask( flag_t & mask, const std::vector< flag_t > & flagMapping ) const
+template< typename FlagField_T, typename... Handlers >
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::translateMask( flag_t & mask, const std::vector< flag_t > & flagMapping ) const
 {
    boundaryHandlers_.get_head().translateMask( mask, flagMapping );
 }
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline CellInterval BoundaryHandlingCollection< FlagField_T, Tuple >::getPackingInterval( stencil::Direction direction,
+template< typename FlagField_T, typename... Handlers >
+inline CellInterval BoundaryHandlingCollection< FlagField_T, Handlers... >::getPackingInterval( stencil::Direction direction,
                                                                                           const uint_t numberOfLayers ) const
 {
    return boundaryHandlers_.get_head().getPackingInterval( direction, numberOfLayers );
@@ -1690,8 +1691,8 @@ inline CellInterval BoundaryHandlingCollection< FlagField_T, Tuple >::getPacking
 
 
 
-template< typename FlagField_T, typename Tuple >
-inline CellInterval BoundaryHandlingCollection< FlagField_T, Tuple >::getUnpackingInterval( stencil::Direction direction,
+template< typename FlagField_T, typename... Handlers >
+inline CellInterval BoundaryHandlingCollection< FlagField_T, Handlers... >::getUnpackingInterval( stencil::Direction direction,
                                                                                             const uint_t numberOfLayers ) const
 {
    return boundaryHandlers_.get_head().getUnpackingInterval( direction, numberOfLayers );
@@ -1699,9 +1700,9 @@ inline CellInterval BoundaryHandlingCollection< FlagField_T, Tuple >::getUnpacki
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail, typename Buffer_T >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::pack( const boost::tuples::cons<Head, Tail> & boundaryHandlers, Buffer_T & buffer,
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::pack( const boost::tuples::cons<Head, Tail> & boundaryHandlers, Buffer_T & buffer,
                                                                     const flag_t mask, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
 {
    boundaryHandlers_.get_head().pack( buffer, mask, x, y, z );
@@ -1710,9 +1711,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::pack( const boost:
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail, typename Buffer_T >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::unpack( boost::tuples::cons<Head, Tail> & boundaryHandlers, Buffer_T & buffer,
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::unpack( boost::tuples::cons<Head, Tail> & boundaryHandlers, Buffer_T & buffer,
                                                                       const flag_t mask, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
 {
    WALBERLA_ASSERT( outerBB_.contains(x,y,z) );
@@ -1731,9 +1732,9 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::unpack( boost::tup
 
 
 
-template< typename FlagField_T, typename Tuple >
+template< typename FlagField_T, typename... Handlers >
 template< typename Head, typename Tail >
-inline void BoundaryHandlingCollection< FlagField_T, Tuple >::toStream( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
+inline void BoundaryHandlingCollection< FlagField_T, Handlers... >::toStream( const boost::tuples::cons<Head, Tail> & boundaryHandlers,
                                                                         std::ostream & os ) const
 {
    os << boundaryHandlers.get_head();
@@ -1748,8 +1749,8 @@ inline void BoundaryHandlingCollection< FlagField_T, Tuple >::toStream( const bo
 
 using boundary::BoundaryHandlingCollection;
 
-template< typename FlagField_T, typename Tuple >
-inline std::ostream & operator<<( std::ostream & os, const BoundaryHandlingCollection< FlagField_T, Tuple > & bhc ) {
+template< typename FlagField_T, typename... Handlers >
+inline std::ostream & operator<<( std::ostream & os, const BoundaryHandlingCollection< FlagField_T, Handlers... > & bhc ) {
 
    bhc.toStream( os );
    return os;
diff --git a/src/lbm/boundary/factories/DefaultBoundaryHandling.h b/src/lbm/boundary/factories/DefaultBoundaryHandling.h
index 3036fa2acb40ca5ba0f4e188c8141fccb728eefa..1e196a12c07d697ec4fb68364011aa66057c01a7 100644
--- a/src/lbm/boundary/factories/DefaultBoundaryHandling.h
+++ b/src/lbm/boundary/factories/DefaultBoundaryHandling.h
@@ -32,8 +32,6 @@
 #include "boundary/BoundaryHandling.h"
 #include "lbm/field/PdfField.h"
 
-#include <boost/tuple/tuple.hpp>
-
 namespace walberla {
 namespace lbm {
 
@@ -82,8 +80,7 @@ public:
    typedef Vector3<real_t>                        Velocity;
    typedef PdfField< LatticeModel >               PdfFieldLM;
 
-   typedef boost::tuple< BcNoSlip, BcFreeSlip, BcSimpleUBB, BcSimpleUBB, BcSimplePressure, BcSimplePressure> BoundaryConditions;
-   typedef walberla::boundary::BoundaryHandling< FlagFieldT, Stencil, BoundaryConditions >                   BoundaryHandling;
+   typedef walberla::boundary::BoundaryHandling< FlagFieldT, Stencil, BcNoSlip, BcFreeSlip, BcSimpleUBB, BcSimpleUBB, BcSimplePressure, BcSimplePressure > BoundaryHandling;
 
    static BlockDataID addBoundaryHandlingToStorage( const shared_ptr< StructuredBlockStorage > & bs, const std::string & identifier,
                                                     BlockDataID flagFieldID, BlockDataID pdfFieldID, const Set< FlagUID > & flagUIDSet,
@@ -165,14 +162,12 @@ DefaultBoundaryHandlingFactory<LatticeModel, FlagFieldT>::operator()( walberla::
       mask = static_cast< flag_t >( mask | flagField->getOrRegisterFlag( *flag ) );
 
    BoundaryHandling * handling = new BoundaryHandling( "default lbm boundary handling", flagField, mask,
-      boost::tuples::make_tuple(
         BcNoSlip        ( getNoSlipBoundaryUID(),    getNoSlip(),    pdfField ),
         BcFreeSlip      ( getFreeSlipBoundaryUID(),  getFreeSlip(),  pdfField, flagField, mask ),
         BcSimpleUBB     ( getVelocity0BoundaryUID(), getVelocity0(), pdfField, velocity0_ ),
         BcSimpleUBB     ( getVelocity1BoundaryUID(), getVelocity1(), pdfField, velocity1_ ),
         BcSimplePressure( getPressure0BoundaryUID(), getPressure0(), pdfField, pressure0_ ),
         BcSimplePressure( getPressure1BoundaryUID(), getPressure1(), pdfField, pressure1_ )
-      )
     );
 
    return handling;
diff --git a/src/lbm/boundary/factories/DefaultBoundaryHandlingCollection.h b/src/lbm/boundary/factories/DefaultBoundaryHandlingCollection.h
index 4bee0c6c165d99ce6c57d746ddef667106e8056f..54aafd3345681bea2ea7776ae6decfd98a0361f3 100644
--- a/src/lbm/boundary/factories/DefaultBoundaryHandlingCollection.h
+++ b/src/lbm/boundary/factories/DefaultBoundaryHandlingCollection.h
@@ -40,11 +40,8 @@ private:
    typedef typename DefaultBoundaryHandlingFactory         < LatticeModel_T,          FlagField_T >::BoundaryHandling    DefaultBoundaryHandling_T;
    typedef typename DefaultDiffusionBoundaryHandlingFactory< DiffusionLatticeModel_T, FlagField_T >::BoundaryHandling_T  DefaultDiffusionBoundaryHandlingFactory_T;
 
-private:
-   typedef boost::tuples::tuple< DefaultBoundaryHandling_T &, DefaultDiffusionBoundaryHandlingFactory_T & >  BoundaryHandlings;
-
 public:
-   typedef BoundaryHandlingCollection< FlagField_T, BoundaryHandlings > BoundaryHandlingCollection_T;
+   typedef BoundaryHandlingCollection< FlagField_T, DefaultBoundaryHandling_T &, DefaultDiffusionBoundaryHandlingFactory_T & > BoundaryHandlingCollection_T;
 
 
 private:
@@ -62,7 +59,7 @@ private:
       DefaultBoundaryHandling_T                 & handling          = * block->getData< DefaultBoundaryHandling_T                 >( handlingID          );
       DefaultDiffusionBoundaryHandlingFactory_T & diffusionHandling = * block->getData< DefaultDiffusionBoundaryHandlingFactory_T >( diffusionHandlingID );
             
-      return new BoundaryHandlingCollection_T( " Diffusion Boundary Handling Collection", flagField, BoundaryHandlings( handling, diffusionHandling ) );
+      return new BoundaryHandlingCollection_T( " Diffusion Boundary Handling Collection", flagField, handling, diffusionHandling );
    }
 
 public:
diff --git a/src/lbm/boundary/factories/DefaultDiffusionBoundaryHandling.h b/src/lbm/boundary/factories/DefaultDiffusionBoundaryHandling.h
index 42e3c7536baaf15361049affebc016eb46f7664d..e32559776caa7bf309cb6d723ecf0d7be605d086 100644
--- a/src/lbm/boundary/factories/DefaultDiffusionBoundaryHandling.h
+++ b/src/lbm/boundary/factories/DefaultDiffusionBoundaryHandling.h
@@ -34,8 +34,6 @@
 
 #include "field/FlagFunctions.h"
 
-#include <boost/tuple/tuple.hpp>
-
 #include <functional>
 
 
@@ -56,12 +54,9 @@ public:
    typedef lbm::NoDiffusion             < LatticeModel_T, flag_t      >  NoDiffusion_T;
    typedef lbm::FreeDiffusion           < LatticeModel_T, FlagField_T >  FreeDiffusion_T;
 
-private:
-   typedef boost::tuples::tuple< DiffusionDirichlet_T, SimpleDiffusionDirichlet_T, NoDiffusion_T, FreeDiffusion_T, SimpleDiffusionDirichlet_T, SimpleDiffusionDirichlet_T >  BoundaryConditions;
-
 public:
-   typedef walberla::boundary::BoundaryHandling< FlagField_T, Stencil, BoundaryConditions > BoundaryHandling_T;
-   typedef walberla::boundary::BoundaryHandling< FlagField_T, Stencil, BoundaryConditions > BoundaryHandling;
+   typedef walberla::boundary::BoundaryHandling< FlagField_T, Stencil, DiffusionDirichlet_T, SimpleDiffusionDirichlet_T, NoDiffusion_T, FreeDiffusion_T, SimpleDiffusionDirichlet_T, SimpleDiffusionDirichlet_T > BoundaryHandling_T;
+   typedef BoundaryHandling_T BoundaryHandling;
 
    const static FlagUID& getDiffusionDirichletFlagUID()        { static FlagUID uid( "DiffusionDirichlet"        ); return uid; }
    const static FlagUID& getSimpleDiffusionDirichletFlagUID()  { static FlagUID uid( "SimpleDiffusionDirichlet"  ); return uid; }
@@ -106,13 +101,13 @@ private:
       for( auto domainFlagUID = domainFlagUIDs.begin(); domainFlagUID != domainFlagUIDs.end(); ++domainFlagUID )
          field::addMask( domainMask, flagField->getOrRegisterFlag( *domainFlagUID ) );
 
-      BoundaryHandling_T * handling = new BoundaryHandling_T( "Diffusion Boundary Handling", flagField, domainMask, boost::tuples::make_tuple( 
+      BoundaryHandling_T * handling = new BoundaryHandling_T( "Diffusion Boundary Handling", flagField, domainMask,
          DiffusionDirichlet_T      ( getDiffusionDirichletBoundaryUID(),        getDiffusionDirichletFlagUID(),        pdfField, flagField ),
          SimpleDiffusionDirichlet_T( getSimpleDiffusionDirichletBoundaryUID(),  getSimpleDiffusionDirichletFlagUID(),  pdfField ),
          NoDiffusion_T             ( getNoDiffusionBoundaryUID(),               getNoDiffusionFlagUID(),               pdfField ),
          FreeDiffusion_T           ( getFreeDiffusionBoundaryUID(),             getFreeDiffusionFlagUID(),             pdfField, flagField, domainMask ),
          SimpleDiffusionDirichlet_T( getSimpleDiffusionDirichletBoundaryUID1(), getSimpleDiffusionDirichletFlagUID1(), pdfField ),
-         SimpleDiffusionDirichlet_T( getSimpleDiffusionDirichletBoundaryUID2(), getSimpleDiffusionDirichletFlagUID2(), pdfField )) );
+         SimpleDiffusionDirichlet_T( getSimpleDiffusionDirichletBoundaryUID2(), getSimpleDiffusionDirichletFlagUID2(), pdfField ) );
 
       if( initFlagUIDs.size() > size_t(0u) )
       {
diff --git a/src/lbm/boundary/factories/ExtendedBoundaryHandlingFactory.h b/src/lbm/boundary/factories/ExtendedBoundaryHandlingFactory.h
index 11ae8af9d3ca08da723e218d37b90796a6e35e85..4c08292b5638039717c3a5f2fa0f6c6e0eee9c68 100644
--- a/src/lbm/boundary/factories/ExtendedBoundaryHandlingFactory.h
+++ b/src/lbm/boundary/factories/ExtendedBoundaryHandlingFactory.h
@@ -33,8 +33,6 @@
 #include "boundary/BoundaryHandling.h"
 #include "lbm/field/PdfField.h"
 
-#include <boost/tuple/tuple.hpp>
-
 namespace walberla {
 namespace lbm {
 
@@ -86,8 +84,7 @@ public:
    typedef Outlet<LatticeModel, FlagFieldT >      BcOutlet;
    typedef Curved<LatticeModel, FlagFieldT >      BcCurved;
 
-   typedef boost::tuple< BcNoSlip, BcFreeSlip, BcPressure, BcUBB, BcOutlet, BcCurved >      BoundaryConditions;
-   typedef walberla::boundary::BoundaryHandling< FlagFieldT, Stencil, BoundaryConditions >  BoundaryHandling;
+   typedef walberla::boundary::BoundaryHandling< FlagFieldT, Stencil, BcNoSlip, BcFreeSlip, BcPressure, BcUBB, BcOutlet, BcCurved >  BoundaryHandling;
 
    static BlockDataID addBoundaryHandlingToStorage( const shared_ptr< StructuredBlockStorage > & bs, const std::string & identifier,
                                                     BlockDataID flagFieldID, BlockDataID pdfFieldID, const Set< FlagUID > & flagUIDSet )
@@ -156,14 +153,12 @@ ExtendedBoundaryHandlingFactory<LatticeModel, FlagFieldT>::operator()( IBlock *
 
 
    BoundaryHandling * handling = new BoundaryHandling( "extended lbm boundary handling", flagField, mask,
-      boost::tuples::make_tuple(
         BcNoSlip    ( getNoSlipBoundaryUID(),   getNoSlip(),   pdfField ),
         BcFreeSlip  ( getFreeSlipBoundaryUID(), getFreeSlip(), pdfField, flagField, mask ),
         BcPressure  ( getPressureBoundaryUID(), getPressure(), pdfField ),
         BcUBB       ( getUBBBoundaryUID(),      getUBB(),      pdfField, flagField, storage->getLevel(*block), block->getAABB() ),
         BcOutlet    ( getOutletBoundaryUID(),   getOutlet(),   pdfField, flagField, mask ),
         BcCurved    ( getCurvedBoundaryUID(),   getCurved(),   pdfField, flagField, mask )
-      )
     );
 
    return handling;
diff --git a/tests/boundary/BoundaryHandling.cpp b/tests/boundary/BoundaryHandling.cpp
index dbc83e1fe2ed0137d412fdcb8e6f8a96f6020505..77213c33d6fb45e04d93644c8451c38e0cd2b09a 100644
--- a/tests/boundary/BoundaryHandling.cpp
+++ b/tests/boundary/BoundaryHandling.cpp
@@ -35,8 +35,6 @@
 
 #include "stencil/D3Q27.h"
 
-#include <boost/tuple/tuple.hpp>
-
 
 namespace walberla {
 
@@ -203,7 +201,7 @@ private:
 // TEST BOUNDARY HANDLING //
 ////////////////////////////
 
-typedef BoundaryHandling< FlagField_T, stencil::D3Q27, boost::tuples::tuple< CopyBoundary, AddBoundary > > TestBoundaryHandling;
+typedef BoundaryHandling< FlagField_T, stencil::D3Q27, CopyBoundary, AddBoundary > TestBoundaryHandling;
 
 
 
@@ -368,7 +366,7 @@ static int main( int argc, char **argv )
    FactorField_T factorField_Ref( xSize, ySize, zSize, gl, uint_c(0) );
 
    TestBoundaryHandling handling( "test boundary handling", &flagField_BH, numeric_cast<flag_t>( domainFlag1 | domainFlag2 ),
-         boost::tuples::make_tuple( CopyBoundary( copy1, copy2, &workField_BH ), AddBoundary( add, &workField_BH ) ) );
+         CopyBoundary( copy1, copy2, &workField_BH ), AddBoundary( add, &workField_BH ) );
 
    // START TESTING / COMPARING
 
diff --git a/tests/lbm/BoundaryHandlingCommunication.cpp b/tests/lbm/BoundaryHandlingCommunication.cpp
index 781b7ca9cfe622768a27d0e9578981c2de7ce16a..6529e69134df4a6ea185864883f86ed4b7bfe9fb 100644
--- a/tests/lbm/BoundaryHandlingCommunication.cpp
+++ b/tests/lbm/BoundaryHandlingCommunication.cpp
@@ -85,9 +85,7 @@ public:
    typedef lbm::NoSlip< LatticeModel_T, flag_t >  NoSlip_T;
    typedef lbm::UBB< LatticeModel_T, flag_t >     UBB_T;
 
-   typedef boost::tuples::tuple< NoSlip_T, UBB_T > BoundaryConditions_T;
-
-   typedef BoundaryHandling< FlagField_T, typename LatticeModel_T::Stencil, BoundaryConditions_T > BoundaryHandling_T;
+   typedef BoundaryHandling< FlagField_T, typename LatticeModel_T::Stencil, NoSlip_T, UBB_T > BoundaryHandling_T;
 
    MyBoundaryHandling( const std::string & id, const BlockDataID & flagField, const BlockDataID & pdfField, const real_t velocity ) :
       id_( id ), flagField_( flagField ), pdfField_( pdfField ), velocity_( velocity ) {}
@@ -119,8 +117,8 @@ MyBoundaryHandling<LatticeModel_T>::operator()( IBlock * const block, const Stru
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * boundaryHandling = new BoundaryHandling_T( std::string("boundary handling ")+id_, flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
-                                       UBB_T( "velocity bounce back", UBB_Flag, pdfField ) ) );
+                                                                   NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
+                                                                   UBB_T( "velocity bounce back", UBB_Flag, pdfField ) );
 
    // closed no slip channel, periodic in x-direction
 
@@ -403,4 +401,4 @@ int main( int argc, char ** argv )
 int main( int argc, char ** argv )
 {
    return walberla::main(argc, argv);
-}
\ No newline at end of file
+}
diff --git a/tests/lbm/SuViscoelasticityTest.cpp b/tests/lbm/SuViscoelasticityTest.cpp
index 3952d7cb48a6faf2cc81c6d716c07bc5db54587a..b29b3678c712c2fe00578e4c310a463bba11dade 100644
--- a/tests/lbm/SuViscoelasticityTest.cpp
+++ b/tests/lbm/SuViscoelasticityTest.cpp
@@ -64,8 +64,7 @@ typedef FlagField< flag_t > FlagField_T;
 const uint_t FieldGhostLayers = 2;
 
 typedef lbm::NoSlip< LatticeModel_T, flag_t> NoSlip_T;
-typedef boost::tuples::tuple< NoSlip_T > BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T> BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T> BoundaryHandling_T;
 
 ///////////
 // FLAGS //
@@ -105,7 +104,7 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
-                                                           boost::tuples::make_tuple(    NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ) ) );
+                                                           NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ) );
 
    const auto noSlip = flagField->getFlag(NoSlip_Flag);
 
diff --git a/tests/lbm/SweepEquivalenceTest.cpp b/tests/lbm/SweepEquivalenceTest.cpp
index e28a00ebbc8295e9716dbea9653141eb0c6f93f4..091033dabda9f72e8cde20c887b95efe48a5b68a 100644
--- a/tests/lbm/SweepEquivalenceTest.cpp
+++ b/tests/lbm/SweepEquivalenceTest.cpp
@@ -94,9 +94,7 @@ public:
    typedef lbm::NoSlip< LatticeModel_T, flag_t >    NoSlip_T;
    typedef lbm::SimpleUBB< LatticeModel_T, flag_t > UBB_T;
 
-   typedef boost::tuples::tuple< NoSlip_T, UBB_T > BoundaryConditions_T;
-
-   typedef BoundaryHandling< FlagField_T, typename LatticeModel_T::Stencil, BoundaryConditions_T > BoundaryHandling_T;
+   typedef BoundaryHandling< FlagField_T, typename LatticeModel_T::Stencil, NoSlip_T, UBB_T > BoundaryHandling_T;
 
    MyBoundaryHandling( const std::string & id, const BlockDataID & flagField, const BlockDataID & pdfField, const real_t velocity ) :
       id_( id ), flagField_( flagField ), pdfField_( pdfField ), velocity_( velocity ) {}
@@ -129,8 +127,8 @@ MyBoundaryHandling<LatticeModel_T>::operator()( IBlock * const block, const Stru
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( std::string("boundary handling ")+id_, flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
-                                       UBB_T( "velocity bounce back", UBB_Flag, pdfField, velocity_, real_c(0), real_c(0) ) ) );
+                                                           NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
+                                                           UBB_T( "velocity bounce back", UBB_Flag, pdfField, velocity_, real_c(0), real_c(0) ) );
 
    // Couette flow -> periodic in x- and y-direction!
 
@@ -867,4 +865,4 @@ int main( int argc, char ** argv )
 int main( int argc, char* argv[] )
 {
   return walberla::main( argc, argv );
-}
\ No newline at end of file
+}
diff --git a/tests/lbm/boundary/BoundaryForce.cpp b/tests/lbm/boundary/BoundaryForce.cpp
index c36db2ec17579532697f6ea9897f372d5da1ecd0..cc1f9e63acf10569164968c64b55d24eb683c03c 100644
--- a/tests/lbm/boundary/BoundaryForce.cpp
+++ b/tests/lbm/boundary/BoundaryForce.cpp
@@ -68,8 +68,7 @@ using FlagField_T = FlagField<flag_t>;
 
 typedef lbm::SimpleUBB< LatticeModel_T, flag_t, false, true >  UBB_Sphere_T;
 typedef lbm::SimpleUBB< LatticeModel_T, flag_t >  UBB_Wall_T;
-typedef boost::tuples::tuple< UBB_Sphere_T, UBB_Wall_T >  BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, UBB_Sphere_T, UBB_Wall_T > BoundaryHandling_T;
 
 
 
@@ -113,8 +112,8 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block ) cons
    const flag_t fluid = flagField->registerFlag( Fluid_Flag ); // register the fluid flag at the flag field
    
    return new BoundaryHandling_T( "boundary handling", flagField, fluid,
-                                 boost::tuples::make_tuple( UBB_Sphere_T( "sphere", UBB_Sphere_Flag, pdfField, Vector3<real_t>() ),
-                                                            UBB_Wall_T(   "wall",   UBB_Wall_Flag,   pdfField, velocity_ ) ) );
+                                  UBB_Sphere_T( "sphere", UBB_Sphere_Flag, pdfField, Vector3<real_t>() ),
+                                  UBB_Wall_T(   "wall",   UBB_Wall_Flag,   pdfField, velocity_ ) );
 }
 
 
diff --git a/tests/lbm/boundary/BoundaryForceCouette.cpp b/tests/lbm/boundary/BoundaryForceCouette.cpp
index 974015b9d5c0d64e705694f28cff9182b182d4e1..01b2038feb3826ad96fa1538a7a0ef0876ef8cb9 100644
--- a/tests/lbm/boundary/BoundaryForceCouette.cpp
+++ b/tests/lbm/boundary/BoundaryForceCouette.cpp
@@ -67,8 +67,7 @@ using FlagField_T = FlagField<flag_t>;
 
 typedef lbm::NoSlip< LatticeModel_T, flag_t, true >  BottomWall_T;
 typedef lbm::SimpleUBB< LatticeModel_T, flag_t, false, true >  TopWall_T;
-typedef boost::tuples::tuple< BottomWall_T, TopWall_T >  BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, BottomWall_T, TopWall_T > BoundaryHandling_T;
 
 
 
@@ -112,8 +111,8 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block ) cons
    const flag_t fluid = flagField->registerFlag( Fluid_Flag ); // register the fluid flag at the flag field
    
    return new BoundaryHandling_T( "boundary handling", flagField, fluid,
-                                 boost::tuples::make_tuple( BottomWall_T( "bottom", BottomWall_Flag, pdfField ),
-                                                            TopWall_T( "top", TopWall_Flag, pdfField, velocity_ ) ) );
+                                  BottomWall_T( "bottom", BottomWall_Flag, pdfField ),
+                                  TopWall_T( "top", TopWall_Flag, pdfField, velocity_ ) );
 }
 
 
diff --git a/tests/lbm/boundary/SimplePABTest.cpp b/tests/lbm/boundary/SimplePABTest.cpp
index 6d5b8c478e11b09873636ee68040b24b91b204a6..a1de4f92ff59a1c018ef7501bb5de55a2231dc2d 100644
--- a/tests/lbm/boundary/SimplePABTest.cpp
+++ b/tests/lbm/boundary/SimplePABTest.cpp
@@ -108,8 +108,7 @@ class BoundaryHandlingCreator
 public:
    typedef lbm::NoSlip<LatticeModel,flag_t>                              MyNoSlip;
    typedef lbm::SimplePAB<LatticeModel,MyFlagField>                      MySimplePAB;
-   typedef boost::tuples::tuple< MyNoSlip, MySimplePAB, MySimplePAB >    BoundaryConditionsT;
-   typedef BoundaryHandling< MyFlagField, Stencil, BoundaryConditionsT > BoundaryHandlingT;
+   typedef BoundaryHandling< MyFlagField, Stencil, MyNoSlip, MySimplePAB, MySimplePAB > BoundaryHandlingT;
 
    BoundaryHandlingCreator( const BlockDataID& flagField, const BlockDataID& pdfField,
                             const real_t leftLatticeDensity, const real_t rightLatticeDensity,
@@ -148,11 +147,10 @@ BoundaryHandlingCreator::BoundaryHandlingT * BoundaryHandlingCreator::operator()
    */
 
    BoundaryHandlingT * handling = new BoundaryHandlingT( "Boundary Handling", flagField, fluidFlag,
-      boost::tuples::make_tuple(
          MyNoSlip( "NoSlip", getNoSlipFlag(), block->getData< PDFField >( pdfField_ ) ),
          MySimplePAB( "SimplePABLeft",  getSimplePABLeftFlag(),  block->getData< PDFField >( pdfField_ ), block->getData< MyFlagField >( flagField_ ), leftLatticeDensity_,  omega_, getFluidFlag(), getNoSlipFlag() ),
          MySimplePAB( "SimplePABRight", getSimplePABRightFlag(), block->getData< PDFField >( pdfField_ ), block->getData< MyFlagField >( flagField_ ), rightLatticeDensity_,  omega_, getFluidFlag(), getNoSlipFlag() )
-      ) );
+      );
 
    CellInterval channelBB(0, 0, 0,
       cell_idx_c(channelWidth_) - 1,  cell_idx_c(channelWidth_) - 1,  cell_idx_c(channelLength_) - 1 );
diff --git a/tests/lbm/refinement/CommunicationEquivalence.cpp b/tests/lbm/refinement/CommunicationEquivalence.cpp
index 150f7894947f95347c820d30b348779303fd5790..7124f5baf8305bd2b59bb442e1215e67da1c681e 100644
--- a/tests/lbm/refinement/CommunicationEquivalence.cpp
+++ b/tests/lbm/refinement/CommunicationEquivalence.cpp
@@ -84,9 +84,7 @@ const uint_t FieldGhostLayers = 4;
 typedef lbm::NoSlip< LatticeModel_T, flag_t >     NoSlip_T;
 typedef lbm::SimpleUBB< LatticeModel_T, flag_t >  UBB_T;
 
-typedef boost::tuples::tuple< NoSlip_T, UBB_T >  BoundaryConditions_T;
-
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, UBB_T > BoundaryHandling_T;
 
 ///////////
 // FLAGS //
@@ -239,8 +237,8 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block ) cons
    const flag_t fluid = flagField->registerFlag( Fluid_Flag );
 
    return new BoundaryHandling_T( "boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
-                                       UBB_T( "velocity bounce back", UBB_Flag, pdfField, topVelocity_, real_c(0), real_c(0) ) ) );
+                                  NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
+                                     UBB_T( "velocity bounce back", UBB_Flag, pdfField, topVelocity_, real_c(0), real_c(0) ) );
 }
 
 
@@ -492,4 +490,4 @@ int main( int argc, char ** argv )
 int main( int argc, char ** argv )
 {
    return walberla::main(argc, argv);
-}
\ No newline at end of file
+}
diff --git a/tests/lbm/refinement/Uniformity.cpp b/tests/lbm/refinement/Uniformity.cpp
index cd6eda7cd83f48dd72b8d8ab162e35668c3b451f..7ff59cbbb954c4444971e72ddcbfa1fae9162e9b 100644
--- a/tests/lbm/refinement/Uniformity.cpp
+++ b/tests/lbm/refinement/Uniformity.cpp
@@ -89,8 +89,7 @@ const uint_t FieldGhostLayers = 4;
 
 // dummy boundary handling
 typedef lbm::NoSlip< LatticeModel_T, flag_t > NoSlip_T;
-using BoundaryConditions_T = boost::tuples::tuple<NoSlip_T>;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T > BoundaryHandling_T;
 
 ///////////
 // FLAGS //
@@ -206,7 +205,7 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block ) cons
    flagField->setWithGhostLayer( fluid );
 
    return new BoundaryHandling_T( "boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "no slip", NoSlip_Flag, pdfField ) ) );
+                                  NoSlip_T( "no slip", NoSlip_Flag, pdfField ) );
 }
 
 
diff --git a/tests/pde/BoundaryTest.cpp b/tests/pde/BoundaryTest.cpp
index 226d597b84a9266d27e879255ef6d6397e587ad1..c232591ada427488f37cf8679e5949703396f5f1 100644
--- a/tests/pde/BoundaryTest.cpp
+++ b/tests/pde/BoundaryTest.cpp
@@ -62,9 +62,7 @@ typedef FlagField < flag_t >   FlagField_T;
 typedef pde::Dirichlet< Stencil_T, flag_t >  Dirichlet_T;
 typedef pde::Neumann< Stencil_T, flag_t >  Neumann_T;
 
-typedef boost::tuples::tuple< Dirichlet_T, Neumann_T >  BoundaryConditions_T;
-
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, Dirichlet_T, Neumann_T > BoundaryHandling_T;
 
 
 const FlagUID  Domain_Flag( "domain" );
@@ -118,9 +116,8 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block ) cons
    // A new boundary handling instance that uses the just fetched flag field is created:
    // Additional, internal flags used by the boundary handling will be stored in this flag field.
    return new BoundaryHandling_T( "boundary handling", flagField, domain,
-         boost::tuples::make_tuple( Dirichlet_T( "dirichlet", Dirichlet_Flag, rhsField, stencilField, adaptStencilField, flagField ) ,
+                                  Dirichlet_T( "dirichlet", Dirichlet_Flag, rhsField, stencilField, adaptStencilField, flagField ) ,
                                     Neumann_T( "neumann", Neumann_Flag, rhsField, stencilField, adaptStencilField, flagField, blockStorage_ )
-         )
    );
 }
 
diff --git a/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp b/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp
index 512f650d9baa06ed758174667a501dd00f2ff225..29b3d0ccab8083d85aba0077b28df6ed7a8c0c22 100644
--- a/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp
+++ b/tests/pe_coupling/discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp
@@ -88,8 +88,7 @@ using FlagField_T = FlagField<flag_t>;
 // boundary handling
 typedef lbm::NoSlip< LatticeModel_T, flag_t >                          NoSlip_T;
 
-using BoundaryConditions_T = boost::tuples::tuple<NoSlip_T>;
-typedef BoundaryHandling<FlagField_T, Stencil_T, BoundaryConditions_T> BoundaryHandling_T;
+typedef BoundaryHandling<FlagField_T, Stencil_T, NoSlip_T> BoundaryHandling_T;
 
 typedef std::tuple<pe::Plane, pe::Sphere> BodyTypeTuple ;
 
@@ -151,7 +150,7 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "Boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ) ) );
+                                                           NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ) );
 
    handling->fillWithDomain( FieldGhostLayers );
    return handling;
diff --git a/tests/pe_coupling/discrete_particle_methods/SphereWallCollisionBehaviorDPM.cpp b/tests/pe_coupling/discrete_particle_methods/SphereWallCollisionBehaviorDPM.cpp
index f1ae8d8c04f5b686129bfa5ade82c5d6bc50c85d..88b5950fdd88a6e0b34538507bb4a5b44cb8b2b5 100644
--- a/tests/pe_coupling/discrete_particle_methods/SphereWallCollisionBehaviorDPM.cpp
+++ b/tests/pe_coupling/discrete_particle_methods/SphereWallCollisionBehaviorDPM.cpp
@@ -87,8 +87,7 @@ using FlagField_T = FlagField<flag_t>;
 // boundary handling
 typedef lbm::NoSlip< LatticeModel_T, flag_t >                          NoSlip_T;
 
-using BoundaryConditions_T = boost::tuples::tuple<NoSlip_T>;
-typedef BoundaryHandling<FlagField_T, Stencil_T, BoundaryConditions_T> BoundaryHandling_T;
+typedef BoundaryHandling<FlagField_T, Stencil_T, NoSlip_T> BoundaryHandling_T;
 
 typedef std::tuple<pe::Plane, pe::Sphere> BodyTypeTuple ;
 
@@ -152,7 +151,7 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "Boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ) ) );
+                                                           NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ) );
 
    const auto noslip = flagField->getFlag( NoSlip_Flag );
 
diff --git a/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp b/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
index e70fe7f8e89f00298930db57617612056feee2fb..cfacb7cd3b58d3af9a12f4ab46047177fc54d1df 100644
--- a/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
@@ -92,8 +92,7 @@ typedef GhostLayerField< pe::BodyID, 1 >  BodyField_T;
 
 // boundary handling
 typedef pe_coupling::SimpleBB< LatticeModel_T, FlagField_T >  MO_T;
-using BoundaryConditions_T = boost::tuples::tuple<MO_T>;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, MO_T > BoundaryHandling_T;
 
 using BodyTypeTuple = std::tuple<pe::Sphere> ;
 
@@ -191,7 +190,7 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "fixed obstacle boundary handling", flagField, fluid,
-                                                           boost::tuples::make_tuple( MO_T (  "MO_BB",  MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+                                                           MO_T (  "MO_BB",  MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) );
 
    handling->fillWithDomain( FieldGhostLayers );
 
diff --git a/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp b/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp
index 59a224fe3fb1de87fad8b71c5f4f5c04d91d3cb2..0edda5942185ca62efd9b8231be0e79db741a8bf 100644
--- a/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/BodyMappingTest.cpp
@@ -74,8 +74,7 @@ typedef GhostLayerField< pe::BodyID, 1 >  BodyField_T;
 // boundary handling
 typedef lbm::NoSlip< LatticeModel_T, flag_t > NoSlip_T;
 typedef pe_coupling::SimpleBB< LatticeModel_T, FlagField_T >  MO_T;
-typedef boost::tuples::tuple< NoSlip_T, MO_T > BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, MO_T > BoundaryHandling_T;
 
 using BodyTypeTuple = std::tuple<pe::Sphere> ;
 
@@ -122,8 +121,8 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "fixed obstacle boundary handling", flagField, fluid,
-                                                           boost::tuples::make_tuple( NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
-                                                                                      MO_T (  "MO_BB",  MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+                                                           NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+                                                           MO_T (  "MO_BB",  MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) );
 
    handling->fillWithDomain( FieldGhostLayers );
 
diff --git a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
index 3d07ede20cb85170b3d6fb3550b8bef909c18941..2f6ac710af30d329c4f6f1b1b4790c30b3a7f55b 100644
--- a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp
@@ -92,8 +92,7 @@ typedef pe_coupling::SimpleBB       < LatticeModel_T, FlagField_T >  MO_BB_T;
 typedef pe_coupling::CurvedLinear   < LatticeModel_T, FlagField_T > MO_CLI_T;
 typedef pe_coupling::CurvedQuadratic< LatticeModel_T, FlagField_T >  MO_MR_T;
 
-typedef boost::tuples::tuple< MO_BB_T, MO_CLI_T, MO_MR_T >               BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, MO_BB_T, MO_CLI_T, MO_MR_T > BoundaryHandling_T;
 
 using BodyTypeTuple = std::tuple<pe::Sphere> ;
 
@@ -175,9 +174,9 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "fixed obstacle boundary handling", flagField, fluid,
-                                                           boost::tuples::make_tuple( MO_BB_T (  "MO_BB",  MO_BB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
-                                                                                     MO_CLI_T ( "MO_CLI", MO_CLI_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
-                                                                                      MO_MR_T (  "MO_MR",  MO_MR_Flag, pdfField, flagField, bodyField, fluid, *storage, *block, pdfFieldPreCol ) ) );
+                                                           MO_BB_T (  "MO_BB",  MO_BB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
+                                                           MO_CLI_T ( "MO_CLI", MO_CLI_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
+                                                           MO_MR_T (  "MO_MR",  MO_MR_Flag, pdfField, flagField, bodyField, fluid, *storage, *block, pdfFieldPreCol ) );
 
    handling->fillWithDomain( FieldGhostLayers );
 
diff --git a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
index ade3a8f05483ef4553125f1c1a91d08cea04f7fc..c60e11186ff6f378f71f73bcc0fbda8e64efdf3e 100644
--- a/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp
@@ -100,8 +100,7 @@ typedef GhostLayerField< pe::BodyID, 1 >  BodyField_T;
 typedef pe_coupling::SimpleBB       < LatticeModel_T, FlagField_T >  MO_BB_T;
 typedef pe_coupling::CurvedLinear   < LatticeModel_T, FlagField_T > MO_CLI_T;
 
-typedef boost::tuples::tuple< MO_BB_T, MO_CLI_T >                        BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, MO_BB_T, MO_CLI_T > BoundaryHandling_T;
 
 using BodyTypeTuple = std::tuple<pe::Sphere> ;
 
@@ -226,8 +225,8 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "fixed obstacle boundary handling", flagField, fluid,
-                                                           boost::tuples::make_tuple( MO_BB_T (  "MO_BB",  MO_BB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
-                                                                                     MO_CLI_T ( "MO_CLI", MO_CLI_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+                                                           MO_BB_T (  "MO_BB",  MO_BB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
+                                                           MO_CLI_T ( "MO_CLI", MO_CLI_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) );
 
    handling->fillWithDomain( FieldGhostLayers );
 
diff --git a/tests/pe_coupling/momentum_exchange_method/GlobalBodyAsBoundaryMEMStaticRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/GlobalBodyAsBoundaryMEMStaticRefinement.cpp
index a2b9a31056d80ef08eff3fa157d917a380f988c4..c19452458880578dc384d808a921a9278e45f842 100644
--- a/tests/pe_coupling/momentum_exchange_method/GlobalBodyAsBoundaryMEMStaticRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/GlobalBodyAsBoundaryMEMStaticRefinement.cpp
@@ -92,8 +92,7 @@ const uint_t FieldGhostLayers = 4;
 // boundary handling
 typedef pe_coupling::SimpleBB< LatticeModel_T, FlagField_T > MO_SBB_T;
 
-using BoundaryConditions_T = boost::tuples::tuple<MO_SBB_T>;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, MO_SBB_T > BoundaryHandling_T;
 
 using BodyTypeTuple = std::tuple<pe::Plane>;
 
@@ -224,7 +223,7 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( MO_SBB_T( "MO_SBB", MO_SBB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+          MO_SBB_T( "MO_SBB", MO_SBB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) );
 
    // boundary conditions are set by mapping the (moving) planes into the domain
 
diff --git a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
index db7d357fa7fe853213a3f21ac2d320ae2ba15a30..502cab1756749f734749c4e9b5de6cb00b22d711 100644
--- a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
@@ -95,9 +95,7 @@ typedef GhostLayerField< pe::BodyID, 1 >  BodyField_T;
 typedef lbm::FreeSlip< LatticeModel_T, FlagField_T>           FreeSlip_T;
 typedef pe_coupling::SimpleBB< LatticeModel_T, FlagField_T > MO_T;
 
-typedef boost::tuples::tuple< FreeSlip_T, MO_T > BoundaryConditions_T;
-
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, FreeSlip_T, MO_T > BoundaryHandling_T;
 
 typedef std::tuple<pe::Sphere, pe::Plane> BodyTypeTuple ;
 
@@ -461,8 +459,8 @@ BoundaryHandling_T * SphSphTestBoundaryHandling::operator()( IBlock * const bloc
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "cf boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( FreeSlip_T( "FreeSlip", FreeSlip_Flag, pdfField, flagField, fluid ),
-                                    MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+                                                           FreeSlip_T( "FreeSlip", FreeSlip_Flag, pdfField, flagField, fluid ),
+                                                           MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) );
 
    const auto freeslip = flagField->getFlag( FreeSlip_Flag );
 
@@ -548,8 +546,8 @@ BoundaryHandling_T * SphWallTestBoundaryHandling::operator()( IBlock * const blo
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "cf boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( FreeSlip_T( "FreeSlip", FreeSlip_Flag, pdfField, flagField, fluid ),
-                                    MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+                                                           FreeSlip_T( "FreeSlip", FreeSlip_Flag, pdfField, flagField, fluid ),
+                                                           MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) );
 
    const auto freeslip = flagField->getFlag( FreeSlip_Flag );
 
diff --git a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
index b130bb7d08f8abf6f26897b3928fd0f4a585894e..380a10f799e4e07c8be8bb8905ed2b7820091384 100644
--- a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
@@ -97,8 +97,7 @@ typedef lbm::NoSlip< LatticeModel_T, flag_t >                 NoSlip_T;
 typedef lbm::SimpleUBB< LatticeModel_T, flag_t >              UBB_T;
 typedef pe_coupling::SimpleBB< LatticeModel_T, FlagField_T > MO_T;
 
-typedef boost::tuples::tuple< NoSlip_T, UBB_T, MO_T > BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, UBB_T, MO_T > BoundaryHandling_T;
 
 typedef std::tuple< pe::Sphere, pe::Plane > BodyTypeTuple;
 
@@ -276,9 +275,9 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "cf boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+                                    NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
                                        UBB_T( "UBB", UBB_Flag, pdfField, velocity_, real_c(0), real_c(0) ),
-                                        MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+                                        MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) );
 
    const auto ubb = flagField->getFlag( UBB_Flag );
 
diff --git a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
index c074a75ee96925e4e083539e126b53f6c0b007d5..73a9d47cf29d4f9e34f6519330665ef5430cf5b1 100644
--- a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
@@ -99,8 +99,7 @@ typedef pe_coupling::SimpleBB       < LatticeModel_T, FlagField_T >  MO_BB_T;
 typedef pe_coupling::CurvedLinear   < LatticeModel_T, FlagField_T > MO_CLI_T;
 typedef pe_coupling::CurvedQuadratic< LatticeModel_T, FlagField_T >  MO_MR_T;
 
-typedef boost::tuples::tuple< NoSlip_T, MO_BB_T, MO_CLI_T, MO_MR_T >     BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, MO_BB_T, MO_CLI_T, MO_MR_T > BoundaryHandling_T;
 
 typedef std::tuple< pe::Sphere, pe::Plane > BodyTypeTuple;
 
@@ -189,10 +188,10 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
-         boost::tuples::make_tuple(    NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+                                       NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
                                        MO_BB_T (  "MO_BB",  MO_BB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
                                        MO_CLI_T( "MO_CLI", MO_CLI_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
-                                       MO_MR_T (  "MO_MR",  MO_MR_Flag, pdfField, flagField, bodyField, fluid, *storage, *block, pdfFieldPreCol ) ) );
+                                       MO_MR_T (  "MO_MR",  MO_MR_Flag, pdfField, flagField, bodyField, fluid, *storage, *block, pdfFieldPreCol ) );
 
    const auto noslip = flagField->getFlag( NoSlip_Flag );
 
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
index c1cc240d498ce6e3ed28d15646107043af8604ad..b8c0b6a0bc0db5bc736bcdc12243a738855169de 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEM.cpp
@@ -95,8 +95,7 @@ typedef lbm::NoSlip< LatticeModel_T, flag_t > NoSlip_T;
 
 typedef pe_coupling::CurvedLinear< LatticeModel_T, FlagField_T > MO_T;
 
-typedef boost::tuples::tuple< NoSlip_T, MO_T > BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, MO_T > BoundaryHandling_T;
 
 typedef std::tuple< pe::Sphere, pe::Plane > BodyTypeTuple;
 
@@ -142,8 +141,8 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
-         boost::tuples::make_tuple(    NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
-                                       MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+                                                           NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+                                                           MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) );
 
    // boundary conditions (no-slip) are set by mapping the planes into the domain
 
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
index a8888df59ef3253447242070939982078bd70cef..9cecbe8b7b631be2c335a43b3a237f8209b0481b 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMDynamicRefinement.cpp
@@ -101,8 +101,7 @@ typedef lbm::NoSlip< LatticeModel_T, flag_t > NoSlip_T;
 
 typedef pe_coupling::CurvedLinear< LatticeModel_T, FlagField_T > MO_T;
 
-typedef boost::tuples::tuple< NoSlip_T, MO_T > BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, MO_T > BoundaryHandling_T;
 
 typedef std::tuple< pe::Sphere, pe::Plane > BodyTypeTuple;
 
@@ -244,8 +243,8 @@ BoundaryHandling_T * MyBoundaryHandling::initialize( IBlock * const block )
    WALBERLA_CHECK_NOT_NULLPTR( blocksPtr );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
-                                                           boost::tuples::make_tuple( NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
-                                                                                      MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *blocksPtr, *block ) ) );
+                                                           NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+                                                           MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *blocksPtr, *block ) );
 
    handling->fillWithDomain( FieldGhostLayers );
 
diff --git a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
index 9e1407c26709a32c455559c12b45bc645a866226..6ba596359a6934ec5292e8395bd644a704267455 100644
--- a/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SettlingSphereMEMStaticRefinement.cpp
@@ -97,8 +97,7 @@ typedef lbm::NoSlip< LatticeModel_T, flag_t > NoSlip_T;
 
 typedef pe_coupling::CurvedLinear< LatticeModel_T, FlagField_T > MO_T;
 
-typedef boost::tuples::tuple< NoSlip_T, MO_T > BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, MO_T > BoundaryHandling_T;
 
 typedef std::tuple< pe::Sphere, pe::Plane > BodyTypeTuple;
 
@@ -230,8 +229,8 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
-         boost::tuples::make_tuple(    NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
-                                       MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+                                                           NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+                                                           MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) );
 
    // boundary conditions (no-slip) are set by mapping the planes into the domain
 
diff --git a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
index 0c8a6ccac228d5ce94d4d9fff0eff338bd5432fc..c13ca3550252c24e951b65c2ba83e42161cb9c09 100644
--- a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
@@ -94,8 +94,7 @@ const uint_t FieldGhostLayers = 1;
 // boundary handling
 typedef pe_coupling::SimpleBB<LatticeModel_T, FlagField_T> MO_BB_T;
 
-using BoundaryConditions_T = boost::tuples::tuple<MO_BB_T>;
-typedef BoundaryHandling<FlagField_T, Stencil_T, BoundaryConditions_T> BoundaryHandling_T;
+typedef BoundaryHandling<FlagField_T, Stencil_T, MO_BB_T> BoundaryHandling_T;
 
 using BodyTypeTuple = std::tuple<pe::Squirmer>;
 
@@ -143,10 +142,8 @@ MyBoundaryHandling::operator()(IBlock *const block, const StructuredBlockStorage
    const auto fluid = flagField->flagExists(Fluid_Flag) ? flagField->getFlag(Fluid_Flag) : flagField->registerFlag(
          Fluid_Flag);
 
-   BoundaryHandling_T *handling = new BoundaryHandling_T("fixed obstacle boundary handling", flagField, fluid,
-                                                         boost::tuples::make_tuple(
-                                                               MO_BB_T("MO_BB", MO_BB_Flag, pdfField, flagField,
-                                                                       bodyField, fluid, *storage, *block)));
+   BoundaryHandling_T *handling = new BoundaryHandling_T( "fixed obstacle boundary handling", flagField, fluid,
+                                                          MO_BB_T("MO_BB", MO_BB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) );
 
    handling->fillWithDomain(FieldGhostLayers);
 
@@ -491,4 +488,4 @@ int main(int argc, char **argv) {
 
 int main( int argc, char **argv ){
    squirmer::main(argc, argv);
-}
\ No newline at end of file
+}
diff --git a/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp b/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
index 7f6e659236f8b76053af10a1ea419677c2a96b2f..67ce8c233fcc0d50a0c579fcf7d76255f9f071b0 100644
--- a/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
@@ -86,8 +86,7 @@ const uint_t FieldGhostLayers = 1;
 
 typedef pe_coupling::CurvedLinear< LatticeModel_T, FlagField_T > MO_T;
 
-using BoundaryConditions_T = boost::tuples::tuple<MO_T>;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, MO_T > BoundaryHandling_T;
 
 typedef std::tuple< pe::Capsule, pe::CylindricalBoundary > BodyTypeTuple;
 
@@ -200,7 +199,7 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "cf boundary handling", flagField, fluid,
-         boost::tuples::make_tuple( MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+                                                           MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) );
 
    CellInterval domainBB = storage->getDomainCellBB();
    domainBB.xMin() -= cell_idx_c( FieldGhostLayers );
diff --git a/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp b/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
index c76d112cd9020919f8f95b195c5d19146866c349..0b649234430aa2820acae8544694316445f0590a 100644
--- a/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp
@@ -92,8 +92,7 @@ typedef pe_coupling::SimpleBB       < LatticeModel_T, FlagField_T >  MO_BB_T;
 typedef pe_coupling::CurvedLinear   < LatticeModel_T, FlagField_T > MO_CLI_T;
 typedef pe_coupling::CurvedQuadratic< LatticeModel_T, FlagField_T >  MO_MR_T;
 
-typedef boost::tuples::tuple< MO_BB_T, MO_CLI_T, MO_MR_T >               BoundaryConditions_T;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, MO_BB_T, MO_CLI_T, MO_MR_T > BoundaryHandling_T;
 
 using BodyTypeTuple = std::tuple<pe::Sphere> ;
 
@@ -158,9 +157,9 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "fixed obstacle boundary handling", flagField, fluid,
-                                                           boost::tuples::make_tuple( MO_BB_T (  "MO_BB",  MO_BB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
-                                                                                     MO_CLI_T ( "MO_CLI", MO_CLI_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
-                                                                                      MO_MR_T (  "MO_MR",  MO_MR_Flag, pdfField, flagField, bodyField, fluid, *storage, *block, pdfFieldPreCol ) ) );
+                                                            MO_BB_T (  "MO_BB",  MO_BB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
+                                                           MO_CLI_T ( "MO_CLI", MO_CLI_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
+                                                            MO_MR_T (  "MO_MR",  MO_MR_Flag, pdfField, flagField, bodyField, fluid, *storage, *block, pdfFieldPreCol ) );
 
    handling->fillWithDomain( FieldGhostLayers );
 
diff --git a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
index 092c0ea36edee47a9ec298c80c638f3a44e93694..b35c41db12a57bcea104e4b6fa81213ca0ca3f14 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/DragForceSpherePSMRefinement.cpp
@@ -99,8 +99,7 @@ typedef std::pair< pe::BodyID, real_t >                              BodyAndVolu
 typedef GhostLayerField< std::vector< BodyAndVolumeFraction_T >, 1 > BodyAndVolumeFractionField_T;
 
 typedef lbm::NoSlip< LatticeModel_T, flag_t > NoSlip_T;
-using BoundaryConditions_T = boost::tuples::tuple<NoSlip_T>;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T > BoundaryHandling_T;
 
 using BodyTypeTuple = std::tuple<pe::Sphere> ;
 
@@ -243,7 +242,7 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "fixed obstacle boundary handling", flagField, fluid,
-                                                           boost::tuples::make_tuple( NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ) ) );
+                                                           NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ) );
 
    handling->fillWithDomain( FieldGhostLayers );
 
diff --git a/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp b/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
index deb8f68963cca2bfa801bfa69e7456cf9aec3f16..fbec373004928d25d7b1d548417e32a2d58f221b 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
@@ -97,8 +97,7 @@ typedef GhostLayerField< std::vector< BodyAndVolumeFraction_T >, 1 > BodyAndVolu
 
 // boundary handling
 typedef lbm::NoSlip< LatticeModel_T, flag_t > NoSlip_T;
-using BoundaryConditions_T = boost::tuples::tuple<NoSlip_T>;
-typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T > BoundaryHandling_T;
 
 typedef std::tuple< pe::Sphere, pe::Plane > BodyTypeTuple;
 
@@ -177,7 +176,7 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
    const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
    BoundaryHandling_T * handling = new BoundaryHandling_T( "boundary handling", flagField, fluid,
-         boost::tuples::make_tuple(    NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ) ) );
+                                                           NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ) );
 
    const auto noslip = flagField->getFlag( NoSlip_Flag );