yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02851
[Branch ~yade-dev/yade/trunk] Rev 1926: 1. Use Cell::trsf instead of Cell::strain (which really meant trsf-Id)
------------------------------------------------------------
revno: 1926
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
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
