yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01358
[svn] r1824 - in trunk: gui/py scripts/test
Author: eudoxos
Date: 2009-06-29 23:06:22 +0200 (Mon, 29 Jun 2009)
New Revision: 1824
Added:
trunk/scripts/test/gts-horse.py
trunk/scripts/test/gts-operators.py
Removed:
trunk/scripts/test/gts.py
Modified:
trunk/gui/py/_packPredicates.cpp
trunk/gui/py/_utils.cpp
trunk/gui/py/pack.py
Log:
1. Make virtual methods on pack.Predicate work in python (inGtsSurface, for instance). Add scripts/test/gts-operators.py to show that.
2. Move scripts/test/gts.py to gts-horse.py
3. Add a few functions to pack for constructing triangulated revolution surfaces from meridians.
Modified: trunk/gui/py/_packPredicates.cpp
===================================================================
--- trunk/gui/py/_packPredicates.cpp 2009-06-29 14:44:29 UTC (rev 1823)
+++ trunk/gui/py/_packPredicates.cpp 2009-06-29 21:06:22 UTC (rev 1824)
@@ -4,6 +4,7 @@
#include<yade/lib-base/Logging.hpp>
#include<yade/lib-base/yadeWm3.hpp>
#include<yade/lib-base/yadeWm3Extra.hpp>
+// #include<yade/gui-py/_utils.hpp> // will be: yade/lib-py/_utils.hpp> at some point
#include<Wm3Vector3.h>
using namespace boost;
@@ -26,63 +27,88 @@
// aux functions
python::tuple vec2tuple(const Vector3r& v){return boost::python::make_tuple(v[0],v[1],v[2]);}
+python::tuple vec2tuple(const Vector2r& v){return boost::python::make_tuple(v[0],v[1]);}
Vector3r tuple2vec(const python::tuple& t){return Vector3r(python::extract<double>(t[0])(),python::extract<double>(t[1])(),python::extract<double>(t[2])());}
+Vector2r tuple2vec2d(const python::tuple& t){return Vector2r(python::extract<double>(t[0])(),python::extract<double>(t[1])());}
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)); }
-class Predicate{
+struct Predicate{
public:
- virtual bool operator() (python::tuple pt,Real pad=0.) const {throw logic_error("Calling virtual operator() of an abstract class Predicate.");}
- virtual python::tuple aabb() const {throw logic_error("Calling virtual aabb() of an abstract class Predicate.");}
+ virtual bool operator() (python::tuple pt,Real pad=0.) const = 0;
+ virtual python::tuple aabb() const = 0;
};
// make the pad parameter optional
BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(PredicateCall_overloads,operator(),1,2);
+/* Since we want to make Predicate::operator() and Predicate::aabb() callable from c++ on python::object
+with the right virtual method resolution, we have to wrap the class in the following way. See
+http://www.boost.org/doc/libs/1_38_0/libs/python/doc/tutorial/doc/html/python/exposing.html for documentation
+on exposing virtual methods.
+
+This makes it possible to derive a python class from Predicate, override its aabb() method, for instance,
+and use it in PredicateUnion, which will call the python implementation of aabb() as it should. This
+approach is used in the inGtsSurface class defined in pack.py.
+
+See scripts/test/gts-operators.py for an example.
+
+NOTE: you still have to call base class ctor in your class' ctor derived in python, e.g.
+super(inGtsSurface,self).__init__() so that virtual methods work as expected.
+*/
+struct PredicateWrap: Predicate, python::wrapper<Predicate>{
+ bool operator()(python::tuple pt, Real pad=0.) const { return this->get_override("__call__")(pt,pad);}
+ python::tuple aabb() const { return this->get_override("aabb")(); }
+};
+// make the pad parameter optional
+BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(PredicateWrapCall_overloads,operator(),1,2);
+
/*********************************************************************************
****************** Boolean operations on predicates ******************************
*********************************************************************************/
+const Predicate& obj2pred(python::object obj){ return python::extract<const Predicate&>(obj)();}
+
class PredicateBoolean: public Predicate{
protected:
- const shared_ptr<Predicate> A,B;
+ const python::object A,B;
public:
- PredicateBoolean(const shared_ptr<Predicate> _A, const shared_ptr<Predicate> _B): A(_A), B(_B){}
- const shared_ptr<Predicate> getA(){ return A;}
- const shared_ptr<Predicate> getB(){ return B;}
+ PredicateBoolean(const python::object _A, const python::object _B): A(_A), B(_B){}
+ const python::object getA(){ return A;}
+ const python::object getB(){ return B;}
};
// http://www.linuxtopia.org/online_books/programming_books/python_programming/python_ch16s03.html
class PredicateUnion: public PredicateBoolean{
public:
- PredicateUnion(const shared_ptr<Predicate> _A, const shared_ptr<Predicate> _B): PredicateBoolean(_A,_B){}
- virtual bool operator()(python::tuple pt,Real pad) const {return (*A)(pt,pad)||(*B)(pt,pad);}
- virtual python::tuple aabb() const { Vector3r minA,maxA,minB,maxB; ttuple2vvec(A->aabb(),minA,maxA); ttuple2vvec(B->aabb(),minB,maxB); return vvec2ttuple(componentMinVector(minA,minB),componentMaxVector(maxA,maxB));}
+ PredicateUnion(const python::object _A, const python::object _B): PredicateBoolean(_A,_B){}
+ virtual bool operator()(python::tuple pt,Real pad) const {return obj2pred(A)(pt,pad)||obj2pred(B)(pt,pad);}
+ virtual python::tuple aabb() const { Vector3r minA,maxA,minB,maxB; ttuple2vvec(obj2pred(A).aabb(),minA,maxA); ttuple2vvec(obj2pred(B).aabb(),minB,maxB); return vvec2ttuple(componentMinVector(minA,minB),componentMaxVector(maxA,maxB));}
};
-PredicateUnion makeUnion(const shared_ptr<Predicate>& A, const shared_ptr<Predicate>& B){ return PredicateUnion(A,B);}
+PredicateUnion makeUnion(const python::object& A, const python::object& B){ return PredicateUnion(A,B);}
class PredicateIntersection: public PredicateBoolean{
public:
- PredicateIntersection(const shared_ptr<Predicate> _A, const shared_ptr<Predicate> _B): PredicateBoolean(_A,_B){}
- virtual bool operator()(python::tuple pt,Real pad) const {return (*A)(pt,pad) && (*B)(pt,pad);}
- virtual python::tuple aabb() const { Vector3r minA,maxA,minB,maxB; ttuple2vvec(A->aabb(),minA,maxA); ttuple2vvec(B->aabb(),minB,maxB); return vvec2ttuple(componentMaxVector(minA,minB),componentMinVector(maxA,maxB));}
+ PredicateIntersection(const python::object _A, const python::object _B): PredicateBoolean(_A,_B){}
+ virtual bool operator()(python::tuple pt,Real pad) const {return obj2pred(A)(pt,pad) && obj2pred(B)(pt,pad);}
+ virtual python::tuple aabb() const { Vector3r minA,maxA,minB,maxB; ttuple2vvec(obj2pred(A).aabb(),minA,maxA); ttuple2vvec(obj2pred(B).aabb(),minB,maxB); return vvec2ttuple(componentMaxVector(minA,minB),componentMinVector(maxA,maxB));}
};
-PredicateIntersection makeIntersection(const shared_ptr<Predicate>& A, const shared_ptr<Predicate>& B){ return PredicateIntersection(A,B);}
+PredicateIntersection makeIntersection(const python::object& A, const python::object& B){ return PredicateIntersection(A,B);}
class PredicateDifference: public PredicateBoolean{
public:
- PredicateDifference(const shared_ptr<Predicate> _A, const shared_ptr<Predicate> _B): PredicateBoolean(_A,_B){}
- virtual bool operator()(python::tuple pt,Real pad) const {return (*A)(pt,pad) && !(*B)(pt,-pad);}
- virtual python::tuple aabb() const { return A->aabb(); }
+ PredicateDifference(const python::object _A, const python::object _B): PredicateBoolean(_A,_B){}
+ virtual bool operator()(python::tuple pt,Real pad) const {return obj2pred(A)(pt,pad) && !obj2pred(B)(pt,-pad);}
+ virtual python::tuple aabb() const { return obj2pred(A).aabb(); }
};
-PredicateDifference makeDifference(const shared_ptr<Predicate>& A, const shared_ptr<Predicate>& B){ return PredicateDifference(A,B);}
+PredicateDifference makeDifference(const python::object& A, const python::object& B){ return PredicateDifference(A,B);}
class PredicateSymmetricDifference: public PredicateBoolean{
public:
- PredicateSymmetricDifference(const shared_ptr<Predicate> _A, const shared_ptr<Predicate> _B): PredicateBoolean(_A,_B){}
- virtual bool operator()(python::tuple pt,Real pad) const {bool inA=(*A)(pt,pad), inB=(*B)(pt,pad); return (inA && !inB) || (!inA && inB);}
- virtual python::tuple aabb() const { Vector3r minA,maxA,minB,maxB; ttuple2vvec(A->aabb(),minA,maxA); ttuple2vvec(B->aabb(),minB,maxB); return vvec2ttuple(componentMinVector(minA,minB),componentMaxVector(maxA,maxB));}
+ PredicateSymmetricDifference(const python::object _A, const python::object _B): PredicateBoolean(_A,_B){}
+ virtual bool operator()(python::tuple pt,Real pad) const {bool inA=obj2pred(A)(pt,pad), inB=obj2pred(B)(pt,pad); return (inA && !inB) || (!inA && inB);}
+ virtual python::tuple aabb() const { Vector3r minA,maxA,minB,maxB; ttuple2vvec(obj2pred(A).aabb(),minA,maxA); ttuple2vvec(obj2pred(B).aabb(),minB,maxB); return vvec2ttuple(componentMinVector(minA,minB),componentMaxVector(maxA,maxB));}
};
-PredicateSymmetricDifference makeSymmetricDifference(const shared_ptr<Predicate>& A, const shared_ptr<Predicate>& B){ return PredicateSymmetricDifference(A,B);}
+PredicateSymmetricDifference makeSymmetricDifference(const python::object& A, const python::object& B){ return PredicateSymmetricDifference(A,B);}
/*********************************************************************************
****************************** Primitive predicates ******************************
@@ -237,7 +263,7 @@
// between both notch planes, closer to the edge than pad (distInPlane<pad)
return false;
}
- // aabb here doesn't make any sense since we are negated. Return just the center point.
+ // This predicate is not bounded, return infinities
python::tuple aabb() const {
Real inf=std::numeric_limits<Real>::infinity();
return vvec2ttuple(Vector3r(-inf,-inf,-inf),Vector3r(inf,inf,inf)); }
@@ -245,22 +271,24 @@
BOOST_PYTHON_MODULE(_packPredicates){
- python::class_<Predicate, shared_ptr<Predicate> >("Predicate")
- .def("__call__",&Predicate::operator(),PredicateCall_overloads(python::args("point","padding"),"Tell whether given point lies within this sphere, still having 'pad' space to the solid boundary"))
- .def("aabb",&Predicate::aabb,"Return minimum and maximum values for AABB")
+ // base predicate class
+ python::class_<PredicateWrap,/* necessary, as methods are pure virtual*/ boost::noncopyable>("Predicate")
+ .def("__call__",python::pure_virtual(&Predicate::operator()))
+ .def("aabb",python::pure_virtual(&Predicate::aabb))
.def("__or__",makeUnion).def("__and__",makeIntersection).def("__sub__",makeDifference).def("__xor__",makeSymmetricDifference);
- python::class_<PredicateBoolean,python::bases<Predicate> >("PredicateBoolean","Boolean operation on 2 predicates (abstract class)",python::no_init)
+ // boolean operations
+ python::class_<PredicateBoolean,python::bases<Predicate>,boost::noncopyable>("PredicateBoolean","Boolean operation on 2 predicates (abstract class)",python::no_init)
.add_property("A",&PredicateBoolean::getA).add_property("B",&PredicateBoolean::getB);
- python::class_<PredicateUnion,python::bases<PredicateBoolean> >("PredicateUnion","Union of 2 predicates",python::init<shared_ptr<Predicate>,shared_ptr<Predicate> >());
- python::class_<PredicateIntersection,python::bases<PredicateBoolean> >("PredicateIntersection","Intersection of 2 predicates",python::init<shared_ptr<Predicate>,shared_ptr<Predicate> >());
- python::class_<PredicateDifference,python::bases<PredicateBoolean> >("PredicateDifference","Difference of 2 predicates",python::init<shared_ptr<Predicate>,shared_ptr<Predicate> >());
- python::class_<PredicateSymmetricDifference,python::bases<PredicateBoolean> >("PredicateSymmetricDifference","SymmetricDifference of 2 predicates",python::init<shared_ptr<Predicate>,shared_ptr<Predicate> >());
+ python::class_<PredicateUnion,python::bases<PredicateBoolean> >("PredicateUnion","Union of 2 predicates",python::init<python::object,python::object>());
+ python::class_<PredicateIntersection,python::bases<PredicateBoolean> >("PredicateIntersection","Intersection of 2 predicates",python::init<python::object,python::object >());
+ python::class_<PredicateDifference,python::bases<PredicateBoolean> >("PredicateDifference","Difference of 2 predicates",python::init<python::object,python::object >());
+ python::class_<PredicateSymmetricDifference,python::bases<PredicateBoolean> >("PredicateSymmetricDifference","SymmetricDifference of 2 predicates",python::init<python::object,python::object >());
+ // primitive predicates
python::class_<inSphere,python::bases<Predicate> >("inSphere","Sphere predicate.",python::init<python::tuple,Real>(python::args("center","radius"),"Ctor taking center (as a 3-tuple) and radius"));
python::class_<inAlignedBox,python::bases<Predicate> >("inAlignedBox","Axis-aligned box predicate",python::init<python::tuple,python::tuple>(python::args("minAABB","maxAABB"),"Ctor taking minumum and maximum points of the box (as 3-tuples)."));
python::class_<inCylinder,python::bases<Predicate> >("inCylinder","Cylinder predicate",python::init<python::tuple,python::tuple,Real>(python::args("centerBottom","centerTop","radius"),"Ctor taking centers of the lateral walls (as 3-tuples) and radius."));
python::class_<inHyperboloid,python::bases<Predicate> >("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)."));
python::class_<inEllipsoid,python::bases<Predicate> >("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)."));
python::class_<notInNotch,python::bases<Predicate> >("notInNotch","Outside of infinite, rectangle-shaped notch predicate",python::init<python::tuple,python::tuple,python::tuple,Real>(python::args("centerPoint","edge","normal","aperture"),"Ctor taking point in the symmetry plane, vector pointing along the edge, plane normal and aperture size.\nThe side inside the notch is edge×normal.\nNormal is made perpendicular to the edge.\nAll vectors are normalized at construction time."));
-
}
Modified: trunk/gui/py/_utils.cpp
===================================================================
--- trunk/gui/py/_utils.cpp 2009-06-29 14:44:29 UTC (rev 1823)
+++ trunk/gui/py/_utils.cpp 2009-06-29 21:06:22 UTC (rev 1824)
@@ -13,8 +13,10 @@
#include<numpy/ndarrayobject.h>
+// #include"_utils.hpp"
+
using namespace boost::python;
#ifdef LOG4CXX
Modified: trunk/gui/py/pack.py
===================================================================
--- trunk/gui/py/pack.py 2009-06-29 14:44:29 UTC (rev 1823)
+++ trunk/gui/py/pack.py 2009-06-29 21:06:22 UTC (rev 1824)
@@ -24,7 +24,9 @@
Note: padding is not supported, warning is given and zero used instead.
"""
def __init__(self,surf):
- import gts
+ # call base class ctor; necessary for virtual methods to work as expected.
+ # see comments in _packPredicates.cpp for struct PredicateWrap.
+ super(inGtsSurface,self).__init__()
if not surf.is_closed(): raise RuntimeError("Surface for inGtsSurface predicate must be closed.")
self.surf=surf
inf=float('inf')
@@ -33,6 +35,7 @@
c=v.coords()
mn,mx=[min(mn[i],c[i]) for i in 0,1,2],[max(mx[i],c[i]) for i in 0,1,2]
self.mn,self.mx=tuple(mn),tuple(mx)
+ import gts
def aabb(self): return self.mn,self.mx
def __call__(self,_pt,pad=0.):
p=gts.Point(*_pt)
@@ -43,8 +46,52 @@
"""Construct facets from given GTS surface. **kw is passed to utils.facet."""
return [utils.facet([v.coords() for v in face.vertices()],**kw) for face in surf]
+def sweptPolylines2gtsSurface(pts,threshold=0,capStart=False,capEnd=False):
+ """Create swept suface (as GTS triangulation) given same-length sequences of points (as 3-tuples).
+ If threshold is given (>0), gts.Surface().cleanup(threshold) will be called before returning.
+ This removes vrtices closer than threshold. Can be used to create closed swept surface (revolved), as
+ we don't check for coincident vertices otherwise.
+ """
+ if not len(set([len(pts1) for pts1 in pts]))==1: raise RuntimeError("Polylines must be all of the same length!")
+ vtxs=[[gts.Vertex(x,y,z) for x,y,z in pts1] for pts1 in pts]
+ sectEdges=[[gts.Edge(vtx[i],vtx[i+1]) for i in xrange(0,len(vtx)-1)] for vtx in vtxs]
+ interSectEdges=[[] for i in range(0,len(vtxs)-1)]
+ for i in range(0,len(vtxs)-1):
+ for j in range(0,len(vtxs[i])):
+ interSectEdges[i].append(gts.Edge(vtxs[i][j],vtxs[i+1][j]))
+ if j<len(vtxs[i])-1: interSectEdges[i].append(gts.Edge(vtxs[i][j],vtxs[i+1][j+1]))
+ surf=gts.Surface()
+ for i in range(0,len(vtxs)-1):
+ for j in range(0,len(vtxs[i])-1):
+ surf.add(gts.Face(interSectEdges[i][2*j+1],sectEdges[i+1][j],interSectEdges[i][2*j]))
+ surf.add(gts.Face(sectEdges[i][j],interSectEdges[i][2*j+2],interSectEdges[i][2*j+1]))
+ def doCap(vtx,edg,start):
+ ret=[]
+ eFan=[edg[0]]+[gts.Edge(vtx[i],vtx[0]) for i in range(2,len(vtx))]
+ for i in range(1,len(edg)):
+ ret+=[gts.Face(eFan[i-1],eFan[i],edg[i]) if start else gts.Face(eFan[i-1],edg[i],eFan[i])]
+ return ret
+ caps=[]
+ if capStart: caps+=doCap(vtxs[0],sectEdges[0],start=True)
+ if capEnd: caps+=doCap(vtxs[-1],sectEdges[-1],start=False)
+ for cap in caps: surf.add(cap)
+ if threshold>0: surf.cleanup(threshold)
+ return surf
+import euclid
+def revolutionSurfaceMeridians(sects,angles,origin=euclid.Vector3(0,0,0),orientation=euclid.Quaternion()):
+ """Revolution surface given sequences of 2d points and sequence of corresponding angles,
+ returning sequences of 3d points representing meridian sections of the revolution surface.
+ The 2d sections are turned around z-axis, but they can be transformed
+ using the origin and orientation arguments to give arbitrary spiral orientation."""
+ import math
+ def toGlobal(x,y,z):
+ return tuple(origin+orientation*(euclid.Vector3(x,y,z)))
+ return [[toGlobal(x2d*math.cos(angles[i]),x2d*math.sin(angles[i]),y2d) for x2d,y2d in sects[i]] for i in range(0,len(sects))]
+
+
+
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."""
Copied: trunk/scripts/test/gts-horse.py (from rev 1822, trunk/scripts/test/gts.py)
===================================================================
--- trunk/scripts/test/gts.py 2009-06-29 08:55:00 UTC (rev 1822)
+++ trunk/scripts/test/gts-horse.py 2009-06-29 21:06:22 UTC (rev 1824)
@@ -0,0 +1,43 @@
+# encoding: utf-8
+# © 2009 Václav Šmilauer <eudoxos@xxxxxxxx>
+"""Script showing how to use GTS to import arbitrary triangulated surface,
+which can further be either filled with packing (if used as predicate) or converted
+to facets representing the surface."""
+
+from yade import pack
+import gts
+
+try:
+ #surf=gts.read(open('horse.gts')); surf.coarsen(1000); surf.write(open('horse.coarse.gts','w'))
+ surf=gts.read(open('horse.coarse.gts'))
+except IOError:
+ print """horse.gts not found, you need to download input data:
+
+ wget http://gts.sourceforge.net/samples/horse.gts.gz
+ gunzip horse.gts.gz
+ """
+ quit()
+print dir(surf)
+if surf.is_closed():
+ pred=pack.inGtsSurface(surf)
+ aabb=pred.aabb()
+ dim0=aabb[1][0]-aabb[0][0]; radius=dim0/30. # get some characteristic dimension, use it for radius
+ O.bodies.append(pack.regularHexa(pred,radius=radius,gap=radius/4.,density=2000))
+ surf.translate(0,0,-(aabb[1][2]-aabb[0][2])) # move surface down so that facets are underneath the falling spheres
+O.bodies.append(pack.gtsSurface2Facets(surf,wire=True))
+
+O.engines=[
+ BexResetter(),
+ BoundingVolumeMetaEngine([InteractingSphere2AABB(),InteractingFacet2AABB(),MetaInteractingGeometry2AABB()]),
+ InsertionSortCollider(),
+ InteractionDispatchers(
+ [ef2_Sphere_Sphere_Dem3DofGeom(),ef2_Facet_Sphere_Dem3DofGeom()],
+ [SimpleElasticRelationships()],
+ [Law2_Dem3Dof_Elastic_Elastic()],
+ ),
+ GravityEngine(gravity=[0,0,-1e4]),
+ NewtonsDampedLaw(damping=.1)
+]
+
+O.dt=1.5*utils.PWaveTimeStep()
+O.saveTmp()
Property changes on: trunk/scripts/test/gts-horse.py
___________________________________________________________________
Name: svn:mergeinfo
+
Added: trunk/scripts/test/gts-operators.py
===================================================================
--- trunk/scripts/test/gts-operators.py 2009-06-29 14:44:29 UTC (rev 1823)
+++ trunk/scripts/test/gts-operators.py 2009-06-29 21:06:22 UTC (rev 1824)
@@ -0,0 +1,32 @@
+""" This file shows 2 ways to fill union of triangulated surfaces:
+You can either use union of 2 inGtsSurface predicates or
+create union of surfaces using GTS calls first and use a single
+isGtsSurface as predicate with the united surface.
+
+Surprisingly, the first variant seems faster.
+
+Note that GTS only moves references to surfaces around, therefore e.g. translating
+surface that is part of the union will move also the part of the united surface.
+Therefore, we use the copy() method for deep copy here.
+"""
+from yade import pack,qt
+import gts
+
+s1=gts.read(open('horse.coarse.gts'))
+s2=gts.Surface(); s2.copy(s1); s2.translate(0.04,0,0)
+O.bodies.append(pack.gtsSurface2Facets(s1,color=(0,1,0))+pack.gtsSurface2Facets(s2,color=(1,0,0)))
+
+s12=gts.Surface(); s12.copy(s1.union(s2)); s12.translate(0,0,.1)
+radius=0.005
+O.bodies.append(pack.gtsSurface2Facets(s12,color=(0,0,1)))
+
+qt.View()
+from time import time
+t0=time()
+O.bodies.append(pack.regularHexa(pack.inGtsSurface(s1) | pack.inGtsSurface(s2),radius,gap=0,color=(0,1,0)))
+t1=time()
+print 'Using predicate union: %gs'%(t1-t0)
+O.bodies.append(pack.regularHexa(pack.inGtsSurface(s12),radius,gap=0.,color=(1,0,0)))
+t2=time()
+print 'Using surface union: %gs'%(t2-t1)
+
Deleted: trunk/scripts/test/gts.py
===================================================================
--- trunk/scripts/test/gts.py 2009-06-29 14:44:29 UTC (rev 1823)
+++ trunk/scripts/test/gts.py 2009-06-29 21:06:22 UTC (rev 1824)
@@ -1,43 +0,0 @@
-# encoding: utf-8
-# © 2009 Václav Šmilauer <eudoxos@xxxxxxxx>
-"""Script showing how to use GTS to import arbitrary triangulated surface,
-which can further be either filled with packing (if used as predicate) or converted
-to facets representing the surface."""
-
-from yade import pack
-import gts
-
-try:
- #surf=gts.read(open('horse.gts')); surf.coarsen(1000); surf.write(open('horse.coarse.gts','w'))
- surf=gts.read(open('horse.coarse.gts'))
-except IOError:
- print """horse.gts not found, you need to download input data:
-
- wget http://gts.sourceforge.net/samples/horse.gts.gz
- gunzip horse.gts.gz
- """
- quit()
-
-if surf.is_closed():
- pred=pack.inGtsSurface(surf)
- aabb=pred.aabb()
- dim0=aabb[1][0]-aabb[0][0]; radius=dim0/30. # get some characteristic dimension, use it for radius
- O.bodies.append(pack.regularHexa(pred,radius=radius,gap=radius/4.,density=2000))
- surf.translate(0,0,-(aabb[1][2]-aabb[0][2])) # move surface down so that facets underneath
-O.bodies.append(pack.gtsSurface2Facets(surf,wire=True))
-
-O.engines=[
- BexResetter(),
- BoundingVolumeMetaEngine([InteractingSphere2AABB(),InteractingFacet2AABB(),MetaInteractingGeometry2AABB()]),
- InsertionSortCollider(),
- InteractionDispatchers(
- [ef2_Sphere_Sphere_Dem3DofGeom(),ef2_Facet_Sphere_Dem3DofGeom()],
- [SimpleElasticRelationships()],
- [Law2_Dem3Dof_Elastic_Elastic()],
- ),
- GravityEngine(gravity=[0,0,-1e4]),
- NewtonsDampedLaw(damping=.1)
-]
-
-O.dt=1.5*utils.PWaveTimeStep()
-O.saveTmp()
