------------------------------------------------------------ revno: 1926 committer: Václav Šmilauer <[email protected]> branch nick: trunk timestamp: Sat 2009-12-26 22:57:37 +0100 message: 1. Use Cell::trsf instead of Cell::strain (which really meant trsf-Id) 2. NewtonIntegrator now precomputes cellTrsfInc for the current step, instead of using Cell::getTrsfInc getting increment from the previous step (!?) 3. Update scripts and other code. modified: core/Cell.hpp pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp pkg/dem/meta/Shop.cpp py/yadeWrapper/yadeWrapper.cpp scripts/test/periodic-shear.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 21:58:25 +0000 +++ core/Cell.hpp 2009-12-26 21:57:37 +0000 @@ -3,18 +3,20 @@ /*! Periodic cell parameters and routines. Usually instantiated as Scene::cell. -The cell has size and shear, which are independent. +The Cell has reference box configuration (*refSize*) which is transformed (using *trsf*) to the current parallelepiped configuration. Using notation from http://en.wikipedia.org/wiki/Finite_strain_theory: + +* Vector3r *refSize* is undeformed cell size (box) +* Matrix3r *trsf* is "deformation gradient tensor" F (written as matrix) +* Matrix3r *velGrad* is ⦠*/ class Cell: public Serializable{ public: - Cell(): refSize(Vector3r(1,1,1)), strain(Matrix3r::ZERO), velGrad(Matrix3r::ZERO){ integrateAndUpdate(0); } + Cell(): refSize(Vector3r(1,1,1)), trsf(Matrix3r::IDENTITY), velGrad(Matrix3r::ZERO){ integrateAndUpdate(0); } - //! size of the cell, projected on axes (for non-zero shear) + //! size of the reference cell Vector3r refSize; - //! Strain matrix (current total strain) - //! Ordering see http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/strain.cfm#matrix - Matrix3r strain; + Matrix3r trsf; Matrix3r velGrad; //Matrix3r Hsize; @@ -22,16 +24,16 @@ const Vector3r& getSize() const { return _size; } //! Return copy of the current size (used only by the python wrapper) Vector3r getSize_copy() const { return _size; } - //! return normal strain component (strain matrix diagonal as Vector3r) - Vector3r getExtensionalStrain() const { return Vector3r(strain[0][0],strain[1][1],strain[2][2]); } + //! stretch ratio Î(n) (http://en.wikipedia.org/wiki/Finite_strain_theory#Stretch_ratio) + Vector3r getStretchRatio() const { return Vector3r(trsf[0][0],trsf[1][1],trsf[2][2]); } //! return matrix containing cosines of shear angles const Matrix3r& getCosMatrix() const {return _cosMatrix;} //! transformation matrix applying puse shear (normal strain removed: ones on the diagonal) - const Matrix3r& getShearTrsfMatrix() const { return _shearTrsfMatrix; } + const Matrix3r& getShearTrsf() const { return _shearTrsf; } //! inverse of getShearTrsfMatrix(). - const Matrix3r& getUnshearTrsfMatrix() const {return _unshearTrsfMatrix;} - //! transformation increment matrix applying arbitrary field - const Matrix3r& getStrainIncrMatrix() const { return _strainInc; } + const Matrix3r& getUnshearTrsf() const {return _unshearTrsf;} + //! transformation increment matrix applying arbitrary field (remove if not used in NewtonIntegrator!) + // const Matrix3r& getTrsfInc() const { return _trsfInc; } /*! return pointer to column-major OpenGL 4x4 matrix applying pure shear. This matrix is suitable as argument for glMultMatrixd. @@ -51,16 +53,16 @@ // caches; private private: - Matrix3r _strainInc; + Matrix3r _trsfInc; Vector3r _size; bool _hasShear; - Matrix3r _shearTrsfMatrix, _unshearTrsfMatrix, _cosMatrix; + Matrix3r _shearTrsf, _unshearTrsf, _cosMatrix; double _glShearTrsfMatrix[16]; void fillGlShearTrsfMatrix(double m[16]){ - m[0]=1; m[4]=strain[0][1]; m[8]=strain[0][2]; m[12]=0; - m[1]=strain[1][0]; m[5]=1; m[9]=strain[1][2]; m[13]=0; - m[2]=strain[2][0]; m[6]=strain[2][1]; m[10]=1; m[14]=0; - m[3]=0; m[7]=0; m[11]=0; m[15]=1; + m[0]=1; m[4]=trsf[0][1]; m[8]=trsf[0][2]; m[12]=0; + m[1]=trsf[1][0]; m[5]=1; m[9]=trsf[1][2]; m[13]=0; + m[2]=trsf[2][0]; m[6]=trsf[2][1]; m[10]=1; m[14]=0; + m[3]=0; m[7]=0; m[11]=0; m[15]=1; } public: @@ -75,26 +77,26 @@ #endif //incremental displacement gradient - _strainInc=dt*velGrad; + _trsfInc=dt*velGrad; // total transformation; M = (Id+G).M = F.M - strain+=_strainInc*_shearTrsfMatrix; + trsf+=_trsfInc*trsf; #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 shear trsf: ones on diagonal + _shearTrsf=trsf; _shearTrsf[0][0]=_shearTrsf[1][1]=_shearTrsf[2][2]=1.; // pure unshear transformation - _unshearTrsfMatrix=_shearTrsfMatrix.Inverse(); + _unshearTrsf=_shearTrsf.Inverse(); // 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); + _hasShear=(trsf[0][1]!=0 || trsf[0][2]!=0 || trsf[1][0]!=0 || trsf[1][2]!=0 || trsf[2][0]!=0 || trsf[2][1]!=0); // current cell size (in units on physical space axes) - _size=refSize+diagMult(getExtensionalStrain(),refSize); + _size=diagMult(getStretchRatio(),refSize); // 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]))); + for(int i=0; i<3; i++) for(int j=0; j<3; j++) _cosMatrix[i][j]=(i==j?0:cos(atan(trsf[i][j]))); } /*! Return point inside periodic cell, even if shear is applied */ @@ -102,9 +104,9 @@ /*! Return point inside periodic cell, even if shear is applied; store cell coordinates in period. */ Vector3r wrapShearedPt(const Vector3r& pt, Vector3<int>& period){ return shearPt(wrapPt(unshearPt(pt),period)); } /*! Apply inverse shear on point; to put it inside (unsheared) periodic cell, apply wrapPt on the returned value. */ - Vector3r unshearPt(const Vector3r& pt){ return _unshearTrsfMatrix*pt; } + Vector3r unshearPt(const Vector3r& pt){ return _unshearTrsf*pt; } //! Apply shear on point. - Vector3r shearPt(const Vector3r& pt){ return _shearTrsfMatrix*pt; } + Vector3r shearPt(const Vector3r& pt){ return _shearTrsf*pt; } /*! Wrap point to inside the periodic cell; don't compute number of periods wrapped */ Vector3r wrapPt(const Vector3r pt)const{ Vector3r ret; for(int i=0;i<3;i++) ret[i]=wrapNum(pt[i],_size[i]); return ret; @@ -122,7 +124,7 @@ Real norm=x/sz; period=(int)floor(norm); return (norm-period)*sz; } void postProcessAttributes(bool deserializing){ if(deserializing) integrateAndUpdate(0); } - REGISTER_ATTRIBUTES(Serializable,(refSize)(strain)(velGrad)); + REGISTER_ATTRIBUTES(Serializable,(refSize)(trsf)(velGrad)); REGISTER_CLASS_AND_BASE(Cell,Serializable); }; REGISTER_SERIALIZABLE(Cell); === modified file 'pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp' --- pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp 2009-12-25 21:58:25 +0000 +++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp 2009-12-26 21:57:37 +0000 @@ -65,6 +65,8 @@ { scene->forces.sync(); Real dt=scene->dt; + // precompute transformation increment; using Cell::getTrsfInc would get increment from the previous step, which is not right... (?) + if(scene->isPeriodic && homotheticCellResize) cellTrsfInc=dt*scene->cell->velGrad; // 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 @@ -163,7 +165,8 @@ state->pos += state->vel*dt + scene->forces.getMove(id); assert(homotheticCellResize>=0 && homotheticCellResize<=2); if(homotheticCellResize>0){ - Vector3r dPos(scene->cell->getStrainIncrMatrix()*scene->cell->wrapShearedPt(state->pos)); + //Vector3r dPos(scene->cell->getTrsfInc()*scene->cell->wrapShearedPt(state->pos)); + Vector3r dPos(cellTrsfInc*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) === modified file 'pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp' --- pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp 2009-12-23 01:27:30 +0000 +++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp 2009-12-26 21:57:37 +0000 @@ -66,6 +66,8 @@ bool exactAsphericalRot; //! Enable artificially moving all bodies with the periodic cell, such that its resizes are isotropic. 0: disabled (default), 1: position update, 2: velocity update. int homotheticCellResize; + //! Store transformation increment for the current step (updated automatically) + Matrix3r cellTrsfInc; #ifdef YADE_OPENMP vector<Real> threadMaxVelocitySq; === modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp' --- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2009-12-25 21:58:25 +0000 +++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2009-12-26 21:57:37 +0000 @@ -97,7 +97,7 @@ void PeriTriaxController::strainStressStiffUpdate(){ // update strain first const Vector3r& cellSize(scene->cell->getSize()); - for(int i=0; i<3; i++) strain[i]=scene->cell->strain[i][i]; + for(int i=0; i<3; i++) strain[i]=scene->cell->trsf[i][i]-1.; // stress and stiffness Vector3r sumForce(Vector3r::ZERO), sumStiff(Vector3r::ZERO), sumLength(Vector3r::ZERO); int n=0; === modified file 'pkg/dem/meta/Shop.cpp' --- pkg/dem/meta/Shop.cpp 2009-12-25 21:58:25 +0000 +++ pkg/dem/meta/Shop.cpp 2009-12-26 21:57:37 +0000 @@ -74,14 +74,14 @@ */ Matrix3r Shop::flipCell(const Matrix3r& _flip){ - Scene* scene=Omega::instance().getScene().get(); const shared_ptr<Cell>& cell(scene->cell); const Matrix3r& strain(cell->strain); + Scene* scene=Omega::instance().getScene().get(); const shared_ptr<Cell>& cell(scene->cell); const Matrix3r& trsf(cell->trsf); Vector3r size=cell->getSize(); Matrix3r flip; if(_flip==Matrix3r::ZERO){ 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])); + flip[i][j]=-floor(.5+trsf[i][j]/(size[j]/size[i])); if(flip[i][j]!=0) hasNonzero=true; } if(!hasNonzero) {LOG_TRACE("No flip necessary."); return Matrix3r::ZERO;} @@ -96,17 +96,17 @@ if(!b) continue; cell->wrapShearedPt(b->state->pos,oldCells[b->getId()]); } - // change cell strain here - Matrix3r strainInc; + // change cell trsf here + Matrix3r trsfInc; 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; } + if(i==j) { if(flip[i][j]!=0) LOG_WARN("Non-zero diagonal term at ["<<i<<","<<j<<"] is meaningless and will be ignored."); trsfInc[i][j]=0; continue; } // make sure non-diagonal entries are "integers" if(flip[i][j]!=double(int(flip[i][j]))) LOG_WARN("Flip matrix entry "<<flip[i][j]<<" at ["<<i<<","<<j<<"] not integer?! (will be rounded)"); - strainInc[i][j]=int(flip[i][j])*size[j]/size[i]; + trsfInc[i][j]=int(flip[i][j])*size[j]/size[i]; } - TRWM3MAT(cell->strain); - TRWM3MAT(strainInc); - cell->strain+=strainInc; + TRWM3MAT(cell->trsf); + TRWM3MAT(trsfInc); + cell->trsf+=trsfInc; cell->postProcessAttributes(/*deserializing*/true); // new cell coords of bodies === modified file 'py/yadeWrapper/yadeWrapper.cpp' --- py/yadeWrapper/yadeWrapper.cpp 2009-12-25 21:58:25 +0000 +++ py/yadeWrapper/yadeWrapper.cpp 2009-12-26 21:57:37 +0000 @@ -765,10 +765,10 @@ python::class_<Cell,shared_ptr<Cell>, python::bases<Serializable>, noncopyable>("Cell",python::no_init) .def_readwrite("refSize",&Cell::refSize) - .def_readwrite("strain",&Cell::strain) + .def_readwrite("trsf",&Cell::trsf) .def_readwrite("velGrad",&Cell::velGrad) //.def_readwrite("Hsize",&Cell::Hsize) - .add_property("extension",&Cell::getExtensionalStrain) + .add_property("stretch",&Cell::getStretchRatio) //.add_property("size",&Cell::getSize,python::return_value_policy<python::return_internal_referece>() ; === modified file 'scripts/test/periodic-shear.py' --- scripts/test/periodic-shear.py 2009-12-25 14:46:48 +0000 +++ scripts/test/periodic-shear.py 2009-12-26 21:57:37 +0000 @@ -19,7 +19,7 @@ ] def doCellFlip(): - flip=1 if O.cell.strain[1,2]<0 else -1; + flip=1 if O.cell.trsf[1,2]<0 else -1; utils.flipCell(Matrix3(0,0,0, 0,0,flip, 0,0,0)) #g=0. @@ -27,7 +27,7 @@ # 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 -O.cell.strain=Matrix3(0,0,0, 0,0,.5, 0,0,0) +O.cell.trsf=Matrix3(1,0,0, 0,1,.5, 0,0,1) O.dt=2e-2*utils.PWaveTimeStep() O.step() O.saveTmp() === modified file 'scripts/test/periodic-triax-velgrad.py' --- scripts/test/periodic-triax-velgrad.py 2009-12-25 21:58:25 +0000 +++ scripts/test/periodic-triax-velgrad.py 2009-12-26 21:57:37 +0000 @@ -7,7 +7,9 @@ #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) +#O.cell.Hsize=Matrix3(0.1,0,0, 0,0.1,0, 0,0,0.1) +O.cell.trsf=Matrix3().IDENTITY; + sp=pack.SpherePack() radius=5e-3 num=sp.makeCloud(Vector3().ZERO,O.cell.refSize,radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
