yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2008-08-29 22:54:03 +0200 (Fri, 29 Aug 2008)
New Revision: 1497

Modified:
   trunk/gui/py/_utils.cpp
   trunk/gui/py/eudoxos.py
Log:
1. Add code for computing histogram of interaction directions in given plane in _utils.
2. This is code is used by yade.eudoxos.plotDirections() to create the actual (nice) figure.



Modified: trunk/gui/py/_utils.cpp
===================================================================
--- trunk/gui/py/_utils.cpp	2008-08-28 21:01:20 UTC (rev 1496)
+++ trunk/gui/py/_utils.cpp	2008-08-29 20:54:03 UTC (rev 1497)
@@ -3,6 +3,8 @@
 #include<yade/core/MetaBody.hpp>
 #include<yade/core/Omega.hpp>
 #include<yade/pkg-common/Sphere.hpp>
+#include<yade/pkg-dem/SpheresContactGeometry.hpp>
+#include<cmath>
 using namespace boost::python;
 
 python::tuple vec2tuple(const Vector3r& v){return boost::python::make_tuple(v[0],v[1],v[2]);};
@@ -49,15 +51,41 @@
 	}
 }
 
-
 Real PWaveTimeStep(){return Shop::PWaveTimeStep();};
 
+
+/* return histogram ([bin1Min,bin2Min,…],[value1,value2,…]) from projections of interactions
+ * to the plane perpendicular to axis∈[0,1,2]; the number of bins can be specified and they cover
+ * the range (0,π), since interactions are bidirecional, hence periodically the same on (π,2π).
+ *
+ * only contacts using SpheresContactGeometry are considered.
+ */
+python::tuple interactionAnglesHistogram(int axis, size_t bins=20){
+	if(axis<0||axis>2) throw invalid_argument("Axis must be from {0,1,2}=x,y,z.");
+	Real binStep=Mathr::PI/bins; int axis2=(axis+1)%3, axis3=(axis+2)%3;
+	vector<Real> cummProj(bins,0.);
+	FOREACH(const shared_ptr<Interaction>& i, *Omega::instance().getRootBody()->transientInteractions){
+		shared_ptr<SpheresContactGeometry> scg=dynamic_pointer_cast<SpheresContactGeometry>(i->interactionGeometry);
+		if(!scg) continue;
+		Vector3r n(scg->normal); n[axis]=0.; Real nLen=n.Length();
+		Real theta=acos(n[axis2]/nLen)*(n[axis3]>0?1:-1); if(theta<0) theta+=Mathr::PI;
+		int binNo=theta/binStep;
+		cummProj[binNo]+=nLen;
+	}
+	python::list val,binMid;
+	for(size_t i=0; i<(size_t)bins; i++){ val.append(cummProj[i]); binMid.append(i*binStep);}
+	return python::make_tuple(binMid,val);
+}
+BOOST_PYTHON_FUNCTION_OVERLOADS(interactionAnglesHistogram_overloads,interactionAnglesHistogram,1,2);
+
+
 BOOST_PYTHON_MODULE(_utils){
 	def("PWaveTimeStep",PWaveTimeStep);
 	def("aabbExtrema",aabbExtrema);
 	def("negPosExtremeIds",negPosExtremeIds,negPosExtremeIds_overloads(args("axis","distFactor")));
 	def("coordsAndDisplacements",coordsAndDisplacements);
 	def("setRefSe3",setRefSe3);
+	def("interactionAnglesHistogram",interactionAnglesHistogram,interactionAnglesHistogram_overloads(args("axis","bins")));
 }
 
 

Modified: trunk/gui/py/eudoxos.py
===================================================================
--- trunk/gui/py/eudoxos.py	2008-08-28 21:01:20 UTC (rev 1496)
+++ trunk/gui/py/eudoxos.py	2008-08-29 20:54:03 UTC (rev 1497)
@@ -4,6 +4,18 @@
 # I doubt there functions will be useful for anyone besides me.
 #
 
+
+def plotDirections(bins=20):
+	"Plot 3 histograms for distribution of interaction directions, in yz,xz and xy planes."
+	import pylab
+	from yade import utils
+	for axis in [0,1,2]:
+		d=utils.interactionAnglesHistogram(axis,bins)
+		fc=[0,0,0]; fc[axis]=1.
+		pylab.subplot(220+axis+1,polar=True);
+		bar(d[0],d[1],width=math.pi/(1.2*bins),fc=fc,alpha=.7,label=['yz','xz','xy'][axis])
+	pylab.show()
+
 def estimatePoissonYoung(principalAxis,stress=0,plot=False):
 	"""Estimate Poisson's ration given the "principal" axis of straining.
 	For every base direction, homogenized strain is computed