yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02708
[Branch ~yade-dev/yade/trunk] Rev 1902: 1. Remove Scene::{cellMin,cellMax}, use just Scene::cellSize; update all related code and scripts
------------------------------------------------------------
revno: 1902
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
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 <eudoxos@xxxxxxxx>
+#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])
