Hello,
I want to modify the scripts c++ JointedCohesiveFroctionalPM.cpp and .hpp to calculate the dissiped energy.
After writing the modifications, i compile yade with "make install" in directory "YADE/sources/build".
When i run the code, i obtain this error message:
*** Error in `/usr/bin/python': free(): invalid pointer: 0x0a0f1718 ***
Abandon (core dumped)
this is the modified JointedCohesiveFroctionalPM.cpp script:
/* LucScholtes2010 */
#include"JointedCohesiveFrictionalPM.hpp"
#include<core/Scene.hpp>
#include<pkg/dem/ScGeom.hpp>
#include<core/Omega.hpp>
YADE_PLUGIN((JCFpmMat)(JCFpmState)(JCFpmPhys)(Ip2_JCFpmMat_JCFpmMat_JCFpmPhys)(Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM));
/********************** Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM ****************************/
CREATE_LOGGER(Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM);
bool Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
const int &id1 = contact->getId1();
const int &id2 = contact->getId2();
ScGeom* geom = static_cast<ScGeom*>(ig.get());
JCFpmPhys* phys = static_cast<JCFpmPhys*>(ip.get());
Body* b1 = Body::byId(id1,scene).get();
Body* b2 = Body::byId(id2,scene).get();
Real Dtensile=phys->FnMax/phys->kn;
string fileCracks = "cracks_"+Key+".txt";
//string
/// Defines the interparticular distance used for computation
Real D = 0;
/*this is for setting the equilibrium distance between all cohesive elements in the first contact detection*/
if ( contact->isFresh(scene) ) {
phys->normalForce = Vector3r::Zero();
phys->shearForce = Vector3r::Zero();
if ((smoothJoint) && (phys->isOnJoint)) {
phys->jointNormal = geom->normal.dot(phys->jointNormal)*phys->jointNormal; //to set the joint normal colinear with the interaction normal
phys->jointNormal.normalize();
phys->initD = std::abs((b1->state->pos - b2->state->pos).dot(phys->jointNormal)); // to set the initial gap as the equilibrium gap
} else {
phys->initD = geom->penetrationDepth;
}
}
if ( smoothJoint && phys->isOnJoint ) {
if ( phys->more || ( phys-> jointCumulativeSliding > (2*min(geom->radius1,geom->radius2)) ) ) {
if (!neverErase) return false;
else {
phys->shearForce = Vector3r::Zero();
phys->normalForce = Vector3r::Zero();
phys->isCohesive =0;
phys->FnMax = 0;
phys->FsMax = 0;
return true; // do we need this? -> yes if it ends the loop (avoid the following calculations)
}
} else {
D = phys->initD - std::abs((b1->state->pos - b2->state->pos).dot(phys->jointNormal));
}
} else {
D = geom->penetrationDepth - phys->initD;
}
phys->crackJointAperture = D<0? -D : 0.; // for DFNFlow
/* Determination of interaction */
if (D < 0) { //spheres do not touch
if (!phys->isCohesive) {
if (!neverErase) return false;
else {
phys->shearForce = Vector3r::Zero();
phys->normalForce = Vector3r::Zero();
phys->isCohesive =0;
phys->FnMax = 0;
phys->FsMax = 0;
return true; // do we need this? not sure -> yes, it ends the loop (avoid the following calculations)
}
}
if ( phys->isCohesive && (phys->FnMax>0) && (std::abs(D)>Dtensile) ) {
// update body state with the number of broken bonds
JCFpmState* st1=dynamic_cast<JCFpmState*>(b1->state.get());
JCFpmState* st2=dynamic_cast<JCFpmState*>(b2->state.get());
st1->tensBreak+=1;
st2->tensBreak+=1;
st1->tensBreakRel+=1.0/st1->noIniLinks;
st2->tensBreakRel+=1.0/st2->noIniLinks;
totalCracksE=totalCracksE+0.5*( ((phys->normalForce.norm()*phys->normalForce.norm())/phys->kn) + ((phys->shearForce.norm()*phys->shearForce.norm())/phys->ks) ); // written by Jabrane Hamdi
// create a text file to record properties of the broken bond (iteration, position, type (tensile), cross section, contact normal orientation and dissipated energy)
if (recordCracks) {
Real scalarNF=phys->normalForce.norm();
Real scalarSF=phys->shearForce.norm();
std::ofstream file (fileCracks.c_str(), !cracksFileExist ? std::ios::trunc : std::ios::app);
if(file.tellp()==0){ file <<"i p0 p1 p2 t s norm0 norm1 norm2 e"<<endl; }
Vector3r crackNormal=Vector3r::Zero();
if ((smoothJoint) && (phys->isOnJoint)) { crackNormal=phys->jointNormal; } else {crackNormal=geom->normal;}
file << boost::lexical_cast<string> ( scene->iter )<<" "<< boost::lexical_cast<string> ( geom->contactPoint[0] ) <<" "<< boost::lexical_cast<string> ( geom->contactPoint[1] ) <<" "<< boost::lexical_cast<string> ( geom->contactPoint[2] ) <<" "<< 0 <<" "<< boost::lexical_cast<string> ( 0.5*(geom->radius1+geom->radius2) ) <<" "<< boost::lexical_cast<string> ( crackNormal[0] ) <<" "<< boost::lexical_cast<string> ( crackNormal[1] ) <<" "<< boost::lexical_cast<string> ( crackNormal[2] ) <<" "<< boost::lexical_cast<string> ( 0.5*( ((scalarNF*scalarNF)/phys->kn) + ((scalarSF*scalarSF)/phys->ks) ) ) << endl;
}
cracksFileExist=true;
/// Timos
if (!neverErase) return false;
else {
phys->shearForce = Vector3r::Zero();
phys->normalForce = Vector3r::Zero();
phys->isCohesive =0;
phys->FnMax = 0;
phys->FsMax = 0;
phys->isBroken = true;
return true; // do we need this? -> yes, it ends the loop (avoid the following calculations)
}
}
}
/* NormalForce */
Real Fn = 0;
Fn = phys->kn*D;
/* ShearForce */
Vector3r& shearForce = phys->shearForce;
Real jointSliding=0;
/* Energy calculated by particles sliding written by Jabrane Hamdi */
totalSlipE=0;
if ( smoothJoint && phys->isOnJoint ) {
/// incremental formulation (OK?)
Vector3r relativeVelocity = (b2->state->vel - b1->state->vel); // angVel are not taken into account as particles on joint don't rotate ????
Vector3r slidingVelocity = relativeVelocity - phys->jointNormal.dot(relativeVelocity)*phys->jointNormal;
Vector3r incrementalSliding = slidingVelocity*scene->dt;
shearForce -= phys->ks*incrementalSliding;
jointSliding = incrementalSliding.norm();
phys->jointCumulativeSliding += jointSliding;
} else {
shearForce = geom->rotate(phys->shearForce);
const Vector3r& incrementalShear = geom->shearIncrement();
shearForce -= phys->ks*incrementalShear;
}
/* Mohr-Coulomb criterion */
Real maxFs = phys->FsMax + Fn*phys->tanFrictionAngle;
Real scalarShearForce = shearForce.norm();
if (scalarShearForce > maxFs) {
if (!phys->isCohesive) {
// why these lines? Can we really have purely normal compressive loading? I guess yes...
if (scalarShearForce != 0) {
Vector3r trialForce=shearForce;
// /// static friction sliding
shearForce*=maxFs/scalarShearForce;
/// energy dissipated by particles sliding written by Jabrane Hamdi
totalSlipE=totalSlipE+((1./phys->ks)*(trialForce-shearForce))/*plastic disp*/ .dot(shearForce)/*active force*/;
// /// dynamic friction sliding on joint - TEST
// if ( phys->isOnJoint ) {
// shearForce*=(Fn*std::tan(5*180/Mathr::PI))/scalarShearForce;
// } else {
// shearForce*=maxFs/scalarShearForce;
// }
/// create a text file to record properties of the slipping contact (iteration, position, type (slip), cross section, contact normal orientation and dissipated energy)
int geometryIndex1 = b1->shape->getClassIndex();
int geometryIndex2 = b2->shape->getClassIndex();
if ( (recordSlips) && (maxFs!=0) && (geometryIndex1==geometryIndex2) ) { // avoid non working forces -> one of the 2 tests should be enough, right
std::ofstream file (fileCracks.c_str(), !cracksFileExist ? std::ios::trunc : std::ios::app);
if(file.tellp()==0){ file <<"i p0 p1 p2 t s norm0 norm1 norm2 e"<<endl; }
Vector3r crackNormal=Vector3r::Zero();
if ((smoothJoint) && (phys->isOnJoint)) { crackNormal=phys->jointNormal; } else {crackNormal=geom->normal;}
file << boost::lexical_cast<string> ( scene->iter )<<" "<< boost::lexical_cast<string> ( geom->contactPoint[0] ) <<" "<< boost::lexical_cast<string> ( geom->contactPoint[1] ) <<" "<< boost::lexical_cast<string> ( geom->contactPoint[2] ) <<" "<< 2 <<" "<< boost::lexical_cast<string> ( 0.5*(geom->radius1+geom->radius2) ) <<" "<< boost::lexical_cast<string> ( crackNormal[0] ) <<" "<< boost::lexical_cast<string> ( crackNormal[1] ) <<" "<< boost::lexical_cast<string> ( crackNormal[2] ) <<" "<< boost::lexical_cast<string> ( ((1./phys->ks)*(trialForce-shearForce))/*plastic disp*/ .dot(shearForce)/*active force*/ ) << endl;
}
cracksFileExist=true;
} else {
shearForce=Vector3r::Zero();
}
if ((smoothJoint) && (phys->isOnJoint)) {phys->dilation=phys->jointCumulativeSliding*phys->tanDilationAngle-D; phys->initD+=(jointSliding*phys->tanDilationAngle);}
}
// take into account shear cracking -> are those lines critical? -> TODO testing with and without
if (phys->isCohesive) {
// update body state with the number of broken bonds
JCFpmState* st1=dynamic_cast<JCFpmState*>(b1->state.get());
JCFpmState* st2=dynamic_cast<JCFpmState*>(b2->state.get());
st1->shearBreak+=1;
st2->shearBreak+=1;
st1->shearBreakRel+=1.0/st1->noIniLinks;
st2->shearBreakRel+=1.0/st2->noIniLinks;
// create a text file to record properties of the broken bond (iteration, position, type (shear), cross section, contact normal orientation and dissipated energy)
if (recordCracks) {
Real scalarNF=phys->normalForce.norm();
Real scalarSF=phys->shearForce.norm();
std::ofstream file (fileCracks.c_str(), !cracksFileExist ? std::ios::trunc : std::ios::app);
if(file.tellp()==0){ file <<"i p0 p1 p2 t s norm0 norm1 norm2 e"<<endl; }
Vector3r crackNormal=Vector3r::Zero();
if ((smoothJoint) && (phys->isOnJoint)) { crackNormal=phys->jointNormal; } else {crackNormal=geom->normal;}
file << boost::lexical_cast<string> ( scene->iter )<<" "<< boost::lexical_cast<string> ( geom->contactPoint[0] ) <<" "<< boost::lexical_cast<string> ( geom->contactPoint[1] ) <<" "<< boost::lexical_cast<string> ( geom->contactPoint[2] ) <<" "<< 1 <<" "<< boost::lexical_cast<string> ( 0.5*(geom->radius1+geom->radius2) ) <<" "<< boost::lexical_cast<string> ( crackNormal[0] ) <<" "<< boost::lexical_cast<string> ( crackNormal[1] ) <<" "<< boost::lexical_cast<string> ( crackNormal[2] ) <<" "<< boost::lexical_cast<string> ( 0.5*( ((scalarNF*scalarNF)/phys->kn) + ((scalarSF*scalarSF)/phys->ks) ) ) << endl;
}
cracksFileExist=true;
// // option 1: set the contact properties to friction if in compression, delete contact if in tension
// phys->isBroken = true;
// phys->isCohesive = 0;
// phys->FnMax = 0;
// phys->FsMax = 0;
// // shearForce *= Fn*phys->tanFrictionAngle/scalarShearForce; // now or at the next timestep?
// if ( D < 0 ) { // spheres do not touch
// if (!neverErase) return false;
// else {
// phys->shearForce = Vector3r::Zero();
// phys->normalForce = Vector3r::Zero();
// return true; // do we need this? not sure -> yes, it ends the loop (avoid the following calculations)
// }
// }
// option 2: delete contact whatever compressive or tensile state
if (!neverErase) return false;
else {
phys->shearForce = Vector3r::Zero();
phys->normalForce = Vector3r::Zero();
phys->isCohesive =0;
phys->FnMax = 0;
phys->FsMax = 0;
phys->isBroken = true;
return true; // do we need this? -> yes, it ends the loop (avoid the following calculations)
}
}
}
/* Apply forces */
if ((smoothJoint) && (phys->isOnJoint)) { phys->normalForce = Fn*phys->jointNormal; } else { phys->normalForce = Fn*geom->normal; }
Vector3r f = phys->normalForce + shearForce;
/// applyForceAtContactPoint computes torque also and, for now, we don't want rotation for particles on joint (some errors in calculation due to specific geometry)
//applyForceAtContactPoint(f, geom->contactPoint, I->getId2(), b2->state->pos, I->getId1(), b1->state->pos, scene);
scene->forces.addForce (id1,-f);
scene->forces.addForce (id2, f);
// simple solution to avoid torque computation for particles interacting on a smooth joint
if ( (phys->isOnJoint)&&(smoothJoint) ) return true;
/// those lines are needed if rootBody->forces.addForce and rootBody->forces.addMoment are used instead of applyForceAtContactPoint -> NOTE need to check for accuracy!!!
scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(-f));
scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(-f));
return true;
}
/********************** Ip2_JCFpmMat_JCFpmMat_JCFpmPhys ****************************/
CREATE_LOGGER(Ip2_JCFpmMat_JCFpmMat_JCFpmPhys);
void Ip2_JCFpmMat_JCFpmMat_JCFpmPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction){
/* avoid updates of interaction if it already exists */
if( interaction->phys ) return;
ScGeom* geom=dynamic_cast<ScGeom*>(interaction->geom.get());
assert(geom);
const shared_ptr<JCFpmMat>& yade1 = YADE_PTR_CAST<JCFpmMat>(b1);
const shared_ptr<JCFpmMat>& yade2 = YADE_PTR_CAST<JCFpmMat>(b2);
JCFpmState* st1=dynamic_cast<JCFpmState*>(Body::byId(interaction->getId1(),scene)->state.get());
JCFpmState* st2=dynamic_cast<JCFpmState*>(Body::byId(interaction->getId2(),scene)->state.get());
shared_ptr<JCFpmPhys> contactPhysics(new JCFpmPhys());
/* From material properties */
Real E1 = yade1->young;
Real E2 = yade2->young;
Real v1 = yade1->poisson;
Real v2 = yade2->poisson;
Real f1 = yade1->frictionAngle;
Real f2 = yade2->frictionAngle;
Real SigT1 = yade1->tensileStrength;
Real SigT2 = yade2->tensileStrength;
Real Coh1 = yade1->cohesion;
Real Coh2 = yade2->cohesion;
/* From interaction geometry */
Real R1= geom->radius1;
Real R2= geom->radius2;
contactPhysics->crossSection = Mathr::PI*pow(min(R1,R2),2);
/* Pass values to JCFpmPhys. In case of a "jointed" interaction, the following values will be replaced by other ones later (in few if(){} blocks)*/
// frictional properties
contactPhysics->kn = 2.*E1*R1*E2*R2/(E1*R1+E2*R2);
if ( (v1==0)&&(v2==0) )
contactPhysics->ks = 0;
else
contactPhysics->ks = 2.*E1*R1*v1*E2*R2*v2/(E1*R1*v1+E2*R2*v2);
contactPhysics->tanFrictionAngle = std::tan(std::min(f1,f2));
// cohesive properties
///to set if the contact is cohesive or not
if ( ((cohesiveTresholdIteration < 0) || (scene->iter < cohesiveTresholdIteration)) && (std::min(SigT1,SigT2)>0 || std::min(Coh1,Coh2)>0) && (yade1->type == yade2->type)){
contactPhysics->isCohesive=true;
st1->noIniLinks++;
st2->noIniLinks++;
}
if ( contactPhysics->isCohesive ) {
contactPhysics->FnMax = std::min(SigT1,SigT2)*contactPhysics->crossSection;
contactPhysics->FsMax = std::min(Coh1,Coh2)*contactPhysics->crossSection;
}
/// +++ Jointed interactions ->NOTE: geom->normal is oriented from 1 to 2 / jointNormal from plane to sphere
if ( st1->onJoint && st2->onJoint )
{
if ( (((st1->jointNormal1.cross(st2->jointNormal1)).norm()<0.1) && (st1->jointNormal1.dot(st2->jointNormal1)<0)) || (((st1->jointNormal1.cross(st2->jointNormal2)).norm()<0.1) && (st1->jointNormal1.dot(st2->jointNormal2)<0)) || (((st1->jointNormal1.cross(st2->jointNormal3)).norm()<0.1) && (st1->jointNormal1.dot(st2->jointNormal3)<0)) )
{
contactPhysics->isOnJoint = true;
contactPhysics->jointNormal = st1->jointNormal1;
}
else if ( (((st1->jointNormal2.cross(st2->jointNormal1)).norm()<0.1) && (st1->jointNormal2.dot(st2->jointNormal1)<0)) || (((st1->jointNormal2.cross(st2->jointNormal2)).norm()<0.1) && (st1->jointNormal2.dot(st2->jointNormal2)<0)) || (((st1->jointNormal2.cross(st2->jointNormal3)).norm()<0.1) && (st1->jointNormal2.dot(st2->jointNormal3)<0)) )
{
contactPhysics->isOnJoint = true;
contactPhysics->jointNormal = st1->jointNormal2;
}
else if ( (((st1->jointNormal3.cross(st2->jointNormal1)).norm()<0.1) && (st1->jointNormal3.dot(st2->jointNormal1)<0)) || (((st1->jointNormal3.cross(st2->jointNormal2)).norm()<0.1) && (st1->jointNormal3.dot(st2->jointNormal2)<0)) || (((st1->jointNormal3.cross(st2->jointNormal3)).norm()<0.1) && (st1->jointNormal3.dot(st2->jointNormal3)<0)) )
{
contactPhysics->isOnJoint = true;
contactPhysics->jointNormal = st1->jointNormal3;
}
else if ( (st1->joint>3 || st2->joint>3) && ( ( ((st1->jointNormal1.cross(st2->jointNormal1)).norm()>0.1) && ((st1->jointNormal1.cross(st2->jointNormal2)).norm()>0.1) && ((st1->jointNormal1.cross(st2->jointNormal3)).norm()>0.1) ) || ( ((st1->jointNormal2.cross(st2->jointNormal1)).norm()>0.1) && ((st1->jointNormal2.cross(st2->jointNormal2)).norm()>0.1) && ((st1->jointNormal2.cross(st2->jointNormal3)).norm()>0.1) ) || ( ((st1->jointNormal3.cross(st2->jointNormal1)).norm()>0.1) && ((st1->jointNormal3.cross(st2->jointNormal2)).norm()>0.1) && ((st1->jointNormal3.cross(st2->jointNormal3)).norm()>0.1) ) ) ) { contactPhysics->isOnJoint = true; contactPhysics->more = true; contactPhysics->jointNormal = geom->normal; }
}
///to specify joint properties
if ( contactPhysics->isOnJoint ) {
Real jf1 = yade1->jointFrictionAngle;
Real jf2 = yade2->jointFrictionAngle;
Real jkn1 = yade1->jointNormalStiffness;
Real jkn2 = yade2->jointNormalStiffness;
Real jks1 = yade1->jointShearStiffness;
Real jks2 = yade2->jointShearStiffness;
Real jdil1 = yade1->jointDilationAngle;
Real jdil2 = yade2->jointDilationAngle;
Real jcoh1 = yade1->jointCohesion;
Real jcoh2 = yade2->jointCohesion;
Real jSigT1 = yade1->jointTensileStrength;
Real jSigT2 = yade2->jointTensileStrength;
contactPhysics->tanFrictionAngle = std::tan(std::min(jf1,jf2));
//contactPhysics->kn = jointNormalStiffness*2.*R1*R2/(R1+R2); // very first expression from Luc
//contactPhysics->kn = (jkn1+jkn2)/2.0*2.*R1*R2/(R1+R2); // after putting jointNormalStiffness in material
contactPhysics->kn = ( jkn1 + jkn2 ) /2.0 * contactPhysics->crossSection; // for a size independant expression
contactPhysics->ks = ( jks1 + jks2 ) /2.0 * contactPhysics->crossSection; // for a size independant expression
contactPhysics->tanDilationAngle = std::tan(std::min(jdil1,jdil2));
///to set if the contact is cohesive or not
if ( ((cohesiveTresholdIteration < 0) || (scene->iter < cohesiveTresholdIteration)) && (std::min(jcoh1,jcoh2)>0 || std::min(jSigT1,jSigT2)>0) ) {
contactPhysics->isCohesive=true;
st1->noIniLinks++;
st2->noIniLinks++;
}
else { contactPhysics->isCohesive=false; contactPhysics->FnMax=0; contactPhysics->FsMax=0; }
if ( contactPhysics->isCohesive ) {
contactPhysics->FnMax = std::min(jSigT1,jSigT2)*contactPhysics->crossSection;
contactPhysics->FsMax = std::min(jcoh1,jcoh2)*contactPhysics->crossSection;
}
}
interaction->phys = contactPhysics;
}
JCFpmPhys::~JCFpmPhys(){}
and this is the modified JointedCohesiveFroctionalPM.hpp script:
/* lucScholtes2010 */
#pragma once
#include<pkg/common/ElastMat.hpp>
#include<pkg/common/Dispatching.hpp>
#include<pkg/common/NormShearPhys.hpp>
#include<pkg/dem/ScGeom.hpp>
/** This class holds information associated with each body state*/
class JCFpmState: public State {
YADE_CLASS_BASE_DOC_ATTRS_CTOR(JCFpmState,State,"JCFpm state information about each body.",
((int,tensBreak,0,,"Number of tensile breakages. [-]"))
((int,shearBreak,0,,"Number of shear breakages. [-]"))
((int,noIniLinks,0,,"Number of initial cohesive interactions. [-]"))
((Real,tensBreakRel,0,,"Relative number (in [0;1], compared with :yref:`noIniLinks<JCFpmState.noIniLinks>`) of tensile breakages. [-]"))
((Real,shearBreakRel,0,,"Relative number (in [0;1], compared with :yref:`noIniLinks<JCFpmState.noIniLinks>`) of shear breakages. [-]"))
((bool,onJoint,false,,"Identifies if the particle is on a joint surface."))
((int,joint,0,,"Indicates the number of joint surfaces to which the particle belongs (0-> no joint, 1->1 joint, etc..). [-]"))
((Vector3r,jointNormal1,Vector3r::Zero(),,"Specifies the normal direction to the joint plane 1. Rk: the ideal here would be to create a vector of vector wich size is defined by the joint integer (as much joint normals as joints). However, it needs to make the pushback function works with python since joint detection is done through a python script. lines 272 to 312 of cpp file should therefore be adapted. [-]"))
((Vector3r,jointNormal2,Vector3r::Zero(),,"Specifies the normal direction to the joint plane 2. [-]"))
((Vector3r,jointNormal3,Vector3r::Zero(),,"Specifies the normal direction to the joint plane 3. [-]"))
,
createIndex();
);
REGISTER_CLASS_INDEX(JCFpmState,State);
};
REGISTER_SERIALIZABLE(JCFpmState);
/** This class holds information associated with each body */
class JCFpmMat: public FrictMat {
public:
virtual shared_ptr<State> newAssocState() const { return shared_ptr<State>(new JCFpmState); }
virtual bool stateTypeOk(State* s) const { return (bool)dynamic_cast<JCFpmState*>(s); }
YADE_CLASS_BASE_DOC_ATTRS_CTOR(JCFpmMat,FrictMat,"Possibly jointed, cohesive frictional material, for use with other JCFpm classes",
((Real,cohesion,0.,,"Defines the maximum admissible tangential force in shear, for Fn=0, in the matrix (:yref:`FsMax<JCFpmPhys.FsMax>` = cohesion * :yref:`crossSection<JCFpmPhys.crossSection>`). [Pa]"))
((Real,jointCohesion,0.,,"Defines the :yref:`maximum admissible tangential force in shear<JCFpmPhys.FsMax>`, for Fn=0, on the joint surface. [Pa]"))
((Real,jointDilationAngle,0,,"Defines the dilatancy of the joint surface (only valid for :yref:`smooth contact logic<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.smoothJoint>`). [rad]"))
((Real,jointFrictionAngle,-1,,"Defines Coulomb friction on the joint surface. [rad]"))
((Real,jointNormalStiffness,0.,,"Defines the normal stiffness on the joint surface. [Pa/m]"))
((Real,jointShearStiffness,0.,,"Defines the shear stiffness on the joint surface. [Pa/m]"))
((Real,jointTensileStrength,0.,,"Defines the :yref:`maximum admissible normal force in traction<JCFpmPhys.FnMax>` on the joint surface. [Pa]"))
((int,type,0,,"If particles of two different types interact, it will be with friction only (no cohesion).[-]"))
((Real,tensileStrength,0.,,"Defines the maximum admissible normal force in traction in the matrix (:yref:`FnMax<JCFpmPhys.FnMax>` = tensileStrength * :yref:`crossSection<JCFpmPhys.crossSection>`). [Pa]"))
,
createIndex();
);
REGISTER_CLASS_INDEX(JCFpmMat,FrictMat);
};
REGISTER_SERIALIZABLE(JCFpmMat);
/** This class holds information associated with each interaction */
class JCFpmPhys: public NormShearPhys {
public:
virtual ~JCFpmPhys();
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(JCFpmPhys,NormShearPhys,"Representation of a single interaction of the JCFpm type, storage for relevant parameters",
((Real,initD,0,,"equilibrium distance for interacting particles. Computed as the interparticular distance at first contact detection."))
((bool,isCohesive,false,,"If false, particles interact in a frictional way. If true, particles are bonded regarding the given :yref:`cohesion<JCFpmMat.cohesion>` and :yref:`tensile strength<JCFpmMat.tensileStrength>` (or their jointed variants)."))
((bool,more,false,,"specifies if the interaction is crossed by more than 3 joints. If true, interaction is deleted (temporary solution)."))
((bool,isOnJoint,false,,"defined as true when both interacting particles are :yref:`on joint<JCFpmState.onJoint>` and are in opposite sides of the joint surface. In this case, mechanical parameters of the interaction are derived from the ''joint...'' material properties of the particles. Furthermore, the normal of the interaction may be re-oriented (see :yref:`Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.smoothJoint`)."))
((Real,tanFrictionAngle,0,,"tangent of Coulomb friction angle for this interaction (auto. computed). [-]"))
((Real,crossSection,0,,"crossSection=pi*Rmin^2. [m2]"))
((Real,FnMax,0,,"positiv value computed from :yref:`tensile strength<JCFpmMat.tensileStrength>` (or joint variant) to define the maximum admissible normal force in traction: Fn >= -FnMax. [N]"))
((Real,FsMax,0,,"computed from :yref:`cohesion<JCFpmMat.cohesion>` (or jointCohesion) to define the maximum admissible tangential force in shear, for Fn=0. [N]"))
((Vector3r,jointNormal,Vector3r::Zero(),,"normal direction to the joint, deduced from e.g. :yref:`<JCFpmState.jointNormal1>`."))
((Real,jointCumulativeSliding,0,,"sliding distance for particles interacting on a joint. Used, when :yref:`<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.smoothJoint>` is true, to take into account dilatancy due to shearing. [-]"))
((Real,tanDilationAngle,0,,"tangent of the angle defining the dilatancy of the joint surface (auto. computed from :yref:`JCFpmMat.jointDilationAngle`). [-]"))
((Real,dilation,0,,"defines the normal displacement in the joint after sliding treshold. [m]"))
((bool,isBroken,false,,"flag for broken interactions"))
((Real,crackJointAperture,0,,"Relative displacement between 2 spheres (in case of a crack it is equivalent of the crack aperture)"))
,
createIndex();
,
);
DECLARE_LOGGER;
REGISTER_CLASS_INDEX(JCFpmPhys,NormShearPhys);
};
REGISTER_SERIALIZABLE(JCFpmPhys);
/** 2d functor creating InteractionPhysics (Ip2) taking JCFpmMat and JCFpmMat of 2 bodies, returning type JCFpmPhys */
class Ip2_JCFpmMat_JCFpmMat_JCFpmPhys: public IPhysFunctor{
public:
virtual void go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction);
FUNCTOR2D(JCFpmMat,JCFpmMat);
DECLARE_LOGGER;
YADE_CLASS_BASE_DOC_ATTRS(Ip2_JCFpmMat_JCFpmMat_JCFpmPhys,IPhysFunctor,"Converts 2 :yref:`JCFpmMat` instances to one :yref:`JCFpmPhys` instance, with corresponding parameters.",
((int,cohesiveTresholdIteration,1,,"should new contacts be cohesive? If strictly negativ, they will in any case. If positiv, they will before this iter, they won't afterward."))
);
};
REGISTER_SERIALIZABLE(Ip2_JCFpmMat_JCFpmMat_JCFpmPhys);
/** 2d functor creating the interaction law (Law2) based on SphereContactGeometry (ScGeom) and JCFpmPhys of 2 bodies, returning type JointedCohesiveFrictionalPM */
class Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM: public LawFunctor{
public:
virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
FUNCTOR2D(ScGeom,JCFpmPhys);
YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM,LawFunctor,"Interaction law for cohesive frictional material, e.g. rock, possibly presenting joint surfaces, that can be mechanically described with a smooth contact logic [Ivars2011]_ (implemented in Yade in [Scholtes2012]_). See examples/jointedCohesiveFrictionalPM for script examples. Joint surface definitions (through stl meshes or direct definition with gts module) are illustrated there.",
((string,Key,"",,"string specifying the name of saved file 'cracks___.txt', when :yref:`recordCracks<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.recordCracks>` is true."))
((bool,cracksFileExist,false,,"if true (and if :yref:`recordCracks<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.recordCracks>`), data are appended to an existing 'cracksKey' text file; otherwise its content is reset."))
((bool,smoothJoint,false,,"if true, interactions of particles belonging to joint surface (:yref:`JCFpmPhys.isOnJoint`) are handled according to a smooth contact logic [Ivars2011]_, [Scholtes2012]_."))
((bool,recordCracks,false,,"if true, data about cohesive interactions that break are stored in a text file cracksKey.txt (see :yref:`Key<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.Key>` and :yref:`cracksFileExist<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.cracksFileExist>`). It contains 9 columns: the break iteration, the 3 coordinates of the contact point, the type (1 means shear break, while 0 corresponds to tensile break), the ''cross section'' (mean radius of the 2 spheres), the 3 coordinates of the contact normal and the energy released."))
((bool,recordSlips,false,,"if true, data about frictional interactions that slip are stored in a text file cracksKey.txt."))
((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene"))
((Real,totalSlipE,0.,,"calculate the sum of the energy dissipated by particles sliding we can get the value in the recorder through interactionLaw.totalSlipE")) // written by Jabrane Hamdi
((Real,totalCracksE,0.,,"calculate the sum of the energy dissipated by particles contact broken we can get the value in the recorder through interactionLaw.totalSlipE")) // written by Jabrane Hamdi
);
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM);
Jabrane
