diff --git a/src/pe_coupling/utility/BodiesForceTorqueContainer.h b/src/pe_coupling/utility/BodiesForceTorqueContainer.h
index a32dc589d99d0a1138b4e5d1cab974338dc0933e..d69b6791b23734526521e10a3b42a6acb8b5f6bf 100644
--- a/src/pe_coupling/utility/BodiesForceTorqueContainer.h
+++ b/src/pe_coupling/utility/BodiesForceTorqueContainer.h
@@ -21,13 +21,13 @@
#pragma once
+#include "blockforest/StructuredBlockForest.h"
#include "core/math/Vector3.h"
-#include "domain_decomposition/StructuredBlockStorage.h"
#include "pe/rigidbody/BodyIterators.h"
#include "pe/synchronization/SyncForces.h"
#include <map>
-#include <vector>
+#include <array>
namespace walberla {
namespace pe_coupling {
@@ -36,9 +36,14 @@ class BodiesForceTorqueContainer
{
public:
- BodiesForceTorqueContainer( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & bodyStorageID )
- : blockStorage_( blockStorage ), bodyStorageID_( bodyStorageID )
- { }
+ typedef std::map< walberla::id_t, std::array<real_t,6> > ForceTorqueStorage_T;
+
+ BodiesForceTorqueContainer( const shared_ptr<StructuredBlockForest> & blockForest, const BlockDataID & bodyStorageID )
+ : blockForest_( blockForest ), bodyStorageID_( bodyStorageID )
+ {
+ // has to be added to the forest (not the storage) to register correctly
+ bodyForceTorqueStorageID_ = blockForest->addBlockData(make_shared<blockforest::AlwaysCreateBlockDataHandling<ForceTorqueStorage_T> >(), "BodiesForceTorqueContainer");
+ }
void operator()()
{
@@ -50,38 +55,18 @@ public:
// clear map
clear();
- // sum up all forces/torques from shadow copies on local body (owner)
- pe::reduceForces( blockStorage_->getBlockStorage(), bodyStorageID_ );
- // send total forces/torques to shadow owners
- pe::distributeForces( blockStorage_->getBlockStorage(), bodyStorageID_ );
-
// (re-)build map
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ 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 = bodyForceTorqueMap_[ bodyIt->getSystemID() ];
+ auto & f = (*bodyForceTorqueStorage)[ bodyIt->getSystemID() ];
- // only add if body has not been added already before (from another block)
- if( f.empty() )
- {
- const auto & force = bodyIt->getForce();
- f.push_back( force[0] );
- f.push_back( force[1] );
- f.push_back( force[2] );
-
- const auto & torque = bodyIt->getTorque();
- f.push_back( torque[0] );
- f.push_back( torque[1] );
- f.push_back( torque[2] );
- }
-
- // reset of force/torque on remote bodies necessary to erase the multiple occurrences of forces/torques on bodies
- // (due to call to distributeForces() before)
- if ( bodyIt->isRemote() )
- {
- bodyIt->resetForceAndTorque();
- }
+ const auto & force = bodyIt->getForce();
+ const auto & torque = bodyIt->getTorque();
+ f = {{force[0], force[1], force[2], torque[0], torque[1], torque[2] }};
}
}
@@ -89,42 +74,52 @@ public:
void setOnBodies()
{
- // owning process sets the force/torque on the bodies
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ // set the force/torque stored in the block-local map onto all bodies
+ for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
{
- for( auto bodyIt = pe::LocalBodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+ auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
+
+ for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
- const auto &f = bodyForceTorqueMap_[bodyIt->getSystemID()];
- WALBERLA_ASSERT( !f.empty(), "When attempting to set force/torque on local body " << bodyIt->getSystemID() << " at position " << bodyIt->getPosition() << ", body was not found in map!");
- bodyIt->addForce ( f[0], f[1], f[2] );
- bodyIt->addTorque( f[3], f[4], f[5] );
+ 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 clear()
{
- bodyForceTorqueMap_.clear();
+ for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
+ {
+ auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
+ bodyForceTorqueStorage->clear();
+ }
}
void swap( BodiesForceTorqueContainer & other )
{
- std::swap( this->bodyForceTorqueMap_, other.bodyForceTorqueMap_ );
+ std::swap( bodyForceTorqueStorageID_, other.bodyForceTorqueStorageID_);
}
private:
- shared_ptr<StructuredBlockStorage> blockStorage_;
+ shared_ptr<StructuredBlockStorage> blockForest_;
const BlockDataID bodyStorageID_;
- std::map< walberla::id_t, std::vector<real_t> > bodyForceTorqueMap_;
+ BlockDataID bodyForceTorqueStorageID_;
};
class BodyContainerSwapper
{
public:
- BodyContainerSwapper( const shared_ptr<BodiesForceTorqueContainer> & cont1,
- const shared_ptr<BodiesForceTorqueContainer> & cont2 )
+ BodyContainerSwapper( const shared_ptr<BodiesForceTorqueContainer> & cont1, const shared_ptr<BodiesForceTorqueContainer> & cont2 )
: cont1_( cont1 ), cont2_( cont2 )
{ }
diff --git a/src/pe_coupling/utility/BodySelectorFunctions.cpp b/src/pe_coupling/utility/BodySelectorFunctions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4c10f95028c47e88b3340aef74840337dc723992
--- /dev/null
+++ b/src/pe_coupling/utility/BodySelectorFunctions.cpp
@@ -0,0 +1,50 @@
+//======================================================================================================================
+//
+// 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 BodySelectorFunctions.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "BodySelectorFunctions.h"
+
+#include "pe/rigidbody/RigidBody.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+bool selectAllBodies(pe::BodyID /*bodyID*/)
+{
+ return true;
+}
+
+bool selectRegularBodies(pe::BodyID bodyID)
+{
+ return !bodyID->hasInfiniteMass() && !bodyID->isGlobal();
+}
+
+bool selectFixedBodies(pe::BodyID bodyID)
+{
+ return bodyID->hasInfiniteMass() && !bodyID->isGlobal();
+}
+
+bool selectGlobalBodies(pe::BodyID bodyID)
+{
+ return bodyID->isGlobal();
+}
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/utility/BodySelectorFunctions.h b/src/pe_coupling/utility/BodySelectorFunctions.h
index 75bc5c93b1e12c5c970c8b0f11556f08e5980d06..a811d416b2fe440d2ff4a0c062b18b7f0823eae2 100644
--- a/src/pe_coupling/utility/BodySelectorFunctions.h
+++ b/src/pe_coupling/utility/BodySelectorFunctions.h
@@ -21,29 +21,18 @@
#pragma once
-#include "pe/Types.h"
+#include "pe/rigidbody/RigidBody.h"
namespace walberla {
namespace pe_coupling {
-bool selectAllBodies(pe::BodyID /*bodyID*/)
-{
- return true;
-}
-
-bool selectRegularBodies(pe::BodyID bodyID)
-{
- return !bodyID->hasInfiniteMass() && !bodyID->isGlobal();
-}
-
-bool selectFixedBodies(pe::BodyID bodyID)
-{
- return bodyID->hasInfiniteMass() && !bodyID->isGlobal();
-}
-
-bool selectGlobalBodies(pe::BodyID bodyID)
-{
- return bodyID->isGlobal();
-}
+bool selectAllBodies(pe::BodyID /*bodyID*/);
+
+bool selectRegularBodies(pe::BodyID bodyID);
+
+bool selectFixedBodies(pe::BodyID bodyID);
+
+bool selectGlobalBodies(pe::BodyID bodyID);
+
} // namespace pe_coupling
} // namespace walberla
diff --git a/src/pe_coupling/utility/TimeStep.h b/src/pe_coupling/utility/TimeStep.h
index d89e0a6e7aff337ba570c252e605fe4d006c8e9c..f32c306142dbe01c904f7e081e8aa49dd2cd3455 100644
--- a/src/pe_coupling/utility/TimeStep.h
+++ b/src/pe_coupling/utility/TimeStep.h
@@ -35,7 +35,7 @@
#include <functional>
#include <map>
-
+#include <array>
namespace walberla {
namespace pe_coupling {
@@ -88,43 +88,25 @@ public:
// 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.
- // sum up all forces/torques from shadow copies on local body (owner)
- pe::reduceForces( blockStorage_->getBlockStorage(), bodyStorageID_ );
- // send total forces/torques to shadow owners
- pe::distributeForces( blockStorage_->getBlockStorage(), bodyStorageID_ );
-
// generate map from all known bodies (process local) to total forces/torques
- std::map< walberla::id_t, std::vector< real_t > > forceMap;
+ // 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
- // Remote bodies are required since during the then following collision response time steps,
- // bodies can change ownership (i.e. a now remote body could become a local body ).
- // Since only the owning process re-sets the force/torque later, remote body values have to be stored as well.
for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
{
- auto & f = forceMap[ bodyIt->getSystemID() ];
+ auto & f = blockLocalForceTorqueMap[ bodyIt->getSystemID() ];
- // only add if body has not been added already before (from another block)
- if( f.empty() )
- {
- const auto & force = bodyIt->getForce();
- f.push_back( force[0] );
- f.push_back( force[1] );
- f.push_back( force[2] );
-
- const auto & torque = bodyIt->getTorque();
- f.push_back( torque[0] );
- f.push_back( torque[1] );
- f.push_back( torque[2] );
- }
+ const auto & force = bodyIt->getForce();
+ const auto & torque = bodyIt->getTorque();
- // reset of force/torque on remote bodies necessary to erase the multiple occurrences of forces/torques on bodies
- // (due to call to distributeForces() before)
- if ( bodyIt->isRemote() )
- {
- bodyIt->resetForceAndTorque();
- }
+ f = {{force[0], force[1], force[2], torque[0], torque[1], torque[2] }};
}
}
@@ -133,18 +115,26 @@ public:
for( uint_t i = 0; i != numberOfSubIterations_; ++i )
{
- // in the first iteration, forces on local bodies are already set by force synchronization
+ // in the first iteration, forces are already set
if( i != 0 )
{
for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
{
- // only owning process sets force/torque on bodies
- for( auto bodyIt = pe::LocalBodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+ 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 = forceMap[ bodyIt->getSystemID() ];
- WALBERLA_ASSERT( !f.empty(), "When attempting to set force/torque on local body " << bodyIt->getSystemID() << " at position " << bodyIt->getPosition() << ", body was not found in map!");
- bodyIt->addForce ( f[0], f[1], f[2] );
- bodyIt->addTorque( f[3], f[4], f[5] );
+
+ 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] );
+ }
}
}
}
diff --git a/tests/pe_coupling/CMakeLists.txt b/tests/pe_coupling/CMakeLists.txt
index 1803b3090ccc1f744a4b1711a81ac3b8732d2828..b53b734b5b16b8c18841938c4839e88e0dd7885f 100644
--- a/tests/pe_coupling/CMakeLists.txt
+++ b/tests/pe_coupling/CMakeLists.txt
@@ -156,3 +156,15 @@ waLBerla_execute_test( NAME SphereWallCollisionBehaviorDPMFuncTest COMMAND $<TAR
waLBerla_compile_test( FILES discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp DEPENDS blockforest pe timeloop )
waLBerla_execute_test( NAME HinderedSettlingDynamicsDPMFuncTest COMMAND $<TARGET_FILE:HinderedSettlingDynamicsDPM> --funcTest PROCESSES 4 LABELS longrun CONFIGURATIONS RelWithDbgInfo )
+
+###################################################################################################
+# Utility tests
+###################################################################################################
+
+waLBerla_compile_test( FILES utility/BodiesForceTorqueContainerTest.cpp DEPENDS blockforest pe timeloop )
+waLBerla_execute_test( NAME BodiesForceTorqueContainerTest COMMAND $<TARGET_FILE:BodiesForceTorqueContainerTest> PROCESSES 1 )
+waLBerla_execute_test( NAME BodiesForceTorqueContainerParallelTest COMMAND $<TARGET_FILE:BodiesForceTorqueContainerTest> PROCESSES 3 )
+
+waLBerla_compile_test( FILES utility/PeSubCyclingTest.cpp DEPENDS blockforest pe timeloop )
+waLBerla_execute_test( NAME PeSubCyclingTest COMMAND $<TARGET_FILE:PeSubCyclingTest> PROCESSES 1 )
+waLBerla_execute_test( NAME PeSubCyclingParallelTest COMMAND $<TARGET_FILE:PeSubCyclingTest> PROCESSES 3 )
\ No newline at end of file
diff --git a/tests/pe_coupling/utility/BodiesForceTorqueContainerTest.cpp b/tests/pe_coupling/utility/BodiesForceTorqueContainerTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ed015b72f6c78f0ecf445ea4e7e8295a036f1d1e
--- /dev/null
+++ b/tests/pe_coupling/utility/BodiesForceTorqueContainerTest.cpp
@@ -0,0 +1,487 @@
+//======================================================================================================================
+//
+// 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 BodiesForceTorqueContainerTest.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/math/all.h"
+
+#include "pe/basic.h"
+#include "pe/utility/DestroyBody.h"
+
+#include <pe_coupling/utility/all.h>
+
+namespace force_torque_container_test
+{
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+
+typedef boost::tuple<pe::Sphere> BodyTypeTuple ;
+
+
+/*!\brief Test cases for the force torque container provided by the coupling module
+ *
+ * Spheres at different positions are created and force(s) and torque(s) are applied.
+ * Then, they are store in the container and reset on the body.
+ * Then, in some cases, the sphere is moved to cross block borders, and the stored forces/torques are reapplied by the container again.
+ * The obtained force/torque values are compared against the originally set (and thus expected) values.
+ */
+//////////
+// MAIN //
+//////////
+int main( int argc, char **argv )
+{
+ debug::enterTestMode();
+
+ mpi::Environment env( argc, argv );
+
+ // uncomment to have logging
+ //logging::Logging::instance()->setLogLevel(logging::Logging::LogLevel::DETAIL);
+
+ const real_t dx = real_t(1);
+ const real_t radius = real_t(5);
+
+ ///////////////////////////
+ // DATA STRUCTURES SETUP //
+ ///////////////////////////
+
+ Vector3<uint_t> blocksPerDirection(uint_t(3), uint_t(1), uint_t(1));
+ Vector3<uint_t> cellsPerBlock(uint_t(20), uint_t(20), uint_t(20));
+ Vector3<bool> periodicity(true, false, false);
+
+ // create fully periodic domain with refined blocks
+ auto blocks = blockforest::createUniformBlockGrid( blocksPerDirection[0], blocksPerDirection[1], blocksPerDirection[2],
+ cellsPerBlock[0], cellsPerBlock[1], cellsPerBlock[2],
+ dx,
+ 0, false, false,
+ periodicity[0], periodicity[1], periodicity[2],
+ false );
+
+
+ // pe body storage
+ pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
+ shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
+ auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "Storage");
+ auto sphereMaterialID = pe::createMaterial( "sphereMat", real_t(1) , real_t(0.3), real_t(0.2), real_t(0.2), real_t(0.24), real_t(200), real_t(200), real_t(0), real_t(0) );
+
+ // pe coupling
+ const real_t overlap = real_t( 1.5 ) * dx;
+ std::function<void(void)> syncCall = std::bind( pe::syncNextNeighbors<BodyTypeTuple>, boost::ref(blocks->getBlockForest()), bodyStorageID, static_cast<WcTimingTree*>(NULL), overlap, false );
+
+ // sphere positions for test scenarios
+ Vector3<real_t> positionInsideBlock(real_t(10), real_t(10), real_t(10));
+ Vector3<real_t> positionAtBlockBorder(real_t(19.5), real_t(10), real_t(10));
+ Vector3<real_t> positionAtBlockBorderUpdated(real_t(20.5), real_t(10), real_t(10));
+
+ Vector3<real_t> positionAtBlockBorder2(real_t(20) + radius + overlap - real_t(0.5), real_t(10), real_t(10));
+ Vector3<real_t> positionAtBlockBorderUpdated2(real_t(20) + radius + overlap + real_t(0.5), real_t(10), real_t(10));
+
+
+ Vector3<real_t> testForce(real_t(2), real_t(1), real_t(0));
+ Vector3<real_t> torqueOffset = Vector3<real_t>(real_t(1), real_t(0), real_t(0));
+
+ pe_coupling::ForceTorqueOnBodiesResetter resetter(blocks, bodyStorageID);
+ shared_ptr<pe_coupling::BodiesForceTorqueContainer> container1 = make_shared<pe_coupling::BodiesForceTorqueContainer>(blocks, bodyStorageID);
+ shared_ptr<pe_coupling::BodiesForceTorqueContainer> container2 = make_shared<pe_coupling::BodiesForceTorqueContainer>(blocks, bodyStorageID);
+ pe_coupling::BodyContainerSwapper containerSwapper(container1, container2);
+
+ //////////////////
+ // Inside block //
+ //////////////////
+ {
+ std::string testIdentifier("Test: sphere inside block");
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+ pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+ positionInsideBlock, radius, sphereMaterialID, false, true, false);
+
+ syncCall();
+
+ uint_t applicationCounter( 0 );
+
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ ++applicationCounter;
+ auto pos = bodyIt->getPosition();
+ bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+ }
+ }
+ mpi::allReduceInplace(applicationCounter, mpi::SUM);
+
+ container1->store();
+
+ resetter();
+
+ containerSwapper();
+
+ container2->setOnBodies();
+
+ Vector3<real_t> expectedForce = applicationCounter * testForce;
+ Vector3<real_t> expectedTorque = applicationCounter * ( torqueOffset % testForce );
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting force: " << expectedForce);
+ WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting torque: " << expectedTorque);
+
+ Vector3<real_t> actingForce(real_t(0));
+ Vector3<real_t> actingTorque(real_t(0));
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ actingForce += bodyIt->getForce();
+ actingTorque += bodyIt->getTorque();
+
+ }
+ }
+
+ mpi::allReduceInplace(actingForce[0], mpi::SUM);
+ mpi::allReduceInplace(actingForce[1], mpi::SUM);
+ mpi::allReduceInplace(actingForce[2], mpi::SUM);
+
+ mpi::allReduceInplace(actingTorque[0], mpi::SUM);
+ mpi::allReduceInplace(actingTorque[1], mpi::SUM);
+ mpi::allReduceInplace(actingTorque[2], mpi::SUM);
+
+ WALBERLA_ROOT_SECTION()
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL(actingForce[0], expectedForce[0], "Mismatch in force0");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingForce[1], expectedForce[1], "Mismatch in force1");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingForce[2], expectedForce[2], "Mismatch in force2");
+
+ WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[0], expectedTorque[0], "Mismatch in torque0");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[1], expectedTorque[1], "Mismatch in torque1");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[2], expectedTorque[2], "Mismatch in torque2");
+ }
+
+ // clean up
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->markForDeletion();
+ }
+ }
+ syncCall();
+
+ container1->clear();
+ container2->clear();
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+ }
+
+ /////////////////////
+ // At block border //
+ /////////////////////
+ {
+ std::string testIdentifier("Test: sphere at block border");
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+ pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+ positionAtBlockBorder, radius, sphereMaterialID, false, true, false);
+
+ syncCall();
+
+ uint_t applicationCounter( 0 );
+
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ ++applicationCounter;
+ auto pos = bodyIt->getPosition();
+ bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+ }
+ }
+ mpi::allReduceInplace(applicationCounter, mpi::SUM);
+
+ container1->store();
+
+ resetter();
+
+ containerSwapper();
+
+ container2->setOnBodies();
+
+ Vector3<real_t> expectedForce = applicationCounter * testForce;
+ Vector3<real_t> expectedTorque = applicationCounter * ( torqueOffset % testForce );
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting force: " << expectedForce);
+ WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting torque: " << expectedTorque);
+
+ Vector3<real_t> actingForce(real_t(0));
+ Vector3<real_t> actingTorque(real_t(0));
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ actingForce += bodyIt->getForce();
+ actingTorque += bodyIt->getTorque();
+ }
+ }
+
+ mpi::allReduceInplace(actingForce[0], mpi::SUM);
+ mpi::allReduceInplace(actingForce[1], mpi::SUM);
+ mpi::allReduceInplace(actingForce[2], mpi::SUM);
+
+ mpi::allReduceInplace(actingTorque[0], mpi::SUM);
+ mpi::allReduceInplace(actingTorque[1], mpi::SUM);
+ mpi::allReduceInplace(actingTorque[2], mpi::SUM);
+
+ WALBERLA_ROOT_SECTION()
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL(actingForce[0], expectedForce[0], "Mismatch in force0");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingForce[1], expectedForce[1], "Mismatch in force1");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingForce[2], expectedForce[2], "Mismatch in force2");
+
+ WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[0], expectedTorque[0], "Mismatch in torque0");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[1], expectedTorque[1], "Mismatch in torque1");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[2], expectedTorque[2], "Mismatch in torque2");
+ }
+
+ // clean up
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->markForDeletion();
+ }
+ }
+ syncCall();
+
+ container1->clear();
+ container2->clear();
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+ }
+
+ ///////////////////////////////////////////
+ // At block border with updated position //
+ ///////////////////////////////////////////
+ {
+ std::string testIdentifier("Test: sphere at block border with updated position");
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+ pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+ positionAtBlockBorder, radius, sphereMaterialID, false, true, false);
+
+ syncCall();
+
+ uint_t applicationCounter( 0 );
+
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ ++applicationCounter;
+ auto pos = bodyIt->getPosition();
+ bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+ }
+ }
+ mpi::allReduceInplace(applicationCounter, mpi::SUM);
+
+ container1->store();
+
+ resetter();
+
+ containerSwapper();
+
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) {
+ for (auto bodyIt = pe::LocalBodyIterator::begin(*blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt) {
+ bodyIt->setPosition(positionAtBlockBorderUpdated);
+ }
+ }
+ syncCall();
+
+ container2->setOnBodies();
+
+ Vector3<real_t> expectedForce = applicationCounter * testForce;
+ Vector3<real_t> expectedTorque = applicationCounter * ( torqueOffset % testForce );
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting force: " << expectedForce);
+ WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting torque: " << expectedTorque);
+
+ Vector3<real_t> actingForce(real_t(0));
+ Vector3<real_t> actingTorque(real_t(0));
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ actingForce += bodyIt->getForce();
+ actingTorque += bodyIt->getTorque();
+ }
+ }
+
+ mpi::allReduceInplace(actingForce[0], mpi::SUM);
+ mpi::allReduceInplace(actingForce[1], mpi::SUM);
+ mpi::allReduceInplace(actingForce[2], mpi::SUM);
+
+ mpi::allReduceInplace(actingTorque[0], mpi::SUM);
+ mpi::allReduceInplace(actingTorque[1], mpi::SUM);
+ mpi::allReduceInplace(actingTorque[2], mpi::SUM);
+
+ WALBERLA_ROOT_SECTION()
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL(actingForce[0], expectedForce[0], "Mismatch in force0");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingForce[1], expectedForce[1], "Mismatch in force1");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingForce[2], expectedForce[2], "Mismatch in force2");
+
+ WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[0], expectedTorque[0], "Mismatch in torque0");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[1], expectedTorque[1], "Mismatch in torque1");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[2], expectedTorque[2], "Mismatch in torque2");
+ }
+
+ // clean up
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->markForDeletion();
+ }
+ }
+ syncCall();
+
+ container1->clear();
+ container2->clear();
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+ }
+
+ /////////////////////////////////////////////
+ // At block border with updated position 2 //
+ /////////////////////////////////////////////
+ {
+ std::string testIdentifier("Test: sphere at block border with updated position 2");
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+ pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+ positionAtBlockBorder2, radius, sphereMaterialID, false, true, false);
+
+ syncCall();
+
+ uint_t applicationCounter( 0 );
+
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ ++applicationCounter;
+ auto pos = bodyIt->getPosition();
+ bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+ }
+ }
+ mpi::allReduceInplace(applicationCounter, mpi::SUM);
+
+ container1->store();
+
+ resetter();
+
+ containerSwapper();
+
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) {
+ for (auto bodyIt = pe::LocalBodyIterator::begin(*blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt) {
+ bodyIt->setPosition(positionAtBlockBorderUpdated2);
+ }
+ }
+ syncCall();
+
+ container2->setOnBodies();
+
+ --applicationCounter; // sphere has vanished from one block
+
+ // in this case, the complete force can no longer be applied onto the body since the part from the block
+ // from which the body vanished is lost
+
+ // HOWEVER: this case will never appear in a coupled simulation since NO FORCE will be acting a body
+ // that is about to vanish from a block since it will not be mapped to the block from which it will vanish
+
+ // If you are expecting a different behavior, you need to change the container mechanism by first all-reducing the
+ // forces/torques to all processes/blocks that know this body such that every process/block has
+ // the full information and then the process/block that owns this body sets the complete force
+
+ Vector3<real_t> expectedForce = applicationCounter * testForce;
+ Vector3<real_t> expectedTorque = applicationCounter * ( torqueOffset % testForce );
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting force: " << expectedForce);
+ WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting torque: " << expectedTorque);
+
+ Vector3<real_t> actingForce(real_t(0));
+ Vector3<real_t> actingTorque(real_t(0));
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ actingForce += bodyIt->getForce();
+ actingTorque += bodyIt->getTorque();
+ }
+ }
+
+ mpi::allReduceInplace(actingForce[0], mpi::SUM);
+ mpi::allReduceInplace(actingForce[1], mpi::SUM);
+ mpi::allReduceInplace(actingForce[2], mpi::SUM);
+
+ mpi::allReduceInplace(actingTorque[0], mpi::SUM);
+ mpi::allReduceInplace(actingTorque[1], mpi::SUM);
+ mpi::allReduceInplace(actingTorque[2], mpi::SUM);
+
+ WALBERLA_ROOT_SECTION()
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL(actingForce[0], expectedForce[0], "Mismatch in force0");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingForce[1], expectedForce[1], "Mismatch in force1");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingForce[2], expectedForce[2], "Mismatch in force2");
+
+ WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[0], expectedTorque[0], "Mismatch in torque0");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[1], expectedTorque[1], "Mismatch in torque1");
+ WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[2], expectedTorque[2], "Mismatch in torque2");
+ }
+
+ // clean up
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->markForDeletion();
+ }
+ }
+ syncCall();
+
+ container1->clear();
+ container2->clear();
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+ }
+
+ return 0;
+
+}
+
+} //namespace force_torque_container_test
+
+int main( int argc, char **argv ){
+ force_torque_container_test::main(argc, argv);
+}
diff --git a/tests/pe_coupling/utility/PeSubCyclingTest.cpp b/tests/pe_coupling/utility/PeSubCyclingTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c21cdff12ee606166f60c42cc53e9a42435c2146
--- /dev/null
+++ b/tests/pe_coupling/utility/PeSubCyclingTest.cpp
@@ -0,0 +1,455 @@
+//======================================================================================================================
+//
+// 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 PeSubCyclingTest.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/math/all.h"
+
+#include "pe/basic.h"
+#include "pe/utility/DestroyBody.h"
+#include "pe/cr/DEM.h"
+
+#include <pe_coupling/utility/all.h>
+
+namespace pe_sub_cycling_test
+{
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+
+typedef boost::tuple<pe::Sphere> BodyTypeTuple ;
+
+/*!\brief test case to check functionality of sub cycling in the pe time step provided by the coupling module
+ *
+ * During this time step, currently acting forces/torques on a body are kept constant.
+ *
+ * This test first computes the resulting position offset as well as the linear and rotational velocity after a
+ * certain force has been applied to a sphere when using several 'real' pe time steps and re-applying the force 10 times.
+ *
+ * It then checks the pe time step sub-cycling functionality of the coupling module by creating same-sized spheres
+ * at different locations and applying the same force. But this time, 10 sub-cycles are carried out.
+ * As a result, the same position offset and velocities must be obtained as in the regular case.
+ *
+ */
+//////////
+// MAIN //
+//////////
+int main( int argc, char **argv )
+{
+ debug::enterTestMode();
+
+ mpi::Environment env( argc, argv );
+
+ // uncomment to have logging
+ //logging::Logging::instance()->setLogLevel(logging::Logging::LogLevel::DETAIL);
+
+ const real_t dx = real_t(1);
+ const real_t radius = real_t(5);
+
+ ///////////////////////////
+ // DATA STRUCTURES SETUP //
+ ///////////////////////////
+
+ Vector3<uint_t> blocksPerDirection(uint_t(3), uint_t(1), uint_t(1));
+ Vector3<uint_t> cellsPerBlock(uint_t(20), uint_t(20), uint_t(20));
+ Vector3<bool> periodicity(true, false, false);
+
+ // create fully periodic domain with refined blocks
+ auto blocks = blockforest::createUniformBlockGrid( blocksPerDirection[0], blocksPerDirection[1], blocksPerDirection[2],
+ cellsPerBlock[0], cellsPerBlock[1], cellsPerBlock[2],
+ dx,
+ 0, false, false,
+ periodicity[0], periodicity[1], periodicity[2],
+ false );
+
+
+ // pe body storage
+ pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
+ shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
+ auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "Storage");
+ auto sphereMaterialID = pe::createMaterial( "sphereMat", real_t(1) , real_t(0.3), real_t(0.2), real_t(0.2), real_t(0.24), real_t(200), real_t(200), real_t(0), real_t(0) );
+
+ auto ccdID = blocks->addBlockData(pe::ccd::createHashGridsDataHandling( globalBodyStorage, bodyStorageID ), "CCD");
+ auto fcdID = blocks->addBlockData(pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, pe::fcd::AnalyticCollideFunctor>(), "FCD");
+ pe::cr::DEM cr(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID, NULL);
+
+ // set up synchronization procedure
+ const real_t overlap = real_t( 1.5 ) * dx;
+ std::function<void(void)> syncCall = std::bind( pe::syncNextNeighbors<BodyTypeTuple>, boost::ref(blocks->getBlockForest()), bodyStorageID, static_cast<WcTimingTree*>(NULL), overlap, false );
+
+
+ // sphere positions for test scenarios
+ Vector3<real_t> positionInsideBlock(real_t(10), real_t(10), real_t(10));
+ Vector3<real_t> positionAtBlockBorder(real_t(19.9), real_t(10), real_t(10));
+ Vector3<real_t> positionAtBlockBorder2(real_t(20) + radius + overlap - real_t(0.1), real_t(10), real_t(10));
+
+ Vector3<real_t> testForce(real_t(2), real_t(1), real_t(0));
+ Vector3<real_t> torqueOffset = Vector3<real_t>(real_t(1), real_t(0), real_t(0));
+
+ uint_t peSubCycles( 10 );
+ real_t dtPe( real_t(10) );
+ real_t dtPeSubCycle = dtPe / real_c(peSubCycles);
+
+ pe_coupling::TimeStep timestep(blocks, bodyStorageID, cr, syncCall, dtPe, peSubCycles);
+
+ // evaluate how far the sphere will travel with a specific force applied which is the reference distance for later
+ // (depends on the chosen time integrator in the DEM and thus can not generally be computed a priori here)
+
+ Vector3<real_t> expectedPosOffset(real_t(0));
+ Vector3<real_t> expectedLinearVel(real_t(0));
+ Vector3<real_t> expectedAngularVel(real_t(0));
+ {
+
+ Vector3<real_t> startPos = positionInsideBlock;
+
+ pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+ startPos, radius, sphereMaterialID, false, true, false);
+
+ for( uint_t t = 0; t < peSubCycles; ++t )
+ {
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->addForceAtPos(testForce, bodyIt->getPosition() + torqueOffset);
+ }
+ }
+ cr.timestep(dtPeSubCycle);
+ syncCall();
+ }
+
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+ {
+ expectedPosOffset = bodyIt->getPosition() - startPos;
+ expectedLinearVel = bodyIt->getLinearVel();
+ expectedAngularVel = bodyIt->getAngularVel();
+ }
+ }
+
+ mpi::allReduceInplace(expectedPosOffset[0], mpi::SUM);
+ mpi::allReduceInplace(expectedPosOffset[1], mpi::SUM);
+ mpi::allReduceInplace(expectedPosOffset[2], mpi::SUM);
+ mpi::allReduceInplace(expectedLinearVel[0], mpi::SUM);
+ mpi::allReduceInplace(expectedLinearVel[1], mpi::SUM);
+ mpi::allReduceInplace(expectedLinearVel[2], mpi::SUM);
+ mpi::allReduceInplace(expectedAngularVel[0], mpi::SUM);
+ mpi::allReduceInplace(expectedAngularVel[1], mpi::SUM);
+ mpi::allReduceInplace(expectedAngularVel[2], mpi::SUM);
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting position offset: " << expectedPosOffset);
+ WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting linear vel: " << expectedLinearVel);
+ WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting angular vel: " << expectedAngularVel);
+
+ // clean up
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->markForDeletion();
+ }
+ }
+ syncCall();
+
+ }
+
+ //////////////////
+ // Inside block //
+ //////////////////
+ {
+ std::string testIdentifier("Test: sphere inside block");
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+
+ Vector3<real_t> startPos = positionInsideBlock;
+
+ pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+ startPos, radius, sphereMaterialID, false, true, false);
+
+ syncCall();
+
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+ {
+ auto pos = bodyIt->getPosition();
+ bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+ }
+ }
+
+ timestep();
+
+
+ Vector3<real_t> curPosOffset(real_t(0));
+ Vector3<real_t> curLinearVel(real_t(0));
+ Vector3<real_t> curAngularVel(real_t(0));
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+ {
+ curPosOffset = bodyIt->getPosition() - startPos;
+ curLinearVel = bodyIt->getLinearVel();
+ curAngularVel = bodyIt->getAngularVel();
+
+ WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[0], expectedPosOffset[0], "Mismatch in posOffset0");
+ WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[1], expectedPosOffset[1], "Mismatch in posOffset1");
+ WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[2], expectedPosOffset[2], "Mismatch in posOffset2");
+
+ WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[0], expectedLinearVel[0], "Mismatch in linearVel0");
+ WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[1], expectedLinearVel[1], "Mismatch in linearVel1");
+ WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[2], expectedLinearVel[2], "Mismatch in linearVel2");
+
+ WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[0], expectedAngularVel[0], "Mismatch in angularVel0");
+ WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[1], expectedAngularVel[1], "Mismatch in angularVel1");
+ WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[2], expectedAngularVel[2], "Mismatch in angularVel2");
+
+
+ }
+ }
+
+ // clean up
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->markForDeletion();
+ }
+ }
+ syncCall();
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+ }
+
+ ///////////////////////
+ // At block border 1 //
+ ///////////////////////
+ {
+ std::string testIdentifier("Test: sphere at block border 1");
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+
+ Vector3<real_t> startPos = positionAtBlockBorder;
+
+ pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+ startPos, radius, sphereMaterialID, false, true, false);
+
+ syncCall();
+
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+ {
+ auto pos = bodyIt->getPosition();
+ bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+ }
+ }
+
+ timestep();
+
+
+ Vector3<real_t> curPosOffset(real_t(0));
+ Vector3<real_t> curLinearVel(real_t(0));
+ Vector3<real_t> curAngularVel(real_t(0));
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+ {
+ curPosOffset = bodyIt->getPosition() - startPos;
+ curLinearVel = bodyIt->getLinearVel();
+ curAngularVel = bodyIt->getAngularVel();
+
+ WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[0], expectedPosOffset[0], "Mismatch in posOffset0");
+ WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[1], expectedPosOffset[1], "Mismatch in posOffset1");
+ WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[2], expectedPosOffset[2], "Mismatch in posOffset2");
+
+ WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[0], expectedLinearVel[0], "Mismatch in linearVel0");
+ WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[1], expectedLinearVel[1], "Mismatch in linearVel1");
+ WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[2], expectedLinearVel[2], "Mismatch in linearVel2");
+
+ WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[0], expectedAngularVel[0], "Mismatch in angularVel0");
+ WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[1], expectedAngularVel[1], "Mismatch in angularVel1");
+ WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[2], expectedAngularVel[2], "Mismatch in angularVel2");
+
+ }
+ }
+
+ // clean up
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->markForDeletion();
+ }
+ }
+ syncCall();
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+ }
+
+ ////////////////////////////
+ // At block border 1, mod //
+ ////////////////////////////
+ {
+ std::string testIdentifier("Test: sphere at block border 1 mod");
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+
+ Vector3<real_t> startPos = positionAtBlockBorder;
+
+ pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+ startPos, radius, sphereMaterialID, false, true, false);
+
+ syncCall();
+
+ // also add on shadow copy, but only half of it to have same total force/torque
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ auto pos = bodyIt->getPosition();
+ bodyIt->addForceAtPos(real_t(0.5)*testForce, pos+torqueOffset);
+ }
+ }
+
+ timestep();
+
+
+ Vector3<real_t> curPosOffset(real_t(0));
+ Vector3<real_t> curLinearVel(real_t(0));
+ Vector3<real_t> curAngularVel(real_t(0));
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+ {
+ curPosOffset = bodyIt->getPosition() - startPos;
+ curLinearVel = bodyIt->getLinearVel();
+ curAngularVel = bodyIt->getAngularVel();
+
+ WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[0], expectedPosOffset[0], "Mismatch in posOffset0");
+ WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[1], expectedPosOffset[1], "Mismatch in posOffset1");
+ WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[2], expectedPosOffset[2], "Mismatch in posOffset2");
+
+ WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[0], expectedLinearVel[0], "Mismatch in linearVel0");
+ WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[1], expectedLinearVel[1], "Mismatch in linearVel1");
+ WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[2], expectedLinearVel[2], "Mismatch in linearVel2");
+
+ WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[0], expectedAngularVel[0], "Mismatch in angularVel0");
+ WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[1], expectedAngularVel[1], "Mismatch in angularVel1");
+ WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[2], expectedAngularVel[2], "Mismatch in angularVel2");
+
+ }
+ }
+
+ // clean up
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->markForDeletion();
+ }
+ }
+ syncCall();
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+ }
+
+ ///////////////////////
+ // At block border 2 //
+ ///////////////////////
+ {
+ std::string testIdentifier("Test: sphere at block border 2");
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+
+ Vector3<real_t> startPos = positionAtBlockBorder2;
+
+ pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+ startPos, radius, sphereMaterialID, false, true, false);
+
+ syncCall();
+
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+ {
+ auto pos = bodyIt->getPosition();
+ bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+ }
+ }
+
+ timestep();
+
+
+ Vector3<real_t> curPosOffset(real_t(0));
+ Vector3<real_t> curLinearVel(real_t(0));
+ Vector3<real_t> curAngularVel(real_t(0));
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+ {
+ curPosOffset = bodyIt->getPosition() - startPos;
+ curLinearVel = bodyIt->getLinearVel();
+ curAngularVel = bodyIt->getAngularVel();
+
+ WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[0], expectedPosOffset[0], "Mismatch in posOffset0");
+ WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[1], expectedPosOffset[1], "Mismatch in posOffset1");
+ WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[2], expectedPosOffset[2], "Mismatch in posOffset2");
+
+ WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[0], expectedLinearVel[0], "Mismatch in linearVel0");
+ WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[1], expectedLinearVel[1], "Mismatch in linearVel1");
+ WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[2], expectedLinearVel[2], "Mismatch in linearVel2");
+
+ WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[0], expectedAngularVel[0], "Mismatch in angularVel0");
+ WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[1], expectedAngularVel[1], "Mismatch in angularVel1");
+ WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[2], expectedAngularVel[2], "Mismatch in angularVel2");
+
+ }
+ }
+
+ // clean up
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->markForDeletion();
+ }
+ }
+ syncCall();
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+ }
+
+ return 0;
+
+}
+
+} //namespace pe_sub_cycling_test
+
+int main( int argc, char **argv ){
+ pe_sub_cycling_test::main(argc, argv);
+}