yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3664
committer: Raphael Maurin <raph_maurin@xxxxxxxxxxx>
timestamp: Tue 2015-06-02 10:51:06 +0200
message:
  getStressProfile: output format modification + add calculation of granular temperature and kinetic stress tensor
modified:
  pkg/dem/Shop.hpp
  pkg/dem/Shop_02.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/Shop.hpp'
--- pkg/dem/Shop.hpp	2015-03-03 19:06:49 +0000
+++ pkg/dem/Shop.hpp	2015-06-02 08:51:06 +0000
@@ -122,7 +122,7 @@
 		static Matrix3r getCapillaryStress(Real volume=0, bool mindlin=false);
 		static Matrix3r stressTensorOfPeriodicCell() { LOG_WARN("Shop::stressTensorOfPeriodicCelli is DEPRECATED: use getStress instead"); return Shop::getStress(); }
 
-		static vector<Matrix3r> getStressProfile(Real volume, int nCell, Real dz, Real zRef, vector<Real> vPartAverageX, vector<Real> vPartAverageY, vector<Real> vPartAverageZ);
+		static py::tuple getStressProfile(Real volume, int nCell, Real dz, Real zRef, vector<Real> vPartAverageX, vector<Real> vPartAverageY, vector<Real> vPartAverageZ);
 
 		//! Compute overall ("macroscopic") stress of periodic cell, returning 2 tensors
 		//! (contribution of normal and shear forces)

=== modified file 'pkg/dem/Shop_02.cpp'
--- pkg/dem/Shop_02.cpp	2015-04-24 06:07:03 +0000
+++ pkg/dem/Shop_02.cpp	2015-06-02 08:51:06 +0000
@@ -373,14 +373,17 @@
 }
 
 
-
-vector<Matrix3r> Shop::getStressProfile(Real volume, int nCell, Real dz, Real zRef, vector<Real> vPartAverageX, vector<Real> vPartAverageY, vector<Real> vPartAverageZ){
+py::tuple Shop::getStressProfile(Real volume, int nCell, Real dz, Real zRef, vector<Real> vPartAverageX, vector<Real> vPartAverageY, vector<Real> vPartAverageZ){
 	int minZ=0;
 	int maxZ=0;
 	Real minPosZ=0.;
 	Real maxPosZ=0.;
 	Scene* scene=Omega::instance().getScene().get();
 	vector<Matrix3r> stressTensorProfile(nCell,Matrix3r::Zero());
+	vector<Matrix3r> kineticStressTensorProfile(nCell,Matrix3r::Zero());
+	vector<Real> granularTemperatureProfile(nCell,0.0);
+	vector<Real> numPart(nCell,0.0);
+
 	const bool isPeriodic = scene->isPeriodic;
 
 	//
@@ -393,8 +396,14 @@
 			Vector3r vFluct = b->state->vel - Vector3r(vPartAverageX[Np],vPartAverageY[Np],vPartAverageZ[Np]); 
 			//Classical dynamical expression of the stress tensor
 			stressTensorProfile[Np]+= -1/volume*b->state->mass*vFluct*vFluct.transpose(); 
+			kineticStressTensorProfile[Np]+= -1/volume*b->state->mass*vFluct*vFluct.transpose(); 
+			granularTemperatureProfile[Np] += 1/3.*(pow(vFluct[0],2) + pow(vFluct[1],2) + pow(vFluct[2],2));
+			numPart[Np]+=1.;
 		}
 	}
+	for(int n=0;n<nCell;n++) {
+		if (numPart[n]>0) granularTemperatureProfile[n]/=numPart[n];
+	}
 	
 	//
 	//Love Weber contribution (same as getStress(), but with different layers)
@@ -467,7 +476,7 @@
 			}
 		}
 	}
-	return stressTensorProfile;
+	return py::make_tuple(stressTensorProfile,kineticStressTensorProfile,granularTemperatureProfile);
 }