yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-06-17 17:47:48 +0200 (Wed, 17 Jun 2009)
New Revision: 1803

Modified:
   trunk/gui/py/_packPredicates.cpp
   trunk/gui/py/pack.py
   trunk/scripts/test/regular-sphere-pack.py
Log:
Anton's new things: axis-aligned ellipsoid predicate and fix in pack.regularHexa.

I added elliptical leg to the sample script and mouth as cylinder. The reason is that
critical timestep depends (also) on the minimum sphere size, therefore the simulation
would run slower...



Modified: trunk/gui/py/_packPredicates.cpp
===================================================================
--- trunk/gui/py/_packPredicates.cpp	2009-06-17 12:49:09 UTC (rev 1802)
+++ trunk/gui/py/_packPredicates.cpp	2009-06-17 15:47:48 UTC (rev 1803)
@@ -110,6 +110,32 @@
 	}
 };
 
+/* Axis-aligned ellipsoid predicate */
+class inEllipsoid{
+	Vector3r c, abc;
+public:
+	inEllipsoid(python::tuple _c, python::tuple _abc)
+		{c=tuple2vec(_c); abc=tuple2vec(_abc);}
+	bool operator()(python::tuple _pt, Real pad=0.){
+		Vector3r pt=tuple2vec(_pt);
+		
+		//Define the ellipsoid X-coordinate of given Y and Z
+		Real x = sqrt((1-pow((pt[1]-c[1]),2)/((abc[1]-pad)*(abc[1]-pad))-pow((pt[2]-c[2]),2)/((abc[2]-pad)*(abc[2]-pad)))*((abc[0]-pad)*(abc[0]-pad)))+c[0]; 
+		Vector3r edgeEllipsoid(x,pt[1],pt[2]); // create a vector of these 3 coordinates
+		
+		if ((pt-c).Length()<=(edgeEllipsoid-c).Length()) { //check whether given coordinates lie inside ellipsoid or not
+			return true;
+		} else {
+			return false;
+		}
+	}
+	python::tuple aabb(){
+		Vector3r& center(c); Vector3r& ABC(abc);
+		return vvec2ttuple(Vector3r(center[0]-ABC[0],center[1]-ABC[1],center[2]-ABC[2]),Vector3r(center[0]+ABC[0],center[1]+ABC[1],center[2]+ABC[2]));
+	}
+};
+
+
 BOOST_PYTHON_MODULE(_packPredicates){
 	boost::python::class_<inSphere>("inSphere","Sphere predicate.",python::init<python::tuple,Real>(python::args("center","radius"),"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");
@@ -119,5 +145,7 @@
 		.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");
 	boost::python::class_<inHyperboloid>("inHyperboloid","Hyperboloid predicate",python::init<python::tuple,python::tuple,Real,Real>(python::args("centerBottom","centerTop","radius","skirt"),"Ctor taking centers of the lateral walls (as 3-tuples), radius at bases and skirt (middle radius)."))
 		.def("__call__",&inHyperboloid::operator(),"Tell whether given point lies within this hyperboloid, still having 'pad' space to the solid boundary\n(not accurate, since distance perpendicular to the axis, not the surface, is taken in account)").def("aabb",&inHyperboloid::aabb,"Return minimum and maximum values for AABB");
+	boost::python::class_<inEllipsoid>("inEllipsoid","Ellipsoid predicate",python::init<python::tuple,python::tuple>(python::args("centerPoint","abc"),"Ctor taking center of the ellipsoid (3-tuple) and its 3 radii (3-tuple)."))
+		.def("__call__",&inEllipsoid::operator(),"Tell whether given point lies within this inEllipsoid").def("aabb",&inEllipsoid::aabb,"Return minimum and maximum values for AABB");
 }
 

Modified: trunk/gui/py/pack.py
===================================================================
--- trunk/gui/py/pack.py	2009-06-17 12:49:09 UTC (rev 1802)
+++ trunk/gui/py/pack.py	2009-06-17 15:47:48 UTC (rev 1803)
@@ -24,7 +24,7 @@
 	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)]
+	ii,jj,kk=[range(0,int(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.

Modified: trunk/scripts/test/regular-sphere-pack.py
===================================================================
--- trunk/scripts/test/regular-sphere-pack.py	2009-06-17 12:49:09 UTC (rev 1802)
+++ trunk/scripts/test/regular-sphere-pack.py	2009-06-17 15:47:48 UTC (rev 1803)
@@ -8,13 +8,14 @@
 	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
+	+pack.regularHexa(pack.inCylinder((-1,2.2,3.3),(1,2.2,3.3),2*rad),radius=rad,gap=gap/3,color=(0.929,0.412,0.412),density=10*rho) #mouth
 )
 
 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.regularOrtho(pack.inEllipsoid((-1,0,-4),(1,1,2)),radius=rad,gap=0,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.inHyperboloid((+2,0,1),(+6,0,1),1,.5),radius=rad,gap=gap,color=(0,0,1),**kw), # right hand
+	pack.regularHexa (pack.inHyperboloid((+2,0,1),(+6,0,0),1,.5),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)