yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2990
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Wed 2012-01-11 22:30:18 +0100
message:
  - improvements in TW interface so that one can save and load "state" files to be used for defining deformations (this will actually be used by E. Ando for analyzing experimental data!)
  - code cleaning 
modified:
  lib/triangulation/KinematicLocalisationAnalyser.cpp
  lib/triangulation/KinematicLocalisationAnalyser.hpp
  lib/triangulation/basicVTKwritter.cpp
  lib/triangulation/basicVTKwritter.hpp
  pkg/dem/TesselationWrapper.cpp
  pkg/dem/TesselationWrapper.hpp


--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription

=== modified file 'lib/triangulation/KinematicLocalisationAnalyser.cpp'
--- lib/triangulation/KinematicLocalisationAnalyser.cpp	2010-08-15 04:19:57 +0000
+++ lib/triangulation/KinematicLocalisationAnalyser.cpp	2012-01-11 21:30:18 +0000
@@ -82,6 +82,7 @@
 	bz2 = usebz2;
 	TS1->from_file(state_file1,/*use bz2?*/ bz2);
 	TS0->from_file(state_file0,/*use bz2?*/ bz2);
+	//FIXME: redundant?
 	Delta_epsilon(3,3) = TS1->eps3 - TS0->eps3;
 	Delta_epsilon(1,1) = TS1->eps1 - TS0->eps1;
 	Delta_epsilon(2,2) = TS1->eps2 - TS0->eps2;
@@ -193,15 +194,22 @@
 	return filteredList;
 }
 
+bool KinematicLocalisationAnalyser::DefToFile(const char* state_file1, const char* state_file0, const char* output_file_name, bool usebz2)
+{
+	consecutive = false;
+	bz2 = usebz2;
+	TS1->from_file(state_file1,/*use bz2?*/ bz2);
+	TS0->from_file(state_file0,/*use bz2?*/ bz2);
+	DefToFile(output_file_name);
+}
+
 bool KinematicLocalisationAnalyser::DefToFile(const char* output_file_name)
 {
 	ComputeParticlesDeformation();
 	Tesselation& Tes = TS1->tesselation();
 	RTriangulation& Tri = Tes.Triangulation();
-		
-	//basicVTKwritter vtk((unsigned int)(Tri.number_of_vertices()), (unsigned int)(Tri.number_of_finite_cells()));
 	basicVTKwritter vtk(n_real_vertices, n_real_cells);
-	vtk.open(output_file_name, "test bidon"); // <- what a pretty comment!
+	vtk.open(output_file_name, "Output file generated by Yade's KinematicLocalisationAnalyser");
 
 	vtk.begin_vertices();
 	RTriangulation::Finite_vertices_iterator  V_it = Tri.finite_vertices_begin();
@@ -822,18 +830,8 @@
 {
 	return output_file;
 }
-// Tenseur3 KinematicLocalisationAnalyser::Grad_u (Point &p1, Point &p2, Point &p3)
-// {
-//  Tenseur3 T;
-//  Vecteur V = (Deplacement(p1)+Deplacement(p2)+Deplacement(p3))/3.00;
-//  Grad_u(p1, p2, p3, V, T);
-//  return T;
-//  //Vecteur V = (Deplacement(p1)+Deplacement(p2)+Deplacement(p3))/3;
-// }
-
-//Vecteur KinematicLocalisationAnalyser::Deplacement (Point &p) {return (p-CGAL::ORIGIN)/100;}
-
-Vecteur KinematicLocalisationAnalyser::Deplacement(Finite_cells_iterator cell, int facet)    // D�p. moyen sur une facette
+
+Vecteur KinematicLocalisationAnalyser::Deplacement(Finite_cells_iterator cell, int facet)  
 {
 	Vecteur v(0.f, 0.f, 0.f);
 	int id;// ident. de la particule
@@ -846,9 +844,9 @@
 			Vecteur meanFieldDisp =Vecteur(TS0->grain(id).sphere.point().x(), TS0->grain(id).sphere.point().y(), TS0->grain(id).sphere.point().z())-fixedPoint;
 			if (1){//fluctuations
 				meanFieldDisp = Vecteur(
-				meanFieldDisp[0]*(TS1->larg-TS0->larg)/TS0->larg,
-				meanFieldDisp[1]*(TS1->haut-TS0->haut)/TS0->haut,
-				meanFieldDisp[2]*(TS1->prof-TS0->prof)/TS0->prof);
+				meanFieldDisp[0]*Delta_epsilon(0,0),
+				meanFieldDisp[1]*Delta_epsilon(1,1),
+				meanFieldDisp[2]*Delta_epsilon(2,2));
 			} else meanFieldDisp=Vecteur(0,0,0);
 			if (consecutive)
 				v = v + TS1->grain(id).translation-meanFieldDisp;
@@ -876,11 +874,7 @@
 							  (cell->vertex(l_vertices[facet][1])->point())) /2.f;
 	Somme(T, V, S);
 }
-// void KinematicLocalisationAnalyser::Grad_u (Point &p1, Point &p2, Point &p3, Vecteur &V, Tenseur3& T) // rotation 1->2->3 orient�e vers l'ext�rieur
-// {
-//  Vecteur S = 0.5*cross_product(p2-p1, p3-p2);
-//      Somme (T, V, S);
-// }
+
 
 void KinematicLocalisationAnalyser::Grad_u(Finite_cells_iterator cell,
 		Tenseur3& T, bool vol_divide)// Calcule le gradient de d�p.
@@ -894,91 +888,6 @@
 	if (vol_divide) T/= Tesselation::Volume(cell);
 }
 
-
-
-
-
-
-// Calcul du tenseur d'orientation des voisins
-// Tenseur_sym3 Orientation_voisins (Tesselation &Tes)
-// {
-//  RTriangulation& T = Tes.Triangulation();
-//  Tenseur3 Tens;
-//  Vecteur v;
-//  long Nv = 0; //nombre de voisins
-//  for (Edge_iterator ed_it = T.edges_begin(); ed_it != T.edges_end(); ed_it++)
-//  {
-//   if (!T.is_infinite(*ed_it))
-//   {
-//    Nv++;
-//    v = T.segment(*ed_it).to_vector()/sqrt(T.segment(*ed_it).squared_length());
-//    for (int i=1; i<4; i++) for (int j=3; j>=i; j--)Tens(i,j) += v[i-1]*v[j-1];
-//   }
-//  }
-//  Tens /= Nv;
-//  return Tens;
-// }
-
-
-
-
-
-
-//
-//;------------------------------------------------
-//;orient_n_Y
-//;calcul de la distribution des orientations de contact
-//;(angle entre n et l'axe y)
-//;param�tres d'entr�e : N_classes (nombre de classes), e_filtre (�paisseur de la bordure supprim�e)
-//;sortie : table 20
-//
-//def orient_n_Z
-//
-//X_min = w_x(waddg) + e_filtre
-//X_max = l + w_x(waddd) - e_filtre
-//Y_min = w_y(waddfa) + e_filtre
-//Y_max = l + w_y(waddfo) - e_filtre
-//Z_min = w_z(waddb) + e_filtre
-//Z_max = h + w_z(waddh) - e_filtre
-//DZ_classes = 1./(N_classes*1.)
-//
-//command
-//table 20 erase
-//table 21 erase
-//endcommand
-//loop temp (0, N_classes)
-// xtable (20, temp+1) = temp * DZ_classes
-//endloop
-//
-//cp=contact_head
-//N_filtre1 = 0
-//n_cont = 0
-//loop while cp#null
-// if c_nforce(cp) # 0
-// n_cont = n_cont +1
-// if filtre1 = 1
-// N_filtre1 = N_filtre1 +1
-//  classe = int(abs(c_zun(cp))/DZ_classes) + 1
-//  ytable (20, classe) = ytable (20, classe) +1
-// endif
-// endif
-//cp=c_next(cp)
-//endloop
-//
-//command
-//print N_filtre1 n_cont
-//set logfile orient.txt
-//set log on ov
-//pr table 20
-//set log off
-//set logfile res.dat
-//endcommand
-//;Normalisation :
-//
-//end
-
-
-
 const vector<Tenseur3>& KinematicLocalisationAnalyser::ComputeParticlesDeformation(void)
 {
 	Tesselation& Tes = TS1->tesselation();
@@ -1047,7 +956,7 @@
 		if (V_it->info().v()) ParticleDeformation[V_it->info().id()]/=V_it->info().v();
 	}
 	grad_u_total_g /= v_total_g;
-	if(0){
+	if(1){
 	cerr << "sym_grad_u_total_g (wrong averaged strain):"<< endl << Tenseur_sym3(grad_u_total_g) << endl;
 	if (v_total) grad_u_total /= v_total;
 	cerr << "Total volume = " << v_total << ", grad_u = " << endl << grad_u_total << endl << "sym_grad_u (true average strain): " << endl << Tenseur_sym3(grad_u_total) << endl;

=== modified file 'lib/triangulation/KinematicLocalisationAnalyser.hpp'
--- lib/triangulation/KinematicLocalisationAnalyser.hpp	2010-08-15 04:19:57 +0000
+++ lib/triangulation/KinematicLocalisationAnalyser.hpp	2012-01-11 21:30:18 +0000
@@ -48,6 +48,7 @@
 		bool DistribsToFile (const char* output_file_name);
 		///Write the averaged deformation on each grain in a file (vertices and cells lists included in the file), no need to call ComputeParticlesDeformation()
 		bool DefToFile (const char* output_file_name = "deformations");
+		bool DefToFile (const char* state_file1, const char* state_file0, const char* output_file_name="deformation.vtk", bool usebz2=false);
 		///Save/Load states using bz2 compression
 		bool bz2;
 		ofstream& ContactDistributionToFile ( ofstream& output_file );

=== modified file 'lib/triangulation/basicVTKwritter.cpp'
--- lib/triangulation/basicVTKwritter.cpp	2010-08-15 04:19:57 +0000
+++ lib/triangulation/basicVTKwritter.cpp	2012-01-11 21:30:18 +0000
@@ -20,6 +20,11 @@
   return true;
 }
 
+bool basicVTKwritter::close()
+{
+  file.close(); return true;
+}
+
 void basicVTKwritter::begin_vertices()
 {
   file << "POINTS " << nbVertices << " float" << endl;
@@ -60,8 +65,8 @@
 {
 switch(pos)
 {
-  case POINT_DATA : file << "POINT_DATA " << nbVertices << endl; break;
-  case CELL_DATA  : file << "CELL_DATA " << nbCells << endl; break;
+  case POINT_DATA : if (!hasPointData) {file << "POINT_DATA " << nbVertices << endl; hasPointData=true;} break;
+  case CELL_DATA  : if (!hasCellData) {file << "CELL_DATA " << nbCells << endl; hasCellData=true;} break;
 }
 
   switch (name)
@@ -77,11 +82,12 @@
   case FLOAT : file << " float"; break;
 }
 
-  if (name == SCALARS) file << " 1" << endl;
-  else                file << endl;
+  if (name == SCALARS) {
+	  file << " 1" << endl;
+	  file << "LOOKUP_TABLE default" << endl;}
+  else file << endl;
+}  
   
-  file << "LOOKUP_TABLE default" << endl;
-}
 
 
 void basicVTKwritter::write_data(float value)

=== modified file 'lib/triangulation/basicVTKwritter.hpp'
--- lib/triangulation/basicVTKwritter.hpp	2010-08-15 04:19:57 +0000
+++ lib/triangulation/basicVTKwritter.hpp	2012-01-11 21:30:18 +0000
@@ -11,10 +11,14 @@
 {
 	std::ofstream file;
 	unsigned int nbVertices, nbCells;
-
-	basicVTKwritter(unsigned int nV, unsigned int nC) : nbVertices(nV),nbCells(nC) { }
+	bool hasPointData;
+	bool hasCellData;
+	
+	
+	basicVTKwritter(unsigned int nV, unsigned int nC) : nbVertices(nV),nbCells(nC),hasPointData(false),hasCellData(false) {}
 	
 	bool open(const char * filename, const char * comment);
+	bool close();
 	
 	void begin_vertices();
 	void write_point(float x, float y, float z);

=== modified file 'pkg/dem/TesselationWrapper.cpp'
--- pkg/dem/TesselationWrapper.cpp	2011-10-21 17:43:08 +0000
+++ pkg/dem/TesselationWrapper.cpp	2012-01-11 21:30:18 +0000
@@ -117,22 +117,14 @@
 	//cerr << " loaded : " << Ng<<", triangulated : "<<TW.n_spheres<<", mean radius = " << TW.mean_radius<<endl;
 }
 
-
-// is this a joke?? #include<limits> std::numeric_limits<double>::infinity();
-//__attribute__((unused)) static double inf = 1e10;
-double pminx=0;
+/*double pminx=0;
 double pminy=0;
 double pminz=0;
 double pmaxx=0;
 double pmaxy=0;
-double pmaxz=0;
+double pmaxz=0;*/
 double thickness = 0;
 
-//  Finite_edges_iterator facet_it; //an edge in a triangulation is a facet in corresponding tesselation, remember...
-//That explain the name.
-
-
-
 TesselationWrapper::~TesselationWrapper() { if (Tes) delete Tes;}
 
 void TesselationWrapper::clear(void)
@@ -170,7 +162,7 @@
 // 	clock.top("Triangulation");
 }
 
-double TesselationWrapper::Volume(unsigned int id) {return (Tes->Max_id() > id) ? Tes->Volume(id) : 0;}
+double TesselationWrapper::Volume(unsigned int id) {return (Tes->Max_id() > id) ? Tes->Volume(id) : -1;}
 
 bool TesselationWrapper::insert(double x, double y, double z, double rad, unsigned int id)
 {
@@ -319,13 +311,31 @@
 
 void TesselationWrapper::setState (bool state){ mma.setState(state ? 2 : 1);}
 
+void TesselationWrapper::loadState (string filename, bool stateNumber, bool bz2){
+	CGT::TriaxialState& TS = stateNumber? *(mma.analyser->TS1) :*( mma.analyser->TS0);
+	TS.from_file(filename.c_str(),bz2);
+}
+
+void TesselationWrapper::saveState (string filename, bool stateNumber, bool bz2){
+	CGT::TriaxialState& TS = stateNumber? *(mma.analyser->TS1) :*( mma.analyser->TS0);
+	TS.to_file(filename.c_str(),bz2);
+}
+
+void TesselationWrapper::defToVtkWithInput (string inputFile1, string inputFile2, string outputFile, bool bz2){
+	mma.analyser->DefToFile(inputFile1.c_str(),inputFile2.c_str(),outputFile.c_str(),bz2);
+}
+
+void TesselationWrapper::defToVtk (string outputFile){
+	mma.analyser->DefToFile(outputFile.c_str());
+}
+
 python::dict TesselationWrapper::getVolPoroDef(bool deformation)
 {
 		Scene* scene=Omega::instance().getScene().get();
 		delete Tes;
 		CGT::TriaxialState* ts;
 		if (deformation){//use the final state to compute volumes
-			mma.analyser->ComputeParticlesDeformation();
+			const vector<CGT::Tenseur3>& def = mma.analyser->ComputeParticlesDeformation();
 			Tes = &mma.analyser->TS1->tesselation();
 			ts = mma.analyser->TS1;
 			}
@@ -335,7 +345,8 @@
 		Pmin=ts->box.base; Pmax=ts->box.sommet;
 		//if (!scene->isPeriodic) AddBoundingPlanes();
 		ComputeVolumes();
-		int bodiesDim = scene->bodies->size();
+		int bodiesDim = Tes->Max_id() + 1; //=scene->bodies->size();
+		cerr<<"bodiesDim="<<bodiesDim<<endl;
 		int dim1[]={bodiesDim};
 		int dim2[]={bodiesDim,9};
 		/// This is the code that needs numpy include

=== modified file 'pkg/dem/TesselationWrapper.hpp'
--- pkg/dem/TesselationWrapper.hpp	2010-11-17 15:19:09 +0000
+++ pkg/dem/TesselationWrapper.hpp	2012-01-11 21:30:18 +0000
@@ -82,6 +82,12 @@
 
 	/// make the current state the initial (0) or final (1) configuration for the definition of displacement increments, use only state=0 if you just want to get only volmumes and porosity
 	void setState (bool state=0);
+	void loadState (string fileName, bool stateNumber=0, bool bz2=false);
+	void saveState (string fileName, bool stateNumber=0, bool bz2=false);
+	/// read two state files and write per-particle deformation to a vtk file. The second variant uses existing states.
+ 	void defToVtkWithInput (string inputFile1, string inputFile2, string outputFile="def.vtk", bool bz2=false);
+	void defToVtk (string outputFile="def.vtk");
+
 	/// return python array containing voronoi volumes, per-particle porosity, and optionaly per-particle deformation, if states 0 and 1 have been assigned
 	boost::python::dict getVolPoroDef(bool deformation);//FIXME ; unexplained crash for now
 
@@ -93,7 +99,7 @@
 	CGT::Finite_edges_iterator facet_it;
 	MicroMacroAnalyser mma;
 
-	YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(TesselationWrapper,GlobalEngine,"Handle the triangulation of spheres in a scene, build tesselation on request, and give access to computed quantities : currently volume and porosity of each Voronoï sphere. Example script :\n tt=TriaxialTest()\ntt.generate('test.xml')\nO.load('test.xml')\nO.run(100) //for unknown reasons, this procedure crashes at iteration 0\nTW=TesselationWrapper()\nTW.triangulate() //compute regular Delaunay triangulation, don't construct tesselation\nTW.computeVolumes() //will silently tesselate the packing\nTW.volume(10) //get volume associated to sphere of id 10",
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(TesselationWrapper,GlobalEngine,"Handle the triangulation of spheres in a scene, build tesselation on request, and give access to computed quantities : currently volume and porosity of each Voronoï sphere. Example script :\n\n tt=TriaxialTest()\n\ntt.generate('test.xml')\n\nO.load('test.xml')\n\nO.run(100)\n\nTW=TesselationWrapper()\n\nTW.triangulate() //compute regular Delaunay triangulation, don't construct tesselation\n\nTW.computeVolumes() //will silently tesselate the packing\n\nTW.volume(10) //get volume associated to sphere of id 10",
 	((unsigned int,n_spheres,0,,"|ycomp|"))
 	,/*ctor*/
   	Tes = new CGT::Tesselation;
@@ -105,7 +111,11 @@
 	,/*py*/
 	.def("triangulate",&TesselationWrapper::insertSceneSpheres,(python::arg("reset")=true),"triangulate spheres of the packing")
  	.def("setState",&TesselationWrapper::setState,(python::arg("state")=0),"Make the current state the initial (0) or final (1) configuration for the definition of displacement increments, use only state=0 if you just want to get only volmumes and porosity.")
+ 	.def("loadState",&TesselationWrapper::loadState,(python::arg("stateNumber")=0),"Load a file with positions to define state 0 or 1.")
+ 	.def("saveState",&TesselationWrapper::saveState,(python::arg("outputFile")="state",python::arg("bz2")=true),"Save a file with positions, can be later reloaded in order to define state 0 or 1.")
  	.def("volume",&TesselationWrapper::Volume,(python::arg("id")=0),"Returns the volume of Voronoi's cell of a sphere.")
+ 	.def("defToVtk",&TesselationWrapper::defToVtk,(python::arg("outputFile")="def.vtk"),"Write local deformations in vtk format from states 0 and 1.")
+ 	.def("defToVtkWithInput",&TesselationWrapper::defToVtkWithInput,(python::arg("input1")="state1",python::arg("input2")="state2",python::arg("outputFile")="def.vtk",python::arg("bz2")=false),"Write local deformations in vtk format from position input files.")
  	.def("computeVolumes",&TesselationWrapper::ComputeVolumes,"Compute volumes of all Voronoi's cells.")
 	.def("getVolPoroDef",&TesselationWrapper::getVolPoroDef,(python::arg("deformation")=false),"Return a table with per-sphere computed quantities. Include deformations on the increment defined by states 0 and 1 if deformation=True (make sure to define states 0 and 1 consistently).")
 	);