diff --git a/src/pe/cr/HCSITS.h b/src/pe/cr/HCSITS.h
index 896f358080c9aed186ecf1d203b503f3285b297c..67056df71da2ccc3b15a3085517e0529c951d502 100644
--- a/src/pe/cr/HCSITS.h
+++ b/src/pe/cr/HCSITS.h
@@ -84,12 +84,12 @@ private:
{
void resize(const size_t n);
- std::vector<Vec3> r1_, r2_;
+ std::vector<Vec3> r1_, r2_; // vector pointing from body1/body2 to the contact position
std::vector<BodyID> body1_, body2_;
- std::vector<Vec3> n_, t_, o_;
- std::vector<real_t> dist_;
- std::vector<real_t> mu_;
- std::vector<Mat3> diag_nto_;
+ std::vector<Vec3> n_, t_, o_; // contact normal and the two other directions perpendicular to the normal
+ std::vector<real_t> dist_; // overlap length, a contact is present if dist_ < 0
+ std::vector<real_t> mu_; // contact friction
+ std::vector<Mat3> diag_nto_;
std::vector<Mat3> diag_nto_inv_;
std::vector<Mat2> diag_to_inv_;
std::vector<real_t> diag_n_inv_;
@@ -102,9 +102,9 @@ public:
enum RelaxationModel {
InelasticFrictionlessContact,
ApproximateInelasticCoulombContactByDecoupling,
-// ApproximateInelasticCoulombContactByOrthogonalProjections,
+ ApproximateInelasticCoulombContactByOrthogonalProjections,
InelasticCoulombContactByDecoupling,
-// InelasticCoulombContactByOrthogonalProjections,
+ InelasticCoulombContactByOrthogonalProjections,
InelasticGeneralizedMaximumDissipationContact
};
//**********************************************************************************************
@@ -136,6 +136,7 @@ public:
inline const std::map<IBlockID::IDType, ContactCache> getContactCache() const { return blockToContactCache_; }
inline real_t getSpeedLimitFactor() const;
inline size_t getMaxIterations() const { return maxIterations_; }
+ inline real_t getOverRelaxationParameter() const { return overRelaxationParam_; }
inline real_t getRelaxationParameter() const { return relaxationParam_; }
inline real_t getErrorReductionParameter() const { return erp_; }
inline RelaxationModel getRelaxationModel() const { return relaxationModel_; }
@@ -145,6 +146,7 @@ public:
//**Set functions*******************************************************************************
/*!\name Set functions */
//@{
+ inline void setOverRelaxationParameter( real_t omega );
inline void setRelaxationParameter( real_t f );
inline void setMaxIterations( size_t n );
inline void setRelaxationModel( RelaxationModel relaxationModel );
@@ -235,6 +237,7 @@ private:
size_t maxSubIterations_; //!< Maximum number of iterations of iterative solvers in the one-contact problem.
real_t abortThreshold_; //!< If L-infinity iterate difference drops below this threshold the iteration is aborted.
RelaxationModel relaxationModel_; //!< The method used to relax unilateral contacts
+ real_t overRelaxationParam_; //!< Parameter specifying the convergence speed for othogonal projection models.
real_t relaxationParam_; //!< Parameter specifying underrelaxation of velocity corrections for boundary bodies.
real_t maximumPenetration_;
size_t numContacts_;
@@ -322,6 +325,28 @@ inline real_t HardContactSemiImplicitTimesteppingSolvers::getSpeedLimitFactor()
//
//=================================================================================================
+//*************************************************************************************************
+/*!\brief Sets the relaxation parameter for boundary bodies.
+ *
+ * \param f The overrelaxation parameter.
+ * \return void
+ *
+ * The overrelaxation parameter \omega is only used when the relaxation model is one of
+ * - ApproximateInelasticCoulombContactByOrthogonalProjections
+ * - InelasticCoulombContactByOrthogonalProjections
+ *
+ * It is used to control the convergence of the model. Large values show faster convergence,
+ * but they can also lead to divergence ("exploding" particles). The default values is 1.0.
+ */
+inline void HardContactSemiImplicitTimesteppingSolvers::setOverRelaxationParameter( real_t omega )
+{
+ WALBERLA_ASSERT_GREATER( omega, 0, "Overrelaxation parameter must be positive." );
+
+ overRelaxationParam_ = omega;
+}
+//*************************************************************************************************
+
+
//*************************************************************************************************
/*!\brief Sets the relaxation parameter for boundary bodies.
*
diff --git a/src/pe/cr/HCSITS.impl.h b/src/pe/cr/HCSITS.impl.h
index ba8fd6aa29d0ec028448b838838406c9d49c50f2..7218b4580519cd623a74976c678a5786eba69184 100644
--- a/src/pe/cr/HCSITS.impl.h
+++ b/src/pe/cr/HCSITS.impl.h
@@ -97,6 +97,7 @@ inline HardContactSemiImplicitTimesteppingSolvers::HardContactSemiImplicitTimest
, maxSubIterations_ ( 20 )
, abortThreshold_ ( real_c(1e-7) )
, relaxationModel_ ( InelasticFrictionlessContact )
+ , overRelaxationParam_( real_c(1.0) )
, relaxationParam_ ( real_c(0.75) )
, maximumPenetration_ ( real_c(0.0) )
, numContacts_ ( 0 )
@@ -388,17 +389,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 );
@@ -1049,55 +1050,68 @@ inline real_t HardContactSemiImplicitTimesteppingSolvers::relaxInelasticCoulombC
bodyCache.dv_[contactCache.body2_[i]->index_] += contactCache.body2_[i]->getInvMass() * contactCache.p_[i];
bodyCache.dw_[contactCache.body2_[i]->index_] += contactCache.body2_[i]->getInvInertia() * ( contactCache.r2_[i] % contactCache.p_[i] );
- // Calculate the relative contact velocity in the global world frame (if no contact reaction is present at contact i)
- Vec3 gdot ( ( bodyCache.v_[contactCache.body1_[i]->index_] + bodyCache.dv_[contactCache.body1_[i]->index_] ) - ( bodyCache.v_[contactCache.body2_[i]->index_] + bodyCache.dv_[contactCache.body2_[i]->index_] ) + ( bodyCache.w_[contactCache.body1_[i]->index_] + bodyCache.dw_[contactCache.body1_[i]->index_] ) % contactCache.r1_[i] - ( bodyCache.w_[contactCache.body2_[i]->index_] + bodyCache.dw_[contactCache.body2_[i]->index_] ) % contactCache.r2_[i] /* + diag_[i] * p */ );
+ // Calculate the relative contact velocity in the global world frame (as if no contact reaction was present at contact i)
+ Vec3 gdot( ( bodyCache.v_[contactCache.body1_[i]->index_] + bodyCache.dv_[contactCache.body1_[i]->index_] )
+ - ( bodyCache.v_[contactCache.body2_[i]->index_] + bodyCache.dv_[contactCache.body2_[i]->index_] )
+ + ( bodyCache.w_[contactCache.body1_[i]->index_] + bodyCache.dw_[contactCache.body1_[i]->index_] ) % contactCache.r1_[i]
+ - ( bodyCache.w_[contactCache.body2_[i]->index_] + bodyCache.dw_[contactCache.body2_[i]->index_] ) % contactCache.r2_[i] /* + diag_[i] * p */ );
// Change from the global world frame to the contact frame
Mat3 contactframe( contactCache.n_[i], contactCache.t_[i], contactCache.o_[i] );
Vec3 gdot_nto( contactframe.getTranspose() * gdot );
- //real_t gdot_n ( trans( contactCache.n_[i] ) * gdot ); // The component of gdot along the contact normal n
- //Vec3 gdot_t ( gdot - gdot_n * contactCache.n_[i] ); // The components of gdot tangential to the contact normal n
- //real_t g_n ( gdot_n * dt /* + trans( contactCache.n_[i] ) * ( contactCache.body1_[i]->getPosition() + contactCache.r1_[i] ) - ( contactCache.body2_[i]->getPosition() + contactCache.r2_[i] ) */ + contactCache.dist_[i] ); // The gap in normal direction
+ //real_t gdot_n( trans( contactCache.n_[i] ) * gdot ); // The component of gdot along the contact normal n
+ //Vec3 gdot_t( gdot - gdot_n * contactCache.n_[i] ); // The components of gdot tangential to the contact normal n
+ //real_t g_n( gdot_n * dt /* + trans( contactCache.n_[i] ) * ( contactCache.body1_[i]->getPosition() + contactCache.r1_[i] ) - ( contactCache.body2_[i]->getPosition() + contactCache.r2_[i] ) */ + contactCache.dist_[i] ); // The gap in normal direction
// The constraint in normal direction is actually a positional constraint but instead of g_n we use g_n/dt equivalently and call it gdot_n
- gdot_nto[0] += ( /* + trans( contactCache.n_[i] ) * ( contactCache.body1_[i]->getPosition() + contactCache.r1_[i] ) - ( contactCache.body2_[i]->getPosition() + contactCache.r2_[i] ) */ + contactCache.dist_[i] ) * dtinv;
+ gdot_nto[0] += ( /* + trans( contactCache.n_[i] ) * ( contactCache.body1_[i]->getPosition() + contactCache.r1_[i] )
+ - ( contactCache.body2_[i]->getPosition() + contactCache.r2_[i] ) */ + contactCache.dist_[i] ) * dtinv;
+
+ // Tasora et al., 2011 starts here
+ // -------------------------------------------------------------------------------------
- const real_t w( 1 ); // w > 0
Vec3 p_cf( contactframe.getTranspose() * contactCache.p_[i] );
if( approximate ) {
- // Calculate next iterate (Anitescu/Tasora).
- p_cf = p_cf - w * ( contactCache.diag_nto_[i] * p_cf + gdot_nto );
+ // Calculate next iterate (Tasora et al., 2011).
+ p_cf = p_cf - overRelaxationParam_ * ( contactCache.diag_nto_[i] * p_cf + gdot_nto );
}
else {
// Calculate next iterate (De Saxce/Feng).
Vec3 tmp( contactCache.diag_nto_[i] * p_cf + gdot_nto );
tmp[0] += std::sqrt( math::sq( tmp[1] ) + math::sq( tmp[2] ) ) * contactCache.mu_[i];
- p_cf = p_cf - w * tmp;
+ p_cf = p_cf - overRelaxationParam_ * tmp;
}
- // Project.
+ // Project on friction cone (Tasora et al., 2011; Fig. 2; Eq. (46)).
real_t flimit( contactCache.mu_[i] * p_cf[0] );
real_t fsq( p_cf[1] * p_cf[1] + p_cf[2] * p_cf[2] );
+
if( p_cf[0] > 0 && fsq < flimit * flimit ) {
- // Unconstrained minimum is in cone leading to a static contact and no projection
- // is necessary.
+ // Unconstrained minimum is INSIDE THE UPPER CONE,
+ // with slope n*mu
+ // leading to a static contact and no projection is necessary.
}
- else if( p_cf[0] < 0 && fsq < p_cf[0] / math::sq( contactCache.mu_[i] ) ) {
- // Unconstrained minimum is in dual cone leading to a separating contact where no contact
- // reaction is present (the unconstrained minimum is projected to the tip of the cone).
-
+ else if( p_cf[0] < 0 && fsq < math::sq( p_cf[0] / contactCache.mu_[i] ) ) {
+ // Unconstrained minimum is INSIDE THE LOWER CONE,
+ // with slope n/mu
+ // leading to a separating contact where no contact reaction is present,
+ // i.e. the unconstrained minimum is projected to the tip of the lower cone.
p_cf = Vec3();
}
else {
- // The contact is dynamic.
- real_t f( std::sqrt( fsq ) );
- p_cf[0] = ( f * contactCache.mu_[i] + p_cf[0] ) / ( math::sq( contactCache.mu_[i] ) + 1 );
+ // Unconstrained minimum is OUTSIDE THE CONE -> Project on cone surface, Eq. (45) in Tasora et al., 2011
+ real_t f( std::sqrt( fsq ) ); // gamma_r
+ p_cf[0] = ( f * contactCache.mu_[i] + p_cf[0] ) / ( math::sq( contactCache.mu_[i] ) + 1 ); // p_cf[0] = gamma_n
+
real_t factor( contactCache.mu_[i] * p_cf[0] / f );
p_cf[1] *= factor;
p_cf[2] *= factor;
}
+ // Tasora et al., 2011 ends here
+ // -------------------------------------------------------------------------------------
+
Vec3 p_wf( contactframe * p_cf );
Vec3 dp( contactCache.p_[i] - p_wf );
delta_max = std::max( delta_max, std::max( std::abs( dp[0] ), std::max( std::abs( dp[1] ), std::abs( dp[2] ) ) ) );