Merge branch 'master' of i10git.cs.fau.de:walberla/walberla

CMakeLists.txt

@@ -1027,31 +1027,22 @@ endif()
 option ( WALBERLA_THREAD_SAFE_LOGGING "Enables/Disables thread-safe logging" ON )
 
 if ( WALBERLA_BUILD_WITH_OPENMP )
-    if ( WALBERLA_CXX_COMPILER_IS_INTEL )
+    if ( WALBERLA_CXX_COMPILER_IS_INTEL AND "${CMAKE_CXX_COMPILER_VERSION}" VERSION_LESS "16.0.3" )
        add_flag ( CMAKE_C_FLAGS   "-openmp" )
        add_flag ( CMAKE_CXX_FLAGS "-openmp" )
-    elseif ( CMAKE_COMPILER_IS_GNUCXX )
-      add_flag ( CMAKE_C_FLAGS   "-fopenmp" )
-      add_flag ( CMAKE_CXX_FLAGS "-fopenmp" )
-    elseif ( WALBERLA_CXX_COMPILER_IS_CLANG )
-       add_flag ( CMAKE_C_FLAGS   "-fopenmp" )
-       add_flag ( CMAKE_CXX_FLAGS "-fopenmp" )
-    elseif ( WALBERLA_CXX_COMPILER_IS_MSVC )
-      add_flag ( CMAKE_C_FLAGS   "/openmp" )
-      add_flag ( CMAKE_CXX_FLAGS "/openmp" )
-    elseif ( WALBERLA_CXX_COMPILER_IS_IBM )
-      add_flag ( CMAKE_C_FLAGS   "-qsmp=omp" )
-      add_flag ( CMAKE_CXX_FLAGS "-qsmp=omp" )
-      # There has been an internal compiler error with the IBM compiler, so WALBERLA_THREAD_SAFE_LOGGING is disabled by default for this compiler
-      set ( WALBERLA_THREAD_SAFE_LOGGING OFF CACHE BOOL "Enables/Disables thread-safe logging" FORCE )
     elseif ( WALBERLA_CXX_COMPILER_IS_NEC )
-      add_flag ( CMAKE_C_FLAGS   "-Popenmp" )
-      add_flag ( CMAKE_CXX_FLAGS "-Popenmp" )
+       add_flag ( CMAKE_C_FLAGS   "-Popenmp" )
+       add_flag ( CMAKE_CXX_FLAGS "-Popenmp" )
+    else()
+       find_package( OpenMP )
+       add_flag ( CMAKE_C_FLAGS   "${OpenMP_C_FLAGS}" )
+       add_flag ( CMAKE_CXX_FLAGS "${OpenMP_CXX_FLAGS}" )
+       list ( APPEND SERVICE_LIBS ${OpenMP_CXX_LIBRARIES} )
     endif()
 else()
     if ( WALBERLA_CXX_COMPILER_IS_CRAY )
-      add_flag ( CMAKE_C_FLAGS   "-h noomp" )
-      add_flag ( CMAKE_CXX_FLAGS "-h noomp" )
+       add_flag ( CMAKE_C_FLAGS   "-h noomp" )
+       add_flag ( CMAKE_CXX_FLAGS "-h noomp" )
     endif()
 endif()
 ############################################################################################################################

src/blockforest/loadbalancing/DynamicCurve.h

@@ -792,8 +792,9 @@ void DynamicCurveBalance< PhantomData_T >::balanceWeighted( const std::vector< s
          long double weight( 0 );
          int numBlocks( 0 );
          while( c < blocks.size() &&
-                std::abs( pWeight - weight - numeric_cast< long double >( allBlocks[ uint_c( blocks[c].first ) ][ blocks[c].second ].second ) ) <=
-                std::abs( pWeight - weight ) &&
+                ( isIdentical(weight, 0.0l) || 
+                  std::abs( pWeight - weight - numeric_cast< long double >( allBlocks[ uint_c( blocks[c].first ) ][ blocks[c].second ].second ) ) <=
+                  std::abs( pWeight - weight ) ) &&
                 numBlocks < maxBlocksPerProcess_ )
          {
             targets[ uint_c( blocks[c].first ) ][ blocks[c].second ] = pid_c(p);

src/blockforest/loadbalancing/StaticCurve.cpp

@@ -123,9 +123,9 @@ uint_t StaticLevelwiseCurveBalanceWeighted::operator()( SetupBlockForest & fores
       {
          const workload_t pWeight = totalWeight / workload_c( numberOfProcesses - p );
          workload_t weight( 0 );
-         while( c < blocksOnLevel.size() &&
+         while( c < blocksOnLevel.size() && ( isIdentical(weight, workload_t(0)) ||
                 std::abs( pWeight - weight - blocksOnLevel[c]->getWorkload() ) <=
-                std::abs( pWeight - weight ) )
+                std::abs( pWeight - weight ) ) )
          {
             blocksOnLevel[c]->assignTargetProcess(p);
 

src/core/Abort.cpp

@@ -63,7 +63,7 @@ void Abort::defaultAbort( const std::string & logMessage, const std::string & ca
 
 
 #ifdef _OPENMP
-   #pragma omp critical (abort)
+   #pragma omp critical (Abort_abort)
    {
 #endif
 

src/core/math/Random.cpp

@@ -43,7 +43,7 @@ std::mt19937 & getGenerator() // std::mt19937_64
 void seedRandomGenerator( const std::mt19937::result_type & seed )
 {
 #ifdef _OPENMP
-   #pragma omp critical (random)
+   #pragma omp critical (Random_random)
 #endif
    internal::getGenerator().seed( seed );
 }

src/core/math/Random.h

@@ -59,7 +59,7 @@ INT intRandom( const INT min, const INT max, std::mt19937 & generator )
 
    INT value;
 #ifdef _OPENMP
-   #pragma omp critical (random)
+   #pragma omp critical (Random_random)
 #endif
    { value = distribution( generator ); }
 
@@ -76,7 +76,7 @@ inline char intRandom<char>( const char min, const char max, std::mt19937 & gene
 
    char value;
 #ifdef _OPENMP
-   #pragma omp critical (random)
+   #pragma omp critical (Random_random)
 #endif
    { value = static_cast<char>( distribution( generator ) ); }
 
@@ -93,7 +93,7 @@ inline unsigned char intRandom<unsigned char>( const unsigned char min, const un
 
    unsigned char value;
 #ifdef _OPENMP
-   #pragma omp critical (random)
+   #pragma omp critical (Random_random)
 #endif
    { value = static_cast<unsigned char>( distribution( generator ) ); }
 
@@ -110,7 +110,7 @@ inline signed char intRandom<signed char>( const signed char min, const signed c
 
    signed char value;
 #ifdef _OPENMP
-   #pragma omp critical (random)
+   #pragma omp critical (Random_random)
 #endif
    { value = static_cast<signed char>( distribution( generator ) ); }
 
@@ -168,7 +168,7 @@ REAL realRandom( const REAL min = REAL(0), const REAL max = REAL(1), std::mt1993
 
    REAL value;
 #ifdef _OPENMP
-   #pragma omp critical (random)
+   #pragma omp critical (Random_random)
 #endif
    { value = distribution( generator ); }
 

src/field/interpolators/NearestNeighborFieldInterpolator.h

@@ -121,7 +121,7 @@ public:
                            *interpolationResultBegin = baseField_(curCell, f);
                            ++interpolationResultBegin;
                         }
-                        break;
+                        return;
                      }
                   }
                }

src/geometry/initializer/BoundaryFromDomainBorder.impl.h

@@ -106,7 +106,7 @@ void BoundaryFromDomainBorder<Handling>::init( stencil::Direction direction,
 
       const cell_idx_t wd = wallDistance;
 
-      CellInterval dBB = blocks_.getDomainCellBB();
+      CellInterval dBB = blocks_.getDomainCellBB( blocks_.getLevel(*blockIt) );
 
       for( uint_t dim = 0; dim< 3; ++dim )
          switch ( stencil::c[dim][direction] )

src/geometry/initializer/BoundarySetter.h

@@ -194,7 +194,10 @@ namespace initializer {
          }
 
          if ( boundaryConfigBlock_ )
+         {
             bcConfig_ = boundaryHandling_->createBoundaryConfiguration( boundaryUID_, boundaryConfigBlock_ );
+            boundaryConfigBlock_ = Config::BlockHandle(); // discard the config block so we don't unnecessarily run createBoundaryConfiguration more than once with the same arguments
+         }
 
          flag_ = boundaryHandling_->getBoundaryMask( boundaryUID_ );
          if ( ! field::isFlag( flag_ ) )

src/lbm/boundary/DynamicUBB.h

@@ -78,13 +78,29 @@ public:
       origin_[0] = aabb.xMin() + real_c(0.5) * dx_[0];
       origin_[1] = aabb.yMin() + real_c(0.5) * dx_[1];
       origin_[2] = aabb.zMin() + real_c(0.5) * dx_[2];
+
+      if( !timeTracker_ )
+      {
+         velocity_( time_ );
+      }
    }
+   DynamicUBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField * const pdfField,
+               const uint_t level, const VelocityFunctor_T & velocity, const AABB & aabb ) :
+      DynamicUBB( boundaryUID, uid, pdfField, nullptr, level, velocity, aabb )
+   {}
 
    shared_ptr< TimeTracker > getTimeTracker() { return timeTracker_; }
 
    void pushFlags( std::vector< FlagUID > & uids ) const { uids.push_back( uid_ ); }
 
-   void beforeBoundaryTreatment() { time_ = timeTracker_->getTime( level_ ); velocity_( time_ ); }
+   void beforeBoundaryTreatment()
+   {
+     if( timeTracker_ )
+     {
+        time_ = timeTracker_->getTime( level_ );
+        velocity_( time_ );
+     }
+   }
    void  afterBoundaryTreatment() const {}
 
    template< typename Buffer_T >
@@ -113,7 +129,7 @@ public:
       WALBERLA_ASSERT_EQUAL( nz, z + cell_idx_c( stencil::cz[ dir ] ) );
       WALBERLA_ASSERT_UNEQUAL( mask & this->mask_, numeric_cast<flag_t>(0) );
       WALBERLA_ASSERT_EQUAL( mask & this->mask_, this->mask_ ); // only true if "this->mask_" only contains one single flag, which is the case for the
-                                                                // current implementation of this boundary condition (SimpleUBB)
+                                                                // current implementation of this boundary condition (DynamicUBB)
 
       const Vector3< real_t > pos( origin_[0] + real_c(nx) * dx_[0],
                                    origin_[1] + real_c(ny) * dx_[1],

src/lbm/boundary/Pressure.h

@@ -211,7 +211,7 @@ inline void Pressure< LatticeModel_T, flag_t >::treatDirection( const cell_idx_t
 
    WALBERLA_ASSERT_UNEQUAL( ( mask & this->mask_ ), numeric_cast<flag_t>(0) );
    WALBERLA_ASSERT_EQUAL( ( mask & this->mask_ ), this->mask_ ); // only true if "this->mask_" only contains one single flag, which is the case for the
-                                                                 // current implementation of this boundary condition (SimplePressure)
+                                                                 // current implementation of this boundary condition (Pressure)
    Vector3<real_t> u = pdfField_->getVelocity(x,y,z);
 
    // result will be streamed to (x,y,z, stencil::inverseDir[d]) during sweep

src/lbm/boundary/all.h

@@ -31,6 +31,7 @@
 #include "NoDiffusion.h"
 #include "NoSlip.h"
 #include "Outlet.h"
+#include "ParserUBB.h"
 #include "Pressure.h"
 #include "SimpleDiffusionDirichlet.h"
 #include "SimplePAB.h"

src/lbm/boundary/factories/ExtendedBoundaryHandlingFactory.h

@@ -24,7 +24,7 @@
 #include "lbm/boundary/FreeSlip.h"
 #include "lbm/boundary/NoSlip.h"
 #include "lbm/boundary/Pressure.h"
-#include "lbm/boundary/UBB.h"
+#include "lbm/boundary/ParserUBB.h"
 #include "lbm/boundary/Outlet.h"
 #include "lbm/boundary/Curved.h"
 
@@ -82,7 +82,7 @@ public:
    typedef NoSlip< LatticeModel, flag_t >         BcNoSlip;
    typedef FreeSlip< LatticeModel, FlagFieldT >   BcFreeSlip;
    typedef Pressure< LatticeModel, flag_t >       BcPressure;
-   typedef UBB<LatticeModel, flag_t>              BcUBB;
+   typedef ParserUBB<LatticeModel, flag_t>        BcUBB;
    typedef Outlet<LatticeModel, FlagFieldT >      BcOutlet;
    typedef Curved<LatticeModel, FlagFieldT >      BcCurved;
 
@@ -144,7 +144,7 @@ ExtendedBoundaryHandlingFactory<LatticeModel, FlagFieldT>::ExtendedBoundaryHandl
 template <typename LatticeModel, typename FlagFieldT >
 typename ExtendedBoundaryHandlingFactory<LatticeModel, FlagFieldT>::BoundaryHandling *
 ExtendedBoundaryHandlingFactory<LatticeModel, FlagFieldT>::operator()( IBlock * const block,
-                                                                      const walberla::StructuredBlockStorage * const /*storage*/ ) const
+                                                                      const walberla::StructuredBlockStorage * const storage ) const
 {
    PdfFieldLM * const pdfField  = block->getData< PdfFieldLM >( pdfField_  );
    FlagFieldT * const flagField = block->getData< FlagFieldT >( flagField_ );
@@ -160,7 +160,7 @@ ExtendedBoundaryHandlingFactory<LatticeModel, FlagFieldT>::operator()( IBlock *
         BcNoSlip    ( getNoSlipBoundaryUID(),   getNoSlip(),   pdfField ),
         BcFreeSlip  ( getFreeSlipBoundaryUID(), getFreeSlip(), pdfField, flagField, mask ),
         BcPressure  ( getPressureBoundaryUID(), getPressure(), pdfField ),
-        BcUBB       ( getUBBBoundaryUID(),      getUBB(),      pdfField ),
+        BcUBB       ( getUBBBoundaryUID(),      getUBB(),      pdfField, flagField, storage->getLevel(*block), block->getAABB() ),
         BcOutlet    ( getOutletBoundaryUID(),   getOutlet(),   pdfField, flagField, mask ),
         BcCurved    ( getCurvedBoundaryUID(),   getCurved(),   pdfField, flagField, mask )
       )

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

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();

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() );
 

tests/fft/FftTest.cpp

@@ -41,7 +41,7 @@ int main (int argc, char** argv)
    {
       Field_T *data_in = block->getData< Field_T >( originalFieldId );
       Field_T *data_out = block->getData< Field_T >( fftFieldId );
-      WALBERLA_FOR_ALL_CELLS_XYZ(data_in, {
+      WALBERLA_FOR_ALL_CELLS_XYZ_OMP(data_in, omp critical, {
          Vector3<real_t> point( real_c(x), real_c(y), real_c(z) );
          blocks->transformBlockLocalToGlobal(point, *block);
          data_in->get(x,y,z) = real_c(std::ranlux48_base(uint_c(point[0])+(uint_c(point[1])*L+uint_c(point[2]))*L)())*real_c(std::pow(2,-48));

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;
 }
 

tests/pe/MinMaxRefinement.cpp

@@ -144,7 +144,7 @@ int main( int argc, char ** argv )
       ccd->reloadBodies();
    }
 
-   WALBERLA_CHECK_EQUAL( blockforest.size(), mpi::MPIManager::instance()->worldRank() == 6 ? 1 : 0);
+   WALBERLA_CHECK_EQUAL( blockforest.size(), mpi::MPIManager::instance()->worldRank() == 0 ? 1 : 0);
    WALBERLA_LOG_DEVEL( infoCollection->size() );
 
    for (unsigned int i = 0; i < 30; ++i)