yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06554
[Branch ~yade-dev/yade/trunk] Rev 2614: - finish implementation of contacts tracking on chained cylinders.
------------------------------------------------------------
revno: 2614
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Wed 2010-12-15 11:11:41 +0100
message:
- finish implementation of contacts tracking on chained cylinders.
- some cleaning and doc update
modified:
pkg/common/Cylinder.cpp
pkg/common/Cylinder.hpp
pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp
pkg/dem/ScGeom.cpp
scripts/test/chained-cylinder-roots.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 'pkg/common/Cylinder.cpp'
--- pkg/common/Cylinder.cpp 2010-12-13 12:11:43 +0000
+++ pkg/common/Cylinder.cpp 2010-12-15 10:11:41 +0000
@@ -5,6 +5,7 @@
#endif
#include<yade/pkg/common/Aabb.hpp>
#include<yade/pkg/dem/FrictPhys.hpp>
+#include "../../lib/base/Math.hpp"
Cylinder::~Cylinder(){}
ChainedCylinder::~ChainedCylinder(){}
@@ -13,7 +14,6 @@
// Ig2_Sphere_ChainedCylinder_CylScGeom::~Ig2_Sphere_ChainedCylinder_CylScGeom() {}
// Ig2_ChainedCylinder_ChainedCylinder_ScGeom6D::~Ig2_ChainedCylinder_ChainedCylinder_ScGeom6D() {}
-
YADE_PLUGIN((Cylinder)(ChainedCylinder)(ChainedState)(CylScGeom)(Ig2_Sphere_ChainedCylinder_CylScGeom)(Ig2_ChainedCylinder_ChainedCylinder_ScGeom6D)(Law2_CylScGeom_FrictPhys_CundallStrack)
#ifdef YADE_OPENGL
(Gl1_Cylinder)(Gl1_ChainedCylinder)
@@ -33,7 +33,6 @@
const State& state1, const State& state2, const Vector3r& shift2, const bool& force,
const shared_ptr<Interaction>& c)
{
-// cerr<<"Ig2_Sphere_ChainedCylinder_CylScGeom::go"<<endl;
const State* sphereSt=YADE_CAST<const State*>(&state1);
const ChainedState* cylinderSt=YADE_CAST<const ChainedState*>(&state2);
ChainedCylinder *cylinder=YADE_CAST<ChainedCylinder*>(cm2.get());
@@ -41,48 +40,94 @@
assert(sphereSt && cylinderSt && cylinder && sphere);
bool isLast = (cylinderSt->chains[cylinderSt->chainNumber].size()==(cylinderSt->rank+1));
bool isNew = !c->geom;
-// if (cylinderSt->chains[cylinderSt->chainNumber].size()==(cylinderSt->rank+1)) {cerr << "last cylinder - ignored"<<endl; return false;}
+
+ shared_ptr<const ChainedState> statePrev;
+ if (cylinderSt->rank>0)
+ statePrev = YADE_PTR_CAST<const ChainedState> (Body::byId(cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank-1],scene)->state);
+
+ shared_ptr<CylScGeom> scm;
+ if (!isNew) {scm=YADE_PTR_CAST<CylScGeom>(c->geom);}
//FIXME : definition of segment in next line breaks periodicity
-// cerr << cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1]<<endl;
-// cerr<<Body::byId(cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1],scene)->state->pos<<endl;
shared_ptr<Body> cylinderNext;
Vector3r segment, branch, direction;
Real length, dist;
+ branch = sphereSt->pos-cylinderSt->pos-shift2;
if (!isLast) {
cylinderNext = Body::byId(cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1],scene);
segment = cylinderNext->state->pos-cylinderSt->pos;
- branch = sphereSt->pos-cylinderSt->pos-shift2;
- if ((segment.dot(branch)>(segment.dot(segment)/*+interactionDetectionFactor*cylinder->radius*/) && isNew)) return false;//position _after_ end of cylinder
+ if (segment.dot(branch)>(segment.dot(segment)/*+interactionDetectionFactor*cylinder->radius*/)) {//position _after_ end of cylinder
+ //FIXME : scm->penetrationDepth=-1 is defined to workaround interactions never being erased when scm->isDuplicate=2 on the true interaction.
+ if (isNew) return false; else if (scm->isDuplicate) {scm->isDuplicate=2; scm->penetrationDepth=-1; return true;}}
length = segment.norm();
direction = segment/length;
dist = direction.dot(branch);
- if ((dist<-interactionDetectionFactor*cylinder->radius) && isNew) return false;
- if (cylinderSt->rank>0){//make sure there is no contact with the previous element in the chain, or consider it a duplicate and continue
- const shared_ptr<Body> cylinderPrev = Body::byId(cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank-1],scene);
- Vector3r branchP = sphereSt->pos-cylinderPrev->state->pos;
- Vector3r segmentP = cylinderSt->pos - cylinderPrev->state->pos;
- if (segmentP.dot(branchP)<segmentP.dot(segmentP)) return false;//will give false on existing interaction, as expected
- }
+ if (dist<-interactionDetectionFactor*cylinder->radius &&
+ branch.squaredNorm() > pow(interactionDetectionFactor*(sphere->radius+cylinder->radius),2)) {
+ if (isNew) return false; else if (scm->isDuplicate) {scm->isDuplicate=2; scm->penetrationDepth=-1; return true;}}
} else {//handle the last node with length=0
segment = Vector3r(0,0,0);
- branch = sphereSt->pos-cylinderSt->pos-shift2;
length = 0;
direction = Vector3r(0,1,0);
- dist = 0;}
+ dist = 0;
+ if (branch.squaredNorm() > interactionDetectionFactor*(sphere->radius+cylinder->radius)) {
+ if (isNew) return false; else if (scm->isDuplicate) {scm->isDuplicate=2; scm->penetrationDepth=-1; return true;}}
+ }
//Check sphere-cylinder distance
Vector3r projectedP = cylinderSt->pos+shift2 + direction*dist;
branch = projectedP-sphereSt->pos;
- if ((branch.squaredNorm()>(pow(sphere->radius+cylinder->radius,2))) && isNew) return false;
-
- shared_ptr<CylScGeom> scm;
- if(!isNew) scm=YADE_PTR_CAST<CylScGeom>(c->geom);
- else { scm=shared_ptr<CylScGeom>(new CylScGeom()); c->geom=scm; }
+ if (branch.squaredNorm()>(pow(sphere->radius+cylinder->radius,2))) {
+ if (isNew) return false; else if (scm->isDuplicate) {scm->isDuplicate=2; scm->penetrationDepth=-1; return true;}}
+
+ if (!isNew) scm->isDuplicate = false;//reset here at each step, and recompute below
+
+ //make sure there is no contact with the previous element in the chain, else consider this one a duplicate and get data from the other contact. two interactions will share the same geometry and physics.
+ if (cylinderSt->rank>0 && dist<=0) {
+ Vector3r branchP = sphereSt->pos - statePrev->pos;
+ Vector3r segmentP = cylinderSt->pos - statePrev->pos;
+ if (segmentP.dot(branchP)<segmentP.dot(segmentP)) {
+ if (isNew) {
+ const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank-1]);
+ assert(intr);//we know there is a contact, so there should be at least a virtual interaction created by collider
+ if (!intr->geom || !intr->phys) return false;
+ else {
+ c->geom = intr->geom;
+ c->phys = intr->phys;
+ scm=YADE_PTR_CAST<CylScGeom>(c->geom);
+// scm->duplicate = intr;
+ scm->isDuplicate = true;
+ isNew = false;
+ return true;}
+ } else scm->isDuplicate=true;
+ scm->trueInt = cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank-1];
+ }
+ }
+ //similarly, make sure there is no contact with the next element in the chain
+ else if (!isLast && dist>length) {
+ if ( (cylinderNext->state->pos-sphereSt->pos).squaredNorm() < pow(sphere->radius+cylinder->radius,2)) {
+ if (isNew) {
+ const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1]);
+ assert(intr);
+ if (!intr->geom || !intr->phys) return false;
+ c->geom = intr->geom;
+ c->phys = intr->phys;
+ scm=YADE_PTR_CAST<CylScGeom>(c->geom);
+// scm->duplicate = intr;
+ scm->isDuplicate = true;
+ isNew = false;
+ return true;
+ } else scm->isDuplicate=true;
+ scm->trueInt = cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1];
+ }
+ }
+
+ //We didn't find any special case, do normal geometry definition
+ if (isNew) { scm=shared_ptr<CylScGeom>(new CylScGeom()); c->geom=scm;}
scm->radius1=sphere->radius;
scm->radius2=cylinder->radius;
- if (!isLast) scm->id3=cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1];
+ if (!isLast && !scm->id3) scm->id3=cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1];
scm->start=cylinderSt->pos+shift2; scm->end=scm->start+segment;
//FIXME : there should be other checks without distanceFactor?
@@ -167,7 +212,6 @@
//bs1->segment used for fast BBs and projections + display
bs1->segment= bchain2.pos-bchain1.pos;
#ifdef YADE_OPENGL
- //bs1->length and s1->chainedOrientation used for display only,
bs1->length=length;
bs1->chainedOrientation.setFromTwoVectors(Vector3r::UnitZ(),bchain1.ori.conjugate()*segt);
#endif
@@ -236,18 +280,11 @@
void Gl1_Cylinder::drawCylinder(bool wire, Real radius, Real length, const Quaternionr& shift) const
{
-// GLERROR;
-/* if (glCylinderList<0) {
- glCylinderList = glGenLists(1);
- glNewList(glCylinderList,GL_COMPILE);*/
glPushMatrix();
GLUquadricObj *quadObj = gluNewQuadric();
gluQuadricDrawStyle(quadObj, (GLenum) (wire ? GLU_SILHOUETTE : GLU_FILL));
gluQuadricNormals(quadObj, (GLenum) GLU_SMOOTH);
gluQuadricOrientation(quadObj, (GLenum) GLU_OUTSIDE);
-// glTranslatef(0.0,0.0,-length*0.5);
- //scaling needs to adapt spheres or they will be elipsoids. They actually seem to disappear when commented glList code is uncommented, the cylinders are displayed correclty.
-// glScalef(1,length,1);
AngleAxisr aa(shift);
glRotatef(aa.angle()*180.0/Mathr::PI,aa.axis()[0],aa.axis()[1],aa.axis()[2]);
gluCylinder(quadObj, radius, radius, length, glutSlices,glutStacks);
@@ -259,43 +296,8 @@
// gluDisk(quadObj,0.0,radius,glutSlices,_loops);
gluDeleteQuadric(quadObj);
glPopMatrix();
-// GLERROR;
-
-// glEndList();}
-// glCallList(glCylinderList);
-
}
-// void Gl1_Cylinder::drawCylinder(bool wire, Real radius, Real length, const Quaternionr& shift) const
-// {
-// // GLERROR;
-// /* if (glCylinderList<0) {
-// glCylinderList = glGenLists(1);
-// glNewList(glCylinderList,GL_COMPILE);*/
-// glPushMatrix();
-// GLUquadricObj *quadObj = gluNewQuadric();
-// gluQuadricDrawStyle(quadObj, (GLenum) (wire ? GLU_SILHOUETTE : GLU_FILL));
-// gluQuadricNormals(quadObj, (GLenum) GLU_SMOOTH);
-// gluQuadricOrientation(quadObj, (GLenum) GLU_OUTSIDE);
-// // glTranslatef(0.0,0.0,-length*0.5);
-// //scaling needs to adapt spheres or they will be elipsoids. They actually seem to disappear when commented glList code is uncommented, the cylinders are displayed correclty.
-// // glScalef(1,length,1);
-// gluCylinder(quadObj, radius, radius, length, glutSlices,glutStacks);
-// gluQuadricOrientation(quadObj, (GLenum) GLU_INSIDE);
-// glutSolidSphere(radius,glutSlices,glutStacks);
-// glTranslatef(0.0,0.0,length);
-// AngleAxisr aa(shift);
-// glRotatef(aa.angle(),aa.axis()[0],aa.axis()[1],aa.axis()[2]);
-// glutSolidSphere(radius,glutSlices,glutStacks);
-// // gluDisk(quadObj,0.0,radius,glutSlices,_loops);
-// gluDeleteQuadric(quadObj);
-// glPopMatrix();
-// // GLERROR;
-//
-// // glEndList();}
-// // glCallList(glCylinderList);
-//
-// }
//!################## BOUNDS FUNCTOR #####################
@@ -315,23 +317,8 @@
}
return;
}
- // adjust box size along axes so that cylinder doesn't stick out of the box even if sheared (i.e. parallelepiped)
-// if(scene->cell->hasShear()) {
-// Vector3r refHalfSize(minkSize);
-// const Vector3r& cos=scene->cell->getCos();
-// for(int i=0; i<3; i++){
-// //cerr<<"cos["<<i<<"]"<<cos[i]<<" ";
-// int i1=(i+1)%3,i2=(i+2)%3;
-// minkSize[i1]+=.5*refHalfSize[i1]*(1/cos[i]-1);
-// minkSize[i2]+=.5*refHalfSize[i2]*(1/cos[i]-1);
-// }
-// }
-// //cerr<<" || "<<halfSize<<endl;
-// aabb->min = scene->cell->unshearPt(se3.position)-minkSize;
-// aabb->max = scene->cell->unshearPt(se3.position)+minkSize;
}
-
void Bo1_ChainedCylinder_Aabb::go(const shared_ptr<Shape>& cm, shared_ptr<Bound>& bv, const Se3r& se3, const Body* b){
Cylinder* cylinder = static_cast<Cylinder*>(cm.get());
if(!bv){ bv=shared_ptr<Bound>(new Aabb); }
@@ -351,7 +338,7 @@
void Law2_CylScGeom_FrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
int id1 = contact->getId1(), id2 = contact->getId2();
- CylScGeom* geom= static_cast<CylScGeom*>(ig.get());
+ CylScGeom* geom= static_cast<CylScGeom*>(ig.get());
FrictPhys* phys = static_cast<FrictPhys*>(ip.get());
if(geom->penetrationDepth <0){
if (neverErase) {
@@ -359,6 +346,13 @@
phys->normalForce = Vector3r::Zero();}
else scene->interactions->requestErase(id1,id2);
return;}
+ if (geom->isDuplicate) {
+ if (id2!=geom->trueInt) {
+// cerr<<"skip duplicate "<<id1<<" "<<id2<<endl;
+ if (geom->isDuplicate==2) scene->interactions->requestErase(id1,id2);
+// cerr<<"erase duplicate "<<id1<<" "<<id2<<endl;
+ return;}
+ }
Real& un=geom->penetrationDepth;
phys->normalForce=phys->kn*std::max(un,(Real) 0)*geom->normal;
=== modified file 'pkg/common/Cylinder.hpp'
--- pkg/common/Cylinder.hpp 2010-12-13 12:11:43 +0000
+++ pkg/common/Cylinder.hpp 2010-12-15 10:11:41 +0000
@@ -52,10 +52,13 @@
public:
/// Emulate a sphere whose position is the projection of sphere's center on cylinder sphere, and with motion linearly interpolated between nodes
State fictiousState;
+// shared_ptr<Interaction> duplicate;
virtual ~CylScGeom ();
YADE_CLASS_BASE_DOC_ATTRS_CTOR(CylScGeom,ScGeom,"Geometry of a cylinder-sphere contact.",
((bool,onNode,false,,"contact on node?"))
+ ((int,isDuplicate,0,,"this flag is turned true (1) automaticaly if the contact is shared between two chained cylinders. A duplicated interaction will be skipped once by the constitutive law, so that only one contact at a time is effective. If isDuplicate=2, it means one of the two duplicates has no longer geometric interaction, and should be erased by the constitutive laws."))
+ ((int,trueInt,-1,,"Defines the body id of the cylinder where the contact is real, when :yref:`CylScGeom::isDuplicate`>0."))
((Vector3r,start,Vector3r::Zero(),,"position of 1st node |yupdate|"))
((Vector3r,end,Vector3r::Zero(),,"position of 2nd node |yupdate|"))
((Body::id_t,id3,0,,"id of next chained cylinder |yupdate|"))
@@ -73,6 +76,7 @@
static unsigned int currentChain;
vector<Body::id_t> barContacts;
vector<Body::id_t> nodeContacts;
+// shared_ptr<ChainedState> statePrev;
virtual ~ChainedState ();
void addToChain(Body::id_t bodyId) {
@@ -80,13 +84,17 @@
chainNumber=currentChain;
rank=chains[currentChain].size();
chains[currentChain].push_back(bodyId);
- bId=bodyId;}
+ bId=bodyId;
+// if (rank>0) statePrev = Body::byId(chains[chainNumber][rank-1],scene)->state;
+ }
void postLoad (ChainedState&){
if (bId<0) return;//state has not been chained yet
if (chains.size()<=currentChain) chains.resize(currentChain+1);
- if (chains[currentChain].size()<=rank) chains.resize(rank+1);
- chains[currentChain][rank]=bId;}
+ if (chains[currentChain].size()<=rank) chains[currentChain].resize(rank+1);
+ chains[currentChain][rank]=bId;
+// if (rank>0) statePrev = Body::byId(chains[chainNumber][rank-1],scene)->state;
+ }
YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ChainedState,State,"State of a chained bodies, containing information on connectivity in order to track contacts jumping over contiguous elements. Chains are 1D lists from which id of chained bodies are retrieved via :yref:rank<ChainedState::rank>` and :yref:chainNumber<ChainedState::chainNumber>`.",
((unsigned int,rank,0,,"rank in the chain"))
=== modified file 'pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp'
--- pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp 2010-11-12 08:03:16 +0000
+++ pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp 2010-12-15 10:11:41 +0000
@@ -17,7 +17,7 @@
const shared_ptr<Material>& b2,
const shared_ptr<Interaction>& interaction);
FUNCTOR2D(FrictMat,FrictMat);
- YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictMat_FrictPhys,IPhysFunctor,"Create a :yref:`FrictPhys` from two :yref:`FrictMats<FrictMat>`. The compliance of one sphere under symetric point loads is defined here as 1/(E.r), with E the stiffness of the sphere and r its radius, and corresponds to a compliance 1/(2.E.r)=1/(E.D) from each contact point. The compliance of the contact itself will be the sum of compliances from each sphere, i.e. 1/(E.D1)+1/(E.D2) in the general case, or 1/(E.r) in the special case of equal sizes. Note that summing compliances corresponds to an harmonic average of stiffnesss, which is how kn is actually computed in the :yref:`Ip2_FrictMat_FrictMat_FrictPhys` functor.\n\nThe shear stiffness ks of one sphere is defined via the material parameter :yref:`ElastMat::poisson`, as ks=poisson*kn, and the resulting shear stiffness of the interaction will be also an harmonic average.\n\nThe friction angle of the contact is defined as the minimum angle of the two materials in contact.\n\nOnly interactions with :yref:`ScGeom` or :yref:`Dem3DofGeom` geometry are meaningfully accepted; run-time typecheck can make this functor unnecessarily slow in general. Such design is problematic in itself, though -- from http://www.mail-archive.com/yade-dev@xxxxxxxxxxxxxxxxxxx/msg02603.html:\n\n\tYou have to suppose some exact type of IGeom in the Ip2 functor, but you don't know anything about it (Ip2 only guarantees you get certain IPhys types, via the dispatch mechanism).\n\n\tThat means, unless you use Ig2 functor producing the desired type, the code will break (crash or whatever). The right behavior would be either to accept any type (what we have now, at least in principle), or really enforce IGeom type of the interation passed to that particular Ip2 functor.",
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictMat_FrictPhys,IPhysFunctor,"Create a :yref:`FrictPhys` from two :yref:`FrictMats<FrictMat>`. The compliance of one sphere under symetric point loads is defined here as 1/(E.r), with E the stiffness of the sphere and r its radius, and corresponds to a compliance 1/(2.E.r)=1/(E.D) from each contact point. The compliance of the contact itself will be the sum of compliances from each sphere, i.e. 1/(E.D1)+1/(E.D2) in the general case, or 1/(E.r) in the special case of equal sizes. Note that summing compliances corresponds to an harmonic average of stiffnesss, which is how kn is actually computed in the :yref:`Ip2_FrictMat_FrictMat_FrictPhys` functor.\n\nThe shear stiffness ks of one sphere is defined via the material parameter :yref:`ElastMat::poisson`, as ks=poisson*kn, and the resulting shear stiffness of the interaction will be also an harmonic average.",
((shared_ptr<MatchMaker>,frictAngle,,,"Instance of :yref:`MatchMaker` determining how to compute interaction's friction angle. If ``None``, minimum value is used."))
);
};
=== modified file 'pkg/dem/ScGeom.cpp'
--- pkg/dem/ScGeom.cpp 2010-12-10 20:20:16 +0000
+++ pkg/dem/ScGeom.cpp 2010-12-15 10:11:41 +0000
@@ -38,7 +38,6 @@
Vector3r ScGeom::getIncidentVel(const State* rbp1, const State* rbp2, Real dt, const Vector3r& shift2, const Vector3r& shiftVel, bool avoidGranularRatcheting){
if(avoidGranularRatcheting){
- //FIXME : put the long comment on the wiki and keep only a small abstract and link here.
/* B.C. Comment :
Giving a short explanation of what we want to avoid :
Numerical ratcheting is best understood considering a small elastic cycle at a contact between two grains : assuming b1 is fixed, impose this displacement to b2 :
=== modified file 'scripts/test/chained-cylinder-roots.py'
--- scripts/test/chained-cylinder-roots.py 2010-12-10 11:25:22 +0000
+++ scripts/test/chained-cylinder-roots.py 2010-12-15 10:11:41 +0000
@@ -13,26 +13,24 @@
frictionAngle1=radians(15)
frictionAngle2=radians(15)
frictionAngle3=radians(15)
-O.dt=5e-5
+O.dt=5e-4
O.materials.append(FrictMat(young=4000000.0,poisson=0.5,frictionAngle=frictionAngle1,density=1600,label='spheremat'))
-O.materials.append(CohFrictMat(young=1.0e6,poisson=0.2,density=2.60e3,frictionAngle=frictionAngle2,label='walllmat'))
+O.materials.append(FrictMat(young=1.0e6,poisson=0.2,density=2.60e3,frictionAngle=frictionAngle2,label='walllmat'))
-Ns=10
+Ns=90
sp=pack.SpherePack()
#cohesive spheres crash because sphere-cylnder functor geneteragets ScGeom3D
#O.materials.append(CohFrictMat(young=1.0e5,poisson=0.03,density=2.60e2,frictionAngle=frictionAngle,label='spheremat'))
-if os.path.exists("cloud4cylinders"):
+if os.path.exists("cloud4cylinders"+`Ns`):
print "loading spheres from file"
- sp.load("cloud4cylinders")
- Ns=0
- for x in sp: Ns=Ns+1 #is there another way?
+ sp.load("cloud4cylinders"+`Ns`)
else:
print "generating spheres"
- Ns=sp.makeCloud(Vector3(-0.3,0.2,-1.0),Vector3(+0.3,+0.5,+1.0),-1,.2,Ns,False,0.9)
- sp.save("cloud4cylinders")
+ Ns=sp.makeCloud(Vector3(-0.3,0.2,-1.0),Vector3(+0.3,+0.5,+1.0),-1,.2,Ns,False,0.8)
+ sp.save("cloud4cylinders"+`Ns`)
O.bodies.append([utils.sphere(center,rad,material='spheremat') for center,rad in sp])
walls=utils.aabbWalls((Vector3(-0.3,-0.15,-1),Vector3(+0.3,+1.0,+1)),thickness=.1,material='walllmat')
@@ -55,7 +53,7 @@
## Apply gravity
GravityEngine(gravity=[1,-9.81,0],label='gravity'),
## Motion equation
- NewtonIntegrator(damping=0.20,label='newton'),
+ NewtonIntegrator(damping=0.05,label='newton'),
]
#Assemble cylinders in sinusoidal shapes
@@ -80,9 +78,4 @@
b.state.blockedDOFs=['x','y','z','rx','ry','rz']
ChainedState.currentChain=ChainedState.currentChain+1
-def outp(id=1):
- for i in O.interactions:
- if i.id1 == 1:
- print i.phys.shearForce
- print i.phys.normalForce
- return i
+O.saveTmp()
