yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06440
[Branch ~yade-dev/yade/trunk] Rev 2590: - Cylinders : handle cylinder-sphere contacts with a fictious state representing interpolated mot...
------------------------------------------------------------
revno: 2590
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Mon 2010-12-06 20:24:20 +0100
message:
- Cylinders : handle cylinder-sphere contacts with a fictious state representing interpolated motion between nodes and passed to ScGeom
- Cylinders : implement CundallStrack for interaction between 3 bodies (1 body + 2 nodes)
- Cylinders : a number of fixes - no more bugs it seems
- Cylinders : handle zero-sized element for terminating chains
- Cylinders : update scripts; add example script chained-cylinder-roots.py, with cylinders, spheres, and boxes
- Ig2_Box_Sphere_ScGeom : small class name fixes ("6D")
- utils.py : adapt cylinders linking
- ElasticContactLaw : don't use Body::byId when possible
- ScGeom6D : put quaternion normalization between #ifdef guards (debug only)
added:
scripts/test/chained-cylinder-roots.py
modified:
pkg/common/Cylinder.cpp
pkg/common/Cylinder.hpp
pkg/dem/CohFrictMat.hpp
pkg/dem/ElasticContactLaw.cpp
pkg/dem/ElasticContactLaw.hpp
pkg/dem/Ig2_Box_Sphere_ScGeom.cpp
pkg/dem/Ig2_Box_Sphere_ScGeom.hpp
pkg/dem/ScGeom.cpp
py/tests/cohesive-chain.py
py/utils.py
scripts/test/chained-cylinder-spring.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-11-25 12:16:32 +0000
+++ pkg/common/Cylinder.cpp 2010-12-06 19:24:20 +0000
@@ -4,6 +4,7 @@
#include<yade/lib/opengl/OpenGLWrapper.hpp>
#endif
#include<yade/pkg/common/Aabb.hpp>
+#include<yade/pkg/dem/FrictPhys.hpp>
Cylinder::~Cylinder(){}
ChainedCylinder::~ChainedCylinder(){}
@@ -13,7 +14,7 @@
// Ig2_ChainedCylinder_ChainedCylinder_ScGeom6D::~Ig2_ChainedCylinder_ChainedCylinder_ScGeom6D() {}
-YADE_PLUGIN((Cylinder)(ChainedCylinder)(ChainedState)(CylScGeom)(Ig2_Sphere_ChainedCylinder_CylScGeom)(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)
#endif
@@ -45,32 +46,36 @@
ChainedCylinder *cylinder=YADE_CAST<ChainedCylinder*>(cm2.get());
Sphere *sphere=YADE_CAST<Sphere*>(cm1.get());
assert(sphereSt && cylinderSt && cylinder && sphere);
- if (cylinderSt->chains[cylinderSt->chainNumber].size()==(cylinderSt->rank+1)) {cerr << "last cylinder - ignored"<<endl; return false;}
-
-// int i1=cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1];
-// cerr << "i1 : "<<i1<< " i2: "<<cylinderSt->rank<<" pos: "<<Body::byId(cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1],scene)->state->pos<< " - "<<cylinderSt->pos<<endl;
+ 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;}
//FIXME : definition of segment in next line breaks periodicity
- Vector3r segment = Body::byId(cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1],scene)->state->pos-cylinderSt->pos;
- Vector3r branch = sphereSt->pos-cylinderSt->pos-shift2;
- bool isNew = !c->geom;
-
- //Check position of projection on cylinder axis
-// if ((segment.dot(branch)>(segment.dot(segment)/*+interactionDetectionFactor*cylinder->radius*/) && !c->isReal())) return (false);//position _after_ end of cylinder
- 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*/) && c->isReal())) cerr<<"pb1"<<endl;
-
- Real length = segment.norm();
- Vector3r direction = segment/length;
- Real dist = direction.dot(branch);
-// if ((dist<-interactionDetectionFactor*cylinder->radius) && !c->isReal()) return (false);//position _before_ start of cylinder
- if ((dist<-interactionDetectionFactor*cylinder->radius) && isNew) return (false);//position _before_ start of cylinder
+// 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;
+ 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
+ length = segment.norm();
+ direction = segment/length;
+ dist = direction.dot(branch);
+ if ((dist<-interactionDetectionFactor*cylinder->radius) && isNew) return false;
+ } 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;}
//Check sphere-cylinder distance
Vector3r projectedP = cylinderSt->pos+shift2 + direction*dist;
branch = projectedP-sphereSt->pos;
-// if ((branch.squaredNorm()>(pow(sphere->radius+cylinder->radius,2))) && !c->isReal()) return (false);
- if ((branch.squaredNorm()>(pow(sphere->radius+cylinder->radius,2))) && isNew) return (false);
+ 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);
@@ -78,47 +83,39 @@
scm->radius1=sphere->radius;
scm->radius2=cylinder->radius;
- scm->id3=cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1];
+ if (!isLast) 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?
- if (dist<0.0) {//We have sphere-node contact
+ if (dist<=0) {//We have sphere-node contact
Vector3r normal=(cylinderSt->pos+shift2)-sphereSt->pos;
Real norm=normal.norm();
+ normal /=norm;
+ scm->relPos=0;
scm->onNode=true; scm->relPos=0;
scm->penetrationDepth=sphere->radius+cylinder->radius-norm;
scm->contactPoint=sphereSt->pos+(sphere->radius-0.5*scm->penetrationDepth)*normal;
- scm->precompute(state1,state2,scene,c,normal/norm,isNew,shift2,true);//use sphere-sphere precompute (a node is a sphere)
+ scm->precompute(state1,state2,scene,c,normal,isNew,shift2,true);//use sphere-sphere precompute (a node is a sphere)
} else {//we have sphere-cylinder contact
- scm->onNode=false; scm->relPos=dist/length;
+ scm->onNode=false;
+ scm->relPos=dist/length;
Real norm=branch.norm();
Vector3r normal=branch/norm;
scm->penetrationDepth= sphere->radius+cylinder->radius-norm;
+
+ // define a virtual sphere at the projected center
+ scm->fictiousState.pos = projectedP;
+ scm->fictiousState.vel = (1-scm->relPos)*cylinderSt->vel + scm->relPos*cylinderNext->state->vel;
+ scm->fictiousState.angVel =
+ ((1-scm->relPos)*cylinderSt->angVel + scm->relPos*cylinderNext->state->angVel).dot(direction)*direction //twist part : interpolated
+ + segment.cross(cylinderNext->state->vel - cylinderSt->vel);// non-twist part : defined from nodes velocities
+
if (dist>length) {
scm->penetrationDepth=sphere->radius+cylinder->radius-(cylinderSt->pos+segment-sphereSt->pos).norm();
//FIXME : handle contact jump on next element
}
- scm->contactPoint = sphereSt->pos+scm->normal*(sphere->radius-0.5*scm->penetrationDepth);
-
- //FIXME : replace the block below with smart use of shift2=shift2 + cylinder_spin.cross(branch) - doesn't compile currently
-
- //precompute stuff manually (sphere-sphere precompute doesn't apply here if we want to avoid ratcheting
- /*
- if(!isNew) {
- scm->orthonormal_axis = scm->normal.cross(normal);
- Real angle = scene->dt*0.5*scm->normal.dot(sphereSt->angVel + cylinderSt->angVel);
- scm->twist_axis = angle*normal;}
- else scm->twist_axis=scm->orthonormal_axis=Vector3r::Zero();
- //Update contact normal
- scm->normal=normal;
- Vector3r c1x = sphere->radius*normal;
- Vector3r c2x = -cylinder->radius*normal;
- Vector3r relativeVelocity = (cylinderSt->vel+cylinderSt->angVel.cross(c2x)) - (sphereSt->vel+sphereSt->angVel.cross(c1x));
- //keep the shear part only
- relativeVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
- scm->shearInc = relativeVelocity*scene->dt;
- */
-
+ scm->contactPoint = sphereSt->pos+normal*(sphere->radius-0.5*scm->penetrationDepth);
+ scm->precompute(state1,scm->fictiousState,scene,c,branch/norm,isNew,shift2,false);//use sphere-sphere precompute (with a virtual sphere)
}
return true;
}
@@ -359,3 +356,61 @@
}
}
+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());
+ FrictPhys* phys = static_cast<FrictPhys*>(ip.get());
+ if(geom->penetrationDepth <0){
+ if (neverErase) {
+ phys->shearForce = Vector3r::Zero();
+ phys->normalForce = Vector3r::Zero();}
+ else scene->interactions->requestErase(id1,id2);
+ return;}
+ Real& un=geom->penetrationDepth;
+ phys->normalForce=phys->kn*std::max(un,(Real) 0)*geom->normal;
+
+ 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 (likely(!scene->trackEnergy)){//Update force but don't compute energy terms (see below))
+ // PFC3d SlipModel, is using friction angle. CoulombCriterion
+ if( shearForce.squaredNorm() > maxFs ){
+ Real ratio = 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
+ if(shearForce.squaredNorm() > maxFs){
+ Real ratio = 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;
+ Real dissip=((1/phys->ks)*(trialForce-shearForce))/*plastic disp*/ .dot(shearForce)/*active force*/;
+ if(dissip>0) scene->energy->add(dissip,"plastDissip",plastDissipIx,/*reset*/false);
+ }
+ // compute elastic energy as well
+ scene->energy->add(0.5*(phys->normalForce.squaredNorm()/phys->kn+phys->shearForce.squaredNorm()/phys->ks),"elastPotential",elastPotentialIx,/*reset at every timestep*/true);
+ }
+ if (!scene->isPeriodic) {
+ Vector3r force = -phys->normalForce-shearForce;
+ scene->forces.addForce(id1,force);
+ scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force));
+ //FIXME : include moment due to axis-contact distance in forces on node
+ Vector3r twist = (geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force);
+ scene->forces.addForce(id2,(geom->relPos-1)*force);
+ scene->forces.addTorque(id2,(1-geom->relPos)*twist);
+ if (geom->relPos) { //else we are on node (or on last node - and id3 is junk)
+ scene->forces.addForce(geom->id3,(-geom->relPos)*force);
+ scene->forces.addTorque(geom->id3,geom->relPos*twist);}
+ }
+// applyForceAtContactPoint(-phys->normalForce-shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
+ else {//FIXME : periodicity not implemented here :
+ Vector3r force = -phys->normalForce-shearForce;
+ scene->forces.addForce(id1,force);
+ scene->forces.addForce(id2,-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));
+ }
+}
=== modified file 'pkg/common/Cylinder.hpp'
--- pkg/common/Cylinder.hpp 2010-11-25 12:16:32 +0000
+++ pkg/common/Cylinder.hpp 2010-12-06 19:24:20 +0000
@@ -50,6 +50,9 @@
class CylScGeom: public ScGeom{
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;
+
virtual ~CylScGeom ();
YADE_CLASS_BASE_DOC_ATTRS_CTOR(CylScGeom,ScGeom,"Geometry of a cylinder-sphere contact.",
((bool,onNode,false,,"contact on node?"))
@@ -76,11 +79,19 @@
if (chains.size()<=currentChain) chains.resize(currentChain+1);
chainNumber=currentChain;
rank=chains[currentChain].size();
- chains[currentChain].push_back(bodyId);}
+ chains[currentChain].push_back(bodyId);
+ bId=bodyId;}
+
+ 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;}
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"))
((unsigned int,chainNumber,0,,"chain id"))
+ ((int,bId,-1,,"id of the body containing - for postLoad operations only"))
,
/* additional initializers */
/* ((pos,se3.position))
@@ -206,6 +217,29 @@
};
REGISTER_SERIALIZABLE(Bo1_ChainedCylinder_Aabb);
+
+class Law2_CylScGeom_FrictPhys_CundallStrack: public LawFunctor{
+ public:
+ //OpenMPAccumulator<Real> plasticDissipation;
+ virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ //Real elasticEnergy ();
+ //Real getPlasticDissipation();
+ //void initPlasticDissipation(Real initVal=0);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_CylScGeom_FrictPhys_CundallStrack,LawFunctor,"Law for linear compression, and Mohr-Coulomb plasticity surface without cohesion.\nThis law implements the classical linear elastic-plastic law from [CundallStrack1979]_ (see also [Pfc3dManual30]_). The normal force is (with the convention of positive tensile forces) $F_n=\\min(k_n u_n, 0)$. The shear force is $F_s=k_s u_s$, the plasticity condition defines the maximum value of the shear force : $F_s^{\\max}=F_n\\tan(\\phi)$, with $\\phi$ the friction angle.\n\n.. note::\n This law uses :yref:`ScGeom`; there is also functionally equivalent :yref:`Law2_Dem3DofGeom_FrictPhys_CundallStrack`, which uses :yref:`Dem3DofGeom` (sphere-box interactions are not implemented for the latest).\n\n.. note::\n This law is well tested in the context of triaxial simulation, and has been used for a number of published results (see e.g. [Scholtes2009b]_ and other papers from the same authors). It is generalised by :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment`, which adds cohesion and moments at contact.",
+ ((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
+ ((bool,traceEnergy,false,Attr::hidden,"Define the total energy dissipated in plastic slips at all contacts."))
+ ((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic dissipation (with O.trackEnergy)"))
+ ((int,elastPotentialIx,-1,(Attr::hidden|Attr::noSave),"Index for elastic potential energy (with O.trackEnergy)"))
+ ,,
+ //.def("elasticEnergy",&Law2_ScGeom_FrictPhys_CundallStrack::elasticEnergy,"Compute and return the total elastic energy in all \"FrictPhys\" contacts")
+ //.def("plasticDissipation",&Law2_ScGeom_FrictPhys_CundallStrack::getPlasticDissipation,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:`Law2_ScGeom_FrictPhys_CundallStrack::traceEnergy` is true.")
+ //.def("initPlasticDissipation",&Law2_ScGeom_FrictPhys_CundallStrack::initPlasticDissipation,"Initialize cummulated plastic dissipation to a value (0 by default).")
+ );
+ FUNCTOR2D(CylScGeom,FrictPhys);
+};
+REGISTER_SERIALIZABLE(Law2_CylScGeom_FrictPhys_CundallStrack);
+
+
// Keep this : Cylinders and ChainedCylinders will have different centers maybe.
// class Bo1_ChainedCylinder_Aabb : public Bo1_Cylinder_Aabb
// {
=== modified file 'pkg/dem/CohFrictMat.hpp'
--- pkg/dem/CohFrictMat.hpp 2010-12-02 19:27:56 +0000
+++ pkg/dem/CohFrictMat.hpp 2010-12-06 19:24:20 +0000
@@ -21,8 +21,8 @@
/// Serialization
YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohFrictMat,FrictMat,"",
((bool,isCohesive,true,,""))
- ((Real,normalCohesion,10000000,,""))
- ((Real,shearCohesion,10000000,,""))
+ ((Real,normalCohesion,0,,""))
+ ((Real,shearCohesion,0,,""))
((bool,momentRotationLaw,false,,"Use bending/twisting moment at contact. The contact will have moments only if both bodies have this flag true. See :yref:`CohFrictPhys::cohesionDisablesFriction` for details.")),
createIndex();
);
=== modified file 'pkg/dem/ElasticContactLaw.cpp'
--- pkg/dem/ElasticContactLaw.cpp 2010-12-06 11:54:38 +0000
+++ pkg/dem/ElasticContactLaw.cpp 2010-12-06 19:24:20 +0000
@@ -57,10 +57,7 @@
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=geom->penetrationDepth;
- TRVAR3(geom->penetrationDepth,de1->se3.position,de2->se3.position);
phys->normalForce=phys->kn*std::max(un,(Real) 0)*geom->normal;
Vector3r& shearForce = geom->rotate(phys->shearForce);
@@ -87,8 +84,10 @@
// compute elastic energy as well
scene->energy->add(0.5*(phys->normalForce.squaredNorm()/phys->kn+phys->shearForce.squaredNorm()/phys->ks),"elastPotential",elastPotentialIx,/*reset at every timestep*/true);
}
- if (!scene->isPeriodic)
- applyForceAtContactPoint(-phys->normalForce-shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
+ if (!scene->isPeriodic && !sphericalBodies) {
+ State* de1 = Body::byId(id1,scene)->state.get();
+ State* de2 = Body::byId(id2,scene)->state.get();
+ 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 = -phys->normalForce-shearForce;
scene->forces.addForce(id1,force);
=== modified file 'pkg/dem/ElasticContactLaw.hpp'
--- pkg/dem/ElasticContactLaw.hpp 2010-12-06 11:54:38 +0000
+++ pkg/dem/ElasticContactLaw.hpp 2010-12-06 19:24:20 +0000
@@ -23,6 +23,7 @@
void initPlasticDissipation(Real initVal=0);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_FrictPhys_CundallStrack,LawFunctor,"Law for linear compression, and Mohr-Coulomb plasticity surface without cohesion.\nThis law implements the classical linear elastic-plastic law from [CundallStrack1979]_ (see also [Pfc3dManual30]_). The normal force is (with the convention of positive tensile forces) $F_n=\\min(k_n u_n, 0)$. The shear force is $F_s=k_s u_s$, the plasticity condition defines the maximum value of the shear force : $F_s^{\\max}=F_n\\tan(\\phi)$, with $\\phi$ the friction angle.\n\n.. note::\n This law uses :yref:`ScGeom`; there is also functionally equivalent :yref:`Law2_Dem3DofGeom_FrictPhys_CundallStrack`, which uses :yref:`Dem3DofGeom` (sphere-box interactions are not implemented for the latest).\n\n.. note::\n This law is well tested in the context of triaxial simulation, and has been used for a number of published results (see e.g. [Scholtes2009b]_ and other papers from the same authors). It is generalised by :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment`, which adds cohesion and moments at contact.",
((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
+ ((bool,sphericalBodies,true,,"If true, compute branch vectors from radii (faster), else use contactPoint-position. Turning this flag true is safe for sphere-sphere contacts and a few other specific cases. It will give wrong values of torques on facets or boxes."))
((bool,traceEnergy,false,,"Define the total energy dissipated in plastic slips at all contacts. This will trace only plastic energy in this law, see O.trackEnergy for a more complete energies tracing"))
((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic dissipation (with O.trackEnergy)"))
((int,elastPotentialIx,-1,(Attr::hidden|Attr::noSave),"Index for elastic potential energy (with O.trackEnergy)"))
=== modified file 'pkg/dem/Ig2_Box_Sphere_ScGeom.cpp'
--- pkg/dem/Ig2_Box_Sphere_ScGeom.cpp 2010-11-07 11:46:20 +0000
+++ pkg/dem/Ig2_Box_Sphere_ScGeom.cpp 2010-12-06 19:24:20 +0000
@@ -41,18 +41,18 @@
Box* obb = static_cast<Box*>(cm1.get());
Sphere* s = static_cast<Sphere*>(cm2.get());
-
+
Vector3r extents = obb->extents;
// FIXME: do we need rotation matrix? Can't quaternion do just fine?
- Matrix3r boxAxisT=se31.orientation.toRotationMatrix();
+ Matrix3r boxAxisT=se31.orientation.toRotationMatrix();
Matrix3r boxAxis = boxAxisT.transpose();
-
+
Vector3r relPos21 = se32.position-se31.position; // relative position of centroids
-
+
// cOnBox_boxLocal is the sphere centroid (in box-local coordinates), but projected onto box if it is outside.
// _boxLocal means that ROTATION is local and origin is in box's origin
- Vector3r cOnBox_boxLocal=boxAxis*relPos21;
+ Vector3r cOnBox_boxLocal=boxAxis*relPos21;
if (cOnBox_boxLocal[0]<-extents[0]){cOnBox_boxLocal[0]=-extents[0]; inside=false; }
if (cOnBox_boxLocal[0]> extents[0]){cOnBox_boxLocal[0]= extents[0]; inside=false; }
@@ -60,7 +60,7 @@
if (cOnBox_boxLocal[1]> extents[1]){cOnBox_boxLocal[1]= extents[1]; inside=false; }
if (cOnBox_boxLocal[2]<-extents[2]){cOnBox_boxLocal[2]=-extents[2]; inside=false; }
if (cOnBox_boxLocal[2]> extents[2]){cOnBox_boxLocal[2]= extents[2]; inside=false; }
-
+
shared_ptr<ScGeom> scm;
if (inside){
// sphere center inside box. find largest `cOnBox_boxLocal' value:
@@ -68,14 +68,14 @@
// where cOnBox_boxLocal is minimal (i.e. where sphere center is closest perpendicularly to the box)
Real minCBoxDist=extents[0]-fabs(cOnBox_boxLocal[0]); int minCBoxDist_index=0;
for (int i=1; i<3; i++){Real tt=extents[i]-fabs(cOnBox_boxLocal[i]); if (tt<minCBoxDist){minCBoxDist=tt; minCBoxDist_index=i;}}
-
+
// contact normal aligned with box edge along largest `cOnBox_boxLocal' value
Vector3r normal_boxLocal = Vector3r(0,0,0);
normal_boxLocal[minCBoxDist_index]=(cOnBox_boxLocal[minCBoxDist_index]>0)?1.0:-1.0;
-
+
normal = boxAxisT*normal_boxLocal;
normal.normalize();
-
+
// se32 is sphere's se3
/*
*
@@ -99,7 +99,7 @@
bool isNew=!c->geom;
if (isNew) scm = shared_ptr<ScGeom>(new ScGeom());
else scm = YADE_PTR_CAST<ScGeom>(c->geom);
-
+
// contact point is in the middle of overlapping volumes
//(in the direction of penetration, which is normal to the box surface closest to sphere center) of overlapping volumes
scm->contactPoint = 0.5*(pt1+pt2);
@@ -127,10 +127,10 @@
* / | cOnBox_sphere
* | × |
* \ c ≡ se32->position
- * \ /
+ * \ /
* ~ /
* ^~~ ~ ~~^
- *
+ *
*/
pt1=cOnBox_box+se31.position;
@@ -138,10 +138,10 @@
cOnBox_sphere.normalize(); // we want only direction in the following
pt2=se32.position+cOnBox_sphere*s->radius;
-
+
bool isNew=!c->geom;
if (isNew) scm = shared_ptr<ScGeom>(new ScGeom());
- else scm = YADE_PTR_CAST<ScGeom>(c->geom);
+ else scm = YADE_PTR_CAST<ScGeom>(c->geom);
scm->contactPoint = 0.5*(pt1+pt2);
scm->penetrationDepth = depth;
scm->radius1 = s->radius;
@@ -172,7 +172,7 @@
#ifdef YADE_DEVIRT_FUNCTORS
bool Ig2_Box_Sphere_ScGeom6D::go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){ throw runtime_error("Do not call Ig2_Box_Sphere_ScGeom6D::go, use getStaticFunctorPtr and call that function instead."); }
-bool Ig2_Box_Sphere_ScGeom::goStatic(IGeomFunctor* self, const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){
+bool Ig2_Box_Sphere_ScGeom6D::goStatic(IGeomFunctor* self, const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){
#else
bool Ig2_Box_Sphere_ScGeom6D::go(
const shared_ptr<Shape>& cm1,
=== modified file 'pkg/dem/Ig2_Box_Sphere_ScGeom.hpp'
--- pkg/dem/Ig2_Box_Sphere_ScGeom.hpp 2010-11-07 11:46:20 +0000
+++ pkg/dem/Ig2_Box_Sphere_ScGeom.hpp 2010-12-06 19:24:20 +0000
@@ -64,7 +64,7 @@
const shared_ptr<Interaction>& c);
#ifdef YADE_DEVIRT_FUNCTORS
- void* getStaticFuncPtr(){ return (void*)&Ig2_Box_Sphere_ScGeom::goStatic; }
+ void* getStaticFuncPtr(){ return (void*)&Ig2_Box_Sphere_ScGeom6D::goStatic; }
static bool goStatic(IGeomFunctor* self, const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& se32, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
#endif
=== modified file 'pkg/dem/ScGeom.cpp'
--- pkg/dem/ScGeom.cpp 2010-11-17 11:25:26 +0000
+++ pkg/dem/ScGeom.cpp 2010-12-06 19:24:20 +0000
@@ -90,7 +90,9 @@
if (isNew) {
initRotations(rbp1,rbp2);
} else {
- rbp1.ori.normalize(); rbp2.ori.normalize(); initialOrientation2.normalize(); initialOrientation1.normalize();
+#ifdef YADE_SCGEOM_DEBUG
+ rbp1.ori.normalize(); rbp2.ori.normalize(); initialOrientation2.normalize(); initialOrientation1.normalize();
+#endif
Quaternionr delta((rbp1.ori * (initialOrientation1.conjugate())) * (initialOrientation2 * (rbp2.ori.conjugate())));
if (creep) delta = delta * twistCreep;
AngleAxisr aa(delta); // axis of rotation - this is the Moment direction UNIT vector; // angle represents the power of resistant ELASTIC moment
=== modified file 'py/tests/cohesive-chain.py'
--- py/tests/cohesive-chain.py 2010-11-15 17:19:37 +0000
+++ py/tests/cohesive-chain.py 2010-12-06 19:24:20 +0000
@@ -15,12 +15,12 @@
from math import *
class TestCohesiveChain(unittest.TestCase):
- # prefix test names with PBC:
+ # prefix test names with PBC:
def setUp(self):
O.reset();
young=1.0e3
poisson=5
- density=2.60e3
+ density=2.60e3
frictionAngle=radians(30)
O.materials.append(CohFrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,normalCohesion=1e13,shearCohesion=1e13,momentRotationLaw=True))
O.dt=1e-3
@@ -44,7 +44,7 @@
omeg=95.0/float(Ne); hy=0.05; hz=0.07;
px=float(i)*(omeg/60.0); py=sin(float(i)*omeg)*hy; pz=cos(float(i)*omeg)*hz;
px2=float(i+1.)*(omeg/60.0); py2=sin(float(i+1.)*omeg)*hy; pz2=cos(float(i+1.)*omeg)*hz;
- O.bodies.append(utils.chainedCylinder(begin=Vector3(pz,py,px), radius=0.005,end=Vector3(pz2,py2,px2),color=Vector3(0.6,0.5,0.5)))
+ utils.chainedCylinder(begin=Vector3(pz,py,px), radius=0.005,end=Vector3(pz2,py2,px2),color=Vector3(0.6,0.5,0.5))
O.bodies[Ne-1].state.blockedDOFs=['x','y','z','rx','ry','rz']
def testMotion(self):
"CohesiveChain: velocity/positions tested in transient dynamics and equilibrium state"
@@ -60,4 +60,3 @@
self.assertTrue(abs(tv1-v1)<abs(tolerance*tv1) and abs(tp1-p1)<abs(tolerance*tp1))
self.assertTrue(abs(tp2-p2)<abs(tolerance*tp2))
#self.assertTrue(abs(tv2-v2)<abs(tolerance*tv2) and abs(tp2-p2)<abs(tolerance*tp2)) #velocity comparison disabled, see comment above
-
\ No newline at end of file
=== modified file 'py/utils.py'
--- py/utils.py 2010-12-05 17:10:06 +0000
+++ py/utils.py 2010-12-06 19:24:20 +0000
@@ -77,7 +77,7 @@
def SpherePWaveTimeStep(radius,density,young):
r"""Compute P-wave critical timestep for a single (presumably representative) sphere, using formula for P-Wave propagation speed $\Delta t_{c}=\frac{r}{\sqrt{E/\rho}}$.
If you want to compute minimum critical timestep for all spheres in the simulation, use :yref:`yade.utils.PWaveTimeStep` instead.
-
+
>>> SpherePWaveTimeStep(1e-3,2400,30e9)
2.8284271247461903e-07
"""
@@ -184,7 +184,7 @@
>>> import random
>>> def matFactory(): return ElastMat(young=1e10*random.random(),density=1e3+1e3*random.random())
- ...
+ ...
>>> s3=utils.sphere([0,2,0],1,material=matFactory)
>>> s4=utils.sphere([1,2,0],1,material=matFactory)
@@ -204,7 +204,7 @@
"""Create box (cuboid) with given parameters.
:param Vector3 extents: half-sizes along x,y,z axes
-
+
See :yref:`yade.utils.sphere`'s documentation for meaning of other parameters."""
b=Body()
b.shape=Box(extents=extents,color=color if color else randomColor(),wire=wire,highlight=highlight)
@@ -219,14 +219,14 @@
def chainedCylinder(begin=Vector3(0,0,0),end=Vector3(1.,0.,0.),radius=0.2,dynamic=True,wire=False,color=None,highlight=False,material=-1,mask=1):
- """Create and chain a MinkCylinder with given parameters. This shape is the Minkowski sum of line and sphere.
+ """Create and connect a chainedCylinder with given parameters. The shape generated by repetaed calls of this function is the Minkowski sum of polyline and sphere.
:param Real radius: radius of sphere in the Minkowski sum.
:param Vector3 begin: first point positioning the line in the Minkowski sum
:param Vector3 last: last point positioning the line in the Minkowski sum
In order to build a correct chain, last point of element of rank N must correspond to first point of element of rank N+1 in the same chain (with some tolerance, since bounding boxes will be used to create connections.
-
+
:return: Body object with the :yref:`ChainedCylinder` :yref:`shape<Body.shape>`.
"""
segment=end-begin
@@ -236,7 +236,7 @@
V=2*(4./3)*math.pi*radius**3
geomInert=(2./5.)*V*radius**2+b.shape.length*b.shape.length*2*(4./3)*math.pi*radius**3
b.state=ChainedState()
- b.state.addToChain(b.id)
+ b.state.addToChain(O.bodies.append(b))
#b.state.blockedDOFs=['rx','ry','rz']
_commonBodySetup(b,V,Vector3(geomInert,geomInert,geomInert),material,resetState=False)
b.state.pos=b.state.refPos=begin
@@ -244,6 +244,7 @@
b.mask=mask
b.bound=Aabb(color=[0,1,0])
b.state.ori=b.state.ori.setFromTwoVectors(Vector3(0.,0.,1.),segment)
+ if (end == begin): b.state.ori = Quaternion((1,0,0),0)
return b
def wall(position,axis,sense=0,color=None,material=-1,mask=1):
@@ -275,7 +276,7 @@
:param bool wire: if ``True``, facets are shown as skeleton; otherwise facets are filled
:param bool noBound: set :yref:`Body.bounded`
:param Vector3-or-None color: color of the facet; random color will be assigned if ``None``.
-
+
See :yref:`yade.utils.sphere`'s documentation for meaning of other parameters."""
b=Body()
center=inscribedCircleCenter(vertices[0],vertices[1],vertices[2])
@@ -293,7 +294,7 @@
Create arbitrarily-aligned box composed of facets, with given center, extents and orientation.
If any of the box dimensions is zero, corresponding facets will not be created. The facets are oriented outwards from the box.
-
+
:Parameters:
`center`: Vector3
center of the created box
@@ -308,8 +309,8 @@
:Returns: list of facets forming the box.
"""
-
-
+
+
#Defense from zero dimensions
if (wallMask>63):
print "wallMask must be 63 or less"
@@ -323,7 +324,7 @@
if (((extents[0]==0) and (extents[1]==0)) or ((extents[0]==0) and (extents[2]==0)) or ((extents[1]==0) and (extents[2]==0))):
raise RuntimeError("Please, specify at least 2 none-zero dimensions in extents!");
# ___________________________
-
+
mn,mx=[-extents[i] for i in 0,1,2],[extents[i] for i in 0,1,2]
def doWall(a,b,c,d):
return [facet((a,b,c),**kw),facet((a,c,d),**kw)]
@@ -344,12 +345,12 @@
if wallMask&16: ret+=doWall(A,B,C,D)
if wallMask&32: ret+=doWall(E,H,G,F)
return ret
-
+
def facetCylinder(center,radius,height,orientation=Quaternion.Identity,segmentsNumber=10,wallMask=7,angleRange=None,closeGap=False,**kw):
"""
Create arbitrarily-aligned cylinder composed of facets, with given center, radius, height and orientation.
Return List of facets forming the cylinder;
-
+
:param Vector3 center: center of the created cylinder
:param float radius: cylinder radius
:param float height: cylinder height
@@ -360,7 +361,7 @@
:param bool closeGap: close range skipped in angleRange with triangular facets at cylinder bases.
:param **kw: (unused keyword arguments) passed to utils.facet;
"""
-
+
# check zero dimentions
if (segmentsNumber<3): raise RuntimeError("The segmentsNumber should be at least 3");
if (height<=0): raise RuntimeError("The height should have the positive value");
@@ -375,7 +376,7 @@
P1=[]; P2=[]
P1.append(Vector3(0,0,-height/2))
P2.append(Vector3(0,0,+height/2))
-
+
for i in anglesInRad:
X=radius*math.cos(i)
Y=radius*math.sin(i)
@@ -384,7 +385,7 @@
for i in range(0,len(P1)):
P1[i]=orientation*P1[i]+center
P2[i]=orientation*P2[i]+center
-
+
ret=[]
for i in range(2,len(P1)):
if wallMask&2:
@@ -401,10 +402,10 @@
pp=[(pts[0][0],pts[0][1],Z),(pts[1][0],pts[1][1],Z),(0,0,Z)]
pp=[orientation*p+center for p in pp]
ret.append(facet(pp,**kw))
-
+
return ret
-
-
+
+
def aabbWalls(extrema=None,thickness=None,oversizeFactor=1.5,**kw):
"""Return 6 boxes that will wrap existing packing as walls from all sides;
extrema are extremal points of the Aabb of the packing (will be calculated if not specified)
@@ -470,7 +471,7 @@
if b.dynamic or not onlyDynamic: b.shape.color=color
def avgNumInteractions(cutoff=0.,skipFree=False):
- r"""Return average numer of interactions per particle, also known as *coordination number* $Z$. This number is defined as
+ r"""Return average numer of interactions per particle, also known as *coordination number* $Z$. This number is defined as
.. math:: Z=2C/N
@@ -508,7 +509,7 @@
def plotDirections(aabb=(),mask=0,bins=20,numHist=True,noShow=False):
"""Plot 3 histograms for distribution of interaction directions, in yz,xz and xy planes and
(optional but default) histogram of number of interactions per body.
-
+
:returns: If *noShow* is ``False``, displays the figure and returns nothing. If *noShow*, the figure object is returned without being displayed (works the same way as :yref:`yade.plot.plot`).
"""
import pylab,math
@@ -616,7 +617,7 @@
def NormalRestitution2DampingRate(en):
r"""Compute the normal damping rate as a function of the normal coefficient of restitution $e_n$. For $e_n\in\langle0,1\rangle$ damping rate equals
-
+
.. math:: -\frac{\log e_n}{\sqrt{e_n^2+\pi^2}}
"""
@@ -652,7 +653,7 @@
class TableParamReader():
"""Class for reading simulation parameters from text file.
-
+
Each parameter is represented by one column, each parameter set by one line. Colums are separated by blanks (no quoting).
First non-empty line contains column titles (without quotes).
@@ -707,16 +708,16 @@
values[l]['description']=dd
descs.add(dd)
self.values=values
-
+
def paramDict(self):
"""Return dictionary containing data from file given to constructor. Keys are line numbers (which might be non-contiguous and refer to real line numbers that one can see in text editors), values are dictionaries mapping parameter names to their values given in the file. The special value '=' has already been interpreted, ``!`` (bangs) (if any) were already removed from column titles, ``description`` column has already been added (if absent)."""
return self.values
-
+
if __name__=="__main__":
tryTable="""head1 important2! !OMP_NUM_THREADS! abcd
1 1.1 1.2 1.3
'a' 'b' 'c' 'd' ### comment
-
+
# empty line
1 = = g
"""
@@ -726,7 +727,7 @@
f.close()
from pprint import *
pprint(TableParamReader(file).paramDict())
-
+
def runningInBatch():
'Tell whether we are running inside the batch or separately.'
import os
@@ -741,12 +742,12 @@
Read parameters from a file and assign them to __builtin__ variables.
The format of the file is as follows (commens starting with # and empty lines allowed)::
-
+
# commented lines allowed anywhere
name1 name2 … # first non-blank line are column headings
# empty line is OK, with or without comment
val1 val2 … # 1st parameter set
- val2 val2 … # 2nd
+ val2 val2 … # 2nd
…
Assigned tags:
@@ -804,7 +805,7 @@
defaults=[]
for k in kw.keys():
dictDefaults[k]=kw[k]
- defaults+=["%s=%s"%(k,kw[k])];
+ defaults+=["%s=%s"%(k,kw[k])];
O.tags['defaultParams']=",".join(defaults)
O.tags['params']=",".join(tagsParams)
dictParams.update(dictDefaults)
=== added file 'scripts/test/chained-cylinder-roots.py'
--- scripts/test/chained-cylinder-roots.py 1970-01-01 00:00:00 +0000
+++ scripts/test/chained-cylinder-roots.py 2010-12-06 19:24:20 +0000
@@ -0,0 +1,120 @@
+#--- bruno.chareyre@xxxxxxxxxxx ---
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+# Experiment beam-like behaviour with chained cylinders + CohFrict connexions
+
+from yade import utils,pack
+import math
+
+young=4.0e6
+poisson=3
+density=1e3
+frictionAngle1=radians(15)
+frictionAngle2=radians(15)
+frictionAngle3=radians(15)
+O.dt=5e-5
+
+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'))
+
+Ns=10
+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"):
+ print "loading spheres from file"
+ sp.load("cloud4cylinders")
+ Ns=0
+ for x in sp: Ns=Ns+1 #is there another way?
+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")
+
+O.bodies.append([utils.sphere(center,rad,material='spheremat') for center,rad in sp])
+
+#Ns=1
+#O.bodies.append(utils.sphere(Vector3(0.,0.4,0.0),0.1,material='spheremat'))
+
+#O.materials.append(FrictMat(young=150e6,poisson=.4,frictionAngle=.2,density=2600,label='wallmat'))#this one crash
+walls=utils.aabbWalls((Vector3(-0.3,-0.15,-1),Vector3(+0.3,+1.0,+1)),thickness=.1,material='walllmat')
+wallIds=O.bodies.append(walls)
+
+O.initializers=[
+ ## Create bounding boxes. They are needed to zoom the 3d view properly before we start the simulation.
+ BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_ChainedCylinder_Aabb(),Bo1_Box_Aabb()])
+]
+
+O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([
+ Bo1_ChainedCylinder_Aabb(),
+ Bo1_Sphere_Aabb(),
+ Bo1_Box_Aabb()
+ ]),
+ InteractionLoop(
+ [Ig2_ChainedCylinder_ChainedCylinder_ScGeom6D(), Ig2_Sphere_ChainedCylinder_CylScGeom(), Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
+ #[Ig2_ChainedCylinder_ChainedCylinder_ScGeom6D()],
+ [Ip2_2xCohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False,label='ipf'),Ip2_FrictMat_FrictMat_FrictPhys()],
+ [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(label='law'),Law2_ScGeom_FrictPhys_CundallStrack(),Law2_CylScGeom_FrictPhys_CundallStrack()]
+ ),
+ ## Apply gravity
+ GravityEngine(gravity=[1,-9.81,0],label='gravity'),
+ ## Motion equation
+ NewtonIntegrator(damping=0.10,label='newton'),
+ PyRunner(iterPeriod=500,command='history()'),
+ #PyRunner(iterPeriod=5000,command='if O.iter<21000 : yade.qt.center()')
+]
+
+#Generate a spiral
+O.materials.append(CohFrictMat(young=young,poisson=poisson,density=3.0*density,frictionAngle=frictionAngle3,normalCohesion=1e40,shearCohesion=1e40,momentRotationLaw=True,label='cylindermat'))
+Ne=30
+dy=0.03
+dx=0.2
+dz=0.2
+Nc=1 #nb. of additional chains
+for j in range(-Nc, Nc+1):
+ dyj = abs(float(j))*dy
+ dxj = abs(float(j))*dx
+ dzj = float(j)*dz
+ for i in range(0, Ne):
+ omega=20/float(Ne); hy=0.0; hz=0.05; hx=1.5;
+ px=float(i)*hx/float(Ne)-0.8+dxj; py=sin(float(i)*omega)*hy+dyj; pz=cos(float(i)*omega)*hz+dzj;
+ px2=float(i+1.)*hx/float(Ne)-0.8+dxj; py2=sin(float(i+1.)*omega)*hy+dyj; pz2=cos(float(i+1.)*omega)*hz+dzj;
+ utils.chainedCylinder(begin=Vector3(pz,py,px), radius=0.02,end=Vector3(pz2,py2,px2),color=Vector3(0.6,0.5,0.5),material='cylindermat')
+ if (i == Ne-1): #close the chain with a node of size 0
+ print "closing chain"
+ b=utils.chainedCylinder(begin=Vector3(pz2,py2,px2), radius=0.02,end=Vector3(pz2,py2,px2),color=Vector3(0.6,0.5,0.5),material='cylindermat')
+ 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
+#from yade import qt
+#qt.View()
+#qt.Controller()
+yade.qt.View();
+
+
+ #plot some results
+#from math import *
+from yade import plot
+plot.plots={'t':('pos1',None,'vel1')}
+def history():
+ plot.addData(pos1=O.bodies[6+Ns].state.pos[1], # potential elastic energy
+ vel1=O.bodies[6+Ns].state.vel[1],
+ t=O.time)
+
+O.run(1,True)
+#ipf.setCohesiononNewContacts=False
+#yade.qt.Renderer().bound=True
+#plot.plot()
+#O.bodies[0].state.angVel=Vector3(0.05,0,0)
+
=== modified file 'scripts/test/chained-cylinder-spring.py'
--- scripts/test/chained-cylinder-spring.py 2010-11-15 09:59:10 +0000
+++ scripts/test/chained-cylinder-spring.py 2010-12-06 19:24:20 +0000
@@ -37,7 +37,7 @@
omega=60.0/float(Ne); hy=0.10; hz=0.15;
px=float(i)*(omega/60.0); py=sin(float(i)*omega)*hy; pz=cos(float(i)*omega)*hz;
px2=float(i+1.)*(omega/60.0); py2=sin(float(i+1.)*omega)*hy; pz2=cos(float(i+1.)*omega)*hz;
- O.bodies.append(utils.chainedCylinder(begin=Vector3(pz,py,px), radius=0.005,end=Vector3(pz2,py2,px2),color=Vector3(0.6,0.5,0.5)))
+ utils.chainedCylinder(begin=Vector3(pz,py,px), radius=0.005,end=Vector3(pz2,py2,px2),color=Vector3(0.6,0.5,0.5))
def outp(id=1):
for i in O.interactions:
Follow ups
