New question #277532 on Yade:
https://answers.launchpad.net/yade/+question/277532
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
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
