Merge branch 'dem' into 'master'

Switchable integrator and contact resolution for DEM

See merge request !32

src/pe/cr/ContactResolvers.h

+#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

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"

src/pe/cr/DEM.cpp → src/pe/cr/DEM.impl.h

@@ -41,13 +41,17 @@ namespace walberla {
 namespace pe {
 namespace cr {
 
-DEM::DEM(    const shared_ptr<BodyStorage>&    globalBodyStorage
+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)
-   : globalBodyStorage_(globalBodyStorage)
+   : integrate_(integrate), resolveContact_(resolveContact)
+   ,globalBodyStorage_(globalBodyStorage)
    , blockStorage_(blockStorage)
    , storageID_(storageID)
    , ccdID_(ccdID)
@@ -60,7 +64,8 @@ DEM::DEM(    const shared_ptr<BodyStorage>&    globalBodyStorage
 
 }
 
-void DEM::timestep( real_t dt )
+template< typename Integrator, typename ContactResolver >
+void DEMSolver<Integrator,ContactResolver>::timestep( real_t dt )
 {
    maxPenetration_          = real_c(0.0);
    numberOfContacts_        = 0;
@@ -83,7 +88,7 @@ void DEM::timestep( real_t dt )
          if (shouldContactBeTreated( &(*cIt), currentBlock.getAABB() ))
          {
             ++numberOfContactsTreated_;
-            resolveContact( &(*cIt) );
+            resolveContact_( &(*cIt) );
          }
       }
 
@@ -119,7 +124,21 @@ void DEM::timestep( real_t dt )
       //      bodyIt->resetNode();
          bodyIt->clearContacts();
 
-         move( *bodyIt, dt );
+         // 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();
@@ -136,147 +155,6 @@ void DEM::timestep( real_t dt )
    }
 }
 
-//=================================================================================================
-//
-//  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

src/pe/cr/Integrators.h

+#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

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"

src/pe/cr/PlainIntegrator.cpp → src/pe/cr/PlainIntegrator.impl.h

@@ -39,11 +39,14 @@ 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)
+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)
@@ -51,7 +54,8 @@ PlainIntegrator::PlainIntegrator(const shared_ptr<BodyStorage>&      globalBodyS
 
 }
 
-void PlainIntegrator::timestep( const real_t dt )
+template< typename Integrator >
+void PlainIntegratorSolver<Integrator>::timestep( const real_t dt )
 {
    if (tt_!=NULL) tt_->start("PlainIntegrator");
    if (tt_!=NULL) tt_->start("Integrate Bodies");
@@ -60,79 +64,29 @@ void PlainIntegrator::timestep( const real_t dt )
       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();
+         // Checking the state of the body
+         WALBERLA_ASSERT( bd->checkInvariants(), "Invalid capsule state detected" );
+         WALBERLA_ASSERT( !bd->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() );
+         // Resetting the contact node and removing all attached contacts
+      //      bd->resetNode();
+         bd->clearContacts();
 
-      // Calculating the angular acceleration by the equation
-      //   R * Iinv * R^T * torque
-      const Vec3 wdot( id->getInvInertia() * id->getTorque() );
+         // Moving the capsule according to the acting forces (don't move a sleeping capsule)
+         if( bd->isAwake() && !bd->hasInfiniteMass() ) {
+            integrate_( *bd, dt, *this );
+         }
 
-      // Updating the linear velocity
-      id->setLinearVel(id->getLinearVel() + vdot * dt);
+         // Resetting the acting forces
+         bd->resetForceAndTorque();
 
-      // 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();
+         // Checking the state of the capsule
+         WALBERLA_ASSERT( bd->checkInvariants(), "Invalid capsule state detected" );
+      }
    }
-
-   // Resetting the acting forces
-   id->resetForceAndTorque();
-
-   // Checking the state of the capsule
-   WALBERLA_ASSERT( id->checkInvariants(), "Invalid capsule state detected" );
+   if (tt_!=NULL) tt_->stop("Integrate Bodies");
+   if (tt_!=NULL) tt_->stop("PlainIntegrator");
 }
-//*************************************************************************************************
 
 } // namespace cr
 } // namespace pe