------------------------------------------------------------ revno: 2525 committer: Václav Šmilauer <eudo...@arcig.cz> 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&>())
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp