yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #00596
[svn] r1494 - in trunk: extra extra/clump gui gui/py
Author: eudoxos
Date: 2008-08-25 09:58:04 +0200 (Mon, 25 Aug 2008)
New Revision: 1494
Added:
trunk/gui/py/_utils.cpp
Modified:
trunk/extra/SimpleScene.cpp
trunk/extra/clump/Shop.cpp
trunk/extra/clump/Shop.hpp
trunk/gui/SConscript
trunk/gui/py/PythonUI_rc.py
trunk/gui/py/eudoxos.py
trunk/gui/py/ipython.py
trunk/gui/py/utils.py
Log:
1. Reimplement some loop-intensive functions from yade.utils in c++ for better performance (as yade._utils, imported into yade.utils automatically - transparent to the user)
2. Add poisson's ratio and Young's modulus estimation in yade.eudoxos based on linear regressions of position?\226?\134?\146displacement mapping along different axes
3. Rename Shop::ElasticWaveTimestepEstimate to Shop::PWaveTimeStep, wrapped in yade.utils
Modified: trunk/extra/SimpleScene.cpp
===================================================================
--- trunk/extra/SimpleScene.cpp 2008-08-22 23:01:31 UTC (rev 1493)
+++ trunk/extra/SimpleScene.cpp 2008-08-25 07:58:04 UTC (rev 1494)
@@ -156,7 +156,7 @@
rootBody->bodies->insert(sphere);
//@
- rootBody->dt=.2*Shop::ElasticWaveTimestepEstimate(rootBody);
+ rootBody->dt=.2*Shop::PWaveTimeStep(rootBody);
//@
return true;
Modified: trunk/extra/clump/Shop.cpp
===================================================================
--- trunk/extra/clump/Shop.cpp 2008-08-22 23:01:31 UTC (rev 1493)
+++ trunk/extra/clump/Shop.cpp 2008-08-25 07:58:04 UTC (rev 1494)
@@ -1071,7 +1071,7 @@
{0.370512,0.415453,0.055970,0.010546},
{-1.,-1.,-1.,-1. } /* sentinel: non-positive radius */
};
-
+#if 0
Real Shop::ElasticWaveTimestepEstimate(shared_ptr<MetaBody> rootBody){
Real minDt=std::numeric_limits<Real>::infinity();
FOREACH(const shared_ptr<Body>& b, *rootBody->bodies){
@@ -1084,6 +1084,7 @@
}
return minDt;
}
+#endif
void Shop::GLDrawLine(Vector3r from, Vector3r to, Vector3r color){
glEnable(GL_LIGHTING); glColor3v(color);
@@ -1109,9 +1110,10 @@
glPopMatrix();
}
-#if 0
-Real Shop::PWaveTimeStep(MetaBody* rb){
- dt=std::numeric_limits<Real>::infinity();
+Real Shop::PWaveTimeStep(shared_ptr<MetaBody> _rb){
+ shared_ptr<MetaBody> rb=_rb;
+ if(!rb)rb=Omega::instance().getRootBody();
+ Real dt=std::numeric_limits<Real>::infinity();
FOREACH(const shared_ptr<Body>& b, *rb->bodies){
if(!b->physicalParameters || !b->geometricalModel) continue;
shared_ptr<ElasticBodyParameters> ebp=dynamic_pointer_cast<ElasticBodyParameters>(b->physicalParameters);
@@ -1122,4 +1124,3 @@
}
return dt;
}
-#endif
Modified: trunk/extra/clump/Shop.hpp
===================================================================
--- trunk/extra/clump/Shop.hpp 2008-08-22 23:01:31 UTC (rev 1493)
+++ trunk/extra/clump/Shop.hpp 2008-08-25 07:58:04 UTC (rev 1494)
@@ -79,9 +79,6 @@
//! Save spheres in the current simulation into a text file
static void saveSpheresToFile(string fileName);
- //! Estimate timestep based on P-wave propagation speed
- static Real ElasticWaveTimestepEstimate(shared_ptr<MetaBody>);
-
static void GLDrawLine(Vector3r from, Vector3r to, Vector3r color=Vector3r(1,1,1));
static void GLDrawArrow(Vector3r from, Vector3r to, Vector3r color=Vector3r(1,1,1));
static void GLDrawText(std::string text, Vector3r pos, Vector3r color=Vector3r(1,1,1));
@@ -101,8 +98,7 @@
static Vector3r& globalStiffness(body_id_t, MetaBody* mb=NULL);
static Vector3r& globalRStiffness(body_id_t, MetaBody* mb=NULL);
};
- #if 0
- Real PWaveTimeStep(MetaBody* rb);
- #endif
+ //! Estimate timestep based on P-wave propagation speed
+ static Real PWaveTimeStep(shared_ptr<MetaBody> rb=shared_ptr<MetaBody>());
};
Modified: trunk/gui/SConscript
===================================================================
--- trunk/gui/SConscript 2008-08-22 23:01:31 UTC (rev 1493)
+++ trunk/gui/SConscript 2008-08-25 07:58:04 UTC (rev 1494)
@@ -52,7 +52,7 @@
# python modules are one level deeper so that you can say: from yade.wrapper import *
env.Install('$PREFIX/lib/yade$SUFFIX/gui/yade',[
env.SharedLibrary('wrapper',['py/yadeControl.cpp'],SHLIBPREFIX='',
- LIBS=env['LIBS']+['yade-base','libboost_python','XMLFormatManager','yade-factory','yade-serialization','Shop',
+ LIBS=env['LIBS']+['libboost_python','XMLFormatManager','yade-factory','yade-serialization','Shop',
'BoundingVolumeMetaEngine',
'GeometricalModelMetaEngine',
'InteractingGeometryMetaEngine',
@@ -63,6 +63,8 @@
'ParallelEngine'
],
),
+ env.SharedLibrary('_utils',['py/_utils.cpp'],SHLIBPREFIX='',
+ LIBS=env['LIBS']+['Shop','libboost_python']),
env.SharedLibrary('log',['py/log.cpp'],SHLIBPREFIX=''),
env.File('__init__.py','py'),
env.File('utils.py','py'),
Modified: trunk/gui/py/PythonUI_rc.py
===================================================================
--- trunk/gui/py/PythonUI_rc.py 2008-08-22 23:01:31 UTC (rev 1493)
+++ trunk/gui/py/PythonUI_rc.py 2008-08-25 07:58:04 UTC (rev 1494)
@@ -62,5 +62,3 @@
import yade.qt
yade.qt.close()
except ImportError: pass
-
-
Added: trunk/gui/py/_utils.cpp
===================================================================
--- trunk/gui/py/_utils.cpp 2008-08-22 23:01:31 UTC (rev 1493)
+++ trunk/gui/py/_utils.cpp 2008-08-25 07:58:04 UTC (rev 1494)
@@ -0,0 +1,63 @@
+#include<yade/extra/Shop.hpp>
+#include<boost/python.hpp>
+#include<yade/core/MetaBody.hpp>
+#include<yade/core/Omega.hpp>
+#include<yade/pkg-common/Sphere.hpp>
+using namespace boost::python;
+
+python::tuple vec2tuple(const Vector3r& v){return boost::python::make_tuple(v[0],v[1],v[2]);};
+
+/* \todo implement groupMask */
+python::tuple aabbExtrema(){
+ Real inf=std::numeric_limits<Real>::infinity();
+ Vector3r minimum(inf,inf,inf),maximum(-inf,-inf,-inf);
+ FOREACH(const shared_ptr<Body>& b, *Omega::instance().getRootBody()->bodies){
+ shared_ptr<Sphere> s=dynamic_pointer_cast<Sphere>(b->geometricalModel); if(!s) continue;
+ Vector3r rrr(s->radius,s->radius,s->radius);
+ minimum=componentMinVector(minimum,b->physicalParameters->se3.position-rrr);
+ maximum=componentMaxVector(maximum,b->physicalParameters->se3.position+rrr);
+ }
+ return python::make_tuple(vec2tuple(minimum),vec2tuple(maximum));
+}
+
+python::tuple negPosExtremeIds(int axis, Real distFactor=1.1){
+ python::tuple extrema=aabbExtrema();
+ Real minCoord=extract<double>(extrema[0][axis])(),maxCoord=extract<double>(extrema[1][axis])();
+ python::list minIds,maxIds;
+ FOREACH(const shared_ptr<Body>& b, *Omega::instance().getRootBody()->bodies){
+ shared_ptr<Sphere> s=dynamic_pointer_cast<Sphere>(b->geometricalModel); if(!s) continue;
+ if(b->physicalParameters->se3.position[axis]-s->radius*distFactor<=minCoord) minIds.append(b->getId());
+ if(b->physicalParameters->se3.position[axis]+s->radius*distFactor>=maxCoord) maxIds.append(b->getId());
+ }
+ return python::make_tuple(minIds,maxIds);
+}
+BOOST_PYTHON_FUNCTION_OVERLOADS(negPosExtremeIds_overloads,negPosExtremeIds,1,2);
+
+python::tuple coordsAndDisplacements(int axis){
+ python::list retCoord,retDispl;
+ FOREACH(const shared_ptr<Body>&b, *Omega::instance().getRootBody()->bodies){
+ retCoord.append(b->physicalParameters->se3.position[axis]);
+ retDispl.append(b->physicalParameters->se3.position[axis]-b->physicalParameters->refSe3.position[axis]);
+ }
+ return python::make_tuple(retCoord,retDispl);
+}
+
+void setRefSe3(){
+ FOREACH(const shared_ptr<Body>& b, *Omega::instance().getRootBody()->bodies){
+ b->physicalParameters->refSe3.position=b->physicalParameters->se3.position;
+ b->physicalParameters->refSe3.orientation=b->physicalParameters->se3.orientation;
+ }
+}
+
+
+Real PWaveTimeStep(){return Shop::PWaveTimeStep();};
+
+BOOST_PYTHON_MODULE(_utils){
+ def("PWaveTimeStep",PWaveTimeStep);
+ def("aabbExtrema",aabbExtrema);
+ def("negPosExtremeIds",negPosExtremeIds,negPosExtremeIds_overloads(args("axis","distFactor")));
+ def("coordsAndDisplacements",coordsAndDisplacements);
+ def("setRefSe3",setRefSe3);
+}
+
+
Modified: trunk/gui/py/eudoxos.py
===================================================================
--- trunk/gui/py/eudoxos.py 2008-08-22 23:01:31 UTC (rev 1493)
+++ trunk/gui/py/eudoxos.py 2008-08-25 07:58:04 UTC (rev 1494)
@@ -4,7 +4,39 @@
# I doubt there functions will be useful for anyone besides me.
#
+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
+ (slope in linear regression on discrete function particle coordinate →
+ → particle displacement in the same direction as returned by
+ utils.coordsAndDisplacements) and, (if axis '0' is the strained
+ axis) the poisson's ratio is given as -½(ε₁+ε₂)/ε₀.
+ Young's modulus is computed as σ/ε₀; if stress σ is not given (default 0),
+ the result is 0.
+ """
+ dd=[] # storage for linear regression parameters
+ import pylab,numpy,stats
+ from yade import utils
+ for axis in [0,1,2]:
+ w,dw=utils.coordsAndDisplacements(axis)
+ l,ll=stats.linregress(w,dw)[0:2] # use only tangent and section
+ dd.append((l,ll,min(w),max(w)))
+ if plot: pylab.plot(w,dw,'.',label='xyz'[axis])
+ if plot:
+ for axis in [0,1,2]:
+ dist=dd[axis][-1]-dd[axis][-2]
+ c=numpy.linspace(dd[axis][-2]-.2*dist,dd[axis][-1]+.2*dist)
+ d=[dd[axis][0]*cc+dd[axis][1] for cc in c]
+ pylab.plot(c,d,label='interp '+'xyz'[axis])
+ pylab.legend()
+ pylab.show()
+ otherAxes=(principalAxis+1)%3,(principalAxis+2)%3
+ avgTransHomogenizedStrain=.5*(dd[otherAxes[0]][0]+dd[otherAxes[1]][0])
+ principalHomogenizedStrain=dd[principalAxis][0]
+ return -avgTransHomogenizedStrain/principalHomogenizedStrain,stress/principalHomogenizedStrain
+
+
def oofemTextExport():
"""Export simulation data in text format
Modified: trunk/gui/py/ipython.py
===================================================================
--- trunk/gui/py/ipython.py 2008-08-22 23:01:31 UTC (rev 1493)
+++ trunk/gui/py/ipython.py 2008-08-25 07:58:04 UTC (rev 1494)
@@ -16,4 +16,4 @@
import IPython.ipapi
ip=IPython.ipapi.get()
ip.set_hook('complete_command',yade_completers,re_key='.*\\b[a-zA-Z0-9_]+\[["\'][^"\'\]]*$')
-
+ip.options.profile='yade'
Modified: trunk/gui/py/utils.py
===================================================================
--- trunk/gui/py/utils.py 2008-08-22 23:01:31 UTC (rev 1493)
+++ trunk/gui/py/utils.py 2008-08-25 07:58:04 UTC (rev 1494)
@@ -4,6 +4,8 @@
#
# 2008 © Václav Šmilauer <eudoxos@xxxxxxxx>
+#from yade._utils import *
+
import math,random
from yade.wrapper import *
try: # use psyco if available
@@ -11,6 +13,10 @@
psyco.full()
except ImportError: pass
+# c++ implementations for performance reasons
+from yade._utils import *
+
+
def sphere(center,radius,density=1,young=30e9,poisson=.3,frictionAngle=0.5236,dynamic=True,wire=False,color=[1,1,1],physParamsClass='BodyMacroParameters'):
"""Create default sphere, with given parameters. Physical properties such as mass and inertia are calculated automatically."""
s=Body()
@@ -18,7 +24,7 @@
s.mold=InteractingGeometry('InteractingSphere',{'radius':radius,'diffuseColor':color})
V=(4./3)*math.pi*radius**3
inert=(2./5.)*V*density*radius**2
- s.phys=PhysicalParameters(physParamsClass,{'se3':[center[0],center[1],center[2],1,0,0,0],'mass':V*density,'inertia':[inert,inert,inert],'young':young,'poisson':poisson,'frictionAngle':frictionAngle})
+ s.phys=PhysicalParameters(physParamsClass,{'se3':[center[0],center[1],center[2],1,0,0,0],'refSe3':[center[0],center[1],center[2],1,0,0,0],'mass':V*density,'inertia':[inert,inert,inert],'young':young,'poisson':poisson,'frictionAngle':frictionAngle})
s.bound=BoundingVolume('AABB',{'diffuseColor':[0,1,0]})
s['isDynamic']=dynamic
return s
@@ -29,29 +35,11 @@
b.shape=GeometricalModel('Box',{'extents':extents,'diffuseColor':color,'wire':wire,'visible':True})
b.mold=InteractingGeometry('InteractingBox',{'extents':extents,'diffuseColor':color})
mass=8*extents[0]*extents[1]*extents[2]*density
- b.phys=PhysicalParameters(physParamsClass,{'se3':[center[0],center[1],center[2],orientation[0],orientation[1],orientation[2],orientation[3]],'mass':mass,'inertia':[mass*4*(extents[1]**2+extents[2]**2),mass*4*(extents[0]**2+extents[2]**2),mass*4*(extents[0]**2+extents[1]**2)],'young':young,'poisson':poisson,'frictionAngle':frictionAngle})
+ b.phys=PhysicalParameters(physParamsClass,{'se3':[center[0],center[1],center[2],orientation[0],orientation[1],orientation[2],orientation[3]],'refSe3':[center[0],center[1],center[2],orientation[0],orientation[1],orientation[2],orientation[3]],'mass':mass,'inertia':[mass*4*(extents[1]**2+extents[2]**2),mass*4*(extents[0]**2+extents[2]**2),mass*4*(extents[0]**2+extents[1]**2)],'young':young,'poisson':poisson,'frictionAngle':frictionAngle})
b.bound=BoundingVolume('AABB',{'diffuseColor':[0,1,0]})
b['isDynamic']=dynamic
return b
-def negPosExtremes(axis,distFactor=1.1):
- """Get 2 lists (negative and positive extremes) of sphere ids that are close to the boundary
- in the direction of requested axis (0=x,1=y,2=z).
-
- distFactor enlarges radius of the sphere, it can be considered 'extremal' even if it doesn't touch the extreme.
- """
- inf=float('inf')
- extremes=[inf,-inf]
- ret=[[],[]]
- o=Omega()
- for b in o.bodies:
- extremes[1]=max(extremes[1],+b.shape['radius']+b.phys['se3'][axis])
- extremes[0]=min(extremes[0],-b.shape['radius']+b.phys['se3'][axis])
- for b in o.bodies:
- if b.phys['se3'][axis]-b.shape['radius']*distFactor<=extremes[0]: ret[0].append(b.id)
- if b.phys['se3'][axis]+b.shape['radius']*distFactor>=extremes[1]: ret[1].append(b.id)
- return ret
-
def aabbWalls(extrema=None,thickness=None,oversizeFactor=1.5,**kw):
"""return 6 walls that will wrap existing packing;
extrema are extremal points of the AABB of the packing (will be calculated if not specified)
@@ -74,23 +62,13 @@
walls[-1].shape['visible']=True
return walls
-def aabbExtrema(consider=lambda body:True):
- """return min and max points of an aabb around spherical packing (non-spheres are ignored)"""
- inf=float('inf')
- minimum,maximum=[inf,inf,inf],[-inf,-inf,-inf]
- o=Omega()
- for b in o.bodies:
- if consider(b) and b.shape.name=='Sphere':
- minimum=[min(minimum[i],b.phys['se3'][i]-b.shape['radius']) for i in range(3)]
- maximum=[max(maximum[i],b.phys['se3'][i]+b.shape['radius']) for i in range(3)]
- return minimum,maximum
-def perpendicularArea(axis,consider=lambda body: True):
+def perpendicularArea(axis):
"""return area perpendicular to given axis (0=x,1=y,2=z) generated by bodies
for which the function consider returns True (defaults to returning True always)
and which is of the type "Sphere"
"""
- ext=aabbExtrema(consider)
+ ext=aabbExtrema()
other=((axis+1)%3,(axis+2)%3)
return (ext[1][other[0]]-ext[0][other[0]])*(ext[1][other[1]]-ext[0][other[1]])
@@ -107,19 +85,6 @@
if onShapes and (b['isDynamic'] or not onlyDynamic): b.shape['diffuseColor']=color
if onMolds and (b['isDynamic'] or not onlyDynamic): b.mold['diffuseColor']=color
-def PWaveTimeStep():
- """Estimate timestep by the elastic wave propagation speed (p-wave).
-
- Multiply the value returned by some safety factor < 1, since shear waves are not taken into account here."""
- o=Omega()
- dt=float('inf')
- for b in o.bodies:
- if not (b.phys and b.shape): continue
- if not (b.phys.has_key('young') and b.shape.has_key('radius')): continue
- density=b.phys['mass']/((4./3.)*math.pi*b.shape['radius']**3)
- thisDt=b.shape['radius']/math.sqrt(b.phys['young']/density)
- dt=min(dt,thisDt)
- return dt
def spheresFromFile(filename,scale=1.,wenjieFormat=False,**kw):
"""Load sphere coordinates from file, create spheres, insert them to the simulation.
@@ -150,3 +115,48 @@
out.write('%g\t%g\t%g\t%g\n'%(b.phys['se3'][0],b.phys['se3'][1],b.phys['se3'][2],b.shape['radius']))
out.close()
+def negPosExtremes(**kw): raise DeprecationWarning("negPosExtremes is deprecated, use negPosExtremalIds instead.")
+
+# reimplemented in _utils in c++ (results verified to be identical)
+if 0:
+ def negPosExtremes(axis,distFactor=1.1):
+ """Get 2 lists (negative and positive extremes) of sphere ids that are close to the boundary
+ in the direction of requested axis (0=x,1=y,2=z).
+
+ distFactor enlarges radius of the sphere, it can be considered 'extremal' even if it doesn't touch the extreme.
+ """
+ inf=float('inf')
+ extremes=[inf,-inf]
+ ret=[[],[]]
+ o=Omega()
+ for b in o.bodies:
+ extremes[1]=max(extremes[1],+b.shape['radius']+b.phys['se3'][axis])
+ extremes[0]=min(extremes[0],-b.shape['radius']+b.phys['se3'][axis])
+ for b in o.bodies:
+ if b.phys['se3'][axis]-b.shape['radius']*distFactor<=extremes[0]: ret[0].append(b.id)
+ if b.phys['se3'][axis]+b.shape['radius']*distFactor>=extremes[1]: ret[1].append(b.id)
+ return ret
+ def aabbExtrema(consider=lambda body:True):
+ """return min and max points of an aabb around spherical packing (non-spheres are ignored)"""
+ inf=float('inf')
+ minimum,maximum=[inf,inf,inf],[-inf,-inf,-inf]
+ o=Omega()
+ for b in o.bodies:
+ if consider(b) and b.shape.name=='Sphere':
+ minimum=[min(minimum[i],b.phys['se3'][i]-b.shape['radius']) for i in range(3)]
+ maximum=[max(maximum[i],b.phys['se3'][i]+b.shape['radius']) for i in range(3)]
+ return minimum,maximum
+ def PWaveTimeStep():
+ """Estimate timestep by the elastic wave propagation speed (p-wave).
+
+ Multiply the value returned by some safety factor < 1, since shear waves are not taken into account here."""
+ o=Omega()
+ dt=float('inf')
+ for b in o.bodies:
+ if not (b.phys and b.shape): continue
+ if not (b.phys.has_key('young') and b.shape.has_key('radius')): continue
+ density=b.phys['mass']/((4./3.)*math.pi*b.shape['radius']**3)
+ thisDt=b.shape['radius']/math.sqrt(b.phys['young']/density)
+ dt=min(dt,thisDt)
+ return dt
+
