yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #04615
[Branch ~yade-dev/yade/trunk] Rev 2257: 1. Add Peri3dController, done by Jan and me during a rainy day
------------------------------------------------------------
revno: 2257
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Wed 2010-05-26 19:49:15 +0200
message:
1. Add Peri3dController, done by Jan and me during a rainy day
2. Implement memoizeDb for pack.randomPeriPack
3. add polar decomposition to Matrix3 in python
4. make the 3d view update even if there are no bodies
added:
scripts/test/peri3d.py
modified:
doc/references.bib
gui/qt3/SimulationController.cpp
pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp
pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp
pkg/dem/meta/ConcretePM.cpp
py/mathWrap/miniEigen.cpp
py/pack/pack.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 'doc/references.bib'
--- doc/references.bib 2010-05-07 11:15:06 +0000
+++ doc/references.bib 2010-05-26 17:49:15 +0000
@@ -105,4 +105,17 @@
publisher = {American Physical Society}
}
+@inproceedings{Kuhl2001,
+ title = {Microplane modelling and particle modelling of cohesive-frictional materials},
+ author = {E. Kuhl and G. A. D'Addetta and M. Leukart and E. Ramm},
+ booktitle = {Continuous and Discontinuous Modelling of Cohesive-Frictional Materials},
+ publisher = {Springer Berlin / Heidelberg},
+ issn = {1616-6361},
+ isbn = {978-3-540-41525-1},
+ volume = {568},
+ series = {Lecture Notes in Physics},
+ year = {2001},
+ pages = {31--46},
+ doi = {10.1007/3-540-44424-6_3},
+}
=== modified file 'gui/qt3/SimulationController.cpp'
--- gui/qt3/SimulationController.cpp 2010-05-24 21:14:57 +0000
+++ gui/qt3/SimulationController.cpp 2010-05-26 17:49:15 +0000
@@ -430,7 +430,7 @@
if(sbRefreshTime->value()!=refreshTime) sbRefreshTime->setValue(refreshTime);
/* enable/disable controls here, dynamically */
- hasSimulation=(scene ? scene->bodies->size()>0 : false );
+ hasSimulation=(scene); // ? scene->bodies->size()>0 : false );
bool isRunning=Omega::instance().isRunning() || syncRunning,
hasTimeStepper=scene->timeStepperPresent(),
usesTimeStepper=scene->timeStepperActive(),
=== modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2010-05-24 19:56:22 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2010-05-26 17:49:15 +0000
@@ -10,7 +10,7 @@
using namespace std;
-YADE_PLUGIN((PeriIsoCompressor)(PeriTriaxController))
+YADE_PLUGIN((PeriIsoCompressor)(PeriTriaxController)(Peri3dController))
CREATE_LOGGER(PeriIsoCompressor);
@@ -245,3 +245,128 @@
}
}
}
+
+#ifndef YADE_WM3
+
+CREATE_LOGGER(Peri3dController);
+void Peri3dController::update(){
+ /* strain
+ polar decomposition of transformation:
+ compute strain tensor from the non-rotational part of trsf, then rotate it back to global coords
+ */
+ Matrix3r rot,nonrot; //nonrot=skew+normal deformation
+ const Matrix3r& trsf=scene->cell->trsf;
+ Eigen::SVD<Matrix3r>(trsf).computeUnitaryPositive(&rot,&nonrot);
+ // FIXME: should be rot*(log(nonrot)), but matrix logarithm is not yet implemented in eigen
+ strain=rot*(nonrot-Matrix3r::Identity());
+ LOG_TRACE("Updated strain value\n"<<strain);
+ /* stress and stiffness
+ */
+ K=Matrix6r::Zero();
+ stress=Matrix3r::Zero();
+ const Real volume=scene->cell->trsf.determinant()*scene->cell->refSize[0]*scene->cell->refSize[1]*scene->cell->refSize[2];
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(I->interactionGeometry.get());
+ NormShearPhys* phys=YADE_CAST<NormShearPhys*>(I->interactionPhysics.get());
+ // not clear whether this should be the reference or the current distance
+ // current: gives consistent results for same configuration with different initial state
+ // reference: does not change stress tensor when the same forces exist on interactions with changing length
+ //const Real& d0=geom->refLength;
+ const Real d0=(geom->se31.position-geom->se32.position).norm();
+ const Vector3r& n=geom->normal;
+ const Vector3r& fT=phys->shearForce;
+ const Real fN=phys->normalForce.dot(n);
+ for(int i=0; i<3; i++) for(int j=0;j<3; j++){
+ stress(i,j)+=d0*(fN*n[i]*n[j]+.5*(fT[i]*n[j]+fT[j]*n[i]));
+ }
+ const Real& kN=phys->kn; const Real& kT=phys->ks;
+ // mapping between 6x6 matrix indices and tensor indices (Voigt notation)
+ const int map[6][6][4]={
+ {{0,0,0,0},{0,0,1,1},{0,0,2,2},{0,0,1,2},{0,0,2,0},{0,0,0,1}},
+ {{8,8,8,8},{1,1,1,1},{1,1,2,2},{1,1,1,2},{1,1,2,0},{1,1,0,1}},
+ {{8,8,8,8},{8,8,8,8},{2,2,2,2},{2,2,1,2},{2,2,2,0},{2,2,0,1}},
+ {{8,8,8,8},{8,8,8,8},{8,8,8,8},{1,2,1,2},{1,2,2,0},{1,2,0,1}},
+ {{8,8,8,8},{8,8,8,8},{8,8,8,8},{8,8,8,8},{2,0,2,0},{2,0,0,1}},
+ {{8,8,8,8},{8,8,8,8},{8,8,8,8},{8,8,8,8},{8,8,8,8},{0,1,0,1}}};
+ const int kron[3][3]={{1,0,0},{0,1,0},{0,0,1}}; // Kronecker delta
+ for(int p=0; p<6; p++) for(int q=p;q<6;q++){
+ int i=map[p][q][0], j=map[p][q][1], k=map[p][q][2], l=map[p][q][3];
+ K(p,q)+=d0*d0*(kN*n[i]*n[j]*n[k]*n[l]+kT*(.25*(n[j]*n[k]*kron[i][l]+n[j]*n[l]*kron[i][k]+n[i]*n[k]*kron[j][l]+n[i]*n[l]*kron[j][k])-n[i]*n[j]*n[k]*n[l]));
+ }
+ }
+ stress/=volume;
+ K/=volume;
+ for(int p=0; p<6; p++)for(int q=p+1;q<6;q++) K(q,p)=K(p,q);
+ LOG_TRACE("Updated stress\n"<<stress);
+ LOG_TRACE("Updated stiffness tensor\n"<<K);
+}
+
+void Peri3dController::action(){
+ // TODO: only call sometimes
+ update();
+ typedef Eigen::Matrix<Real,Eigen::Dynamic,Eigen::Dynamic> MatrixXr;
+ typedef Eigen::Matrix<Real,Eigen::Dynamic,1> VectorXr;
+ /*
+ sigma = K * eps
+
+ decompose the stiffness matrix depending on what is prescribed
+
+ | sigma_u | = | Kuu Kup | * | eps_u |
+ | sigma_p | | Kpu Kpp | | eps_p |
+
+ then
+
+ eps_u=(Kuu^-1)*(sigma_u-Kup eps_p)
+
+ (we replace all eps and sigma by their rates in the computation)
+ */
+ int numP=0,numU=0; // number of prescribed and unprescribed strain components
+ for(int i=0;i<6;i++) if(stressMask&(1<<i))numU++;
+ numP=6-numU;
+ // sub-matrices of the stiffness matrix
+ MatrixXr Kuu(numU,numU), Kup(numU,numP);
+ // epsP, epsU, sigU are _rates_
+ VectorXr epsP(numP), epsU(numU), sigU(numU);
+ // conversion from Voigt indices to 3x3 tensor indices (upper-triangle part)
+ const int mapI[6]={0,1,2,1,2,0};
+ const int mapJ[6]={0,1,2,2,0,1};
+ int jU=0,jP=0;
+ for(int j=0; j<6; j++){
+ // prescribed stress, i.e. un-prescribed strain
+ if(stressMask&(1<<j)){
+ sigU[jU]=goal(mapI[j],mapJ[j]);
+ int iU=0;
+ for(int i=0;i<6;i++){ if((stressMask&(1<<i))) Kuu(iU++,jU)=K(i,j);}
+ jU++;
+ } else {
+ epsP[jP]=(j<3?1.:2.)*goal(mapI[j],mapJ[j]); /* multiply tensor shear by 2, see http://en.wikiversity.org/wiki/Introduction_to_Elasticity/Constitutive_relations */
+ int iP=0;
+ for(int i=0;i<6;i++){ if((stressMask&(1<<i))) Kup(iP++,jP)=K(i,j);}
+ jP++;
+ }
+ }
+ assert(jU==numU); assert(jP==numP);
+ LOG_TRACE("Kuu=\n"<<Kuu<<"\nKup=\n"<<Kup<<"\nepsP=\n"<<epsP<<"\nsigU=\n"<<sigU);
+ // if Kuu is (nearly) singular, the goal strain is inf and the corresponding velGrad component then NaN
+ // in such case, sanitize it with some random value (irrelevant which one)
+ // FIXME: find a better way for this than determinant (expensive for larger matrices)
+ if(Kuu.rows()*Kuu.cols()>0 && Kuu.determinant()<1e-20){ Kuu+=MatrixXr(Kuu.cols(),Kuu.rows()).setIdentity(); LOG_TRACE("Kuu after sanitization\n"<<Kuu); }
+ epsU=Kuu.inverse()*(sigU-Kup*epsP);
+ // assemble strain tensor (as matrix), from prescribed values epsP and computed values epsU
+ jU=0; jP=0;
+ Matrix3r eps;
+ for(int j=0; j<6; j++){
+ if(stressMask&(1<<j)) eps(mapI[j],mapJ[j])=epsU[jU++];
+ else eps(mapI[j],mapJ[j])=(j<3?1.:.5)*epsP[jP++]; /* multiply shear components back by 1/2 when converting from Voigt vector back to tensor */
+ }
+ eps(2,1)=eps(1,2); eps(0,2)=eps(2,0); eps(1,0)=eps(0,1);
+ Matrix3r& velGrad=scene->cell->velGrad;
+ // rate of (goal strain - current strain)
+ velGrad=(eps-strain)/scene->dt;
+ Real mx=max(abs(velGrad.minCoeff()),abs(velGrad.maxCoeff()));
+ if(mx>abs(maxStrainRate)) velGrad*=abs(maxStrainRate)/mx;
+ LOG_TRACE("epsU=\n"<<epsU<<"\neps=\n"<<"\nabs(maxCoeff)="<<mx<<"\nvelGrad=\n"<<velGrad);
+ // TODO: check unbalanced force and run some hook when the goal state is achieved
+}
+#endif
=== modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp 2010-04-25 13:18:11 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.hpp 2010-05-26 17:49:15 +0000
@@ -70,3 +70,31 @@
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(PeriTriaxController);
+#ifndef YADE_WM3
+
+#include<Eigen/SVD>
+class Peri3dController: public BoundaryController{
+ public:
+ typedef Eigen::Matrix<Real,6,6> Matrix6r;
+ typedef Eigen::Matrix<Real,6,1> Vector6r;
+
+ // stiffness matrix, updated automatically
+ Matrix6r K;
+
+ virtual void action();
+ void update();
+ YADE_CLASS_BASE_DOC_ATTRS(Peri3dController,BoundaryController,"Experimental controller of full strain/stress tensors on periodic cell. Stress and strain tensors are computed using formulas derived in [Kuhl2001]_, in particular equations (33) and (35).",
+ ((Matrix3r,strain,Matrix3r::Zero(),"Current deformation tensor |yupdate|"))
+ ((Matrix3r,stress,Matrix3r::Zero(),"Current stress tensor |yupdate|"))
+ ((Matrix3r,goal,Matrix3r::Zero(),"Goal state."))
+ ((int,stressMask,((void)"all strains",0),"mask determining whether components of :yref:`goal<Peri3dController.goal>` are strain (0) or stress (1). The order is 00,11,22,12,02,01 from the least significant bit. (e.g. 0b000011 is stress 00 and stress 11)."))
+ ((Real,maxStrainRate,1,"Maximum absolute value of strain rate (both normal and shear components of :yref:`Cell.velGrad`)"))
+ // not yet used
+ //((Real,currUnbalanced,NaN,"current unbalanced force |yupdate|"))
+ //((Real,maxUnbalanced,1e-4,"Maximum unbalanced force"))
+ );
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Peri3dController);
+
+#endif
=== modified file 'pkg/dem/meta/ConcretePM.cpp'
--- pkg/dem/meta/ConcretePM.cpp 2010-05-13 20:19:38 +0000
+++ pkg/dem/meta/ConcretePM.cpp 2010-05-26 17:49:15 +0000
@@ -133,12 +133,7 @@
#ifdef YADE_CPM_FULL_MODEL_AVAILABLE
- // assure interoperability with older version
- #define SquaredLength squaredNorm
- #define Length norm
- #include"../../../../brefcom-mm.hh"
- #undef SquaredLength
- #undef Length
+ #include"../../../../brefcom-mm.hh"
#endif
=== modified file 'py/mathWrap/miniEigen.cpp'
--- py/mathWrap/miniEigen.cpp 2010-05-11 10:25:49 +0000
+++ py/mathWrap/miniEigen.cpp 2010-05-26 17:49:15 +0000
@@ -112,6 +112,10 @@
static Vector3i Vector3i_cross(const Vector3i& self, const Vector3i& v){ return self.cross(v); }
static bool Quaternionr__eq__(const Quaternionr& q1, const Quaternionr& q2){ return q1==q2; }
static bool Quaternionr__neq__(const Quaternionr& q1, const Quaternionr& q2){ return q1!=q2; }
+#ifndef YADE_WM3
+ #include<Eigen/SVD>
+ static bp::tuple Matrix3r_polarDecomposition(const Matrix3r& self){ Matrix3r unitary,positive; Eigen::SVD<Matrix3r>(self).computeUnitaryPositive(&unitary,&positive); return bp::make_tuple(unitary,positive); }
+#endif
#define WM3_COMPAT
@@ -230,6 +234,10 @@
.def("determinant",&Matrix3r::determinant)
.def("inverse",&Matrix3r_inverse)
.def("transpose",&Matrix3r_transpose)
+ #ifndef YADE_WM3
+ .def("polarDecomposition",&Matrix3r_polarDecomposition)
+ #endif
+
//
.def("__neg__",&Matrix3r__neg__)
.def("__add__",&Matrix3r__add__Matrix3r).def("__iadd__",&Matrix3r__iadd__Matrix3r)
=== modified file 'py/pack/pack.py'
--- py/pack/pack.py 2010-05-14 10:01:37 +0000
+++ py/pack/pack.py 2010-05-26 17:49:15 +0000
@@ -221,6 +221,71 @@
if predicate(s[0],s[1]): ret+=[utils.sphere(s[0],radius=s[1],**kw)]
return ret
+def _memoizePacking(memoizeDb,sp,radius,rRelFuzz,wantPeri,fullDim):
+ if not memoizeDb: return
+ import cPickle,sqlite3,time,os
+ if os.path.exists(memoizeDb):
+ conn=sqlite3.connect(memoizeDb)
+ else:
+ conn=sqlite3.connect(memoizeDb)
+ c=conn.cursor()
+ c.execute('create table packings (radius real, rRelFuzz real, dimx real, dimy real, dimz real, N integer, timestamp real, periodic integer, pack blob)')
+ c=conn.cursor()
+ packBlob=buffer(cPickle.dumps(sp.toList_pointsAsTuples(),cPickle.HIGHEST_PROTOCOL))
+ packDim=sp.cellSize if wantPeri else fullDim
+ c.execute('insert into packings values (?,?,?,?,?,?,?,?,?)',(radius,rRelFuzz,packDim[0],packDim[1],packDim[2],len(sp),time.time(),wantPeri,packBlob,))
+ c.close()
+ conn.commit()
+ print "Packing saved to the database",memoizeDb
+
+def _getMemoizedPacking(memoizeDb,radius,rRelFuzz,x1,y1,z1,fullDim,wantPeri,fillPeriodic,spheresInCell,memoDbg=False):
+ """Return suitable SpherePack read from *memoizeDb* if found, None otherwise.
+
+ :param fillPeriodic: whether to fill fullDim by repeating periodic packing
+ :param wantPeri: only consider periodic packings
+ """
+ import os,os.path,sqlite3,time,cPickle,sys
+ if memoDbg:
+ def memoDbgMsg(s): print s
+ else:
+ def memoDbgMsg(s): pass
+ if not memoizeDb or not os.path.exists(memoizeDb):
+ if memoizeDb: memoDbgMsg("Database %s does not exist."%memoizeDb)
+ return None
+ # find suitable packing and return it directly
+ conn=sqlite3.connect(memoizeDb); c=conn.cursor();
+ try:
+ c.execute('select radius,rRelFuzz,dimx,dimy,dimz,N,timestamp,periodic from packings order by N')
+ except sqlite3.OperationalError:
+ raise RuntimeError("ERROR: database `"+memoizeDb+"' not compatible with randomDensePack (corrupt, deprecated format or not a db created by randomDensePack)")
+ for row in c:
+ R,rDev,X,Y,Z,NN,timestamp,isPeri=row[0:8]; scale=radius/R
+ rDev*=scale; X*=scale; Y*=scale; Z*=scale
+ memoDbgMsg("Considering packing (radius=%g±%g,N=%g,dim=%g×%g×%g,%s,scale=%g), created %s"%(R,.5*rDev,NN,X,Y,Z,"periodic" if isPeri else "non-periodic",scale,time.asctime(time.gmtime(timestamp))))
+ if not isPeri and wantPeri: memoDbgMsg("REJECT: is not periodic, which is requested."); continue
+ if (rRelFuzz==0 and rDev!=0) or (rRelFuzz!=0 and rDev==0) or (rRelFuzz!=0 and abs((rDev-rRelFuzz)/rRelFuzz)>1e-2): memoDbgMsg("REJECT: radius fuzz differs too much (%g, %g desired)"%(rDev,rRelFuzz)); continue # radius fuzz differs too much
+ if isPeri and wantPeri:
+ if spheresInCell>NN: memoDbgMsg("REJECT: Number of spheres in the packing too small"); continue
+ if abs((y1/x1)/(Y/X)-1)>0.3 or abs((z1/x1)/(Z/X)-1)>0.3: memoDbgMsg("REJECT: proportions (y/x=%g, z/x=%g) differ too much from what is desired (%g, %g)."%(Y/X,Z/X,y1/x1,z1/x1)); continue
+ else:
+ if (X<fullDim[0] or Y<fullDim[1] or Z<fullDim[2]): memoDbgMsg("REJECT: not large enough"); continue # not large enough
+ memoDbgMsg("ACCEPTED");
+ print "Found suitable packing in %s (radius=%g±%g,N=%g,dim=%g×%g×%g,%s,scale=%g), created %s"%(memoizeDb,R,rDev,NN,X,Y,Z,"periodic" if isPeri else "non-periodic",scale,time.asctime(time.gmtime(timestamp)))
+ c.execute('select pack from packings where timestamp=?',(timestamp,))
+ sp=SpherePack(cPickle.loads(str(c.fetchone()[0])))
+ if isPeri and wantPeri:
+ sp.cellSize=(X,Y,Z);
+ if fillPeriodic: sp.cellFill(Vector3(fullDim[0],fullDim[1],fullDim[2]));
+ sp.scale(scale);
+ #sp.cellSize=(0,0,0) # resetting cellSize avoids warning when rotating
+ return sp
+ #if orientation: sp.rotate(*orientation.toAxisAngle())
+ #return filterSpherePack(predicate,sp,material=material)
+ #print "No suitable packing in database found, running",'PERIODIC compression' if wantPeri else 'triaxial'
+ #sys.stdout.flush()
+
+
+
def randomDensePack(predicate,radius,material=-1,dim=None,cropLayers=0,rRelFuzz=0.,spheresInCell=0,memoizeDb=None,useOBB=True,memoDbg=False,color=None):
"""Generator of random dense packing with given geometry properties, using TriaxialTest (aperiodic)
or PeriIsoCompressor (periodic). The priodicity depens on whether the spheresInCell parameter is given.
@@ -274,37 +339,13 @@
maxR=radius*(1+rRelFuzz)
x1=max(x1,8*maxR); y1=max(y1,8*maxR); z1=max(z1,8*maxR); vol1=x1*y1*z1
N100*=vol1/vol0 # volume might have been increased, increase number of spheres to keep porosity the same
- if(memoizeDb and os.path.exists(memoizeDb)):
- if memoDbg:
- def memoDbgMsg(s): print s
- else:
- def memoDbgMsg(s): pass
- # find suitable packing and return it directly
- conn=sqlite3.connect(memoizeDb); c=conn.cursor();
- try:
- c.execute('select radius,rRelFuzz,dimx,dimy,dimz,N,timestamp,periodic from packings order by N')
- except sqlite3.OperationalError:
- raise RuntimeError("ERROR: database `"+memoizeDb+"' not compatible with randomDensePack (corrupt, deprecated format or not a db created by randomDensePack)")
- for row in c:
- R,rDev,X,Y,Z,NN,timestamp,isPeri=row[0:8]; scale=radius/R
- rDev*=scale; X*=scale; Y*=scale; Z*=scale
- memoDbgMsg("Considering packing (radius=%g±%g,N=%g,dim=%g×%g×%g,%s,scale=%g), created %s"%(R,.5*rDev,NN,X,Y,Z,"periodic" if isPeri else "non-periodic",scale,time.asctime(time.gmtime(timestamp))))
- if (rRelFuzz==0 and rDev!=0) or (rRelFuzz==0 and rDev!=0) or (rRelFuzz!=0 and abs((rDev-rRelFuzz)/rRelFuzz)>1e-2): memoDbgMsg("REJECT: radius fuzz differs too much (%g, %g desired)"%(rDev,rRelFuzz)); continue # radius fuzz differs too much
- if isPeri and wantPeri:
- if spheresInCell>NN: memoDbgMsg("REJECT: Number of spheres in the packing too small"); continue
- if abs((y1/x1)/(Y/X)-1)>0.3 or abs((z1/x1)/(Z/X)-1)>0.3: memoDbgMsg("REJECT: proportions (y/x=%g, z/x=%g) differ too much from what is desired (%g, %g)."%(Y/X,Z/X,y1/x1,z1/x1)); continue
- else:
- if (X<fullDim[0] or Y<fullDim[1] or Z<fullDim[2]): memoDbgMsg("REJECT: not large enough"); continue # not large enough
- memoDbgMsg("ACCEPTED");
- print "Found suitable packing in %s (radius=%g±%g,N=%g,dim=%g×%g×%g,%s,scale=%g), created %s"%(memoizeDb,R,rDev,NN,X,Y,Z,"periodic" if isPeri else "non-periodic",scale,time.asctime(time.gmtime(timestamp)))
- c.execute('select pack from packings where timestamp=?',(timestamp,))
- sp=SpherePack(cPickle.loads(str(c.fetchone()[0])))
- if isPeri and wantPeri:
- sp.cellSize=(X,Y,Z); sp.cellFill(Vector3(fullDim[0],fullDim[1],fullDim[2])); sp.cellSize=(0,0,0) # resetting cellSize avoids warning when rotating
- sp.scale(scale);
- if orientation: sp.rotate(*orientation.toAxisAngle())
+ sp=_getMemoizedPacking(memoizeDb,radius,rRelFuzz,x1,y1,z1,fullDim,wantPeri,fillPeriodic=True,spheresInCell=spheresInCell,memoDbg=False)
+ if sp:
+ if orientation:
+ sp.cellSize=(0,0,0) # resetting cellSize avoids warning when rotating
+ sp.rotate(*orientation.toAxisAngle())
return filterSpherePack(predicate,sp,material=material)
- print "No suitable packing in database found, running",'PERIODIC compression' if wantPeri else 'triaxial'
+ else: print "No suitable packing in database found, running",'PERIODIC compression' if wantPeri else 'triaxial'
sys.stdout.flush()
O.switchScene(); O.resetThisScene() ### !!
if wantPeri:
@@ -337,27 +378,14 @@
O.run(); O.wait()
sp=SpherePack(); sp.fromSimulation()
O.switchScene() ### !!
- if(memoizeDb):
- if os.path.exists(memoizeDb):
- conn=sqlite3.connect(memoizeDb)
- else:
- conn=sqlite3.connect(memoizeDb)
- c=conn.cursor()
- c.execute('create table packings (radius real, rRelFuzz real, dimx real, dimy real, dimz real, N integer, timestamp real, periodic integer, pack blob)')
- c=conn.cursor()
- packBlob=buffer(cPickle.dumps(sp.toList_pointsAsTuples(),cPickle.HIGHEST_PROTOCOL))
- packDim=sp.cellSize if wantPeri else fullDim
- c.execute('insert into packings values (?,?,?,?,?,?,?,?,?)',(radius,rRelFuzz,packDim[0],packDim[1],packDim[2],len(sp),time.time(),wantPeri,packBlob,))
- c.close()
- conn.commit()
- print "Packing saved to the database",memoizeDb
+ _memoizePacking(memoizeDb,sp,radius,rRelFuzz,wantPeri,fullDim)
if wantPeri: sp.cellFill(Vector3(fullDim[0],fullDim[1],fullDim[2]))
if orientation:
sp.cellSize=(0,0,0); # reset periodicity to avoid warning when rotating periodic packing
sp.rotate(*orientation.toAxisAngle())
return filterSpherePack(predicate,sp,material=material,color=color)
-def randomPeriPack(radius,rRelFuzz,initSize):
+def randomPeriPack(radius,rRelFuzz,initSize,memoizeDb=None):
"""Generate periodic dense packing. EXPERIMENTAL, you at your own risk.
A cell of initSize is stuffed with as many spheres as possible (ignore the warning from SpherePack::makeCloud about
@@ -374,6 +402,8 @@
"""
from math import pi
+ sp=_getMemoizedPacking(memoizeDb,radius,rRelFuzz,initSize[0],initSize[1],initSize[2],fullDim=Vector3(0,0,0),wantPeri=True,fillPeriodic=False,spheresInCell=-1,memoDbg=True)
+ if sp: return sp
O.switchScene(); O.resetThisScene()
sp=SpherePack()
O.periodic=True
@@ -387,6 +417,7 @@
O.run(); O.wait()
ret=SpherePack()
ret.fromSimulation()
+ _memoizePacking(memoizeDb,sp,radius,rRelFuzz,wantPeri=True,fullDim=Vector3(0,0,0)) # fullDim unused
O.switchScene()
return ret
=== added file 'scripts/test/peri3d.py'
--- scripts/test/peri3d.py 1970-01-01 00:00:00 +0000
+++ scripts/test/peri3d.py 2010-05-26 17:49:15 +0000
@@ -0,0 +1,29 @@
+from yade import pack,log
+sp=pack.randomPeriPack(.5,.5,Vector3(10,10,10),memoizeDb='/tmp/packDb.sqlite')
+O.periodic=True
+O.cell.refSize=sp.cellSize
+O.bodies.append([utils.sphere(c,r) for c,r in sp])
+O.dt=utils.PWaveTimeStep()
+log.setLevel('Peri3dController',log.TRACE)
+O.engines=[
+ ForceResetter(),
+ BoundDispatcher([Bo1_Sphere_Aabb()]),
+ InsertionSortCollider(),
+ InteractionDispatchers(
+ [Ig2_Sphere_Sphere_Dem3DofGeom()],
+ [Ip2_FrictMat_FrictMat_FrictPhys()],
+ [Law2_Dem3DofGeom_FrictPhys_Basic()]
+ ),
+ NewtonIntegrator(homotheticCellResize=1,damping=.4),
+ #Peri3dController(goal=Matrix3(-.1,0,0, 0,-.1,0, 0,0,-.1),stressMask=0,maxStrainRate=1e-1,label='ctrl'),
+ #Peri3dController(goal=Matrix3(0,.2,0, 0,0,0, 0,0,0),stressMask=0,maxStrainRate=1,label='ctrl'),
+ Peri3dController(goal=Matrix3(-1e3,.1,0, 0,.1,0, 0,0,0),stressMask=0b0001,maxStrainRate=1,label='ctrl'),
+ PeriodicPythonRunner(iterPeriod=10,command='addData()')
+]
+O.step()
+from yade import plot
+def addData():
+ plot.addData(sxx=ctrl.stress[0,0],syy=ctrl.stress[1,1],szz=ctrl.stress[2,2],exx=ctrl.strain[0,0],eyy=ctrl.strain[1,1],ezz=ctrl.strain[2,2],t=O.time)
+
+plot.plots={'exx':('sxx'),}
+
