yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01360
[svn] r1826 - in trunk: . gui gui/py lib pkg/dem/Engine/EngineUnit pkg/dem/PreProcessor scripts/test
Author: eudoxos
Date: 2009-06-30 22:10:30 +0200 (Tue, 30 Jun 2009)
New Revision: 1826
Modified:
trunk/SConstruct
trunk/gui/SConscript
trunk/gui/py/_packPredicates.cpp
trunk/gui/py/pack.py
trunk/lib/SConscript
trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingBox4InteractionOfMyTetrahedron.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingMyTetrahedron4InteractionOfMyTetrahedron.cpp
trunk/pkg/dem/PreProcessor/TetrahedronsTest.cpp
trunk/scripts/test/gts-operators.py
Log:
1. Fix crashers in InteractingMyTetrahedron* classes (introduced by me when updating interaction logic)
2. Reimplement inGtsSurface in c/c++ which makes it faster due to pygts interface limitations (addressed on pygts mailing list meanwhile), but is dirty programming.
3. Set CCOMSTR to correspond to CXXCOMSTR in SConstruct (to show nice line when compiling c source file, for pygts)
4. Set good-looking timestep in the TetrahedronsTest generator
Modified: trunk/SConstruct
===================================================================
--- trunk/SConstruct 2009-06-29 22:09:04 UTC (rev 1825)
+++ trunk/SConstruct 2009-06-30 20:10:30 UTC (rev 1826)
@@ -368,9 +368,10 @@
if env['pretty']:
## http://www.scons.org/wiki/HidingCommandLinesInOutput
env.Replace(CXXCOMSTR='C ${SOURCES}', # → ${TARGET.file}')
- SHCXXCOMSTR='C ${SOURCES}', #→ ${TARGET.file}')
- SHLINKCOMSTR='L ${TARGET.file}', # → ${TARGET.file}')
- LINKCOMSTR='L ${TARGET.file}', # → ${TARGET.file}')
+ CCOMSTR='C ${SOURCES}',
+ SHCXXCOMSTR='C ${SOURCES}',
+ SHLINKCOMSTR='L ${TARGET.file}',
+ LINKCOMSTR='L ${TARGET.file}',
INSTALLSTR='⇒ $TARGET',
QT_UICCOMSTR='U ${SOURCES}',
QT_MOCCOMSTR='M ${SOURCES}')
Modified: trunk/gui/SConscript
===================================================================
--- trunk/gui/SConscript 2009-06-29 22:09:04 UTC (rev 1825)
+++ trunk/gui/SConscript 2009-06-30 20:10:30 UTC (rev 1826)
@@ -73,9 +73,12 @@
'Clump'
],
),
- env.SharedLibrary('_utils',['py/_utils.cpp'],SHLIBPREFIX='',
- LIBS=env['LIBS']+['Shop','ConcretePM']),
- env.SharedLibrary('_packPredicates',['py/_packPredicates.cpp'],SHLIBPREFIX=''),
+ env.SharedLibrary('_utils',['py/_utils.cpp'],SHLIBPREFIX='',LIBS=env['LIBS']+['Shop','ConcretePM']),
+ env.SharedLibrary('_packPredicates',['py/_packPredicates.cpp'],SHLIBPREFIX='',
+ # if we compile with GTS, link to the python module, as inGtsSurface uses some of its symbols.
+ # because the module doesn't have the lib- suffix, we put it directly to SHLINKFLAGS
+ # using the -l: syntax (see man ld) and declare the dependency below
+ SHLINKFLAGS=env['SHLINKFLAGS']+(['-l:$PREFIX/lib/yade$SUFFIX/py/gts/_gts.so'] if 'GTS' in env['features'] else [])),
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'),
@@ -90,5 +93,7 @@
env.File('pack.py','py'),
])
if os.path.exists('../../brefcom-mm.hh'): Depends('py/_eudoxos.cpp','../../brefcom-mm.hh')
+ # see comments for _packPredicates above
+ if 'GTS' in env['features']: env.Depends('_packPredicates.so','$PREFIX/lib/yade$SUFFIX/py/gts/_gts.so')
Modified: trunk/gui/py/_packPredicates.cpp
===================================================================
--- trunk/gui/py/_packPredicates.cpp 2009-06-29 22:09:04 UTC (rev 1825)
+++ trunk/gui/py/_packPredicates.cpp 2009-06-30 20:10:30 UTC (rev 1826)
@@ -269,7 +269,61 @@
return vvec2ttuple(Vector3r(-inf,-inf,-inf),Vector3r(inf,inf,inf)); }
};
+#ifdef YADE_GTS
+extern "C" {
+#include<yade/lib-py/pygts.h>
+}
+/* Helper function for inGtsSurface::aabb() */
+static void vertex_aabb(GtsVertex *vertex, pair<Vector3r,Vector3r> *bb)
+{
+ GtsPoint *_p=GTS_POINT(vertex);
+ Vector3r p(_p->x,_p->y,_p->z);
+ bb->first=componentMinVector(bb->first,p);
+ bb->second=componentMaxVector(bb->second,p);
+}
+/*
+This class plays tricks getting aroung pyGTS to get GTS objects and cache bb tree to speed
+up point inclusion tests. For this reason, we have to link with _gts.so (see corresponding
+SConscript file), which is at the same time the python module.
+*/
+class inGtsSurface: public Predicate{
+ python::object pySurf; // to hold the reference so that surf is valid
+ GtsSurface *surf;
+ bool is_open, noPad, noPadWarned;
+ GNode* tree;
+public:
+ inGtsSurface(python::object _surf, bool _noPad=false): pySurf(_surf), noPad(_noPad), noPadWarned(false) {
+ if(!pygts_surface_check(_surf.ptr())) throw invalid_argument("Ctor must receive a gts.Surface() instance.");
+ surf=PYGTS_SURFACE_AS_GTS_SURFACE(PYGTS_SURFACE(_surf.ptr()));
+ if(!gts_surface_is_closed(surf)) throw invalid_argument("Surface is not closed.");
+ is_open=gts_surface_volume(surf)<0.;
+ if((tree=gts_bb_tree_surface(surf))==NULL) throw runtime_error("Could not create GTree.");
+ }
+ ~inGtsSurface(){g_node_destroy(tree);}
+ python::tuple aabb() const {
+ Real inf=std::numeric_limits<Real>::infinity();
+ pair<Vector3r,Vector3r> bb; bb.first=Vector3r(inf,inf,inf); bb.second=Vector3r(-inf,-inf,-inf);
+ gts_surface_foreach_vertex(surf,(GtsFunc)vertex_aabb,&bb);
+ return vvec2ttuple(bb.first,bb.second);
+ }
+ bool ptCheck(Vector3r pt) const{
+ GtsPoint gp; gp.x=pt[0]; gp.y=pt[1]; gp.z=pt[2];
+ return (bool)gts_point_is_inside_surface(&gp,tree,is_open);
+ }
+ bool operator()(python::tuple _pt, Real pad=0.) const {
+ Vector3r pt=tuple2vec(_pt);
+ if(noPad){
+ if(pad!=0. && noPadWarned) LOG_WARN("inGtsSurface constructed with noPad; requested non-zero pad set to zero.");
+ return ptCheck(pt);
+ }
+ return ptCheck(pt) && ptCheck(pt-Vector3r(pad,0,0)) && ptCheck(pt+Vector3r(pad,0,0)) && ptCheck(pt-Vector3r(0,pad,0))&& ptCheck(pt+Vector3r(0,pad,0)) && ptCheck(pt-Vector3r(0,0,pad))&& ptCheck(pt+Vector3r(0,0,pad));
+ }
+};
+
+#endif
+
+
BOOST_PYTHON_MODULE(_packPredicates){
// base predicate class
python::class_<PredicateWrap,/* necessary, as methods are pure virtual*/ boost::noncopyable>("Predicate")
@@ -290,5 +344,8 @@
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."));
+ #ifdef YADE_GTS
+ python::class_<inGtsSurface,python::bases<Predicate> >("inGtsSurface","GTS surface predicate",python::init<python::object,python::optional<bool> >(python::args("surface","noPad"),"Ctor taking a gts.Surface() instance, which must not be modified during instance lifetime.\nThe optional noPad can disable padding (if set to True), which speeds up calls several times.\nNote: padding checks inclusion of 6 points along +- cardinal directions in the pad distance from given point, which is not exact."));
+ #endif
}
Modified: trunk/gui/py/pack.py
===================================================================
--- trunk/gui/py/pack.py 2009-06-29 22:09:04 UTC (rev 1825)
+++ trunk/gui/py/pack.py 2009-06-30 20:10:30 UTC (rev 1826)
@@ -13,9 +13,16 @@
# make c++ predicates available in this module
from _packPredicates import *
+class inGtsSurface_py(Predicate):
+ """This class was re-implemented in c++, but should stay here to serve as reference for implementing
+ Predicates in pure python code. C++ allows us to play dirty tricks in GTS which are not accessible
+ through pygts itself; the performance penalty of pygts comes from fact that if constructs and destructs
+ bb tree for the surface at every invocation of gts.Point().is_inside(). That is cached in the c++ code,
+ provided that the surface is not manipulated with during lifetime of the object (user's responsibility).
-class inGtsSurface(Predicate):
- """Predicate for GTS surfaces. Constructed using an already existing surfaces, which must be closed.
+ ---
+
+ Predicate for GTS surfaces. Constructed using an already existing surfaces, which must be closed.
import gts
surf=gts.read(open('horse.gts'))
@@ -25,13 +32,13 @@
must be enabled in the ctor by saying doSlowPad=True. If it is not enabled and pad is not zero,
warning is issued.
"""
- def __init__(self,surf,doSlowPad=False):
+ def __init__(self,surf,noPad=False):
# 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
- self.doSlowPad=doSlowPad
+ self.noPad=noPad
inf=float('inf')
mn,mx=[inf,inf,inf],[-inf,-inf,-inf]
for v in surf.vertices():
@@ -42,8 +49,8 @@
def aabb(self): return self.mn,self.mx
def __call__(self,_pt,pad=0.):
p=gts.Point(*_pt)
- if not self.doSlowPad:
- if pad!=0: warnings.warn("Pad distance not supported for GTS surfaces, using 0 instead.")
+ if self.noPad:
+ if pad!=0: warnings.warn("Padding disabled in ctor, using 0 instead.")
return p.is_inside(self.surf)
pp=[gts.Point(_pt[0]-pad,_pt[1],_pt[2]),gts.Point(_pt[0]+pad,_pt[1],_pt[2]),gts.Point(_pt[0],_pt[1]-pad,_pt[2]),gts.Point(_pt[0],_pt[1]+pad,_pt[2]),gts.Point(_pt[0],_pt[1],_pt[2]-pad),gts.Point(_pt[0],_pt[1],_pt[2]+pad)]
return p.is_inside(self.surf) and pp[0].is_inside(self.surf) and pp[1].is_inside(self.surf) and pp[2].is_inside(self.surf) and pp[3].is_inside(self.surf) and pp[4].is_inside(self.surf) and pp[5].is_inside(self.surf)
Modified: trunk/lib/SConscript
===================================================================
--- trunk/lib/SConscript 2009-06-29 22:09:04 UTC (rev 1825)
+++ trunk/lib/SConscript 2009-06-30 20:10:30 UTC (rev 1826)
@@ -79,7 +79,12 @@
env.File('py/pygts-0.3.1/__init__.py'),
env.File('py/pygts-0.3.1/pygts.py')
])
+ #env.Install('$PREFIX/lib/yade$SUFFIX/lib',[
+ # # the same, but to link the pack module with
+ # env.SharedLibrary('_gts',['py/pygts-0.3.1/cleanup.c','py/pygts-0.3.1/edge.c','py/pygts-0.3.1/face.c','py/pygts-0.3.1/object.c','py/pygts-0.3.1/point.c','py/pygts-0.3.1/pygts.c','py/pygts-0.3.1/segment.c','py/pygts-0.3.1/surface.c','py/pygts-0.3.1/triangle.c','py/pygts-0.3.1/vertex.c'],CPPDEFINES=env['CPPDEFINES']+['PYGTS_HAS_NUMPY']),
+ #])
+
env.Install('$PREFIX/lib/yade$SUFFIX/lib',[
env.SharedLibrary('yade-base',
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingBox4InteractionOfMyTetrahedron.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingBox4InteractionOfMyTetrahedron.cpp 2009-06-29 22:09:04 UTC (rev 1825)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingBox4InteractionOfMyTetrahedron.cpp 2009-06-30 20:10:30 UTC (rev 1826)
@@ -74,7 +74,7 @@
shared_ptr<InteractionOfMyTetrahedron> imt;
// depending whether it's a new interaction: create new one, or use the existing one.
- if (c->interactionGeometry)
+ if (!c->interactionGeometry)
imt = shared_ptr<InteractionOfMyTetrahedron>(new InteractionOfMyTetrahedron());
else
imt = YADE_PTR_CAST<InteractionOfMyTetrahedron>(c->interactionGeometry);
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingMyTetrahedron4InteractionOfMyTetrahedron.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingMyTetrahedron4InteractionOfMyTetrahedron.cpp 2009-06-29 22:09:04 UTC (rev 1825)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingMyTetrahedron4InteractionOfMyTetrahedron.cpp 2009-06-30 20:10:30 UTC (rev 1826)
@@ -36,10 +36,10 @@
shared_ptr<InteractionOfMyTetrahedron> imt;
// depending whether it's a new interaction: create new one, or use the existing one.
- if (c->interactionGeometry)
+ if (!c->interactionGeometry)
imt = shared_ptr<InteractionOfMyTetrahedron>(new InteractionOfMyTetrahedron());
else
- imt = YADE_PTR_CAST<InteractionOfMyTetrahedron>(c->interactionGeometry);
+ imt = YADE_PTR_CAST<InteractionOfMyTetrahedron>(c->interactionGeometry);
bool isInteracting = false;
for(int i=0 ; i<4 ; ++i )
@@ -66,7 +66,6 @@
isInteracting = true;
}
-
c->interactionGeometry = imt;
return isInteracting;
}
Modified: trunk/pkg/dem/PreProcessor/TetrahedronsTest.cpp
===================================================================
--- trunk/pkg/dem/PreProcessor/TetrahedronsTest.cpp 2009-06-29 22:09:04 UTC (rev 1825)
+++ trunk/pkg/dem/PreProcessor/TetrahedronsTest.cpp 2009-06-30 20:10:30 UTC (rev 1826)
@@ -153,8 +153,10 @@
setProgress(current++/all);
}
}
+
+ rootBody->dt=1e-2;
- message="foo bar "+boost::lexical_cast<std::string>(42);
+ message="";
return true;
}
Modified: trunk/scripts/test/gts-operators.py
===================================================================
--- trunk/scripts/test/gts-operators.py 2009-06-29 22:09:04 UTC (rev 1825)
+++ trunk/scripts/test/gts-operators.py 2009-06-30 20:10:30 UTC (rev 1826)
@@ -3,7 +3,9 @@
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.
+The disadvantage of the predicate union | is that each sphere must fit whole in one
+surface or another: with padding, several points on the sphere are tested. Therefore,
+areas near both surfaces' boundary will not be filled at all.
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.
@@ -17,7 +19,7 @@
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
+radius=0.002
O.bodies.append(pack.gtsSurface2Facets(s12,color=(0,0,1)))
qt.View()
@@ -29,4 +31,3 @@
O.bodies.append(pack.regularHexa(pack.inGtsSurface(s12),radius,gap=0.,color=(1,0,0)))
t2=time()
print 'Using surface union: %gs'%(t2-t1)
-
