------------------------------------------------------------ revno: 1902 committer: Václav Šmilauer <[email protected]> branch nick: trunk timestamp: Tue 2009-12-15 10:41:07 +0100 message: 1. Remove Scene::{cellMin,cellMax}, use just Scene::cellSize; update all related code and scripts 2. Change python interface: Omega.periodicCell=Vector3 (only size, no more min/max) 3. Add (not yet used) core/Cell.hpp header. Will be accessible as Scene::cell, once used. added: core/Cell.hpp modified: core/Interaction.hpp core/Scene.cpp core/Scene.hpp examples/concrete/periodic.py examples/concrete/uniax.py gui/qt3/GLViewer.cpp pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp pkg/common/Engine/GlobalEngine/InsertionSortCollider.hpp pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp pkg/common/RenderingEngine/OpenGLRenderingEngine.hpp pkg/dem/DataClass/SpherePack.cpp pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp py/pack.py py/yadeWrapper/yadeWrapper.cpp scripts/test/pack-cloud.py scripts/test/periodic-compress.py scripts/test/periodic-grow.py scripts/test/periodic-simple.py scripts/test/periodic-triax.py
-- 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.
=== added file 'core/Cell.hpp' --- core/Cell.hpp 1970-01-01 00:00:00 +0000 +++ core/Cell.hpp 2009-12-15 09:41:07 +0000 @@ -0,0 +1,52 @@ +// 2009 © Václav Å milauer <[email protected]> +#include<yade/lib-base/yadeWm3Extra.hpp> + +/*! Periodic cell parameters and routines. Usually instantiated as Scene::cell. + +The cell has size and skew, which are independent. + +*/ +class Cell{ + public: + Vector3r size; + Vector3r skew; + Vector3r refSize; + + // caches + Vector3r _skewAngle; + Vector3r _skewSine; + Matrix3r _skewTrsf; + Matrix3r _unskewTrsf; + // should be called before every step + void updateCache(){ + for(int i=0; i<3; i++) { + _skewAngle[i]=atan(skew[i]); + _skewSine[i]=sin(_skewAngle); + } + _skewTrsf=Matrix3r(1,skew[1],skew[2],skew[0],1,skew[2],skew[0],skew[1],1); + _unskewTrsf=_skewTrst.Inverse(); + } + + /*! Apply inverse skew on point; to put it inside (unskewed) periodic cell, apply wrapPt on the returned value. */ + Vector3r unskewPt(Vector3r pt){ return _unskewTrsf*pt; } + /*! Wrap number to interval 0â¦sz */ + Real wrapNum(const Real& x, const Real& sz){ + Real norm=x/mx; return (norm-floor(norm))*sz; + } + /*! Wrap number to interval 0â¦sz; store how many intervals were wrapped in period */ + Real wrapNum(const Real& x, const Real& sz, int& period){ + Real norm=x/mx; period=(int)floor(norm); return (norm-period[i])*sz[i]; + } + /*! Wrap point to inside the periodic cell; don't compute number of periods wrapped */ + Vector3r wrapPt(const Vector3r pt){ + Vector3r ret; for(int i=0;i<3;i++) ret[i]=wrapNum(pt[i],size[i]); return ret; + } + /*! Wrap point to inside the periodic cell; period will contain by how many cells it was wrapped. */ + Vector3r wrapPt(const Vector3r pt, Vector3<int>& period){ + Vector3r ret; for(int i=0; i<3; i++){ ret[i]=wrapNum(pt[i],size[i],period[i]); } return ret; + } + + REGISTER_ATTRIBUTES(Serializable,(size)(skew)); + REGISTER_CLASS_AND_BASE(Cell,Serializable); +}; +REGISTER_SERIALIZABLE(Cell); === modified file 'core/Interaction.hpp' --- core/Interaction.hpp 2009-12-04 23:07:34 +0000 +++ core/Interaction.hpp 2009-12-15 09:41:07 +0000 @@ -38,7 +38,7 @@ //! NOTE : TriangulationCollider needs this (nothing else) bool isNeighbor; - /*! relative distance between bodies, given in (Scene::cellMax-Scene::cellMin) units + /*! relative distance between bodies, given in (Scene::cellSize) units Position of id1 must be incremented by that distance so that there is spatial interaction NOTE (tricky): cellDist must survive Interaction::reset(), it is only initialized in ctor === modified file 'core/Scene.cpp' --- core/Scene.cpp 2009-12-11 18:10:19 +0000 +++ core/Scene.cpp 2009-12-15 09:41:07 +0000 @@ -44,7 +44,7 @@ dt=1e-8; selectedBody=-1; isPeriodic=false; - cellMin=cellMax=Vector3r::ZERO; + cellSize=Vector3r::ZERO; // FIXME: move SceneShape to core and create it here right away // shape=shared_ptr<Shape>(new SceneShape); === modified file 'core/Scene.hpp' --- core/Scene.hpp 2009-12-11 18:10:19 +0000 +++ core/Scene.hpp 2009-12-15 09:41:07 +0000 @@ -65,7 +65,7 @@ long stopAtIteration; Real stopAtVirtTime; Real stopAtRealTime; - Vector3r cellMin, cellMax; + Vector3r cellSize; bool isPeriodic; bool needsInitializers; @@ -87,8 +87,7 @@ (currentIteration) (simulationTime) (stopAtIteration) - (cellMin) - (cellMax) + (cellSize) (isPeriodic) ); REGISTER_CLASS_AND_BASE(Scene,Serializable); === modified file 'examples/concrete/periodic.py' --- examples/concrete/periodic.py 2009-12-11 18:10:19 +0000 +++ examples/concrete/periodic.py 2009-12-15 09:41:07 +0000 @@ -80,7 +80,7 @@ import numpy avgRadius=numpy.average([r for c,r in sp]) O.bodies.append([utils.sphere(c,r,color=utils.randomColor()) for c,r in sp]) -O.periodicCell=(Vector3(0,0,0),sp.cellSize) +O.periodicCell=sp.cellSize axis=2 ax1=(axis+1)%3 ax2=(axis+2)%3 @@ -162,7 +162,6 @@ # FIXME: only temporary, should be .5 minMaxRatio=0.5 if mode=='tension' else 0.5 if extremum==0: return - cellSize=O.periodicCell[1]-O.periodicCell[0] print O.tags['id'],mode,strainer['strain'][axis],sigma[-1] #print 'strain',strainer['strain'],'stress',strainer['stress'] import sys; sys.stdout.flush() === modified file 'examples/concrete/uniax.py' --- examples/concrete/uniax.py 2009-12-11 18:10:19 +0000 +++ examples/concrete/uniax.py 2009-12-15 09:41:07 +0000 @@ -65,7 +65,7 @@ # z-oriented hyperboloid, length 20cm, diameter 10cm, skirt 8cm # using spheres 7mm of diameter concreteId=O.materials.append(CpmMat(young=young,frictionAngle=frictionAngle,poisson=poisson,density=4800)) -spheres=pack.randomDensePack(pack.inHyperboloid((0,0,-.5*specimenLength),(0,0,.5*specimenLength),.25*specimenLength,.2*specimenLength),spheresInCell=2000,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',material=concreteId) +spheres=pack.randomDensePack(pack.inHyperboloid((0,0,-.5*specimenLength),(0,0,.5*specimenLength),.25*specimenLength,.17*specimenLength),spheresInCell=2000,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',material=concreteId) O.bodies.append(spheres) bb=utils.uniaxialTestFeatures() negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area'] === modified file 'gui/qt3/GLViewer.cpp' --- gui/qt3/GLViewer.cpp 2009-12-04 23:46:45 +0000 +++ gui/qt3/GLViewer.cpp 2009-12-15 09:41:07 +0000 @@ -331,8 +331,8 @@ void GLViewer::centerPeriodic(){ Scene* rb=Omega::instance().getScene().get(); assert(rb->isPeriodic); - Vector3r center=.5*(rb->cellMin+rb->cellMax); - Vector3r halfSize=.5*(rb->cellMax-rb->cellMin); + Vector3r center=.5*rb->cellSize; + Vector3r halfSize=.5*rb->cellSize; float radius=std::max(halfSize[0],std::max(halfSize[1],halfSize[2])); LOG_DEBUG("Periodic scene center="<<center<<", halfSize="<<halfSize<<", radius="<<radius); setSceneCenter(qglviewer::Vec(center[0],center[1],center[2])); === modified file 'pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp' --- pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2009-12-11 07:27:29 +0000 +++ pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2009-12-15 09:41:07 +0000 @@ -37,7 +37,7 @@ geomDispatcher->updateScenePtr(); physDispatcher->updateScenePtr(); lawDispatcher->updateScenePtr(); - Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cellMax-scene->cellMin; + Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cellSize; bool removeUnseenIntrs=(scene->interactions->iterColliderLastRun>=0 && scene->interactions->iterColliderLastRun==scene->currentIteration); #ifdef YADE_OPENMP const long size=scene->interactions->size(); === modified file 'pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp' --- pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp 2009-12-09 16:43:25 +0000 +++ pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp 2009-12-15 09:41:07 +0000 @@ -48,7 +48,7 @@ updateScenePtr(); shared_ptr<BodyContainer>& bodies = scene->bodies; - Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cellMax-scene->cellMin; + Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cellSize; bool removeUnseenIntrs=(scene->interactions->iterColliderLastRun>=0 && scene->interactions->iterColliderLastRun==scene->currentIteration); #ifdef YADE_OPENMP const long size=scene->interactions->size(); === modified file 'pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp' --- pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp 2009-12-09 17:11:51 +0000 +++ pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp 2009-12-15 09:41:07 +0000 @@ -214,7 +214,7 @@ BBj[i].coord=((BBj[i].flags.hasBB=((bool)bv)) ? (BBj[i].flags.isMin ? bv->min[j] : bv->max[j]) : (b->state->pos[j])) - (periodic ? BBj.cellDim*BBj[i].period : 0.); } else { BBj[i].flags.hasBB=false; /* for vanished body, keep the coordinate as-is, to minimize inversions. */ } // if initializing periodic, shift coords & record the period into BBj[i].period - if(doInitSort && periodic) BBj[i].coord=cellWrap(BBj[i].coord,BBj.cellMin,BBj.cellMax,BBj[i].period); + if(doInitSort && periodic) BBj[i].coord=cellWrap(BBj[i].coord,0,BBj.cellDim,BBj[i].period); } } // for each body, copy its minima and maxima, for quick checks of overlaps later @@ -296,7 +296,7 @@ return x0+(xNorm-period)*(x1-x0); } -// return coordinate wrapped to x0â¦x1, relative to x0; don't care about period +// return coordinate wrapped to 0â¦x1, relative to x0; don't care about period Real InsertionSortCollider::cellWrapRel(const Real x, const Real x0, const Real x1){ Real xNorm=(x-x0)/(x1-x0); return (xNorm-floor(xNorm))*(x1-x0); @@ -309,7 +309,7 @@ const long i=v.norm(_i); const long i_1=v.norm(i-1); //switch period of (i) if the coord is below the lower edge cooridnate-wise and just above the split - if(i==loIdx && v[i].coord<v.cellMin){ v[i].period-=1; v[i].coord+=v.cellDim; loIdx=v.norm(loIdx+1); } + if(i==loIdx && v[i].coord<0){ v[i].period-=1; v[i].coord+=v.cellDim; loIdx=v.norm(loIdx+1); } // coordinate of v[i] used to check inversions // if crossing the split, adjust by cellDim; // if we get below the loIdx however, the v[i].coord will have been adjusted already, no need to do that here @@ -323,7 +323,7 @@ long j1=v.norm(j+1); // OK, now if many bodies move at the same pace through the cell and at one point, there is inversion, // this can happen without any side-effects - if (false && v[j].coord>v.cellMax+v.cellDim){ + if (false && v[j].coord>2*v.cellDim){ // this condition is not strictly necessary, but the loop of insertionSort would have to run more times. // Since size of particle is required to be < .5*cellDim, this would mean simulation explosion anyway LOG_FATAL("Body #"<<v[j].id<<" going faster than 1 cell in one step? Not handled."); @@ -332,7 +332,7 @@ Bounds& vNew(v[j1]); // elt at j+1 being overwritten by the one at j and adjusted vNew=v[j]; // inversions close the the split need special care - if(j==loIdx && vi.coord<v.cellMin) { vi.period-=1; vi.coord+=v.cellDim; loIdx=v.norm(loIdx+1); } + if(j==loIdx && vi.coord<0) { vi.period-=1; vi.coord+=v.cellDim; loIdx=v.norm(loIdx+1); } else if(j1==loIdx) { vNew.period+=1; vNew.coord-=v.cellDim; loIdx=v.norm(loIdx-1); } if(doCollide && viHasBB && v[j].flags.hasBB){ if(vi.id==vNew.id){ // BUG!! @@ -399,7 +399,7 @@ assert(periodic); assert(id1!=id2); // programming error, or weird bodies (too large?) for(int axis=0; axis<3; axis++){ - Real dim=rb->cellMax[axis]-rb->cellMin[axis]; + Real dim=rb->cellSize[axis]; // too big bodies in interaction assert(maxima[3*id1+axis]-minima[3*id1+axis]<.99*dim); assert(maxima[3*id2+axis]-minima[3*id2+axis]<.99*dim); // find body of which when taken as period start will make the gap smaller === modified file 'pkg/common/Engine/GlobalEngine/InsertionSortCollider.hpp' --- pkg/common/Engine/GlobalEngine/InsertionSortCollider.hpp 2009-12-09 17:11:51 +0000 +++ pkg/common/Engine/GlobalEngine/InsertionSortCollider.hpp 2009-12-15 09:41:07 +0000 @@ -53,7 +53,7 @@ Periodicity control =================== c++: - Scene::isPeriodic, Scene::cellMin, Scene::cellMax + Scene::isPeriodic, Scene::cellSize python: O.periodicCell=((0,0,0),(10,10,10) # activates periodic boundary O.periodicCell=() # deactivates periodic boundary @@ -148,7 +148,7 @@ // axis set in the ctor int axis; std::vector<Bounds> vec; - Real cellMin, cellMax, cellDim; + Real cellDim; // cache vector size(), update at every step in action() long size; // index of the lowest coordinate element, before which the container wraps @@ -159,7 +159,7 @@ void updatePeriodicity(Scene* rb){ assert(rb->isPeriodic); assert(axis>=0 && axis <=2); - cellMin=rb->cellMin[axis]; cellMax=rb->cellMax[axis]; cellDim=cellMax-cellMin; + cellDim=rb->cellSize[axis]; } // normalize given index to the right range (wraps around) long norm(long i) const { if(i<0) i+=size; long ret=i%size; assert(ret>=0 && ret<size); return ret;} === modified file 'pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp' --- pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp 2009-12-14 18:49:06 +0000 +++ pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp 2009-12-15 09:41:07 +0000 @@ -111,13 +111,13 @@ /* mostly copied from PeriodicInsertionSortCollider FIXME: common implementation somewhere */ -Real OpenGLRenderingEngine::wrapCell(const Real x, const Real x0, const Real x1){ - Real xNorm=(x-x0)/(x1-x0); - return x0+(xNorm-floor(xNorm))*(x1-x0); +Real OpenGLRenderingEngine::wrapCell(const Real x, const Real x1){ + Real xNorm=(x)/(x1); + return (xNorm-floor(xNorm))*(x1); } Vector3r OpenGLRenderingEngine::wrapCellPt(const Vector3r& pt, Scene* rb){ if(!rb->isPeriodic) return pt; - return Vector3r(wrapCell(pt[0],rb->cellMin[0],rb->cellMax[0]),wrapCell(pt[1],rb->cellMin[1],rb->cellMax[1]),wrapCell(pt[2],rb->cellMin[2],rb->cellMax[2])); + return Vector3r(wrapCell(pt[0],rb->cellSize[0]),wrapCell(pt[1],rb->cellSize[1]),wrapCell(pt[2],rb->cellSize[2])); } void OpenGLRenderingEngine::setBodiesDispInfo(const shared_ptr<Scene>& rootBody){ @@ -148,8 +148,8 @@ if(!rootBody->isPeriodic) return; glPushMatrix(); glColor3v(Vector3r(1,1,0)); - Vector3r cent=.5*(rootBody->cellMin+rootBody->cellMax); Vector3r size=rootBody->cellMax-rootBody->cellMin; - glTranslate(cent[0],cent[1],cent[2]); glScale(size[0],size[1],size[2]); + Vector3r cent=.5*rootBody->cellSize; + glTranslate(cent[0],cent[1],cent[2]); glScale(rootBody->cellSize[0],rootBody->cellSize[1],rootBody->cellSize[2]); glutWireCube(1); glPopMatrix(); } @@ -373,7 +373,7 @@ // if the body goes over the cell margin, draw it in positions where the bbox overlaps with the cell in wire // precondition: pos is inside the cell. if(b->bound && rootBody->isPeriodic){ - const Vector3r& cellMin(rootBody->cellMin); const Vector3r& cellMax(rootBody->cellMax); Vector3r cellSize=cellMax-cellMin; + const Vector3r& cellSize(rootBody->cellSize); // traverse all periodic cells around the body, to see if any of them touches Vector3r halfSize=b->bound->max-b->bound->min; halfSize*=.5; Vector3r pmin,pmax; @@ -382,9 +382,9 @@ if(i[0]==0 && i[1]==0 && i[2]==0) continue; // middle; already rendered above Vector3r pt=pos+Vector3r(cellSize[0]*i[0],cellSize[1]*i[1],cellSize[2]*i[2]); pmin=pt-halfSize; pmax=pt+halfSize; - if(pmin[0]<=cellMax[0] && pmax[0]>=cellMin[0] && - pmin[1]<=cellMax[1] && pmax[1]>=cellMin[1] && - pmin[2]<=cellMax[2] && pmax[2]>=cellMin[2]) { + if(pmin[0]<=cellSize[0] && pmax[0]>=0 && + pmin[1]<=cellSize[1] && pmax[1]>=0 && + pmin[2]<=cellSize[2] && pmax[2]>=0) { glPushMatrix(); glTranslatev(pt); glRotatef(angle*Mathr::RAD_TO_DEG,axis[0],axis[1],axis[2]); === modified file 'pkg/common/RenderingEngine/OpenGLRenderingEngine.hpp' --- pkg/common/RenderingEngine/OpenGLRenderingEngine.hpp 2009-12-11 18:10:19 +0000 +++ pkg/common/RenderingEngine/OpenGLRenderingEngine.hpp 2009-12-15 09:41:07 +0000 @@ -39,8 +39,8 @@ Real normSaw(Real t, Real period){ Real xi=(t-period*((int)(t/period)))/period; /* normalized value, (0-1⪠*/ return (xi<.5?2*xi:2-2*xi); } Real normSquare(Real t, Real period){ Real xi=(t-period*((int)(t/period)))/period; /* normalized value, (0-1⪠*/ return (xi<.5?0:1); } - //! wrap number to interval x0â¦x1 - Real wrapCell(const Real x, const Real x0, const Real x1); + //! wrap number to interval 0â¦x1 + Real wrapCell(const Real x, const Real x1); //! wrap point to inside Scene's cell (identity if !Scene::isPeriodic) Vector3r wrapCellPt(const Vector3r& pt, Scene* rb); void drawPeriodicCell(Scene*); === modified file 'pkg/dem/DataClass/SpherePack.cpp' --- pkg/dem/DataClass/SpherePack.cpp 2009-12-13 20:30:13 +0000 +++ pkg/dem/DataClass/SpherePack.cpp 2009-12-15 09:41:07 +0000 @@ -75,7 +75,7 @@ if(!intSph) continue; pack.push_back(Sph(b->state->pos,intSph->radius)); } - if(rootBody->isPeriodic) { cellSize=rootBody->cellMax-rootBody->cellMin; } + if(rootBody->isPeriodic) { cellSize=rootBody->cellSize; } } long SpherePack::makeCloud(Vector3r mn, Vector3r mx, Real rMean, Real rRelFuzz, size_t num, bool periodic){ === modified file 'pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp' --- pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp 2009-12-09 17:11:51 +0000 +++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp 2009-12-15 09:41:07 +0000 @@ -66,9 +66,8 @@ // account for motion of the periodic boundary, if we remember its last position // its velocity will count as max velocity of bodies // otherwise the collider might not run if only the cell were changing without any particle motion - if(scene->isPeriodic && (prevCellMax!=scene->cellMax || prevCellMin!=scene->cellMin)){ cellChanged=true; maxVelocitySq=(prevCellMax-prevCellMin-(scene->cellMax-scene->cellMin)).SquaredLength()/pow(dt,2); } + if(scene->isPeriodic && (prevCellSize!=scene->cellSize)){ cellChanged=true; maxVelocitySq=(prevCellSize-scene->cellSize).SquaredLength()/pow(dt,2); } else { maxVelocitySq=0; cellChanged=false; } - prevCellSize=prevCellMax-prevCellMin; cellSize=scene->cellMax-scene->cellMin; haveBins=(bool)velocityBins; if(haveBins) velocityBins->binVelSqInitialize(maxVelocitySq); @@ -152,7 +151,7 @@ FOREACH(const Real& thrMaxVSq, threadMaxVelocitySq) { maxVelocitySq=max(maxVelocitySq,thrMaxVSq); } #endif if(haveBins) velocityBins->binVelSqFinalize(); - if(scene->isPeriodic) { prevCellMax=scene->cellMax; prevCellMin=scene->cellMin; } + if(scene->isPeriodic) { prevCellSize=scene->cellSize; } } inline void NewtonIntegrator::leapfrogTranslate(Scene* scene, State* state, const body_id_t& id, const Real& dt ) @@ -162,8 +161,9 @@ // which is then used to compute new position x=ξ(xâ-xâ)+xâ. (per component) // Then we update either velocity by (x-x')/Ît (homotheticCellResize==1) // or position by (x-x') (homotheticCellResize==2) - Vector3r dPos; - if(cellChanged && homotheticCellResize){ for(int i=0; i<3; i++) dPos[i]=((state->pos[i]-prevCellMin[i])/(prevCellSize[i]))*cellSize[i]+scene->cellMin[i]-state->pos[i]; } + // FIXME: wrap state->pos first, then apply the shift; otherwise result might be garbage + Vector3r dPos(Vector3r::ZERO); // initialize to avoid warning; find way to avoid it in a better way + if(cellChanged && homotheticCellResize){ for(int i=0; i<3; i++) dPos[i]=(state->pos[i]/prevCellSize[i])*scene->cellSize[i]-state->pos[i]; } blockTranslateDOFs(state->blockedDOFs, state->accel); state->vel+=dt*state->accel; if(homotheticCellResize==1) state->vel+=dPos/dt; state->pos += state->vel*dt + scene->bex.getMove(id); if(homotheticCellResize==2) state->pos+=dPos; === modified file 'pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp' --- pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp 2009-12-09 17:11:51 +0000 +++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp 2009-12-15 09:41:07 +0000 @@ -50,10 +50,8 @@ Quaternionr DotQ(const Vector3r& angVel, const Quaternionr& Q); inline void blockTranslateDOFs(unsigned blockedDOFs, Vector3r& v); inline void blockRotateDOFs(unsigned blockedDOFs, Vector3r& v); - // cell corners from previous step, used to detect change, find max velocity and update positions if linearCellResize enabled - Vector3r prevCellMax, prevCellMin; - // cache sizes, recomputed before the loop at every step - Vector3r prevCellSize, cellSize; + // cell size from previous step, used to detect change, find max velocity and update positions if linearCellResize enabled + Vector3r prevCellSize; // whether the cell has changed from the previous step bool cellChanged; === modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp' --- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2009-12-06 17:18:16 +0000 +++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2009-12-15 09:41:07 +0000 @@ -33,7 +33,7 @@ } if(maxDisplPerStep<0) maxDisplPerStep=1e-2*charLen; // this should be tuned somehow⦠const long& step=rb->currentIteration; - Vector3r cellSize=rb->cellMax-rb->cellMin; //unused: Real cellVolume=cellSize[0]*cellSize[1]*cellSize[2]; + Vector3r cellSize=rb->cellSize; //unused: Real cellVolume=cellSize[0]*cellSize[1]*cellSize[2]; Vector3r cellArea=Vector3r(cellSize[1]*cellSize[2],cellSize[0]*cellSize[2],cellSize[0]*cellSize[1]); Real minSize=min(cellSize[0],min(cellSize[1],cellSize[2])), maxSize=max(cellSize[0],max(cellSize[1],cellSize[2])); if(minSize<2.1*maxSpan){ throw runtime_error("Minimum cell size is smaller than 2.1*span_of_the_biggest_body! (periodic collider requirement)"); } @@ -42,7 +42,7 @@ sigma=-Vector3r(sumForces[0]/cellArea[0],sumForces[1]/cellArea[1],sumForces[2]/cellArea[2]); LOG_TRACE("Updated sigma="<<sigma<<", avgStiffness="<<avgStiffness); } - Real sigmaGoal=stresses[state]; assert(sigmaGoal>0); + Real sigmaGoal=stresses[state]; assert(sigmaGoal<0); // expansion of cell in this step (absolute length) Vector3r cellGrow(Vector3r::ZERO); // is the stress condition satisfied in all directions? @@ -72,7 +72,7 @@ } } TRVAR4(cellGrow,sigma,sigmaGoal,avgStiffness); - rb->cellMin-=.5*cellGrow; rb->cellMax+=.5*cellGrow; + rb->cellSize+=cellGrow; // handle state transitions if(allStressesOK){ if((step%globalUpdateInt)==0) currUnbalanced=Shop::unbalancedForce(/*useMaxForce=*/false,rb); @@ -91,7 +91,7 @@ void PeriTriaxController::strainStressStiffUpdate(){ // update strain first - const Vector3r cellSize(scene->cellMax-scene->cellMin); + const Vector3r& cellSize(scene->cellSize); for(int i=0; i<3; i++) strain[i]=(cellSize[i]-refSize[i])/refSize[i]; // stress and stiffness Vector3r sumForce(Vector3r::ZERO), sumStiff(Vector3r::ZERO), sumLength(Vector3r::ZERO); @@ -122,7 +122,7 @@ void PeriTriaxController::action(Scene* scene){ if(!scene->isPeriodic){ throw runtime_error("PeriTriaxController run on aperiodic simulation."); } - Vector3r cellSize=scene->cellMax-scene->cellMin; + const Vector3r& cellSize=scene->cellSize; Vector3r cellArea=Vector3r(cellSize[1]*cellSize[2],cellSize[0]*cellSize[2],cellSize[0]*cellSize[1]); // initial updates if(refSize[0]<0) refSize=cellSize; @@ -182,7 +182,7 @@ if(stiff[axis]>0) stress[axis]+=(cellGrow[axis]/refSize[axis])*(stiff[axis]/cellArea[axis]); //-bogusPoisson*(cellGrow[ax1]/refSize[ax1])*(stiff[ax1]/cellArea[ax1])-bogusPoisson*(cellGrow[ax2]/refSize[ax2])*(stiff[ax2]/cellArea[ax2]); } // change cell size now - scene->cellMax+=cellGrow; + scene->cellSize+=cellGrow; strainRate=cellGrow/scene->dt; if(allOk){ === modified file 'py/pack.py' --- py/pack.py 2009-12-09 17:11:51 +0000 +++ py/pack.py 2009-12-15 09:41:07 +0000 @@ -304,10 +304,10 @@ if wantPeri: # x1,y1,z1 already computed above sp=SpherePack() - O.periodicCell=((0,0,0),(x1,y1,z1)) + O.periodicCell=Vector3(x1,y1,z1) #print cloudPorosity,beta,gamma,N100,x1,y1,z1,O.periodicCell #print x1,y1,z1,radius,rRelFuzz - num=sp.makeCloud(O.periodicCell[0],O.periodicCell[1],radius,rRelFuzz,spheresInCell,True) + num=sp.makeCloud(Vector3().ZERO,O.periodicCell,radius,rRelFuzz,spheresInCell,True) O.engines=[BexResetter(),BoundDispatcher([InteractingSphere2AABB()]),InsertionSortCollider(nBins=5,sweepLength=.05*radius),InteractionDispatchers([ef2_Sphere_Sphere_Dem3DofGeom()],[SimpleElasticRelationships()],[Law2_Dem3Dof_Elastic_Elastic()]),PeriIsoCompressor(charLen=radius/5.,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=5,keepProportions=True),NewtonIntegrator(damping=.6)] O.materials.append(GranularMat(young=30e9,frictionAngle=.5,poisson=.3,density=1e3)) for s in sp: O.bodies.append(utils.sphere(s[0],s[1])) @@ -373,8 +373,8 @@ from math import pi O.switchScene(); O.resetThisScene() sp=SpherePack() - O.periodicCell=((0,0,0),Vector3(initSize)) - sp.makeCloud(O.periodicCell[0],O.periodicCell[1],radius,rRelFuzz,int(initSize[0]*initSize[1]*initSize[2]/((4/3.)*pi*radius**3)),True) + O.periodicCell=Vector3(initSize) + sp.makeCloud(Vector3().ZERO,O.periodicCell,radius,rRelFuzz,int(initSize[0]*initSize[1]*initSize[2]/((4/3.)*pi*radius**3)),True) O.engines=[BexResetter(),BoundDispatcher([InteractingSphere2AABB()]),InsertionSortCollider(nBins=2,sweepLength=.05*radius),InteractionDispatchers([Ig2_Sphere_Sphere_Dem3DofGeom()],[SimpleElasticRelationships()],[Law2_Dem3Dof_Elastic_Elastic()]),PeriIsoCompressor(charLen=2*radius,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=20,keepProportions=True),NewtonIntegrator(damping=.8)] O.materials.append(GranularMat(young=30e9,frictionAngle=.1,poisson=.3,density=1e3)) for s in sp: O.bodies.append(utils.sphere(s[0],s[1])) === modified file 'py/yadeWrapper/yadeWrapper.cpp' --- py/yadeWrapper/yadeWrapper.cpp 2009-12-13 22:51:22 +0000 +++ py/yadeWrapper/yadeWrapper.cpp 2009-12-15 09:41:07 +0000 @@ -521,15 +521,13 @@ oa << YadeSimulation; } #endif - void periodicCell_set(python::tuple& t){ - if(python::len(t)==2){ OMEGA.getScene()->cellMin=python::extract<Vector3r>(t[0]); OMEGA.getScene()->cellMax=python::extract<Vector3r>(t[1]); OMEGA.getScene()->isPeriodic=true; } - else if (python::len(t)==0) {OMEGA.getScene()->isPeriodic=false; } - else throw invalid_argument("Must pass either 2 Vector3's to set periodic cell, or () to disable periodicity (got "+lexical_cast<string>(python::len(t))+" instead)."); + void periodicCell_set(Vector3r& size){ + if(size[0]>0 && size[1]>0 && size[2]>0) { OMEGA.getScene()->cellSize=size; OMEGA.getScene()->isPeriodic=true; } + else {OMEGA.getScene()->isPeriodic=false; } } - python::tuple periodicCell_get(){ - vector<Vector3r> ret; - if(OMEGA.getScene()->isPeriodic){ return python::make_tuple(OMEGA.getScene()->cellMin,OMEGA.getScene()->cellMax); } - return python::make_tuple(); + Vector3r periodicCell_get(){ + if(OMEGA.getScene()->isPeriodic) return OMEGA.getScene()->cellSize; + return Vector3r(-1,-1,-1); } void disableGdb(){ signal(SIGSEGV,SIG_DFL); === modified file 'scripts/test/pack-cloud.py' --- scripts/test/pack-cloud.py 2009-12-04 23:27:26 +0000 +++ scripts/test/pack-cloud.py 2009-12-15 09:41:07 +0000 @@ -1,8 +1,8 @@ """ Generate random periodic sphere packing using SpherePack::makeCloud """ from yade import pack p=pack.SpherePack() -O.periodicCell=(0,0,0),(10,10,10) -print p.makeCloud(O.periodicCell[0],O.periodicCell[1],.5,.5,1200,True) +O.periodicCell=Vector3(10,10,10) +print p.makeCloud(Vector3().ZERO,O.periodicCell[1],.5,.5,1200,True) for s in p: O.bodies.append(utils.sphere(s[0],s[1])) O.runEngine(BoundDispatcher([InteractingSphere2AABB()])) === modified file 'scripts/test/periodic-compress.py' --- scripts/test/periodic-compress.py 2009-12-11 18:10:19 +0000 +++ scripts/test/periodic-compress.py 2009-12-15 09:41:07 +0000 @@ -1,7 +1,7 @@ -O.periodicCell=(0,0,0),(20,20,10) +O.periodicCell=Vector3(20,20,10) from yade import pack,log,timing p=pack.SpherePack() -p.makeCloud(O.periodicCell[0],O.periodicCell[1],1,.5,700,True) +p.makeCloud(Vector3().ZERO,O.periodicCell,1,.5,700,True) for sph in p: O.bodies.append(utils.sphere(sph[0],sph[1])) === modified file 'scripts/test/periodic-grow.py' --- scripts/test/periodic-grow.py 2009-12-11 18:10:19 +0000 +++ scripts/test/periodic-grow.py 2009-12-15 09:41:07 +0000 @@ -19,7 +19,7 @@ O.bodies.append(utils.sphere(Vector3(10*random.random(),10*random.random(),10*random.random()),.5+random.random())) cubeSize=20 # absolute positioning of the cell is not important -O.periodicCell=((-.5*cubeSize,-.5*cubeSize,0),(.5*cubeSize,.5*cubeSize,cubeSize)) +O.periodicCell=Vector3(cubeSize,cubeSize,cubeSize) O.dt=utils.PWaveTimeStep() O.saveTmp() from yade import qt @@ -28,12 +28,10 @@ O.run(200,True) for i in range(0,250): O.run(200,True) - mn,mx=O.periodicCell - step=(mx-mn); step=Vector3(.002*step[0],.002*step[1],.002*step[2]) - O.periodicCell=mn+step,mx-step + O.periodicCell=O.periodicCell*.998 if (i%10==0): F,stiff=utils.totalForceInVolume() - dim=mx-mn; A=Vector3(dim[1]*dim[2],dim[0]*dim[2],dim[0]*dim[1]) + dim=O.periodicCell; A=Vector3(dim[1]*dim[2],dim[0]*dim[2],dim[0]*dim[1]) avgStress=sum([F[i]/A[i] for i in 0,1,2])/3. print 'strain',(cubeSize-dim[0])/cubeSize,'avg. stress ',avgStress,'unbalanced ',utils.unbalancedForce() #O.timingEnabled=True; timing.reset(); O.run(200000,True); timing.stats() === modified file 'scripts/test/periodic-simple.py' --- scripts/test/periodic-simple.py 2009-12-11 18:10:19 +0000 +++ scripts/test/periodic-simple.py 2009-12-15 09:41:07 +0000 @@ -24,7 +24,7 @@ O.bodies.append(utils.sphere([0,2,5.5],2)) O.bodies.appendClumped([utils.sphere([0,4,8],.8),utils.sphere([0,5,7],.6)]) # sets up the periodic cell -O.periodicCell=((-5,-5,0),(5,5,10)) +O.periodicCell=Vector3(10,10,10) O.dt=.1*utils.PWaveTimeStep() O.saveTmp() from yade import qt === modified file 'scripts/test/periodic-triax.py' --- scripts/test/periodic-triax.py 2009-12-09 17:11:51 +0000 +++ scripts/test/periodic-triax.py 2009-12-15 09:41:07 +0000 @@ -4,10 +4,10 @@ from yade import * from yade import pack,log,qt #log.setLevel('PeriTriaxController',log.DEBUG) -O.periodicCell=(0,0,0),(.1,.1,.1) +O.periodicCell=Vector3(.1,.1,.1) sp=pack.SpherePack() radius=5e-3 -num=sp.makeCloud(O.periodicCell[0],O.periodicCell[1],radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic +num=sp.makeCloud(Vector3().ZERO,O.periodicCell,radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic O.bodies.append([utils.sphere(s[0],s[1]) for s in sp])
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
