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
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
