diff --git a/src/mesh/boundary/BoundarySetup.h b/src/mesh/boundary/BoundarySetup.h
index 1512fe618a89281707dae40e6c6465b01fb0258d..bdc0fa5098aa56b33407784e8f35b8f7ec95c125 100644
--- a/src/mesh/boundary/BoundarySetup.h
+++ b/src/mesh/boundary/BoundarySetup.h
@@ -66,6 +66,11 @@ public:
template<typename FlagField_T>
void setFlag( const BlockDataID flagFieldID, field::FlagUID flagUID, Location boundaryLocation );
+ // set flag with FlagUID wherever boundaryUID is set
+ template< typename FlagField_T, typename BoundaryFunction, typename Stencil = stencil::D3Q27 >
+ void setBoundaryFlag(const BlockDataID flagFieldID, field::FlagUID flagUID, boundary::BoundaryUID boundaryUID,
+ const BoundaryFunction& boundaryFunction, Location boundaryLocation);
+
void writeVTKVoxelfile( const std::string & identifier = "voxelization", bool writeGhostLayers = false, const std::string & baseFolder = std::string("vtk_out"),
const std::string & executionFolder = std::string("voxelization") );
@@ -83,7 +88,7 @@ private:
shared_ptr< BlockDataID > voxelizationFieldId_;
DistanceFunction distanceFunction_; /// function providing the squared signed distance to an object
uint_t numGhostLayers_;
- size_t cellVectorChunkSize_; /// Number of boundary cells which are setup simultaneously
+ size_t cellVectorChunkSize_; /// Number of boundary cells which are setup simultaneously
};
@@ -105,19 +110,15 @@ void BoundarySetup::setDomainCells( const BlockDataID boundaryHandlingId, const
std::vector<Cell> domainCells;
- for( cell_idx_t z = -cell_idx_c( numGhostLayers_ ); z < cell_idx_c( voxelizationField->zSize() + numGhostLayers_ ); ++z )
- for( cell_idx_t y = -cell_idx_c( numGhostLayers_ ); y < cell_idx_c( voxelizationField->ySize() + numGhostLayers_ ); ++y )
- 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.emplace_back( x, y, z );
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(voxelizationField, {
+ if (voxelizationField->get(x, y, z) == domainValue) { domainCells.emplace_back(x, y, z); }
- if( domainCells.size() > cellVectorChunkSize_ )
- {
- boundaryHandling->setDomain( domainCells.begin(), domainCells.end() );
- domainCells.clear();
- }
- }
+ if (domainCells.size() > cellVectorChunkSize_)
+ {
+ boundaryHandling->setDomain(domainCells.begin(), domainCells.end());
+ domainCells.clear();
+ }
+ })
boundaryHandling->setDomain( domainCells.begin(), domainCells.end() );
}
@@ -140,17 +141,69 @@ void BoundarySetup::setFlag( const BlockDataID flagFieldID, field::FlagUID flagU
const uint8_t domainValue = boundaryLocation == INSIDE ? uint8_t(0) : uint8_t(1);
- for( cell_idx_t z = -cell_idx_c( numGhostLayers_ ); z < cell_idx_c( voxelizationField->zSize() + numGhostLayers_ ); ++z )
- for( cell_idx_t y = -cell_idx_c( numGhostLayers_ ); y < cell_idx_c( voxelizationField->ySize() + numGhostLayers_ ); ++y )
- for( cell_idx_t x = -cell_idx_c( numGhostLayers_ ); x < cell_idx_c( voxelizationField->xSize() + numGhostLayers_ ); ++x )
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(voxelizationField, {
+ if (voxelizationField->get(x, y, z) != domainValue) { flagField->addFlag(x, y, z, flag); }
+ })
+ }
+}
+
+template< typename FlagField_T, typename BoundaryFunction, typename Stencil >
+void BoundarySetup::setBoundaryFlag(const BlockDataID flagFieldID, field::FlagUID flagUID,
+ boundary::BoundaryUID boundaryUID, const BoundaryFunction& boundaryFunction,
+ Location boundaryLocation)
+{
+ for (auto& block : *structuredBlockStorage_)
+ {
+ FlagField_T* flagField = block.getData< FlagField_T >(flagFieldID);
+ WALBERLA_CHECK_NOT_NULLPTR(flagField, "flagFieldID invalid!");
+ auto flag = flagField->getFlag(flagUID);
+
+ 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!");
+
+ // get where the (fluid) domain is located (on the inside/outside of the mesh)
+ const uint8_t domainValue = boundaryLocation == INSIDE ? uint8_t(0) : uint8_t(1);
+
+ // get the cell interval covering the whole block
+ const CellInterval blockCi = voxelizationField->xyzSizeWithGhostLayer();
+
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(voxelizationField, {
+ // only boundaries are set here, i.e., skip regions inside the (fluid) domain
+ if (voxelizationField->get(x, y, z) == domainValue) { continue; }
+
+ for (auto dirIt = Stencil::beginNoCenter(); dirIt != Stencil::end(); ++dirIt)
+ {
+ const Cell cell(x, y, z);
+ const Cell neighborCell = cell + *dirIt;
+
+ // this cell (x,y,z) is considered a boundary cell if the neighboring cell (x,y,z)+*dirIt
+ // - is still located on the same block AND
+ // - belongs to the (fluid) domain, i.e., not to the boundary
+ if (blockCi.contains(neighborCell) && voxelizationField->get(neighborCell) == domainValue)
{
- if( voxelizationField->get( x, y, z ) != domainValue )
+ const Vector3< real_t > cellCenter = structuredBlockStorage_->getBlockLocalCellCenter(block, cell);
+ const BoundaryInfo& bi = boundaryFunction(cellCenter);
+
+ if (boundaryUID == bi.getUid())
+ {
+ // clear existing flags
+ flagField->removeMask(x, y, z, flagField->get(x, y, z));
+
+ // add new flag
flagField->addFlag(x, y, z, flag);
+ }
+
+ break;
}
+ }
+ })
}
}
-
template< typename BoundaryHandlingType, typename BoundaryFunction, typename Stencil >
void BoundarySetup::setBoundaries( const BlockDataID boundaryHandlingID, const BoundaryFunction & boundaryFunction, Location boundaryLocation )
{
@@ -168,33 +221,27 @@ void BoundarySetup::setBoundaries( const BlockDataID boundaryHandlingID, const B
const CellInterval blockCi = voxelizationField->xyzSizeWithGhostLayer();
- for( cell_idx_t z = -cell_idx_c( numGhostLayers_ ); z < cell_idx_c( voxelizationField->zSize() + numGhostLayers_ ); ++z )
- for( cell_idx_t y = -cell_idx_c( numGhostLayers_ ); y < cell_idx_c( voxelizationField->ySize() + numGhostLayers_ ); ++y )
- for( cell_idx_t x = -cell_idx_c( numGhostLayers_ ); x < cell_idx_c( voxelizationField->xSize() + numGhostLayers_ ); ++x )
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(voxelizationField, {
+ if (voxelizationField->get(x, y, z) == domainValue) { continue; }
+
+ for (auto dirIt = Stencil::beginNoCenter(); dirIt != Stencil::end(); ++dirIt)
+ {
+ const Cell cell(x, y, z);
+ const Cell neighborCell = cell + *dirIt;
+ if (blockCi.contains(neighborCell) && voxelizationField->get(neighborCell) == domainValue)
{
- if( voxelizationField->get( x, y, z ) == domainValue )
- continue;
+ const Vector3< real_t > cellCenter = structuredBlockStorage_->getBlockLocalCellCenter(block, cell);
+ const BoundaryInfo& bi = boundaryFunction(cellCenter);
+ const auto boundaryMask = boundaryHandling->getBoundaryMask(bi.getUid());
- for(auto dirIt = Stencil::beginNoCenter(); dirIt != Stencil::end(); ++dirIt)
- {
- const Cell cell( x, y, z );
- const Cell neighborCell = cell + *dirIt;
- if( blockCi.contains( neighborCell ) && voxelizationField->get( neighborCell ) == domainValue )
- {
- const Vector3< real_t > cellCenter = structuredBlockStorage_->getBlockLocalCellCenter( block, cell );
- const BoundaryInfo & bi = boundaryFunction( cellCenter );
- const auto boundaryMask = boundaryHandling->getBoundaryMask( bi.getUid() );
-
- boundaryHandling->setBoundary( boundaryMask, cell.x(), cell.y(), cell.z(), *( bi.getConfig() ) );
-
- break;
- }
- }
+ boundaryHandling->setBoundary(boundaryMask, cell.x(), cell.y(), cell.z(), *(bi.getConfig()));
+
+ break;
}
+ }
+ })
}
}
-
-
} // namespace mesh
} // namespace walberla
\ No newline at end of file