yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01294
[svn] r1794 - in trunk: gui gui/py scripts/test
Author: eudoxos
Date: 2009-06-13 17:09:51 +0200 (Sat, 13 Jun 2009)
New Revision: 1794
Added:
trunk/gui/py/_packPredicates.cpp
trunk/gui/py/pack.py
Modified:
trunk/gui/SConscript
trunk/gui/py/PythonUI.cpp
trunk/gui/py/PythonUI_rc.py
trunk/scripts/test/regular-sphere-pack.py
Log:
1. Move all python modules to PREFIX/lib/yadeSUFFIX/py/yade instead of */gui/yade
2. Add pack module for generating regular packings (thanks to Anton Gladky for his ideas & code)
3. Add _packPredicates module for solid inclusion testing.
4. Update regular-sphere-pack.py to test new packing features.
Modified: trunk/gui/SConscript
===================================================================
--- trunk/gui/SConscript 2009-06-09 11:23:06 UTC (rev 1793)
+++ trunk/gui/SConscript 2009-06-13 15:09:51 UTC (rev 1794)
@@ -49,7 +49,7 @@
env.File('PythonUI_rc.py','py'),
])
if 'qt3' not in env['exclude']:
- env.Install('$PREFIX/lib/yade$SUFFIX/gui/yade',[
+ env.Install('$PREFIX/lib/yade$SUFFIX/py/yade',[
env.SharedLibrary('_qt',['qt3/QtGUI-python.cpp'],SHLIBPREFIX='',LIBS=env['LIBS']+['QtGUI'],CPPPATH=env['CPPPATH']+[env['buildDir']+'/gui/qt3']), # CPPPATH is for files generated by moc which are indirectly included
env.File('qt.py','qt3'),
])
@@ -57,7 +57,7 @@
env.InstallAs('$PREFIX/bin/yade$SUFFIX-multi',env.File('yade-multi','py'))
# python modules are one level deeper so that you can say: from yade.wrapper import *
- env.Install('$PREFIX/lib/yade$SUFFIX/gui/yade',[
+ env.Install('$PREFIX/lib/yade$SUFFIX/py/yade',[
env.SharedLibrary('wrapper',['py/yadeControl.cpp'],SHLIBPREFIX='',
LIBS=env['LIBS']+['XMLFormatManager','yade-factory','yade-serialization','Shop',
'BoundingVolumeMetaEngine',
@@ -75,6 +75,7 @@
),
env.SharedLibrary('_utils',['py/_utils.cpp'],SHLIBPREFIX='',
LIBS=env['LIBS']+['Shop','ConcretePM']),
+ env.SharedLibrary('_packPredicates',['py/_packPredicates.cpp'],SHLIBPREFIX=''),
env.SharedLibrary('_eudoxos',['py/_eudoxos.cpp'],SHLIBPREFIX='',CXXFLAGS=env['CXXFLAGS']+([] if not os.path.exists('../../brefcom-mm.hh') else ['-include','../brefcom-mm.hh']),LIBS=env['LIBS']+['Shop','ConcretePM']),
env.SharedLibrary('log',['py/log.cpp'],SHLIBPREFIX=''),
env.File('__init__.py','py'),
@@ -86,6 +87,7 @@
env.File('linterpolation.py','py'),
env.File('timing.py','py'),
env.File('PythonTCPServer.py','py'),
+ env.File('pack.py','py'),
])
if os.path.exists('../../brefcom-mm.hh'): Depends('py/_eudoxos.cpp','../../brefcom-mm.hh')
Modified: trunk/gui/py/PythonUI.cpp
===================================================================
--- trunk/gui/py/PythonUI.cpp 2009-06-09 11:23:06 UTC (rev 1793)
+++ trunk/gui/py/PythonUI.cpp 2009-06-13 15:09:51 UTC (rev 1794)
@@ -79,7 +79,7 @@
PyGILState_STATE pyState = PyGILState_Ensure();
LOG_DEBUG("Got Global Interpreter Lock, good.");
/* import yade (for startUI()) and yade.runtime (initially empty) namespaces */
- PyRun_SimpleString("import sys; sys.path.insert(0,'" PREFIX "/lib/yade" SUFFIX "/gui')");
+ PyRun_SimpleString("import sys; sys.path.insert(0,'" PREFIX "/lib/yade" SUFFIX "/py')");
PyRun_SimpleString("import yade");
PyRun_SimpleString("from __future__ import division");
Modified: trunk/gui/py/PythonUI_rc.py
===================================================================
--- trunk/gui/py/PythonUI_rc.py 2009-06-09 11:23:06 UTC (rev 1793)
+++ trunk/gui/py/PythonUI_rc.py 2009-06-13 15:09:51 UTC (rev 1794)
@@ -122,7 +122,7 @@
| | (_| | |_| | __/ | |__| (_) | | | \__ \ (_) | | __/
|_|\__,_|____/ \___| \____\___/|_| |_|___/\___/|_|\___|
""",exit_msg='Bye.'
- ,rc_override={'execfile':[runtime.prefix+'/lib/yade'+runtime.suffix+'/gui/yade/ipython.py']})
+ ,rc_override={'execfile':[runtime.prefix+'/lib/yade'+runtime.suffix+'/py/yade/ipython.py']})
ipshell()
# save history -- a workaround for atexit handlers not being run (why?)
Added: trunk/gui/py/_packPredicates.cpp
===================================================================
--- trunk/gui/py/_packPredicates.cpp 2009-06-09 11:23:06 UTC (rev 1793)
+++ trunk/gui/py/_packPredicates.cpp 2009-06-13 15:09:51 UTC (rev 1794)
@@ -0,0 +1,93 @@
+// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
+#include<boost/python.hpp>
+#include<yade/extra/boost_python_len.hpp>
+#include<yade/lib-base/Logging.hpp>
+#include<yade/lib-base/yadeWm3.hpp>
+#include<Wm3Vector3.h>
+
+using namespace boost;
+using namespace std;
+#ifdef LOG4CXX
+ log4cxx::LoggerPtr logger=log4cxx::Logger::getLogger("yade.predicates");
+#endif
+
+/*
+This file contains various predicates that say whether a given point is within the solid,
+or, not closer than "pad" to its boundary, if pad is nonzero
+Besides the (point,pad) operator, each predicate defines aabb() method that returns
+(min,max) tuple defining minimum and maximum point of axis-aligned bounding box
+for the predicate.
+
+These classes are primarily used for yade.pack.* functions creating packings.
+See scripts/test/regular-sphere-pack.py for an example.
+
+*/
+
+// aux functions
+python::tuple vec2tuple(const Vector3r& v){return boost::python::make_tuple(v[0],v[1],v[2]);}
+Vector3r tuple2vec(const python::tuple& t){return Vector3r(python::extract<double>(t[0])(),python::extract<double>(t[1])(),python::extract<double>(t[2])());}
+void ttuple2vvec(const python::tuple& t, Vector3r& v1, Vector3r& v2){ v1=tuple2vec(python::extract<python::tuple>(t[0])()); v2=tuple2vec(python::extract<python::tuple>(t[1])()); }
+python::tuple vvec2ttuple(const Vector3r&v1, const Vector3r&v2){ return python::make_tuple(vec2tuple(v1),vec2tuple(v2)); }
+
+
+/*! Sphere predicate */
+class inSphere {
+ Vector3r center; Real radius;
+public:
+ inSphere(python::tuple _center, Real _radius){center=tuple2vec(_center); radius=_radius;}
+ bool operator()(python::tuple _pt, Real pad=0.){
+ Vector3r pt=tuple2vec(_pt);
+ return ((pt-center).Length()-pad<=radius-pad);
+ }
+ python::tuple aabb(){return vvec2ttuple(Vector3r(center[0]-radius,center[1]-radius,center[2]-radius),Vector3r(center[0]+radius,center[1]+radius,center[2]+radius));}
+};
+
+/* Axis-aligned box predicate */
+class inAlignedBox{
+ Vector3r mn, mx;
+public:
+ inAlignedBox(python::tuple _mn, python::tuple _mx){mn=tuple2vec(_mn); mx=tuple2vec(_mx);}
+ bool operator()(python::tuple _pt, Real pad=0.){
+ Vector3r pt=tuple2vec(_pt);
+ return
+ mn[0]+pad<=pt[0] && mx[0]-pad>=pt[0] &&
+ mn[1]+pad<=pt[1] && mx[1]-pad>=pt[1] &&
+ mn[2]+pad<=pt[2] && mx[2]-pad>=pt[2];
+ }
+ python::tuple aabb(){ return vvec2ttuple(mn,mx); }
+};
+
+/* Arbitrarily oriented cylinder predicate */
+class inCylinder{
+ Vector3r c1,c2,c12; Real radius,ht;
+public:
+ inCylinder(python::tuple _c1, python::tuple _c2, Real _radius){c1=tuple2vec(_c1); c2=tuple2vec(_c2); c12=c2-c1; radius=_radius; ht=c12.Length(); }
+ bool operator()(python::tuple _pt, Real pad=0.){
+ Vector3r pt=tuple2vec(_pt);
+ Real u=(pt.Dot(c12)-c1.Dot(c12))/(ht*ht); // normalized coordinate along the c1--c2 axis
+ if((u*ht<0+pad) || (u*ht>ht-pad)) return false; // out of cylinder along the axis
+ Real axisDist=((pt-c1).Cross(pt-c2)).Length()/ht;
+ if(axisDist>radius-pad) return false;
+ return true;
+ }
+ python::tuple aabb(){
+ // see http://www.gamedev.net/community/forums/topic.asp?topic_id=338522&forum_id=20&gforum_id=0 for the algorithm
+ Vector3r& A(c1); Vector3r& B(c2);
+ Vector3r k(
+ sqrt((pow(A[1]-B[1],2)+pow(A[2]-B[2],2)))/ht,
+ sqrt((pow(A[0]-B[0],2)+pow(A[2]-B[2],2)))/ht,
+ sqrt((pow(A[0]-B[0],2)+pow(A[1]-B[1],2)))/ht);
+ Vector3r mn(min(A[0],B[0]),min(A[1],B[1]),min(A[2],B[2])), mx(max(A[0],B[0]),max(A[1],B[1]),max(A[2],B[2]));
+ return vvec2ttuple(mn-radius*k,mx+radius*k);
+ }
+};
+
+BOOST_PYTHON_MODULE(_packPredicates){
+ boost::python::class_<inSphere>("inSphere","Sphere predicate.",python::init<python::tuple,Real>("Ctor taking center (as a 3-tuple) and radius"))
+ .def("__call__",&inSphere::operator(),"Tell whether given point lies within this sphere, still having 'pad' space to the solid boundary").def("aabb",&inSphere::aabb,"Return minimum and maximum values for AABB");
+ boost::python::class_<inAlignedBox>("inAlignedBox","Axis-aligned box predicate",python::init<python::tuple,python::tuple>("Ctor taking minumum and maximum points of the box (as 3-tuples)."))
+ .def("__call__",&inAlignedBox::operator(),"Tell whether given point lies within this box, still having 'pad' space to the solid boundary").def("aabb",&inAlignedBox::aabb,"Return minimum and maximum values for AABB");
+ boost::python::class_<inCylinder>("inCylinder","Cylinder predicate",python::init<python::tuple,python::tuple,Real>("Ctor taking centers of the lateral walls (as 3-tuples) and radius."))
+ .def("__call__",&inCylinder::operator(),"Tell whether given point lies within this cylinder, still having 'pad' space to the solid boundary").def("aabb",&inCylinder::aabb,"Return minimum and maximum values for AABB");
+}
+
Added: trunk/gui/py/pack.py
===================================================================
--- trunk/gui/py/pack.py 2009-06-09 11:23:06 UTC (rev 1793)
+++ trunk/gui/py/pack.py 2009-06-13 15:09:51 UTC (rev 1794)
@@ -0,0 +1,34 @@
+
+from numpy import arange; from math import sqrt
+import itertools
+from yade import utils
+
+# make predicates available from yade.pack
+from _packPredicates import *
+
+def regularOrtho(predicate,radius,gap,**kw):
+ """Return set of spheres in regular orthogonal grid, clipped inside solid given by predicate.
+ Created spheres will have given radius and will be separated by gap space."""
+ ret=[]
+ mn,mx=predicate.aabb()
+ xx,yy,zz=[arange(mn[i]+radius,mx[i]-radius,2*radius+gap) for i in 0,1,2]
+ for xyz in itertools.product(xx,yy,zz):
+ if predicate(xyz,radius): ret+=[utils.sphere(xyz,radius=radius,**kw)]
+ return ret
+
+def regularHexa(predicate,radius,gap,**kw):
+ """Return set of spheres in regular hexagonal grid, clipped inside solid given by predicate.
+ Created spheres will have given radius and will be separated by gap space."""
+ ret=[]
+ a=2*radius+gap
+ h=a*sqrt(3)/2.
+ mn,mx=predicate.aabb()
+ dim=[mx[i]-mn[i] for i in 0,1,2]
+ ii,jj,kk=[range(0,int(2*dim[0]/a)+1),range(0,int(dim[1]/h)+1),range(0,int(dim[2]/h)+1)]
+ for i,j,k in itertools.product(ii,jj,kk):
+ x,y,z=mn[0]+radius+i*a,mn[1]+radius+j*h,mn[2]+radius+k*h
+ if j%2==0: x+= a/2. if k%2==0 else -a/2.
+ if k%2!=0: x+=a/2.; y+=h/2.
+ if predicate((x,y,z),radius): ret+=[utils.sphere((x,y,z),radius=radius,**kw)]
+ return ret
+
Modified: trunk/scripts/test/regular-sphere-pack.py
===================================================================
--- trunk/scripts/test/regular-sphere-pack.py 2009-06-09 11:23:06 UTC (rev 1793)
+++ trunk/scripts/test/regular-sphere-pack.py 2009-06-13 15:09:51 UTC (rev 1794)
@@ -1,3 +1,43 @@
-O.bodies.append(utils.regularSphereOrthoPack([0,0,0],extents=[2,2,2],radius=.1,gap=.1,color=(1,0,1)))
-O.bodies.append(utils.regularSphereOrthoPack([0,0,4],extents=2,radius=.1,gap=.1,color=(0,1,0),velocity=[0,10,0]))
+from yade import pack
+rad,gap=.15,.02
+rho=1e3
+kw={'density':rho}
+
+O.bodies.append(
+ pack.regularHexa(pack.inSphere((0,0,4),2),radius=rad,gap=gap,color=(0,1,0),density=10*rho) # head
+ +[utils.sphere((.8,1.9,5),radius=.2,color=(.6,.6,.6),density=rho),utils.sphere((-.8,1.9,5),radius=.2,color=(.6,.6,.6),density=rho),utils.sphere((0,2.4,4),radius=.4,color=(1,0,0),density=rho)] # eyes and nose
+ +[utils.facet([(12,0,-6),(0,12,-6,),(-12,-12,-6)],dynamic=False)] # ground
+)
+
+for part in [
+ pack.regularHexa (pack.inAlignedBox((-2,-2,-2),(2,2,2)),radius=1.5*rad,gap=2*gap,color=(1,0,1),**kw), # body,
+ pack.regularOrtho(pack.inCylinder((-1,0,-2),(-1,0,-6),1),radius=rad,gap=gap,color=(0,1,1),**kw), # left leg
+ pack.regularHexa (pack.inCylinder((+1,1,-2.5),(0,3,-5),1),radius=rad,gap=gap,color=(0,1,1),**kw), # right leg
+ pack.regularHexa (pack.inCylinder((+2,0,1),(+6,0,1),1),radius=rad,gap=gap,color=(0,0,1),**kw), # right hand
+ pack.regularOrtho(pack.inCylinder((-2,0,2),(-5,0,4),1),radius=rad,gap=gap,color=(0,0,1),**kw) # left hand
+ ]: O.bodies.appendClumped(part)
+
+
+
+try:
+ from yade import qt
+ qt.Controller()
+except ImportError: pass
+
+O.engines=[
+ BexResetter(),
+ BoundingVolumeMetaEngine([InteractingSphere2AABB(),InteractingFacet2AABB(),MetaInteractingGeometry2AABB()]),
+ InsertionSortCollider(),
+ InteractionDispatchers(
+ [ef2_Sphere_Sphere_Dem3DofGeom(),ef2_Facet_Sphere_Dem3DofGeom()],
+ [SimpleElasticRelationships()],
+ [ef2_Dem3Dof_Elastic_ElasticLaw()],
+ ),
+ GravityEngine(gravity=(1e-2,1e-2,-1000)),
+ NewtonsDampedLaw(damping=.1)
+]
+# we don't care about physical accuracy here, over-critical step is fine as long as the simulation doesn't explode
+O.dt=2*utils.PWaveTimeStep()
+O.saveTmp('init')
+#for b in O.bodies:
