From 6c5f354eaf58758249c1d19cf106581edddf375a Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Thu, 18 May 2017 14:18:13 +0200
Subject: [PATCH] added speed limiter to HCSITS
---
src/pe/cr/HCSITS.h | 43 +++++++++++++++++++++++++++++++++++++++
src/pe/cr/HCSITS.impl.h | 45 ++++++++++++++++++++++++++---------------
tests/pe/HCSITS.cpp | 23 +++++++++++++++++++++
3 files changed, 95 insertions(+), 16 deletions(-)
diff --git a/src/pe/cr/HCSITS.h b/src/pe/cr/HCSITS.h
index 41f863e37..0af5c6fa9 100644
--- a/src/pe/cr/HCSITS.h
+++ b/src/pe/cr/HCSITS.h
@@ -134,6 +134,7 @@ public:
virtual inline size_t getNumberOfContacts() const WALBERLA_OVERRIDE;
virtual inline size_t getNumberOfContactsTreated() const WALBERLA_OVERRIDE;
inline const std::map<IBlockID::IDType, ContactCache> getContactCache() const { return blockToContactCache_; }
+ inline real_t getSpeedLimitFactor() const;
//@}
//**********************************************************************************************
@@ -145,6 +146,7 @@ public:
inline void setRelaxationModel( RelaxationModel relaxationModel );
inline void setErrorReductionParameter( real_t erp );
inline void setAbortThreshold( real_t threshold );
+ inline void setSpeedLimiter( bool active, const real_t speedLimitFactor = real_t(0.0) );
//@}
//**********************************************************************************************
@@ -153,6 +155,7 @@ public:
//@{
inline bool isSyncRequired() const;
inline bool isSyncRequiredLocally() const;
+ inline bool isSpeedLimiterActive() const;
//@}
//**********************************************************************************************
@@ -235,6 +238,9 @@ private:
size_t numContacts_;
size_t numContactsTreated_;
+ bool speedLimiterActive_; //!< is the speed limiter active?
+ real_t speedLimitFactor_; //!< what multiple of boundingbox edge length is the body allowed to travel in one timestep
+
//**********************************************************************************************
/*! \cond WALBERLA_INTERNAL */
/*!\brief Functor for comparing the system ID of two bodies.
@@ -299,6 +305,11 @@ inline size_t HardContactSemiImplicitTimesteppingSolvers::getNumberOfContactsTre
{
return numContactsTreated_;
}
+
+inline real_t HardContactSemiImplicitTimesteppingSolvers::getSpeedLimitFactor() const
+{
+ return speedLimitFactor_;
+}
//*************************************************************************************************
@@ -393,6 +404,23 @@ inline void HardContactSemiImplicitTimesteppingSolvers::setAbortThreshold( real_
+//*************************************************************************************************
+/*!\brief Activates/Deactivates the speed limiter and sets the limit
+*
+* \param active activate/deactivate speed limtier
+* \param speedLimitFactor size of bounding box will be multiplied by this factor to get the maximal distance a body is allowed to travel within one timestep
+* \return void
+*/
+inline void HardContactSemiImplicitTimesteppingSolvers::setSpeedLimiter( bool active, const real_t speedLimitFactor )
+{
+ speedLimiterActive_ = active;
+ speedLimitFactor_ = speedLimitFactor;
+}
+//*************************************************************************************************
+
+
+
+
//=================================================================================================
//
// QUERY FUNCTIONS
@@ -453,6 +481,21 @@ inline bool HardContactSemiImplicitTimesteppingSolvers::isSyncRequiredLocally()
}
//*************************************************************************************************
+
+
+
+//*************************************************************************************************
+/*!\brief Returns if speed limiter is currently active and working.
+ *
+ * \return status of the speed limiter
+ */
+inline bool HardContactSemiImplicitTimesteppingSolvers::isSpeedLimiterActive() const
+{
+ return speedLimiterActive_;
+}
+//*************************************************************************************************
+
+
} // namespace cr
} // namespace pe
diff --git a/src/pe/cr/HCSITS.impl.h b/src/pe/cr/HCSITS.impl.h
index d2c8b9ba7..8be3ec9fd 100644
--- a/src/pe/cr/HCSITS.impl.h
+++ b/src/pe/cr/HCSITS.impl.h
@@ -102,6 +102,8 @@ inline HardContactSemiImplicitTimesteppingSolvers::HardContactSemiImplicitTimest
, maximumPenetration_ ( real_c(0.0) )
, numContacts_ ( 0 )
, numContactsTreated_( 0)
+ , speedLimiterActive_( false )
+ , speedLimitFactor_ ( real_c(1.0) )
, requireSync_ ( false )
{
// Logging the successful setup of the collision system
@@ -364,17 +366,17 @@ inline void HardContactSemiImplicitTimesteppingSolvers::timestep( const real_t d
delta_max = relaxApproximateInelasticCoulombContactsByDecoupling( dtinv, contactCache, bodyCache );
break;
-// case ApproximateInelasticCoulombContactByOrthogonalProjections:
-// delta_max = relaxInelasticCoulombContactsByOrthogonalProjections( dtinv, true, contactCache, bodyCache );
-// break;
+ // case ApproximateInelasticCoulombContactByOrthogonalProjections:
+ // delta_max = relaxInelasticCoulombContactsByOrthogonalProjections( dtinv, true, contactCache, bodyCache );
+ // break;
case InelasticCoulombContactByDecoupling:
delta_max = relaxInelasticCoulombContactsByDecoupling( dtinv, contactCache, bodyCache );
break;
-// case InelasticCoulombContactByOrthogonalProjections:
-// delta_max = relaxInelasticCoulombContactsByOrthogonalProjections( dtinv, false, contactCache, bodyCache );
-// break;
+ // case InelasticCoulombContactByOrthogonalProjections:
+ // delta_max = relaxInelasticCoulombContactsByOrthogonalProjections( dtinv, false, contactCache, bodyCache );
+ // break;
case InelasticGeneralizedMaximumDissipationContact:
delta_max = relaxInelasticGeneralizedMaximumDissipationContacts( dtinv, contactCache, bodyCache );
@@ -455,10 +457,10 @@ inline void HardContactSemiImplicitTimesteppingSolvers::timestep( const real_t d
BodyStorage& shadowStorage = (*storage)[1];
for( auto body = ConcatIterator<BodyStorage::Iterator>
- (localStorage.begin(),
- localStorage.end(),
- shadowStorage.begin(),
- shadowStorage.end()); body != ConcatIterator<BodyStorage::Iterator>(); ++body )
+ (localStorage.begin(),
+ localStorage.end(),
+ shadowStorage.begin(),
+ shadowStorage.end()); body != ConcatIterator<BodyStorage::Iterator>(); ++body )
{
if (!body->isCommunicating())
{
@@ -1665,7 +1667,7 @@ inline void HardContactSemiImplicitTimesteppingSolvers::synchronizeVelocities( )
if (tt_ != NULL) tt_->stop("Velocity Sync Update Processing");
if (tt_ != NULL) tt_->start("Velocity Sync Globals");
-/*
+ /*
{
size_t i;
std::vector<real_t> reductionBuffer( globalBodyStorage_->size() * 6, real_c(0) );
@@ -1775,6 +1777,17 @@ inline void HardContactSemiImplicitTimesteppingSolvers::integratePositions( Body
if( body->isAwake() )
{
+ if ( isSpeedLimiterActive() )
+ {
+ const auto speed = v.length();
+ const auto aabb = body->getAABB();
+ const auto edge = std::min(aabb.xSize(), std::min(aabb.ySize(), aabb.zSize()));
+ if (speed * dt > edge * getSpeedLimitFactor() )
+ {
+ v = v * (edge * getSpeedLimitFactor() / dt / speed );
+ }
+ }
+
// Calculating the translational displacement
body->setPosition( body->getPosition() + v * dt );
@@ -1821,19 +1834,19 @@ void configure( const Config::BlockHandle& config, pe::cr::HCSITS& cr )
cr::HCSITS::RelaxationModel HCSITSRelaxationModel;
if (HCSITSRelaxationModelStr == "InelasticFrictionlessContact")
{
- HCSITSRelaxationModel = cr::HCSITS::InelasticFrictionlessContact;
+ HCSITSRelaxationModel = cr::HCSITS::InelasticFrictionlessContact;
} else if (HCSITSRelaxationModelStr == "ApproximateInelasticCoulombContactByDecoupling")
{
- HCSITSRelaxationModel = cr::HCSITS::ApproximateInelasticCoulombContactByDecoupling;
+ HCSITSRelaxationModel = cr::HCSITS::ApproximateInelasticCoulombContactByDecoupling;
} else if (HCSITSRelaxationModelStr == "InelasticCoulombContactByDecoupling")
{
- HCSITSRelaxationModel = cr::HCSITS::InelasticCoulombContactByDecoupling;
+ HCSITSRelaxationModel = cr::HCSITS::InelasticCoulombContactByDecoupling;
} else if (HCSITSRelaxationModelStr == "InelasticGeneralizedMaximumDissipationContact")
{
- HCSITSRelaxationModel = cr::HCSITS::InelasticGeneralizedMaximumDissipationContact;
+ HCSITSRelaxationModel = cr::HCSITS::InelasticGeneralizedMaximumDissipationContact;
} else
{
- WALBERLA_ABORT("Unknown HCSITSRelaxationModel: " << HCSITSRelaxationModelStr);
+ WALBERLA_ABORT("Unknown HCSITSRelaxationModel: " << HCSITSRelaxationModelStr);
}
Vec3 globalLinearAcceleration = config.getParameter<Vec3>("globalLinearAcceleration", Vec3(0, 0, 0));
diff --git a/tests/pe/HCSITS.cpp b/tests/pe/HCSITS.cpp
index ae0472195..93faeca36 100644
--- a/tests/pe/HCSITS.cpp
+++ b/tests/pe/HCSITS.cpp
@@ -77,6 +77,25 @@ void normalReactionTest(cr::HCSITS& cr, SphereID sp)
WALBERLA_CHECK_FLOAT_EQUAL( sp->getLinearVel(), Vec3(0,0,0) );
}
+void speedLimiterTest(cr::HCSITS& cr, SphereID sp)
+{
+ cr.setErrorReductionParameter( real_t(1.0) );
+
+ sp->setPosition( Vec3(5,5,6) );
+ sp->setLinearVel( Vec3(0,0,0) );
+ cr.setSpeedLimiter( true, real_t(0.2) );
+ cr.timestep( real_c( real_t(1.0) ) );
+ WALBERLA_CHECK_FLOAT_EQUAL( sp->getPosition() , Vec3(5,5,real_t(6.1)) );
+ WALBERLA_CHECK_FLOAT_EQUAL( sp->getLinearVel(), Vec3(0,0,real_t(0.1)) );
+
+ sp->setPosition( Vec3(5,5,5.5) );
+ sp->setLinearVel( Vec3(0,0,0) );
+ cr.setSpeedLimiter( true, real_t(0.2) );
+ cr.timestep( real_c( real_t(1.0) ) );
+ WALBERLA_CHECK_FLOAT_EQUAL( sp->getPosition() , Vec3(5,5,real_t(5.94)) );
+ WALBERLA_CHECK_FLOAT_EQUAL( sp->getLinearVel(), Vec3(0,0,real_t(0.44)) );
+}
+
int main( int argc, char** argv )
{
walberla::debug::enterTestMode();
@@ -135,5 +154,9 @@ int main( int argc, char** argv )
cr.setRelaxationModel( cr::HardContactSemiImplicitTimesteppingSolvers::InelasticGeneralizedMaximumDissipationContact );
normalReactionTest(cr, sp);
+ WALBERLA_LOG_PROGRESS("SpeedLimiter Test: InelasticFrictionlessContact");
+ cr.setRelaxationModel( cr::HardContactSemiImplicitTimesteppingSolvers::InelasticFrictionlessContact );
+ speedLimiterTest(cr, sp);
+
return EXIT_SUCCESS;
}
--
GitLab