diff --git a/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.cpp b/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9821dbb80e8e9565c889ee61be52b83a1a733e14
--- /dev/null
+++ b/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.cpp
@@ -0,0 +1,49 @@
+//======================================================================================================================
+//
+// 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 ExtrapolationDirectionFinder.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+
+#include "ExtrapolationDirectionFinder.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+
+void SphereNormalExtrapolationDirectionFinder
+::getDirection( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, Vector3<cell_idx_t> & extrapolationDirection ) const
+{
+ BodyField_T * bodyField = block->getData< BodyField_T >( bodyFieldID_ );
+
+ WALBERLA_ASSERT_NOT_NULLPTR( bodyField );
+ WALBERLA_ASSERT_NOT_NULLPTR( (*bodyField)(x,y,z) );
+
+ real_t cx, cy, cz;
+ blockStorage_->getBlockLocalCellCenter( *block, Cell(x,y,z), cx, cy, cz );
+
+ Vector3<real_t> bodyCenterPosition = (*bodyField)(x,y,z)->getPosition();
+ WALBERLA_ASSERT( !math::isnan(bodyCenterPosition) );
+
+ Vector3<real_t> bodyNormal( cx - bodyCenterPosition[0], cy - bodyCenterPosition[1], cz - bodyCenterPosition[2] );
+
+ findCorrespondingLatticeDirection< stencil::D3Q27 >( bodyNormal, extrapolationDirection );
+}
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.h b/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.h
index e23d38be85023d5d2976b62d0a995cb09a517de5..7e9968d83b16e7070931e6be79f7e79a92efba53 100644
--- a/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.h
+++ b/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.h
@@ -134,24 +134,5 @@ private:
};
-void SphereNormalExtrapolationDirectionFinder
-::getDirection( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, Vector3<cell_idx_t> & extrapolationDirection ) const
-{
- BodyField_T * bodyField = block->getData< BodyField_T >( bodyFieldID_ );
-
- WALBERLA_ASSERT_NOT_NULLPTR( bodyField );
- WALBERLA_ASSERT_NOT_NULLPTR( (*bodyField)(x,y,z) );
-
- real_t cx, cy, cz;
- blockStorage_->getBlockLocalCellCenter( *block, Cell(x,y,z), cx, cy, cz );
-
- Vector3<real_t> bodyCenterPosition = (*bodyField)(x,y,z)->getPosition();
- WALBERLA_ASSERT( !math::isnan(bodyCenterPosition) );
-
- Vector3<real_t> bodyNormal( cx - bodyCenterPosition[0], cy - bodyCenterPosition[1], cz - bodyCenterPosition[2] );
-
- findCorrespondingLatticeDirection< stencil::D3Q27 >( bodyNormal, extrapolationDirection );
-}
-
} // namespace pe_coupling
} // namespace walberla
diff --git a/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.cpp b/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..75ea1ed2f3f6870152d54fc0b921e1219fcc8a4a
--- /dev/null
+++ b/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.cpp
@@ -0,0 +1,235 @@
+//======================================================================================================================
+//
+// 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 BodyAndVolumeFractionMapping.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "BodyAndVolumeFractionMapping.h"
+
+#include "geometry/bodies/BodyOverlapFunctions.h"
+
+#include "pe/rigidbody/BodyIterators.h"
+
+#include "pe_coupling/geometry/PeOverlapFraction.h"
+#include "pe_coupling/mapping/BodyBBMapping.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+
+void mapPSMBodyAndVolumeFraction( const pe::BodyID body, IBlock & block, StructuredBlockStorage & blockStorage,
+ const BlockDataID bodyAndVolumeFractionFieldID )
+{
+ typedef std::pair< pe::BodyID, real_t > BodyAndVolumeFraction_T;
+ typedef GhostLayerField< std::vector< BodyAndVolumeFraction_T >, 1 > BodyAndVolumeFractionField_T;
+
+ BodyAndVolumeFractionField_T * bodyAndVolumeFractionField = block.getData< BodyAndVolumeFractionField_T >( bodyAndVolumeFractionFieldID );
+ WALBERLA_ASSERT_NOT_NULLPTR( bodyAndVolumeFractionField );
+
+ // get bounding box of body
+ CellInterval cellBB = getCellBB( body, block, blockStorage, bodyAndVolumeFractionField->nrOfGhostLayers() );
+
+ uint_t level = blockStorage.getLevel( block );
+ Vector3<real_t> dxVec(blockStorage.dx(level), blockStorage.dy(level), blockStorage.dz(level));
+
+ for( auto cellIt = cellBB.begin(); cellIt != cellBB.end(); ++cellIt )
+ {
+ Cell cell( *cellIt );
+
+ // get the cell's center
+ Vector3<real_t> cellCenter;
+ cellCenter = blockStorage.getBlockLocalCellCenter( block, cell );
+
+ const real_t fraction = overlapFractionPe( *body, cellCenter, dxVec );
+
+ // if the cell intersected with the body, store a pointer to that body and the corresponding volume fraction in the field
+ if( fraction > real_t(0) )
+ {
+ bodyAndVolumeFractionField->get(cell).push_back( BodyAndVolumeFraction_T( body, fraction ) );
+ }
+ }
+}
+
+
+void BodyAndVolumeFractionMapping::initialize()
+{
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ BodyAndVolumeFractionField_T * bodyAndVolumeFractionField = blockIt->getData< BodyAndVolumeFractionField_T >( bodyAndVolumeFractionFieldID_ );
+
+ if( updatedBodyAndVolumeFractionField_ == NULL )
+ {
+ // hold internally an identical field for swapping
+ updatedBodyAndVolumeFractionField_ = shared_ptr<BodyAndVolumeFractionField_T>( bodyAndVolumeFractionField->cloneUninitialized() );
+
+ auto xyzFieldSize = updatedBodyAndVolumeFractionField_->xyzSize();
+ for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
+ {
+ (updatedBodyAndVolumeFractionField_->get(*fieldIt)).clear();
+ }
+ }
+
+ // clear the field
+ auto xyzFieldSize = bodyAndVolumeFractionField->xyzSize();
+ for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
+ {
+ (bodyAndVolumeFractionField->get(*fieldIt)).clear();
+ }
+
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
+ {
+ mapPSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, *blockStorage_, bodyAndVolumeFractionFieldID_ );
+ lastUpdatedPositionMap_.insert( std::pair< walberla::id_t, Vector3< real_t > >( bodyIt->getSystemID(), bodyIt->getPosition() ) );
+ }
+ }
+
+ for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
+ {
+ if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
+ {
+ mapPSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, *blockStorage_, bodyAndVolumeFractionFieldID_);
+ lastUpdatedPositionMap_.insert( std::pair< walberla::id_t, Vector3< real_t > >( bodyIt->getSystemID(), bodyIt->getPosition() ) );
+ }
+ }
+
+ }
+}
+
+void BodyAndVolumeFractionMapping::update()
+{
+ std::map< walberla::id_t, Vector3< real_t > > tempLastUpdatedPositionMap;
+
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ BodyAndVolumeFractionField_T * bodyAndVolumeFractionField = blockIt->getData< BodyAndVolumeFractionField_T >( bodyAndVolumeFractionFieldID_ );
+
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
+ {
+ updatePSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, bodyAndVolumeFractionField, tempLastUpdatedPositionMap);
+ }
+ }
+
+ for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
+ {
+ if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
+ {
+ updatePSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, bodyAndVolumeFractionField, tempLastUpdatedPositionMap);
+ }
+ }
+
+ bodyAndVolumeFractionField->swapDataPointers( *updatedBodyAndVolumeFractionField_ );
+
+ auto xyzFieldSize = updatedBodyAndVolumeFractionField_->xyzSize();
+ for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
+ {
+ (updatedBodyAndVolumeFractionField_->get(*fieldIt)).clear();
+ }
+ }
+
+ lastUpdatedPositionMap_ = tempLastUpdatedPositionMap;
+}
+
+void BodyAndVolumeFractionMapping::updatePSMBodyAndVolumeFraction( pe::BodyID body, IBlock & block,
+ BodyAndVolumeFractionField_T * oldBodyAndVolumeFractionField,
+ std::map< walberla::id_t, Vector3< real_t > > & tempLastUpdatedPositionMap )
+{
+
+ WALBERLA_ASSERT_NOT_NULLPTR( oldBodyAndVolumeFractionField );
+
+ // estimate traveled distance since last volume fraction update
+ real_t traveledSquaredDistance( real_t(0) );
+ auto mapBodyIt = lastUpdatedPositionMap_.find( body->getSystemID() );
+ if( mapBodyIt != lastUpdatedPositionMap_.end() )
+ {
+ // body found and traveled distance can be estimated
+ Vector3<real_t> distanceVec = body->getPosition() - mapBodyIt->second;
+ traveledSquaredDistance = distanceVec.sqrLength();
+ } else
+ {
+ // body was not found in map -> is a new body, so no information is known
+ traveledSquaredDistance = std::numeric_limits<real_t>::max();
+ }
+
+ // get bounding box of body
+ CellInterval cellBB = getCellBB( body, block, *blockStorage_, oldBodyAndVolumeFractionField->nrOfGhostLayers() );
+
+ uint_t level = blockStorage_->getLevel( block );
+ Vector3<real_t> dxVec(blockStorage_->dx(level), blockStorage_->dy(level), blockStorage_->dz(level));
+
+ // if body has not moved (specified by some epsilon), just reuse old fraction values
+ if( body->getLinearVel().sqrLength() < velocityUpdatingEpsilonSquared_ &&
+ body->getAngularVel().sqrLength() < velocityUpdatingEpsilonSquared_ &&
+ traveledSquaredDistance < positionUpdatingEpsilonSquared_ )
+ {
+ for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
+ {
+ for( cell_idx_t y = cellBB.yMin(); y <= cellBB.yMax(); ++y )
+ {
+ for( cell_idx_t x = cellBB.xMin(); x <= cellBB.xMax(); ++x )
+ {
+
+ auto oldVec = oldBodyAndVolumeFractionField->get(x,y,z);
+ for( auto pairIt = oldVec.begin(); pairIt != oldVec.end(); ++pairIt )
+ {
+ if( pairIt->first == body )
+ {
+ updatedBodyAndVolumeFractionField_->get(x,y,z).push_back(*pairIt);
+ break;
+ }
+ }
+ }
+ }
+ }
+ tempLastUpdatedPositionMap.insert( std::pair< walberla::id_t, Vector3< real_t > >( body->getSystemID(), mapBodyIt->second ) );
+
+ } else
+ {
+ // else body has moved significantly or is new to this block, thus the volume fraction has to be calculated
+ for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
+ {
+ for( cell_idx_t y = cellBB.yMin(); y <= cellBB.yMax(); ++y )
+ {
+ for( cell_idx_t x = cellBB.xMin(); x <= cellBB.xMax(); ++x )
+ {
+ // get the cell's center
+ Vector3<real_t> cellCenter;
+ cellCenter = blockStorage_->getBlockLocalCellCenter( block, Cell(x,y,z) );
+
+ const real_t fraction = overlapFractionPe( *body, cellCenter, dxVec, superSamplingDepth_ );
+
+ // if the cell intersected with the body, store a pointer to that body and the corresponding volume fraction in the field
+ if( fraction > real_t(0) )
+ {
+ updatedBodyAndVolumeFractionField_->get(x,y,z).push_back( BodyAndVolumeFraction_T( body, fraction ) );
+ }
+ }
+ }
+ }
+
+ tempLastUpdatedPositionMap.insert( std::pair< walberla::id_t, Vector3< real_t > >( body->getSystemID(), body->getPosition() ) );
+ }
+}
+
+
+} // namespace pe_coupling
+} // namespace walberla
+
diff --git a/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h b/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h
index 5f59d9754cb9e3edd911d1127d912244f5bdfae2..59fb85590fd9db424523e3f2c6cf6db0d9d84bef 100644
--- a/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h
+++ b/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h
@@ -22,16 +22,8 @@
#pragma once
#include "domain_decomposition/StructuredBlockStorage.h"
-
#include "field/GhostLayerField.h"
-
-#include "geometry/bodies/BodyOverlapFunctions.h"
-
-#include "pe/rigidbody/BodyIterators.h"
#include "pe/Types.h"
-
-#include "pe_coupling/geometry/PeOverlapFraction.h"
-#include "pe_coupling/mapping/BodyBBMapping.h"
#include "pe_coupling/utility/BodySelectorFunctions.h"
#include <functional>
@@ -49,37 +41,7 @@ namespace pe_coupling {
*
*/
void mapPSMBodyAndVolumeFraction( const pe::BodyID body, IBlock & block, StructuredBlockStorage & blockStorage,
- const BlockDataID bodyAndVolumeFractionFieldID )
-{
- typedef std::pair< pe::BodyID, real_t > BodyAndVolumeFraction_T;
- typedef GhostLayerField< std::vector< BodyAndVolumeFraction_T >, 1 > BodyAndVolumeFractionField_T;
-
- BodyAndVolumeFractionField_T * bodyAndVolumeFractionField = block.getData< BodyAndVolumeFractionField_T >( bodyAndVolumeFractionFieldID );
- WALBERLA_ASSERT_NOT_NULLPTR( bodyAndVolumeFractionField );
-
- // get bounding box of body
- CellInterval cellBB = getCellBB( body, block, blockStorage, bodyAndVolumeFractionField->nrOfGhostLayers() );
-
- uint_t level = blockStorage.getLevel( block );
- Vector3<real_t> dxVec(blockStorage.dx(level), blockStorage.dy(level), blockStorage.dz(level));
-
- for( auto cellIt = cellBB.begin(); cellIt != cellBB.end(); ++cellIt )
- {
- Cell cell( *cellIt );
-
- // get the cell's center
- Vector3<real_t> cellCenter;
- cellCenter = blockStorage.getBlockLocalCellCenter( block, cell );
-
- const real_t fraction = overlapFractionPe( *body, cellCenter, dxVec );
-
- // if the cell intersected with the body, store a pointer to that body and the corresponding volume fraction in the field
- if( fraction > real_t(0) )
- {
- bodyAndVolumeFractionField->get(cell).push_back( BodyAndVolumeFraction_T( body, fraction ) );
- }
- }
-}
+ const BlockDataID bodyAndVolumeFractionFieldID );
/*!\brief Mapping class that can be used inside the timeloop to update the volume fractions and body mappings
*
@@ -154,168 +116,6 @@ private:
const uint_t superSamplingDepth_;
};
-void BodyAndVolumeFractionMapping::initialize()
-{
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- BodyAndVolumeFractionField_T * bodyAndVolumeFractionField = blockIt->getData< BodyAndVolumeFractionField_T >( bodyAndVolumeFractionFieldID_ );
-
- if( updatedBodyAndVolumeFractionField_ == NULL )
- {
- // hold internally an identical field for swapping
- updatedBodyAndVolumeFractionField_ = shared_ptr<BodyAndVolumeFractionField_T>( bodyAndVolumeFractionField->cloneUninitialized() );
-
- auto xyzFieldSize = updatedBodyAndVolumeFractionField_->xyzSize();
- for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
- {
- (updatedBodyAndVolumeFractionField_->get(*fieldIt)).clear();
- }
- }
-
- // clear the field
- auto xyzFieldSize = bodyAndVolumeFractionField->xyzSize();
- for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
- {
- (bodyAndVolumeFractionField->get(*fieldIt)).clear();
- }
-
- for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
- {
- mapPSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, *blockStorage_, bodyAndVolumeFractionFieldID_ );
- lastUpdatedPositionMap_.insert( std::pair< walberla::id_t, Vector3< real_t > >( bodyIt->getSystemID(), bodyIt->getPosition() ) );
- }
- }
-
- for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
- {
- if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
- {
- mapPSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, *blockStorage_, bodyAndVolumeFractionFieldID_);
- lastUpdatedPositionMap_.insert( std::pair< walberla::id_t, Vector3< real_t > >( bodyIt->getSystemID(), bodyIt->getPosition() ) );
- }
- }
-
- }
-}
-
-void BodyAndVolumeFractionMapping::update()
-{
- std::map< walberla::id_t, Vector3< real_t > > tempLastUpdatedPositionMap;
-
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- BodyAndVolumeFractionField_T * bodyAndVolumeFractionField = blockIt->getData< BodyAndVolumeFractionField_T >( bodyAndVolumeFractionFieldID_ );
-
- for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
- {
- updatePSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, bodyAndVolumeFractionField, tempLastUpdatedPositionMap);
- }
- }
-
- for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
- {
- if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
- {
- updatePSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, bodyAndVolumeFractionField, tempLastUpdatedPositionMap);
- }
- }
-
- bodyAndVolumeFractionField->swapDataPointers( *updatedBodyAndVolumeFractionField_ );
-
- auto xyzFieldSize = updatedBodyAndVolumeFractionField_->xyzSize();
- for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
- {
- (updatedBodyAndVolumeFractionField_->get(*fieldIt)).clear();
- }
- }
-
- lastUpdatedPositionMap_ = tempLastUpdatedPositionMap;
-}
-
-void BodyAndVolumeFractionMapping::updatePSMBodyAndVolumeFraction( pe::BodyID body, IBlock & block,
- BodyAndVolumeFractionField_T * oldBodyAndVolumeFractionField,
- std::map< walberla::id_t, Vector3< real_t > > & tempLastUpdatedPositionMap )
-{
-
- WALBERLA_ASSERT_NOT_NULLPTR( oldBodyAndVolumeFractionField );
-
- // estimate traveled distance since last volume fraction update
- real_t traveledSquaredDistance( real_t(0) );
- auto mapBodyIt = lastUpdatedPositionMap_.find( body->getSystemID() );
- if( mapBodyIt != lastUpdatedPositionMap_.end() )
- {
- // body found and traveled distance can be estimated
- Vector3<real_t> distanceVec = body->getPosition() - mapBodyIt->second;
- traveledSquaredDistance = distanceVec.sqrLength();
- } else
- {
- // body was not found in map -> is a new body, so no information is known
- traveledSquaredDistance = std::numeric_limits<real_t>::max();
- }
-
- // get bounding box of body
- CellInterval cellBB = getCellBB( body, block, *blockStorage_, oldBodyAndVolumeFractionField->nrOfGhostLayers() );
-
- uint_t level = blockStorage_->getLevel( block );
- Vector3<real_t> dxVec(blockStorage_->dx(level), blockStorage_->dy(level), blockStorage_->dz(level));
-
- // if body has not moved (specified by some epsilon), just reuse old fraction values
- if( body->getLinearVel().sqrLength() < velocityUpdatingEpsilonSquared_ &&
- body->getAngularVel().sqrLength() < velocityUpdatingEpsilonSquared_ &&
- traveledSquaredDistance < positionUpdatingEpsilonSquared_ )
- {
- for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
- {
- for( cell_idx_t y = cellBB.yMin(); y <= cellBB.yMax(); ++y )
- {
- for( cell_idx_t x = cellBB.xMin(); x <= cellBB.xMax(); ++x )
- {
-
- auto oldVec = oldBodyAndVolumeFractionField->get(x,y,z);
- for( auto pairIt = oldVec.begin(); pairIt != oldVec.end(); ++pairIt )
- {
- if( pairIt->first == body )
- {
- updatedBodyAndVolumeFractionField_->get(x,y,z).push_back(*pairIt);
- break;
- }
- }
- }
- }
- }
- tempLastUpdatedPositionMap.insert( std::pair< walberla::id_t, Vector3< real_t > >( body->getSystemID(), mapBodyIt->second ) );
-
- } else
- {
- // else body has moved significantly or is new to this block, thus the volume fraction has to be calculated
- for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
- {
- for( cell_idx_t y = cellBB.yMin(); y <= cellBB.yMax(); ++y )
- {
- for( cell_idx_t x = cellBB.xMin(); x <= cellBB.xMax(); ++x )
- {
- // get the cell's center
- Vector3<real_t> cellCenter;
- cellCenter = blockStorage_->getBlockLocalCellCenter( block, Cell(x,y,z) );
-
- const real_t fraction = overlapFractionPe( *body, cellCenter, dxVec, superSamplingDepth_ );
-
- // if the cell intersected with the body, store a pointer to that body and the corresponding volume fraction in the field
- if( fraction > real_t(0) )
- {
- updatedBodyAndVolumeFractionField_->get(x,y,z).push_back( BodyAndVolumeFraction_T( body, fraction ) );
- }
- }
- }
- }
-
- tempLastUpdatedPositionMap.insert( std::pair< walberla::id_t, Vector3< real_t > >( body->getSystemID(), body->getPosition() ) );
- }
-}
} // namespace pe_coupling
diff --git a/src/pe_coupling/utility/BodiesForceTorqueContainer.cpp b/src/pe_coupling/utility/BodiesForceTorqueContainer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bed1c267c35770e232638d6b423e283cf3bc9172
--- /dev/null
+++ b/src/pe_coupling/utility/BodiesForceTorqueContainer.cpp
@@ -0,0 +1,97 @@
+//======================================================================================================================
+//
+// 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 BodiesForceTorqueContainer.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "BodiesForceTorqueContainer.h"
+
+#include "blockforest/StructuredBlockForest.h"
+#include "core/math/Vector3.h"
+#include "pe/rigidbody/BodyIterators.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+
+
+void BodiesForceTorqueContainer::store()
+{
+ // clear map
+ clear();
+
+ // (re-)build map
+ for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
+ {
+ auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
+
+ for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ auto & f = (*bodyForceTorqueStorage)[ bodyIt->getSystemID() ];
+
+ const auto & force = bodyIt->getForce();
+ const auto & torque = bodyIt->getTorque();
+ f = {{force[0], force[1], force[2], torque[0], torque[1], torque[2] }};
+ }
+ }
+}
+
+void BodiesForceTorqueContainer::setOnBodies()
+{
+ // set the force/torque stored in the block-local map onto all bodies
+ for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
+ {
+ auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
+
+ for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ const auto f = bodyForceTorqueStorage->find( bodyIt->getSystemID() );
+
+ if( f != bodyForceTorqueStorage->end() )
+ {
+ const auto & ftValues = f->second;
+ bodyIt->addForce ( ftValues[0], ftValues[1], ftValues[2] );
+ bodyIt->addTorque( ftValues[3], ftValues[4], ftValues[5] );
+ }
+ // else: new body has arrived that was not known before
+ }
+ }
+}
+
+void BodiesForceTorqueContainer::clear()
+{
+ for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
+ {
+ auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
+ bodyForceTorqueStorage->clear();
+ }
+}
+
+void BodiesForceTorqueContainer::swap( BodiesForceTorqueContainer & other )
+{
+ std::swap( bodyForceTorqueStorageID_, other.bodyForceTorqueStorageID_);
+}
+
+
+
+
+
+
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/utility/BodiesForceTorqueContainer.h b/src/pe_coupling/utility/BodiesForceTorqueContainer.h
index d69b6791b23734526521e10a3b42a6acb8b5f6bf..fe367b5e095075f7882c19f9cf213d793be40351 100644
--- a/src/pe_coupling/utility/BodiesForceTorqueContainer.h
+++ b/src/pe_coupling/utility/BodiesForceTorqueContainer.h
@@ -22,9 +22,6 @@
#pragma once
#include "blockforest/StructuredBlockForest.h"
-#include "core/math/Vector3.h"
-#include "pe/rigidbody/BodyIterators.h"
-#include "pe/synchronization/SyncForces.h"
#include <map>
#include <array>
@@ -50,63 +47,13 @@ public:
store();
}
- void store()
- {
- // clear map
- clear();
-
- // (re-)build map
- for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
- {
- auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
-
- for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- auto & f = (*bodyForceTorqueStorage)[ bodyIt->getSystemID() ];
-
- const auto & force = bodyIt->getForce();
- const auto & torque = bodyIt->getTorque();
- f = {{force[0], force[1], force[2], torque[0], torque[1], torque[2] }};
- }
- }
-
- }
+ void store();
- void setOnBodies()
- {
- // set the force/torque stored in the block-local map onto all bodies
- for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
- {
- auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
-
- for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- const auto f = bodyForceTorqueStorage->find( bodyIt->getSystemID() );
-
- if( f != bodyForceTorqueStorage->end() )
- {
- const auto & ftValues = f->second;
- bodyIt->addForce ( ftValues[0], ftValues[1], ftValues[2] );
- bodyIt->addTorque( ftValues[3], ftValues[4], ftValues[5] );
- }
- // else: new body has arrived that was not known before
- }
- }
- }
+ void setOnBodies();
- void clear()
- {
- for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
- {
- auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
- bodyForceTorqueStorage->clear();
- }
- }
+ void clear();
- void swap( BodiesForceTorqueContainer & other )
- {
- std::swap( bodyForceTorqueStorageID_, other.bodyForceTorqueStorageID_);
- }
+ void swap( BodiesForceTorqueContainer & other );
private:
diff --git a/src/pe_coupling/utility/ForceOnBodiesAdder.cpp b/src/pe_coupling/utility/ForceOnBodiesAdder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4c7c17f8376763af77004bedb8cd83186515460e
--- /dev/null
+++ b/src/pe_coupling/utility/ForceOnBodiesAdder.cpp
@@ -0,0 +1,49 @@
+//======================================================================================================================
+//
+// 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 ForceOnBodiesAdder.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "ForceOnBodiesAdder.h"
+
+#include "core/math/Vector3.h"
+#include "domain_decomposition/StructuredBlockStorage.h"
+#include "pe/rigidbody/BodyIterators.h"
+
+
+namespace walberla {
+namespace pe_coupling {
+
+void ForceOnBodiesAdder::operator()()
+{
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->addForce ( force_ );
+ }
+ }
+}
+
+void ForceOnBodiesAdder::updateForce( const Vector3<real_t> & newForce )
+{
+ force_ = newForce;
+}
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/utility/ForceOnBodiesAdder.h b/src/pe_coupling/utility/ForceOnBodiesAdder.h
index 2d6e3f4ead588cc83b1fbc147e76fed8d4a1d38c..b0b7056d52aa54498197317eb48f4fdf294c4b0b 100644
--- a/src/pe_coupling/utility/ForceOnBodiesAdder.h
+++ b/src/pe_coupling/utility/ForceOnBodiesAdder.h
@@ -23,8 +23,6 @@
#include "core/math/Vector3.h"
#include "domain_decomposition/StructuredBlockStorage.h"
-#include "pe/rigidbody/BodyIterators.h"
-
namespace walberla {
namespace pe_coupling {
@@ -39,21 +37,9 @@ public:
{ }
// set a constant force on all (only local, to avoid force duplication) bodies
- void operator()()
- {
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
- {
- bodyIt->addForce ( force_ );
- }
- }
- }
+ void operator()();
- void updateForce( const Vector3<real_t> & newForce )
- {
- force_ = newForce;
- }
+ void updateForce( const Vector3<real_t> & newForce );
private:
diff --git a/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.cpp b/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e3301650a01ce5b2a5a6cfc7f517576891ddfaf4
--- /dev/null
+++ b/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.cpp
@@ -0,0 +1,41 @@
+//======================================================================================================================
+//
+// 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 ForceTorqueOnBodiesResetter.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "ForceTorqueOnBodiesResetter.h"
+
+#include "pe/rigidbody/BodyIterators.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+void ForceTorqueOnBodiesResetter::operator()()
+{
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->resetForceAndTorque();
+ }
+ }
+}
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.h b/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.h
index c42e96b0c9cd6bdbac95688654aa9d4b1b979e17..a86368f83fdaf26b4f2c586ac90f178042d567f6 100644
--- a/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.h
+++ b/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.h
@@ -21,10 +21,7 @@
#pragma once
-#include "core/math/Vector3.h"
#include "domain_decomposition/StructuredBlockStorage.h"
-#include "pe/rigidbody/BodyIterators.h"
-
namespace walberla {
namespace pe_coupling {
@@ -38,17 +35,7 @@ public:
{ }
// resets forces and torques on all (local and remote) bodies
- void operator()()
- {
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- bodyIt->resetForceAndTorque();
- }
- }
- }
-
+ void operator()();
private:
diff --git a/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.cpp b/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..677c5788d7dea105e6952d911d9d8942deb0cbce
--- /dev/null
+++ b/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.cpp
@@ -0,0 +1,55 @@
+//======================================================================================================================
+//
+// 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 ForceTorqueOnBodiesScaler.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "ForceTorqueOnBodiesScaler.h"
+
+#include "core/math/Vector3.h"
+#include "pe/rigidbody/BodyIterators.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+void ForceTorqueOnBodiesScaler::operator()()
+{
+ Vector3<real_t> force(real_t(0));
+ Vector3<real_t> torque(real_t(0));
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ force = scalingFactor_ * bodyIt->getForce();
+ torque = scalingFactor_ * bodyIt->getTorque();
+
+ bodyIt->resetForceAndTorque();
+
+ bodyIt->setForce ( force );
+ bodyIt->setTorque( torque );
+ }
+ }
+}
+
+void ForceTorqueOnBodiesScaler::resetScalingFactor( const real_t newScalingFactor )
+{
+ scalingFactor_ = newScalingFactor;
+}
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.h b/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.h
index ead115a69ea789ec629c600d396b7c2028d6fa7e..f903e8ddc1606112f60dbf19a872379077a4c325 100644
--- a/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.h
+++ b/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.h
@@ -21,15 +21,12 @@
#pragma once
-#include "core/math/Vector3.h"
#include "domain_decomposition/StructuredBlockStorage.h"
-#include "pe/rigidbody/BodyIterators.h"
-
namespace walberla {
namespace pe_coupling {
-// scales force/torquew on all bodies (local and remote) by a constant scalar value
+// scales force/torque on all bodies (local and remote) by a constant scalar value
// can e.g. be used to average the force/torque over two time steps
class ForceTorqueOnBodiesScaler
{
@@ -41,29 +38,9 @@ public:
{ }
// resets forces and torques on all (local and remote) bodies
- void operator()()
- {
- Vector3<real_t> force(real_t(0));
- Vector3<real_t> torque(real_t(0));
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- force = scalingFactor_ * bodyIt->getForce();
- torque = scalingFactor_ * bodyIt->getTorque();
-
- bodyIt->resetForceAndTorque();
-
- bodyIt->setForce ( force );
- bodyIt->setTorque( torque );
- }
- }
- }
+ void operator()();
- void resetScalingFactor( const real_t newScalingFactor )
- {
- scalingFactor_ = newScalingFactor;
- }
+ void resetScalingFactor( const real_t newScalingFactor );
private:
diff --git a/src/pe_coupling/utility/LubricationCorrection.cpp b/src/pe_coupling/utility/LubricationCorrection.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..eaf1952149fc4c35d7b6b45b23549c3516557f76
--- /dev/null
+++ b/src/pe_coupling/utility/LubricationCorrection.cpp
@@ -0,0 +1,278 @@
+//======================================================================================================================
+//
+// 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 LubricationCorrection.cpp
+//! \ingroup pe_coupling
+//! \author Kristina Pickl <kristina.pickl@fau.de>
+//! \author Dominik Bartuschat
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "LubricationCorrection.h"
+
+#include "core/logging/all.h"
+#include "core/debug/Debug.h"
+
+#include "pe/rigidbody/BodyIterators.h"
+#include "pe/utility/Distance.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+void LubricationCorrection::operator ()()
+{
+ WALBERLA_LOG_PROGRESS( "Calculating Lubrication Force" );
+
+ for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt)
+ {
+ // loop over all rigid bodies
+ for( auto body1It = pe::BodyIterator::begin( *blockIt, bodyStorageID_ ); body1It != pe::BodyIterator::end(); ++body1It )
+ {
+ // lubrication forces for spheres
+ if ( body1It->getTypeID() == pe::Sphere::getStaticTypeID() )
+ {
+ pe::SphereID sphereI = static_cast<pe::SphereID> ( body1It.getBodyID() );
+
+ auto copyBody1It = body1It;
+ // loop over all rigid bodies after current body1 to avoid double forces
+ for( auto body2It = (++copyBody1It); body2It != pe::BodyIterator::end(); ++body2It )
+ {
+ // sphere-sphere lubrication
+ if ( body2It->getTypeID() == pe::Sphere::getStaticTypeID() )
+ {
+ pe::SphereID sphereJ = static_cast<pe::SphereID>( body2It.getBodyID() );
+ treatLubricationSphrSphr( sphereI, sphereJ, blockIt->getAABB() );
+ }
+ }
+ }
+ }
+
+ // lubrication correction for local bodies with global bodies (for example sphere-plane)
+ for( auto body1It = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID_ ); body1It != pe::LocalBodyIterator::end(); ++body1It )
+ {
+ if ( body1It->getTypeID() == pe::Sphere::getStaticTypeID() )
+ {
+ pe::SphereID sphereI = static_cast<pe::SphereID> ( body1It.getBodyID() );
+
+ for (auto body2It = globalBodyStorage_->begin(); body2It != globalBodyStorage_->end(); ++body2It)
+ {
+ if ( body2It->getTypeID() == pe::Plane::getStaticTypeID() )
+ {
+ // sphere-plane lubrication
+ pe::PlaneID planeJ = static_cast<pe::PlaneID>( body2It.getBodyID() );
+ treatLubricationSphrPlane( sphereI, planeJ );
+ } else if ( body2It->getTypeID() == pe::Sphere::getStaticTypeID() )
+ {
+ // sphere-sphere lubrication
+ pe::SphereID sphereJ = static_cast<pe::SphereID>( body2It.getBodyID() );
+ treatLubricationSphrSphr( sphereI, sphereJ, blockIt->getAABB() );
+ }
+ }
+ }
+ }
+ }
+}
+//*****************************************************************************************************************************************
+
+
+
+//////////////////////
+// Helper Functions //
+//////////////////////
+
+void LubricationCorrection::treatLubricationSphrSphr( const pe::SphereID sphereI, const pe::SphereID sphereJ, const math::AABB & blockAABB )
+{
+
+ WALBERLA_ASSERT_UNEQUAL( sphereI->getSystemID(), sphereJ->getSystemID() );
+
+ real_t gap = pe::getSurfaceDistance( sphereI, sphereJ );
+
+ if ( gap > cutOffDistance_ || gap < real_t(0) )
+ {
+ WALBERLA_LOG_DETAIL("gap " << gap << " larger than cutOff " << cutOffDistance_ << " - ignoring pair");
+ return;
+ }
+
+ if ( gap < minimalGapSize_ )
+ {
+ WALBERLA_LOG_DETAIL("gap " << gap << " smaller than minimal gap " << minimalGapSize_ << " - using minimal gap");
+ gap = minimalGapSize_;
+ }
+
+ const pe::Vec3 &posSphereI = sphereI->getPosition();
+ const pe::Vec3 &posSphereJ = sphereJ->getPosition();
+ pe::Vec3 fLub(0);
+
+ // compute (global) coordinate between spheres' centers of gravity
+ pe::Vec3 midPoint( (posSphereI + posSphereJ ) * real_c(0.5) );
+
+ // Let process on which midPoint lies do the lubrication correction
+ // or the local process of sphereI if sphereJ is global
+ if ( blockAABB.contains(midPoint) || sphereJ->isGlobal() )
+ {
+ fLub = compLubricationSphrSphr(gap, sphereI, sphereJ);
+ sphereI->addForce( fLub);
+ sphereJ->addForce(-fLub);
+
+ WALBERLA_LOG_DETAIL( "Lubrication force on sphere " << sphereI->getID() << " from sphere " << sphereJ->getID() << " is:" << fLub);
+ WALBERLA_LOG_DETAIL( "Lubrication force on sphere " << sphereJ->getID() << " from sphere " << sphereI->getID() << " is:" << -fLub << "\n");
+ }
+
+}
+
+void LubricationCorrection::treatLubricationSphrPlane( const pe::SphereID sphereI, const pe::ConstPlaneID planeJ )
+{
+
+ real_t gap = pe::getSurfaceDistance( sphereI, planeJ );
+
+ if ( gap > cutOffDistance_ || gap < real_t(0) )
+ {
+ WALBERLA_LOG_DETAIL("gap " << gap << " larger than cutOff " << cutOffDistance_ << " - ignoring pair");
+ return;
+ }
+
+ if ( gap < minimalGapSize_ )
+ {
+ WALBERLA_LOG_DETAIL("gap " << gap << " smaller than minimal gap " << minimalGapSize_ << " - using minimal gap");
+ gap = minimalGapSize_;
+ }
+
+ pe::Vec3 fLub = compLubricationSphrPlane( gap, sphereI, planeJ);
+
+ WALBERLA_LOG_DETAIL( "Lubrication force on sphere " << sphereI->getID() << " to plane with id " << planeJ->getID() << " is:" << fLub << std::endl );
+ sphereI->addForce( fLub );
+
+}
+
+
+
+//*****************************************************************************************************************************************
+/*! \brief Computes lubrication correction force between spheres.
+ * \ingroup pe_coupling
+ *
+ * Lubrication correction according to Ladd and Verberg, 2001
+ * ("Lattice-Boltzmann Simulations of Particle-Fluid Suspensions")
+ *
+ * Note: Verified quantitatively by computation in spreadsheet
+ * and qualitatively by considering direction of force for example setup.
+ */
+//*****************************************************************************************************************************************
+pe::Vec3 LubricationCorrection::compLubricationSphrSphr( real_t gap, const pe::SphereID sphereI, const pe::SphereID sphereJ) const
+{
+ const pe::Vec3 &posSphereI = sphereI->getPosition();
+ const pe::Vec3 &posSphereJ = sphereJ->getPosition();
+
+ const pe::Vec3 &tmpVec = posSphereJ - posSphereI;
+ const pe::Vec3 &rIJ = tmpVec.getNormalized();
+
+ real_t radiusSphereI = sphereI->getRadius();
+ real_t radiusSphereJ = sphereJ->getRadius();
+
+ pe::Vec3 velDiff(sphereI->getLinearVel() - sphereJ->getLinearVel());
+
+ real_t length = velDiff * rIJ;
+
+ real_t radiiSQR = ( radiusSphereI * radiusSphereJ ) * ( radiusSphereI * radiusSphereJ );
+ real_t radiiSumSQR = ( radiusSphereI + radiusSphereJ ) * ( radiusSphereI + radiusSphereJ );
+
+ pe::Vec3 fLub = ( -real_t(6) * dynamicViscosity_ * walberla::math::PI * radiiSQR / radiiSumSQR * ( real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
+
+ WALBERLA_LOG_DETAIL_SECTION()
+ {
+ std::stringstream ss;
+ ss << "Sphere I: \n uid:" << sphereI->getID() << "\n";
+ ss << "vel: " << sphereI->getLinearVel() << "\n";
+ ss << "rad: " << radiusSphereI << "\n";
+ ss << "pos: " << posSphereI << "\n\n";
+
+ ss << "Sphere J: \n uid:" << sphereJ->getID() << "\n";
+ ss << "vel: " << sphereJ->getLinearVel() << "\n";
+ ss << "rad: " << radiusSphereJ << "\n";
+ ss << "pos: " << posSphereJ << "\n\n";
+
+ real_t distance = gap + radiusSphereI + radiusSphereJ;
+ ss << "distance: " << distance << "\n";
+ ss << "viscosity: " << dynamicViscosity_ << "\n";
+
+ ss << "gap: " << gap << "\n";
+ ss << "cutOff: " << cutOffDistance_ << "\n";
+ ss << "velDiff " << velDiff << "\n";
+ ss << "rIJ " << rIJ << "\n\n";
+
+ ss << "Resulting lubrication force: " << fLub << "\n";
+
+ WALBERLA_LOG_DETAIL( ss.str() );
+ }
+
+ return fLub;
+}
+//*****************************************************************************************************************************************
+
+
+
+//*****************************************************************************************************************************************
+/*! \brief Computes lubrication correction force between sphere and wall.
+ * \ingroup pe_coupling
+ *
+ * Lubrication correction according to Ladd and Verberg, 2001
+ * ("Lattice-Boltzmann Simulations of Particle-Fluid Suspensions" )
+ *
+ * Note: Verified quantitatively by computation in spreadsheet
+ * and qualitatively by considering direction of force for example setup.
+ */
+//*****************************************************************************************************************************************
+pe::Vec3 LubricationCorrection::compLubricationSphrPlane( real_t gap, const pe::SphereID sphereI, const pe::ConstPlaneID planeJ) const
+{
+ const pe::Vec3 &posSphereI( sphereI->getPosition() );
+ real_t radiusSphereI = sphereI->getRadius();
+
+ const pe::Vec3 &planeNormal( planeJ->getNormal() ); // took negative of normal from sphere to plane (plane's normal) - sign cancels out anyway
+ pe::Vec3 rIJ( planeNormal.getNormalized() ); // for safety reasons, normalize normal
+
+ real_t length = sphereI->getLinearVel() * rIJ;
+
+ real_t radiiSQR = radiusSphereI * radiusSphereI;
+
+ pe::Vec3 fLub( -real_t(6) * dynamicViscosity_ * walberla::math::PI * radiiSQR * (real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
+
+ WALBERLA_LOG_DETAIL_SECTION() {
+ std::stringstream ss;
+ ss << "Sphere I: \n uid:" << sphereI->getID() << "\n";
+ ss << "vel: " << sphereI->getLinearVel() << "\n";
+ ss << "rad: " << radiusSphereI << "\n";
+ ss << "pos: " << posSphereI << "\n\n";
+
+ real_t distance = gap + radiusSphereI;
+ ss << "distance: " << distance << "\n";
+ ss << "viscosity: " << dynamicViscosity_ << "\n";
+
+ ss << "gap: " << gap << "\n";
+ ss << "cutOff: " << cutOffDistance_ << "\n";
+ ss << "velDiff " << sphereI->getLinearVel() << "\n";
+ ss << "rIJ " << -rIJ << "\n\n";
+
+ ss << "Resulting lubrication force: " << fLub << "\n";
+
+ WALBERLA_LOG_DETAIL( ss.str() );
+ }
+
+ return fLub;
+}
+
+
+} // pe_coupling
+} // walberla
diff --git a/src/pe_coupling/utility/LubricationCorrection.h b/src/pe_coupling/utility/LubricationCorrection.h
index 4a5afa830d09352ff53ae87c69e88c96c7533eae..b4f34181505f1a42472222875c72f74596f62162 100644
--- a/src/pe_coupling/utility/LubricationCorrection.h
+++ b/src/pe_coupling/utility/LubricationCorrection.h
@@ -24,15 +24,11 @@
#pragma once
-#include "core/logging/all.h"
-#include "core/debug/Debug.h"
-#include "domain_decomposition/BlockStorage.h"
-#include "lbm/field/PdfField.h"
-#include "pe/rigidbody/BodyIterators.h"
+#include "domain_decomposition/StructuredBlockStorage.h"
+
#include "pe/rigidbody/Plane.h"
#include "pe/rigidbody/Sphere.h"
-#include "pe/utility/Distance.h"
namespace walberla {
namespace pe_coupling {
@@ -67,7 +63,6 @@ private:
// member variables
shared_ptr<StructuredBlockStorage> blockStorage_;
shared_ptr<pe::BodyStorage> globalBodyStorage_;
- const BlockDataID pdfFieldID_;
const BlockDataID bodyStorageID_;
real_t dynamicViscosity_;
@@ -76,247 +71,5 @@ private:
}; // class LubricationCorrection
-
-void LubricationCorrection::operator ()()
-{
- WALBERLA_LOG_PROGRESS( "Calculating Lubrication Force" );
-
- for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt)
- {
- // loop over all rigid bodies
- for( auto body1It = pe::BodyIterator::begin( *blockIt, bodyStorageID_ ); body1It != pe::BodyIterator::end(); ++body1It )
- {
- // lubrication forces for spheres
- if ( body1It->getTypeID() == pe::Sphere::getStaticTypeID() )
- {
- pe::SphereID sphereI = static_cast<pe::SphereID> ( body1It.getBodyID() );
-
- auto copyBody1It = body1It;
- // loop over all rigid bodies after current body1 to avoid double forces
- for( auto body2It = (++copyBody1It); body2It != pe::BodyIterator::end(); ++body2It )
- {
- // sphere-sphere lubrication
- if ( body2It->getTypeID() == pe::Sphere::getStaticTypeID() )
- {
- pe::SphereID sphereJ = static_cast<pe::SphereID>( body2It.getBodyID() );
- treatLubricationSphrSphr( sphereI, sphereJ, blockIt->getAABB() );
- }
- }
- }
- }
-
- // lubrication correction for local bodies with global bodies (for example sphere-plane)
- for( auto body1It = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID_ ); body1It != pe::LocalBodyIterator::end(); ++body1It )
- {
- if ( body1It->getTypeID() == pe::Sphere::getStaticTypeID() )
- {
- pe::SphereID sphereI = static_cast<pe::SphereID> ( body1It.getBodyID() );
-
- for (auto body2It = globalBodyStorage_->begin(); body2It != globalBodyStorage_->end(); ++body2It)
- {
- if ( body2It->getTypeID() == pe::Plane::getStaticTypeID() )
- {
- // sphere-plane lubrication
- pe::PlaneID planeJ = static_cast<pe::PlaneID>( body2It.getBodyID() );
- treatLubricationSphrPlane( sphereI, planeJ );
- } else if ( body2It->getTypeID() == pe::Sphere::getStaticTypeID() )
- {
- // sphere-sphere lubrication
- pe::SphereID sphereJ = static_cast<pe::SphereID>( body2It.getBodyID() );
- treatLubricationSphrSphr( sphereI, sphereJ, blockIt->getAABB() );
- }
- }
- }
- }
- }
-}
-//*****************************************************************************************************************************************
-
-
-
-//////////////////////
-// Helper Functions //
-//////////////////////
-
-void LubricationCorrection::treatLubricationSphrSphr( const pe::SphereID sphereI, const pe::SphereID sphereJ, const math::AABB & blockAABB )
-{
-
- WALBERLA_ASSERT_UNEQUAL( sphereI->getSystemID(), sphereJ->getSystemID() );
-
- real_t gap = pe::getSurfaceDistance( sphereI, sphereJ );
-
- if ( gap > cutOffDistance_ || gap < real_t(0) )
- {
- WALBERLA_LOG_DETAIL("gap " << gap << " larger than cutOff " << cutOffDistance_ << " - ignoring pair");
- return;
- }
-
- if ( gap < minimalGapSize_ )
- {
- WALBERLA_LOG_DETAIL("gap " << gap << " smaller than minimal gap " << minimalGapSize_ << " - using minimal gap");
- gap = minimalGapSize_;
- }
-
- const pe::Vec3 &posSphereI = sphereI->getPosition();
- const pe::Vec3 &posSphereJ = sphereJ->getPosition();
- pe::Vec3 fLub(0);
-
- // compute (global) coordinate between spheres' centers of gravity
- pe::Vec3 midPoint( (posSphereI + posSphereJ ) * real_c(0.5) );
-
- // Let process on which midPoint lies do the lubrication correction
- // or the local process of sphereI if sphereJ is global
- if ( blockAABB.contains(midPoint) || sphereJ->isGlobal() )
- {
- fLub = compLubricationSphrSphr(gap, sphereI, sphereJ);
- sphereI->addForce( fLub);
- sphereJ->addForce(-fLub);
-
- WALBERLA_LOG_DETAIL( "Lubrication force on sphere " << sphereI->getID() << " from sphere " << sphereJ->getID() << " is:" << fLub);
- WALBERLA_LOG_DETAIL( "Lubrication force on sphere " << sphereJ->getID() << " from sphere " << sphereI->getID() << " is:" << -fLub << "\n");
- }
-
-}
-
-void LubricationCorrection::treatLubricationSphrPlane( const pe::SphereID sphereI, const pe::ConstPlaneID planeJ )
-{
-
- real_t gap = pe::getSurfaceDistance( sphereI, planeJ );
-
- if ( gap > cutOffDistance_ || gap < real_t(0) )
- {
- WALBERLA_LOG_DETAIL("gap " << gap << " larger than cutOff " << cutOffDistance_ << " - ignoring pair");
- return;
- }
-
- if ( gap < minimalGapSize_ )
- {
- WALBERLA_LOG_DETAIL("gap " << gap << " smaller than minimal gap " << minimalGapSize_ << " - using minimal gap");
- gap = minimalGapSize_;
- }
-
- pe::Vec3 fLub = compLubricationSphrPlane( gap, sphereI, planeJ);
-
- WALBERLA_LOG_DETAIL( "Lubrication force on sphere " << sphereI->getID() << " to plane with id " << planeJ->getID() << " is:" << fLub << std::endl );
- sphereI->addForce( fLub );
-
-}
-
-
-
-//*****************************************************************************************************************************************
-/*! \brief Computes lubrication correction force between spheres.
- * \ingroup pe_coupling
- *
- * Lubrication correction according to Ladd and Verberg, 2001
- * ("Lattice-Boltzmann Simulations of Particle-Fluid Suspensions")
- *
- * Note: Verified quantitatively by computation in spreadsheet
- * and qualitatively by considering direction of force for example setup.
- */
-//*****************************************************************************************************************************************
-pe::Vec3 LubricationCorrection::compLubricationSphrSphr( real_t gap, const pe::SphereID sphereI, const pe::SphereID sphereJ) const
-{
- const pe::Vec3 &posSphereI = sphereI->getPosition();
- const pe::Vec3 &posSphereJ = sphereJ->getPosition();
-
- const pe::Vec3 &tmpVec = posSphereJ - posSphereI;
- const pe::Vec3 &rIJ = tmpVec.getNormalized();
-
- real_t radiusSphereI = sphereI->getRadius();
- real_t radiusSphereJ = sphereJ->getRadius();
-
- pe::Vec3 velDiff(sphereI->getLinearVel() - sphereJ->getLinearVel());
-
- real_t length = velDiff * rIJ;
-
- real_t radiiSQR = ( radiusSphereI * radiusSphereJ ) * ( radiusSphereI * radiusSphereJ );
- real_t radiiSumSQR = ( radiusSphereI + radiusSphereJ ) * ( radiusSphereI + radiusSphereJ );
-
- pe::Vec3 fLub = ( -real_t(6) * dynamicViscosity_ * walberla::math::PI * radiiSQR / radiiSumSQR * ( real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
-
- WALBERLA_LOG_DETAIL_SECTION()
- {
- std::stringstream ss;
- ss << "Sphere I: \n uid:" << sphereI->getID() << "\n";
- ss << "vel: " << sphereI->getLinearVel() << "\n";
- ss << "rad: " << radiusSphereI << "\n";
- ss << "pos: " << posSphereI << "\n\n";
-
- ss << "Sphere J: \n uid:" << sphereJ->getID() << "\n";
- ss << "vel: " << sphereJ->getLinearVel() << "\n";
- ss << "rad: " << radiusSphereJ << "\n";
- ss << "pos: " << posSphereJ << "\n\n";
-
- real_t distance = gap + radiusSphereI + radiusSphereJ;
- ss << "distance: " << distance << "\n";
- ss << "viscosity: " << dynamicViscosity_ << "\n";
-
- ss << "gap: " << gap << "\n";
- ss << "cutOff: " << cutOffDistance_ << "\n";
- ss << "velDiff " << velDiff << "\n";
- ss << "rIJ " << rIJ << "\n\n";
-
- ss << "Resulting lubrication force: " << fLub << "\n";
-
- WALBERLA_LOG_DETAIL( ss.str() );
- }
-
- return fLub;
-}
-//*****************************************************************************************************************************************
-
-
-
-//*****************************************************************************************************************************************
-/*! \brief Computes lubrication correction force between sphere and wall.
- * \ingroup pe_coupling
- *
- * Lubrication correction according to Ladd and Verberg, 2001
- * ("Lattice-Boltzmann Simulations of Particle-Fluid Suspensions" )
- *
- * Note: Verified quantitatively by computation in spreadsheet
- * and qualitatively by considering direction of force for example setup.
- */
-//*****************************************************************************************************************************************
-pe::Vec3 LubricationCorrection::compLubricationSphrPlane( real_t gap, const pe::SphereID sphereI, const pe::ConstPlaneID planeJ) const
-{
- const pe::Vec3 &posSphereI( sphereI->getPosition() );
- real_t radiusSphereI = sphereI->getRadius();
-
- const pe::Vec3 &planeNormal( planeJ->getNormal() ); // took negative of normal from sphere to plane (plane's normal) - sign cancels out anyway
- pe::Vec3 rIJ( planeNormal.getNormalized() ); // for safety reasons, normalize normal
-
- real_t length = sphereI->getLinearVel() * rIJ;
-
- real_t radiiSQR = radiusSphereI * radiusSphereI;
-
- pe::Vec3 fLub( -real_t(6) * dynamicViscosity_ * walberla::math::PI * radiiSQR * (real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
-
- WALBERLA_LOG_DETAIL_SECTION() {
- std::stringstream ss;
- ss << "Sphere I: \n uid:" << sphereI->getID() << "\n";
- ss << "vel: " << sphereI->getLinearVel() << "\n";
- ss << "rad: " << radiusSphereI << "\n";
- ss << "pos: " << posSphereI << "\n\n";
-
- real_t distance = gap + radiusSphereI;
- ss << "distance: " << distance << "\n";
- ss << "viscosity: " << dynamicViscosity_ << "\n";
-
- ss << "gap: " << gap << "\n";
- ss << "cutOff: " << cutOffDistance_ << "\n";
- ss << "velDiff " << sphereI->getLinearVel() << "\n";
- ss << "rIJ " << -rIJ << "\n\n";
-
- ss << "Resulting lubrication force: " << fLub << "\n";
-
- WALBERLA_LOG_DETAIL( ss.str() );
- }
-
- return fLub;
-}
-
-
} // pe_coupling
} // walberla
diff --git a/src/pe_coupling/utility/TimeStep.cpp b/src/pe_coupling/utility/TimeStep.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..79a9d85d03946047bb91c2bcfda7d40c7a44aaea
--- /dev/null
+++ b/src/pe_coupling/utility/TimeStep.cpp
@@ -0,0 +1,112 @@
+//======================================================================================================================
+//
+// 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 TimeStep.cpp
+//! \ingroup pe_coupling
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "TimeStep.h"
+
+#include "pe/rigidbody/BodyIterators.h"
+
+#include <map>
+#include <array>
+
+namespace walberla {
+namespace pe_coupling {
+
+void TimeStep::operator()()
+{
+ if( numberOfSubIterations_ == 1 )
+ {
+ forceEvaluationFunc_();
+
+ collisionResponse_.timestep( timeStepSize_ );
+ synchronizeFunc_();
+ }
+ else
+ {
+ // during the intermediate time steps of the collision response, the currently acting forces
+ // (interaction forces, gravitational force, ...) have to remain constant.
+ // Since they are reset after the call to collision response, they have to be stored explicitly before.
+ // Then they are set again after each intermediate step.
+
+ // generate map from all known bodies (process local) to total forces/torques
+ // this has to be done on a block-local basis, since the same body could reside on several blocks from this process
+ typedef domain_decomposition::IBlockID::IDType BlockID_T;
+ std::map< BlockID_T, std::map< walberla::id_t, std::array< real_t, 6 > > > forceTorqueMap;
+
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ BlockID_T blockID = blockIt->getId().getID();
+ auto& blockLocalForceTorqueMap = forceTorqueMap[blockID];
+
+ // iterate over local and remote bodies and store force/torque in map
+ for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ auto & f = blockLocalForceTorqueMap[ bodyIt->getSystemID() ];
+
+ const auto & force = bodyIt->getForce();
+ const auto & torque = bodyIt->getTorque();
+
+ f = {{force[0], force[1], force[2], torque[0], torque[1], torque[2] }};
+ }
+ }
+
+ // perform pe time steps
+ const real_t subTimeStepSize = timeStepSize_ / real_c( numberOfSubIterations_ );
+ for( uint_t i = 0; i != numberOfSubIterations_; ++i )
+ {
+
+ // in the first iteration, forces are already set
+ if( i != 0 )
+ {
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ BlockID_T blockID = blockIt->getId().getID();
+ auto& blockLocalForceTorqueMap = forceTorqueMap[blockID];
+
+ // re-set stored force/torque on bodies
+ for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+
+ const auto f = blockLocalForceTorqueMap.find( bodyIt->getSystemID() );
+
+ if( f != blockLocalForceTorqueMap.end() )
+ {
+ const auto & ftValues = f->second;
+ bodyIt->addForce ( ftValues[0], ftValues[1], ftValues[2] );
+ bodyIt->addTorque( ftValues[3], ftValues[4], ftValues[5] );
+ }
+ }
+ }
+ }
+
+ // evaluate forces (e.g. lubrication forces)
+ forceEvaluationFunc_();
+
+ collisionResponse_.timestep( subTimeStepSize );
+ synchronizeFunc_();
+ }
+ }
+}
+
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/utility/TimeStep.h b/src/pe_coupling/utility/TimeStep.h
index f32c306142dbe01c904f7e081e8aa49dd2cd3455..283282b6893e9b6307b7dcb77b2ddc734dee9523 100644
--- a/src/pe_coupling/utility/TimeStep.h
+++ b/src/pe_coupling/utility/TimeStep.h
@@ -23,20 +23,12 @@
#pragma once
-#include "core/debug/Debug.h"
#include "core/timing/TimingTree.h"
-
-#include "domain_decomposition/BlockStorage.h"
-
+#include "domain_decomposition/StructuredBlockStorage.h"
#include "pe/cr/ICR.h"
-#include "pe/rigidbody/BodyIterators.h"
-#include "pe/synchronization/SyncForces.h"
#include <functional>
-#include <map>
-#include <array>
-
namespace walberla {
namespace pe_coupling {
@@ -72,83 +64,9 @@ public:
, forceEvaluationFunc_( forceEvaluationFunc )
{}
- void operator()()
- {
- if( numberOfSubIterations_ == 1 )
- {
- forceEvaluationFunc_();
-
- collisionResponse_.timestep( timeStepSize_ );
- synchronizeFunc_();
- }
- else
- {
- // during the intermediate time steps of the collision response, the currently acting forces
- // (interaction forces, gravitational force, ...) have to remain constant.
- // Since they are reset after the call to collision response, they have to be stored explicitly before.
- // Then they are set again after each intermediate step.
-
- // generate map from all known bodies (process local) to total forces/torques
- // this has to be done on a block-local basis, since the same body could reside on several blocks from this process
- typedef domain_decomposition::IBlockID::IDType BlockID_T;
- std::map< BlockID_T, std::map< walberla::id_t, std::array< real_t, 6 > > > forceTorqueMap;
-
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- BlockID_T blockID = blockIt->getId().getID();
- auto& blockLocalForceTorqueMap = forceTorqueMap[blockID];
-
- // iterate over local and remote bodies and store force/torque in map
- for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- auto & f = blockLocalForceTorqueMap[ bodyIt->getSystemID() ];
-
- const auto & force = bodyIt->getForce();
- const auto & torque = bodyIt->getTorque();
-
- f = {{force[0], force[1], force[2], torque[0], torque[1], torque[2] }};
- }
- }
-
- // perform pe time steps
- const real_t subTimeStepSize = timeStepSize_ / real_c( numberOfSubIterations_ );
- for( uint_t i = 0; i != numberOfSubIterations_; ++i )
- {
-
- // in the first iteration, forces are already set
- if( i != 0 )
- {
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- BlockID_T blockID = blockIt->getId().getID();
- auto& blockLocalForceTorqueMap = forceTorqueMap[blockID];
-
- // re-set stored force/torque on bodies
- for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
-
- const auto f = blockLocalForceTorqueMap.find( bodyIt->getSystemID() );
-
- if( f != blockLocalForceTorqueMap.end() )
- {
- const auto & ftValues = f->second;
- bodyIt->addForce ( ftValues[0], ftValues[1], ftValues[2] );
- bodyIt->addTorque( ftValues[3], ftValues[4], ftValues[5] );
- }
- }
- }
- }
-
- // evaluate forces (e.g. lubrication forces)
- forceEvaluationFunc_();
-
- collisionResponse_.timestep( subTimeStepSize );
- synchronizeFunc_();
- }
- }
- }
+ void operator()();
-protected:
+private:
const real_t timeStepSize_;
const uint_t numberOfSubIterations_;