diff --git a/CMakeLists.txt b/CMakeLists.txt
index 85125dcaf83ccf81a7ef108eb38f92b4b3ee89ed..0ee1e578c7a963ec1f5f40c907fce8aa647b1b6b 100644
--- a/CMakeLists.txt
+++ b/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()
############################################################################################################################
diff --git a/src/blockforest/loadbalancing/DynamicCurve.h b/src/blockforest/loadbalancing/DynamicCurve.h
index d9767ecf055e8b1106bbcd2746381c5be4b89db1..d72b60726f0b1dab41e8e5708051b1bc7b238a22 100644
--- a/src/blockforest/loadbalancing/DynamicCurve.h
+++ b/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);
diff --git a/src/blockforest/loadbalancing/StaticCurve.cpp b/src/blockforest/loadbalancing/StaticCurve.cpp
index 9bb9c184dcc24b10cfd1bc93f388d59523911301..e90ec3bc10ac72a04ebbfe9ce6609552c4276158 100644
--- a/src/blockforest/loadbalancing/StaticCurve.cpp
+++ b/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);
diff --git a/src/core/Abort.cpp b/src/core/Abort.cpp
index 910dcb209f178dcdd1595477f8339eb08726d06a..1ff3ad2cb045f2132518f6a854fd79f48053a1a8 100644
--- a/src/core/Abort.cpp
+++ b/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
diff --git a/src/core/math/Random.cpp b/src/core/math/Random.cpp
index 72707f208e5053cba44daed005e2685f86333e8b..5acdc01d34fbe236fcb94fb98fa952a92564eaa2 100644
--- a/src/core/math/Random.cpp
+++ b/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 );
}
diff --git a/src/core/math/Random.h b/src/core/math/Random.h
index ec61aa51d65b76f3a10c243086d08a48eb166590..d63d2e7ede6338606219e743d7fc44444ca23010 100644
--- a/src/core/math/Random.h
+++ b/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 ); }
diff --git a/src/field/interpolators/NearestNeighborFieldInterpolator.h b/src/field/interpolators/NearestNeighborFieldInterpolator.h
index 45e9110de725173af89a4aed20e1d6290b5278a5..975c55edb4c52d0ea197391a88b271588bbf7729 100644
--- a/src/field/interpolators/NearestNeighborFieldInterpolator.h
+++ b/src/field/interpolators/NearestNeighborFieldInterpolator.h
@@ -121,7 +121,7 @@ public:
*interpolationResultBegin = baseField_(curCell, f);
++interpolationResultBegin;
}
- break;
+ return;
}
}
}
diff --git a/src/geometry/initializer/BoundaryFromDomainBorder.impl.h b/src/geometry/initializer/BoundaryFromDomainBorder.impl.h
index 9984fcaca34a6a4da0ca99d5e385d7b7469275da..dca93f84d01ce55b345fe585b1b34263856249f8 100644
--- a/src/geometry/initializer/BoundaryFromDomainBorder.impl.h
+++ b/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] )
diff --git a/src/geometry/initializer/BoundarySetter.h b/src/geometry/initializer/BoundarySetter.h
index 256ec58444530bcce20502a4541cab8b3a1b951c..b9fe341228dca130ada4f2f8647d98cfcc1f9b3d 100644
--- a/src/geometry/initializer/BoundarySetter.h
+++ b/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_ ) )
diff --git a/src/lbm/boundary/DynamicUBB.h b/src/lbm/boundary/DynamicUBB.h
index 03f4af7c8fc980b210eff4576bf64f0a7d893010..67425aa05cc721a13cfe746e000520bbf6052fdd 100644
--- a/src/lbm/boundary/DynamicUBB.h
+++ b/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],
diff --git a/src/lbm/boundary/ParserUBB.h b/src/lbm/boundary/ParserUBB.h
new file mode 100644
index 0000000000000000000000000000000000000000..82bfff5e2ba37d7366f1bf1318a1992e664a7602
--- /dev/null
+++ b/src/lbm/boundary/ParserUBB.h
@@ -0,0 +1,546 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file ParserUBB.h
+//! \ingroup lbm
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "lbm/field/DensityAndVelocity.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/refinement/TimeTracker.h"
+
+#include "boundary/Boundary.h"
+
+#include "core/DataTypes.h"
+#include "core/cell/CellInterval.h"
+#include "core/config/Config.h"
+#include "core/debug/Debug.h"
+#include "core/math/Vector3.h"
+#include "core/math/ParserOMP.h"
+
+#include "field/FlagUID.h"
+
+#include "stencil/Directions.h"
+
+#include <vector>
+
+
+
+namespace walberla {
+namespace lbm {
+
+
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce = false >
+class ParserUBB : public Boundary<flag_t>
+{
+ typedef lbm::PdfField< LatticeModel_T > PDFField;
+ typedef typename LatticeModel_T::Stencil Stencil;
+
+public:
+
+ static const bool threadsafe = true;
+
+ class Parser : public BoundaryConfiguration
+ {
+ public:
+ inline Parser( const Config::BlockHandle & config );
+ inline Parser( std::array< std::string, 3 > & equations );
+ Vector3< real_t > operator()( const Vector3< real_t > & x, const real_t t ) const;
+ Vector3< real_t > operator()( const Vector3< real_t > & x ) const;
+ bool isTimeDependent() const { return timeDependent_; }
+ const std::array< std::string, 3 > & equations() const { return equations_; }
+
+ private:
+ std::array< math::FunctionParserOMP, 3 > parsers_;
+ std::array< std::string, 3 > equations_;
+ bool timeDependent_;
+ }; // class Parser
+
+ // constant velocity class is for API compatibility with UBB
+ class Velocity : public BoundaryConfiguration {
+ public:
+ Velocity( const Vector3< real_t > & _velocity ) : velocity_( _velocity ) {}
+ Velocity( const real_t _x, const real_t _y, const real_t _z ) : velocity_(_x,_y,_z) {}
+ inline Velocity( const Config::BlockHandle & config );
+
+ const Vector3< real_t > & velocity() const { return velocity_; }
+
+ const real_t & x() const { return velocity_[0]; }
+ const real_t & y() const { return velocity_[1]; }
+ const real_t & z() const { return velocity_[2]; }
+
+ Vector3< real_t > & velocity() { return velocity_; }
+
+ real_t & x() { return velocity_[0]; }
+ real_t & y() { return velocity_[1]; }
+ real_t & z() { return velocity_[2]; }
+
+ private:
+ Vector3< real_t > velocity_;
+ }; // class Velocity
+
+ typedef GhostLayerField< shared_ptr<Parser>, 1 > ParserField;
+ typedef GhostLayerField< Vector3<real_t>, 1 > VelocityField;
+
+ static shared_ptr<Parser> createConfiguration( const Config::BlockHandle & config )
+ { return make_shared<Parser>( config ); }
+
+
+
+ ParserUBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField * const pdfField,
+ FlagField<flag_t> * const flagField, const shared_ptr< TimeTracker > & timeTracker,
+ const uint_t level, const AABB & aabb );
+ ParserUBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField * const pdfField,
+ FlagField<flag_t> * const flagField, const uint_t level, const AABB & aabb );
+
+ shared_ptr< TimeTracker > getTimeTracker() { return timeTracker_; }
+
+ void pushFlags( std::vector< FlagUID > & uids ) const { uids.push_back( uid_ ); }
+
+ void beforeBoundaryTreatment()
+ {
+ if( timeTracker_ )
+ time_ = timeTracker_->getTime( level_ );
+ }
+ void afterBoundaryTreatment() const {}
+
+ template< typename Buffer_T >
+ inline void packCell( Buffer_T &, const cell_idx_t, const cell_idx_t, const cell_idx_t ) const;
+
+ template< typename Buffer_T >
+ inline void registerCell( Buffer_T &, const flag_t, const cell_idx_t, const cell_idx_t, const cell_idx_t );
+
+ inline void registerCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const BoundaryConfiguration & parser );
+ inline void registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & parser );
+ template< typename CellIterator >
+ inline void registerCells( const flag_t, const CellIterator & begin, const CellIterator & end, const BoundaryConfiguration & parser );
+
+ inline void unregisterCell( const flag_t, const cell_idx_t, const cell_idx_t, const cell_idx_t ) const;
+
+#ifndef NDEBUG
+ inline void treatDirection( 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, const flag_t mask );
+#else
+ inline void treatDirection( 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, const flag_t /*mask*/ );
+#endif
+
+ inline Vector3<real_t> getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const;
+ inline Vector3<real_t> getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z, real_t t ) const;
+
+private:
+
+ const FlagUID uid_;
+
+ PDFField * const pdfField_;
+
+ Vector3< real_t > origin_;
+ Vector3< real_t > dx_;
+
+ // required to keep track of the simulation time
+ shared_ptr< TimeTracker > timeTracker_; // -> addPostBoundaryHandlingVoidFunction (when used with refinement time step)
+ real_t time_;
+ uint_t level_;
+
+ shared_ptr<ParserField> parserField_;
+ shared_ptr<VelocityField> velocityField_;
+
+}; // class ParserUBB
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
+inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::Parser( const Config::BlockHandle & config )
+: parsers_(), equations_(), timeDependent_( false )
+{
+ if( !config )
+ return;
+
+ if( config.isDefined( "x" ) )
+ {
+ equations_[0] = config.getParameter<std::string>( "x" );
+ parsers_[0].parse( equations_[0] );
+ if( parsers_[0].symbolExists( "t" ) )
+ timeDependent_ = true;
+ }
+ if( config.isDefined( "y" ) )
+ {
+ equations_[1] = config.getParameter<std::string>( "y" );
+ parsers_[1].parse( equations_[1] );
+ if( parsers_[1].symbolExists( "t" ) )
+ timeDependent_ = true;
+ }
+ if( config.isDefined( "z" ) )
+ {
+ equations_[2] = config.getParameter<std::string>( "z" );
+ parsers_[2].parse( equations_[2] );
+ if( parsers_[2].symbolExists( "t" ) )
+ timeDependent_ = true;
+ }
+}
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
+inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::Parser( std::array< std::string, 3 > & equations )
+: parsers_(), equations_( equations ), timeDependent_( false )
+{
+ if( equations_[0].length() > 0 )
+ {
+ parsers_[0].parse( equations_[0] );
+ if( parsers_[0].symbolExists( "t" ) )
+ timeDependent_ = true;
+ }
+ if( equations_[1].length() > 0 )
+ {
+ parsers_[1].parse( equations_[1] );
+ if( parsers_[1].symbolExists( "t" ) )
+ timeDependent_ = true;
+ }
+ if( equations_[2].length() > 0 )
+ {
+ parsers_[2].parse( equations_[2] );
+ if( parsers_[2].symbolExists( "t" ) )
+ timeDependent_ = true;
+ }
+}
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+inline ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::Velocity::Velocity( const Config::BlockHandle & config )
+{
+ velocity_[0] = ( config && config.isDefined( "x" ) ) ? config.getParameter<real_t>( "x" ) : real_c(0.0);
+ velocity_[1] = ( config && config.isDefined( "y" ) ) ? config.getParameter<real_t>( "y" ) : real_c(0.0);
+ velocity_[2] = ( config && config.isDefined( "z" ) ) ? config.getParameter<real_t>( "z" ) : real_c(0.0);
+}
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
+Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::operator()( const Vector3< real_t > & x, const real_t t ) const
+{
+ std::map< std::string, double > symbols;
+ symbols["x"] = x[0];
+ symbols["y"] = x[1];
+ symbols["z"] = x[2];
+ symbols["t"] = t;
+
+ Vector3< real_t > v;
+ v[0] = parsers_[0].evaluate( symbols );
+ v[1] = parsers_[1].evaluate( symbols );
+ v[2] = parsers_[2].evaluate( symbols );
+ return v;
+}
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
+Vector3< real_t > ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::Parser::operator()( const Vector3< real_t > & x ) const
+{
+ WALBERLA_ASSERT( !timeDependent_ );
+
+ std::map< std::string, double > symbols;
+ symbols["x"] = x[0];
+ symbols["y"] = x[1];
+ symbols["z"] = x[2];
+
+ Vector3< real_t > v;
+ v[0] = parsers_[0].evaluate( symbols );
+ v[1] = parsers_[1].evaluate( symbols );
+ v[2] = parsers_[2].evaluate( symbols );
+ return v;
+}
+
+
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
+inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::ParserUBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField * const pdfField,
+ FlagField<flag_t> * const flagField, const shared_ptr< TimeTracker > & timeTracker,
+ const uint_t level, const AABB & aabb )
+ : Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField ), timeTracker_( timeTracker ), time_( real_t(0) ), level_( level )
+{
+ WALBERLA_ASSERT_NOT_NULLPTR( pdfField_ );
+ dx_[0] = aabb.xSize() / real_c( pdfField_->xSize() );
+ dx_[1] = aabb.ySize() / real_c( pdfField_->ySize() );
+ dx_[2] = aabb.zSize() / real_c( pdfField_->zSize() );
+ 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(flagField != NULL)
+ {
+ parserField_ = make_shared<ParserField> ( pdfField->xSize(), pdfField->ySize(), pdfField->zSize(), flagField->nrOfGhostLayers(), field::zyxf );
+ velocityField_ = make_shared<VelocityField>( pdfField->xSize(), pdfField->ySize(), pdfField->zSize(), flagField->nrOfGhostLayers(), field::zyxf );
+ }
+ else
+ {
+ parserField_ = make_shared<ParserField> ( pdfField->xSize(), pdfField->ySize(), pdfField->zSize(), pdfField->nrOfGhostLayers(), field::zyxf );
+ velocityField_ = make_shared<VelocityField>( pdfField->xSize(), pdfField->ySize(), pdfField->zSize(), pdfField->nrOfGhostLayers(), field::zyxf );
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+template< typename Buffer_T >
+inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::packCell( Buffer_T & buffer, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
+{
+ if( parserField_->get( x, y, z ) )
+ {
+ auto & eqs = parserField_->get( x, y, z )->equations();
+ buffer << true << eqs[0] << eqs[1] << eqs[2];
+ }
+ else
+ {
+ buffer << false << velocityField_->get( x, y, z );
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+template< typename Buffer_T >
+inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCell( Buffer_T & buffer, const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
+{
+ bool isparser;
+ buffer >> isparser;
+ if( isparser )
+ {
+ std::array< std::string, 3> eqs;
+ buffer >> eqs[0] >> eqs[1] >> eqs[2];
+
+ auto p = make_shared<Parser>(eqs);
+ parserField_->get( x, y, z ) = p;
+ }
+ else
+ {
+ buffer >> velocityField_->get( x, y, z );
+ parserField_->get( x, y, z ) = nullptr;
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
+inline ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::ParserUBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField * const pdfField,
+ FlagField<flag_t> * const flagField, const uint_t level, const AABB & aabb )
+ : ParserUBB( boundaryUID, uid, pdfField, flagField, nullptr, level, aabb )
+{}
+
+
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
+ const BoundaryConfiguration & parser )
+{
+ WALBERLA_ASSERT_EQUAL( dynamic_cast< const Parser * >( &parser ), &parser );
+ WALBERLA_ASSERT_NOT_NULLPTR( parserField_ );
+
+ if( auto v = dynamic_cast< const Velocity * >( &parser ) )
+ {
+ velocityField_->get( x, y, z ) = v->velocity();
+ parserField_->get( x, y, z ) = nullptr;
+ return;
+ }
+
+ auto & p = dynamic_cast< const Parser & >( parser );
+
+ if( p.isTimeDependent() )
+ {
+ parserField_->get( x, y, z ) = make_shared<Parser>(p);
+ }
+ else
+ {
+ const Vector3< real_t > pos( origin_[0] + real_c(x) * dx_[0],
+ origin_[1] + real_c(y) * dx_[1],
+ origin_[2] + real_c(z) * dx_[2] );
+ velocityField_->get( x, y, z ) = p(pos);
+ parserField_->get( x, y, z ) = nullptr;
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCells( const flag_t, const CellInterval & cells, const BoundaryConfiguration & parser )
+{
+ WALBERLA_ASSERT_EQUAL( dynamic_cast< const Parser * >( &parser ), &parser );
+ WALBERLA_ASSERT_NOT_NULLPTR( parserField_ );
+
+ if( auto v = dynamic_cast< const Velocity * >( &parser ) )
+ {
+ for( auto cell = parserField_->beginSliceXYZ( cells ); cell != parserField_->end(); ++cell )
+ {
+ velocityField_->get( cell.x(), cell.y(), cell.z() ) = v->velocity();
+ *cell = nullptr;
+ }
+ return;
+ }
+
+ auto & p = dynamic_cast< const Parser & >( parser );
+
+ if( p.isTimeDependent() )
+ {
+ auto shared_p = make_shared<Parser>(p);
+ for( auto cell = parserField_->beginSliceXYZ( cells ); cell != parserField_->end(); ++cell )
+ *cell = shared_p;
+ }
+ else
+ {
+ for( auto cell = parserField_->beginSliceXYZ( cells ); cell != parserField_->end(); ++cell )
+ {
+ const Vector3< real_t > pos( origin_[0] + real_c(cell.x()) * dx_[0],
+ origin_[1] + real_c(cell.y()) * dx_[1],
+ origin_[2] + real_c(cell.z()) * dx_[2] );
+ velocityField_->get( cell.x(), cell.y(), cell.z() ) = p(pos);
+ *cell = nullptr;
+ }
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+template< typename CellIterator >
+inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::registerCells( const flag_t, const CellIterator & begin, const CellIterator & end,
+ const BoundaryConfiguration & parser )
+{
+ WALBERLA_ASSERT_EQUAL( dynamic_cast< const Parser * >( &parser ), &parser );
+ WALBERLA_ASSERT_NOT_NULLPTR( parserField_ );
+
+ if( auto v = dynamic_cast< const Velocity * >( &parser ) )
+ {
+ for( auto cell = begin; cell != end; ++cell )
+ {
+ velocityField_->get( cell.x(), cell.y(), cell.z() ) = v->velocity();
+ parserField_->get( cell->x(), cell->y(), cell->z() ) = nullptr;
+ }
+ return;
+ }
+
+ auto & p = dynamic_cast< const Parser & >( parser );
+
+ if( p.isTimeDependent() )
+ {
+ auto shared_p = make_shared<Parser>(p);
+ for( auto cell = begin; cell != end; ++cell )
+ {
+ parserField_->get( cell->x(), cell->y(), cell->z() ) = shared_p;
+ }
+ }
+ else
+ {
+ for( auto cell = begin; cell != end; ++cell )
+ {
+ const Vector3< real_t > pos( origin_[0] + real_c(cell->x()) * dx_[0],
+ origin_[1] + real_c(cell->y()) * dx_[1],
+ origin_[2] + real_c(cell->z()) * dx_[2] );
+ velocityField_->get( cell->x(), cell->y(), cell->z() ) = p(pos);
+ parserField_->get( cell->x(), cell->y(), cell->z() ) = nullptr;
+ }
+ }
+}
+
+
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+inline void ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::unregisterCell( const flag_t, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
+{
+ parserField_->get(x,y,z) = nullptr;
+}
+
+
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce>
+#ifndef NDEBUG
+inline void ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::ParserUBB::treatDirection( 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, const flag_t mask )
+#else
+inline void ParserUBB<LatticeModel_T, flag_t, AdaptVelocityToExternalForce>::ParserUBB::treatDirection( 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, const flag_t /*mask*/ )
+#endif
+ {
+ WALBERLA_ASSERT_EQUAL( nx, x + cell_idx_c( stencil::cx[ dir ] ) );
+ WALBERLA_ASSERT_EQUAL( ny, y + cell_idx_c( stencil::cy[ dir ] ) );
+ 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 (ParserUBB)
+
+ Vector3<real_t> velocity;
+ if( parserField_->get(nx,ny,nz) )
+ {
+ WALBERLA_ASSERT_NOT_NULLPTR( getTimeTracker(), "A TimeTracker is needed for time-dependent equations" );
+ const Vector3< real_t > pos( origin_[0] + real_c(nx) * dx_[0],
+ origin_[1] + real_c(ny) * dx_[1],
+ origin_[2] + real_c(nz) * dx_[2] );
+ velocity = (*parserField_->get(nx,ny,nz))( pos, time_ );
+ }
+ else
+ {
+ velocity = velocityField_->get(nx,ny,nz);
+ }
+
+ if( LatticeModel_T::compressible )
+ {
+ const auto density = pdfField_->getDensity(x,y,z);
+ const auto vel = AdaptVelocityToExternalForce ? lbm::internal::AdaptVelocityToForce<LatticeModel_T>::get( x, y, z, pdfField_->latticeModel(), velocity, density ) :
+ velocity;
+
+ pdfField_->get( nx, ny, nz, Stencil::invDirIdx(dir) ) = pdfField_->get( x, y, z, Stencil::idx[dir] ) -
+ ( real_c(6) * density * real_c(LatticeModel_T::w[ Stencil::idx[dir] ]) *
+ ( real_c(stencil::cx[ dir ]) * vel[0] +
+ real_c(stencil::cy[ dir ]) * vel[1] +
+ real_c(stencil::cz[ dir ]) * vel[2] ) );
+ }
+ else
+ {
+ const auto vel = AdaptVelocityToExternalForce ? lbm::internal::AdaptVelocityToForce<LatticeModel_T>::get( x, y, z, pdfField_->latticeModel(), velocity, real_t(1) ) :
+ velocity;
+
+ pdfField_->get( nx, ny, nz, Stencil::invDirIdx(dir) ) = pdfField_->get( x, y, z, Stencil::idx[dir] ) -
+ ( real_c(6) * real_c(LatticeModel_T::w[ Stencil::idx[dir] ]) *
+ ( real_c(stencil::cx[ dir ]) * vel[0] +
+ real_c(stencil::cy[ dir ]) * vel[1] +
+ real_c(stencil::cz[ dir ]) * vel[2] ) );
+ }
+ }
+
+
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+inline Vector3<real_t> ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z ) const
+{
+ return getValue( x, y, z, time_ );
+}
+
+template< typename LatticeModel_T, typename flag_t, bool AdaptVelocityToExternalForce >
+inline Vector3<real_t> ParserUBB< LatticeModel_T, flag_t, AdaptVelocityToExternalForce >::getValue( const cell_idx_t x, cell_idx_t y, cell_idx_t z, real_t t ) const
+{
+ Vector3<real_t> velocity;
+ if( parserField_->get(x,y,z) )
+ {
+ WALBERLA_ASSERT_NOT_NULLPTR( getTimeTracker(), "A TimeTracker is needed for time-dependent equations" );
+ const Vector3< real_t > pos( origin_[0] + real_c(x) * dx_[0],
+ origin_[1] + real_c(y) * dx_[1],
+ origin_[2] + real_c(z) * dx_[2] );
+ velocity = (*parserField_->get(x,y,z))( pos, t );
+ }
+ else
+ {
+ velocity = velocityField_->get(x,y,z);
+ }
+ return velocity;
+}
+
+
+} // namespace lbm
+} // namespace walberla
diff --git a/src/lbm/boundary/Pressure.h b/src/lbm/boundary/Pressure.h
index a637b14d17133c33369a546163961a75e8f465a6..ae8afb05372001a7062c72df7e53a645e566e1b5 100644
--- a/src/lbm/boundary/Pressure.h
+++ b/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
diff --git a/src/lbm/boundary/all.h b/src/lbm/boundary/all.h
index 158a9f887948684e90d55437748ff481b52e22a7..31709c60e0e538dddaaf537fe0d661fb2b4833c1 100644
--- a/src/lbm/boundary/all.h
+++ b/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"
diff --git a/src/lbm/boundary/factories/ExtendedBoundaryHandlingFactory.h b/src/lbm/boundary/factories/ExtendedBoundaryHandlingFactory.h
index b12fdc4bb5ed9046d01ca2733e7e9ba66503d1c6..11ae8af9d3ca08da723e218d37b90796a6e35e85 100644
--- a/src/lbm/boundary/factories/ExtendedBoundaryHandlingFactory.h
+++ b/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 )
)
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/fft/FftTest.cpp b/tests/fft/FftTest.cpp
index 58d80c486de68e7c33d9e2a14263fd08cdd0eeb7..99e934dbfca25d21ad8fb7883634ac6d591c6bd4 100644
--- a/tests/fft/FftTest.cpp
+++ b/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));
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;
}
diff --git a/tests/pe/MinMaxRefinement.cpp b/tests/pe/MinMaxRefinement.cpp
index 458052c0e4f1853e3ab03b8f45e91c040915699e..3c13d4cfbc7fc6aa974a9cb6f0688e98d813965c 100644
--- a/tests/pe/MinMaxRefinement.cpp
+++ b/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)