← Back to team overview

yade-dev team mailing list archive

[Branch ~yade-pkg/yade/git-trunk] Rev 3699: Merge branch 'master' of github.com:yade/trunk

 

Merge authors:
  Christian Jakob (jakob-ifgt)
  jduriez (jduriez)
------------------------------------------------------------
revno: 3699 [merge]
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Fri 2013-10-04 17:30:32 +0200
message:
  Merge branch 'master' of github.com:yade/trunk
  
  Conflicts:
  	pkg/dem/VTKRecorder.cpp
  	pkg/dem/VTKRecorder.hpp
modified:
  doc/references.bib
  doc/sphinx/references.bib
  doc/sphinx/user.rst
  examples/clumps/replaceByClumps-example.py
  pkg/dem/JointedCohesiveFrictionalPM.cpp
  pkg/dem/JointedCohesiveFrictionalPM.hpp
  pkg/dem/VTKRecorder.cpp
  pkg/dem/VTKRecorder.hpp
  py/wrapper/yadeWrapper.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 'doc/references.bib'
--- doc/references.bib	2013-09-23 17:17:34 +0000
+++ doc/references.bib	2013-10-04 12:35:58 +0000
@@ -619,3 +619,14 @@
 	url={http://stacks.iop.org/0295-5075/94/i=5/a=50004},
 	year={2011}
 }
+
+@article{Ivars2011,
+title = "The synthetic rock mass approach for jointed rock mass modelling ",
+journal = "International Journal of Rock Mechanics and Mining Sciences ",
+volume = "48",
+number = "2",
+pages = "219 - 244",
+year = "2011",
+doi = "10.1016/j.ijrmms.2010.11.014",
+author = "Diego Mas Ivars and Matthew E. Pierce and Caroline Darcel and Juan Reyes-Montes and David O. Potyondy and R. Paul Young and Peter A. Cundall"
+}

=== modified file 'doc/sphinx/references.bib'
--- doc/sphinx/references.bib	2013-08-28 10:38:02 +0000
+++ doc/sphinx/references.bib	2013-10-04 12:35:58 +0000
@@ -79,6 +79,17 @@
 	url={ http://geo.hmg.inpg.fr/frederic/Discrete_Element_Group_FVD.html }
 }
 
+@article{Ivars2011,
+title = "The synthetic rock mass approach for jointed rock mass modelling ",
+journal = "International Journal of Rock Mechanics and Mining Sciences ",
+volume = "48",
+number = "2",
+pages = "219 - 244",
+year = "2011",
+doi = "10.1016/j.ijrmms.2010.11.014",
+author = "Diego Mas Ivars and Matthew E. Pierce and Caroline Darcel and Juan Reyes-Montes and David O. Potyondy and R. Paul Young and Peter A. Cundall"
+}
+
 @article{Jerier2009,
 	journal={Granular Matter},
 	title={A geometric algorithm based on tetrahedral meshes to generate a dense polydisperse sphere packing},

=== modified file 'doc/sphinx/user.rst'
--- doc/sphinx/user.rst	2013-10-02 15:32:00 +0000
+++ doc/sphinx/user.rst	2013-10-04 13:41:43 +0000
@@ -173,7 +173,7 @@
 
 	Yade [4]: O.bodies.clump(bodyList)
 
-	Yade [5]: RC=O.bodies.getRoundness([]) # empty list [] is needed if no body should be excluded
+	Yade [5]: RC=O.bodies.getRoundness()
 
 	Yade [3]: print RC
 	
@@ -268,7 +268,7 @@
 #. :yref:`yade.pack.gtsSurface2Facets` function can create the triangulated surface (from :yref:`Facet` particles) in the simulation itself, as shown in the funnel example. (Triangulated surface can also be imported directly from a STL file using :yref:`yade.ymport.stl`.)
 #. :yref:`yade._packPredicates.inGtsSurface` predicate can be created, using the surface as boundary representation of the enclosed volume.
 
-The :ysrc:`scripts/gts-horse/gts-horse.py` (img. img-horse_) shows both possibilities; first, a GTS surface is imported::
+The :ysrc:`examples/gts-horse/gts-horse.py` (img. img-horse_) shows both possibilities; first, a GTS surface is imported::
 
 	import gts
 	surf=gts.read(open('horse.coarse.gts'))
@@ -495,7 +495,7 @@
 Individual interactions on demand
 ----------------------------------
 
-It is possible to create an interaction between a pair of particles independently of collision detection using :yref:`createInteraction<yade.utils.createInteraction>`. This function looks for and uses matching ``Ig2`` and ``Ip2`` functors. Interaction will be created regardless of distance between given particles (by passing a special parameter to the ``Ig2`` functor to force creation of the interaction even without any geometrical contact). Appropriate constitutive law should be used to avoid deletion of the interaction at the next simulation step.
+It is possible to create an interaction between a pair of particles independently of collision detection using :yref:`createInteraction<yade._utils.createInteraction>`. This function looks for and uses matching ``Ig2`` and ``Ip2`` functors. Interaction will be created regardless of distance between given particles (by passing a special parameter to the ``Ig2`` functor to force creation of the interaction even without any geometrical contact). Appropriate constitutive law should be used to avoid deletion of the interaction at the next simulation step.
 
 .. ipython::
 
@@ -523,7 +523,7 @@
 	# created by functors in InteractionLoop
 	Yade [2]: i.geom, i.phys
 
-This method will be rather slow if many interaction are to be created (the functor lookup will be repeated for each of them). In such case, ask on yade-dev@xxxxxxxxxxxxxxxxxxx to have the :yref:`createInteraction<yade.utils.createInteraction>` function accept list of pairs id's as well.
+This method will be rather slow if many interaction are to be created (the functor lookup will be repeated for each of them). In such case, ask on yade-dev@xxxxxxxxxxxxxxxxxxx to have the :yref:`createInteraction<yade._utils.createInteraction>` function accept list of pairs id's as well.
 
 Base engines
 =============
@@ -919,7 +919,7 @@
 
 A special engine :yref:`PyRunner` can be used to periodically call python code, specified via the ``command`` parameter. Periodicity can be controlled by specifying computation time (``realPeriod``), virutal time (``virtPeriod``) or iteration number (``iterPeriod``).
 
-For instance, to print kinetic energy (using :yref:`kineticEnergy<yade.utils.kineticEnergy>`) every 5 seconds, this engine will be put to ``O.engines``::
+For instance, to print kinetic energy (using :yref:`kineticEnergy<yade._utils.kineticEnergy>`) every 5 seconds, this engine will be put to ``O.engines``::
 	PyRunner(command="print 'kinetic energy',kineticEnergy()",realPeriod=5)
 
 For running more complex commands, it is convenient to define an external function and only call it from within the engine. Since the ``command`` is run in the script's namespace, functions defined within scripts can be called. Let us print information on interaction between bodies 0 and 1 periodically::
@@ -1138,7 +1138,7 @@
 
 Frequently, decisions have to be made based on evolution of the simulation itself, which is not yet known. In such case, a function checking some specific condition is called periodically; if the condition is satisfied, ``O.pause`` or other functions can be called to stop the stimulation. See documentation for :yref:`Omega.run`, :yref:`Omega.pause`, :yref:`Omega.step`, :yref:`Omega.stopAtIter` for details.
 
-For simulations that seek static equilibrium, the :yref:`unbalancedForce<yade.utils.unbalancedForce>` can provide a useful metrics (see its documentation for details); for a desired value of ``1e-2`` or less, for instance, we can use::
+For simulations that seek static equilibrium, the :yref:`unbalancedForce<yade._utils.unbalancedForce>` can provide a useful metrics (see its documentation for details); for a desired value of ``1e-2`` or less, for instance, we can use::
 
 	
 	def checkUnbalanced():
@@ -1555,7 +1555,7 @@
 
 Micro-stress
 ------------
-Stress fields can be generated by combining the volume returned by TesselationWrapper to per-particle stress given by :yref:`bodyStressTensors<yade.utils.bodyStressTensors>`. Since the stress $\sigma$ from bodyStressTensor implies a division by the volume $V_b$ of the solid particle, one has to re-normalize it in order to obtain the micro-stress as defined in [Catalano2013a]_ (equation 39 therein), i.e. $\overline{\sigma}^k = \sigma^k \times V_b^k / V_{\sigma}^k$ where $V_{\sigma}^k$ is the volume assigned to particle $k$ in the tesselation. For instance:
+Stress fields can be generated by combining the volume returned by TesselationWrapper to per-particle stress given by :yref:`bodyStressTensors<yade._utils.bodyStressTensors>`. Since the stress $\sigma$ from bodyStressTensor implies a division by the volume $V_b$ of the solid particle, one has to re-normalize it in order to obtain the micro-stress as defined in [Catalano2013a]_ (equation 39 therein), i.e. $\overline{\sigma}^k = \sigma^k \times V_b^k / V_{\sigma}^k$ where $V_{\sigma}^k$ is the volume assigned to particle $k$ in the tesselation. For instance:
 
 .. code-block:: python
 

=== modified file 'examples/clumps/replaceByClumps-example.py'
--- examples/clumps/replaceByClumps-example.py	2013-04-16 12:37:01 +0000
+++ examples/clumps/replaceByClumps-example.py	2013-10-04 13:41:43 +0000
@@ -33,7 +33,7 @@
 O.bodies.append([sphere(c,r,material=Mat) for c,r in sp])
 
 print len(sp),' particles generated.'
-print 'Roundness coefficient without clumps is: ',O.bodies.getRoundness([])			#give an empty list [] if no body should be excluded
+print 'Roundness coefficient without clumps is: ',O.bodies.getRoundness()
 
 
 #### show how to use makeClumpTemplate():
@@ -73,7 +73,7 @@
 	if b.isStandalone:
 		standaloneList.append(b.id)
 
-print 'Roundness coefficient for spheres and clumps is: ',O.bodies.getRoundness([])			#give an empty list [] if no body should be excluded
+print 'Roundness coefficient for spheres and clumps is: ',O.bodies.getRoundness()
 print 'Roundness coefficient just for clumps is: ',O.bodies.getRoundness(standaloneList)
 
 O.dt=1e-6

=== modified file 'pkg/dem/JointedCohesiveFrictionalPM.cpp'
--- pkg/dem/JointedCohesiveFrictionalPM.cpp	2013-09-19 08:50:04 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.cpp	2013-10-04 12:35:58 +0000
@@ -60,6 +60,8 @@
 	    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;
 	    
     	    // create a text file to record properties of the broken bond (iteration, position, type (tensile), cross section and contact normal orientation)
 	    if (recordCracks){
@@ -118,6 +120,8 @@
 	    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)
 	    if (recordCracks){
@@ -195,7 +199,13 @@
 	
 	// cohesive properties
 	///to set if the contact is cohesive or not
-	if ((scene->iter < cohesiveTresholdIteration) && (tensileStrength>0 || cohesion>0) && (yade1->type == yade2->type)){ contactPhysics->isCohesive=true; }
+	if ((scene->iter < cohesiveTresholdIteration) && (tensileStrength>0 || cohesion>0) && (yade1->type == yade2->type)){ 
+	  contactPhysics->isCohesive=true;
+// 	  JCFpmState* st1=dynamic_cast<JCFpmState*>(Body::byId(interaction->getId1(),scene)->state.get());
+// 	  st1->noIniLinks++;
+// 	  JCFpmState* st2=dynamic_cast<JCFpmState*>(Body::byId(interaction->getId2(),scene)->state.get());
+// 	  st2->noIniLinks++;
+	}
 	
 	if ( contactPhysics->isCohesive ) {
 	  contactPhysics->FnMax = tensileStrength*contactPhysics->crossSection;
@@ -221,7 +231,13 @@
 			contactPhysics->tanDilationAngle = std::tan(contactPhysics->dilationAngle);
 		  
 			///to set if the contact is cohesive or not
-			if ((scene->iter < cohesiveTresholdIteration) && (jointCohesion>0 || jointTensileStrength>0)) { contactPhysics->isCohesive=true; } 
+			if ((scene->iter < cohesiveTresholdIteration) && (jointCohesion>0 || jointTensileStrength>0)) {
+			  contactPhysics->isCohesive=true;
+// 			  JCFpmState* st1=dynamic_cast<JCFpmState*>(Body::byId(interaction->getId1(),scene)->state.get());
+// 			  st1->noIniLinks++;
+// 			  JCFpmState* st2=dynamic_cast<JCFpmState*>(Body::byId(interaction->getId2(),scene)->state.get());
+// 			  st2->noIniLinks++;
+			} 
 			else { contactPhysics->isCohesive=false; contactPhysics->FnMax=0; contactPhysics->FsMax=0; }
 		  
 			if ( contactPhysics->isCohesive ) {

=== modified file 'pkg/dem/JointedCohesiveFrictionalPM.hpp'
--- pkg/dem/JointedCohesiveFrictionalPM.hpp	2013-09-26 09:10:15 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.hpp	2013-10-04 12:35:58 +0000
@@ -100,8 +100,8 @@
 		virtual void 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 jointed frictional material. Basically, this law adds the possibility to define joint surfaces into a cohesive frictional material as defined by :yref:`Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM`. Joint surfaces can be defined in a preprocessing phase through .stl meshes (see ref for details of the procedure).",
-			((bool,smoothJoint,false,,"if true, particles belonging to joint surface have a smooth contact logic."))
+		YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM,LawFunctor,"Interaction law for cohesive frictional material, e.g. rock, possibly presenting joint surfaces. Joint surfaces can be defined in a preprocessing phase through .stl meshes (see ref for details of the procedure), and can be mechanically described with a smooth contact logic [Ivars2011]_ (implemented in Yade in [Scholtes2012]_).",
+			((bool,smoothJoint,false,,"if true, particles belonging to joint surface have a smooth contact logic [Ivars2011]_, [Scholtes2012]_."))
 			((bool,recordCracks,false,,"if true a text file cracks.txt will be created (iteration, position, type (tensile or shear), cross section and contact normal)."))
 			((bool,cracksFileExist,false,,"If true, text file already exists and its content won't be reset."))
 		);

=== modified file 'pkg/dem/VTKRecorder.cpp'
--- pkg/dem/VTKRecorder.cpp	2013-10-04 15:25:30 +0000
+++ pkg/dem/VTKRecorder.cpp	2013-10-04 15:30:32 +0000
@@ -28,6 +28,7 @@
 #include<yade/pkg/common/Box.hpp>
 #include<yade/pkg/dem/ConcretePM.hpp>
 #include<yade/pkg/dem/WirePM.hpp>
+#include<yade/pkg/dem/JointedCohesiveFrictionalPM.hpp>
 #include<yade/pkg/dem/Shop.hpp>
 
 
@@ -66,12 +67,17 @@
 		else if((rec=="clumpids") || (rec=="clumpId")) recActive[REC_CLUMPID]=true;
 		else if(rec=="materialId") recActive[REC_MATERIALID]=true;
 		else if(rec=="stress") recActive[REC_STRESS]=true;
+		else if(rec=="jcfpm") recActive[REC_JCFPM]=true;
+		else if(rec=="cracks") recActive[REC_CRACKS]=true;
 		else if(rec=="pericell" && scene->isPeriodic) recActive[REC_PERICELL]=true;
-		else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: all, spheres, velocity, facets, boxes, color, stress, cpm, wpm, intr, id, clumpId, materialId, pericell). Ignored.");
+		else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: all, spheres, velocity, facets, boxes, color, stress, cpm, wpm, intr, id, clumpId, materialId, jcfpm, cracks, pericell). Ignored.");
 	}
 	// cpm needs interactions
 	if(recActive[REC_CPM]) recActive[REC_INTR]=true;
 	
+	// jcfpm needs interactions
+	if(recActive[REC_JCFPM]) recActive[REC_INTR]=true;
+
 	// wpm needs interactions
 	if(recActive[REC_WPM]) recActive[REC_INTR]=true;
 	
@@ -205,6 +211,52 @@
 	cpmStress->SetNumberOfComponents(9);
 	cpmStress->SetName("cpmStress");
 
+	// extras for JCFpm
+	vtkSmartPointer<vtkFloatArray> damage = vtkSmartPointer<vtkFloatArray>::New();
+	damage->SetNumberOfComponents(1);;
+	damage->SetName("damage");
+	vtkSmartPointer<vtkFloatArray> damageRel = vtkSmartPointer<vtkFloatArray>::New();
+	damageRel->SetNumberOfComponents(1);;
+	damageRel->SetName("damageRel");
+	vtkSmartPointer<vtkFloatArray> intrIsCohesive = vtkSmartPointer<vtkFloatArray>::New();
+	intrIsCohesive->SetNumberOfComponents(1);
+	intrIsCohesive->SetName("isCohesive");
+	vtkSmartPointer<vtkFloatArray> intrIsOnJoint = vtkSmartPointer<vtkFloatArray>::New();
+	intrIsOnJoint->SetNumberOfComponents(1);
+	intrIsOnJoint->SetName("isOnJoint");
+	
+	// extras for cracks
+	vtkSmartPointer<vtkPoints> crackPos = vtkSmartPointer<vtkPoints>::New();
+	vtkSmartPointer<vtkCellArray> crackCells = vtkSmartPointer<vtkCellArray>::New();
+	vtkSmartPointer<vtkFloatArray> iter = vtkSmartPointer<vtkFloatArray>::New();
+	iter->SetNumberOfComponents(1);
+	iter->SetName("iter");
+	vtkSmartPointer<vtkFloatArray> crackType = vtkSmartPointer<vtkFloatArray>::New();
+	crackType->SetNumberOfComponents(1);
+	crackType->SetName("crackType");
+	vtkSmartPointer<vtkFloatArray> crackSize = vtkSmartPointer<vtkFloatArray>::New();
+	crackSize->SetNumberOfComponents(1);
+	crackSize->SetName("crackSize");
+	vtkSmartPointer<vtkFloatArray> crackNorm = vtkSmartPointer<vtkFloatArray>::New();
+	crackNorm->SetNumberOfComponents(3);
+	crackNorm->SetName("crackNorm");
+
+// 	// the same for newly created cracks
+// 	vtkSmartPointer<vtkPoints> crackPosNew = vtkSmartPointer<vtkPoints>::New();
+// 	vtkSmartPointer<vtkCellArray> crackCellsNew = vtkSmartPointer<vtkCellArray>::New();
+// 	vtkSmartPointer<vtkFloatArray> iterNew = vtkSmartPointer<vtkFloatArray>::New();
+// 	iterNew->SetNumberOfComponents(1);
+// 	iterNew->SetName("iter");
+// 	vtkSmartPointer<vtkFloatArray> crackTypeNew = vtkSmartPointer<vtkFloatArray>::New();
+// 	crackTypeNew->SetNumberOfComponents(1);
+// 	crackTypeNew->SetName("crackType");
+// 	vtkSmartPointer<vtkFloatArray> crackSizeNew = vtkSmartPointer<vtkFloatArray>::New();
+// 	crackSizeNew->SetNumberOfComponents(1);
+// 	crackSizeNew->SetName("crackSize");
+// 	vtkSmartPointer<vtkFloatArray> crackNormNew = vtkSmartPointer<vtkFloatArray>::New();
+// 	crackNormNew->SetNumberOfComponents(3);
+// 	crackNormNew->SetName("crackNorm");
+	
 	// extras for WireMatPM
 	vtkSmartPointer<vtkFloatArray> wpmNormalForce = vtkSmartPointer<vtkFloatArray>::New();
 	wpmNormalForce->SetNumberOfComponents(1);
@@ -290,6 +342,13 @@
 						wpmLimitFactor->InsertNextValue(NaN);
 					}
 				}
+				else if (recActive[REC_JCFPM]){
+					const JCFpmPhys* jcfpmphys = YADE_CAST<JCFpmPhys*>(I->phys.get());
+					intrIsCohesive->InsertNextValue(jcfpmphys->isCohesive);
+					intrIsOnJoint->InsertNextValue(jcfpmphys->isOnJoint);
+					intrForceN->InsertNextValue(fn);
+				}
+
 				else {
 					intrForceN->InsertNextValue(fn);
 				}
@@ -351,6 +410,11 @@
 					cpmStress->InsertNextTupleValue(s);
 				}
 				
+				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);
+				}
+
 				if (recActive[REC_MATERIALID]) spheresMaterialId->InsertNextValue(b->material->id);
 				continue;
 			}
@@ -496,6 +560,9 @@
 			spheresUg->GetPointData()->AddArray(cpmStress);
 		}
 
+		if (recActive[REC_JCFPM]) {
+			spheresUg->GetPointData()->AddArray(damage);
+		}
 		if (recActive[REC_MATERIALID]) spheresUg->GetPointData()->AddArray(spheresMaterialId);
 
 		#ifdef YADE_VTK_MULTIBLOCK
@@ -565,6 +632,10 @@
 		intrPd->SetLines(intrCells);
 		intrPd->GetCellData()->AddArray(intrForceN);
 		intrPd->GetCellData()->AddArray(intrAbsForceT);
+		if (recActive[REC_JCFPM]) { 
+			intrPd->GetCellData()->AddArray(intrIsCohesive);
+			intrPd->GetCellData()->AddArray(intrIsOnJoint);
+		}
 		if (recActive[REC_WPM]){
 			intrPd->GetCellData()->AddArray(wpmNormalForce);
 			intrPd->GetCellData()->AddArray(wpmLimitFactor);
@@ -599,6 +670,70 @@
 			writer->Write();
 		}
 	}
+
+	if (recActive[REC_CRACKS]) {
+		std::ifstream file ("cracks.txt",std::ios::in);
+		vtkSmartPointer<vtkUnstructuredGrid> crackUg = vtkSmartPointer<vtkUnstructuredGrid>::New();
+		vtkSmartPointer<vtkUnstructuredGrid> crackUgNew = vtkSmartPointer<vtkUnstructuredGrid>::New();
+		
+		 if(file){
+			 while ( !file.eof() ){
+				std::string line;
+				Real i,p0,p1,p2,t,s,n0,n1,n2;
+				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;
+					vtkIdType pid[1];
+					pid[0] = crackPos->InsertNextPoint(p0, p1, p2);
+					crackCells->InsertNextCell(1,pid);
+					crackType->InsertNextValue(t);
+					crackSize->InsertNextValue(s);
+					iter->InsertNextValue(i);
+					float n[3] = { n0, n1, n2 };
+					crackNorm->InsertNextTupleValue(n);
+// 					if (i > scene->iter - iterPeriod)
+// 					{
+// 					  pid[0] = crackPosNew->InsertNextPoint(p0, p1, p2);
+// 					  crackCellsNew->InsertNextCell(1,pid);
+// 					  crackTypeNew->InsertNextValue(t);
+// 					  crackSizeNew->InsertNextValue(s);
+// 					  iterNew->InsertNextValue(i);
+// 					  crackNormNew->InsertNextTupleValue(n);
+// 					}
+					  
+				}
+			 }
+			 file.close();
+		} 
+
+		crackUg->SetPoints(crackPos);
+		crackUg->SetCells(VTK_VERTEX, crackCells);
+		crackUg->GetPointData()->AddArray(iter);
+		crackUg->GetPointData()->AddArray(crackType);
+		crackUg->GetPointData()->AddArray(crackSize);
+		crackUg->GetPointData()->AddArray(crackNorm); //orientation of 2D glyphs does not match this direction (some work to do in order to have the good orientation) 
+		
+		/*crackUgNew->SetPoints(crackPosNew);
+		crackUgNew->SetCells(VTK_VERTEX, crackCellsNew);
+		crackUgNew->GetPointData()->AddArray(iterNew);
+		crackUgNew->GetPointData()->AddArray(crackTypeNew);
+		crackUgNew->GetPointData()->AddArray(crackSizeNew);
+		crackUgNew->GetPointData()->AddArray(crackNormNew); //same remark...*/
+	
+		vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
+		if(compress) writer->SetCompressor(compressor);
+		if(ascii) writer->SetDataModeToAscii();
+		string fn=fileName+"cracks."+lexical_cast<string>(scene->iter)+".vtu";
+		writer->SetFileName(fn.c_str());
+		writer->SetInput(crackUg);
+		writer->Write();
+		
+// 		fn=fileName+"newcracks."+lexical_cast<string>(scene->iter)+".vtu";
+// 		writer->SetFileName(fn.c_str());
+// 		writer->SetInput(crackUgNew);
+// 		writer->Write();
+	}
+
 	#ifdef YADE_VTK_MULTIBLOCK
 		if(multiblock){
 			vtkSmartPointer<vtkMultiBlockDataSet> multiblockDataset = vtkSmartPointer<vtkMultiBlockDataSet>::New();
@@ -615,8 +750,6 @@
 			writer->Write();	
 		}
 	#endif
-
-
 };
 
 void VTKRecorder::addWallVTK (vtkSmartPointer<vtkQuad>& boxes, vtkSmartPointer<vtkPoints>& boxesPos, Vector3r& W1, Vector3r& W2, Vector3r& W3, Vector3r& W4) {

=== modified file 'pkg/dem/VTKRecorder.hpp'
--- pkg/dem/VTKRecorder.hpp	2013-10-04 15:25:30 +0000
+++ pkg/dem/VTKRecorder.hpp	2013-10-04 15:30:32 +0000
@@ -11,7 +11,7 @@
 
 class VTKRecorder: public PeriodicEngine {
 	public:
-		enum {REC_SPHERES=0,REC_FACETS,REC_BOXES,REC_COLORS,REC_MASS,REC_CPM,REC_INTR,REC_VELOCITY,REC_ID,REC_CLUMPID,REC_SENTINEL,REC_MATERIALID,REC_STRESS,REC_MASK,REC_RPM,REC_WPM,REC_PERICELL};
+  enum {REC_SPHERES=0,REC_FACETS,REC_BOXES,REC_COLORS,REC_MASS,REC_CPM,REC_INTR,REC_VELOCITY,REC_ID,REC_CLUMPID,REC_SENTINEL,REC_MATERIALID,REC_STRESS,REC_MASK,REC_RPM,REC_JCFPM,REC_CRACKS,REC_WPM,REC_PERICELL};
 		virtual void action();
 		void addWallVTK (vtkSmartPointer<vtkQuad>& boxes, vtkSmartPointer<vtkPoints>& boxesPos, Vector3r& W1, Vector3r& W2, Vector3r& W3, Vector3r& W4);
 	YADE_CLASS_BASE_DOC_ATTRS_CTOR(VTKRecorder,PeriodicEngine,"Engine recording snapshots of simulation into series of \\*.vtu files, readable by VTK-based postprocessing programs such as Paraview. Both bodies (spheres and facets) and interactions can be recorded, with various vector/scalar quantities that are defined on them.\n\n:yref:`PeriodicEngine.initRun` is initialized to ``True`` automatically.",
@@ -23,7 +23,7 @@
 			((bool,multiblock,false,,"Use multi-block (``.vtm``) files to store data, rather than separate ``.vtu`` files."))
 		#endif
 		((string,fileName,"",,"Base file name; it will be appended with {spheres,intrs,facets}-243100.vtu (unless *multiblock* is ``True``) depending on active recorders and step number (243100 in this case). It can contain slashes, but the directory must exist already."))
-		((vector<string>,recorders,vector<string>(1,string("all")),,"List of active recorders (as strings). ``all`` (the default value) enables all base and generic recorders.\n\n.. admonition:: Base recorders\n\n\tBase recorders save the geometry (unstructured grids) on which other data is defined. They are implicitly activated by many of the other recorders. Each of them creates a new file (or a block, if :yref:`multiblock <VTKRecorder.multiblock>` is set).\n\n\t``spheres``\n\t\tSaves positions and radii (``radii``) of :yref:`spherical<Sphere>` particles.\n\t``facets``\n\t\tSave :yref:`facets<Facet>` positions (vertices).\n\t``boxes``\n\t\tSave :yref:`boxes<Box>` positions (edges).\n\t``intr``\n\t\tStore interactions as lines between nodes at respective particles positions. Additionally stores magnitude of normal (``forceN``) and shear (``absForceT``) forces on interactions (the :yref:`geom<Interaction.geom> must be of type :yref:`NormShearPhys`). \n\n.. admonition:: Generic recorders\n\n\tGeneric recorders do not depend on specific model being used and save commonly useful data.\n\n\t``id``\n\t\tSaves id's (field ``id``) of spheres; active only if ``spheres`` is active.\n\t``mass``\n\t\tSaves masses (field ``mass``) of spheres; active only if ``spheres`` is active.\n\t``clumpId``\n\t\tSaves id's of clumps to which each sphere belongs (field ``clumpId``); active only if ``spheres`` is active.\n\t``colors``\n\t\tSaves colors of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``color``); only active if ``spheres`` or ``facets`` are activated.\n\t``mask``\n\t\tSaves groupMasks of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``mask``); only active if ``spheres`` or ``facets`` are activated.\n\t``materialId``\n\t\tSaves materialID of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>`; only active if ``spheres`` or ``facets`` are activated.\n\t``velocity``\n\t\tSaves linear and angular velocities of spherical particles as Vector3 and length(fields ``linVelVec``, ``linVelLen`` and ``angVelVec``, ``angVelLen`` respectively``); only effective with ``spheres``.\n\t``stress``\n\t\tSaves stresses of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>`  as Vector3 and length; only active if ``spheres`` or ``facets`` are activated.\n\n.. admonition:: Specific recorders\n\n\tThe following should only be activated in appropriate cases, otherwise crashes can occur due to violation of type presuppositions.\n\n\t``cpm``\n\t\tSaves data pertaining to the :yref:`concrete model<Law2_ScGeom_CpmPhys_Cpm>`: ``cpmDamage`` (normalized residual strength averaged on particle), ``cpmStress`` (stress on particle); ``intr`` is activated automatically by ``cpm``\n\t``wpm``\n\t\tSaves data pertaining to the :yref:`wire particle model<Law2_ScGeom_WirePhys_WirePM>`: ``wpmForceNFactor`` shows the loading factor for the wire, e.g. normal force divided by threshold normal force.\n\n"))
+		((vector<string>,recorders,vector<string>(1,string("all")),,"List of active recorders (as strings). ``all`` (the default value) enables all base and generic recorders.\n\n.. admonition:: Base recorders\n\n\tBase recorders save the geometry (unstructured grids) on which other data is defined. They are implicitly activated by many of the other recorders. Each of them creates a new file (or a block, if :yref:`multiblock <VTKRecorder.multiblock>` is set).\n\n\t``spheres``\n\t\tSaves positions and radii (``radii``) of :yref:`spherical<Sphere>` particles.\n\t``facets``\n\t\tSave :yref:`facets<Facet>` positions (vertices).\n\t``boxes``\n\t\tSave :yref:`boxes<Box>` positions (edges).\n\t``intr``\n\t\tStore interactions as lines between nodes at respective particles positions. Additionally stores magnitude of normal (``forceN``) and shear (``absForceT``) forces on interactions (the :yref:`geom<Interaction.geom> must be of type :yref:`NormShearPhys`). \n\n.. admonition:: Generic recorders\n\n\tGeneric recorders do not depend on specific model being used and save commonly useful data.\n\n\t``id``\n\t\tSaves id's (field ``id``) of spheres; active only if ``spheres`` is active.\n\t``mass``\n\t\tSaves masses (field ``mass``) of spheres; active only if ``spheres`` is active.\n\t``clumpId``\n\t\tSaves id's of clumps to which each sphere belongs (field ``clumpId``); active only if ``spheres`` is active.\n\t``colors``\n\t\tSaves colors of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``color``); only active if ``spheres`` or ``facets`` are activated.\n\t``mask``\n\t\tSaves groupMasks of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``mask``); only active if ``spheres`` or ``facets`` are activated.\n\t``materialId``\n\t\tSaves materialID of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>`; only active if ``spheres`` or ``facets`` are activated.\n\t``velocity``\n\t\tSaves linear and angular velocities of spherical particles as Vector3 and length(fields ``linVelVec``, ``linVelLen`` and ``angVelVec``, ``angVelLen`` respectively``); only effective with ``spheres``.\n\t``stress``\n\t\tSaves stresses of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>`  as Vector3 and length; only active if ``spheres`` or ``facets`` are activated.\n\n.. admonition:: Specific recorders\n\n\tThe following should only be activated in appropriate cases, otherwise crashes can occur due to violation of type presuppositions.\n\n\t``cpm``\n\t\tSaves data pertaining to the :yref:`concrete model<Law2_ScGeom_CpmPhys_Cpm>`: ``cpmDamage`` (normalized residual strength averaged on particle), ``cpmStress`` (stress on particle); ``intr`` is activated automatically by ``cpm``\n\t``wpm``\n\t\tSaves data pertaining to the :yref:`wire particle model<Law2_ScGeom_WirePhys_WirePM>`: ``wpmForceNFactor`` shows the loading factor for the wire, e.g. normal force divided by threshold normal force.\n\t``jcfpm``\n\t\tSaves data pertaining to the :yref:`rock (smooth)-jointed model<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM>`: ``damage`` is defined by :yref:`JCFpmState.tensBreak` + :yref:`JCFpmState.shearBreak`; ``intr`` is activated automatically by ``jcfpm``, and :yref:`on joint<JCFpmPhys.isOnJoint>` or :yref:`cohesive<JCFpmPhys.isCohesive>` interactions can be vizualized.\n\t``cracks``\n\t\tSaves other data pertaining to the :yref:`rock model<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM>`: ``cracks`` shows locations where cohesive bonds failed during the simulation, with their types (0/1  for tensile/shear breakages), their sizes (0.5*(R1+R2)), and their normal directions.\n\n"))
 		((int,mask,0,,"If mask defined, only bodies with corresponding groupMask will be exported. If 0, all bodies will be exported.")),
 		/*ctor*/
 		initRun=true;

=== modified file 'py/wrapper/yadeWrapper.cpp'
--- py/wrapper/yadeWrapper.cpp	2013-10-04 11:47:37 +0000
+++ py/wrapper/yadeWrapper.cpp	2013-10-04 13:41:43 +0000
@@ -865,7 +865,7 @@
 		.def("addToClump",&pyBodyContainer::addToClump,(python::arg("discretization")=15,python::arg("integrateInertia")=true),"Add body b (or a list of bodies) to an existing clump c. c must be clump and b may not be a clump member of c. Clump masses and inertia are adapted automatically (for details see :yref:`clump()<BodyContainer.clump>`).\n\nSee :ysrc:`examples/clumps/addToClump-example.py` for an example script.\n\n.. note:: If b is a clump itself, then all members will be added to c and b will be deleted. If b is a clump member of clump d, then all members from d will be added to c and d will be deleted. If you need to add just clump member b, :yref:`release<BodyContainer.releaseFromClump>` this member from d first.")
 		.def("releaseFromClump",&pyBodyContainer::releaseFromClump,(python::arg("discretization")=15,python::arg("integrateInertia")=true),"Release body b from clump c. b must be a clump member of c. Clump masses and inertia are adapted automatically (for details see :yref:`clump()<BodyContainer.clump>`).\n\nSee :ysrc:`examples/clumps/releaseFromClump-example.py` for an example script.\n\n.. note:: If c contains only 2 members b will not be released and a warning will appear. In this case clump c should be :yref:`erased<BodyContainer.erase>`.")
 		.def("replaceByClumps",&pyBodyContainer::replaceByClumps,(python::arg("discretization")=15,python::arg("integrateInertia")=true),"Replace spheres by clumps using a list of clump templates and a list of amounts; returns a list of tuples: ``[(clumpId1,[memberId1,memberId2,...]),(clumpId2,[memberId1,memberId2,...]),...]``. A new clump will have the same volume as the sphere, that was replaced. Clump masses and inertia are adapted automatically (for details see :yref:`clump()<BodyContainer.clump>`). \n\n\t *O.bodies.replaceByClumps( [utils.clumpTemplate([1,1],[.5,.5])] , [.9] ) #will replace 90 % of all standalone spheres by 'dyads'*\n\nSee :ysrc:`examples/clumps/replaceByClumps-example.py` for an example script.")
-		.def("getRoundness",&pyBodyContainer::getRoundness,"Returns roundness coefficient RC = R2/R1. R1 is the theoretical radius of a sphere, with same volume as clump. R2 is the minimum radius of a sphere, that imbeds clump. If just spheres are present RC = 1. If clumps are present 0 < RC < 1. Bodies can be excluded from the calculation by giving a list of ids: *O.bodies.getRoundness([ids])*.\n\nSee :ysrc:`examples/clumps/replaceByClumps-example.py` for an example script.")
+		.def("getRoundness",&pyBodyContainer::getRoundness,(python::arg("excludeList")=python::list()),"Returns roundness coefficient RC = R2/R1. R1 is the theoretical radius of a sphere, with same volume as clump. R2 is the minimum radius of a sphere, that imbeds clump. If just spheres are present RC = 1. If clumps are present 0 < RC < 1. Bodies can be excluded from the calculation by giving a list of ids: *O.bodies.getRoundness([ids])*.\n\nSee :ysrc:`examples/clumps/replaceByClumps-example.py` for an example script.")
 		.def("clear", &pyBodyContainer::clear,"Remove all bodies (interactions not checked)")
 		.def("erase", &pyBodyContainer::erase,"Erase body with the given id; all interaction will be deleted by InteractionLoop in the next step.")
 		.def("replace",&pyBodyContainer::replace);


Follow ups