yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #13112
[Branch ~yade-pkg/yade/git-trunk] Rev 4040: update of the JCFPM particle model with new or adapted functionalities to record number of cracks...
------------------------------------------------------------
revno: 4040
committer: luc scholtes <lscholtes63@xxxxxxxxx>
timestamp: Thu 2017-05-11 13:42:34 +0200
message:
update of the JCFPM particle model with new or adapted functionalities to record number of cracks and associated released energy. An example script will be provided.
modified:
pkg/dem/JointedCohesiveFrictionalPM.cpp
pkg/dem/JointedCohesiveFrictionalPM.hpp
pkg/dem/VTKRecorder.cpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'pkg/dem/JointedCohesiveFrictionalPM.cpp'
--- pkg/dem/JointedCohesiveFrictionalPM.cpp 2014-11-05 15:08:12 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.cpp 2017-05-11 11:42:34 +0000
@@ -27,7 +27,7 @@
/// 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*/
+ /*this is for setting the equilibrium distance between all cohesive elements at the first contact detection*/
if ( contact->isFresh(scene) ) {
phys->normalForce = Vector3r::Zero();
phys->shearForce = Vector3r::Zero();
@@ -49,7 +49,7 @@
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)
+ return true;
}
} else {
D = phys->initD - std::abs((b1->state->pos - b2->state->pos).dot(phys->jointNormal));
@@ -70,29 +70,34 @@
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)
+ return true;
}
}
if ( phys->isCohesive && (phys->FnMax>0) && (std::abs(D)>Dtensile) ) {
- // update body state with the number of broken bonds
+ nbTensCracks++;
+ // update body state with the number of broken bonds -> do we really need that?
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;
+ st1->nbBrokenBonds++;
+ st2->nbBrokenBonds++;
+ st1->damageIndex+=1.0/st1->nbInitBonds;
+ st2->damageIndex+=1.0/st2->nbInitBonds;
+
+ Real scalarNF=phys->normalForce.norm();
+ Real scalarSF=phys->shearForce.norm();
+ totalTensCracksE+=0.5*( ((scalarNF*scalarNF)/phys->kn) + ((scalarSF*scalarSF)/phys->ks) );
- // create a text file to record properties of the broken bond (iteration, position, type (tensile), cross section and contact normal orientation)
if (recordCracks){
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"<<endl; }
+ if(file.tellp()==0){ file <<"iter time p0 p1 p2 type size norm0 norm1 norm2 nrg"<<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] ) << endl;
+ file << boost::lexical_cast<string> ( scene->iter ) << " " << boost::lexical_cast<string> ( scene->time ) <<" "<< 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;
+
/// Timos
if (!neverErase) return false;
else {
@@ -102,9 +107,8 @@
phys->FnMax = 0;
phys->FsMax = 0;
phys->isBroken = true;
- return true; // do we need this? not sure -> yes, it ends the loop (avoid the following calculations)
+ return true;
}
-// return true; // do we need this? no
}
}
@@ -145,42 +149,60 @@
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) {
+//
+// nbSlips++;
+// totalSlipE+=((1./phys->ks)*(trialForce-shearForce))/*plastic disp*/.dot(shearForce)/*active force*/;
+//
+// if ( (recordSlips) && (maxFs!=0) ) {
+// std::ofstream file (fileCracks.c_str(), !cracksFileExist ? std::ios::trunc : std::ios::app);
+// if(file.tellp()==0){ file <<"iter time p0 p1 p2 type size norm0 norm1 norm2 nrg"<<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> ( scene->time ) <<" "<< 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> ( ((1./phys->ks)*(trialForce-shearForce)).dot(shearForce) ) << endl;
+// }
+// cracksFileExist=true;
+// }
+
if ( phys->isCohesive ) {
- // update body state with the number of broken bonds
+ nbShearCracks++;
+ // update body state with the number of broken bonds -> do we really need that?
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 and contact normal orientation)
+ st1->nbBrokenBonds++;
+ st2->nbBrokenBonds++;
+ st1->damageIndex+=1.0/st1->nbInitBonds;
+ st2->damageIndex+=1.0/st2->nbInitBonds;
+
+ Real scalarNF=phys->normalForce.norm();
+ Real scalarSF=phys->shearForce.norm();
+ totalShearCracksE+=0.5*( ((scalarNF*scalarNF)/phys->kn) + ((scalarSF*scalarSF)/phys->ks) );
+
if (recordCracks){
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"<<endl; }
+ if(file.tellp()==0){ file <<"iter time p0 p1 p2 type size norm0 norm1 norm2 nrg"<<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] ) << endl;
+ file << boost::lexical_cast<string> ( scene->iter ) << " " << boost::lexical_cast<string> ( scene->time ) <<" "<< 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> ( 0.5*( ((scalarNF*scalarNF)/phys->kn) + ((scalarSF*scalarSF)/phys->ks) ) ) <<endl;
}
cracksFileExist=true;
-
+
// 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?
+// 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)
+ return true;
}
}
-// return true; // do we need this one? no
}
}
@@ -252,8 +274,8 @@
///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++;
+ st1->nbInitBonds++;
+ st2->nbInitBonds++;
}
if ( contactPhysics->isCohesive ) {
@@ -308,15 +330,15 @@
///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++;
+ contactPhysics->isCohesive=true;
+ st1->nbInitBonds++;
+ st2->nbInitBonds++;
}
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;
+ contactPhysics->FnMax = std::min(jSigT1,jSigT2)*contactPhysics->crossSection;
+ contactPhysics->FsMax = std::min(jcoh1,jcoh2)*contactPhysics->crossSection;
}
}
interaction->phys = contactPhysics;
=== modified file 'pkg/dem/JointedCohesiveFrictionalPM.hpp'
--- pkg/dem/JointedCohesiveFrictionalPM.hpp 2016-01-13 18:56:48 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.hpp 2017-05-11 11:42:34 +0000
@@ -10,17 +10,15 @@
/** 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. [-]"))
+ ((int,nbInitBonds,0,,"Number of initial bonds. [-]"))
+ ((int,nbBrokenBonds,0,,"Number of broken bonds. [-]"))
+ ((Real,damageIndex,0,,"Ratio of broken bonds over initial bonds. [-]"))
((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);
@@ -34,15 +32,15 @@
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",
+ ((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]"))
((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,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]"))
((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();
);
@@ -57,9 +55,11 @@
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,isBroken,false,,"flag for broken interactions"))
((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`)."))
+ ((bool,isOnSlot,false,,"defined as true when interaction is located in the perforation slot (surface)."))
((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]"))
@@ -68,7 +68,6 @@
((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();
@@ -86,8 +85,7 @@
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. See :yref:`JCFpmMat` and [Duriez2016]_ for details",
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_JCFpmMat_JCFpmMat_JCFpmPhys,IPhysFunctor,"Converts 2 :yref:`JCFpmMat` instances to one :yref:`JCFpmPhys` instance, with corresponding parameters. See :yref:`JCFpmMat` and [Duriez2016]_ for details",
((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."))
);
@@ -101,11 +99,18 @@
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 interactions that lose their cohesive feature 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) and the 3 coordinates of the contact normal."))
((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene"))
+ ((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."))
+ ((string,Key,"",,"string specifying the name of saved file 'cracks___.txt', when :yref:`recordCracks<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.recordCracks>` is true."))
+ ((bool,recordCracks,false,,"if true, data about interactions that lose their cohesive feature are stored in the 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) and the 3 coordinates of the contact normal."))
+ ((int,nbTensCracks,0,,"number of tensile microcracks."))
+ ((int,nbShearCracks,0,,"number of shear microcracks."))
+ ((Real,totalTensCracksE,0.,,"calculate the overall energy dissipated by interparticle microcracking in tension."))
+ ((Real,totalShearCracksE,0.,,"calculate the overall energy dissipated by interparticle microcracking in shear."))
+// ((bool,recordSlips,false,,"if true, data about frictional interactions that slip are stored in a text file cracksKey.txt."))
+// ((int,nbSlips,0,,"number of slips."))
+// ((Real,totalSlipE,0.,,"calculate the overall energy dissipated by interparticle friction."))
);
DECLARE_LOGGER;
};
=== modified file 'pkg/dem/VTKRecorder.cpp'
--- pkg/dem/VTKRecorder.cpp 2016-05-13 13:07:22 +0000
+++ pkg/dem/VTKRecorder.cpp 2017-05-11 11:42:34 +0000
@@ -102,7 +102,6 @@
// liquid control needs interactions
if(recActive[REC_LIQ]) recActive[REC_INTR]=true;
-
// spheres
vtkSmartPointer<vtkPoints> spheresPos = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> spheresCells = vtkSmartPointer<vtkCellArray>::New();
@@ -341,12 +340,12 @@
cpmStress->SetName("cpmStress");
// extras for JCFpm
- vtkSmartPointer<vtkDoubleArray> damage = vtkSmartPointer<vtkDoubleArray>::New();
- damage->SetNumberOfComponents(1);
- damage->SetName("damage");
- vtkSmartPointer<vtkDoubleArray> damageRel = vtkSmartPointer<vtkDoubleArray>::New();
- damageRel->SetNumberOfComponents(1);
- damageRel->SetName("damageRel");
+ vtkSmartPointer<vtkDoubleArray> nbCracks = vtkSmartPointer<vtkDoubleArray>::New();
+ nbCracks->SetNumberOfComponents(1);
+ nbCracks->SetName("nbCracks");
+ vtkSmartPointer<vtkDoubleArray> jcfpmDamage = vtkSmartPointer<vtkDoubleArray>::New();
+ jcfpmDamage->SetNumberOfComponents(1);
+ jcfpmDamage->SetName("damage");
vtkSmartPointer<vtkDoubleArray> intrIsCohesive = vtkSmartPointer<vtkDoubleArray>::New();
intrIsCohesive->SetNumberOfComponents(1);
intrIsCohesive->SetName("isCohesive");
@@ -357,18 +356,24 @@
// extras for cracks
vtkSmartPointer<vtkPoints> crackPos = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> crackCells = vtkSmartPointer<vtkCellArray>::New();
- vtkSmartPointer<vtkDoubleArray> iter = vtkSmartPointer<vtkDoubleArray>::New();
- iter->SetNumberOfComponents(1);
- iter->SetName("iter");
+ vtkSmartPointer<vtkDoubleArray> crackIter = vtkSmartPointer<vtkDoubleArray>::New();
+ crackIter->SetNumberOfComponents(1);
+ crackIter->SetName("iter");
+ vtkSmartPointer<vtkDoubleArray> crackTime = vtkSmartPointer<vtkDoubleArray>::New();
+ crackTime->SetNumberOfComponents(1);
+ crackTime->SetName("time");
vtkSmartPointer<vtkDoubleArray> crackType = vtkSmartPointer<vtkDoubleArray>::New();
crackType->SetNumberOfComponents(1);
- crackType->SetName("crackType");
+ crackType->SetName("type");
vtkSmartPointer<vtkDoubleArray> crackSize = vtkSmartPointer<vtkDoubleArray>::New();
crackSize->SetNumberOfComponents(1);
- crackSize->SetName("crackSize");
+ crackSize->SetName("size");
vtkSmartPointer<vtkDoubleArray> crackNorm = vtkSmartPointer<vtkDoubleArray>::New();
crackNorm->SetNumberOfComponents(3);
- crackNorm->SetName("crackNorm");
+ crackNorm->SetName("norm");
+ vtkSmartPointer<vtkDoubleArray> crackNrg = vtkSmartPointer<vtkDoubleArray>::New();
+ crackNrg->SetNumberOfComponents(1);
+ crackNrg->SetName("nrg");
#ifdef YADE_LIQMIGRATION
vtkSmartPointer<vtkDoubleArray> liqVol = vtkSmartPointer<vtkDoubleArray>::New();
@@ -624,8 +629,8 @@
}
if (recActive[REC_JCFPM]){
- damage->InsertNextValue(YADE_PTR_CAST<JCFpmState>(b->state)->tensBreak + YADE_PTR_CAST<JCFpmState>(b->state)->shearBreak);
- damageRel->InsertNextValue(YADE_PTR_CAST<JCFpmState>(b->state)->tensBreakRel + YADE_PTR_CAST<JCFpmState>(b->state)->shearBreakRel);
+ nbCracks->InsertNextValue(YADE_PTR_CAST<JCFpmState>(b->state)->nbBrokenBonds);
+ jcfpmDamage->InsertNextValue(YADE_PTR_CAST<JCFpmState>(b->state)->damageIndex);
}
if (recActive[REC_COORDNUMBER]){
spheresCoordNumb->InsertNextValue(b->coordNumber());
@@ -847,7 +852,8 @@
}
if (recActive[REC_JCFPM]) {
- spheresUg->GetPointData()->AddArray(damage);
+ spheresUg->GetPointData()->AddArray(nbCracks);
+ spheresUg->GetPointData()->AddArray(jcfpmDamage);
}
if (recActive[REC_BSTRESS])
{
@@ -1014,17 +1020,19 @@
if(file){
while ( !file.eof() ){
std::string line;
- Real i,p0,p1,p2,t,s,n0,n1,n2;
+ Real iter,time,p0,p1,p2,type,size,n0,n1,n2,nrg;
while ( std::getline(file, line)) {/* writes into string "line", a line of file "file". To go along diff. lines*/
- file >> i >> p0 >> p1 >> p2 >> t >> s >> n0 >> n1 >> n2;
+ file >> iter >> time >> p0 >> p1 >> p2 >> type >> size >> n0 >> n1 >> n2 >> nrg;
vtkIdType pid[1];
pid[0] = crackPos->InsertNextPoint(p0, p1, p2);
crackCells->InsertNextCell(1,pid);
- crackType->InsertNextValue(t);
- crackSize->InsertNextValue(s);
- iter->InsertNextValue(i);
+ crackIter->InsertNextValue(iter);
+ crackTime->InsertNextValue(time);
+ crackType->InsertNextValue(type);
+ crackSize->InsertNextValue(size);
Real n[3] = { n0, n1, n2 };
crackNorm->InsertNextTupleValue(n);
+ crackNrg->InsertNextValue(nrg);
}
}
file.close();
@@ -1032,10 +1040,12 @@
//
crackUg->SetPoints(crackPos);
crackUg->SetCells(VTK_VERTEX, crackCells);
- crackUg->GetPointData()->AddArray(iter);
+ crackUg->GetPointData()->AddArray(crackIter);
+ crackUg->GetPointData()->AddArray(crackTime);
crackUg->GetPointData()->AddArray(crackType);
crackUg->GetPointData()->AddArray(crackSize);
crackUg->GetPointData()->AddArray(crackNorm); //see https://www.mail-archive.com/paraview@xxxxxxxxxxxx/msg08166.html to obtain Paraview 2D glyphs conforming to this normal
+ crackUg->GetPointData()->AddArray(crackNrg);
vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
if(compress) writer->SetCompressor(compressor);
Follow ups
