yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2160
committer: Bruno Chareyre <bchareyre@r1arduina>
branch nick: trunk
timestamp: Mon 2010-04-19 15:30:50 +0200
message:
  - Some fixes in triangulation. Python array still crashing in TW. 
modified:
  pkg/dem/Engine/GlobalEngine/MicroMacroAnalyser.cpp
  pkg/dem/Engine/GlobalEngine/TesselationWrapper.cpp
  pkg/dem/Engine/GlobalEngine/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 'pkg/dem/Engine/GlobalEngine/MicroMacroAnalyser.cpp'
--- pkg/dem/Engine/GlobalEngine/MicroMacroAnalyser.cpp	2010-04-18 17:40:36 +0000
+++ pkg/dem/Engine/GlobalEngine/MicroMacroAnalyser.cpp	2010-04-19 13:30:50 +0000
@@ -77,7 +77,7 @@
 void MicroMacroAnalyser::setState(unsigned int state, bool save_states, bool computeIncrement)
 {
 	LOG_INFO("MicroMacroAnalyser::setState");
-	CGT::TriaxialState& TS = makeState(state, false);
+	CGT::TriaxialState& TS = makeState(state);
 	if (state == 2) {
 		analyser->Delta_epsilon(3,3) = analyser->TS1->eps3 - analyser->TS0->eps3;
 		analyser->Delta_epsilon(1,1) = analyser->TS1->eps1 - analyser->TS0->eps1;
@@ -144,7 +144,6 @@
 	}
 	TS.mean_radius /= Ng;//rayon moyen
 	LOG_INFO(" loaded : " << Ng << " grains with mean radius = " << TS.mean_radius);
-		
 	if (fictiousVtx.size()>=6){//boxes found, simulate them with big spheres
 		CGT::Point& Pmin = TS.box.base; CGT::Point& Pmax = TS.box.sommet; 
 		Real FAR = 1e4;
@@ -159,7 +158,6 @@
 		TS.grains[fictiousVtx[5]].sphere =
 		CGT::Sphere(CGT::Point(0.5*(Pmin.x()+Pmax.x()),0.5*(Pmax.y()+Pmin.y()),Pmax.z()+FAR*(Pmax.y()-Pmin.y())),FAR*(Pmax.y()-Pmin.y()));
 	}
-
 	InteractionContainer::iterator ii    = scene->interactions->begin();
 	InteractionContainer::iterator iiEnd = scene->interactions->end();
 	for (; ii!=iiEnd ; ++ii) {
@@ -194,22 +192,31 @@
 			c->frictional_work = 0;
 		}
 	}
-	//Save various parameters
-	TS.wszzh = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_top][1];//find_parameter("wszzh=", Statefile);
-	TS.wsxxd = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_right][0];//find_parameter("wsxxd=", Statefile);
-	TS.wsyyfa = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_front][2];//find_parameter("wsyyfa=", Statefile);
-	TS.eps3 = triaxialCompressionEngine->strain[2];//find_parameter("eps3=", Statefile);
-	TS.eps1 = triaxialCompressionEngine->strain[0];//find_parameter("eps1=", Statefile);
-	TS.eps2 = triaxialCompressionEngine->strain[1];//find_parameter("eps2=", Statefile);
-	TS.haut = triaxialCompressionEngine->height;//find_parameter("haut=", Statefile);
-	TS.larg = triaxialCompressionEngine->width;//find_parameter("larg=", Statefile);
-	TS.prof = triaxialCompressionEngine->depth;//find_parameter("prof=", Statefile);
-	TS.porom = analyser->ComputeMacroPorosity();//find_parameter("porom=", Statefile);
-	TS.ratio_f = triaxialCompressionEngine-> ComputeUnbalancedForce(scene);  //find_parameter("ratio_f=", Statefile);
-	if (filename!=NULL) {
-		//ostringstream oss;
-		TS.to_file(filename);
+	//Save various parameters if triaxialCompressionEngine is defined
+	if (!triaxialCompressionEngine) {
+		vector<shared_ptr<Engine> >::iterator itFirst = scene->engines.begin();
+		vector<shared_ptr<Engine> >::iterator itLast = scene->engines.end();
+		for (;itFirst!=itLast; ++itFirst) {
+			if ((*itFirst)->getClassName() == "TriaxialCompressionEngine") {
+				LOG_DEBUG("stress controller engine found");
+				triaxialCompressionEngine =  YADE_PTR_CAST<TriaxialCompressionEngine> (*itFirst);}
+		}
+		if (!triaxialCompressionEngine) LOG_ERROR("stress controller engine not found");}
+	
+	if (triaxialCompressionEngine) {
+		TS.wszzh = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_top][1];
+		TS.wsxxd = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_right][0];
+		TS.wsyyfa = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_front][2];
+		TS.eps3 = triaxialCompressionEngine->strain[2];//find_parameter("eps3=", Statefile);
+		TS.eps1 = triaxialCompressionEngine->strain[0];//find_parameter("eps1=", Statefile);
+		TS.eps2 = triaxialCompressionEngine->strain[1];//find_parameter("eps2=", Statefile);
+		TS.haut = triaxialCompressionEngine->height;//find_parameter("haut=", Statefile);
+		TS.larg = triaxialCompressionEngine->width;//find_parameter("larg=", Statefile);
+		TS.prof = triaxialCompressionEngine->depth;//find_parameter("prof=", Statefile);
+		TS.porom = 0/*analyser->ComputeMacroPorosity() crasher?*/;//find_parameter("porom=", Statefile);
+		TS.ratio_f = triaxialCompressionEngine-> ComputeUnbalancedForce(scene);  //find_parameter("ratio_f=", Statefile);
 	}
+	if (filename!=NULL) TS.to_file(filename);
 	return TS;
 }
 

=== modified file 'pkg/dem/Engine/GlobalEngine/TesselationWrapper.cpp'
--- pkg/dem/Engine/GlobalEngine/TesselationWrapper.cpp	2010-04-19 10:18:42 +0000
+++ pkg/dem/Engine/GlobalEngine/TesselationWrapper.cpp	2010-04-19 13:30:50 +0000
@@ -112,11 +112,8 @@
 }
 
 
-// static CGT::Point Pmin;
-// static CGT::Point Pmax;
-
 // is this a joke?? #include<limits> std::numeric_limits<double>::infinity();
-__attribute__((unused)) static double inf = 1e10; 
+//__attribute__((unused)) static double inf = 1e10; 
 double pminx=0;
 double pminy=0;
 double pminz=0;
@@ -316,14 +313,18 @@
 python::dict TesselationWrapper::getVolPoroDef(bool deformation)
 {
 		Scene* scene=Omega::instance().getScene().get();
-		CGT::Tesselation* pTes;
+		delete Tes;
+		CGT::TriaxialState* ts;
 		if (deformation){//use the final state to compute volumes
 			mma.analyser->ComputeParticlesDeformation();
-			pTes = &mma.analyser->TS1->tesselation();}
-		else pTes = &mma.analyser->TS0->tesselation();//no reason to use the final state if we don't want to compute deformations, keep using the initial
-		CGT::Tesselation& Tes = *pTes;
-		CGT::RTriangulation& Tri = Tes.Triangulation();
-		if (!scene->isPeriodic) AddBoundingPlanes();
+			Tes = &mma.analyser->TS1->tesselation();
+			ts = mma.analyser->TS1;
+			}
+		else {	Tes = &mma.analyser->TS0->tesselation();//no reason to use the final state if we don't want to compute deformations, keep using the initial
+			ts = mma.analyser->TS0;}
+		CGT::RTriangulation& Tri = Tes->Triangulation();
+		Pmin=ts->box.base; Pmax=ts->box.sommet;
+		//if (!scene->isPeriodic) AddBoundingPlanes();
 		ComputeVolumes();
 		int bodiesDim = scene->bodies->size();
 		int dim1[]={bodiesDim};

=== modified file 'pkg/dem/Engine/GlobalEngine/TesselationWrapper.hpp'
--- pkg/dem/Engine/GlobalEngine/TesselationWrapper.hpp	2010-04-18 17:40:36 +0000
+++ pkg/dem/Engine/GlobalEngine/TesselationWrapper.hpp	2010-04-19 13:30:50 +0000
@@ -78,22 +78,7 @@
 	/// 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);	
 	/// 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);//{
-// 		Scene* scene=Omega::instance().getScene().get();
-// 		int dim1[]={scene->bodies->size()};
-// 		int dim2[]={scene->bodies->size(),3};
-// 		numpy_boost<Real,1> mass(dim1);
-// 		numpy_boost<Real,2> vel(dim2);
-// 		FOREACH(const shared_ptr<Body>& b, *scene->bodies){
-// 			if(!b) continue;
-// 			mass[b->getId()]=b->state->mass;
-// 			VECTOR3R_TO_NUMPY(vel[b->getId()],b->state->vel);
-// 		}
-// 		python::dict ret;
-// 		ret["mass"]=mass;
-// 		ret["vel"]=vel;
-// 		return ret;
-// 	}
+	boost::python::dict getVolPoroDef(bool deformation);//FIXME ; unexplained crash for now
 
 	
 public:
@@ -111,7 +96,10 @@
 	facet_begin = Tes->Triangulation().finite_edges_begin();
 	facet_end = Tes->Triangulation().finite_edges_end();
 	facet_it = Tes->Triangulation().finite_edges_begin();
+	inf=1e10;
 	,/*py*/
+ 	.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("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).")
 	);	
 	DECLARE_LOGGER;
 };