diff --git a/src/pe/Materials.cpp b/src/pe/Materials.cpp
index 1bdc0326bb166b0590b6b33c39efc85a18ca2a1a..e7619d68e9fca921dd780d2059b79c5d96cbd39b 100644
--- a/src/pe/Materials.cpp
+++ b/src/pe/Materials.cpp
@@ -160,7 +160,8 @@ bool Material::activateMaterials()
*/
MaterialID createMaterial( const std::string& name, real_t density, real_t cor,
real_t csf, real_t cdf, real_t poisson, real_t young,
- real_t stiffness, real_t dampingN, real_t dampingT, real_t surfEnergy )
+ real_t stiffness, real_t dampingN, real_t dampingT,
+ real_t surfEnergy, real_t yieldStrength )
{
using M = Material;
@@ -207,8 +208,12 @@ MaterialID createMaterial( const std::string& name, real_t density, real_t cor,
if( surfEnergy < real_c(0) )
throw std::invalid_argument( "Invalid surface energy" );
+ if( yieldStrength < real_c(0) )
+ throw std::invalid_argument( "Invalid yield strength" );
+
// Registering the new material
- M::materials_.push_back( Material( name, density, cor, csf, cdf, poisson, young, stiffness, dampingN, dampingT, surfEnergy ) );
+ M::materials_.push_back( Material( name, density, cor, csf, cdf, poisson, young, stiffness, dampingN, dampingT,
+ surfEnergy, yieldStrength ) );
const MaterialID mat( M::materials_.size()-1 );
// Updating the restitution table, the static friction table and the dynamic friction table
@@ -255,7 +260,7 @@ MaterialID createMaterial( const std::string& name, real_t density, real_t cor,
// Note that the material has to be created on all processes in an MPI parallel simulation.
*/
MaterialID createMaterial( real_t density, real_t cor, real_t csf, real_t cdf, real_t poisson, real_t young,
- real_t stiffness, real_t dampingN, real_t dampingT, real_t surfEnergy )
+ real_t stiffness, real_t dampingN, real_t dampingT, real_t surfEnergy, real_t yieldStrength )
{
std::ostringstream sstr;
@@ -266,7 +271,8 @@ MaterialID createMaterial( real_t density, real_t cor, real_t csf, real_t cdf, r
sstr << "Material" << ++Material::anonymousMaterials_;
} while( Material::find( sstr.str() ) != invalid_material );
- return createMaterial( sstr.str(), density, cor, csf, cdf, poisson, young, stiffness, dampingN, dampingT, surfEnergy );
+ return createMaterial( sstr.str(), density, cor, csf, cdf, poisson, young, stiffness, dampingN, dampingT,
+ surfEnergy, yieldStrength );
}
//*************************************************************************************************
@@ -327,9 +333,10 @@ std::string Material::toString( const MaterialID& v )
<< " stiffness = " << Material::getStiffness(v) << "\n"
<< " dampingN = " << Material::getDampingN(v) << "\n"
<< " dampingT = " << Material::getDampingT(v) << "\n"
- << " surface energy = " << Material::getSurfaceEnergy(v);
+ << " surface energy = " << Material::getSurfaceEnergy(v) << "\n"
+ << " yield strength = " << Material::getYieldStrength(v);
return ss.str();
}
} // namespace pe
-} // namespace walberla
+} // namespace walberla
diff --git a/src/pe/Materials.h b/src/pe/Materials.h
index 59a46e8f5f56353b08e1c93d2ca4ce8af4984de1..06e36abcfa8a886ae51db1d07414c7ea46274aa8 100644
--- a/src/pe/Materials.h
+++ b/src/pe/Materials.h
@@ -106,7 +106,8 @@ public:
//@{
explicit inline Material( const std::string& name, real_t density, real_t cor,
real_t csf, real_t cdf, real_t poisson, real_t young, real_t stiffness,
- real_t dampingN, real_t dampingT, real_t surfEnergy = real_t(0) );
+ real_t dampingN, real_t dampingT,
+ real_t surfEnergy = real_t(0), real_t yieldStrength = real_t(1E8) );
// No explicitly declared copy constructor.
//@}
//**********************************************************************************************
@@ -133,6 +134,7 @@ public:
inline real_t getDampingN() const;
inline real_t getDampingT() const;
inline real_t getSurfaceEnergy() const;
+ inline real_t getYieldStrength() const;
static MaterialID find( const std::string& name );
static std::vector<MaterialID> findPrefix( const std::string& prefix );
@@ -154,6 +156,7 @@ public:
static inline real_t getDampingT( MaterialID material );
static inline real_t getDampingT( MaterialID material1, MaterialID material2 );
static inline real_t getSurfaceEnergy( MaterialID material );
+ static inline real_t getYieldStrength( MaterialID material );
static std::string toString( const MaterialID& v );
//@}
//**********************************************************************************************
@@ -236,7 +239,8 @@ private:
/*!< Friction counteracts the tangential relative velocity and thus can be
modelled as a viscous damper with a limited damping force. The viscous
damping coefficient corresponds to this parameter.*/
- real_t surf_; //!< The surface energy at the contact region in mJ/m².
+ real_t surf_; //!< The surface energy at the contact region in \f$ mJ/m^2 \f$.
+ real_t yieldStrength_; //!< The elastic limit or yield strength of the material in \f$ N/m^2 \f$.
static Materials materials_; //!< Vector for the registered materials.
static MatN corTable_; //!< Table for the coefficients of restitution.
@@ -251,10 +255,12 @@ private:
/*! \cond internal */
friend MaterialID createMaterial( const std::string& name, real_t density, real_t cor,
real_t csf, real_t cdf, real_t poisson, real_t young,
- real_t stiffness, real_t dampingN, real_t dampingT, real_t surfEnergy = real_t(1) );
+ real_t stiffness, real_t dampingN, real_t dampingT,
+ real_t surfEnergy = real_t(1), real_t yieldStrength = real_t(1E8) );
friend MaterialID createMaterial( real_t density, real_t cor, real_t csf, real_t cdf,
real_t poisson, real_t young,
- real_t stiffness, real_t dampingN, real_t dampingT, real_t surfEnergy = real_t(1) );
+ real_t stiffness, real_t dampingN, real_t dampingT,
+ real_t surfEnergy = real_t(1), real_t yieldStrength = real_t(1E8) );
/*! \endcond */
//**********************************************************************************************
};
@@ -274,10 +280,12 @@ private:
* \param stiffness The stiffness in normal direction of the material's contact region.
* \param dampingN The damping coefficient in normal direction of the material's contact region.
* \param dampingT The damping coefficient in tangential direction of the material's contact region.
+ * \param surfEnergy The surface energy at the contact region in \f$ mJ/m^2 \f$.
+ * \param yieldStrength The yield strength of the material in \f$ N/m^2 \f$.
*/
inline Material::Material( const std::string& name, real_t density, real_t cor,
real_t csf, real_t cdf, real_t poisson, real_t young,
- real_t stiffness, real_t dampingN, real_t dampingT, real_t surfEnergy )
+ real_t stiffness, real_t dampingN, real_t dampingT, real_t surfEnergy, real_t yieldStrength )
: name_ ( name ) // The name of the material
, density_ ( density ) // The density of the material
, restitution_( cor ) // The coefficient of restitution of the material
@@ -289,6 +297,7 @@ inline Material::Material( const std::string& name, real_t density, real_t cor,
, dampingN_ ( dampingN ) // The damping coefficient in normal direction of the material's contact region.
, dampingT_ ( dampingT ) // The damping coefficient in tangential direction of the material's contact region.
, surf_ ( surfEnergy ) // The surface energy at the contact region.
+ , yieldStrength_( yieldStrength ) // The yield strength of the material.
{}
//*************************************************************************************************
@@ -425,6 +434,18 @@ inline real_t Material::getSurfaceEnergy() const
//*************************************************************************************************
+//*************************************************************************************************
+/*!\brief Returns the yield strength of the material in \f$ N/m^2 \f$.
+ *
+ * \return The yield strength of the material.
+ */
+inline real_t Material::getYieldStrength() const
+{
+ return yieldStrength_;
+}
+//*************************************************************************************************
+
+
//*************************************************************************************************
/*!\brief Returns the name of the given material.
*
@@ -726,6 +747,20 @@ inline real_t Material::getSurfaceEnergy( MaterialID material )
//*************************************************************************************************
+//*************************************************************************************************
+/*!\brief Returns the yield strength of the material in \f$ N/m^2 \f$.
+ *
+ * \param material The material to be queried.
+ * \return The yield strength of the given material.
+ */
+inline real_t Material::getYieldStrength( MaterialID material )
+{
+ WALBERLA_ASSERT( material < materials_.size(), "Invalid material ID" );
+ return materials_[material].getYieldStrength();
+}
+//*************************************************************************************************
+
+
//*************************************************************************************************
/*!\brief Setting the coefficient of restitution between material \a material1 and \a material2.
*
@@ -1114,9 +1149,11 @@ const MaterialID invalid_material = static_cast<MaterialID>( -1 );
/*!\name Material functions */
//@{
WALBERLA_PUBLIC MaterialID createMaterial( const std::string& name, real_t density, real_t cor,
- real_t csf, real_t cdf, real_t poisson, real_t young, real_t stiffness, real_t dampingN, real_t dampingT, real_t surfEnergy );
+ real_t csf, real_t cdf, real_t poisson, real_t young, real_t stiffness, real_t dampingN, real_t dampingT,
+ real_t surfEnergy, real_t yieldStrength );
WALBERLA_PUBLIC MaterialID createMaterial( real_t density, real_t cor, real_t csf, real_t cdf,
- real_t poisson, real_t young, real_t stiffness, real_t dampingN, real_t dampingT, real_t surfEnergy );
+ real_t poisson, real_t young, real_t stiffness, real_t dampingN, real_t dampingT,
+ real_t surfEnergy, real_t yieldStrength );
//@}
//*************************************************************************************************
diff --git a/src/pe/contact/ContactFunctions.h b/src/pe/contact/ContactFunctions.h
index 1e328b628504afb0c0bf3b46c0bd6ac9a29b7773..f753b1f2e67eb2237f2ae406bb7a35068bb37441 100644
--- a/src/pe/contact/ContactFunctions.h
+++ b/src/pe/contact/ContactFunctions.h
@@ -37,6 +37,7 @@ inline real_t getDampingT(ConstContactID c);
// inline real_t getCorParam(ConstContactID c);
inline real_t getFriction(ConstContactID c);
inline real_t getSurfaceEnergy(ConstContactID c);
+inline real_t getYieldStrength(ConstContactID c);
} // namespace pe
} // namespace walberla
diff --git a/src/pe/contact/ContactFunctions.impl.h b/src/pe/contact/ContactFunctions.impl.h
index 673946214006c41c70c7c776d03bd55cb4ce7ecd..4641e63681096ecd9dc73b32b0d0a0c07dbc5675 100644
--- a/src/pe/contact/ContactFunctions.impl.h
+++ b/src/pe/contact/ContactFunctions.impl.h
@@ -166,5 +166,10 @@ inline real_t getSurfaceEnergy( ConstContactID c )
return math::max( Material::getSurfaceEnergy( c->getBody1()->getMaterial() ), Material::getSurfaceEnergy( c->getBody2()->getMaterial() ) );
}
+inline real_t getYieldStrength( ConstContactID c )
+{
+ return math::min( Material::getYieldStrength( c->getBody1()->getMaterial() ), Material::getYieldStrength( c->getBody2()->getMaterial() ) );
+}
+
} // namespace pe
} // namespace walberla
diff --git a/src/pe/cr/HCSITS.h b/src/pe/cr/HCSITS.h
index 2099644a1a42764ef398a0a066ba18bda0628c71..e1113588e405452ee95bf26ecd85c2cdd360429c 100644
--- a/src/pe/cr/HCSITS.h
+++ b/src/pe/cr/HCSITS.h
@@ -73,7 +73,7 @@ private:
// bodies
struct BodyCache
{
- void resize(const size_t n);
+ void resize( size_t n );
std::vector<Vec3> v_, w_, dv_, dw_;
};
@@ -82,19 +82,20 @@ private:
// contacts
struct ContactCache
{
- void resize(const size_t n);
+ void resize( size_t n );
std::vector<Vec3> r1_, r2_; // vector pointing from body1/body2 to the contact position
- std::vector<BodyID> body1_, body2_;
+ std::vector<BodyID> body1_, body2_; // the two bodies involved into the contact
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<real_t> surf_; // surface energy
- std::vector<Mat3> diag_nto_;
- std::vector<Mat3> diag_nto_inv_;
- std::vector<Mat2> diag_to_inv_;
- std::vector<real_t> diag_n_inv_;
- std::vector<Vec3> p_;
+ std::vector<real_t> k_; // yield strength
+ std::vector<Mat3> diag_nto_; // transforms impulse (n/t/o) to velocity in the contact frame
+ std::vector<Mat3> diag_nto_inv_; // transforms velocity to impulse (n/t/o) in the contact frame
+ std::vector<Mat2> diag_to_inv_; // transforms velocity to tangential impulses (t/o) in the world frame
+ std::vector<real_t> diag_n_inv_; // transforms velocity to normal impulse (n) in the world frame
+ std::vector<Vec3> p_; // the contact impulse
};
std::map<IBlockID::IDType, ContactCache> blockToContactCache_;
@@ -118,7 +119,7 @@ public:
domain_decomposition::BlockDataID storageID,
domain_decomposition::BlockDataID ccdID,
domain_decomposition::BlockDataID fcdID,
- WcTimingTree* tt = NULL );
+ WcTimingTree* tt = nullptr );
//@}
//**********************************************************************************************
//**Destructor**********************************************************************************
@@ -153,7 +154,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) );
+ inline void setSpeedLimiter( bool active, real_t speedLimitFactor = real_t(0.0) );
//@}
//**********************************************************************************************
@@ -168,9 +169,9 @@ public:
/// forwards to timestep
/// Convenience operator to make class a functor.
- void operator()(const real_t dt) { timestep(dt); }
+ void operator()( real_t dt ) { timestep(dt); }
/// Advances the simulation dt seconds.
- void timestep( const real_t dt );
+ void timestep( real_t dt );
private:
@@ -217,7 +218,7 @@ private:
//**Utility functions***************************************************************************
/*!\name Utility functions */
//@{
- bool checkUpdateFlags();
+ bool checkUpdateFlags();
//@}
//**********************************************************************************************
diff --git a/src/pe/cr/HCSITS.impl.h b/src/pe/cr/HCSITS.impl.h
index f5eb4d6c49191dabc478a661235367bb8647b80c..9388cd22976cf23defece2e4eb3e6502842231d3 100644
--- a/src/pe/cr/HCSITS.impl.h
+++ b/src/pe/cr/HCSITS.impl.h
@@ -30,6 +30,7 @@
#include "pe/ccd/ICCD.h"
#include "pe/fcd/IFCD.h"
#include "pe/contact/ContactFunctions.h"
+#include "pe/cr/YieldSurface.h"
#include "core/math/Constants.h"
#include "core/math/Limits.h"
@@ -43,7 +44,7 @@ namespace walberla {
namespace pe {
namespace cr {
-inline void HardContactSemiImplicitTimesteppingSolvers::BodyCache::resize(const size_t n)
+inline void HardContactSemiImplicitTimesteppingSolvers::BodyCache::resize( const size_t n )
{
v_.resize(n);
w_.resize(n);
@@ -51,7 +52,7 @@ inline void HardContactSemiImplicitTimesteppingSolvers::BodyCache::resize(const
dw_.resize(n);
}
-inline void HardContactSemiImplicitTimesteppingSolvers::ContactCache::resize(const size_t n)
+inline void HardContactSemiImplicitTimesteppingSolvers::ContactCache::resize( const size_t n )
{
r1_.resize(n);
r2_.resize(n);
@@ -63,6 +64,7 @@ inline void HardContactSemiImplicitTimesteppingSolvers::ContactCache::resize(con
dist_.resize(n);
mu_.resize(n);
surf_.resize(n);
+ k_.resize(n);
diag_nto_.resize(n);
diag_nto_inv_.resize(n);
diag_to_inv_.resize(n);
@@ -260,6 +262,7 @@ inline void HardContactSemiImplicitTimesteppingSolvers::timestep( const real_t d
contactCache.mu_[j] = getFriction(c);
contactCache.surf_[j] = getSurfaceEnergy(c);
+ contactCache.k_[i] = getYieldStrength(c);
Mat3 diag = -( math::skewSymCrossProduct(contactCache.r1_[j], math::skewSymCrossProduct(b1->getInvInertia(), contactCache.r1_[j]))
+ math::skewSymCrossProduct(contactCache.r2_[j], math::skewSymCrossProduct(b2->getInvInertia(), contactCache.r2_[j])));
@@ -540,26 +543,25 @@ inline real_t HardContactSemiImplicitTimesteppingSolvers::relaxInelasticFriction
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 );
// 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;
if( gdot_nto[0] >= 0 ) {
- // Contact is separating if no contact reaction is present at contact i.
-
+ // Contact is separating if no contact reaction is present at contact i
delta_max = std::max( delta_max, std::max( std::abs( contactCache.p_[i][0] ), std::max( std::abs( contactCache.p_[i][1] ), std::abs( contactCache.p_[i][2] ) ) ) );
+ // Reset impulse to zero.
contactCache.p_[i] = Vec3();
-
- // No need to apply zero impulse.
}
else {
// Contact is persisting.
@@ -1086,35 +1088,40 @@ inline real_t HardContactSemiImplicitTimesteppingSolvers::relaxInelasticCoulombC
p_cf = p_cf - overRelaxationParam_ * tmp;
}
- // TODO: calculate cohesion from contact area
- const real_t cohesion = contactCache.surf_[i];
- p_cf[0] += cohesion;
+ // Cohesion -->
+ real_t c( 0 );
+ if( contactCache.surf_[i] > real_t(0) )
+ c = contactCache.surf_[i] * contactCache.dist_[i] * contactCache.dist_[i];
+ // <-- Cohesion
- // 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] );
+// // 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 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 < 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( cohesion, 0, 0 );
+// }
+// else {
+// // 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;
+// }
- if( p_cf[0] > 0 && fsq < flimit * flimit ) {
- // 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 < 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( cohesion, 0, 0 );
- }
- else {
- // 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;
- }
+ CoulombOrowan cone( contactCache.mu_[i], contactCache.k_[i], c );
+ cone.project( p_cf );
// Tasora et al., 2011 ends here
// -------------------------------------------------------------------------------------
diff --git a/src/pe/cr/YieldSurface.h b/src/pe/cr/YieldSurface.h
new file mode 100644
index 0000000000000000000000000000000000000000..e9260b56ae885a41fd984c5a01382f042df4a141
--- /dev/null
+++ b/src/pe/cr/YieldSurface.h
@@ -0,0 +1,249 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file HCSITS.h
+//! \author Tobias Preclik
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//! \brief Header file for the hard contact solver
+//
+//======================================================================================================================
+
+#pragma once
+
+//*************************************************************************************************
+// Includes
+//*************************************************************************************************
+
+#include "pe/Types.h"
+
+namespace walberla {
+namespace pe {
+namespace cr {
+
+
+//=================================================================================================
+//
+// CLASS DEFINITIONS
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/**
+ * \ingroup pe
+ * \brief Abstract base class for yield surfaces.
+ */
+class YieldSurface
+{
+ /**
+ * \brief Checks if the contact reaction is inside the yield surface.
+ * \param p The contact reaction in contact frame coordinates.
+ * \return A bool indicating whether the contact reaction is inside the yield surface.
+ */
+ virtual bool contains( const Vec3 & p ) = 0;
+
+ /**
+ * \brief Handles the projection of the contact reaction on the yield surface.
+ * \param p The contact reaction in contact frame coordinates.
+ */
+ virtual void project( Vec3 & p ) = 0;
+};
+
+
+/**
+ * \ingroup pe
+ * \brief The Coulomb friction cone with optional cohesion.
+ */
+class Coulomb : public YieldSurface
+{
+public:
+ explicit Coulomb( real_t mu, real_t c = real_t(0) )
+ : mu_( mu ), c_( c ) { }
+
+ bool contains( const Vec3 & p ) WALBERLA_OVERRIDE
+ {
+ // Unconstrained minimum is INSIDE THE UPPER CONE,
+ // with slope n*mu and minimum -c.
+ return ( p[0] > -c_ && p[1] < mu_ * p[0] );
+ }
+
+ void project( Vec3 & p ) WALBERLA_OVERRIDE
+ {
+ // Project on friction cone (Tasora et al., 2011; Fig. 2; Eq. (46)).
+ real_t pt_sq( p[1] * p[1] + p[2] * p[2] );
+ real_t pt_lim( mu_ * p[0] );
+
+ if( p[0] > -c_ && pt_sq < pt_lim * pt_lim )
+ {
+ // Unconstrained minimum is INSIDE THE UPPER CONE,
+ // with slope n*mu
+ // leading to a static contact and no projection is necessary.
+ }
+ else if( p[0] < -c_ && pt_sq < math::sq( p[0] / mu_ ) )
+ {
+ // 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 = Vec3( -c_, 0, 0 );
+ }
+ else
+ {
+ // Unconstrained minimum is OUTSIDE THE CONE -> Project on cone surface, Eq. (45) in Tasora et al., 2011
+ real_t pt( std::sqrt( pt_sq ) ); // gamma_r
+ p[0] = ( pt * mu_ + p[0] ) / ( math::sq( mu_ ) + 1 ); // p_cf[0] = gamma_n
+ real_t factor( mu_ * p[0] / pt );
+ p[1] *= factor;
+ p[2] *= factor;
+ }
+ }
+
+private:
+ const real_t mu_;
+ const real_t c_;
+}; // class Coulomb
+
+
+/**
+ * \ingroup pe
+ * \brief The Tresca law yield surface.
+ */
+class Tresca : public YieldSurface
+{
+public:
+ explicit Tresca( const real_t pt_cf_max )
+ : pt_cf_max_( pt_cf_max ) { }
+
+ bool contains( const Vec3 & p ) WALBERLA_OVERRIDE
+ {
+ return p[1] < pt_cf_max_;
+ }
+
+ void project( Vec3 & p ) WALBERLA_OVERRIDE
+ {
+ // Ratio of (squared) tangential reaction and (squared) tangential limit
+ const real_t factor = ( p[1] * p[1] + p[2] * p[2] ) / pt_cf_max_ * pt_cf_max_;
+
+ // Actual reaction exceeds the limit
+ if( factor > real_t(1) )
+ {
+ p[1] /= factor;
+ p[2] /= factor;
+ }
+ }
+
+private:
+ const real_t pt_cf_max_;
+}; // class Tresca
+
+
+///**
+// * \ingroup pe
+// * \brief The Shaw law yield surface with optional cohesion.
+// */
+//class Shaw : public YieldSurface
+//{
+//public:
+// Shaw( real_t k, real_t c = real_t(0) )
+// : k_( k ), c_( c ) { }
+//
+// bool contains( const Vec3 & p_cf ) WALBERLA_OVERRIDE
+// {
+//
+// }
+//
+// void project( Vec3 & p_cf ) WALBERLA_OVERRIDE
+// {
+//
+// }
+//
+//
+//private:
+// const real_t k_;
+// const real_t c_;
+//}; // class Shaw
+
+
+/**
+ * \ingroup pe
+ * \brief The Coulomb-Orowan yield surface with optional cohesion.
+ */
+class CoulombOrowan : public YieldSurface
+{
+public:
+
+ /**
+ * \brief The constructor.
+ * \param mu The coulomb friction coefficient \in [0, 1.0].
+ * \param k The elastic limit or yield strength contact reaction of the material.
+ * \param c The optional cohesion of the material.
+ */
+ CoulombOrowan( real_t mu, real_t k, real_t c = real_t(0) )
+ : mu_( mu ), k_( k ), c_( c ), pnc_( k / mu ) { }
+
+ bool contains( const Vec3 & p ) WALBERLA_OVERRIDE
+ {
+ return ( p[0] > -c_ && p[1] < math::min( p[0] * mu_, k_ ) );
+ }
+
+ void project( Vec3 & p ) WALBERLA_OVERRIDE
+ {
+ // Project on friction cone (Tasora et al., 2011; Fig. 2; Eq. (46)).
+ real_t pt_sq(p[1] * p[1] + p[2] * p[2] );
+ real_t pt_lim( math::min(mu_ * p[0], k_ ) );
+
+ if(p[0] >= pnc_ && pt_sq > k_ * k_ )
+ {
+ // OUTSIDE THE UPPER CONE ABOVE THE CRITICAL NORMAL REACTION
+ // where k is the limit tangential reaction
+ const real_t factor = k_ * k_ / pt_sq;
+ p[1] *= factor;
+ p[2] *= factor;
+ }
+ else if(p[0] > -c_ && pt_sq < pt_lim * pt_lim )
+ {
+ // INSIDE THE UPPER CONE,
+ // with slope n*mu
+ // leading to a static contact and no projection is necessary.
+ }
+ else if(p[0] < -c_ && pt_sq < math::sq(p[0] / mu_ ) )
+ {
+ // 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 = Vec3(-c_, 0, 0 );
+ }
+ else
+ {
+ // OUTSIDE THE UPPER CONE BELOW -> Project on cone surface, Eq. (45) in Tasora et al., 2011
+ real_t pt( std::sqrt( pt_sq ) ); // gamma_r
+ p[0] = (pt * mu_ + p[0] ) / (math::sq(mu_ ) + 1 ); // p[0] = gamma_n
+ const real_t factor( mu_ * p[0] / pt );
+ p[1] *= factor;
+ p[2] *= factor;
+ }
+ }
+
+
+private:
+ const real_t mu_;
+ const real_t k_;
+ const real_t c_;
+ const real_t pnc_;
+}; // class CoulombOrowan
+
+} // namespace cr
+} // namespace pe
+} // namespace walberla