yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02844
[Branch ~yade-dev/yade/trunk] Rev 1925: 1. Add cell strain rate integration to the Scene loop.
------------------------------------------------------------
revno: 1925
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Fri 2009-12-25 22:58:25 +0100
message:
1. Add cell strain rate integration to the Scene loop.
2. Remove Cell::Hsize and friends, remove VELGRAD; unify both approaches.
3. Update some scripts, fix Shop::flipCell (demo in scripts/test/periodic-triax-velgrad.py)
modified:
core/Cell.hpp
core/Scene.cpp
pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp
pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp
pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp
pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp
pkg/dem/meta/Shop.cpp
py/_utils.cpp
py/yadeWrapper/yadeWrapper.cpp
scripts/test/periodic-compress.py
scripts/test/periodic-grow.py
scripts/test/periodic-triax-velgrad.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/Cell.hpp'
--- core/Cell.hpp 2009-12-25 15:10:35 +0000
+++ core/Cell.hpp 2009-12-25 21:58:25 +0000
@@ -1,24 +1,6 @@
// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
// #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;
- class Serializable{};
-#endif
-// end yade compatibility
-
-
-#define VELGRAD
-
/*! Periodic cell parameters and routines. Usually instantiated as Scene::cell.
The cell has size and shear, which are independent.
@@ -26,33 +8,15 @@
*/
class Cell: public Serializable{
public:
- Cell(): refSize(Vector3r(1,1,1)), strain(Matrix3r::ZERO)
-#ifdef VELGRAD
- , velGrad(Matrix3r::ZERO)
- , Hsize(Matrix3r::ZERO)
- ,_shearTrsfMatrix(Matrix3r::IDENTITY)
- { updateCache(0); }
-#else
- { updateCache(); }
-#endif
+ Cell(): refSize(Vector3r(1,1,1)), strain(Matrix3r::ZERO), velGrad(Matrix3r::ZERO){ integrateAndUpdate(0); }
//! size of the cell, projected on axes (for non-zero shear)
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;
+ //! Strain matrix (current total strain)
+ //! Ordering see http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/strain.cfm#matrix
Matrix3r strain;
-#ifdef VELGRAD
Matrix3r velGrad;
- Matrix3r Hsize;
- Matrix3r _shearIncrt;
-#endif
-
- //! reference values of size and shear (for rendering, mainly)
- // Vector3r refShear, refCenter;
+ //Matrix3r Hsize;
//! Get current size (refSize × normal strain)
const Vector3r& getSize() const { return _size; }
@@ -66,10 +30,8 @@
const Matrix3r& getShearTrsfMatrix() const { return _shearTrsfMatrix; }
//! inverse of getShearTrsfMatrix().
const Matrix3r& getUnshearTrsfMatrix() const {return _unshearTrsfMatrix;}
-#ifdef VELGRAD
//! transformation increment matrix applying arbitrary field
- const Matrix3r& getShearIncrMatrix() const { return _shearIncrt; }
-#endif
+ const Matrix3r& getStrainIncrMatrix() const { return _strainInc; }
/*! return pointer to column-major OpenGL 4x4 matrix applying pure shear. This matrix is suitable as argument for glMultMatrixd.
@@ -87,13 +49,12 @@
//! Whether any shear (non-diagonal) component of the strain matrix is nonzero.
bool hasShear() const {return _hasShear; }
- // caches
+ // caches; private
private:
+ Matrix3r _strainInc;
Vector3r _size;
bool _hasShear;
Matrix3r _shearTrsfMatrix, _unshearTrsfMatrix, _cosMatrix;
- // Vector3r _shearAngle, _shearSin, _shearCos;
- //Matrix3r _shearTrsf, _unshearTrsf;
double _glShearTrsfMatrix[16];
void fillGlShearTrsfMatrix(double m[16]){
m[0]=1; m[4]=strain[0][1]; m[8]=strain[0][2]; m[12]=0;
@@ -104,56 +65,38 @@
public:
- // should be called before every step
-#ifdef VELGRAD
- void updateCache(double dt){
- //initialize Hsize for "lazy" default scripts, after that Hsize is free to change
- if (refSize[0]!=1 && Hsize[0][0]==0) {Hsize[0][0]=refSize[0]; Hsize[1][1]=refSize[1]; Hsize[2][2]=refSize[2];}
- //incremental disp gradient
- _shearIncrt=dt*velGrad;
- //update Hsize (redundant with total transformation perhaps)
- Hsize=Hsize+_shearIncrt*Hsize;
- //total transformation
- _shearTrsfMatrix = _shearTrsfMatrix + _shearIncrt*_shearTrsfMatrix;// M = (Id+G).M = F.M
- //compatibility with Vaclav's code :
- _size[0]=Hsize[0][0]; _size[1]=Hsize[1][1]; _size[2]=Hsize[2][2];
- _hasShear = true;
+ //! "integrate" velGrad, update cached values used by public getter
+ void integrateAndUpdate(Real dt){
+ #if 0
+ //initialize Hsize for "lazy" default scripts, after that Hsize is free to change
+ if (refSize[0]!=1 && Hsize[0][0]==0) {Hsize[0][0]=refSize[0]; Hsize[1][1]=refSize[1]; Hsize[2][2]=refSize[2];}
+ //update Hsize (redundant with total transformation perhaps)
+ Hsize=Hsize+_shearInc*Hsize;
+ #endif
+
+ //incremental displacement gradient
+ _strainInc=dt*velGrad;
+ // total transformation; M = (Id+G).M = F.M
+ strain+=_strainInc*_shearTrsfMatrix;
+ #if 0
+ cerr<<"velGrad "; for(int i=0; i<3; i++) for(int j=0; j<3; j++) cerr<<velGrad[i][j]<<endl;
+ cerr<<"_strainInc "; for(int i=0; i<3; i++) for(int j=0; j<3; j++) cerr<<_strainInc[i][j]<<endl;
+ cerr<<"strain "; for(int i=0; i<3; i++) for(int j=0; j<3; j++) cerr<<strain[i][j]<<endl;
+ #endif
+ // pure shear transformation: strain with 1s on the diagonal
+ _shearTrsfMatrix=strain; _shearTrsfMatrix[0][0]=_shearTrsfMatrix[1][1]=_shearTrsfMatrix[2][2]=1.;
+ // pure unshear transformation
_unshearTrsfMatrix=_shearTrsfMatrix.Inverse();
- strain=_shearTrsfMatrix;
-#else
- void updateCache(){
+ // some parts can branch depending on presence/absence of shear
+ _hasShear=!(strain[0][1]==0 && strain[0][2]==0 && strain[1][0]==0 && strain[1][2]==0 && strain[2][0]==0 && strain[2][1]==0);
+ // current cell size (in units on physical space axes)
_size=refSize+diagMult(getExtensionalStrain(),refSize);
- _hasShear=!(strain[0][1]==0 && strain[0][2]==0 && strain[1][0]==0 && strain[1][2]==0 && strain[2][0]==0 && strain[2][1]==0);
- _shearTrsfMatrix=strain;
- _shearTrsfMatrix[0][0]=_shearTrsfMatrix[1][1]=_shearTrsfMatrix[2][2]=1.;
- _unshearTrsfMatrix=_shearTrsfMatrix.Inverse();
-#endif
+ // OpenGL shear matrix (used frequently)
fillGlShearTrsfMatrix(_glShearTrsfMatrix);
+ // precompute cosines of angles, used for enlarge factor when computing Aabb's
for(int i=0; i<3; i++) for(int j=0; j<3; j++) _cosMatrix[i][j]=(i==j?0:cos(atan(strain[i][j])));
}
-// #ifdef VELGRAD
-// void updateCache(){ updateCache(0);}
-// #endif
-
- // doesn't seem to be really useful
- #if 0
- /*! Prepare OpenGL matrix with current shear and given translation and scale.
-
- This matrix is ccumulated product of this sequence:
-
- glTranslatev(translation);
- glScalev(scale);
- glMultMatrixd(scene->cell->_glShearMatrix);
-
- */
- void glTrsfMatrix(double m[16], const Vector3r& translation, const Vector3r& scale){
- m[0]=scale[0]; m[4]=scale[2]*shear[2]; m[8]=scale[1]*shear[1]; m[12]=translation[0];
- m[1]=scale[2]*shear[2]; m[5]=scale[1]; m[9]=scale[0]*shear[0]; m[13]=translation[1];
- m[2]=scale[1]*shear[1]; m[6]=shear[0]; m[10]=scale[2]; m[14]=translation[2];
- m[3]=0; m[7]=0; m[11]=0; m[15]=1;
- }
- #endif
/*! Return point inside periodic cell, even if shear is applied */
Vector3r wrapShearedPt(const Vector3r& pt){ return shearPt(wrapPt(unshearPt(pt))); }
/*! Return point inside periodic cell, even if shear is applied; store cell coordinates in period. */
@@ -170,7 +113,6 @@
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;
}
-
/*! Wrap number to interval 0…sz */
static Real wrapNum(const Real& x, const Real& sz){
Real norm=x/sz; return (norm-floor(norm))*sz;
@@ -179,16 +121,8 @@
static Real wrapNum(const Real& x, const Real& sz, int& period){
Real norm=x/sz; period=(int)floor(norm); return (norm-period)*sz;
}
-#ifdef VELGRAD
- void postProcessAttributes(bool deserializing){ if(deserializing) updateCache(0); }
-#else
- void postProcessAttributes(bool deserializing){ if(deserializing) updateCache(); }
-#endif
- REGISTER_ATTRIBUTES(Serializable,(refSize)(strain)
-#ifdef VELGRAD
- (velGrad)(Hsize)
-#endif
- );
+ void postProcessAttributes(bool deserializing){ if(deserializing) integrateAndUpdate(0); }
+ REGISTER_ATTRIBUTES(Serializable,(refSize)(strain)(velGrad));
REGISTER_CLASS_AND_BASE(Cell,Serializable);
};
REGISTER_SERIALIZABLE(Cell);
=== modified file 'core/Scene.cpp'
--- core/Scene.cpp 2009-12-25 14:46:48 +0000
+++ core/Scene.cpp 2009-12-25 21:58:25 +0000
@@ -80,12 +80,7 @@
forces.resize(bodies->size());
needsInitializers=false;
}
- // update cell data
-#ifndef VELGRAD
- if(isPeriodic) cell->updateCache();
-#else
- if(isPeriodic) cell->updateCache(dt);
-#endif
+ if(isPeriodic) cell->integrateAndUpdate(dt);
//forces.reset(); // uncomment if ForceResetter 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 'pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp'
--- pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp 2009-12-21 22:19:11 +0000
+++ pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp 2009-12-25 21:58:25 +0000
@@ -382,8 +382,10 @@
if(pmin[0]<=cellSize[0] && pmax[0]>=0 &&
pmin[1]<=cellSize[1] && pmax[1]>=0 &&
pmin[2]<=cellSize[2] && pmax[2]>=0) {
+ Vector3r pt=scene->cell->shearPt(pos2);
+ if(pointClipped(pt)) continue;
glPushMatrix();
- glTranslatev(scene->cell->shearPt(pos2));
+ glTranslatev(pt);
glRotatef(angle*Mathr::RAD_TO_DEG,axis[0],axis[1],axis[2]);
shapeDispatcher(b->shape,b->state,/*Body_wire*/ true, viewInfo);
glPopMatrix();
=== modified file 'pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp'
--- pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp 2009-12-25 14:46:48 +0000
+++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp 2009-12-25 21:58:25 +0000
@@ -63,11 +63,6 @@
void NewtonIntegrator::action(Scene*)
{
- // only temporary
-#ifndef VELGRAD
- if(homotheticCellResize) throw runtime_error("Homothetic resizing not yet implemented with new, core/Cell.hpp based cell (extension+shear).");
-#endif
-
scene->forces.sync();
Real dt=scene->dt;
// account for motion of the periodic boundary, if we remember its last position
@@ -163,26 +158,18 @@
inline void NewtonIntegrator::leapfrogTranslate(Scene* scene, State* state, const body_id_t& id, const Real& dt )
{
- // for homothetic resize of the cell (if enabled), compute the position difference of the homothetic transformation
- // the term ξ=(x'-xâ‚€')/(xâ‚'-xâ‚€') is normalized cell coordinate (' meaning at previous step),
- // 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)
- // 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
-#ifndef VELGRAD
- if(cellChanged && homotheticCellResize){ for(int i=0; i<3; i++) dPos[i]=(state->pos[i]/prevCellSize[i])*scene->cell->getSize()[i]-state->pos[i]; }
-#else
- if(homotheticCellResize){ dPos = scene->cell->getShearIncrMatrix()*scene->cell->wrapShearedPt(state->pos);
- }
-#endif
blockTranslateDOFs(state->blockedDOFs, state->accel);
state->vel+=dt*state->accel;
state->pos += state->vel*dt + scene->forces.getMove(id);
- //apply cell deformation
- if(homotheticCellResize>=1) state->pos+=dPos;
- //update velocity for usage in rate dependant equations (e.g. viscous law) FIXME : it is not recommended to do that because it impacts the dynamics (this modified velocity will be used as reference in the next time-step)
- if(homotheticCellResize==2) state->vel+=dPos/dt;
+ assert(homotheticCellResize>=0 && homotheticCellResize<=2);
+ if(homotheticCellResize>0){
+ Vector3r dPos(scene->cell->getStrainIncrMatrix()*scene->cell->wrapShearedPt(state->pos));
+ // apply cell deformation
+ if(homotheticCellResize>=1) state->pos+=dPos;
+ // update velocity for usage in rate dependant equations (e.g. viscous law)
+ // FIXME : it is not recommended to do that because it impacts the dynamics (this modified velocity will be used as reference in the next time-step)
+ if(homotheticCellResize==2) state->vel+=dPos/dt;
+ }
}
inline void NewtonIntegrator::leapfrogSphericalRotate(Scene* scene, State* state, const body_id_t& id, const Real& dt )
=== modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2009-12-23 12:21:15 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2009-12-25 21:58:25 +0000
@@ -74,12 +74,9 @@
}
}
TRVAR4(cellGrow,sigma,sigmaGoal,avgStiffness);
-#ifdef VELGRAD
- for(int axis=0; axis<3; axis++){
- if (scene->dt!=0) scene->cell->velGrad[axis][axis]=cellGrow[axis]/(scene->dt*scene->cell->refSize[axis]);
- }
-#endif
- scene->cell->refSize+=cellGrow;
+ assert(scene->dt>0);
+ for(int axis=0; axis<3; axis++){ scene->cell->velGrad[axis][axis]=cellGrow[axis]/(scene->dt*scene->cell->refSize[axis]); }
+ // scene->cell->refSize+=cellGrow;
// handle state transitions
if(allStressesOK){
@@ -100,7 +97,7 @@
void PeriTriaxController::strainStressStiffUpdate(){
// update strain first
const Vector3r& cellSize(scene->cell->getSize());
- for(int i=0; i<3; i++) strain[i]=(cellSize[i]-refSize[i])/refSize[i];
+ for(int i=0; i<3; i++) strain[i]=scene->cell->strain[i][i];
// stress and stiffness
Vector3r sumForce(Vector3r::ZERO), sumStiff(Vector3r::ZERO), sumLength(Vector3r::ZERO);
int n=0;
@@ -128,12 +125,12 @@
CREATE_LOGGER(PeriTriaxController);
-void PeriTriaxController::action(Scene* scene){
+void PeriTriaxController::action(Scene*){
if(!scene->isPeriodic){ throw runtime_error("PeriTriaxController run on aperiodic simulation."); }
const Vector3r& cellSize=scene->cell->getSize();
Vector3r cellArea=Vector3r(cellSize[1]*cellSize[2],cellSize[0]*cellSize[2],cellSize[0]*cellSize[1]);
// initial updates
- if(refSize[0]<0) refSize=cellSize;
+ const Vector3r& refSize=scene->cell->refSize;
if(maxBodySpan[0]<=0){
FOREACH(const shared_ptr<Body>& b,*scene->bodies){
if(!b || !b->bound) continue;
@@ -182,18 +179,21 @@
}
}
}
+ assert(scene->dt>0.);
// update stress and strain
for(int axis=0; axis<3; axis++){
-#ifdef VELGRAD
- if (scene->dt!=0) scene->cell->velGrad[axis][axis]=cellGrow[axis]/(scene->dt*refSize[axis]);
-#endif
- strain[axis]+=cellGrow[axis]/refSize[axis];
+ // either prescribe velocity gradient
+ scene->cell->velGrad[axis][axis]=cellGrow[axis]/(scene->dt*refSize[axis]);
+ // or strain increment (but NOT both)
+ // strain[axis]+=cellGrow[axis]/refSize[axis];
+
// take in account something like poisson's effect here…
//Real bogusPoisson=0.25; int ax1=(axis+1)%3,ax2=(axis+2)%3;
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->cell->refSize+=cellGrow;
+ // scene->cell->refSize+=cellGrow;
+
strainRate=cellGrow/scene->dt;
if(allOk){
=== modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp 2009-12-09 17:11:51 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp 2009-12-25 21:58:25 +0000
@@ -65,8 +65,6 @@
int globUpdate;
//! python command to be run when the desired state is reached
string doneHook;
- //! reference cell size (set automatically; if refSize[0]<0 (initial default), the current size is referenced)
- Vector3r refSize;
//! maximum body dimension (set automatically)
Vector3r maxBodySpan;
@@ -86,8 +84,8 @@
void action(Scene*);
void strainStressStiffUpdate();
- PeriTriaxController(): reversedForces(false),goal(Vector3r::ZERO),stressMask(0),maxStrainRate(Vector3r(1,1,1)),maxUnbalanced(1e-4),absStressTol(1e3),relStressTol(3e-5),growDamping(.25),globUpdate(5),refSize(Vector3r(-1,-1,-1)),maxBodySpan(Vector3r(-1,-1,-1)),stress(Vector3r::ZERO),strain(Vector3r::ZERO),strainRate(Vector3r::ZERO),stiff(Vector3r::ZERO),currUnbalanced(-1),prevGrow(Vector3r::ZERO){}
- REGISTER_ATTRIBUTES(GlobalEngine,(reversedForces)(goal)(stressMask)(maxStrainRate)(maxUnbalanced)(absStressTol)(relStressTol)(growDamping)(globUpdate)(doneHook)(refSize)(stress)(strain)(strainRate)(stiff));
+ PeriTriaxController(): reversedForces(false),goal(Vector3r::ZERO),stressMask(0),maxStrainRate(Vector3r(1,1,1)),maxUnbalanced(1e-4),absStressTol(1e3),relStressTol(3e-5),growDamping(.25),globUpdate(5),maxBodySpan(Vector3r(-1,-1,-1)),stress(Vector3r::ZERO),strain(Vector3r::ZERO),strainRate(Vector3r::ZERO),stiff(Vector3r::ZERO),currUnbalanced(-1),prevGrow(Vector3r::ZERO){}
+ REGISTER_ATTRIBUTES(GlobalEngine,(reversedForces)(goal)(stressMask)(maxStrainRate)(maxUnbalanced)(absStressTol)(relStressTol)(growDamping)(globUpdate)(doneHook)(stress)(strain)(strainRate)(stiff));
DECLARE_LOGGER;
REGISTER_CLASS_AND_BASE(PeriTriaxController,GlobalEngine);
};
=== modified file 'pkg/dem/meta/Shop.cpp'
--- pkg/dem/meta/Shop.cpp 2009-12-25 14:46:48 +0000
+++ pkg/dem/meta/Shop.cpp 2009-12-25 21:58:25 +0000
@@ -70,15 +70,22 @@
/*! Flip periodic cell by given number of cells.
+Still broken, some interactions are missed. Should be checked.
*/
Matrix3r Shop::flipCell(const Matrix3r& _flip){
- Scene* scene=Omega::instance().getScene().get(); const shared_ptr<Cell>& cell(scene->cell);
+ Scene* scene=Omega::instance().getScene().get(); const shared_ptr<Cell>& cell(scene->cell); const Matrix3r& strain(cell->strain);
+ Vector3r size=cell->getSize();
Matrix3r flip;
if(_flip==Matrix3r::ZERO){
- LOG_ERROR("Computing optimal cell flip not yet implemented, no flipping done!");
- return Matrix3r::ZERO;
- // flip=optimal matrix …
+ bool hasNonzero=false;
+ for(int i=0; i<3; i++) for(int j=0; j<3; j++) {
+ if(i==j){ flip[i][j]=0; continue; }
+ flip[i][j]=-floor(.5+strain[i][j]/(size[j]/size[i]));
+ if(flip[i][j]!=0) hasNonzero=true;
+ }
+ if(!hasNonzero) {LOG_TRACE("No flip necessary."); return Matrix3r::ZERO;}
+ LOG_DEBUG("Computed flip matrix: upper "<<flip[0][1]<<","<<flip[0][2]<<","<<flip[1][2]<<"; lower "<<flip[1][0]<<","<<flip[2][0]<<","<<flip[2][1]);
} else {
flip=_flip;
}
@@ -90,7 +97,6 @@
}
// change cell strain here
- Vector3r size=cell->getSize();
Matrix3r strainInc;
for(int i=0; i<3; i++) for(int j=0; j<3; j++){
if(i==j) { if(flip[i][j]!=0) LOG_WARN("Non-zero diagonal term at ["<<i<<","<<j<<"] is meaningless and will be ignored."); strainInc[i][j]=0; continue; }
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2009-12-25 14:46:48 +0000
+++ py/_utils.cpp 2009-12-25 21:58:25 +0000
@@ -456,7 +456,7 @@
def("wireAll",wireAll);
def("wireNone",wireNone);
def("wireNoSpheres",wireNoSpheres);
- def("flipCell",&Shop::flipCell);
+ def("flipCell",&Shop::flipCell,(python::arg("flip")=Matrix3r::ZERO));
}
=== modified file 'py/yadeWrapper/yadeWrapper.cpp'
--- py/yadeWrapper/yadeWrapper.cpp 2009-12-25 15:28:07 +0000
+++ py/yadeWrapper/yadeWrapper.cpp 2009-12-25 21:58:25 +0000
@@ -766,10 +766,8 @@
python::class_<Cell,shared_ptr<Cell>, python::bases<Serializable>, noncopyable>("Cell",python::no_init)
.def_readwrite("refSize",&Cell::refSize)
.def_readwrite("strain",&Cell::strain)
-#ifdef VELGRAD
.def_readwrite("velGrad",&Cell::velGrad)
- .def_readwrite("Hsize",&Cell::Hsize)
-#endif
+ //.def_readwrite("Hsize",&Cell::Hsize)
.add_property("extension",&Cell::getExtensionalStrain)
//.add_property("size",&Cell::getSize,python::return_value_policy<python::return_internal_referece>()
;
=== modified file 'scripts/test/periodic-compress.py'
--- scripts/test/periodic-compress.py 2009-12-25 14:46:48 +0000
+++ scripts/test/periodic-compress.py 2009-12-25 21:58:25 +0000
@@ -7,7 +7,7 @@
for sph in p:
O.bodies.append(utils.sphere(sph[0],sph[1]))
-#log.setLevel("PeriIsoCompressor",log.TRACE)
+log.setLevel("PeriIsoCompressor",log.DEBUG)
O.timingEnabled=True
O.engines=[
ForceResetter(),
@@ -19,7 +19,7 @@
[Law2_Dem3Dof_Elastic_Elastic()],
),
PeriIsoCompressor(charLen=.5,stresses=[-50e9,-1e8],doneHook="print 'FINISHED'; O.pause() ",keepProportions=True),
- NewtonIntegrator(damping=.4)
+ NewtonIntegrator(damping=.4,homotheticCellResize=1)
]
O.dt=utils.PWaveTimeStep()
O.saveTmp()
@@ -28,6 +28,9 @@
O.run()
O.wait()
timing.stats()
+#while True:
+# O.step()
+
# now take that packing and pad some larger volume with it
#sp=pack.SpherePack()
=== modified file 'scripts/test/periodic-grow.py'
--- scripts/test/periodic-grow.py 2009-12-25 14:46:48 +0000
+++ scripts/test/periodic-grow.py 2009-12-25 21:58:25 +0000
@@ -12,7 +12,7 @@
[SimpleElasticRelationships()],
[Law2_Dem3Dof_Elastic_Elastic()],
),
- NewtonIntegrator(damping=.6)
+ NewtonIntegrator(damping=.6,homotheticCellResize=1)
]
import random
for i in xrange(250):
@@ -25,14 +25,13 @@
O.saveTmp()
from yade import qt
qt.Controller(); qt.View()
-step=.01
O.run(200,True)
-for i in range(0,250):
- O.run(200,True)
- O.cell.refSize=O.cell.refSize*.998
- if (i%10==0):
- F,stiff=utils.totalForceInVolume()
- dim=O.cell.refSize; 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()
+rate=-1e-4*cubeSize/(O.dt*200)*Matrix3().IDENTITY
+O.cell['velGrad']=rate
+for i in range(0,25):
+ O.run(2000,True)
+ F,stiff=utils.totalForceInVolume()
+ dim=O.cell.refSize; 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-triax-velgrad.py'
--- scripts/test/periodic-triax-velgrad.py 2009-12-25 14:46:48 +0000
+++ scripts/test/periodic-triax-velgrad.py 2009-12-25 21:58:25 +0000
@@ -4,6 +4,7 @@
from yade import *
from yade import pack,log,qt
#log.setLevel('PeriTriaxController',log.DEBUG)
+#log.setLevel('Shop',log.TRACE)
O.periodic=True
O.cell.refSize=Vector3(.1,.1,.1)
O.cell.Hsize=Matrix3(0.1,0,0, 0,0.1,0, 0,0,0.1)
@@ -24,13 +25,14 @@
),
PeriTriaxController(goal=[-1e5,-1e5,0],stressMask=3,globUpdate=5,maxStrainRate=[1.,1.,1.],doneHook='triaxDone()',label='triax'),
NewtonIntegrator(damping=.6, homotheticCellResize=1),
+ PeriodicPythonRunner(command='utils.flipCell()',iterPeriod=1000), # broken for larger strains?
]
-O.dt=0.1*utils.PWaveTimeStep()
+O.dt=0.5*utils.PWaveTimeStep()
O.run(1)
qt.View()
O.cell.velGrad=Matrix3(0,5,0,0,0,0, 0,0,-5)
O.run();
-
+rrr=qt.Renderer(); rrr['intrAllWire'],rrr['Body_interacting_geom']=True,False
phase=0
def triaxDone():
