Skip to content
Snippets Groups Projects
Commit 2150f2b7 authored by Michael Kuron's avatar Michael Kuron :mortar_board:
Browse files

Switchable integrator and contact resolution for DEM

parent 5af6b820
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/pe/cr/ContactResolvers.h 0 → 100644
+ 59
− 0
View file @ 2150f2b7
#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
+ 31
− 10
View file @ 2150f2b7
......@@ -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
+ 160
− 0
View file @ 2150f2b7
......@@ -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 0 → 100644
+ 108
− 0
View file @ 2150f2b7
#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
+ 25
− 8
View file @ 2150f2b7
......@@ -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
+ 93
− 0
View file @ 2150f2b7
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment