Commit c6267218 authored by Tobias Schruff's avatar Tobias Schruff Committed by Sebastian Eibl
Browse files

Fixed small bug in pe's HCSITS solver and added overrelaxation parameter

parent 8c0a300e
......@@ -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.
*
......
......@@ -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] ) ) ) );
......
Markdown is supported
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