Commit 25f35ef4 authored by Christoph Rettinger's avatar Christoph Rettinger

Merge branch 'pe_coupling_refinement_fixes' into 'master'

Pe coupling refinement fixes

See merge request walberla/walberla!67
parents 26b7f120 87eb7592
Pipeline #6337 passed with stages
in 80 minutes and 58 seconds
add_subdirectory( CouetteFlow )
add_subdirectory( ForcesOnSphereNearPlaneInShearFlow )
add_subdirectory( NonUniformGrid )
add_subdirectory( MotionSingleHeavySphere )
add_subdirectory( PoiseuilleChannel )
......
waLBerla_add_executable ( NAME ForcesOnSphereNearPlaneInShearFlow
DEPENDS blockforest boundary core domain_decomposition field lbm pe pe_coupling postprocessing timeloop vtk )
\ No newline at end of file
......@@ -900,9 +900,9 @@ int main( int argc, char **argv )
//initialization of the PDFs inside the particles, important for PSM M3
if( psmVariant == PSMVariant::SC1W1 || psmVariant == PSMVariant::SC2W1 || psmVariant == PSMVariant::SC3W1 )
pe_coupling::initializeDomainForPSM< LatticeModel_T, 1 >( blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
pe_coupling::initializeDomainForPSM< LatticeModel_T, 1 >( *blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
else
pe_coupling::initializeDomainForPSM< LatticeModel_T, 2 >( blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
pe_coupling::initializeDomainForPSM< LatticeModel_T, 2 >( *blocks, pdfFieldID, bodyAndVolumeFractionFieldID );
}
else
......@@ -912,11 +912,11 @@ int main( int argc, char **argv )
if( memVariant == MEMVariant::CLI )
{
pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_CLI_Flag );
pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_CLI_Flag );
}else if ( memVariant == MEMVariant::MR ){
pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_MR_Flag );
pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_MR_Flag );
}else{
pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_BB_Flag );
pe_coupling::mapMovingBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MEM_BB_Flag );
}
}
......
......@@ -30,8 +30,8 @@ real_t intersectionRatioBisection( const Body & body,
const Vector3<real_t> & direction,
const real_t epsilon )
{
WALBERLA_ASSERT( !geometry::contains( body, fluidPoint ) );
WALBERLA_ASSERT( geometry::contains( body, fluidPoint + direction ) );
WALBERLA_ASSERT( !geometry::contains( body, fluidPoint ), "fluid point: " << fluidPoint );
WALBERLA_ASSERT( geometry::contains( body, fluidPoint + direction ), "fluidPoint + direction: " << fluidPoint + direction );
const real_t sqEpsilon = epsilon * epsilon;
const real_t sqDirectionLength = direction.sqrLength();
......
......@@ -658,7 +658,7 @@ void PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::communicateLocalCoa
#ifndef NDEBUG
if( boundaryHandling->isDomain(sx,sy,sz) )
WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy, sz, idx ) ) );
WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy, sz, idx ) ), sx << ", " << sy << ", " << sz << ", " << idx << " coarse sender block = " << coarseSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
#endif
}
rx += cell_idx_t(2);
......@@ -697,7 +697,7 @@ void PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::communicateLocalCoa
#ifndef NDEBUG
if( boundaryHandling->isDomain(sx,sy,sz) )
WALBERLA_ASSERT( !math::isnan( value ) );
WALBERLA_ASSERT( !math::isnan( value ), "value at " << sx << ", " << sy << ", " << sz << ", " << idx << " coarse sender block = " << coarseSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
#endif
rf->get( rx, ry, rz, idx ) = value;
......@@ -779,16 +779,16 @@ void PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::packDataFineToCoars
#ifndef NDEBUG
if( boundaryHandling->isDomain(x,y,z) )
{
WALBERLA_ASSERT( !math::isnan( field->get( x, y, z, idx ) ) );
WALBERLA_ASSERT( !math::isnan( field->get( x + cell_idx_t(1), y, z, idx ) ) );
WALBERLA_ASSERT( !math::isnan( field->get( x, y + cell_idx_t(1), z, idx ) ) );
WALBERLA_ASSERT( !math::isnan( field->get( x + cell_idx_t(1), y + cell_idx_t(1), z, idx ) ) );
WALBERLA_ASSERT( !math::isnan( field->get( x, y, z, idx ) ), x << ", " << y << ", " << z << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
WALBERLA_ASSERT( !math::isnan( field->get( x + cell_idx_t(1), y, z, idx ) ), x + cell_idx_t(1) << ", " << y << ", " << z << ", " << idx <<" fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
WALBERLA_ASSERT( !math::isnan( field->get( x, y + cell_idx_t(1), z, idx ) ), x << ", " << y + cell_idx_t(1) << ", " << z << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
WALBERLA_ASSERT( !math::isnan( field->get( x + cell_idx_t(1), y + cell_idx_t(1), z, idx ) ), x + cell_idx_t(1) << ", " << y + cell_idx_t(1) << ", " << z << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
if( Stencil::D == uint_t(3) )
{
WALBERLA_ASSERT( !math::isnan( field->get( x, y, z + cell_idx_t(1), idx ) ) );
WALBERLA_ASSERT( !math::isnan( field->get( x + cell_idx_t(1), y, z + cell_idx_t(1), idx ) ) );
WALBERLA_ASSERT( !math::isnan( field->get( x, y + cell_idx_t(1), z + cell_idx_t(1), idx ) ) );
WALBERLA_ASSERT( !math::isnan( field->get( x + cell_idx_t(1), y + cell_idx_t(1), z + cell_idx_t(1), idx ) ) );
WALBERLA_ASSERT( !math::isnan( field->get( x, y, z + cell_idx_t(1), idx ) ), x << ", " << y << ", " << z + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
WALBERLA_ASSERT( !math::isnan( field->get( x + cell_idx_t(1), y, z + cell_idx_t(1), idx ) ), x + cell_idx_t(1) << ", " << y << ", " << z + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
WALBERLA_ASSERT( !math::isnan( field->get( x, y + cell_idx_t(1), z + cell_idx_t(1), idx ) ), x << ", " << y + cell_idx_t(1) << ", " << z + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
WALBERLA_ASSERT( !math::isnan( field->get( x + cell_idx_t(1), y + cell_idx_t(1), z + cell_idx_t(1), idx ) ), x + cell_idx_t(1) << ", " << y + cell_idx_t(1) << ", " << z + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> " );
}
}
#endif
......@@ -925,16 +925,16 @@ void PdfFieldPackInfo< LatticeModel_T, BoundaryHandling_T >::communicateLocalFin
#ifndef NDEBUG
if( boundaryHandling->isDomain(sx,sy,sz) )
{
WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy, sz, idx ) ) );
WALBERLA_ASSERT( !math::isnan( sf->get( sx + cell_idx_t(1), sy, sz, idx ) ) );
WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy + cell_idx_t(1), sz, idx ) ) );
WALBERLA_ASSERT( !math::isnan( sf->get( sx + cell_idx_t(1), sy + cell_idx_t(1), sz, idx ) ) );
WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy, sz, idx ) ), sx << ", " << sy << ", " << sz << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
WALBERLA_ASSERT( !math::isnan( sf->get( sx + cell_idx_t(1), sy, sz, idx ) ), sx + cell_idx_t(1) << ", " << sy << ", " << sz << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy + cell_idx_t(1), sz, idx ) ), sx << ", " << sy + cell_idx_t(1) << ", " << sz << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
WALBERLA_ASSERT( !math::isnan( sf->get( sx + cell_idx_t(1), sy + cell_idx_t(1), sz, idx ) ), sx + cell_idx_t(1) << ", " << sy + cell_idx_t(1) << ", " << sz << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
if( Stencil::D == uint_t(3) )
{
WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy, sz + cell_idx_t(1), idx ) ) );
WALBERLA_ASSERT( !math::isnan( sf->get( sx + cell_idx_t(1), sy, sz + cell_idx_t(1), idx ) ) );
WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy + cell_idx_t(1), sz + cell_idx_t(1), idx ) ) );
WALBERLA_ASSERT( !math::isnan( sf->get( sx + cell_idx_t(1), sy + cell_idx_t(1), sz + cell_idx_t(1), idx ) ) );
WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy, sz + cell_idx_t(1), idx ) ), sx << ", " << sy << ", " << sz + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
WALBERLA_ASSERT( !math::isnan( sf->get( sx + cell_idx_t(1), sy, sz + cell_idx_t(1), idx ) ), sx + cell_idx_t(1) << ", " << sy << ", " << sz + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
WALBERLA_ASSERT( !math::isnan( sf->get( sx, sy + cell_idx_t(1), sz + cell_idx_t(1), idx ) ), sx << ", " << sy + cell_idx_t(1) << ", " << sz + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
WALBERLA_ASSERT( !math::isnan( sf->get( sx + cell_idx_t(1), sy + cell_idx_t(1), sz + cell_idx_t(1), idx ) ), sx + cell_idx_t(1) << ", " << sy + cell_idx_t(1) << ", " << sz + cell_idx_t(1) << ", " << idx << " fine sender block = " << fineSender->getId() << " in dir <" << stencil::cx[dir] << "," << stencil::cy[dir] << "," << stencil::cz[dir] << "> ");
}
}
#endif
......
//======================================================================================================================
//
// 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 RefinementFunctorWrapper.h
//! \ingroup lbm
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#pragma once
#include "domain_decomposition/BlockStorage.h"
#include <boost/function.hpp>
namespace walberla {
namespace lbm {
namespace refinement {
class FunctorWrapper {
public:
FunctorWrapper(boost::function<void()> fct)
: fct_(fct) {
}
void operator()(const uint_t /*level*/, const uint_t /*executionCounter*/) {
fct_();
}
private:
boost::function<void(void)> fct_;
};
class SweepAsFunctorWrapper {
public:
SweepAsFunctorWrapper( boost::function<void(IBlock * )> fct, const shared_ptr <StructuredBlockStorage> &blockStorage )
: fct_(fct), blockStorage_(blockStorage) {
}
void operator()(const uint_t level, const uint_t /*executionCounter*/) {
for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt) {
if (blockStorage_->getLevel(*blockIt) != level) continue;
fct_(blockIt.get());
}
}
private:
boost::function<void(IBlock * )> fct_;
shared_ptr <StructuredBlockStorage> blockStorage_;
};
} // namespace refinement
} // namespace lbm
} // namespace walberla
......@@ -1591,91 +1591,130 @@ void TimeStep< LatticeModel_T, Sweep_T, BoundaryHandling_T >::recursiveStep( con
std::vector< Block * > blocks = selectedBlocks( level );
WALBERLA_LOG_DETAIL("Starting recursive step with level " << level << " and execution count " << executionCount);
if( asynchronousCommunication_ && level != coarsestLevel )
{
WALBERLA_LOG_DETAIL("Start communication coarse to fine, initiated by fine level " << level );
startCommunicationCoarseToFine( level ); // [start] explosion (initiated by fine level, involves "level" and "level-1")
}
WALBERLA_LOG_DETAIL("Colliding on level " << level);
collide( blocks, level, executionCount1st );
if( asynchronousCommunication_ )
{
WALBERLA_LOG_DETAIL("Start communication equal level, initiated by level " << level );
startCommunicationEqualLevel( level ); // [start] equal level communication
}
if( level != finestLevel )
{
WALBERLA_LOG_DETAIL("Calling recursive step with level " << level + uint_t(1) << " and execution count " << executionCount1st );
recursiveStep( level + uint_t(1), executionCount1st );
if( asynchronousCommunication_ )
startCommunicationFineToCoarse( level + uint_t(1) ); // [start] coalescence (initiated by coarse level)
if( asynchronousCommunication_ ) {
WALBERLA_LOG_DETAIL("Start communication fine to coarse, initiated by coarse level " << level );
startCommunicationFineToCoarse(level + uint_t(1)); // [start] coalescence (initiated by coarse level)
}
}
if( level != coarsestLevel )
{
if( !asynchronousCommunication_ )
startCommunicationCoarseToFine( level ); // [start] explosion (initiated by fine level, involves "level" and "level-1")
if( !asynchronousCommunication_ ) {
WALBERLA_LOG_DETAIL("Start communication coarse to fine, initiated by fine level " << level );
startCommunicationCoarseToFine(level); // [start] explosion (initiated by fine level, involves "level" and "level-1")
}
WALBERLA_LOG_DETAIL("End communication coarse to fine, initiated by fine level " << level );
endCommunicationCoarseToFine( level ); // [end] explosion (initiated by fine level, involves "level" and "level-1")
WALBERLA_LOG_DETAIL("Perform linear explosion on level " << level );
performLinearExplosion( blocks, level );
}
if( !asynchronousCommunication_ )
startCommunicationEqualLevel( level ); // [start] equal level communication
if( !asynchronousCommunication_ ) {
WALBERLA_LOG_DETAIL("Start communication equal level, initiated by level " << level );
startCommunicationEqualLevel(level); // [start] equal level communication
}
WALBERLA_LOG_DETAIL("End communication equal level, initiated by level " << level );
endCommunicationEqualLevel( level ); // [end] equal level communication
// performLinearExplosion( blocks, level ); // if equal level neighbor values are needed, linear explosion should be performed here
if( level == finestLevel && level != coarsestLevel )
{
WALBERLA_LOG_DETAIL("Stream + collide on level " << level );
streamCollide( blocks, level, executionCount1st );
}
else
{
WALBERLA_LOG_DETAIL("Stream on level " << level );
stream( blocks, level, executionCount1st );
if( level != finestLevel )
{
if( !asynchronousCommunication_ )
startCommunicationFineToCoarse( level + uint_t(1) ); // [start] coalescence (initiated by coarse level)
if( !asynchronousCommunication_ ) {
WALBERLA_LOG_DETAIL("Start communication fine to coarse, initiated by coarse level " << level );
startCommunicationFineToCoarse(level + uint_t(1)); // [start] coalescence (initiated by coarse level)
}
WALBERLA_LOG_DETAIL("End communication fine to coarse, initiated by coarse level " << level );
endCommunicationFineToCoarse( level + uint_t(1) ); // [end] coalescence (initiated by coarse level)
}
WALBERLA_LOG_DETAIL("Finish stream on level " << level );
finishStream( blocks, level, executionCount1st );
if( level == coarsestLevel )
{
WALBERLA_LOG_DETAIL("End recursive step on level " << level );
return;
}
WALBERLA_LOG_DETAIL("Colliding on level " << level);
collide( blocks, level, executionCount2nd );
}
if( asynchronousCommunication_ )
startCommunicationEqualLevel( level ); // [start] equal level communication
if( asynchronousCommunication_ ) {
WALBERLA_LOG_DETAIL("Start communication equal level, initiated by level " << level );
startCommunicationEqualLevel(level); // [start] equal level communication
}
if( level != finestLevel )
{
WALBERLA_LOG_DETAIL("Calling recursive step with level " << level + uint_t(1) << " and execution count " << executionCount2nd );
recursiveStep( level + uint_t(1), executionCount2nd );
if( asynchronousCommunication_ )
startCommunicationFineToCoarse( level + uint_t(1) ); // [start] coalescence (initiated by coarse level)
if( asynchronousCommunication_ ) {
WALBERLA_LOG_DETAIL("Start communication fine to coarse, initiated by coarse level " << level );
startCommunicationFineToCoarse(level + uint_t(1)); // [start] coalescence (initiated by coarse level)
}
}
if( !asynchronousCommunication_ )
startCommunicationEqualLevel( level ); // [start] equal level communication
if( !asynchronousCommunication_ ) {
WALBERLA_LOG_DETAIL("Start communication equal level, initiated by level " << level );
startCommunicationEqualLevel(level); // [start] equal level communication
}
WALBERLA_LOG_DETAIL("End communication equal level, initiated by level " << level );
endCommunicationEqualLevel( level ); // [end] equal level communication
WALBERLA_LOG_DETAIL("Stream on level " << level );
stream( blocks, level, executionCount2nd );
if( level != finestLevel )
{
if( !asynchronousCommunication_ )
{
WALBERLA_LOG_DETAIL("Start communication fine to coarse, initiated by coarse level " << level );
startCommunicationFineToCoarse( level + uint_t(1) ); // [start] coalescence (initiated by coarse level)
}
WALBERLA_LOG_DETAIL("End communication fine to coarse, initiated by coarse level " << level );
endCommunicationFineToCoarse( level + uint_t(1) ); // [end] coalescence (initiated by coarse level)
}
finishStream( blocks, level, executionCount2nd );
}
WALBERLA_LOG_DETAIL("Finish stream on level " << level );
finishStream( blocks, level, executionCount2nd );
WALBERLA_LOG_DETAIL("End recursive step on level " << level );
}
......
......@@ -25,6 +25,7 @@
#include "BoundarySetup.h"
#include "PdfFieldPackInfo.h"
#include "PdfFieldSyncPackInfo.h"
#include "RefinementFunctorWrapper.h"
#include "TimeStep.h"
#include "TimeStepPdfPackInfo.h"
#include "TimeTracker.h"
......@@ -28,7 +28,7 @@
namespace walberla {
namespace pe_coupling {
CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, const shared_ptr<StructuredBlockStorage> & blockStorage,
CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, StructuredBlockStorage & blockStorage,
const uint_t numberOfGhostLayersToInclude )
{
......@@ -38,16 +38,33 @@ CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, const
if( body->isFinite() )
{
blockStorage->getCellBBFromAABB( cellBB, body->getAABB(), blockStorage->getLevel(block) );
} else
blockStorage.getCellBBFromAABB( cellBB, body->getAABB(), blockStorage.getLevel(block) );
}
else
{
blockStorage->getCellBBFromAABB( cellBB, body->getAABB().getIntersection( blockStorage->getDomain() ), blockStorage->getLevel(block) );
// if body is infinite (global), its AABB is also infinite which then requires special treatment
auto level = blockStorage.getLevel(block);
const real_t dx = blockStorage.dx(level);
const real_t dy = blockStorage.dy(level);
const real_t dz = blockStorage.dz(level);
Vector3<real_t> aabbExtensionByGhostLayers(real_c(numberOfGhostLayersToInclude) * dx,
real_c(numberOfGhostLayersToInclude) * dy,
real_c(numberOfGhostLayersToInclude) * dz);
auto extendedBlockAABB = blockStorage.getAABB(block.getId()).getExtended( aabbExtensionByGhostLayers );
// intersect the infinte (global) body with the block AABB, extended by its ghost layers
// then determine the cell bounding box of the intersection
blockStorage.getCellBBFromAABB( cellBB, body->getAABB().getIntersection( extendedBlockAABB ), level );
// if infinte body does not intersect with the extended block AABB, return an empty interval
if( cellBB.empty() ) return CellInterval();
}
cellBB.xMin() -= cell_idx_t(1); cellBB.yMin() -= cell_idx_t(1); cellBB.zMin() -= cell_idx_t(1);
cellBB.xMax() += cell_idx_t(1); cellBB.yMax() += cell_idx_t(1); cellBB.zMax() += cell_idx_t(1);
CellInterval blockBB = blockStorage->getBlockCellBB( block );
CellInterval blockBB = blockStorage.getBlockCellBB( block );
cell_idx_t layers = cell_idx_c( numberOfGhostLayersToInclude );
......@@ -56,12 +73,10 @@ CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, const
cellBB.intersect( blockBB );
blockStorage->transformGlobalToBlockLocalCellInterval( cellBB, block );
blockStorage.transformGlobalToBlockLocalCellInterval( cellBB, block );
return cellBB;
}
} // namespace pe_coupling
} // namespace walberla
......@@ -32,7 +32,7 @@ namespace walberla {
namespace pe_coupling {
CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, const shared_ptr<StructuredBlockStorage> & blockStorage,
CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, StructuredBlockStorage & blockStorage,
const uint_t numberOfGhostLayersToInclude );
......
......@@ -37,24 +37,16 @@
namespace walberla {
namespace pe_coupling {
// general mapping functions for a given single body on a given single block
template< typename BoundaryHandling_T >
void mapBody( const pe::BodyID & body, IBlock & block, const shared_ptr<StructuredBlockStorage> & blockStorage,
const BlockDataID & boundaryHandlingID, const FlagUID & obstacle,
const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
void mapBody( const pe::BodyID & body, IBlock & block, StructuredBlockStorage & blockStorage,
BoundaryHandling_T * boundaryHandlingPtr, const FlagUID & obstacle )
{
WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage->getBlockStorage()) );
WALBERLA_UNUSED(moBodiesOnly); // undo when other coupling algorithms are available
WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage.getBlockStorage()) );
WALBERLA_ASSERT_NOT_NULLPTR( boundaryHandlingPtr );
if( ( fixedBodiesOnly && !body->isFixed() ) /*|| ( moBodiesOnly && !isMOBody( body ) )*/ )
return;
BoundaryHandling_T * boundaryHandling = block.getData< BoundaryHandling_T >( boundaryHandlingID );
auto * flagField = boundaryHandling->getFlagField();
auto * flagField = boundaryHandlingPtr->getFlagField();
WALBERLA_ASSERT_NOT_NULLPTR( boundaryHandling );
WALBERLA_ASSERT_NOT_NULLPTR( flagField );
WALBERLA_ASSERT( flagField->flagExists( obstacle ) );
......@@ -62,10 +54,12 @@ void mapBody( const pe::BodyID & body, IBlock & block, const shared_ptr<Structur
CellInterval cellBB = getCellBB( body, block, blockStorage, flagField->nrOfGhostLayers() );
Vector3<real_t> startCellCenter = blockStorage->getBlockLocalCellCenter( block, cellBB.min() );
const real_t dx = blockStorage->dx( blockStorage->getLevel(block) );
const real_t dy = blockStorage->dy( blockStorage->getLevel(block) );
const real_t dz = blockStorage->dz( blockStorage->getLevel(block) );
if( cellBB.empty() ) return;
Vector3<real_t> startCellCenter = blockStorage.getBlockLocalCellCenter( block, cellBB.min() );
const real_t dx = blockStorage.dx( blockStorage.getLevel(block) );
const real_t dy = blockStorage.dy( blockStorage.getLevel(block) );
const real_t dz = blockStorage.dz( blockStorage.getLevel(block) );
real_t cz = startCellCenter[2];
for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
......@@ -77,58 +71,119 @@ void mapBody( const pe::BodyID & body, IBlock & block, const shared_ptr<Structur
for( cell_idx_t x = cellBB.xMin(); x <= cellBB.xMax(); ++x )
{
if( body->containsPoint(cx,cy,cz) )
boundaryHandling->forceBoundary( obstacleFlag, x, y, z );
boundaryHandlingPtr->forceBoundary( obstacleFlag, x, y, z );
cx += dx;
}
cy += dy;
}
cz += dz;
}
}
template< typename BoundaryHandling_T >
void mapBody( const pe::BodyID & body, IBlock & block, StructuredBlockStorage & blockStorage,
const BlockDataID & boundaryHandlingID, const FlagUID & obstacle )
{
WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage.getBlockStorage()) );
BoundaryHandling_T * boundaryHandling = block.getData< BoundaryHandling_T >( boundaryHandlingID );
mapBody(body, block, blockStorage, boundaryHandling, obstacle );
}
// mapping function to map all bodies from the body storage - with certain properties - onto all blocks
template< typename BoundaryHandling_T >
void mapBodies( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID, const BlockDataID & bodyStorageID, const FlagUID & obstacle,
void mapBodies( StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
const BlockDataID & bodyStorageID, const FlagUID & obstacle,
const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
{
for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
{
for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle, fixedBodiesOnly, moBodiesOnly );
{
WALBERLA_UNUSED(moBodiesOnly); // undo when other coupling algorithms are available
if( fixedBodiesOnly && bodyIt->isFixed() )
continue;
mapBody<BoundaryHandling_T>( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle );
}
}
}
// mapping function to map all global bodies - with certain properties - onto all blocks
template< typename BoundaryHandling_T >
void mapGlobalBodies( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID, pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle,
void mapGlobalBodies( StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle,
const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
{
for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
{
for( auto bodyIt = globalBodyStorage.begin(); bodyIt != globalBodyStorage.end(); ++bodyIt)
for( auto bodyIt = globalBodyStorage.begin(); bodyIt != globalBodyStorage.end(); ++bodyIt )
{
mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle, fixedBodiesOnly, moBodiesOnly );
WALBERLA_UNUSED(moBodiesOnly); // undo when other coupling algorithms are available
if( fixedBodiesOnly && bodyIt->isFixed() )
continue;
mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle );
}
}
}
// mapping function to map a given single global body onto all blocks
template< typename BoundaryHandling_T >
void mapGlobalBody( const id_t globalBodySystemID,
const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID, pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle,
const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
StructuredBlockStorage & blockStorage, const BlockDataID & boundaryHandlingID,
pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle )
{
for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
for( auto blockIt = blockStorage.begin(); blockIt != blockStorage.end(); ++blockIt )
{
auto bodyIt = globalBodyStorage.find( globalBodySystemID );
if( bodyIt != globalBodyStorage.end() )
{
mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle, fixedBodiesOnly, moBodiesOnly );
mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle );
}
}
}
// mapping function to map all global bodies - with certain properties - onto a given single block
template< typename BoundaryHandling_T >
void mapGlobalBodiesOnBlock( IBlock & block,
StructuredBlockStorage & blockStorage, BoundaryHandling_T * boundaryHandlingPtr,
pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle,
const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
{
for( auto bodyIt = globalBodyStorage.begin(); bodyIt != globalBodyStorage.end(); ++bodyIt)
{
WALBERLA_UNUSED(moBodiesOnly); // undo when other coupling algorithms are available
if( fixedBodiesOnly && bodyIt->isFixed() )
continue;
mapBody< BoundaryHandling_T >( *bodyIt, block, blockStorage, boundaryHandlingPtr, obstacle );
}
}
// mapping function to map a given single global body onto a given single block
template< typename BoundaryHandling_T >
void mapGlobalBodyOnBlock( const id_t globalBodySystemID, IBlock & block,
StructuredBlockStorage & blockStorage, BoundaryHandling_T * boundaryHandlingPtr,
pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle )
{
auto bodyIt = globalBodyStorage.find( globalBodySystemID );
if( bodyIt != globalBodyStorage.end() )
{
mapBody< BoundaryHandling_T >( *bodyIt, block, blockStorage, boundaryHandlingPtr, obstacle );
}
}
} // namespace pe_coupling
} // namespace walberla
......@@ -108,7 +108,7 @@ void BodyMapping< BoundaryHandling_T >::operator()( IBlock * const block ) const
// policy: every body manages only its own flags
CellInterval cellBB = getCellBB( *bodyIt, *block, blockStorage_, flagField->nrOfGhostLayers() );
CellInterval cellBB = getCellBB( *bodyIt, *block, *blockStorage_, flagField->nrOfGhostLayers() );
Vector3<real_t> startCellCenter = blockStorage_->getBlockLocalCellCenter( *block, cellBB.min() );
const real_t dx = blockStorage_->dx( blockStorage_->getLevel(*block) );
......@@ -174,12 +174,12 @@ void BodyMapping< BoundaryHandling_T >::operator()( IBlock * const block ) const
////////////////////////////
template< typename BoundaryHandling_T >
void mapMovingBody( const pe::BodyID body, IBlock & block, const shared_ptr<StructuredBlockStorage> & blockStorage,
void mapMovingBody( const pe::BodyID body, IBlock & block, StructuredBlockStorage & blockStorage,
const BlockDataID & boundaryHandlingID, const BlockDataID & bodyFieldID, const FlagUID & obstacle )
{
typedef Field< pe::BodyID, 1 > BodyField_T;
WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage->getBlockStorage()) );
WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage.getBlockStorage()) );
if( body->isFixed() /*|| !body->isFinite()*/ )
return;
......@@ -198,10 +198,12 @@ void mapMovingBody( const pe::BodyID body, IBlock & block, const shared_ptr<Stru