yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05654
[Branch ~yade-dev/yade/trunk] Rev 2429: 1. The mindlin law now prevents granular ratcheting by default.
------------------------------------------------------------
revno: 2429
committer: Chiara Modenese <chia@engs-018373>
branch nick: trunk
timestamp: Thu 2010-09-09 15:02:49 +0100
message:
1. The mindlin law now prevents granular ratcheting by default.
2. Ig2_Sphere_Sphere_ScGeom has now avoidGranularRatcheting bool (true by default, which was the hardcoded value till now as well); add documentation from the source to the real doc, add a few references on that.
3. SpherePack no longer saves vectors as tuples; this also fixes memoization issues with pack.random{Dense,Peri}Pack
4. Add Law2_ScGeom_MindlinPhys_HertzWithLinearShear used for testing differences between linear and non-linear law.
modified:
doc/references.bib
doc/yade-conferences.bib
pkg/dem/DataClass/SpherePack.cpp
pkg/dem/DataClass/SpherePack.hpp
pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp
pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.hpp
pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp
pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp
pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp
py/pack/_packSpheres.cpp
py/pack/pack.py
py/plot.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-07-16 18:12:30 +0000
+++ doc/references.bib 2010-09-09 14:02:49 +0000
@@ -86,6 +86,27 @@
url = {http://www.comphys.ethz.ch/hans/p/334.pdf}
}
+@article{McNamara2008,
+ title={Microscopic origin of granular ratcheting},
+ author={S. McNamara and R. GarcÃa-Rojo and H. J. Herrmann},
+ journal={Physical Review E},
+ year={2008},
+ number={3},
+ volume={77},
+ doi={11.1103/PhysRevE.77.031304}
+}
+
+@InProceedings{GarciaRojo2004,
+ author = {GarcÃa-Rojo, R. and McNamara, S. and Herrmann, H. J.},
+ title = {Discrete element methods for the micro-mechanical investigation of granular ratcheting},
+ booktitle = {Proceedings ECCOMAS 2004},
+ year = {2004},
+ address = {Jyvaskyla},
+ url = {http://www.ica1.uni-stuttgart.de/publications/2004/GMH04}
+}
+
+
+
@book{Allen1989,
author = {Allen, M. P. and Tildesley, D. J.},
title = {Computer simulation of liquids},
=== modified file 'doc/yade-conferences.bib'
--- doc/yade-conferences.bib 2010-07-16 18:12:30 +0000
+++ doc/yade-conferences.bib 2010-09-09 14:02:49 +0000
@@ -5,7 +5,7 @@
booktitle={Proceedings of the International Conference on Modelling and Simulation 2010, Prague},
year={2010},
month={June},
- publisher={preprint},
+ url={https://yade-dem.org/w/images/6/64/Stransky2010-Macroscopic-elastic-properties-of-particle-models.pdf}
}
=== modified file 'pkg/dem/DataClass/SpherePack.cpp'
--- pkg/dem/DataClass/SpherePack.cpp 2010-09-09 09:51:23 +0000
+++ pkg/dem/DataClass/SpherePack.cpp 2010-09-09 14:02:49 +0000
@@ -28,6 +28,11 @@
const python::tuple& t=python::extract<python::tuple>(l[i]);
python::extract<Vector3r> vec(t[0]);
if(vec.check()) { pack.push_back(Sph(vec(),python::extract<double>(t[1]))); continue; }
+ #if 0
+ // compatibility block
+ python::extract<python::tuple> tup(t[0]);
+ if(tup.check()) { pack.push_back(Sph(Vector3r(python::extract<double>(tup[0]),python::extract<double>(tup[1]),python::extract<double>(tup[2]))),python::extract<double>(t[1])); continue; }
+ #endif
PyErr_SetString(PyExc_TypeError, "List elements must be (tuple or Vector3, float)!");
python::throw_error_already_set();
}
@@ -38,12 +43,14 @@
FOREACH(const Sph& s, pack) ret.append(python::make_tuple(s.c,s.r));
return ret;
};
-
+#if 0
python::list SpherePack::toList_pointsAsTuples() const {
python::list ret;
FOREACH(const Sph& s, pack) ret.append(python::make_tuple(s.c[0],s.c[1],s.c[2],s.r));
return ret;
};
+#endif
+
void SpherePack::fromFile(string file) {
typedef pair<Vector3r,Real> pairVector3rReal;
=== modified file 'pkg/dem/DataClass/SpherePack.hpp'
--- pkg/dem/DataClass/SpherePack.hpp 2010-09-02 10:09:05 +0000
+++ pkg/dem/DataClass/SpherePack.hpp 2010-09-09 14:02:49 +0000
@@ -47,7 +47,9 @@
// I/O
void fromList(const python::list& l);
python::list toList() const;
- python::list toList_pointsAsTuples() const;
+ #if 0
+ python::list toList_pointsAsTuples() const;
+ #endif
void fromFile(const string file);
void toFile(const string file) const;
void fromSimulation();
=== modified file 'pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp 2010-07-18 19:21:08 +0000
+++ pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp 2010-09-09 14:02:49 +0000
@@ -39,7 +39,7 @@
scm->penetrationDepth=penetrationDepth;
scm->radius1=s1->radius;
scm->radius2=s2->radius;
- scm->precompute(state1,state2,scene,c,normal,isNew,true);
+ scm->precompute(state1,state2,scene,c,normal,isNew,avoidGranularRatcheting);
return true;
}
return false;
=== modified file 'pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.hpp'
--- pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.hpp 2010-08-24 12:54:14 +0000
+++ pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.hpp 2010-09-09 14:02:49 +0000
@@ -20,6 +20,17 @@
#endif
YADE_CLASS_BASE_DOC_ATTRS(Ig2_Sphere_Sphere_ScGeom,InteractionGeometryFunctor,"Create/update a :yref:`ScGeom` instance representing intersection of two :yref:`Spheres<Sphere>`.",
((Real,interactionDetectionFactor,1,,"Enlarge both radii by this factor (if >1), to permit creation of distant interactions.\n\nInteractionGeometry will be computed when interactionDetectionFactor*(rad1+rad2) > distance.\n\n.. note::\n\t This parameter is functionally coupled with :yref:`Bo1_Sphere_Aabb::aabbEnlargeFactor`, which will create larger bounding boxes and should be of the same value.\n\n.. warning::\n\tFunctionally equal class :yref:`Ig2_Sphere_Sphere_Dem3DofGeom` (which creates :yref:`Dem3DofGeom` rather than :yref:`ScGeom`) calls this parameter :yref:`distFactor<Ig2_Sphere_Sphere_Dem3DofGeom::distFactor>`, but its semantics is *different* in some aspects."))
+ ((bool,avoidGranularRatcheting,true,,"Granular ratcheting is mentioned in [GarciaRojo2004]_, [Alonso2004]_, [McNamara2008]_.\n\n"
+ "Unfortunately, published papers tend to focus on the \"good\" ratcheting, i.e. irreversible deformations due to the intrinsic nature of frictional granular materials, while a 2004 talk of McNamara in Paris clearly mentioned a possible \"bad\" ratcheting, purely linked with the formulation of the contact laws in what he called \"molecular dynamics\" (i.e. Cundall's model, as opposed to \"contact dynamics\" from Moreau and Jean).\n\n"
+ "The bad ratcheting is best understood considering this small elastic cycle at a contact between two grains: assuming b1 is fixed, impose this displacement to b2:\n\n"
+ "#. translation *dx* in the normal direction\n"
+ "#. rotation *a*\n"
+ "#. translation *-dx* (back to the initial position)\n"
+ "#. rotation *-a* (back to the initial orientation)\n\n\n"
+ "If the branch vector used to define the relative shear in rotation×branch is not constant (typically if it is defined from the vector center→contactPoint), then the shear displacement at the end of this cycle is not zero: rotations *a* and *-a* are multiplied by branches of different lengths.\n\n"
+ "It results in a finite contact force at the end of the cycle even though the positions and orientations are unchanged, in total contradiction with the elastic nature of the problem. It could also be seen as an *inconsistent energy creation or loss*. Given that DEM simulations tend to generate oscillations around equilibrium (damped mass-spring), it can have a significant impact on the evolution of the packings, resulting for instance in slow creep in iterations under constant load.\n\n"
+ "The solution to avoid that is quite simple in the case of linear-elastic laws: use a constant branch vector, which is what this functor does by default."
+ ))
);
FUNCTOR2D(Sphere,Sphere);
// needed for the dispatcher, even if it is symmetric
=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-07-18 19:21:08 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-09-09 14:02:49 +0000
@@ -49,57 +49,53 @@
void Law2_ScGeom_FrictPhys_Basic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact){
int id1 = contact->getId1(), id2 = contact->getId2();
- ScGeom* currentContactGeometry= static_cast<ScGeom*>(ig.get());
- FrictPhys* currentContactPhysics = static_cast<FrictPhys*>(ip.get());
- if(currentContactGeometry->penetrationDepth <0){
+ ScGeom* geom= static_cast<ScGeom*>(ig.get());
+ FrictPhys* phys = static_cast<FrictPhys*>(ip.get());
+ if(geom->penetrationDepth <0){
if (neverErase) {
- currentContactPhysics->shearForce = Vector3r::Zero();
- currentContactPhysics->normalForce = Vector3r::Zero();}
+ phys->shearForce = Vector3r::Zero();
+ phys->normalForce = Vector3r::Zero();}
else scene->interactions->requestErase(id1,id2);
return;}
State* de1 = Body::byId(id1,scene)->state.get();
State* de2 = Body::byId(id2,scene)->state.get();
- Real& un=currentContactGeometry->penetrationDepth;
- TRVAR3(currentContactGeometry->penetrationDepth,de1->se3.position,de2->se3.position);
- currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal;
+ Real& un=geom->penetrationDepth;
+ TRVAR3(geom->penetrationDepth,de1->se3.position,de2->se3.position);
+ phys->normalForce=phys->kn*std::max(un,(Real) 0)*geom->normal;
- Vector3r& shearForce = currentContactGeometry->rotate(currentContactPhysics->shearForce);
- const Vector3r& shearDisp = currentContactGeometry->shearIncrement();
+ Vector3r& shearForce = geom->rotate(phys->shearForce);
+ const Vector3r& shearDisp = geom->shearIncrement();
+ shearForce -= phys->ks*shearDisp;
+ Real maxFs = phys->normalForce.squaredNorm()*std::pow(phys->tangensOfFrictionAngle,2);
if (!traceEnergy){//Update force but don't compute energy terms (see below))
- shearForce -= currentContactPhysics->ks*shearDisp;
// PFC3d SlipModel, is using friction angle. CoulombCriterion
- Real maxFs = currentContactPhysics->normalForce.squaredNorm()*
- std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
if( shearForce.squaredNorm() > maxFs ){
Real ratio = Mathr::Sqrt(maxFs) / shearForce.norm();
shearForce *= ratio;}
} else {
//almost the same with additional Vector3r instanciated for energy tracing, duplicated block to make sure there is no cost for the instanciation of the vector when traceEnergy==false
- shearForce -= currentContactPhysics->ks*shearDisp;
- Real maxFs = currentContactPhysics->normalForce.squaredNorm()*
- std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
if( shearForce.squaredNorm() > maxFs ){
Real ratio = Mathr::Sqrt(maxFs) / shearForce.norm();
Vector3r trialForce=shearForce;//store prev force for definition of plastic slip
//define the plastic work input and increment the total plastic energy dissipated
shearForce *= ratio;
plasticDissipation +=
- ((1/currentContactPhysics->ks)*(trialForce-shearForce))//plastic disp.
+ ((1/phys->ks)*(trialForce-shearForce))//plastic disp.
.dot(shearForce);//active force
}
}
if (!scene->isPeriodic)
- applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
+ applyForceAtContactPoint(-phys->normalForce-shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
else {//we need to use correct branches in the periodic case, the following apply for spheres only
- Vector3r force = -currentContactPhysics->normalForce-shearForce;
+ Vector3r force = -phys->normalForce-shearForce;
scene->forces.addForce(id1,force);
scene->forces.addForce(id2,-force);
- scene->forces.addTorque(id1,(currentContactGeometry->radius1-0.5*currentContactGeometry->penetrationDepth)* currentContactGeometry->normal.cross(force));
- scene->forces.addTorque(id2,(currentContactGeometry->radius2-0.5*currentContactGeometry->penetrationDepth)* currentContactGeometry->normal.cross(force));
+ scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force));
+ scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force));
}
- //FIXME : make sure currentContactPhysics->prevNormal is not used anywhere, remove it from physics and remove this line :
- currentContactPhysics->prevNormal = currentContactGeometry->normal;
+ //FIXME : make sure phys->prevNormal is not used anywhere, remove it from physics and remove this line :
+ phys->prevNormal = geom->normal;
}
// same as elasticContactLaw, but using Dem3DofGeom
=== modified file 'pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp'
--- pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp 2010-09-09 09:51:23 +0000
+++ pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp 2010-09-09 14:02:49 +0000
@@ -8,7 +8,7 @@
#define pi 3.14159265
-YADE_PLUGIN((MindlinPhys)(Ip2_FrictMat_FrictMat_MindlinPhys)(Law2_ScGeom_MindlinPhys_Mindlin));
+YADE_PLUGIN((MindlinPhys)(Ip2_FrictMat_FrictMat_MindlinPhys)(Law2_ScGeom_MindlinPhys_HertzWithLinearShear)(Law2_ScGeom_MindlinPhys_Mindlin));
Real Law2_ScGeom_MindlinPhys_Mindlin::Real0=0;
Real Law2_ScGeom_MindlinPhys_Mindlin::getfrictionDissipation() {return (Real) frictionDissipation;}
@@ -99,14 +99,51 @@
return normEnergy;
}
+void Law2_ScGeom_MindlinPhys_HertzWithLinearShear::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact){
+ Body::id_t id1(contact->getId1()), id2(contact->getId2());
+ ScGeom* geom = static_cast<ScGeom*>(ig.get());
+ MindlinPhys* phys=static_cast<MindlinPhys*>(ip.get());
+ const Real uN=geom->penetrationDepth;
+ if (uN<0) {scene->interactions->requestErase(id1,id2); return;}
+ // normal force
+ Real Fn=phys->kno*pow(uN,3/2.);
+ phys->normalForce=Fn*geom->normal;
+ //phys->kn=3./2.*phys->kno*std::pow(uN,0.5); // update stiffness, not needed
+
+ // shear force
+ Vector3r& Fs=geom->rotate(phys->shearForce);
+ Real ks= nonLin>0 ? phys->kso*std::pow(uN,0.5) : phys->kso;
+ Vector3r shearIncrement;
+ if(nonLin>1){
+ State *de1=Body::byId(id1,scene)->state.get(), *de2=Body::byId(id2,scene)->state.get();
+ Vector3r shiftVel=scene->isPeriodic ? Vector3r(scene->cell->velGrad*scene->cell->Hsize*contact->cellDist.cast<Real>()) : Vector3r::Zero();
+ Vector3r incidentV = geom->getIncidentVel(de1, de2, scene->dt, shiftVel, /*preventGranularRatcheting*/ nonLin>2 );
+ Vector3r incidentVn = geom->normal.dot(incidentV)*geom->normal; // contact normal velocity
+ Vector3r incidentVs = incidentV-incidentVn; // contact shear velocity
+ shearIncrement=incidentVs*scene->dt;
+ } else { shearIncrement=geom->shearIncrement(); }
+ Fs-=ks*shearIncrement;
+ // Mohr-Coulomb slip
+ Real maxFs2=pow(Fn,2)*pow(phys->tangensOfFrictionAngle,2);
+ if(Fs.squaredNorm()>maxFs2) Fs*=sqrt(maxFs2)/Fs.norm();
+
+ // apply forces
+ Vector3r f=-phys->normalForce-phys->shearForce; /* should be a reference returned by geom->rotate */ assert(phys->shearForce==Fs);
+ scene->forces.addForce(id1,f);
+ scene->forces.addForce(id2,-f);
+ scene->forces.addTorque(id1,(geom->radius1-.5*geom->penetrationDepth)*geom->normal.cross(f));
+ scene->forces.addTorque(id2,(geom->radius2-.5*geom->penetrationDepth)*geom->normal.cross(f));
+}
+
+
/******************** Law2_ScGeom_MindlinPhys_Mindlin *********/
CREATE_LOGGER(Law2_ScGeom_MindlinPhys_Mindlin);
void Law2_ScGeom_MindlinPhys_Mindlin::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact){
const Real& dt = scene->dt; // get time step
- int id1 = contact->getId1(); // get id body 1
- int id2 = contact->getId2(); // get id body 2
+ Body::id_t id1 = contact->getId1(); // get id body 1
+ Body::id_t id2 = contact->getId2(); // get id body 2
State* de1 = Body::byId(id1,scene)->state.get();
State* de2 = Body::byId(id2,scene)->state.get();
@@ -271,7 +308,7 @@
/****************/
if (!scene->isPeriodic)
- applyForceAtContactPoint(-phys->normalForce - phys->shearForce, scg->contactPoint , id1, de1->se3.position, id2, de2->se3.position);
+ applyForceAtContactPoint(-phys->normalForce - phys->shearForce, scg->contactPoint , id1, de1->se3.position, id2, de2->se3.position);
else { // in scg we do not wrap particles positions, hence "applyForceAtContactPoint" cannot be used
Vector3r force = -phys->normalForce - phys->shearForce;
scene->forces.addForce(id1,force);
=== modified file 'pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp'
--- pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp 2010-09-09 09:51:23 +0000
+++ pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp 2010-09-09 14:02:49 +0000
@@ -64,6 +64,16 @@
};
REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_MindlinPhys);
+class Law2_ScGeom_MindlinPhys_HertzWithLinearShear: public LawFunctor{
+ public:
+ virtual void go(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*);
+ FUNCTOR2D(ScGeom,MindlinPhys);
+ YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_MindlinPhys_HertzWithLinearShear,LawFunctor,
+ "Constitutive law for the Hertz formulation (using :yref:`MindlinPhys.kno`) and linear beahvior in shear (using :yref:`MindlinPhys.kso` for stiffness and :yref:`FrictPhys.tangensOfFrictionAngle`). \n\n.. note:: No viscosity or damping. If you need those, look at :yref:`Law2_ScGeom_MindlinPhys_Mindlin`, which also includes non-linear Mindlin shear.",
+ ((int,nonLin,0,,"Shear force nonlinearity (the value determines how many features of the non-linearity are taken in account). 1: ks as in HM 2: shearElastic increment computed as in HM 3. granular ratcheting disabled."))
+ );
+};
+REGISTER_SERIALIZABLE(Law2_ScGeom_MindlinPhys_HertzWithLinearShear);
/******************** Law2_ScGeom_MindlinPhys_Mindlin *********/
class Law2_ScGeom_MindlinPhys_Mindlin: public LawFunctor{
@@ -85,7 +95,7 @@
FUNCTOR2D(ScGeom,MindlinPhys);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_MindlinPhys_Mindlin,LawFunctor,"Constitutive law for the Hertz-Mindlin formulation. It includes non linear elasticity in the normal direction as predicted by Hertz for two non-conforming elastic contact bodies. In the shear direction, instead, it reseambles the simplified case without slip discussed in Mindlin's paper, where a linear relationship between shear force and tangential displacement is provided. Finally, the Mohr-Coulomb criterion is employed to established the maximum friction force which can be developed at the contact. Moreover, it is also possible to include the effect of linear viscous damping through the definition of the parameters $\\beta_{n}$ and $\\beta_{s}$.",
- ((bool,preventGranularRatcheting,false,,"bool to avoid granular ratcheting"))
+ ((bool,preventGranularRatcheting,true,,"bool to avoid granular ratcheting"))
((bool,includeAdhesion,false,,"bool to include the adhesion force following the DMT formulation. If true, also the normal elastic energy takes into account the adhesion effect."))
((bool,useDamping,false,,"bool to include contact damping"))
((bool,calcEnergy,false,,"bool to calculate energy terms (shear potential energy, dissipation of energy due to friction and dissipation of energy due to normal and tangential damping)"))
=== modified file 'py/pack/_packSpheres.cpp'
--- py/pack/_packSpheres.cpp 2010-09-09 09:51:23 +0000
+++ py/pack/_packSpheres.cpp 2010-09-09 14:02:49 +0000
@@ -9,7 +9,9 @@
python::class_<SpherePack>("SpherePack","Set of spheres represented as centers and radii. This class is returned by :yref:`yade.pack.randomDensePack`, :yref:`yade.pack.randomPeriPack` and others. The object supports iteration over spheres, as in \n\n\t>>> sp=SpherePack()\n\t>>> for center,radius in sp: print center,radius\n\n\t>>> for sphere in sp: print sphere[0],sphere[1] ## same, but without unpacking the tuple automatically\n\n\t>>> for i in range(0,len(sp)): print sp[i][0], sp[i][1] ## same, but accessing spheres by index\n",python::init<python::optional<python::list> >(python::args("list"),"Empty constructor, optionally taking list [ ((cx,cy,cz),r), … ] for initial data." ))
.def("add",&SpherePack::add,"Add single sphere to packing, given center as 3-tuple and radius")
.def("toList",&SpherePack::toList,"Return packing data as python list.")
+ #if 0
.def("toList_pointsAsTuples",&SpherePack::toList_pointsAsTuples,"Return packing data as python list, but using only pure-python data types (3-tuples instead of Vector3) (for pickling with cPickle)")
+ #endif
.def("fromList",&SpherePack::fromList,"Make packing from given list, same format as for constructor. Discards current data.")
.def("load",&SpherePack::fromFile,(python::arg("fileName")),"Load packing from external text file (current data will be discarded).")
.def("save",&SpherePack::toFile,(python::arg("fileName")),"Save packing to external text file (will be overwritten).")
=== modified file 'py/pack/pack.py'
--- py/pack/pack.py 2010-07-01 10:31:49 +0000
+++ py/pack/pack.py 2010-09-09 14:02:49 +0000
@@ -283,7 +283,7 @@
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))
+ packBlob=buffer(cPickle.dumps(sp.toList(),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()
@@ -441,8 +441,7 @@
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
- not being able to add any more spheres), then we run periodic compression with PeriIsoCompressor, just like with
+ A cell of initSize is stuffed with as many spheres as possible, then we run periodic compression with PeriIsoCompressor, just like with
randomDensePack.
:param radius: mean sphere radius
@@ -450,9 +449,6 @@
:param initSize: initial size of the periodic cell.
:return: SpherePack object, which also contains periodicity information.
-
- :todo: memoization in db; what criteria??
-
"""
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)
@@ -462,7 +458,7 @@
O.periodic=True
O.cell.refSize=Vector3(initSize)
sp.makeCloud(Vector3().Zero,O.cell.refSize,radius,rRelFuzz,-1,True)
- O.engines=[ForceResetter(),BoundDispatcher([Bo1_Sphere_Aabb()]),InsertionSortCollider(nBins=2,sweepLength=.05*radius),InteractionDispatchers([Ig2_Sphere_Sphere_Dem3DofGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_Dem3DofGeom_FrictPhys_Basic()]),PeriIsoCompressor(charLen=2*radius,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=20,keepProportions=True),NewtonIntegrator(damping=.8)]
+ O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb()],nBins=2,sweepLength=.05*radius),InteractionDispatchers([Ig2_Sphere_Sphere_Dem3DofGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_Dem3DofGeom_FrictPhys_Basic()]),PeriIsoCompressor(charLen=2*radius,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=20,keepProportions=True),NewtonIntegrator(damping=.8)]
O.materials.append(FrictMat(young=30e9,frictionAngle=.1,poisson=.3,density=1e3))
for s in sp: O.bodies.append(utils.sphere(s[0],s[1]))
O.dt=utils.PWaveTimeStep()
=== modified file 'py/plot.py'
--- py/plot.py 2010-09-09 09:51:23 +0000
+++ py/plot.py 2010-09-09 14:02:49 +0000
@@ -182,6 +182,8 @@
pylab.legend([xlateLabel(_p[0]) for _p in plots_p_y2],loc='upper right')
pylab.rcParams['lines.color']=origLinesColor
pylab.ylabel((', '.join([xlateLabel(_p[0]) for _p in plots_p_y2])) if p not in xylabels or len(xylabels[p])<3 or not xylabels[p][2] else xylabels[p][2])
+ ## should be repeated for y2 axis, but in that case the 10^.. label goes to y1 axis (bug in matplotlib, it seems)
+ # if scientific: pylab.ticklabel_format(style='sci',scilimits=(0,0),axis='both')
if 'title' in O.tags.keys(): pylab.title(O.tags['title'])
