yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11466
[Branch ~yade-pkg/yade/git-trunk] Rev 3461: DFNFlow unblock cells as fractures reach them + additional attributes in JCFPM
------------------------------------------------------------
revno: 3461
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
timestamp: Thu 2014-10-09 23:03:29 +0200
message:
DFNFlow unblock cells as fractures reach them + additional attributes in JCFPM
modified:
pkg/dem/JointedCohesiveFrictionalPM.cpp
pkg/dem/JointedCohesiveFrictionalPM.hpp
pkg/pfv/DFNFlow.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-07-18 18:18:50 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.cpp 2014-10-09 21:03:29 +0000
@@ -41,9 +41,18 @@
}
if ( smoothJoint && phys->isOnJoint ) {
- if ( phys->more || ( phys->jointCumulativeSliding > (2*min(geom->radius1,geom->radius2)) ) ) {
- return false;
- } else {
+ 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;
+ }
+ }
+ else {
D = phys->initD - std::abs((b1->state->pos - b2->state->pos).dot(phys->jointNormal));
}
} else {
@@ -51,10 +60,17 @@
}
/* Determination of interaction */
- if (D < 0) { //spheres do not "touch" !
- if ( !phys->isCohesive )
- {
- return false;
+ 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;
+ }
}
if ( phys->isCohesive && (phys->FnMax>0) && (std::abs(D)>Dtensile) ) {
@@ -76,11 +92,22 @@
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;
}
cracksFileExist=true;
-
- // delete contact
- return false;
+ /// Timos
+ if (!neverErase) return false;
+ else
+ {
+ phys->shearForce = Vector3r::Zero();
+ phys->normalForce = Vector3r::Zero();
+ phys->isCohesive =0;
+ phys->FnMax = 0;
+ phys->FsMax = 0;
+ phys->interactionIsCracked = 1;
+ }
+ return true;
}
- }
+ }
+// phys->crackJointAperture = D;
+ phys->crackJointAperture = -D;
/* NormalForce */
Real Fn = 0;
@@ -141,15 +168,18 @@
cracksFileExist=true;
// delete contact if in tension, set the contact properties to friction if in compression
- if ( D < 0 )
- {
- return false;
- }
- else
- {
- phys->FnMax = 0;
- phys->FsMax = 0;
- phys->isCohesive=false;
+ if ( D < 0 ) {
+ if (!neverErase) return false;
+ else
+ {
+ phys->shearForce = Vector3r::Zero();
+ phys->normalForce = Vector3r::Zero();
+ phys->isCohesive =0;
+ phys->FnMax = 0;
+ phys->FsMax = 0;
+ phys->interactionIsCracked=1;
+ return true;
+ }
}
}
}
=== modified file 'pkg/dem/JointedCohesiveFrictionalPM.hpp'
--- pkg/dem/JointedCohesiveFrictionalPM.hpp 2014-07-24 11:35:36 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.hpp 2014-10-09 21:03:29 +0000
@@ -59,7 +59,7 @@
((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`)."))
+ ((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, and normal of the interaction is re-oriented (see also :yref:`Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM`)."))
((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,6 +68,8 @@
((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,interactionIsCracked,0,,"flag for cracked interactions"))
+ ((Real,crackJointAperture,0,,"Relative displacement between 2 spheres (in case of a crack it is equivalent of the crack aperture)"))
,
createIndex();
,
@@ -103,6 +105,7 @@
((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"))
);
DECLARE_LOGGER;
};
=== modified file 'pkg/pfv/DFNFlow.cpp'
--- pkg/pfv/DFNFlow.cpp 2014-06-20 20:10:04 +0000
+++ pkg/pfv/DFNFlow.cpp 2014-10-09 21:03:29 +0000
@@ -7,6 +7,9 @@
*************************************************************************/
#ifdef FLOW_ENGINE
+#include<yade/pkg/dem/JointedCohesiveFrictionalPM.hpp>
+// #include<yade/pkg/dem/ScGeom.hpp>
+
//keep this #ifdef for commited versions unless you really have stable version that should be compiled by default
//it will save compilation time for everyone else
//when you want it compiled, you can pass -DDFNFLOW to cmake, or just uncomment the following line
@@ -19,7 +22,11 @@
{
public:
Real anotherVariable;
- void anotherFunction() {};
+ bool crack;
+// bool preExistingJoint;
+// void anotherFunction() {};
+// DFNCellInfo() : FlowCellInfo(),crack(false) {}
+ DFNCellInfo() : crack(false) {}
};
@@ -28,7 +35,97 @@
//same here if needed
};
-typedef TemplateFlowEngine_DFNFlowEngineT<DFNCellInfo,DFNVertexInfo> DFNFlowEngineT;
+typedef CGT::_Tesselation<CGT::TriangulationTypes<DFNVertexInfo,DFNCellInfo> > DFNTesselation;
+
+#ifdef LINSOLV
+class DFNBoundingSphere : public CGT::FlowBoundingSphereLinSolv<DFNTesselation>
+#else
+class DFNBoundingSphere : public CGT::FlowBoundingSphere<DFNTesselation>
+#endif
+{
+public:
+
+ void saveVtk(const char* folder)
+ {
+// CGT::FlowBoundingSphere<DFNTesselation>::saveVtk(folder)
+ RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
+ static unsigned int number = 0;
+ char filename[80];
+ mkdir(folder, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+ sprintf(filename,"%s/out_%d.vtk",folder,number++);
+ int firstReal=-1;
+
+ //count fictious vertices and cells
+ vtkInfiniteVertices=vtkInfiniteCells=0;
+ FiniteCellsIterator cellEnd = Tri.finite_cells_end();
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
+ bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal() && cell->vertex(3)->info().isReal();
+ if (!isDrawable) vtkInfiniteCells+=1;
+ }
+ for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
+ if (!v->info().isReal()) vtkInfiniteVertices+=1;
+ else if (firstReal==-1) firstReal=vtkInfiniteVertices;}
+
+ basicVTKwritter vtkfile((unsigned int) Tri.number_of_vertices()-vtkInfiniteVertices, (unsigned int) Tri.number_of_finite_cells()-vtkInfiniteCells);
+
+ vtkfile.open(filename,"test");
+
+ vtkfile.begin_vertices();
+ double x,y,z;
+ for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
+ if (v->info().isReal()){
+ x = (double)(v->point().point()[0]);
+ y = (double)(v->point().point()[1]);
+ z = (double)(v->point().point()[2]);
+ vtkfile.write_point(x,y,z);}
+ }
+ vtkfile.end_vertices();
+
+ vtkfile.begin_cells();
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
+ bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal() && cell->vertex(3)->info().isReal();
+ if (isDrawable){vtkfile.write_cell(cell->vertex(0)->info().id()-firstReal, cell->vertex(1)->info().id()-firstReal, cell->vertex(2)->info().id()-firstReal, cell->vertex(3)->info().id()-firstReal);}
+ }
+ vtkfile.end_cells();
+
+ if (permeabilityMap){
+ vtkfile.begin_data("Permeability",CELL_DATA,SCALARS,FLOAT);
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
+ bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal() && cell->vertex(3)->info().isReal();
+ if (isDrawable){vtkfile.write_data(cell->info().s);}
+ }
+ vtkfile.end_data();}
+ else{
+ vtkfile.begin_data("Pressure",CELL_DATA,SCALARS,FLOAT);
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
+ bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal() && cell->vertex(3)->info().isReal();
+ if (isDrawable){vtkfile.write_data(cell->info().p());}
+ }
+ vtkfile.end_data();}
+
+ if (1){
+ averageRelativeCellVelocity();
+ vtkfile.begin_data("Velocity",CELL_DATA,VECTORS,FLOAT);
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
+ bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal() && cell->vertex(3)->info().isReal();
+ if (isDrawable){vtkfile.write_data(cell->info().averageVelocity()[0],cell->info().averageVelocity()[1],cell->info().averageVelocity()[2]);}
+ }
+ vtkfile.end_data();}
+// /// Check this one, cell info()->cracked not defined yet
+ if(1){
+// // trickPermeability();
+ vtkfile.begin_data("fracturedCells",CELL_DATA,SCALARS,FLOAT);
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
+ bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal() && cell->vertex(3)->info().isReal();
+ if (isDrawable){vtkfile.write_data(cell->info().crack);}
+ }
+ vtkfile.end_data();}
+ }
+ // e.g. vtk recorders
+};
+
+
+typedef TemplateFlowEngine_DFNFlowEngineT<DFNCellInfo,DFNVertexInfo, DFNTesselation> DFNFlowEngineT;
REGISTER_SERIALIZABLE(DFNFlowEngineT);
YADE_PLUGIN((DFNFlowEngineT));
@@ -36,15 +133,22 @@
{
public :
void trickPermeability();
- void trickPermeability (RTriangulation::Facet_circulator& facet,Real somethingBig);
- void trickPermeability (RTriangulation::Finite_edges_iterator& edge,Real somethingBig);
+ void trickPermeability (RTriangulation::Facet_circulator& facet,Real aperture, Real residualAperture);
+ void trickPermeability (RTriangulation::Finite_edges_iterator& edge,Real aperture, Real residualAperture);
void setPositionsBuffer(bool current);
+// void computeTotalFractureArea(Real totalFracureArea,bool printFractureTotalArea);/// Trying to get fracture's surface
+ Real totalFracureArea; /// Trying to get fracture's surface
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(DFNFlowEngine,DFNFlowEngineT,"documentation here",
- ((Real, myNewAttribute, 0,,"useless example"))
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(DFNFlowEngine,DFNFlowEngineT,"This is an enhancement of the FlowEngine for intact and fractured rocks that takes into acount pre-existing discontinuities and bond breakage between particles. The local conductivity around the broken link is calculated according to parallel plates model",
+ ((Real, myNewAttribute, 0,,"useless example, Input values to be implemented: newCrackResidualPermeability, jointResidualPermeability, crackOpeningFactor"))
((bool, updatePositions, false,,"update particles positions when rebuilding the mesh (experimental)"))
- ,/*DFNFlowEngineT()*/,
- ,
+ ((Real, jointResidual, 0,,"Calibration parameter for residual aperture of joints"))
+ ((bool, printFractureTotalArea, 0,,"The final fracture area computed through the network")) /// Trying to get fracture's surface
+ ,
+ ,
+ ,
+// .def("computeTotalFractureArea",&DFNFlowEngineT::computeTotalFractureArea," Compute and print the total fracture area of the network") /// Trying to get fracture's surface
+// .def("trickPermeability",&DFNFlowEngineT::trickPermeability,"measure the mean trickPermeability in the period")
)
DECLARE_LOGGER;
};
@@ -70,38 +174,86 @@
}
}
-
-void DFNFlowEngine::trickPermeability (RTriangulation::Facet_circulator& facet, Real somethingBig)
+void DFNFlowEngine::trickPermeability (RTriangulation::Facet_circulator& facet, Real aperture, Real residualAperture)
{
const RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
const CellHandle& cell1 = facet->first;
const CellHandle& cell2 = facet->first->neighbor(facet->second);
if ( Tri.is_infinite(cell1) || Tri.is_infinite(cell2)) cerr<<"Infinite cell found in trickPermeability, should be handled somehow, maybe"<<endl;
- cell1->info().kNorm()[facet->second] = somethingBig;
- cell2->info().kNorm()[Tri.mirror_index(cell1, facet->second)] = somethingBig;
+ cell1->info().kNorm()[facet->second] =pow((aperture+residualAperture),3) / (12 * viscosity);
+ cell2->info().kNorm()[Tri.mirror_index(cell1, facet->second)] =pow((aperture+residualAperture),3) / (12 * viscosity);
+ //For vtk recorder:
+ cell1->info().crack= 1;
+ cell2->info().crack= 1;
+ cell2->info().blocked = cell1->info().blocked = cell2->info().Pcondition = cell1->info().Pcondition = false;//those ones will be included in the flow problem
+ Point& CellCentre1 = cell1->info(); /// Trying to get fracture's surface
+ Point& CellCentre2 = cell2->info(); /// Trying to get fracture's surface
+ CVector networkFractureLength = CellCentre1 - CellCentre2; /// Trying to get fracture's surface
+ double networkFractureDistance = sqrt(networkFractureLength.squared_length()); /// Trying to get fracture's surface
+ Real networkFractureArea = pow(networkFractureDistance,2); /// Trying to get fracture's surface
+ totalFracureArea += networkFractureArea; /// Trying to get fracture's surface
+// cout <<" ------------------ The total surface area up to here is --------------------" << totalFracureArea << endl;
+// printFractureTotalArea = totalFracureArea; /// Trying to get fracture's surface
}
-void DFNFlowEngine::trickPermeability(RTriangulation::Finite_edges_iterator& edge, Real somethingBig)
+void DFNFlowEngine::trickPermeability(RTriangulation::Finite_edges_iterator& edge, Real aperture, Real residualAperture)
{
const RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
RTriangulation::Facet_circulator facet1 = Tri.incident_facets(*edge);
RTriangulation::Facet_circulator facet0=facet1++;
- trickPermeability(facet0,somethingBig);
- while ( facet1!=facet0 ) {trickPermeability(facet1,somethingBig); facet1++;}
+ trickPermeability(facet0, aperture,residualAperture);
+ while ( facet1!=facet0 ) {trickPermeability(facet1, aperture, residualAperture); facet1++;}
+ /// Needs the fracture surface for this edge?
+// double edgeArea = solver->T[solver->currentTes].computeVFacetArea(edge); cout<<"edge area="<<edgeArea<<endl;
}
void DFNFlowEngine::trickPermeability()
{
- Real somethingBig=10;//is that big??
const RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
//We want to change permeability perpendicular to the 10th edge, let's say.
//in the end this function should have a loop on all edges I guess
+ const JCFpmPhys* jcfpmphys;
+ const shared_ptr<InteractionContainer> interactions = scene->interactions;
+ int numberOfCrackedOrJoinedInteractions = 0; /// DEBUG
+ Real SumOfApertures = 0.; /// DEBUG
+ Real AverageAperture =0; /// DEBUG
+ Real totalFracureArea=0; /// Trying to get fracture's surface
+// const shared_ptr<IGeom>& ig;
+// const ScGeom* geom; // = static_cast<ScGeom*>(ig.get());
FiniteEdgesIterator edge = Tri.finite_edges_begin();
- for(int k=10; k>0;--k) edge++;
- trickPermeability(edge,somethingBig);
+ for( ; edge!= Tri.finite_edges_end(); ++edge) {
+ const VertexInfo& vi1=(edge->first)->vertex(edge->second)->info();
+ const VertexInfo& vi2=(edge->first)->vertex(edge->third)->info();
+ const shared_ptr<Interaction>& interaction=interactions->find( vi1.id(),vi2.id() );
+ if (interaction && interaction->isReal()) {
+ numberOfCrackedOrJoinedInteractions +=1; /// DEBUG
+
+ jcfpmphys = YADE_CAST<JCFpmPhys*>(interaction->phys.get());
+// if ( (jcfpmphys->interactionIsCracked || jcfpmphys->isOnJoint) ) {Real aperture=jcfpmphys->dilation; Real residualAperture = jcfpmphys->isOnJoint? jointResidual : 0; trickPermeability(edge,aperture, residualAperture);};
+// if ( (jcfpmphys->interactionIsCracked || jcfpmphys->isOnJoint) ) {Real aperture=(geom->penetrationDepth - jcfpmphys->initD); Real residualAperture = jcfpmphys->isOnJoint? jointResidual : 0; trickPermeability(edge,aperture, residualAperture);};
+ if ( (jcfpmphys->interactionIsCracked || jcfpmphys->isOnJoint) ) {
+ Real residualAperture = jcfpmphys->isOnJoint? jointResidual : 1e-5;
+// Real aperture = jcfpmphys->crackJointAperture;
+// Real aperture = jcfpmphys->crackJointAperture > 0? jcfpmphys->crackJointAperture : 0.0000000001;
+ Real aperture = (jcfpmphys->crackJointAperture + residualAperture) > 1e-4? jcfpmphys->crackJointAperture : 1e-4 ;
+// cout<<"crackJointperture = " << aperture <<endl;
+// cout<<"initial Aperture = " << residualAperture <<endl;
+ SumOfApertures += aperture; /// DEBUG
+ trickPermeability(edge,aperture, residualAperture);};
+ // TEST LUC PULL COMMIT
+ }
+ }
+ AverageAperture = SumOfApertures/numberOfCrackedOrJoinedInteractions; /// DEBUG
+// cout << " Average aperture in joint ( -D ) = " << AverageAperture << endl; /// DEBUG
}
+// Real DFNFlowEngine::computeTotalFractureArea(totalFracureArea,printFractureTotalArea) /// Trying to get fracture's surface
+// {
+// if (printFractureTotalArea >0) {
+// cout<< " The total fracture area computed from the Network is: " << totalFracureArea <<endl;
+// }
+// }
#endif //DFNFLOW
#endif //FLOWENGINE
Follow ups
