yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02746
[Branch ~yade-dev/yade/trunk] Rev 1906: 1. Use Scene::cell instead of Scene::cellSize; update all code
------------------------------------------------------------
revno: 1906
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Thu 2009-12-17 08:13:18 +0100
message:
1. Use Scene::cell instead of Scene::cellSize; update all code
2. Have O.cellSize (instead of O.periodicCell; add wrapper with warning) and add O.cellShear
3. Add scripts/test/periodic-shear.py (purely for showing the display things, no interactions etc yet!)
added:
scripts/test/periodic-shear.py
modified:
core/Body.cpp
core/Body.hpp
core/Cell.hpp
core/Scene.cpp
core/Scene.hpp
examples/concrete/periodic.py
gui/qt3/GLViewer.cpp
pkg/common/DataClass/InteractingGeometry/Sphere.hpp
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/dem/DataClass/SpherePack.cpp
pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp
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.
=== modified file 'core/Body.cpp'
--- core/Body.cpp 2009-12-11 07:27:29 +0000
+++ core/Body.cpp 2009-12-17 07:13:18 +0000
@@ -25,6 +25,6 @@
#endif
// we must initialize id = 0, otherwise BodyContainer will crash.
-Body::Body(): id(Body::ID_NONE),groupMask(1),clumpId(ID_NONE), state(shared_ptr<State>(new State)), interactingGeometry(shape), boundingVolume(bound){}
-Body::Body(body_id_t newId, int newGroup): id(newId), groupMask(newGroup), clumpId(ID_NONE), state(shared_ptr<State>(new State)), interactingGeometry(shape), boundingVolume(bound){}
+Body::Body(): id(Body::ID_NONE),groupMask(1),clumpId(ID_NONE), state(shared_ptr<State>(new State)){}
+Body::Body(body_id_t newId, int newGroup): id(newId), groupMask(newGroup), clumpId(ID_NONE), state(shared_ptr<State>(new State)){}
=== modified file 'core/Body.hpp'
--- core/Body.hpp 2009-12-13 20:30:13 +0000
+++ core/Body.hpp 2009-12-17 07:13:18 +0000
@@ -91,10 +91,6 @@
/// Bound is used for quick detection of potential interactions, that can be: AABB, K-Dop
shared_ptr<Bound> bound;
- // for backwards-compatibility, will be removed in the future
- __attribute__ ((deprecated)) shared_ptr<Shape>& interactingGeometry;
- __attribute__ ((deprecated)) shared_ptr<Bound>& boundingVolume;
-
/*! isDynamic is true if the state of the body is not modified by a kinematicEngine.
* It is useful for example for collision detection : if two colliding bodies are only
* kinematic then it is useless to modelise their contact */
=== modified file 'core/Cell.hpp'
--- core/Cell.hpp 2009-12-15 09:41:07 +0000
+++ core/Cell.hpp 2009-12-17 07:13:18 +0000
@@ -1,52 +1,74 @@
// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
-#include<yade/lib-base/yadeWm3Extra.hpp>
+// #include<yade/lib-base/yadeWm3Extra.hpp>
+
+// begin yade compatibility
+#ifndef PREFIX
+ #include<yade/lib-miniWm3/Wm3Vector3.h>
+ #include<yade/lib-miniWm3/Wm3Matrix3.h>
+ #define REGISTER_ATTRIBUTES(a,b)
+ #define REGISTER_SERIALIZABLE(a)
+ #define REGISTER_CLASS_AND_BASE(a,b)
+ typedef Wm3::Vector3<double> Vector3r;
+ using Wm3::Vector3;
+ typedef Wm3::Matrix3<double> Matrix3r;
+ typedef double Real;
+#endif
+// end yade compatibility
/*! Periodic cell parameters and routines. Usually instantiated as Scene::cell.
-The cell has size and skew, which are independent.
+The cell has size and shear, which are independent.
*/
-class Cell{
+class Cell: public Serializable{
public:
+ Cell(): size(Vector3r(1,1,1)), shear(Vector3r::ZERO){ updateCache(); }
+ //! size of the cell, projected on axes (for non-zero shear)
Vector3r size;
- Vector3r skew;
- Vector3r refSize;
+ //! 3 non-diagonal components of the shear matrix, ordered by axis normal to the shear plane,
+ //! i.e. (εyz=εzy,εxz=εzx,εxy=εyx). Shear of the cell is always symmetric (transforms
+ //! a box into a parallelepiped).
+ //! See http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/strain.cfm#matrix
+ //! for details
+ Vector3r shear;
// caches
- Vector3r _skewAngle;
- Vector3r _skewSine;
- Matrix3r _skewTrsf;
- Matrix3r _unskewTrsf;
+ Vector3r _shearAngle;
+ Vector3r _shearSine;
+ Matrix3r _shearTrsf;
+ Matrix3r _unshearTrsf;
// should be called before every step
void updateCache(){
for(int i=0; i<3; i++) {
- _skewAngle[i]=atan(skew[i]);
- _skewSine[i]=sin(_skewAngle);
+ _shearAngle[i]=atan(shear[i]);
+ _shearSine[i]=sin(_shearAngle[i]);
}
- _skewTrsf=Matrix3r(1,skew[1],skew[2],skew[0],1,skew[2],skew[0],skew[1],1);
- _unskewTrsf=_skewTrst.Inverse();
+ _shearTrsf=Matrix3r(1,shear[2],shear[1],shear[2],1,shear[0],shear[1],shear[0],1);
+ _unshearTrsf=_shearTrsf.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; }
+ /*! Apply inverse shear on point; to put it inside (unsheared) periodic cell, apply wrapPt on the returned value. */
+ Vector3r unshearPt(const Vector3r& pt){ return _unshearTrsf*pt; }
+ //! Apply shear on point.
+ Vector3r shearPt(const Vector3r& pt){ return _shearTrsf*pt; }
/*! Wrap number to interval 0…sz */
- Real wrapNum(const Real& x, const Real& sz){
- Real norm=x/mx; return (norm-floor(norm))*sz;
+ Real wrapNum(const Real& x, const Real& sz)const{
+ Real norm=x/sz; 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];
+ Real wrapNum(const Real& x, const Real& sz, int& period)const{
+ Real norm=x/sz; period=(int)floor(norm); return (norm-period)*sz;
}
/*! Wrap point to inside the periodic cell; don't compute number of periods wrapped */
- Vector3r wrapPt(const Vector3r pt){
+ Vector3r wrapPt(const Vector3r pt)const{
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 wrapPt(const Vector3r pt, Vector3<int>& period)const{
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_ATTRIBUTES(Serializable,(size)(shear));
REGISTER_CLASS_AND_BASE(Cell,Serializable);
};
REGISTER_SERIALIZABLE(Cell);
=== modified file 'core/Scene.cpp'
--- core/Scene.cpp 2009-12-15 09:41:07 +0000
+++ core/Scene.cpp 2009-12-17 07:13:18 +0000
@@ -33,8 +33,7 @@
// should be elsewhere, probably
bool TimingInfo::enabled=false;
-Scene::Scene():
- bodies(new BodyVector), interactions(new InteractionVecMap){
+Scene::Scene(): bodies(new BodyVector), interactions(new InteractionVecMap){
needsInitializers=true;
currentIteration=0;
simulationTime=0;
@@ -44,9 +43,6 @@
dt=1e-8;
selectedBody=-1;
isPeriodic=false;
- cellSize=Vector3r::ZERO;
- // FIXME: move SceneShape to core and create it here right away
- // shape=shared_ptr<Shape>(new SceneShape);
// fill default tags
struct passwd* pw;
@@ -83,6 +79,8 @@
bex.resize(bodies->size());
needsInitializers=false;
}
+ // update cell data
+ if(isPeriodic) cell.updateCache();
//bex.reset(); // uncomment if PhysicalActionContainerReseter is removed
bool TimingInfo_enabled=TimingInfo::enabled; // cache the value, so that when it is changed inside the step, the engine that was just running doesn't get bogus values
TimingInfo::delta last=TimingInfo::getNow(); // actually does something only if TimingInfo::enabled, no need to put the condition here
=== modified file 'core/Scene.hpp'
--- core/Scene.hpp 2009-12-15 09:41:07 +0000
+++ core/Scene.hpp 2009-12-17 07:13:18 +0000
@@ -11,6 +11,7 @@
#pragma once
#include"Body.hpp"
+#include"Cell.hpp"
#include"BodyContainer.hpp"
#include"Engine.hpp"
#include"Material.hpp"
@@ -43,10 +44,12 @@
BexContainer bex;
- vector<shared_ptr<Serializable> > miscParams; // will set static parameters during deserialization (primarily for GLDraw functors which otherwise have no attribute access)
- //! tags like mp3 tags: author, date, version, description etc.
+ //! information on periodicity; only should be used if Scene::isPeriodic
+ Cell cell;
+ //! store for arbitrary Serializable objects; will set static parameters during deserialization (primarily for GLDraw functors which otherwise have no attribute access)
+ vector<shared_ptr<Serializable> > miscParams; //! tags like mp3 tags: author, date, version, description etc.
list<string> tags;
- //! "hash maps" of display parameters
+ //! "hash maps" of display parameters (since yade::serialization has no support for maps, emulate it via vector of strings in format key=value)
vector<shared_ptr<DisplayParameters> > dispParams;
shared_ptr<GroupRelationData> grpRelationData;
@@ -65,7 +68,6 @@
long stopAtIteration;
Real stopAtVirtTime;
Real stopAtRealTime;
- Vector3r cellSize;
bool isPeriodic;
bool needsInitializers;
@@ -87,8 +89,8 @@
(currentIteration)
(simulationTime)
(stopAtIteration)
- (cellSize)
(isPeriodic)
+ (cell)
);
REGISTER_CLASS_AND_BASE(Scene,Serializable);
};
=== modified file 'examples/concrete/periodic.py'
--- examples/concrete/periodic.py 2009-12-15 09:41:07 +0000
+++ examples/concrete/periodic.py 2009-12-17 07:13:18 +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=sp.cellSize
+O.cellSize=sp.cellSize
axis=2
ax1=(axis+1)%3
ax2=(axis+2)%3
=== modified file 'gui/qt3/GLViewer.cpp'
--- gui/qt3/GLViewer.cpp 2009-12-15 17:47:47 +0000
+++ gui/qt3/GLViewer.cpp 2009-12-17 07:13:18 +0000
@@ -336,10 +336,10 @@
}
/* Center the scene such that periodic cell is contained in the view */
void GLViewer::centerPeriodic(){
- Scene* rb=Omega::instance().getScene().get();
- assert(rb->isPeriodic);
- Vector3r center=.5*rb->cellSize;
- Vector3r halfSize=.5*rb->cellSize;
+ Scene* scene=Omega::instance().getScene().get();
+ assert(scene->isPeriodic);
+ Vector3r center=.5*scene->cell.size;
+ Vector3r halfSize=.5*scene->cell.size;
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]));
@@ -355,16 +355,16 @@
* "central" (where most bodies is) part very small or even invisible.
*/
void GLViewer::centerMedianQuartile(){
- Scene* rb=Omega::instance().getScene().get();
- if(rb->isPeriodic){ centerPeriodic(); return; }
- long nBodies=rb->bodies->size();
+ Scene* scene=Omega::instance().getScene().get();
+ if(scene->isPeriodic){ centerPeriodic(); return; }
+ long nBodies=scene->bodies->size();
if(nBodies<4) {
LOG_INFO("Less than 4 bodies, median makes no sense; calling centerScene() instead.");
return centerScene();
}
std::vector<Real> coords[3];
for(int i=0;i<3;i++)coords[i].reserve(nBodies);
- FOREACH(const shared_ptr<Body>& b, *rb->bodies){
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
for(int i=0; i<3; i++) coords[i].push_back(b->state->pos[i]);
}
Vector3r median,interQuart;
@@ -630,7 +630,10 @@
}
if(drawGridXYZ[0] || drawGridXYZ[1] || drawGridXYZ[2]){
glColor3v(Vector3r(1,1,0));
- QGLViewer::drawText(x,y,std::string("grid: "+boost::lexical_cast<std::string>(gridStep))+(grid_subdivision?(std::string(" ,subdivision: "+boost::lexical_cast<std::string>(gridStep*0.1))):(std::string(""))));
+ ostringstream oss;
+ oss<<"grid: "<<setprecision(4)<<gridStep;
+ if(grid_subdivision) oss<<" (minor "<<setprecision(4)<<gridStep*.1<<")";
+ QGLViewer::drawText(x,y,oss.str());
y-=lineHt;
}
}
=== modified file 'pkg/common/DataClass/InteractingGeometry/Sphere.hpp'
--- pkg/common/DataClass/InteractingGeometry/Sphere.hpp 2009-12-13 20:30:13 +0000
+++ pkg/common/DataClass/InteractingGeometry/Sphere.hpp 2009-12-17 07:13:18 +0000
@@ -16,7 +16,7 @@
Real radius;
Sphere();
- Sphere(Real _radius): radius(_radius){}
+ Sphere(Real _radius): radius(_radius){ createIndex(); }
virtual ~Sphere ();
REGISTER_ATTRIBUTES(Shape,(radius));
=== modified file 'pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp'
--- pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2009-12-15 09:41:07 +0000
+++ pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2009-12-17 07:13:18 +0000
@@ -37,7 +37,7 @@
geomDispatcher->updateScenePtr();
physDispatcher->updateScenePtr();
lawDispatcher->updateScenePtr();
- Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cellSize;
+ Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cell.size;
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-15 09:41:07 +0000
+++ pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp 2009-12-17 07:13:18 +0000
@@ -48,7 +48,7 @@
updateScenePtr();
shared_ptr<BodyContainer>& bodies = scene->bodies;
- Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cellSize;
+ Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cell.size;
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-15 09:41:07 +0000
+++ pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp 2009-12-17 07:13:18 +0000
@@ -147,7 +147,7 @@
// update periodicity
assert(BB[0].axis==0); assert(BB[1].axis==1); assert(BB[2].axis==2);
- if(periodic) for(int i=0; i<3; i++) BB[i].updatePeriodicity(rb);
+ if(periodic) for(int i=0; i<3; i++) BB[i].updatePeriodicity(scene);
#ifdef COLLIDE_STRIDED
// get the BoundDispatcher and turn it off; we will call it ourselves
@@ -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->cellSize[axis];
+ Real dim=rb->cell.size[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-15 09:41:07 +0000
+++ pkg/common/Engine/GlobalEngine/InsertionSortCollider.hpp 2009-12-17 07:13:18 +0000
@@ -156,10 +156,10 @@
Bounds& operator[](long idx){ assert(idx<size && idx>=0); return vec[idx]; }
const Bounds& operator[](long idx) const { assert(idx<size && idx>=0); return vec[idx]; }
// update number of bodies, periodic properties and size from Scene
- void updatePeriodicity(Scene* rb){
- assert(rb->isPeriodic);
+ void updatePeriodicity(Scene* scene){
+ assert(scene->isPeriodic);
assert(axis>=0 && axis <=2);
- cellDim=rb->cellSize[axis];
+ cellDim=scene->cell.size[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-15 09:41:07 +0000
+++ pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp 2009-12-17 07:13:18 +0000
@@ -117,7 +117,7 @@
}
Vector3r OpenGLRenderingEngine::wrapCellPt(const Vector3r& pt, Scene* rb){
if(!rb->isPeriodic) return pt;
- return Vector3r(wrapCell(pt[0],rb->cellSize[0]),wrapCell(pt[1],rb->cellSize[1]),wrapCell(pt[2],rb->cellSize[2]));
+ return Vector3r(wrapCell(pt[0],rb->cell.size[0]),wrapCell(pt[1],rb->cell.size[1]),wrapCell(pt[2],rb->cell.size[2]));
}
void OpenGLRenderingEngine::setBodiesDispInfo(const shared_ptr<Scene>& rootBody){
@@ -144,12 +144,20 @@
}
// draw periodic cell, if active
-void OpenGLRenderingEngine::drawPeriodicCell(Scene* rootBody){
- if(!rootBody->isPeriodic) return;
+void OpenGLRenderingEngine::drawPeriodicCell(Scene* scene){
+ if(!scene->isPeriodic) return;
+ const Vector3r& shear(scene->cell.shear);
+ // shear center (moves when sheared)
+ Vector3r cent=scene->cell._shearTrsf*(.5*scene->cell.size);
+ // see http://www.songho.ca/opengl/gl_transform.html#matrix
+ GLdouble cellMat[16]={
+ scene->cell.size[0],shear[2],shear[1],cent[0],
+ shear[2],scene->cell.size[1],shear[0],cent[1],
+ shear[1],shear[0],scene->cell.size[2],cent[2],
+ 0,0,0,1};
+ glColor3v(Vector3r(1,1,0));
glPushMatrix();
- glColor3v(Vector3r(1,1,0));
- Vector3r cent=.5*rootBody->cellSize;
- glTranslate(cent[0],cent[1],cent[2]); glScale(rootBody->cellSize[0],rootBody->cellSize[1],rootBody->cellSize[2]);
+ glMultTransposeMatrixd(cellMat);
glutWireCube(1);
glPopMatrix();
}
@@ -373,7 +381,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& cellSize(rootBody->cellSize);
+ const Vector3r& cellSize(rootBody->cell.size);
// 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;
=== modified file 'pkg/dem/DataClass/SpherePack.cpp'
--- pkg/dem/DataClass/SpherePack.cpp 2009-12-15 09:41:07 +0000
+++ pkg/dem/DataClass/SpherePack.cpp 2009-12-17 07:13:18 +0000
@@ -69,13 +69,13 @@
void SpherePack::fromSimulation() {
pack.clear();
- Scene* rootBody=Omega::instance().getScene().get();
- FOREACH(const shared_ptr<Body>& b, *rootBody->bodies){
+ Scene* scene=Omega::instance().getScene().get();
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
shared_ptr<Sphere> intSph=dynamic_pointer_cast<Sphere>(b->shape);
if(!intSph) continue;
pack.push_back(Sph(b->state->pos,intSph->radius));
}
- if(rootBody->isPeriodic) { cellSize=rootBody->cellSize; }
+ if(scene->isPeriodic) { cellSize=scene->cell.size; }
}
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-15 09:41:07 +0000
+++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp 2009-12-17 07:13:18 +0000
@@ -66,7 +66,7 @@
// 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 && (prevCellSize!=scene->cellSize)){ cellChanged=true; maxVelocitySq=(prevCellSize-scene->cellSize).SquaredLength()/pow(dt,2); }
+ if(scene->isPeriodic && (prevCellSize!=scene->cell.size)){ cellChanged=true; maxVelocitySq=(prevCellSize-scene->cell.size).SquaredLength()/pow(dt,2); }
else { maxVelocitySq=0; cellChanged=false; }
haveBins=(bool)velocityBins;
if(haveBins) velocityBins->binVelSqInitialize(maxVelocitySq);
@@ -151,7 +151,7 @@
FOREACH(const Real& thrMaxVSq, threadMaxVelocitySq) { maxVelocitySq=max(maxVelocitySq,thrMaxVSq); }
#endif
if(haveBins) velocityBins->binVelSqFinalize();
- if(scene->isPeriodic) { prevCellSize=scene->cellSize; }
+ if(scene->isPeriodic) { prevCellSize=scene->cell.size; }
}
inline void NewtonIntegrator::leapfrogTranslate(Scene* scene, State* state, const body_id_t& id, const Real& dt )
@@ -163,7 +163,7 @@
// or position by (x-x') (homotheticCellResize==2)
// 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]; }
+ if(cellChanged && homotheticCellResize){ for(int i=0; i<3; i++) dPos[i]=(state->pos[i]/prevCellSize[i])*scene->cell.size[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/PeriIsoCompressor.cpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2009-12-15 09:41:07 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2009-12-17 07:13:18 +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->cellSize; //unused: Real cellVolume=cellSize[0]*cellSize[1]*cellSize[2];
+ Vector3r cellSize=rb->cell.size; //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)"); }
@@ -72,7 +72,7 @@
}
}
TRVAR4(cellGrow,sigma,sigmaGoal,avgStiffness);
- rb->cellSize+=cellGrow;
+ rb->cell.size+=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->cellSize);
+ const Vector3r& cellSize(scene->cell.size);
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."); }
- const Vector3r& cellSize=scene->cellSize;
+ const Vector3r& cellSize=scene->cell.size;
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->cellSize+=cellGrow;
+ scene->cell.size+=cellGrow;
strainRate=cellGrow/scene->dt;
if(allOk){
=== modified file 'py/pack.py'
--- py/pack.py 2009-12-15 09:41:07 +0000
+++ py/pack.py 2009-12-17 07:13:18 +0000
@@ -304,10 +304,10 @@
if wantPeri:
# x1,y1,z1 already computed above
sp=SpherePack()
- O.periodicCell=Vector3(x1,y1,z1)
- #print cloudPorosity,beta,gamma,N100,x1,y1,z1,O.periodicCell
+ O.cellSize=Vector3(x1,y1,z1)
+ #print cloudPorosity,beta,gamma,N100,x1,y1,z1,O.cellSize
#print x1,y1,z1,radius,rRelFuzz
- num=sp.makeCloud(Vector3().ZERO,O.periodicCell,radius,rRelFuzz,spheresInCell,True)
+ num=sp.makeCloud(Vector3().ZERO,O.cellSize,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=Vector3(initSize)
- sp.makeCloud(Vector3().ZERO,O.periodicCell,radius,rRelFuzz,int(initSize[0]*initSize[1]*initSize[2]/((4/3.)*pi*radius**3)),True)
+ O.cellSize=Vector3(initSize)
+ sp.makeCloud(Vector3().ZERO,O.cellSize,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-15 09:41:07 +0000
+++ py/yadeWrapper/yadeWrapper.cpp 2009-12-17 07:13:18 +0000
@@ -521,14 +521,20 @@
oa << YadeSimulation;
}
#endif
- void periodicCell_set(Vector3r& size){
- if(size[0]>0 && size[1]>0 && size[2]>0) { OMEGA.getScene()->cellSize=size; OMEGA.getScene()->isPeriodic=true; }
+ void cellSize_set_depr(const Vector3r& size){ LOG_WARN("O.periodicCell is deprecated (and will be removed), use O.cellSize instead."); cellSize_set(size); }
+ Vector3r cellSize_get_depr(){ LOG_WARN("O.periodicCell is deprecated (and will be removed), use O.cellSize instead."); return cellSize_get(); }
+
+ void cellSize_set(const Vector3r& size){
+ if(size[0]>0 && size[1]>0 && size[2]>0) { OMEGA.getScene()->cell.size=size; OMEGA.getScene()->isPeriodic=true; }
else {OMEGA.getScene()->isPeriodic=false; }
+ OMEGA.getScene()->cell.updateCache();
}
- Vector3r periodicCell_get(){
- if(OMEGA.getScene()->isPeriodic) return OMEGA.getScene()->cellSize;
+ Vector3r cellSize_get(){
+ if(OMEGA.getScene()->isPeriodic) return OMEGA.getScene()->cell.size;
return Vector3r(-1,-1,-1);
}
+ void cellShear_set(const Vector3r& s){ OMEGA.getScene()->cell.shear=s; /* DO NOT FORGET!! */ OMEGA.getScene()->cell.updateCache(); }
+ Vector3r cellShear_get(){ return OMEGA.getScene()->cell.shear; }
void disableGdb(){
signal(SIGSEGV,SIG_DFL);
signal(SIGABRT,SIG_DFL);
@@ -733,7 +739,9 @@
.add_property("timingEnabled",&pyOmega::timingEnabled_get,&pyOmega::timingEnabled_set,"Globally enable/disable timing services (see documentation of yade.timing).")
.add_property("bexSyncCount",&pyOmega::bexSyncCount_get,&pyOmega::bexSyncCount_set,"Counter for number of syncs in BexContainer, for profiling purposes.")
.add_property("numThreads",&pyOmega::numThreads_get /* ,&pyOmega::numThreads_set*/ ,"Get maximum number of threads openMP can use.")
- .add_property("periodicCell",&pyOmega::periodicCell_get,&pyOmega::periodicCell_set, "Get/set periodic cell minimum and maximum (tuple of 2 Vector3's), or () for no periodicity.")
+ .add_property("periodicCell",&pyOmega::cellSize_get_depr,&pyOmega::cellSize_set_depr,"Deprecated alias for cellSize")
+ .add_property("cellSize",&pyOmega::cellSize_get,&pyOmega::cellSize_set, "Get/set periodic cell minimum and maximum (tuple of 2 Vector3's), or () for no periodicity.")
+ .add_property("cellShear",&pyOmega::cellShear_get,&pyOmega::cellShear_set,"Get/set periodic cell shear.")
.def("exitNoBacktrace",&pyOmega::exitNoBacktrace,omega_exitNoBacktrace_overloads("Disable SEGV handler and exit."))
.def("disableGdb",&pyOmega::disableGdb,"Revert SEGV and ABRT handlers to system defaults.")
#ifdef YADE_BOOST_SERIALIZATION
=== modified file 'scripts/test/pack-cloud.py'
--- scripts/test/pack-cloud.py 2009-12-15 09:41:07 +0000
+++ scripts/test/pack-cloud.py 2009-12-17 07:13:18 +0000
@@ -1,8 +1,8 @@
""" Generate random periodic sphere packing using SpherePack::makeCloud """
from yade import pack
p=pack.SpherePack()
-O.periodicCell=Vector3(10,10,10)
-print p.makeCloud(Vector3().ZERO,O.periodicCell[1],.5,.5,1200,True)
+O.cellSize=Vector3(10,10,10)
+print p.makeCloud(Vector3().ZERO,O.cellSize[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-15 09:41:07 +0000
+++ scripts/test/periodic-compress.py 2009-12-17 07:13:18 +0000
@@ -1,7 +1,7 @@
-O.periodicCell=Vector3(20,20,10)
+O.cellSize=Vector3(20,20,10)
from yade import pack,log,timing
p=pack.SpherePack()
-p.makeCloud(Vector3().ZERO,O.periodicCell,1,.5,700,True)
+p.makeCloud(Vector3().ZERO,O.cellSize,1,.5,700,True)
for sph in p:
O.bodies.append(utils.sphere(sph[0],sph[1]))
@@ -21,7 +21,7 @@
]
O.dt=utils.PWaveTimeStep()
O.saveTmp()
-print O.periodicCell
+print O.cellSize
from yade import qt; qt.Controller(); qt.View()
O.run()
O.wait()
=== modified file 'scripts/test/periodic-grow.py'
--- scripts/test/periodic-grow.py 2009-12-15 09:41:07 +0000
+++ scripts/test/periodic-grow.py 2009-12-17 07:13:18 +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=Vector3(cubeSize,cubeSize,cubeSize)
+O.cellSize=Vector3(cubeSize,cubeSize,cubeSize)
O.dt=utils.PWaveTimeStep()
O.saveTmp()
from yade import qt
@@ -28,10 +28,10 @@
O.run(200,True)
for i in range(0,250):
O.run(200,True)
- O.periodicCell=O.periodicCell*.998
+ O.cellSize=O.cellSize*.998
if (i%10==0):
F,stiff=utils.totalForceInVolume()
- dim=O.periodicCell; A=Vector3(dim[1]*dim[2],dim[0]*dim[2],dim[0]*dim[1])
+ dim=O.cellSize; 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()
=== added file 'scripts/test/periodic-shear.py'
--- scripts/test/periodic-shear.py 1970-01-01 00:00:00 +0000
+++ scripts/test/periodic-shear.py 2009-12-17 07:13:18 +0000
@@ -0,0 +1,11 @@
+O.cellSize=Vector3(1,1,1)
+O.bodies.append(utils.sphere([0,0,0],.05))
+g=0.
+import yade.qt,time
+v=yade.qt.View()
+v.axes=True
+v.grid=(True,True,True)
+while True:
+ O.cellShear=Vector3(.2*sin(g),.2*cos(pi*g),.2*sin(2*g)+.2*cos(3*g))
+ time.sleep(0.001)
+ g+=1e-3
=== modified file 'scripts/test/periodic-simple.py'
--- scripts/test/periodic-simple.py 2009-12-15 09:41:07 +0000
+++ scripts/test/periodic-simple.py 2009-12-17 07:13:18 +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=Vector3(10,10,10)
+O.cellSize=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-15 09:41:07 +0000
+++ scripts/test/periodic-triax.py 2009-12-17 07:13:18 +0000
@@ -4,10 +4,10 @@
from yade import *
from yade import pack,log,qt
#log.setLevel('PeriTriaxController',log.DEBUG)
-O.periodicCell=Vector3(.1,.1,.1)
+O.cellSize=Vector3(.1,.1,.1)
sp=pack.SpherePack()
radius=5e-3
-num=sp.makeCloud(Vector3().ZERO,O.periodicCell,radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic
+num=sp.makeCloud(Vector3().ZERO,O.cellSize,radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic
O.bodies.append([utils.sphere(s[0],s[1]) for s in sp])
Follow ups
