Author: eudoxos
Date: 2009-06-04 16:17:27 +0200 (Thu, 04 Jun 2009)
New Revision: 1789
Modified:
trunk/pkg/dem/ConcretePM.cpp
trunk/pkg/dem/ConcretePM.hpp
Log:
Fixes and cleanups in the CPM code.
Modified: trunk/pkg/dem/ConcretePM.cpp
===================================================================
--- trunk/pkg/dem/ConcretePM.cpp 2009-06-02 22:02:25 UTC (rev 1788)
+++ trunk/pkg/dem/ConcretePM.cpp 2009-06-04 14:17:27 UTC (rev 1789)
@@ -118,54 +118,59 @@
/********************** Law2_Dem3DofGeom_CpmPhys_Cpm
****************************/
-Real Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2=1.; /* deactivated if >
0 */
-Real Law2_Dem3DofGeom_CpmPhys_Cpm::yieldLogSpeed=1.;
-Real Law2_Dem3DofGeom_CpmPhys_Cpm::yieldEllipseShift=0.;
+/// yield surface parameters
+int Law2_Dem3DofGeom_CpmPhys_Cpm::yieldSurfType=2;
+Real Law2_Dem3DofGeom_CpmPhys_Cpm::yieldLogSpeed=.1;
+Real Law2_Dem3DofGeom_CpmPhys_Cpm::yieldEllipseShift=0.; // deprecated
+/// compressive plasticity parameters
+// approximates confinement -20MPa precisely, -100MPa a little over, -200 and
-400 are OK (secant)
+Real Law2_Dem3DofGeom_CpmPhys_Cpm::epsSoft=-3e-3; // deactivated if >=0
+Real Law2_Dem3DofGeom_CpmPhys_Cpm::relKnSoft=.25;
+// >=1. to deactivate (never delete damaged contacts)
Real Law2_Dem3DofGeom_CpmPhys_Cpm::omegaThreshold=1.;
-Real Law2_Dem3DofGeom_CpmPhys_Cpm::epsSoft=0.;
-Real Law2_Dem3DofGeom_CpmPhys_Cpm::relKnSoft=.5;
void Law2_Dem3DofGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom,
shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody){
- //timingDeltas->start();
Dem3DofGeom* contGeom=static_cast<Dem3DofGeom*>(_geom.get());
CpmPhys* BC=static_cast<CpmPhys*>(_phys.get());
// shorthands
- Real& epsN(BC->epsN); Real& epsNPl(BC->epsNPl); Vector3r&
epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); const
Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const
Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real&
crossSection(BC->crossSection); const Real&
omegaThreshold(Law2_Dem3DofGeom_CpmPhys_Cpm::omegaThreshold); const Real&
epsCrackOnset(BC->epsCrackOnset); Real&
relResidualStrength(BC->relResidualStrength); const Real&
dt=Omega::instance().getTimeStep(); const Real& epsFracture(BC->epsFracture);
const bool& neverDamage(BC->neverDamage); const Real& dmgTau(BC->dmgTau); const
Real& plTau(BC->plTau); const bool& isCohesive(BC->isCohesive);
- Real& omega(BC->omega); Real& sigmaN(BC->sigmaN); Vector3r&
sigmaT(BC->sigmaT); Real& Fn(BC->Fn); Vector3r& Fs(BC->Fs); // for python access
- const Real& yieldLogSpeed(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldLogSpeed);
const int& yieldSurfType(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldSurfType);
- const Real&
yieldEllipseShift(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldEllipseShift);
- const Real& epsSoft(Law2_Dem3DofGeom_CpmPhys_Cpm::epsSoft);
- const Real& relKnSoft(Law2_Dem3DofGeom_CpmPhys_Cpm::relKnSoft);
+ Real& epsN(BC->epsN); Real& epsNPl(BC->epsNPl); Vector3r&
epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); const
Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const
Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real&
crossSection(BC->crossSection); const Real&
omegaThreshold(Law2_Dem3DofGeom_CpmPhys_Cpm::omegaThreshold); const Real&
epsCrackOnset(BC->epsCrackOnset); Real&
relResidualStrength(BC->relResidualStrength); const Real&
dt=Omega::instance().getTimeStep(); const Real& epsFracture(BC->epsFracture);
const bool& neverDamage(BC->neverDamage); const Real& dmgTau(BC->dmgTau); const
Real& plTau(BC->plTau); const bool& isCohesive(BC->isCohesive);
+ Real& omega(BC->omega); Real& sigmaN(BC->sigmaN); Vector3r&
sigmaT(BC->sigmaT); Real& Fn(BC->Fn); Vector3r& Fs(BC->Fs); // for python access
+ const Real&
yieldLogSpeed(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldLogSpeed); const int&
yieldSurfType(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldSurfType);
+ const Real&
yieldEllipseShift(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldEllipseShift);
+ const Real& epsSoft(Law2_Dem3DofGeom_CpmPhys_Cpm::epsSoft);
+ const Real& relKnSoft(Law2_Dem3DofGeom_CpmPhys_Cpm::relKnSoft);
- #define YADE_VERIFY(condition) if(!(condition)){LOG_FATAL("Verification
`"<<#condition<<"' failed!"); throw;}
-
- #define NNAN(a) YADE_VERIFY(!isnan(a));
- #define NNANV(v) YADE_VERIFY(!isnan(v[0])); assert(!isnan(v[1]));
assert(!isnan(v[2]));
- //timingDeltas->checkpoint("setup");
- // if(contGeom->refR1<0) contGeom->refLength=contGeom->refR2; // make
facet-sphere contact always at equilibrium when touching exactly (and not the
initial distance)
epsN=contGeom->strainN(); epsT=contGeom->strainT();
- #ifdef YADE_DEBUG
- if(isnan(epsN)){
- LOG_FATAL("refLength="<<contGeom->refLength<<";
pos1="<<contGeom->se31.position<<"; pos2="<<contGeom->se32.position<<";
displacementN="<<contGeom->displacementN());
+
+ // debugging
+ #define YADE_VERIFY(condition)
if(!(condition)){LOG_FATAL("Verification `"<<#condition<<"' failed!"); throw;}
+ #define NNAN(a) YADE_VERIFY(!isnan(a));
+ #define NNANV(v) YADE_VERIFY(!isnan(v[0]));
assert(!isnan(v[1])); assert(!isnan(v[2]));
+ #ifdef YADE_DEBUG
+ if(isnan(epsN)){
+ LOG_FATAL("refLength="<<contGeom->refLength<<";
pos1="<<contGeom->se31.position<<"; pos2="<<contGeom->se32.position<<";
displacementN="<<contGeom->displacementN());
throw runtime_error("!! epsN==NaN !!");
- }
+ }
+ NNAN(epsN); NNANV(epsT);
+ #endif
NNAN(epsN); NNANV(epsT);
- #endif
- // already in SpheresContactGeometry:
- // contGeom->relocateContactPoints(); // allow very large mutual
rotations
- if(logStrain && epsN<0){
- Real epsN0=epsN;
- epsN=log(epsN0+1); epsT*=epsN/epsN0;
- }
- NNAN(epsN); NNANV(epsT);
- //timingDeltas->checkpoint("geom");
- epsN+=BC->isoPrestress/E;
+ // constitutive law
#ifdef CPM_MATERIAL_MODEL
+ // complicated version
+ if(epsSoft>=0) epsN+=BC->isoPrestress/E;
+ else{ // take softening into account for the prestress
+ Real sigmaSoft=E*epsSoft;
+ if(BC->isoPrestress>=sigmaSoft)
epsN+=BC->isoPrestress/E; // on the non-softened branch yet
+ // otherwise take the regular and softened branches
separately (different moduli)
+ else
epsN+=sigmaSoft/E+(BC->isoPrestress-sigmaSoft)/(E*relKnSoft);
+ }
CPM_MATERIAL_MODEL
#else
+ // simplified public model
+ epsN+=BC->isoPrestress/E;
// very simplified version of the constitutive law
kappaD=max(max(0.,epsN),kappaD); // internal variable, max
positive strain (non-decreasing)
omega=isCohesive?funcG(kappaD,epsCrackOnset,epsFracture,neverDamage):1.; //
damage variable (non-decreasing, as funcG is also non-decreasing)
@@ -179,18 +184,11 @@
relResidualStrength=isCohesive?(kappaD<epsCrackOnset?1.:(1-omega)*(kappaD)/epsCrackOnset):0;
#endif
sigmaN-=BC->isoPrestress;
- if(contGeom->refR1<0 &&
Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2<=0 &&
epsN<Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2){
- /* move Body2 (the sphere) so that minStrain is satisfied */
-
rootBody->bex.addMove(I->getId2(),contGeom->normal*(Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2-epsN)*contGeom->refLength);
- LOG_TRACE("Moving by
"<<contGeom->normal*(Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2-epsN)*contGeom->refLength);
- }
+
NNAN(kappaD); NNAN(epsCrackOnset); NNAN(epsFracture); NNAN(omega);
NNAN(sigmaN); NNANV(sigmaT); NNAN(crossSection);
- //timingDeltas->checkpoint("material");
- //const int watch1=6300, watch2=6299;
- //#define SHOW(a) if((I->getId1()==watch1 && I->getId2()==watch2) ||
(I->getId2()==watch1 && I->getId1()==watch2)) cerr<<__FILE__<<":"<<__LINE__<<"
"<<a<<endl;
- //SHOW("epsN"<<epsN);
+ // handle broken contacts
if(epsN>0. && ((isCohesive && omega>omegaThreshold) || !isCohesive)){
rootBody->interactions->requestErase(I->getId1(),I->getId2());
if(isCohesive){
@@ -206,7 +204,6 @@
Fs=sigmaT*crossSection; BC->shearForce=Fs;
applyForceAtContactPoint(BC->normalForce+BC->shearForce,
contGeom->contactPoint, I->getId1(), contGeom->se31.position, I->getId2(),
contGeom->se32.position, rootBody);
- //timingDeltas->checkpoint("rest");
}
Modified: trunk/pkg/dem/ConcretePM.hpp
===================================================================
--- trunk/pkg/dem/ConcretePM.hpp 2009-06-02 22:02:25 UTC (rev 1788)
+++ trunk/pkg/dem/ConcretePM.hpp 2009-06-04 14:17:27 UTC (rev 1789)
@@ -248,26 +248,23 @@
if(kappaD<epsCrackOnset || neverDamage) return 0;
return
1.-(epsCrackOnset/kappaD)*exp(-(kappaD-epsCrackOnset)/epsFracture);
}
- bool logStrain;
//! yield function: 0: mohr-coulomb (original); 1: parabolic;
2: logarithmic, 3: log+lin_tension, 4: elliptic, 5: elliptic+log
- int yieldSurfType;
+ static int yieldSurfType;
//! scaling in the logarithmic yield surface (should be <1 for
realistic results; >=0 for meaningful results)
static Real yieldLogSpeed;
//! horizontal scaling of the ellipse (shifts on the +x axis as
interactions with +y are given)
static Real yieldEllipseShift;
//! damage after which the contact disappears (<1), since omega
reaches 1 only for strain →+∞
static Real omegaThreshold;
- //! HACK: limit strain on some contacts by moving body #2 in
the contact; only if refR1<0 (facet); deactivated if > 0
- static Real minStrain_moveBody2;
//! Strain at which softening in compression starts (set to 0.
(default) or positive value to deactivate)
static Real epsSoft;
//! Relative rigidity of the softening branch in compression
(0=perfect elastic-plastic, 1=no softening, >1=hardening)
static Real relKnSoft;
- Law2_Dem3DofGeom_CpmPhys_Cpm(): logStrain(false),
yieldSurfType(0) { /*timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);*/
}
+ Law2_Dem3DofGeom_CpmPhys_Cpm() { }
void go(shared_ptr<InteractionGeometry>& _geom,
shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody);
FUNCTOR2D(Dem3DofGeom,CpmPhys);
REGISTER_CLASS_AND_BASE(Law2_Dem3DofGeom_CpmPhys_Cpm,ConstitutiveLaw);
-
REGISTER_ATTRIBUTES(ConstitutiveLaw,(logStrain)(yieldSurfType)(yieldLogSpeed)(yieldEllipseShift)(minStrain_moveBody2)(omegaThreshold)(epsSoft)(relKnSoft));
+
REGISTER_ATTRIBUTES(ConstitutiveLaw,(yieldSurfType)(yieldLogSpeed)(yieldEllipseShift)(omegaThreshold)(epsSoft)(relKnSoft));
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(Law2_Dem3DofGeom_CpmPhys_Cpm);
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
