------------------------------------------------------------ revno: 1925 committer: Václav Šmilauer <[email protected]> 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 <[email protected]> // #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():
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
