diff --git a/src/mesh/blockforest/BlockForestInitialization.h b/src/mesh/blockforest/BlockForestInitialization.h
index 6f9018beb5617cbb98718e8c285a529eac9ede9e..1c39d1c2ec8c340e53647e36a20577c131fdbcd4 100644
--- a/src/mesh/blockforest/BlockForestInitialization.h
+++ b/src/mesh/blockforest/BlockForestInitialization.h
@@ -21,6 +21,10 @@
 
 #pragma once
 
+#include "BlockExclusion.h"
+#include "mesh/MeshOperations.h"
+#include "mesh/distance_octree/DistanceOctree.h"
+
 #include "blockforest/StructuredBlockForest.h"
 #include "blockforest/SetupBlockForest.h"
 #include "blockforest/loadbalancing/StaticCurve.h"
@@ -110,5 +114,37 @@ private:
 };
 
 
+template< typename MeshType >
+shared_ptr<StructuredBlockForest> createStructuredBlockStorageInsideMesh( const shared_ptr< mesh::DistanceOctree< MeshType > > & distanceOctree, const real_t dx, const Vector3<uint_t> & blockSize )
+{
+   ComplexGeometryStructuredBlockforestCreator creator( distanceOctree->getAABB(), Vector3<real_t>(dx), makeExcludeMeshExterior( distanceOctree, dx ) );
+   return creator.createStructuredBlockForest( blockSize );
+}
+
+
+template< typename MeshType >
+shared_ptr<StructuredBlockForest> createStructuredBlockStorageOutsideMesh( const AABB & aabb, const shared_ptr< mesh::DistanceOctree< MeshType > > & distanceOctree, const real_t dx, const Vector3<uint_t> & blockSize )
+{
+   ComplexGeometryStructuredBlockforestCreator creator( aabb, Vector3<real_t>(dx), makeExcludeMeshInterior( distanceOctree, dx ) );
+   return creator.createStructuredBlockForest( blockSize );
+}
+
+
+template< typename MeshType >
+shared_ptr<StructuredBlockForest> createStructuredBlockStorageInsideMesh( const shared_ptr< mesh::DistanceOctree< MeshType > > & distanceOctree, const real_t dx, const uint_t targetNumRootBlocks )
+{
+   ComplexGeometryStructuredBlockforestCreator creator( distanceOctree->getAABB(), Vector3<real_t>(dx), makeExcludeMeshExterior( distanceOctree, dx ) );
+   return creator.createStructuredBlockForest( targetNumRootBlocks );
+}
+
+
+template< typename MeshType >
+shared_ptr<StructuredBlockForest> createStructuredBlockStorageOutsideMesh( const AABB & aabb, const shared_ptr< mesh::DistanceOctree< MeshType > > & distanceOctree, const real_t dx, const uint_t targetNumRootBlocks )
+{
+   ComplexGeometryStructuredBlockforestCreator creator( aabb, Vector3<real_t>(dx), makeExcludeMeshInterior( distanceOctree, dx ) );
+   return creator.createStructuredBlockForest( targetNumRootBlocks );
+}
+
+
 } // namespace mesh
 } // namespace walberla
\ No newline at end of file
diff --git a/src/mesh/boundary/BoundarySetup.cpp b/src/mesh/boundary/BoundarySetup.cpp
index b8390b6b71ecf287821ab4f2694f5b6fcc628aea..f0094a972fd390c1fb8cf1dc96a8e9dbb85c0a70 100644
--- a/src/mesh/boundary/BoundarySetup.cpp
+++ b/src/mesh/boundary/BoundarySetup.cpp
@@ -66,9 +66,9 @@ void BoundarySetup::allocateOrResetVoxelizationField()
 {
    if( voxelizationFieldId_ )
    {
-      for( auto blockIt = structuredBlockStorage_->begin(); blockIt != structuredBlockStorage_->end(); ++blockIt )
+      for( auto & block : *structuredBlockStorage_ )
       {
-         VoxelizationField * voxelizationField = blockIt->getData< VoxelizationField >( *voxelizationFieldId_ );
+         VoxelizationField * voxelizationField = block.getData< VoxelizationField >( *voxelizationFieldId_ );
          voxelizationField->setWithGhostLayer( uint8_t(0) );
       }
    }
@@ -93,15 +93,15 @@ void BoundarySetup::voxelize()
 {
    allocateOrResetVoxelizationField();
 
-   for( auto blockIt = structuredBlockStorage_->begin(); blockIt != structuredBlockStorage_->end(); ++blockIt )
+   for( auto & block : *structuredBlockStorage_ )
    {
-      VoxelizationField * voxelizationField = blockIt->getData< VoxelizationField >( *voxelizationFieldId_ );
+      VoxelizationField * voxelizationField = block.getData< VoxelizationField >( *voxelizationFieldId_ );
 
       WALBERLA_ASSERT_NOT_NULLPTR( voxelizationField );
       WALBERLA_ASSERT_EQUAL( numGhostLayers_, voxelizationField->nrOfGhostLayers() );
 
       CellInterval blockCi = voxelizationField->xyzSizeWithGhostLayer();
-      structuredBlockStorage_->transformBlockLocalToGlobalCellInterval( blockCi, *blockIt );
+      structuredBlockStorage_->transformBlockLocalToGlobalCellInterval( blockCi, block );
 
       std::queue< CellInterval > ciQueue;
       ciQueue.push( blockCi );
@@ -112,7 +112,7 @@ void BoundarySetup::voxelize()
 
          WALBERLA_ASSERT( !curCi.empty(), "Cell Interval: " << curCi );
 
-         AABB curAABB = structuredBlockStorage_->getAABBFromCellBB( curCi, structuredBlockStorage_->getLevel( *blockIt ) );
+         AABB curAABB = structuredBlockStorage_->getAABBFromCellBB( curCi, structuredBlockStorage_->getLevel( block ) );
 
          WALBERLA_ASSERT( !curAABB.empty(), "AABB: " << curAABB );
 
@@ -125,7 +125,7 @@ void BoundarySetup::voxelize()
             if( ( sqSignedDistance < real_t(0) ) )
             {
                Cell localCell;
-               structuredBlockStorage_->transformGlobalToBlockLocalCell( localCell, *blockIt, curCi.min() );
+               structuredBlockStorage_->transformGlobalToBlockLocalCell( localCell, block, curCi.min() );
                voxelizationField->get( localCell ) = uint8_t(1);
             }
 
@@ -140,7 +140,7 @@ void BoundarySetup::voxelize()
          {
             // clearly the cell interval is fully covered by the mesh
             CellInterval localCi;
-            structuredBlockStorage_->transformGlobalToBlockLocalCellInterval( localCi, *blockIt, curCi );
+            structuredBlockStorage_->transformGlobalToBlockLocalCellInterval( localCi, block, curCi );
             std::fill( voxelizationField->beginSliceXYZ( localCi ), voxelizationField->end(), uint8_t(1) );
 
             ciQueue.pop();
diff --git a/src/mesh/boundary/BoundarySetup.h b/src/mesh/boundary/BoundarySetup.h
index 08840ed0ae24502e7248d9831351b955d15f0d30..2a54d688e65ba10b9b658fe0cdf24d649d4b714c 100644
--- a/src/mesh/boundary/BoundarySetup.h
+++ b/src/mesh/boundary/BoundarySetup.h
@@ -82,9 +82,9 @@ private:
 
    shared_ptr< StructuredBlockStorage >       structuredBlockStorage_;
    shared_ptr< BlockDataID >                  voxelizationFieldId_;
-   DistanceFunction                           distanceFunction_;  // function providing the squared signed distance to an object
+   DistanceFunction                           distanceFunction_;  /// function providing the squared signed distance to an object
    uint_t                                     numGhostLayers_;
-   size_t                                     cellVectorChunkSize_;
+   size_t                                     cellVectorChunkSize_; /// Number of boundary cells which are setup simultaneously 
 };
 
 
@@ -92,12 +92,12 @@ private:
 template< typename BoundaryHandlingType >
 void BoundarySetup::setDomainCells( const BlockDataID boundaryHandlingId, const Location domainLocation )
 {
-   for( auto blockIt = structuredBlockStorage_->begin(); blockIt != structuredBlockStorage_->end(); ++blockIt )
+   for( auto & block : *structuredBlockStorage_ )
    {
-      BoundaryHandlingType * boundaryHandling  = blockIt->getData< BoundaryHandlingType >( boundaryHandlingId  );
+      BoundaryHandlingType * boundaryHandling  = block.getData< BoundaryHandlingType >( boundaryHandlingId  );
       WALBERLA_CHECK_NOT_NULLPTR( boundaryHandling, "boundaryHandlingId invalid!" );
 
-      const VoxelizationField * voxelizationField = blockIt->getData< VoxelizationField >( *voxelizationFieldId_ );
+      const VoxelizationField * voxelizationField = block.getData< VoxelizationField >( *voxelizationFieldId_ );
       WALBERLA_ASSERT_NOT_NULLPTR( voxelizationField );
       WALBERLA_CHECK_LESS_EQUAL( numGhostLayers_, boundaryHandling->getFlagField()->nrOfGhostLayers(), "You want to use mesh boundary setup with " \
                                  << numGhostLayers_ << " but your flag field has only " << boundaryHandling->getFlagField()->nrOfGhostLayers() << " ghost layers!" );
@@ -111,7 +111,7 @@ void BoundarySetup::setDomainCells( const BlockDataID boundaryHandlingId, const
             for( cell_idx_t x = -cell_idx_c( numGhostLayers_ ); x < cell_idx_c( voxelizationField->xSize() + numGhostLayers_ ); ++x )
             {
                if( voxelizationField->get( x, y, z ) == domainValue )
-                  domainCells.push_back( Cell(x,y,z) );
+                  domainCells.emplace_back( x, y, z );
 
                if( domainCells.size() > cellVectorChunkSize_ )
                {
@@ -128,13 +128,13 @@ void BoundarySetup::setDomainCells( const BlockDataID boundaryHandlingId, const
 template<typename FlagField_T>
 void BoundarySetup::setFlag( const BlockDataID flagFieldID, field::FlagUID flagUID, Location boundaryLocation )
 {
-   for( auto blockIt = structuredBlockStorage_->begin(); blockIt != structuredBlockStorage_->end(); ++blockIt )
+   for( auto & block : *structuredBlockStorage_ )
    {
-      FlagField_T * flagField  = blockIt->getData< FlagField_T >( flagFieldID );
+      FlagField_T * flagField  = block.getData< FlagField_T >( flagFieldID );
       WALBERLA_CHECK_NOT_NULLPTR( flagField, "flagFieldID invalid!" );
       auto flag = flagField->getFlag(flagUID);
 
-      const VoxelizationField * voxelizationField = blockIt->getData< VoxelizationField >( *voxelizationFieldId_ );
+      const VoxelizationField * voxelizationField = block.getData< VoxelizationField >( *voxelizationFieldId_ );
       WALBERLA_ASSERT_NOT_NULLPTR( voxelizationField );
       WALBERLA_CHECK_LESS_EQUAL( numGhostLayers_, flagField->nrOfGhostLayers(), "You want to use mesh boundary setup with " \
                                  << numGhostLayers_ << " but your flag field has only " << flagField->nrOfGhostLayers() << " ghost layers!" );
@@ -155,12 +155,12 @@ void BoundarySetup::setFlag( const BlockDataID flagFieldID, field::FlagUID flagU
 template< typename BoundaryHandlingType, typename BoundaryFunction, typename Stencil >
 void BoundarySetup::setBoundaries( const BlockDataID boundaryHandlingID, const BoundaryFunction & boundaryFunction, Location boundaryLocation )
 {
-   for( auto blockIt = structuredBlockStorage_->begin(); blockIt != structuredBlockStorage_->end(); ++blockIt )
+   for( auto & block : *structuredBlockStorage_ )
    {
-      BoundaryHandlingType * boundaryHandling  = blockIt->getData< BoundaryHandlingType >( boundaryHandlingID  );
+      BoundaryHandlingType * boundaryHandling  = block.getData< BoundaryHandlingType >( boundaryHandlingID  );
       WALBERLA_CHECK_NOT_NULLPTR( boundaryHandling, "boundaryHandlingId invalid!" );
 
-      const VoxelizationField * voxelizationField = blockIt->getData< VoxelizationField    >( *voxelizationFieldId_ );
+      const VoxelizationField * voxelizationField = block.getData< VoxelizationField >( *voxelizationFieldId_ );
       WALBERLA_ASSERT_NOT_NULLPTR( voxelizationField );
       WALBERLA_CHECK_LESS_EQUAL( numGhostLayers_, boundaryHandling->getFlagField()->nrOfGhostLayers(), "You want to use mesh boundary setup with " \
                                  << numGhostLayers_ << " but your flag field has only " << boundaryHandling->getFlagField()->nrOfGhostLayers() << " ghost layers!" );
@@ -182,7 +182,7 @@ void BoundarySetup::setBoundaries( const BlockDataID boundaryHandlingID, const B
                   const Cell neighborCell = cell + *dirIt;
                   if( blockCi.contains( neighborCell ) && voxelizationField->get( neighborCell ) == domainValue )
                   {
-                     const Vector3< real_t > cellCenter = structuredBlockStorage_->getBlockLocalCellCenter( *blockIt, cell );
+                     const Vector3< real_t > cellCenter = structuredBlockStorage_->getBlockLocalCellCenter( block, cell );
                      const BoundaryInfo & bi = boundaryFunction( cellCenter );
                      const auto boundaryMask = boundaryHandling->getBoundaryMask( bi.getUid() );
 
diff --git a/tests/mesh/MeshInitilizationTest.cpp b/tests/mesh/MeshInitilizationTest.cpp
index cc98a3e5150eedb4a57374cee41d6ec88fae8827..494f212fb68abf7c99e626196168fa156f542092 100644
--- a/tests/mesh/MeshInitilizationTest.cpp
+++ b/tests/mesh/MeshInitilizationTest.cpp
@@ -102,7 +102,7 @@ void test( const std::string & meshFile, const uint_t numProcesses, const uint_t
    WALBERLA_LOG_INFO_ON_ROOT( "Creating SBF with StaticLevelwiseCurveBalanceWeighted Partitioner" );
    bfc.setTargetProcessAssignmentFunction( blockforest::StaticLevelwiseCurveBalanceWeighted() );
    auto sbf_default = bfc.createSetupBlockForest( Vector3<uint_t>(64,64,64), numProcesses );
-   sbf_default->writeVTKOutput("sbf_default");
+   //sbf_default->writeVTKOutput("sbf_default");
    WALBERLA_LOG_INFO_ON_ROOT( sbf_default->toString() );
 
    return;
@@ -137,6 +137,37 @@ void test( const std::string & meshFile, const uint_t numProcesses, const uint_t
 #endif
 }
 
+template< typename MeshType >
+void testHelperFunctions( const std::string & meshFile, const uint_t numTotalBlocks )
+{
+   auto mesh = make_shared<MeshType>();
+   mesh::readAndBroadcast( meshFile, *mesh);
+   auto triDist = make_shared< mesh::TriangleDistance<MeshType> >( mesh );
+   auto distanceOctree = make_shared< DistanceOctree< MeshType > >( triDist );
+
+   const real_t meshVolume  = real_c( computeVolume( *mesh ) );
+   const real_t blockVolume = meshVolume / real_c( numTotalBlocks );
+   static const real_t cellsPersBlock = real_t(1000);
+   const real_t cellVolume = blockVolume / cellsPersBlock;
+   const real_t dx = std::pow( cellVolume, real_t(1) / real_t(3) );
+
+   WALBERLA_LOG_INFO_ON_ROOT( "Creating SBF with createStructuredBlockStorageInsideMesh with block size" );
+   auto sbf0 = mesh::createStructuredBlockStorageInsideMesh( distanceOctree, dx, numTotalBlocks );
+   
+   WALBERLA_LOG_INFO_ON_ROOT( "Creating SBF with createStructuredBlockStorageInsideMesh with block size" );
+   Vector3<uint_t> blockSize( sbf0->getNumberOfXCells(), sbf0->getNumberOfYCells(), sbf0->getNumberOfZCells() );
+   auto sbf1 = mesh::createStructuredBlockStorageInsideMesh( distanceOctree, dx, blockSize );
+
+   auto exteriorAabb = computeAABB( *mesh ).getScaled( real_t(2) );
+
+   WALBERLA_LOG_INFO_ON_ROOT( "Creating SBF with createStructuredBlockStorageInsideMesh with block size" );
+   auto sbf2 = mesh::createStructuredBlockStorageOutsideMesh( exteriorAabb, distanceOctree, dx, numTotalBlocks );
+
+   WALBERLA_LOG_INFO_ON_ROOT( "Creating SBF with createStructuredBlockStorageInsideMesh with block size" );
+   blockSize = Vector3<uint_t>( sbf2->getNumberOfXCells(), sbf2->getNumberOfYCells(), sbf2->getNumberOfZCells() );
+   auto sbf3 = mesh::createStructuredBlockStorageOutsideMesh( exteriorAabb, distanceOctree, dx, blockSize );
+}
+
 int main( int argc, char * argv[] )
 {
    debug::enterTestMode();
@@ -155,6 +186,8 @@ int main( int argc, char * argv[] )
    //test< mesh::FloatTriangleMesh >( meshFile, numProcesses, numTotalBlocks );
    //test< mesh::PythonTriangleMesh >( meshFile, numProcesses, numTotalBlocks );
 
+   testHelperFunctions< mesh::TriangleMesh >( meshFile, numTotalBlocks );
+
    return EXIT_SUCCESS;
 }