diff --git a/src/pe/cr/ContactResolvers.h b/src/pe/cr/ContactResolvers.h
new file mode 100644
index 0000000000000000000000000000000000000000..d90b6783559c68628ca8f0d91f8b5288a7129d56
--- /dev/null
+++ b/src/pe/cr/ContactResolvers.h
@@ -0,0 +1,59 @@
+#pragma once
+
+#include "pe/Types.h"
+#include "ICR.h"
+#include "pe/contact/Contact.h"
+#include "pe/contact/ContactFunctions.h"
+
+namespace walberla {
+namespace pe {
+namespace cr {
+
+class ResolveContactSpringDashpotHaffWerner {
+public:
+ void operator()( ContactID c ) const
+ {
+ WALBERLA_LOG_DETAIL( "resolving contact: " << c->getID() );
+ BodyID b1( c->getBody1()->getTopSuperBody() );
+ BodyID b2( c->getBody2()->getTopSuperBody() );
+
+ // Global position of contact
+ const Vec3 gpos( c->getPosition() );
+
+ // The absolute value of the penetration length
+ real_t delta( -c->getDistance() );
+
+ // Calculating the relative velocity in normal and tangential direction
+ // The negative signs result from the different definition of the relative
+ // normal velocity of the pe (see Contact::getType)
+ const real_t relVelN( -c->getNormalRelVel() );
+ const Vec3 relVel ( -c->getRelVel() );
+ const Vec3 relVelT( relVel - ( relVelN * c->getNormal() ) );
+
+ real_t fNabs( 0 );
+ Vec3 fN;
+
+ // Calculating the normal force based on a linear spring-dashpot force model
+ fNabs = getStiffness(c) * delta + getDampingN(c) * relVelN;
+ if( fNabs < real_c(0) ) fNabs = real_c(0);
+ fN = fNabs * c->getNormal();
+
+ // Calculating the tangential force based on the model by Haff and Werner
+ const real_t fTabs( std::min( getDampingT(c) * relVelT.length(), getFriction(c) * fNabs ) );
+ const Vec3 fT ( fTabs * relVelT.getNormalizedOrZero() );
+
+ // Add normal force at contact point
+ b1->addForceAtPos( fN, gpos );
+ b2->addForceAtPos( -fN, gpos );
+
+ // Add tangential force at contact point
+ b1->addForceAtPos( fT, gpos );
+ b2->addForceAtPos( -fT, gpos );
+
+ WALBERLA_LOG_DETAIL("nForce: " << fN << "\ntForce: " << fT);
+ }
+};
+
+} // namespace cr
+} // namespace pe
+} // namespace walberla
diff --git a/src/pe/cr/DEM.cpp b/src/pe/cr/DEM.cpp
deleted file mode 100644
index e4915c755efb69c49e0d9499db4e98f5ca1f90ad..0000000000000000000000000000000000000000
--- a/src/pe/cr/DEM.cpp
+++ /dev/null
@@ -1,282 +0,0 @@
-//======================================================================================================================
-//
-// 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 DEM.cpp
-//! \author Klaus Iglberger
-//! \author Sebastian Eibl <sebastian.eibl@fau.de>
-//! \brief Source file for the DEM solver
-//
-//======================================================================================================================
-
-
-//*************************************************************************************************
-// Includes
-//*************************************************************************************************
-
-#include "pe/cr/DEM.h"
-#include "pe/ccd/ICCD.h"
-#include "pe/fcd/IFCD.h"
-#include "pe/bg/IBG.h"
-#include "pe/rigidbody/BodyStorage.h"
-#include "pe/rigidbody/RigidBody.h"
-#include "pe/contact/Contact.h"
-#include "pe/contact/ContactFunctions.h"
-#include "pe/synchronization/SyncForces.h"
-
-#include "core/logging/all.h"
-
-namespace walberla {
-namespace pe {
-namespace cr {
-
-DEM::DEM( const shared_ptr<BodyStorage>& globalBodyStorage
- , const shared_ptr<BlockStorage>& blockStorage
- , domain_decomposition::BlockDataID storageID
- , domain_decomposition::BlockDataID ccdID
- , domain_decomposition::BlockDataID fcdID
- , WcTimingTree* tt)
- : globalBodyStorage_(globalBodyStorage)
- , blockStorage_(blockStorage)
- , storageID_(storageID)
- , ccdID_(ccdID)
- , fcdID_(fcdID)
- , tt_(tt)
- , maxPenetration_(0)
- , numberOfContacts_(0)
- , numberOfContactsTreated_(0)
-{
-
-}
-
-void DEM::timestep( real_t dt )
-{
- maxPenetration_ = real_c(0.0);
- numberOfContacts_ = 0;
- numberOfContactsTreated_ = 0;
-
- for (auto it = blockStorage_->begin(); it != blockStorage_->end(); ++it){
- IBlock & currentBlock = *it;
-
- ccd::ICCD* ccd = currentBlock.getData< ccd::ICCD >( ccdID_ );
- fcd::IFCD* fcd = currentBlock.getData< fcd::IFCD >( fcdID_ );
- ccd->generatePossibleContacts( tt_ );
- if (tt_ != NULL) tt_->start("FCD");
- Contacts& cont = fcd->generateContacts( ccd->getPossibleContacts() );
- if (tt_ != NULL) tt_->stop("FCD");
-
- for (auto cIt = cont.begin(); cIt != cont.end(); ++cIt){
- const real_t overlap( -cIt->getDistance() );
- if( overlap > maxPenetration_ )
- maxPenetration_ = overlap;
- if (shouldContactBeTreated( &(*cIt), currentBlock.getAABB() ))
- {
- ++numberOfContactsTreated_;
- resolveContact( &(*cIt) );
- }
- }
-
- numberOfContacts_ += cont.size();
-
- cont.clear();
- }
-
-// if (numContacts > 0)
-// WALBERLA_LOG_DEVEL_ON_ROOT("#Contacts: " << numContacts << "." );
-
- if (tt_ != NULL) tt_->start("ForceSync");
- reduceForces( *blockStorage_, storageID_, *globalBodyStorage_);
- if (tt_ != NULL) tt_->stop("ForceSync");
-
- for (auto it = blockStorage_->begin(); it != blockStorage_->end(); ++it)
- {
- IBlock & currentBlock = *it;
- Storage * storage = currentBlock.getData< Storage >( storageID_ );
- BodyStorage& localStorage = (*storage)[0];
- BodyStorage& shadowStorage = (*storage)[1];
-
- // Updating the positions and velocities of all locally owned and global rigid bodies (but not shadow copies).
- WALBERLA_LOG_DETAIL( "Time integration starts...\n" );
-
- if (tt_ != NULL) tt_->start("Integration");
-
- for( auto bodyIt = localStorage.begin(); bodyIt != localStorage.end(); ++bodyIt )
- {
- WALBERLA_LOG_DETAIL( "Time integration of body with system id " << bodyIt->getSystemID());// << "\n" << *bodyIt );
-
- // Resetting the contact node and removing all attached contacts
- // bodyIt->resetNode();
- bodyIt->clearContacts();
-
- move( *bodyIt, dt );
-
- // Resetting the acting forces
- bodyIt->resetForceAndTorque();
- }
-
- if (tt_ != NULL) tt_->stop("Integration");
-
- // Reset forces of shadow copies
- for( auto bodyIt = shadowStorage.begin(); bodyIt != shadowStorage.end(); ++bodyIt ) {
- bodyIt->clearContacts();
- bodyIt->resetForceAndTorque();
- }
-
- }
-}
-
-//=================================================================================================
-//
-// SOLVER FUNCTIONS
-//
-//=================================================================================================
-
-//*************************************************************************************************
-/*!\brief Resolves the given colliding contact.
- *
- * \param c The colliding contact.
- * \return The overlap of the contact.
- *
- * TODO
- */
-void DEM::resolveContact( ContactID c ) const
-{
- WALBERLA_LOG_DETAIL( "resolving contact: " << c->getID() );
- BodyID b1( c->getBody1()->getTopSuperBody() );
- BodyID b2( c->getBody2()->getTopSuperBody() );
-
- // Global position of contact
- const Vec3 gpos( c->getPosition() );
-
- // The absolute value of the penetration length
- real_t delta( -c->getDistance() );
-
- // Calculating the relative velocity in normal and tangential direction
- // The negative signs result from the different definition of the relative
- // normal velocity of the pe (see Contact::getType)
- const real_t relVelN( -c->getNormalRelVel() );
- const Vec3 relVel ( -c->getRelVel() );
- const Vec3 relVelT( relVel - ( relVelN * c->getNormal() ) );
-
- real_t fNabs( 0 );
- Vec3 fN;
-
- // Calculating the normal force based on the non-linear extended Hertz model
- // This force model is only applied in case of a sphere/sphere collision, since only
- // then an effective radius can be computed.
-// if( dem::forceModel == dem::hertz && c->hasEffectiveRadius() )
-// {
-// const real_t alpha ( 1.5 );
-// const real_t beta ( 0.5 );
-// const real_t Reff ( c->getEffectiveRadius() );
-// const real_t Eeff ( c->getEffectiveYoungModulus() );
-// const real_t k ( ( real_c(4)/real_c(3) ) * Eeff * sqrt( Reff ) );
-
-// fNabs = k*std::pow( delta, alpha ) + c->getDampingN()*relVelN*std::pow( delta, beta );
-// if( fNabs < real_c(0) ) fNabs = real_c(0);
-// fN = fNabs * c->getNormal();
-// }
-
-// // Calculating the normal force based on a linear spring-dashpot force model
-// else
- {
- fNabs = getStiffness(c) * delta + getDampingN(c) * relVelN;
- if( fNabs < real_c(0) ) fNabs = real_c(0);
- fN = fNabs * c->getNormal();
- }
-
- // Calculating the tangential force based on the model by Haff and Werner
- const real_t fTabs( std::min( getDampingT(c) * relVelT.length(), getFriction(c) * fNabs ) );
- const Vec3 fT ( fTabs * relVelT.getNormalizedOrZero() );
-
- // Add normal force at contact point
- b1->addForceAtPos( fN, gpos );
- b2->addForceAtPos( -fN, gpos );
-
- // Add tangential force at contact point
- b1->addForceAtPos( fT, gpos );
- b2->addForceAtPos( -fT, gpos );
-
- WALBERLA_LOG_DETAIL("nForce: " << fN << "\ntForce: " << fT);
-}
-//*************************************************************************************************
-
-//=================================================================================================
-//
-// SIMULATION FUNCTIONS
-//
-//=================================================================================================
-
-//*************************************************************************************************
-/*!\brief Calculation of a time step of \a dt.
- *
- * \param dt Time step size.
- * \return void
- *
- * Calculating one single time step of size \a dt for the capsule. The global position, the
- * linear and angular velocity and the orientation of the capsule are changed depending on
- * the acting forces and the current velocities.
- */
-void DEM::move( BodyID id, real_t dt )
-{
- // Checking the state of the body
- WALBERLA_ASSERT( id->checkInvariants(), "Invalid body state detected" );
- WALBERLA_ASSERT( !id->hasSuperBody(), "Invalid superordinate body detected" );
-
- // Moving the capsule according to the acting forces (don't move a sleeping capsule)
- if( id->isAwake() && !id->hasInfiniteMass() )
- {
- // Calculating the linear acceleration by the equation
- // force * m^(-1) + gravity
- const Vec3 vdot( id->getForce() * id->getInvMass() + getGlobalLinearAcceleration() );
-
- // Calculating the angular acceleration by the equation
- // R * Iinv * R^T * torque
- const Vec3 wdot( id->getInvInertia() * id->getTorque() );
-
- // Calculating the translational displacement
- id->setPosition( id->getPosition() + id->getLinearVel() * dt + 0.5 * vdot * dt * dt );
-
- // Calculating the rotation angle
- const Vec3 phi( id->getAngularVel() * dt + 0.5 * wdot * dt * dt);
-
- // Calculating the new orientation
- if (!floatIsEqual(phi.length(), 0))
- id->rotate( Quat( phi, phi.length() ) );
- WALBERLA_ASSERT_FLOAT_EQUAL( id->getRotation().getDeterminant(), real_c(1 ), "Corrupted rotation matrix determinant" );
-
- // Updating the linear velocity
- id->setLinearVel(id->getLinearVel() + vdot * dt);
-
- // Updating the angular velocity
- id->setAngularVel( id->getAngularVel() + wdot * dt );
-
- // Setting the axis-aligned bounding box
- id->calcBoundingBox();
-
- // Calculating the current motion of the capsule
- id->calcMotion();
- }
-
- // Resetting the acting forces
- id->resetForceAndTorque();
-
- // Checking the state of the rigid body
- WALBERLA_ASSERT( id->checkInvariants(), "Invalid capsule state detected" );
-}
-//*************************************************************************************************
-
-} // namespace cr
-} // namespace pe
-} // namespace walberla
diff --git a/src/pe/cr/DEM.h b/src/pe/cr/DEM.h
index 100d921ceb2f46b2f2095ec378c7709c67afe22a..f5269b06fa53a1fae90b286532429f0720507c5e 100644
--- a/src/pe/cr/DEM.h
+++ b/src/pe/cr/DEM.h
@@ -28,6 +28,8 @@
//*************************************************************************************************
#include "ICR.h"
+#include "Integrators.h"
+#include "ContactResolvers.h"
#include "pe/Types.h"
#include "domain_decomposition/BlockStorage.h"
@@ -39,15 +41,17 @@ namespace cr {
/**
* \ingroup pe
*/
-class DEM : public ICR
+template< typename Integrator, typename ContactResolver >
+class DEMSolver : public ICR
{
public:
- DEM( const shared_ptr<BodyStorage>& globalBodyStorage
- , const shared_ptr<BlockStorage>& blockStorage
- , domain_decomposition::BlockDataID storageID
- , domain_decomposition::BlockDataID ccdID
- , domain_decomposition::BlockDataID fcdID
- , WcTimingTree* tt = NULL);
+ DEMSolver( const Integrator & integrate, const ContactResolver & resolveContact
+ , const shared_ptr<BodyStorage>& globalBodyStorage
+ , const shared_ptr<BlockStorage>& blockStorage
+ , domain_decomposition::BlockDataID storageID
+ , domain_decomposition::BlockDataID ccdID
+ , domain_decomposition::BlockDataID fcdID
+ , WcTimingTree* tt = NULL);
/// forwards to timestep
/// Convenience operator to make class a functor.
@@ -59,9 +63,8 @@ public:
virtual inline size_t getNumberOfContacts() const WALBERLA_OVERRIDE { return numberOfContacts_; }
virtual inline size_t getNumberOfContactsTreated() const WALBERLA_OVERRIDE { return numberOfContactsTreated_; }
private:
- void resolveContact( ContactID c ) const;
- void move( BodyID id, real_t dt );
-
+ const Integrator integrate_;
+ const ContactResolver resolveContact_;
shared_ptr<BodyStorage> globalBodyStorage_;
shared_ptr<BlockStorage> blockStorage_;
domain_decomposition::BlockDataID storageID_;
@@ -74,6 +77,24 @@ private:
size_t numberOfContactsTreated_;
};
+class DEM : public DEMSolver<IntegrateImplictEuler, ResolveContactSpringDashpotHaffWerner>
+{
+public:
+ DEM( const shared_ptr<BodyStorage>& globalBodyStorage
+ , const shared_ptr<BlockStorage>& blockStorage
+ , domain_decomposition::BlockDataID storageID
+ , domain_decomposition::BlockDataID ccdID
+ , domain_decomposition::BlockDataID fcdID
+ , WcTimingTree* tt = NULL)
+ : DEMSolver<IntegrateImplictEuler, ResolveContactSpringDashpotHaffWerner>(
+ IntegrateImplictEuler(), ResolveContactSpringDashpotHaffWerner(),
+ globalBodyStorage, blockStorage, storageID, ccdID, fcdID, tt )
+ {
+ }
+};
+
} // namespace cr
} // namespace pe
} // namespace walberla
+
+#include "DEM.impl.h"
diff --git a/src/pe/cr/DEM.impl.h b/src/pe/cr/DEM.impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..3a910f72cacde157d4d50cbbba716b43e32f5aed
--- /dev/null
+++ b/src/pe/cr/DEM.impl.h
@@ -0,0 +1,160 @@
+//======================================================================================================================
+//
+// 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 DEM.cpp
+//! \author Klaus Iglberger
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//! \brief Source file for the DEM solver
+//
+//======================================================================================================================
+
+
+//*************************************************************************************************
+// Includes
+//*************************************************************************************************
+
+#include "pe/cr/DEM.h"
+#include "pe/ccd/ICCD.h"
+#include "pe/fcd/IFCD.h"
+#include "pe/bg/IBG.h"
+#include "pe/rigidbody/BodyStorage.h"
+#include "pe/rigidbody/RigidBody.h"
+#include "pe/contact/Contact.h"
+#include "pe/contact/ContactFunctions.h"
+#include "pe/synchronization/SyncForces.h"
+
+#include "core/logging/all.h"
+
+namespace walberla {
+namespace pe {
+namespace cr {
+
+template< typename Integrator, typename ContactResolver >
+DEMSolver<Integrator,ContactResolver>::DEMSolver(
+ const Integrator & integrate, const ContactResolver & resolveContact
+ , const shared_ptr<BodyStorage>& globalBodyStorage
+ , const shared_ptr<BlockStorage>& blockStorage
+ , domain_decomposition::BlockDataID storageID
+ , domain_decomposition::BlockDataID ccdID
+ , domain_decomposition::BlockDataID fcdID
+ , WcTimingTree* tt)
+ : integrate_(integrate), resolveContact_(resolveContact)
+ ,globalBodyStorage_(globalBodyStorage)
+ , blockStorage_(blockStorage)
+ , storageID_(storageID)
+ , ccdID_(ccdID)
+ , fcdID_(fcdID)
+ , tt_(tt)
+ , maxPenetration_(0)
+ , numberOfContacts_(0)
+ , numberOfContactsTreated_(0)
+{
+
+}
+
+template< typename Integrator, typename ContactResolver >
+void DEMSolver<Integrator,ContactResolver>::timestep( real_t dt )
+{
+ maxPenetration_ = real_c(0.0);
+ numberOfContacts_ = 0;
+ numberOfContactsTreated_ = 0;
+
+ for (auto it = blockStorage_->begin(); it != blockStorage_->end(); ++it){
+ IBlock & currentBlock = *it;
+
+ ccd::ICCD* ccd = currentBlock.getData< ccd::ICCD >( ccdID_ );
+ fcd::IFCD* fcd = currentBlock.getData< fcd::IFCD >( fcdID_ );
+ ccd->generatePossibleContacts( tt_ );
+ if (tt_ != NULL) tt_->start("FCD");
+ Contacts& cont = fcd->generateContacts( ccd->getPossibleContacts() );
+ if (tt_ != NULL) tt_->stop("FCD");
+
+ for (auto cIt = cont.begin(); cIt != cont.end(); ++cIt){
+ const real_t overlap( -cIt->getDistance() );
+ if( overlap > maxPenetration_ )
+ maxPenetration_ = overlap;
+ if (shouldContactBeTreated( &(*cIt), currentBlock.getAABB() ))
+ {
+ ++numberOfContactsTreated_;
+ resolveContact_( &(*cIt) );
+ }
+ }
+
+ numberOfContacts_ += cont.size();
+
+ cont.clear();
+ }
+
+// if (numContacts > 0)
+// WALBERLA_LOG_DEVEL_ON_ROOT("#Contacts: " << numContacts << "." );
+
+ if (tt_ != NULL) tt_->start("ForceSync");
+ reduceForces( *blockStorage_, storageID_, *globalBodyStorage_);
+ if (tt_ != NULL) tt_->stop("ForceSync");
+
+ for (auto it = blockStorage_->begin(); it != blockStorage_->end(); ++it)
+ {
+ IBlock & currentBlock = *it;
+ Storage * storage = currentBlock.getData< Storage >( storageID_ );
+ BodyStorage& localStorage = (*storage)[0];
+ BodyStorage& shadowStorage = (*storage)[1];
+
+ // Updating the positions and velocities of all locally owned and global rigid bodies (but not shadow copies).
+ WALBERLA_LOG_DETAIL( "Time integration starts...\n" );
+
+ if (tt_ != NULL) tt_->start("Integration");
+
+ for( auto bodyIt = localStorage.begin(); bodyIt != localStorage.end(); ++bodyIt )
+ {
+ WALBERLA_LOG_DETAIL( "Time integration of body with system id " << bodyIt->getSystemID());// << "\n" << *bodyIt );
+
+ // Resetting the contact node and removing all attached contacts
+ // bodyIt->resetNode();
+ bodyIt->clearContacts();
+
+ // Checking the state of the body
+ WALBERLA_ASSERT( bodyIt->checkInvariants(), "Invalid body state detected" );
+ WALBERLA_ASSERT( !bodyIt->hasSuperBody(), "Invalid superordinate body detected" );
+
+ // Moving the capsule according to the acting forces (don't move a sleeping capsule)
+ if( bodyIt->isAwake() && !bodyIt->hasInfiniteMass() )
+ {
+ integrate_( *bodyIt, dt, *this );
+ }
+
+ // Resetting the acting forces
+ bodyIt->resetForceAndTorque();
+
+ // Checking the state of the rigid body
+ WALBERLA_ASSERT( bodyIt->checkInvariants(), "Invalid capsule state detected" );
+
+ // Resetting the acting forces
+ bodyIt->resetForceAndTorque();
+ }
+
+ if (tt_ != NULL) tt_->stop("Integration");
+
+ // Reset forces of shadow copies
+ for( auto bodyIt = shadowStorage.begin(); bodyIt != shadowStorage.end(); ++bodyIt ) {
+ bodyIt->clearContacts();
+ bodyIt->resetForceAndTorque();
+ }
+
+ }
+}
+
+} // namespace cr
+} // namespace pe
+} // namespace walberla
diff --git a/src/pe/cr/Integrators.h b/src/pe/cr/Integrators.h
new file mode 100644
index 0000000000000000000000000000000000000000..7675caf6d51f7ce27b6f015e68bafdb14654ac9d
--- /dev/null
+++ b/src/pe/cr/Integrators.h
@@ -0,0 +1,108 @@
+#pragma once
+
+#include "pe/Types.h"
+#include "ICR.h"
+
+namespace walberla {
+namespace pe {
+namespace cr {
+
+//*************************************************************************************************
+/*!\brief Integrate the trajectory of one body using implict Euler.
+*
+* \param id Body ID.
+* \param dt Time step size.
+* \param solver Solver.
+* \return void
+*
+* The implicit Euler algorithm, also known as backward Euler, is used. It is a first-order
+* integrator that does conserves energy (i.e. it is symplectic.)
+*/
+class IntegrateImplictEuler {
+public:
+ void operator()( BodyID id, real_t dt, ICR & solver ) const
+ {
+ // Calculating the linear acceleration by the equation
+ // force * m^(-1) + gravity
+ const Vec3 vdot( id->getForce() * id->getInvMass() + solver.getGlobalLinearAcceleration() );
+
+ // Calculating the angular acceleration by the equation
+ // R * Iinv * R^T * torque
+ const Vec3 wdot( id->getInvInertia() * id->getTorque() );
+
+ // Updating the linear velocity
+ id->setLinearVel( id->getLinearVel() + vdot * dt );
+
+ // Updating the angular velocity
+ id->setAngularVel( id->getAngularVel() + wdot * dt );
+
+ // Calculating the translational displacement
+ id->setPosition( id->getPosition() + id->getLinearVel() * dt );
+
+ // Calculating the rotation angle
+ const Vec3 phi( id->getAngularVel() * dt );
+
+ // Calculating the new orientation
+ if (!floatIsEqual(phi.length(), 0))
+ id->rotate( Quat( phi, phi.length() ) );
+ WALBERLA_ASSERT_FLOAT_EQUAL( id->getRotation().getDeterminant(), real_c(1 ), "Corrupted rotation matrix determinant" );
+
+ // Setting the axis-aligned bounding box
+ id->calcBoundingBox();
+
+ // Calculating the current motion of the capsule
+ id->calcMotion();
+ }
+};
+
+//*************************************************************************************************
+/*!\brief Integrate the trajectory of one body using explicit Euler.
+*
+* \param id Body ID.
+* \param dt Time step size.
+* \param solver Solver.
+* \return void
+*
+* The explicit Euler algorithm, also known as forward Euler, is used. It is a first-order
+* integrator that does not conserve energy (i.e. it is not symplectic.)
+*/
+class IntegrateExplictEuler {
+public:
+ void operator()( BodyID id, real_t dt, ICR & solver ) const
+ {
+ // Calculating the linear acceleration by the equation
+ // force * m^(-1) + gravity
+ const Vec3 vdot( id->getForce() * id->getInvMass() + solver.getGlobalLinearAcceleration() );
+
+ // Calculating the angular acceleration by the equation
+ // R * Iinv * R^T * torque
+ const Vec3 wdot( id->getInvInertia() * id->getTorque() );
+
+ // Calculating the translational displacement
+ id->setPosition( id->getPosition() + id->getLinearVel() * dt + 0.5 * vdot * dt * dt );
+
+ // Calculating the rotation angle
+ const Vec3 phi( id->getAngularVel() * dt + 0.5 * wdot * dt * dt);
+
+ // Calculating the new orientation
+ if (!floatIsEqual(phi.length(), 0))
+ id->rotate( Quat( phi, phi.length() ) );
+ WALBERLA_ASSERT_FLOAT_EQUAL( id->getRotation().getDeterminant(), real_c(1 ), "Corrupted rotation matrix determinant" );
+
+ // Updating the linear velocity
+ id->setLinearVel( id->getLinearVel() + vdot * dt );
+
+ // Updating the angular velocity
+ id->setAngularVel( id->getAngularVel() + wdot * dt );
+
+ // Setting the axis-aligned bounding box
+ id->calcBoundingBox();
+
+ // Calculating the current motion of the capsule
+ id->calcMotion();
+ }
+};
+
+} // namespace cr
+} // namespace pe
+} // namespace walberla
diff --git a/src/pe/cr/PlainIntegrator.cpp b/src/pe/cr/PlainIntegrator.cpp
deleted file mode 100644
index 0267e80225ef58785794268c2e9f29e2fdb0f928..0000000000000000000000000000000000000000
--- a/src/pe/cr/PlainIntegrator.cpp
+++ /dev/null
@@ -1,139 +0,0 @@
-//======================================================================================================================
-//
-// 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 PlainIntegrator.h
-//! \author Klaus Iglberger
-//! \author Sebastian Eibl <sebastian.eibl@fau.de>
-//! \brief Source file for the Plain Integrator solver
-//
-//======================================================================================================================
-
-
-//*************************************************************************************************
-// Includes
-//*************************************************************************************************
-
-#include "pe/cr/PlainIntegrator.h"
-#include "pe/ccd/ICCD.h"
-#include "pe/fcd/IFCD.h"
-#include "pe/bg/IBG.h"
-#include "pe/rigidbody/BodyStorage.h"
-#include "pe/rigidbody/RigidBody.h"
-
-#include "core/logging/all.h"
-#include "core/timing/TimingTree.h"
-
-namespace walberla {
-namespace pe {
-namespace cr {
-
-PlainIntegrator::PlainIntegrator(const shared_ptr<BodyStorage>& globalBodyStorage,
- const shared_ptr<BlockStorage>& blockStorage,
- domain_decomposition::BlockDataID storageID,
- WcTimingTree* tt)
- : globalBodyStorage_(globalBodyStorage)
- , blockStorage_(blockStorage)
- , storageID_(storageID)
- , tt_(tt)
-{
-
-}
-
-void PlainIntegrator::timestep( const real_t dt )
-{
- if (tt_!=NULL) tt_->start("PlainIntegrator");
- if (tt_!=NULL) tt_->start("Integrate Bodies");
- for (auto it = blockStorage_->begin(); it != blockStorage_->end(); ++it){
- IBlock & currentBlock = *it;
- Storage * storage = currentBlock.getData< Storage >( storageID_ );
- BodyStorage& localStorage = (*storage)[0];
- for (auto bd = localStorage.begin(); bd != localStorage.end(); ++bd){
- move(*bd, dt);
- }
- }
- if (tt_!=NULL) tt_->stop("Integrate Bodies");
- if (tt_!=NULL) tt_->stop("PlainIntegrator");
-}
-
-//=================================================================================================
-//
-// SIMULATION FUNCTIONS
-//
-//=================================================================================================
-
-//*************************************************************************************************
-/*!\brief Calculation of a time step of \a dt.
- *
- * \param dt Time step size.
- * \return void
- *
- * Calculating one single time step of size \a dt for the capsule. The global position, the
- * linear and angular velocity and the orientation of the capsule are changed depending on
- * the acting forces and the current velocities.
- */
-void PlainIntegrator::move( BodyID id, real_t dt )
-{
- // Checking the state of the body
- WALBERLA_ASSERT( id->checkInvariants(), "Invalid capsule state detected" );
- WALBERLA_ASSERT( !id->hasSuperBody(), "Invalid superordinate body detected" );
-
- // Resetting the contact node and removing all attached contacts
-// id->resetNode();
- id->clearContacts();
-
- // Moving the capsule according to the acting forces (don't move a sleeping capsule)
- if( id->isAwake() && !id->hasInfiniteMass() ) {
- // Calculating the linear acceleration by the equation
- // force * m^(-1) + gravity
- const Vec3 vdot( id->getForce() * id->getInvMass() + getGlobalLinearAcceleration() );
-
- // Calculating the angular acceleration by the equation
- // R * Iinv * R^T * torque
- const Vec3 wdot( id->getInvInertia() * id->getTorque() );
-
- // Updating the linear velocity
- id->setLinearVel(id->getLinearVel() + vdot * dt);
-
- // Updating the angular velocity
- id->setAngularVel( id->getAngularVel() + wdot * dt );
-
- // Calculating the translational displacement
- id->setPosition( id->getPosition() + id->getLinearVel() * dt );
-
- // Calculating the rotation angle
- const Vec3 phi( id->getAngularVel() * dt );
-
- // Calculating the new orientation
- id->rotate( Quat( phi, phi.length() ) );
- WALBERLA_ASSERT( realIsEqual( id->getRotation().getDeterminant(), real_c(1) ), "Corrupted rotation matrix determinant" );
-
- // Setting the axis-aligned bounding box
- id->calcBoundingBox();
-
- // Calculating the current motion of the capsule
- id->calcMotion();
- }
-
- // Resetting the acting forces
- id->resetForceAndTorque();
-
- // Checking the state of the capsule
- WALBERLA_ASSERT( id->checkInvariants(), "Invalid capsule state detected" );
-}
-//*************************************************************************************************
-
-} // namespace cr
-} // namespace pe
-} // namespace walberla
diff --git a/src/pe/cr/PlainIntegrator.h b/src/pe/cr/PlainIntegrator.h
index 378af07f64685bed667819b0a2541d9b498ab7f5..271847f5d3d9b601452553d08d4a85bae8a4eb71 100644
--- a/src/pe/cr/PlainIntegrator.h
+++ b/src/pe/cr/PlainIntegrator.h
@@ -29,6 +29,7 @@
#include "ICR.h"
#include "pe/Types.h"
+#include "Integrators.h"
#include "domain_decomposition/BlockStorage.h"
@@ -36,14 +37,16 @@ namespace walberla {
namespace pe {
namespace cr {
-class PlainIntegrator : public ICR
- , private NonCopyable
+template< typename Integrator >
+class PlainIntegratorSolver : public ICR
+ , private NonCopyable
{
public:
- PlainIntegrator( const shared_ptr<BodyStorage>& globalBodyStorage,
- const shared_ptr<BlockStorage>& blockStorage,
- domain_decomposition::BlockDataID storageID,
- WcTimingTree* tt = NULL );
+ PlainIntegratorSolver( const Integrator & integrate,
+ const shared_ptr<BodyStorage>& globalBodyStorage,
+ const shared_ptr<BlockStorage>& blockStorage,
+ domain_decomposition::BlockDataID storageID,
+ WcTimingTree* tt = NULL );
/// forwards to timestep
/// Convenience operator to make class a functor.
@@ -51,14 +54,28 @@ public:
/// Advances the simulation dt seconds.
void timestep( const real_t dt );
private:
- void move( BodyID id, real_t dt );
-
+ const Integrator integrate_;
shared_ptr<BodyStorage> globalBodyStorage_;
shared_ptr<BlockStorage> blockStorage_;
domain_decomposition::BlockDataID storageID_;
WcTimingTree* tt_;
};
+
+class PlainIntegrator : public PlainIntegratorSolver<IntegrateImplictEuler>
+{
+public:
+ PlainIntegrator( const shared_ptr<BodyStorage>& globalBodyStorage
+ , const shared_ptr<BlockStorage>& blockStorage
+ , domain_decomposition::BlockDataID storageID
+ , WcTimingTree* tt = NULL)
+ : PlainIntegratorSolver<IntegrateImplictEuler>( IntegrateImplictEuler(), globalBodyStorage, blockStorage,
+ storageID, tt )
+ {
+ }
+};
} // namespace cr
} // namespace pe
} // namespace walberla
+
+#include "PlainIntegrator.impl.h"
diff --git a/src/pe/cr/PlainIntegrator.impl.h b/src/pe/cr/PlainIntegrator.impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..7b9c4ce45101e774ed936da68e81caab088f1f95
--- /dev/null
+++ b/src/pe/cr/PlainIntegrator.impl.h
@@ -0,0 +1,93 @@
+//======================================================================================================================
+//
+// 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 PlainIntegrator.h
+//! \author Klaus Iglberger
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//! \brief Source file for the Plain Integrator solver
+//
+//======================================================================================================================
+
+
+//*************************************************************************************************
+// Includes
+//*************************************************************************************************
+
+#include "pe/cr/PlainIntegrator.h"
+#include "pe/ccd/ICCD.h"
+#include "pe/fcd/IFCD.h"
+#include "pe/bg/IBG.h"
+#include "pe/rigidbody/BodyStorage.h"
+#include "pe/rigidbody/RigidBody.h"
+
+#include "core/logging/all.h"
+#include "core/timing/TimingTree.h"
+
+namespace walberla {
+namespace pe {
+namespace cr {
+
+template< typename Integrator >
+PlainIntegratorSolver<Integrator>::PlainIntegratorSolver( const Integrator & integrate,
+ const shared_ptr<BodyStorage>& globalBodyStorage,
+ const shared_ptr<BlockStorage>& blockStorage,
+ domain_decomposition::BlockDataID storageID,
+ WcTimingTree* tt)
+ : integrate_(integrate)
+ , globalBodyStorage_(globalBodyStorage)
+ , blockStorage_(blockStorage)
+ , storageID_(storageID)
+ , tt_(tt)
+{
+
+}
+
+template< typename Integrator >
+void PlainIntegratorSolver<Integrator>::timestep( const real_t dt )
+{
+ if (tt_!=NULL) tt_->start("PlainIntegrator");
+ if (tt_!=NULL) tt_->start("Integrate Bodies");
+ for (auto it = blockStorage_->begin(); it != blockStorage_->end(); ++it){
+ IBlock & currentBlock = *it;
+ Storage * storage = currentBlock.getData< Storage >( storageID_ );
+ BodyStorage& localStorage = (*storage)[0];
+ for (auto bd = localStorage.begin(); bd != localStorage.end(); ++bd){
+ // Checking the state of the body
+ WALBERLA_ASSERT( bd->checkInvariants(), "Invalid capsule state detected" );
+ WALBERLA_ASSERT( !bd->hasSuperBody(), "Invalid superordinate body detected" );
+
+ // Resetting the contact node and removing all attached contacts
+ // bd->resetNode();
+ bd->clearContacts();
+
+ // Moving the capsule according to the acting forces (don't move a sleeping capsule)
+ if( bd->isAwake() && !bd->hasInfiniteMass() ) {
+ integrate_( *bd, dt, *this );
+ }
+
+ // Resetting the acting forces
+ bd->resetForceAndTorque();
+
+ // Checking the state of the capsule
+ WALBERLA_ASSERT( bd->checkInvariants(), "Invalid capsule state detected" );
+ }
+ }
+ if (tt_!=NULL) tt_->stop("Integrate Bodies");
+ if (tt_!=NULL) tt_->stop("PlainIntegrator");
+}
+
+} // namespace cr
+} // namespace pe
+} // namespace walberla