yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06050
[Branch ~yade-dev/yade/trunk] Rev 2525: 1. Fix intiialization problem in Shop::kineticEnergy (thanks, Anton)
------------------------------------------------------------
revno: 2525
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Tue 2010-11-02 13:02:13 +0100
message:
1. Fix intiialization problem in Shop::kineticEnergy (thanks, Anton)
2. Several enhancements of the post2d module
3. Fix a few crashers with deleted particles.
4. Improve utils.facetCylinder
modified:
doc/sphinx/upload
gui/qt4/GLViewer.cpp
gui/qt4/__init__.py
lib/smoothing/WeightedAverage2d.hpp
pkg/common/OpenGLRenderer.cpp
pkg/dem/ElasticContactLaw.hpp
pkg/dem/Shop.cpp
py/WeightedAverage2d.cpp
py/post2d.py
py/utils.py
py/wrapper/yadeWrapper.cpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'doc/sphinx/upload'
--- doc/sphinx/upload 2010-09-10 11:33:21 +0000
+++ doc/sphinx/upload 2010-11-02 12:02:13 +0000
@@ -1,2 +1,2 @@
#!/bin/sh
-rsync -rv _build/html/ gamma:www/NoBackup/yade-sphinx
+rsync -rv _build/html/ ggamma:www/NoBackup/yade-sphinx
=== modified file 'gui/qt4/GLViewer.cpp'
--- gui/qt4/GLViewer.cpp 2010-09-27 17:47:59 +0000
+++ gui/qt4/GLViewer.cpp 2010-11-02 12:02:13 +0000
@@ -426,6 +426,7 @@
Real inf=std::numeric_limits<Real>::infinity();
min=Vector3r(inf,inf,inf); max=Vector3r(-inf,-inf,-inf);
FOREACH(const shared_ptr<Body>& b, *rb->bodies){
+ if(!b) continue;
max=max.cwise().max(b->state->pos);
min=min.cwise().min(b->state->pos);
}
=== modified file 'gui/qt4/__init__.py'
--- gui/qt4/__init__.py 2010-10-29 10:12:44 +0000
+++ gui/qt4/__init__.py 2010-11-02 12:02:13 +0000
@@ -253,7 +253,7 @@
subStepInfo="<br><small>sub %d/%d [%s]</small>"%(subStep,len(O.engines),subStepInfo)
self.subStepCheckbox.setChecked(O.subStepping) # might have been changed async
if stopAtIter<=iter:
- self.realTimeLabel.setText('%02d:%02d:%02d'%(rt//3600,rt//60,rt%60))
+ self.realTimeLabel.setText('%02d:%02d:%02d'%(rt//3600,(rt%3600)//60,rt%60))
self.iterLabel.setText('#%ld, %.1f/s %s'%(iter,self.iterPerSec,subStepInfo))
else:
e=int((stopAtIter-iter)*self.iterPerSec)
=== modified file 'lib/smoothing/WeightedAverage2d.hpp'
--- lib/smoothing/WeightedAverage2d.hpp 2010-08-24 12:54:14 +0000
+++ lib/smoothing/WeightedAverage2d.hpp 2010-11-02 12:02:13 +0000
@@ -14,6 +14,8 @@
#include<boost/python.hpp>
#include<yade/extra/boost_python_len.hpp>
+#include<boost/math/distributions/normal.hpp>
+
#ifndef FOREACH
#define FOREACH BOOST_FOREACH
#endif
@@ -32,6 +34,7 @@
public:
Vector2r getLo(){return lo;}
Vector2r getHi(){return hi;}
+ Vector2r getCellSize(){ return cellSizes; }
vector<vector<vector<T> > > grid;
/* construct grid from lower left corner, upper right corner and number of cells in each direction */
GridContainer(Vector2r _lo, Vector2r _hi, Vector2i _nCells): lo(_lo), hi(_hi), nCells(_nCells){
@@ -95,21 +98,29 @@
template<typename T, typename Tvalue>
struct WeightedAverage{
const shared_ptr<GridContainer<T> > grid;
+ Real weightedSupportArea; // must be computed by derived class
WeightedAverage(const shared_ptr<GridContainer<T> >& _grid):grid(_grid){};
virtual Vector2r getPosition(const T&)=0;
virtual Real getWeight(const Vector2r& refPt, const T&)=0;
virtual Tvalue getValue(const T&)=0;
virtual vector<Vector2i> filterCells(const Vector2r& refPt)=0;
Tvalue computeAverage(const Vector2r& refPt){
+ Real sumValues, sumWeights; sumValuesWeights(refPt,sumValues,sumWeights);
+ return sumValues/sumWeights;
+ }
+ Tvalue computeAvgPerUnitArea(const Vector2r& refPt){
+ Real sumValues, sumWeights; sumValuesWeights(refPt,sumValues,sumWeights);
+ return sumValues/weightedSupportArea;
+ }
+ void sumValuesWeights(const Vector2r& refPt, Real& sumValues, Real& sumWeights){
vector<Vector2i> filtered=filterCells(refPt);
- Real sumValues=0, sumWeights=0;
+ sumValues=sumWeights=0;
FOREACH(Vector2i cell, filtered){
FOREACH(const T& element, grid->grid[cell[0]][cell[1]]){
Real weight=getWeight(refPt,element);
sumValues+=weight*getValue(element); sumWeights+=weight;
}
}
- return sumValues/sumWeights;
}
};
@@ -122,13 +133,25 @@
/* Symmetric Gaussian Distribution Average with scalar values
*/
struct SGDA_Scalar2d: public WeightedAverage<Scalar2d,Real> {
- Real stDev, relThreshold;
- SGDA_Scalar2d(const shared_ptr<GridContainer<Scalar2d> >& _grid, Real _stDev, Real _relThreshold=3): WeightedAverage<Scalar2d,Real>(_grid), stDev(_stDev), relThreshold(_relThreshold){}
+ Real stDev, relThreshold; boost::math::normal_distribution<Real> distrib;
+ SGDA_Scalar2d(const shared_ptr<GridContainer<Scalar2d> >& _grid, Real _stDev, Real _relThreshold=3): WeightedAverage<Scalar2d,Real>(_grid), stDev(_stDev), relThreshold(_relThreshold), distrib(0,stDev) {
+ // approximate probability density function between -stDev*relThreshold,stDev*relThreshold
+ // it is enough to get Φ(-stDev*relThreshold) and subtract twice from 1 (symmetry)
+ // http://en.wikipedia.org/wiki/Normal_distribution#Numerical_approximations_for_the_normal_CDF
+ // using Abramowitz & Stegun approximation, which has error less that 7.5e-8 (fine for us)
+
+ // FIXME: algorithm not correct, as it takes 1d quantile, while we would need PDF for symmetric 2d gaussian!
+ Real clippedQuantile=boost::math::cdf(distrib,-stDev*relThreshold);
+ Real area=M_PI*pow(relThreshold*stDev,2); // area of the support
+ weightedSupportArea=(1-2*clippedQuantile)*area;
+
+ }
virtual Real getWeight(const Vector2r& meanPt, const Scalar2d& e){
Vector2r pos=getPosition(e);
- Real rSq=pow(meanPt[0]-pos[0],2)+pow(meanPt[1]-pos[1],2);
+ Real rSq=(meanPt-pos).squaredNorm(); //pow(meanPt[0]-pos[0],2)+pow(meanPt[1]-pos[1],2);
if(rSq>pow(relThreshold*stDev,2)) return 0.; // discard points further than relThreshold*stDev, by default 3*stDev
- return (1./(stDev*sqrt(2*M_PI)))*exp(-rSq/(2*stDev*stDev));
+ //return (1./(stDev*sqrt(2*M_PI)))*exp(-rSq/(2*stDev*stDev));
+ return boost::math::pdf(distrib,sqrt(rSq));
}
vector<Vector2i> filterCells(const Vector2r& refPt){return WeightedAverage<Scalar2d,Real>::grid->circleFilter(refPt,stDev*relThreshold);}
Real getValue(const Scalar2d& dp){return (Real)dp.val;}
@@ -166,7 +189,8 @@
return false;
}
bool addPt(Real val, python::tuple pos){Scalar2d d; d.pos=tuple2vec2r(pos); if(ptIsClipped(d.pos)) return false; d.val=val; sgda->grid->add(d,d.pos); return true; }
- Real avg(python::tuple _pt){Vector2r pt=tuple2vec2r(_pt); if(ptIsClipped(pt)) return std::numeric_limits<Real>::quiet_NaN(); return sgda->computeAverage(pt);}
+ Real avg(Vector2r pt){ if(ptIsClipped(pt)) return std::numeric_limits<Real>::quiet_NaN(); return sgda->computeAverage(pt);}
+ Real avgPerUnitArea(Vector2r pt){ if(ptIsClipped(pt)) return std::numeric_limits<Real>::quiet_NaN(); return sgda->computeAvgPerUnitArea(pt); }
Real stDev_get(){return sgda->stDev;} void stDev_set(Real s){sgda->stDev=s;}
Real relThreshold_get(){return sgda->relThreshold;} void relThreshold_set(Real rt){sgda->relThreshold=rt;}
python::tuple aabb_get(){return python::make_tuple(sgda->grid->getLo(),sgda->grid->getHi());}
@@ -204,5 +228,11 @@
}
return python::make_tuple(x,y,val);
}
+ Vector2i nCells_get(){ return sgda->grid->getSize(); }
+ int cellNum(const Vector2i& cell){ return sgda->grid->grid[cell[0]][cell[1]].size(); }
+ Real cellSum(const Vector2i& cell){ Real sum=0; FOREACH(const Scalar2d& v, sgda->grid->grid[cell[0]][cell[1]]) sum+=v.val; return sum; }
+ Real cellAvg(const Vector2i& cell){ return cellSum(cell)/cellNum(cell); }
+ Real cellArea(){ Vector2r sz=sgda->grid->getCellSize(); return sz[0]*sz[1]; }
+ Vector2r cellDim(){ return sgda->grid->getCellSize(); }
};
=== modified file 'pkg/common/OpenGLRenderer.cpp'
--- pkg/common/OpenGLRenderer.cpp 2010-10-13 16:23:08 +0000
+++ pkg/common/OpenGLRenderer.cpp 2010-11-02 12:02:13 +0000
@@ -209,6 +209,7 @@
void OpenGLRenderer::renderAllInteractionsWire(){
FOREACH(const shared_ptr<Interaction>& i, *scene->interactions){
+ if(!i->functorCache.geomExists) continue;
glColor3v(i->isReal()? Vector3r(0,1,0) : Vector3r(.5,0,1));
Vector3r p1=Body::byId(i->getId1(),scene)->state->pos;
const Vector3r& size=scene->cell->getSize();
=== modified file 'pkg/dem/ElasticContactLaw.hpp'
--- pkg/dem/ElasticContactLaw.hpp 2010-10-29 10:12:44 +0000
+++ pkg/dem/ElasticContactLaw.hpp 2010-11-02 12:02:13 +0000
@@ -21,7 +21,7 @@
//Real elasticEnergy ();
//Real getPlasticDissipation();
//void initPlasticDissipation(Real initVal=0);
- YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_FrictPhys_CundallStrack,LawFunctor,"Law for linear compression, and Mohr-Coulomb plasticity surface without cohesion.\nThis law implements the classical linear elastic-plastic law from [CundallStrack1979]_ (see also [Pfc3dManual30]_). The normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, 0)$. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force : $F_s^{max}=F_n*tan(\\phi)$, with $\\phi$ the friction angle.\n\n.. note::\n This law uses :yref:`ScGeom`; there is also functionally equivalent :yref:`Law2_Dem3DofGeom_FrictPhys_CundallStrack`, which uses :yref:`Dem3DofGeom` (sphere-box interactions are not implemented for the latest).\n\n.. note::\n This law is well tested in the context of triaxial simulation, and has been used for a number of published results (see e.g. [Scholtes2009b]_ and other papers from the same authors). It is generalised by :yref:`Law2_ScGeom_CohFrictPhys_CohesionMoment`, which adds cohesion and moments at contact.",
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_FrictPhys_CundallStrack,LawFunctor,"Law for linear compression, and Mohr-Coulomb plasticity surface without cohesion.\nThis law implements the classical linear elastic-plastic law from [CundallStrack1979]_ (see also [Pfc3dManual30]_). The normal force is (with the convention of positive tensile forces) $F_n=\\min(k_n u_n, 0)$. The shear force is $F_s=k_s u_s$, the plasticity condition defines the maximum value of the shear force : $F_s^{\\max}=F_n\\tan(\\phi)$, with $\\phi$ the friction angle.\n\n.. note::\n This law uses :yref:`ScGeom`; there is also functionally equivalent :yref:`Law2_Dem3DofGeom_FrictPhys_CundallStrack`, which uses :yref:`Dem3DofGeom` (sphere-box interactions are not implemented for the latest).\n\n.. note::\n This law is well tested in the context of triaxial simulation, and has been used for a number of published results (see e.g. [Scholtes2009b]_ and other papers from the same authors). It is generalised by :yref:`Law2_ScGeom_CohFrictPhys_CohesionMoment`, which adds cohesion and moments at contact.",
((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
((bool,traceEnergy,false,Attr::hidden,"Define the total energy dissipated in plastic slips at all contacts."))
((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic dissipation (with O.trackEnergy)"))
=== modified file 'pkg/dem/Shop.cpp'
--- pkg/dem/Shop.cpp 2010-10-26 13:41:30 +0000
+++ pkg/dem/Shop.cpp 2010-11-02 12:02:13 +0000
@@ -201,7 +201,7 @@
Matrix3r T(state->ori);
// the tensorial expression http://en.wikipedia.org/wiki/Moment_of_inertia#Moment_of_inertia_tensor
// inertia tensor rotation from http://www.kwon3d.com/theory/moi/triten.html
- Matrix3r mI; mI.diagonal()=state->inertia;
+ Matrix3r mI; mI<<state->inertia[0],0,0, 0,state->inertia[1],0, 0,0,state->inertia[2];
E+=.5*state->angVel.transpose().dot((T.transpose()*mI*T)*state->angVel);
}
else { E+=state->angVel.dot(state->inertia.cwise()*state->angVel);}
@@ -465,7 +465,7 @@
shared_ptr<Scene> rb=(_rb?_rb:Omega::instance().getScene());
Real dt=std::numeric_limits<Real>::infinity();
FOREACH(const shared_ptr<Body>& b, *rb->bodies){
- if(!b->material || !b->shape) continue;
+ if(!b || !b->material || !b->shape) continue;
shared_ptr<ElastMat> ebp=dynamic_pointer_cast<ElastMat>(b->material);
shared_ptr<Sphere> s=dynamic_pointer_cast<Sphere>(b->shape);
if(!ebp || !s) continue;
=== modified file 'py/WeightedAverage2d.cpp'
--- py/WeightedAverage2d.cpp 2010-01-04 10:34:53 +0000
+++ py/WeightedAverage2d.cpp 2010-11-02 12:02:13 +0000
@@ -20,11 +20,18 @@
boost::python::class_<pyGaussAverage>("GaussAverage",python::init<python::tuple,python::tuple,python::tuple,Real,python::optional<Real> >(python::args("min","max","nCells","stDev","relThreshold"),"Create empty container for data, which can be added using add and later retrieved using avg."))
.def("add",&pyGaussAverage::addPt)
.def("avg",&pyGaussAverage::avg)
+ .def("avgPerUnitArea",&pyGaussAverage::avgPerUnitArea)
+ .def("cellNum",&pyGaussAverage::cellNum)
+ .def("cellSum",&pyGaussAverage::cellSum)
+ .def("cellAvg",&pyGaussAverage::cellAvg)
.add_property("stDev",&pyGaussAverage::stDev_get,&pyGaussAverage::stDev_set)
.add_property("relThreshold",&pyGaussAverage::relThreshold_get,&pyGaussAverage::relThreshold_set)
.add_property("clips",&pyGaussAverage::clips_get,&pyGaussAverage::clips_set)
.add_property("data",&pyGaussAverage::data_get)
.add_property("aabb",&pyGaussAverage::aabb_get)
+ .add_property("nCells",&pyGaussAverage::nCells_get)
+ .add_property("cellArea",&pyGaussAverage::cellArea)
+ .add_property("cellDim",&pyGaussAverage::cellDim)
;
};
=== modified file 'py/post2d.py'
--- py/post2d.py 2010-10-30 17:08:51 +0000
+++ py/post2d.py 2010-11-02 12:02:13 +0000
@@ -59,6 +59,7 @@
"""
from yade.wrapper import *
+from miniEigen import *
class Flatten:
"""Abstract class for converting 3d point into 2d. Used by post2d.data2d."""
@@ -149,7 +150,7 @@
def normal(self,pos,vec):
return vec[self.axis]
-def data(extractor,flattener,intr=False,onlyDynamic=True,stDev=None,relThreshold=3.,div=(50,50),margin=(0,0),radius=1):
+def data(extractor,flattener,intr=False,onlyDynamic=True,stDev=None,relThreshold=3.,perArea=0,div=(50,50),margin=(0,0),radius=1):
"""Filter all bodies/interactions, project them to 2d and extract required scalar value;
return either discrete array of positions and values, or smoothed data, depending on whether the stDev
value is specified.
@@ -163,6 +164,7 @@
:param bool onlyDynamic: skip all non-dynamic bodies
:param float/None stDev: standard deviation for averaging, enables smoothing; ``None`` (default) means raw mode, where discrete points are returned
:param float relThreshold: threshold for the gaussian weight function relative to stDev (smooth mode only)
+ :param int perArea: if 1, compute weightedSum/weightedArea rather than weighted average (weightedSum/sumWeights); the first is useful to compute average stress; if 2, compute averages on subdivision elements, not using weight function
:param (int,int) div: number of cells for the gaussian grid (smooth mode only)
:param (float,float) margin: x,y margins around bounding box for data (smooth mode only)
:param float/callable radius: Fallback value for radius (for raw plotting) for non-spherical bodies or interactions; if a callable, receives body/interaction and returns radius
@@ -214,14 +216,29 @@
xxx,yyy=[numpy.arange(llo[i]+.5*step[i],hhi[i],step[i]) for i in [0,1]]
ddd=numpy.zeros((len(yyy),len(xxx)),float)
ddd2=numpy.zeros((len(yyy),len(xxx)),float)
+ # set the type of average we are going to use
+ if perArea==0:
+ def compAvg(gauss,coord): return float(gauss.avg(coord))
+ elif perArea==1:
+ def compAvg(gauss,coord,cellCoord): return gauss.avgPerUnitArea(coord)
+ elif perArea==2:
+ def compAvg(gauss,coord,cellCoord):
+ s=gauss.cellSum(cellCoord);
+ return (s/gauss.cellArea) if s>0 else float('nan')
+ elif perArea==3:
+ def compAvg(gauss,coord,cellCoord):
+ s=gauss.cellSum(cellCoord);
+ return s if s>0 else float('nan')
+ else: raise RuntimeError('Invalid value of *perArea*, must be one of 0,1,2,3.')
+ #
for cx in range(0,div[0]):
for cy in range(0,div[1]):
- ddd[cy,cx]=float(ga.avg((xxx[cx],yyy[cy])))
- if nDim>1: ddd2[cy,cx]=float(ga2.avg((xxx[cx],yyy[cy])))
- if nDim==1: return {'type':'smoothScalar','x':xxx,'y':yyy,'val':ddd,'bbox':(llo,hhi)}
- else: return {'type':'smoothVector','x':xxx,'y':yyy,'valX':ddd,'valY':ddd2,'bbox':(llo,hhi)}
+ ddd[cy,cx]=compAvg(ga,(xxx[cx],yyy[cy]),(cx,cy))
+ if nDim>1: ddd2[cy,cx]=compAvg(ga2,(xxx[cx],yyy[cy]),(cx,cy))
+ if nDim==1: return {'type':'smoothScalar','x':xxx,'y':yyy,'val':ddd,'bbox':(llo,hhi),'perArea':perArea,'grid':ga}
+ else: return {'type':'smoothVector','x':xxx,'y':yyy,'valX':ddd,'valY':ddd2,'bbox':(llo,hhi),'grid':ga,'grid2':ga2}
-def plot(data,axes=None,alpha=.5,clabel=True,cbar=False,**kw):
+def plot(data,axes=None,alpha=.5,clabel=True,cbar=False,aspect='equal',**kw):
"""Given output from post2d.data, plot the scalar as discrete or smooth plot.
For raw discrete data, plot filled circles with radii of particles, colored by the scalar value.
@@ -232,7 +249,7 @@
:param axes: matplotlib.axes\ instance where the figure will be plotted; if None, will be created from scratch.
:param data: value returned by :yref:`yade.post2d.data`
- :param bool clabel: show contour labels (smooth mode only)
+ :param bool clabel: show contour labels (smooth mode only), or annotate cells with numbers inside (with perArea==2)
:param bool cbar: show colorbar (equivalent to calling pylab.colorbar(mappable) on the returned mappable)
:return: tuple of ``(axes,mappable)``; mappable can be used in further calls to pylab.colorbar.
@@ -251,17 +268,23 @@
axes.add_collection(coll)
axes.set_xlim(bb.xmin,bb.xmax); axes.set_ylim(bb.ymin,bb.ymax)
if cbar: axes.get_figure().colorbar(coll)
- axes.grid(True); axes.set_aspect('equal')
+ axes.grid(True); axes.set_aspect(aspect)
return axes,coll
elif data['type']=='smoothScalar':
loHi=data['bbox']
- img=axes.imshow(data['val'],extent=(loHi[0][0],loHi[1][0],loHi[0][1],loHi[1][1]),origin='lower',aspect='equal',**kw)
- ct=axes.contour(data['x'],data['y'],data['val'],colors='k',origin='lower',extend='both')
+ if data['perArea'] in (0,1):
+ img=axes.imshow(data['val'],extent=(loHi[0][0],loHi[1][0],loHi[0][1],loHi[1][1]),origin='lower',aspect=aspect,**kw)
+ ct=axes.contour(data['x'],data['y'],data['val'],colors='k',origin='lower',extend='both')
+ if clabel: axes.clabel(ct,inline=1,fontsize=10)
+ else:
+ img=axes.imshow(data['val'],extent=(loHi[0][0],loHi[1][0],loHi[0][1],loHi[1][1]),origin='lower',aspect=aspect,interpolation='nearest',**kw)
+ xStep=(data['x'][1]-data['x'][0]) if len(data['x'])>1 else 0
+ for y,valLine in zip(data['y'],data['val']):
+ for x,val in zip(data['x'],valLine): axes.text(x-.4*xStep,y,('-' if math.isnan(val) else '%5g'%val),size=4)
axes.update_datalim(loHi)
- if clabel: axes.clabel(ct,inline=1,fontsize=10)
axes.set_xlim(loHi[0][0],loHi[1][0]); axes.set_ylim(loHi[0][1],loHi[1][1])
if cbar: axes.get_figure().colorbar(img)
- axes.grid(True); axes.set_aspect('equal')
+ axes.grid(True if data['perArea'] in (0,1) else False); axes.set_aspect(aspect)
return axes,img
elif data['type'] in ('rawVector','smoothVector'):
import numpy
@@ -273,7 +296,7 @@
#axes.update_datalim(loHi)
axes.set_xlim(loHi[0][0],loHi[1][0]); axes.set_ylim(loHi[0][1],loHi[1][1])
if cbar: axes.get_figure().colorbar(coll)
- axes.grid(True); axes.set_aspect('equal')
+ axes.grid(True); axes.set_aspect(aspect)
return axes,quiv
=== modified file 'py/utils.py'
--- py/utils.py 2010-10-30 17:08:51 +0000
+++ py/utils.py 2010-11-02 12:02:13 +0000
@@ -345,7 +345,7 @@
if wallMask&32: ret+=doWall(E,H,G,F)
return ret
-def facetCylinder(center,radius,height,orientation=Quaternion.Identity,segmentsNumber=10,wallMask=7,angleRange=2.0*math.pi,**kw):
+def facetCylinder(center,radius,height,orientation=Quaternion.Identity,segmentsNumber=10,wallMask=7,angleRange=None,closeGap=False,**kw):
"""
Create arbitrarily-aligned cylinder composed of facets, with given center, radius, height and orientation.
Return List of facets forming the cylinder;
@@ -356,7 +356,8 @@
:param Quaternion orientation: orientation of the cylinder; the reference orientation has axis along the $+x$ axis.
:param int segmentsNumber: number of edges on the cylinder surface (>=5)
:param bitmask wallMask: determines which walls will be created, in the order up (1), down (2), side (4). The numbers are ANDed; the default 7 means to create all walls
- :param float angleRange: allows to create only part of cylinder, 2.0*math.pi means the whole cylinder, 1.0*math.pi - the half etc;
+ :param (θmin,Θmax) angleRange: allows to create only part of cylinder by specifying range of angles; if ``None``, (0,2*pi) is assumed.
+ :param bool closeGap: close range skipped in angleRange with triangular facets at cylinder bases.
:param **kw: (unused keyword arguments) passed to utils.facet;
"""
@@ -366,7 +367,11 @@
if (radius<=0): raise RuntimeError("The radius should have the positive value");
import numpy
- anglesInRad = numpy.linspace(0, angleRange, segmentsNumber+1, endpoint=True)
+ if angleRange==None: angleRange=(0,2*math.pi)
+ if isinstance(angleRange,float):
+ print u'WARNING: utils.facetCylinder,angleRange should be (Θmin,Θmax), not just Θmax (one number), update your code.'
+ angleRange=(0,angleRange)
+ anglesInRad = numpy.linspace(angleRange[0], angleRange[1], segmentsNumber+1, endpoint=True)
P1=[]; P2=[]
P1.append(Vector3(0,0,-height/2))
P2.append(Vector3(0,0,+height/2))
@@ -389,6 +394,14 @@
if wallMask&4:
ret.append(facet((P1[i],P2[i],P2[i-1]),**kw))
ret.append(facet((P2[i-1],P1[i-1],P1[i]),**kw))
+ if closeGap and (angleRange[0]%(2*math.pi))!=(angleRange[1]%(2*math.pi)): # some part is skipped
+ pts=[(radius*math.cos(angleRange[i]),radius*math.sin(angleRange[i])) for i in (0,1)]
+ for Z in -height/2,height/2:
+ #print (pts[0][0],pts[0][1],Z),(pts[1][0],pts[1][1],Z),(0,0,Z)
+ pp=[(pts[0][0],pts[0][1],Z),(pts[1][0],pts[1][1],Z),(0,0,Z)]
+ pp=[orientation*p+center for p in pp]
+ ret.append(facet(pp,**kw))
+
return ret
=== modified file 'py/wrapper/yadeWrapper.cpp'
--- py/wrapper/yadeWrapper.cpp 2010-10-29 10:12:44 +0000
+++ py/wrapper/yadeWrapper.cpp 2010-11-02 12:02:13 +0000
@@ -213,6 +213,7 @@
bool serializeSorted_get(){return proxee->serializeSorted;}
void serializeSorted_set(bool ss){proxee->serializeSorted=ss;}
void eraseNonReal(){ proxee->eraseNonReal(); }
+ void erase(Body::id_t id1, Body::id_t id2){ proxee->requestErase(id1,id2); }
};
class pyForceContainer{
@@ -589,6 +590,7 @@
.def("withBody",&pyInteractionContainer::withBody,"Return list of real interactions of given body.")
.def("withBodyAll",&pyInteractionContainer::withBodyAll,"Return list of all (real as well as non-real) interactions of given body.")
.def("eraseNonReal",&pyInteractionContainer::eraseNonReal,"Erase all interactions that are not :yref:`real <InteractionContainer.isReal>`.")
+ .def("erase",&pyInteractionContainer::erase,"Erase one interaction, given by id1, id2 (internally, ``requestErase`` is called -- the interaction might still exist as potential, if the :yref:`Collider` decides so).")
.add_property("serializeSorted",&pyInteractionContainer::serializeSorted_get,&pyInteractionContainer::serializeSorted_set)
.def("clear",&pyInteractionContainer::clear,"Remove all interactions");
python::class_<pyInteractionIterator>("InteractionIterator",python::init<pyInteractionIterator&>())
