yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05432
[Branch ~yade-dev/yade/trunk] Rev 2375: 1. Add the possibility of custom constructors to python (not a hack as it used to be)
------------------------------------------------------------
revno: 2375
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Sun 2010-07-18 13:59:57 +0200
message:
1. Add the possibility of custom constructors to python (not a hack as it used to be)
2. Add cached functors to explicitly created interactions (not saved, though)
3. Improve the ScGeom-based Cpm law
4. Change interface of StepDisplacer to specify mov and rot separately (old interface with deltaSe3 works, but warns)
5. Add scripts/test/cpm-dem3dof-scgeom.py to show some differences between the 2 geometry formulations.
added:
scripts/test/cpm-dem3dof-scgeom.py
modified:
doc/sphinx/prog.rst
examples/simple-scene/simple-scene-plot.py
lib/serialization/Serializable.hpp
pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp
pkg/common/Engine/Dispatcher/InteractionDispatchers.hpp
pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp
pkg/common/Engine/Dispatcher/InteractionPhysicsDispatcher.cpp
pkg/common/Engine/PartialEngine/StepDisplacer.cpp
pkg/common/Engine/PartialEngine/StepDisplacer.hpp
pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp
pkg/dem/meta/ConcretePM.cpp
pkg/dem/meta/Shop.cpp
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/sphinx/prog.rst'
--- doc/sphinx/prog.rst 2010-07-14 19:13:18 +0000
+++ doc/sphinx/prog.rst 2010-07-18 11:59:57 +0000
@@ -655,7 +655,7 @@
Python attribute access
^^^^^^^^^^^^^^^^^^^^^^^^
-The macro :ref:`YADE_CLASS_BASE_DOC` introduced above is (behind the scenes) also used to create functions for accessing attributes from Python. As already noted, set of serialized attributes and set of attributes accessible from Python are identical. Besides attribute access and the ``[]`` operator access, these wrapper classes imitate also other functionality of regular python dictionaries:
+The macro :ref:`YADE_CLASS_BASE_DOC` introduced above is (behind the scenes) also used to create functions for accessing attributes from Python. As already noted, set of serialized attributes and set of attributes accessible from Python are identical. Besides attribute access, these wrapper classes imitate also some functionality of regular python dictionaries:
.. ipython::
@@ -814,6 +814,16 @@
};
+Special python constructors
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+The Python wrapper automatically create constructor that takes keyword (named) arguments corresponding to instance attributes; those attributes are set to values provided in the constructor. In some cases, more flexibility is desired (such as :yref:`InteractionDispatchers`, which takes 3 lists of functors). For such cases, you can override the function ``Serializable::pyHandleCustomCtorArgs``, which can arbitrarily modify the new (already existing) instance. It should modify in-place arguments given to it, as they will be passed further down to the routine which sets attribute values. In such cases, you should document the constructor::
+
+ .. admonition:: Special constructor
+
+ Constructs from lists of …
+
+which then appears in the documentation similar to :yref:`InteractionDispatchers`.
+
Static attributes
^^^^^^^^^^^^^^^^^^^
=== modified file 'examples/simple-scene/simple-scene-plot.py'
--- examples/simple-scene/simple-scene-plot.py 2010-07-06 11:31:29 +0000
+++ examples/simple-scene/simple-scene-plot.py 2010-07-18 11:59:57 +0000
@@ -1,5 +1,7 @@
#!/usr/bin/python
# -*- coding: utf-8 -*-
+import matplotlib
+matplotlib.use('TkAgg')
O.initializers=[
BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()])
=== modified file 'lib/serialization/Serializable.hpp'
--- lib/serialization/Serializable.hpp 2010-07-14 19:13:18 +0000
+++ lib/serialization/Serializable.hpp 2010-07-18 11:59:57 +0000
@@ -271,26 +271,13 @@
REGISTER_SERIALIZABLE_DESCRIPTOR(name,sname,SerializableTypes::CUSTOM_CLASS,isFundamental);
-
-
// helper functions
template <typename T>
-shared_ptr<T> Serializable_ctor_kwAttrs(const python::tuple& t, const python::dict& d){
- if(python::len(t)>1) throw runtime_error("Zero or one (and not more) non-keyword string argument required");
- string clss;
- if(python::len(t)==1){
- python::extract<string> clss_(t[0]); if(!clss_.check()) throw runtime_error("First argument (if given) must be a string.");
- clss=clss_();
- cerr<<"WARN: Constructing class using the Serializable('"<<clss<<"') syntax is deprecated. Use directly "<<clss<<"() instead."<<endl;
- }
+shared_ptr<T> Serializable_ctor_kwAttrs(python::tuple& t, python::dict& d){
shared_ptr<T> instance;
- if(clss.empty()){ instance=shared_ptr<T>(new T); }
- else{
- shared_ptr<Factorable> instance0=ClassFactory::instance().createShared(clss);
- if(!instance0) throw runtime_error("Invalid class `"+clss+"' (not created by ClassFactory).");
- instance=dynamic_pointer_cast<T>(instance0);
- if(!instance) throw runtime_error("Invalid class `"+clss+"' (unable to cast to typeid `"+typeid(T).name()+"')");
- }
+ instance=shared_ptr<T>(new T);
+ instance->pyHandleCustomCtorArgs(t,d); // can change t and d in-place
+ if(python::len(t)>0) throw runtime_error("Zero (not "+lexical_cast<string>(python::len(t))+") non-keyword constructor arguments required [in Serializable_ctor_kwAttrs; Serializable::pyHandleCustomCtorArgs might had changed it after your call].");
instance->pyUpdateAttrs(d);
return instance;
}
@@ -350,6 +337,9 @@
virtual bool checkPyClassRegistersItself(const std::string& thisClassName) const;
// perform class registration; overridden in all classes
virtual void pyRegisterClass(boost::python::object _scope);
+ // perform any manipulation of arbitrary constructor arguments coming from python, manipulating them in-place;
+ // the remainder is passed to the Serializable_ctor_kwAttrs of the respective class (note: args must be empty)
+ virtual void pyHandleCustomCtorArgs(boost::python::tuple& args, boost::python::dict& kw){ return; }
//! update attributes from dictionary
void pyUpdateAttrs(const boost::python::dict& d);
=== modified file 'pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp'
--- pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2010-07-17 17:37:54 +0000
+++ pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2010-07-18 11:59:57 +0000
@@ -3,17 +3,22 @@
YADE_PLUGIN((InteractionDispatchers));
CREATE_LOGGER(InteractionDispatchers);
-//! Pseudo-ctor for InteractionDispatchers, using lists of functors (might be turned into real ctor, perhaps)
-shared_ptr<InteractionDispatchers> InteractionDispatchers_ctor_lists(const std::vector<shared_ptr<InteractionGeometryFunctor> >& gff, const std::vector<shared_ptr<InteractionPhysicsFunctor> >& pff, const std::vector<shared_ptr<LawFunctor> >& cff /*, const std::vector<shared_ptr<IntrCallback> >& cbs*/){
- shared_ptr<InteractionDispatchers> instance(new InteractionDispatchers);
- FOREACH(shared_ptr<InteractionGeometryFunctor> gf, gff) instance->geomDispatcher->add(gf);
- FOREACH(shared_ptr<InteractionPhysicsFunctor> pf, pff) instance->physDispatcher->add(pf);
- FOREACH(shared_ptr<LawFunctor> cf, cff) instance->lawDispatcher->add(cf);
- // FOREACH(shared_ptr<IntrCallback> cb, cbs) instance->callbacks.push_back(cb);
- return instance;
+void InteractionDispatchers::pyHandleCustomCtorArgs(python::tuple& t, python::dict& d){
+ if(python::len(t)==0) return; // nothing to do
+ if(python::len(t)!=3) throw invalid_argument("Exactly 3 lists of functors must be given");
+ // parse custom arguments (3 lists) and do in-place modification of args
+ typedef std::vector<shared_ptr<InteractionGeometryFunctor> > vecGeom;
+ typedef std::vector<shared_ptr<InteractionPhysicsFunctor> > vecPhys;
+ typedef std::vector<shared_ptr<LawFunctor> > vecLaw;
+ vecGeom vg=python::extract<vecGeom>(t[0])();
+ vecPhys vp=python::extract<vecPhys>(t[1])();
+ vecLaw vl=python::extract<vecLaw>(t[2])();
+ FOREACH(shared_ptr<InteractionGeometryFunctor> gf, vg) this->geomDispatcher->add(gf);
+ FOREACH(shared_ptr<InteractionPhysicsFunctor> pf, vp) this->physDispatcher->add(pf);
+ FOREACH(shared_ptr<LawFunctor> cf, vl) this->lawDispatcher->add(cf);
+ t=python::tuple(); // empty the args; not sure if this is OK, as there is some refcounting in raw_constructor code
}
-
// #define IDISP_TIMING
#ifdef IDISP_TIMING
=== modified file 'pkg/common/Engine/Dispatcher/InteractionDispatchers.hpp'
--- pkg/common/Engine/Dispatcher/InteractionDispatchers.hpp 2010-05-06 10:11:20 +0000
+++ pkg/common/Engine/Dispatcher/InteractionDispatchers.hpp 2010-07-18 11:59:57 +0000
@@ -6,9 +6,6 @@
#include<yade/pkg-common/InteractionPhysicsDispatcher.hpp>
#include<yade/pkg-common/LawDispatcher.hpp>
-class InteractionDispatchers;
-shared_ptr<InteractionDispatchers> InteractionDispatchers_ctor_lists(const std::vector<shared_ptr<InteractionGeometryFunctor> >& gff, const std::vector<shared_ptr<InteractionPhysicsFunctor> >& pff, const std::vector<shared_ptr<LawFunctor> >& cff /*, const std::vector<shared_ptr<IntrCallback> >& cb=std::vector<shared_ptr<IntrCallback> >()*/ );
-
class InteractionDispatchers: public GlobalEngine {
bool alreadyWarnedNoCollider;
typedef std::pair<body_id_t, body_id_t> idPair;
@@ -21,8 +18,9 @@
void eraseAfterLoop(body_id_t id1,body_id_t id2){ eraseAfterLoopIds.push_back(idPair(id1,id2)); }
#endif
public:
+ virtual void pyHandleCustomCtorArgs(python::tuple& t, python::dict& d);
virtual void action();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(InteractionDispatchers,GlobalEngine,"Unified dispatcher for handling interaction loop at every step, for parallel performance reasons.",
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(InteractionDispatchers,GlobalEngine,"Unified dispatcher for handling interaction loop at every step, for parallel performance reasons.\n\n.. admonition:: Special constructor\n\n\tConstructs from 3 lists of :yref:`Ig2<InteractionGeometryFunctor>`, :yref:`Ip2<InteractionPhysicsFunctor>`, :yref:`Law<LawFunctor>` functors respectively; they will be passed to interal dispatchers, which you might retrieve. (NOT YET DONE: Optionally, list of :yref:`IntrCallbacks<IntrCallback>` can be provided as fourth argument.)",
((shared_ptr<InteractionGeometryDispatcher>,geomDispatcher,new InteractionGeometryDispatcher,"[will be overridden]"))
((shared_ptr<InteractionPhysicsDispatcher>,physDispatcher,new InteractionPhysicsDispatcher,"[will be overridden]"))
((shared_ptr<LawDispatcher>,lawDispatcher,new LawDispatcher,"[will be overridden]"))
@@ -37,9 +35,6 @@
#endif
,
/*py*/
- .def("__init__",python::make_constructor(InteractionDispatchers_ctor_lists),
- // (python::arg("geomFunctors"),python::arg("physFunctors"),python::arg("lawFunctors") /*,python::arg("callbacks")=std::vector<shared_ptr<IntrCallback> >() */),
- "Construct from lists :yref:`Ig2<InteractionGeometryFunctor>`, :yref:`Ip2<InteractionPhysicsFunctor>`, :yref:`Law<LawFunctor>` functors respectively; they will be passed to interal dispatchers, which you might retrieve. (NOT YET DONE: Optionally, list of :yref:`IntrCallbacks<IntrCallback>` can be provided as fourth argument.)")
.def_readonly("geomDispatcher",&InteractionDispatchers::geomDispatcher,"InteractionGeometryDispatcher object that is used for dispatch.")
.def_readonly("physDispatcher",&InteractionDispatchers::physDispatcher,"InteractionPhysicsDispatcher object used for dispatch.")
.def_readonly("lawDispatcher",&InteractionDispatchers::lawDispatcher,"LawDispatcher object used for dispatch.");
=== modified file 'pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp'
--- pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp 2010-07-06 11:47:43 +0000
+++ pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp 2010-07-18 11:59:57 +0000
@@ -27,15 +27,25 @@
shared_ptr<Interaction> InteractionGeometryDispatcher::explicitAction(const shared_ptr<Body>& b1, const shared_ptr<Body>& b2, bool force){
updateScenePtr();
if(force){
+ #ifdef YADE_DEVIRT_FUNCTORS
+ throw logic_error("InteractionGeometryDispatcher::explicitAction not supported with the devirt-functors feature (yet)");
+ #endif
assert(b1->shape && b2->shape);
- shared_ptr<Interaction> i(new Interaction(b1->getId(),b2->getId()));
- bool op=operator()(b1->shape,b2->shape,*b1->state,*b2->state,/*shift2*/Vector3r::Zero(),/*force*/true,i);
- if(!op) throw runtime_error("InteractionGeometryDispatcher::explicitAction could not dispatch for given types ("+b1->shape->getClassName()+","+b2->shape->getClassName()+") or the dispatchee returned false (it was asked to force creation of InteractionGeometry; that would a bug).");
- return i;
+ shared_ptr<Interaction> I(new Interaction(b1->getId(),b2->getId()));
+ // FIXME: this code is more or less duplicated from InteractionDispatchers :-(
+ bool swap=false;
+ I->functorCache.geom=getFunctor2D(b1->shape,b2->shape,swap);
+ if(!I->functorCache.geom) throw invalid_argument("InteractionGeometryDispatcher::explicitAction could not dispatch for given types ("+b1->shape->getClassName()+","+b2->shape->getClassName()+").");
+ if(swap){I->swapOrder();}
+ const shared_ptr<Body>& b1=Body::byId(I->getId1(),scene);
+ const shared_ptr<Body>& b2=Body::byId(I->getId2(),scene);
+ bool succ=I->functorCache.geom->go(b1->shape,b2->shape,*b1->state,*b2->state,/*shift2*/Vector3r::Zero(),/*force*/true,I);
+ if(!succ) throw logic_error("Functor "+I->functorCache.geom->getClassName()+"::go returned false, even if asked to force InteractionGeometry creation. Please report bug.");
+ return I;
} else {
- shared_ptr<Interaction> interaction(new Interaction(b1->getId(),b2->getId()));
- b1->shape && b2->shape && operator()( b1->shape , b2->shape , *b1->state , *b2->state , Vector3r::Zero(), /*force*/ false, interaction );
- return interaction;
+ shared_ptr<Interaction> I(new Interaction(b1->getId(),b2->getId()));
+ b1->shape && b2->shape && operator()( b1->shape , b2->shape , *b1->state , *b2->state , Vector3r::Zero(), /*force*/ false, I);
+ return I;
}
}
=== modified file 'pkg/common/Engine/Dispatcher/InteractionPhysicsDispatcher.cpp'
--- pkg/common/Engine/Dispatcher/InteractionPhysicsDispatcher.cpp 2010-03-20 12:40:44 +0000
+++ pkg/common/Engine/Dispatcher/InteractionPhysicsDispatcher.cpp 2010-07-18 11:59:57 +0000
@@ -14,12 +14,17 @@
*
* The interaction must be real (needed?).
*/
-void InteractionPhysicsDispatcher::explicitAction(shared_ptr<Material>& pp1, shared_ptr<Material>& pp2, shared_ptr<Interaction>& i){
+void InteractionPhysicsDispatcher::explicitAction(shared_ptr<Material>& pp1, shared_ptr<Material>& pp2, shared_ptr<Interaction>& I){
// should we throw instead of asserting?
//assert(i->isReal());
updateScenePtr();
- if(!i->interactionGeometry) throw runtime_error(string(__FILE__)+": explicitAction received interaction without interactionGeometry.");
- operator()(pp1,pp2,i);
+ if(!I->interactionGeometry) throw invalid_argument(string(__FILE__)+": explicitAction received interaction without interactionGeometry.");
+ if(!I->functorCache.phys){
+ bool dummy;
+ I->functorCache.phys=getFunctor2D(pp1,pp2,dummy);
+ if(!I->functorCache.phys) throw invalid_argument("InteractionPhysicsDispatcher::explicitAction did not find a suitable dispatch for types "+pp1->getClassName()+" and "+pp2->getClassName());
+ I->functorCache.phys->go(pp1,pp2,I);
+ }
}
=== modified file 'pkg/common/Engine/PartialEngine/StepDisplacer.cpp'
--- pkg/common/Engine/PartialEngine/StepDisplacer.cpp 2010-06-08 22:25:00 +0000
+++ pkg/common/Engine/PartialEngine/StepDisplacer.cpp 2010-07-18 11:59:57 +0000
@@ -6,18 +6,24 @@
YADE_PLUGIN((StepDisplacer));
void StepDisplacer::action(){
+ if(!isnan(deltaSe3.position[0])){
+ LOG_WARN("Using StepDisplacer::deltaSe3 is deprecated, use StepDisplacer.mov and StepDisplacer.rot instead (setting automatically now).");
+ mov=deltaSe3.position;
+ rot=deltaSe3.orientation;
+ deltaSe3.position=Vector3r(NaN,NaN,NaN); // to detect next manual change of deltaSe3
+ }
FOREACH(body_id_t id, subscribedBodies){
const shared_ptr<Body>& b=Body::byId(id,scene);
if(setVelocities){
const Real& dt=scene->dt;
- b->state->vel=deltaSe3.position/dt;
- AngleAxisr aa(deltaSe3.orientation); aa.axis().normalize();
+ b->state->vel=mov/dt;
+ AngleAxisr aa(rot); aa.axis().normalize();
b->state->angVel=aa.axis()*aa.angle()/dt;
LOG_DEBUG("Angular velocity set to "<<aa.axis()*aa.angle()/dt<<". Axis="<<aa.axis()<<", angle="<<aa.angle());
}
if(!setVelocities || (setVelocities && !b->isDynamic())){
- b->state->pos+=deltaSe3.position;
- b->state->ori=deltaSe3.orientation*b->state->ori;
+ b->state->pos+=mov;
+ b->state->ori=rot*b->state->ori;
}
}
}
=== modified file 'pkg/common/Engine/PartialEngine/StepDisplacer.hpp'
--- pkg/common/Engine/PartialEngine/StepDisplacer.hpp 2010-06-08 21:22:07 +0000
+++ pkg/common/Engine/PartialEngine/StepDisplacer.hpp 2010-07-18 11:59:57 +0000
@@ -7,8 +7,10 @@
public:
virtual void action();
YADE_CLASS_BASE_DOC_ATTRS(StepDisplacer,PartialEngine,"Apply generalized displacement (displacement or rotation) stepwise on subscribed bodies.",
- ((Se3r,deltaSe3,Se3r(Vector3r::Zero(),Quaternionr::Identity()),"Difference of position/orientation that will be added. Position is added to current :yref:`State::pos` using vector addition, orientation to :yref:`State::ori` using quaternion multiplication (rotation composition)."))
- ((bool,setVelocities,false,"If true, velocity and angularVelocity are modified in such a way that over one iteration (dt), the body will have prescribed se3 jump. In this case, se3 itself is not updated for :yref:`dynamic<Body::dynamic>` bodies, since they would have the delta applied twice (here and in the :yref:`integrator<NewtonIntegrator>`). For non-dynamic bodies however, se3 *is* still updated."))
+ ((Se3r,deltaSe3,Se3r(Vector3r(NaN,NaN,NaN),Quaternionr::Identity()),"|ydeprec| Set moc/rot directly instead; keps only for backwards compat. If the 0th component of the vector is not NaN, then it was updated by the user, warning is issued and mov and rot are updated automatically."))
+ ((Vector3r,mov,Vector3r::Zero(),"Linear displacement step to be applied per iteration, by addition to :yref:`State.pos`."))
+ ((Quaternionr,rot,Quaternionr::Identity(),"Rotation step to be applied per iteration (via rotation composition with :yref:`State.ori`)."))
+ ((bool,setVelocities,false,"If true, velocity and angularVelocity are modified in such a way that over one iteration (dt), the body will have the prescribed jump. In this case, :yref:`position<State.pos>` and :yref:`orientation<State.ori>` itself is not updated for :yref:`dynamic<Body::dynamic>` bodies, since they would have the delta applied twice (here and in the :yref:`integrator<NewtonIntegrator>`). For non-dynamic bodies however, they *are* still updated."))
);
DECLARE_LOGGER;
};
=== modified file 'pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp'
--- pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp 2010-07-15 14:26:48 +0000
+++ pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp 2010-07-18 11:59:57 +0000
@@ -42,14 +42,20 @@
Vector3r strainT(){return displacementT()/refLength;}
Real slipToStrainTMax(Real strainTMax){return slipToDisplacementTMax(strainTMax*refLength)/refLength;}
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(Dem3DofGeom,GenericSpheresContact,
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Dem3DofGeom,GenericSpheresContact,
"Abstract base class for representing contact geometry of 2 elements that has 3 degrees of freedom: normal (1 component) and shear (Vector3r, but in plane perpendicular to the normal).",
((Real,refLength,,"some length used to convert displacements to strains. |ycomp|"))
((Vector3r,contactPoint,,"some reference point for the interaction (usually in the middle). |ycomp|"))
((bool,logCompression,false,"make strain go to -∞ for length going to zero (false by default)."))
((Se3r,se31,,"Copy of body #1 se3 (needed to compute torque from the contact, strains etc). |yupdate|"))
((Se3r,se32,,"Copy of body #2 se3. |yupdate|")),
- createIndex()
+ createIndex();,
+ .def("displacementN",&Dem3DofGeom::displacementN)
+ .def("displacementT",&Dem3DofGeom::displacementT)
+ .def("strainN",&Dem3DofGeom::strainN)
+ .def("strainT",&Dem3DofGeom::strainT)
+ .def("slipToDisplacementTMax",&Dem3DofGeom::slipToDisplacementTMax)
+ .def("slipToStrainTMax",&Dem3DofGeom::slipToStrainTMax)
);
REGISTER_CLASS_INDEX(Dem3DofGeom,InteractionGeometry);
};
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-07-17 17:37:54 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-07-18 11:59:57 +0000
@@ -11,8 +11,6 @@
CREATE_LOGGER(Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM);
void Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact){
- const Real& dt = scene->dt;
-
ScGeom* geom = static_cast<ScGeom*>(ig.get());
CFpmPhys* phys = static_cast<CFpmPhys*>(ip.get());
const int &id1 = contact->getId1();
=== modified file 'pkg/dem/meta/ConcretePM.cpp'
--- pkg/dem/meta/ConcretePM.cpp 2010-07-17 17:37:54 +0000
+++ pkg/dem/meta/ConcretePM.cpp 2010-07-18 11:59:57 +0000
@@ -136,6 +136,8 @@
#include"../../../../brefcom-mm.hh"
#endif
+#undef CPM_MATERIAL_MODEL
+
void Law2_Dem3DofGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I){
Dem3DofGeom* contGeom=static_cast<Dem3DofGeom*>(_geom.get());
@@ -193,7 +195,9 @@
sigmaT=G*epsT; // trial stress
Real yieldSigmaT=max((Real)0.,undamagedCohesion*(1-omega)-sigmaN*tanFrictionAngle); // Mohr-Coulomb law with damage
if(sigmaT.squaredNorm()>yieldSigmaT*yieldSigmaT){
- sigmaT*=yieldSigmaT/sigmaT.norm(); // stress return
+ Real scale=yieldSigmaT/sigmaT.norm();
+ sigmaT*=scale; // stress return
+ epsT*=scale;
epsPlSum+=yieldSigmaT*contGeom->slipToStrainTMax(yieldSigmaT/G); // adjust strain
}
relResidualStrength=isCohesive?(kappaD<epsCrackOnset?1.:(1-omega)*(kappaD)/epsCrackOnset):0;
@@ -253,7 +257,8 @@
Real& epsN(BC->epsN);
Vector3r& epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); const Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real& crossSection(BC->crossSection); const Real& omegaThreshold(Law2_ScGeom_CpmPhys_Cpm::omegaThreshold); const Real& epsCrackOnset(BC->epsCrackOnset); Real& relResidualStrength(BC->relResidualStrength); const Real& epsFracture(BC->epsFracture); const bool& neverDamage(BC->neverDamage); Real& omega(BC->omega); Real& sigmaN(BC->sigmaN); Vector3r& sigmaT(BC->sigmaT); Real& Fn(BC->Fn); Vector3r& Fs(BC->Fs); // for python access
const bool& isCohesive(BC->isCohesive);
- epsN=contGeom->penetrationDepth/BC->refLength;
+ // FIXME: penetrationDepth does not account for interactionDetectionFactor!
+ epsN=(pos2-pos1).norm()/BC->refLength-1;
// FIXME: sign?
epsT=contGeom->rotate(epsT);
@@ -389,7 +394,7 @@
shared_ptr<CpmPhys> phys=dynamic_pointer_cast<CpmPhys>(I->interactionPhysics);
if(!phys) continue;
const body_id_t id1=I->getId1(), id2=I->getId2();
- Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(I->interactionGeometry.get());
+ GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(I->interactionGeometry.get());
Vector3r normalStress=((1./phys->crossSection)*geom->normal.dot(phys->normalForce))*geom->normal;
bodyStats[id1].sigma+=normalStress; bodyStats[id2].sigma+=normalStress;
=== modified file 'pkg/dem/meta/Shop.cpp'
--- pkg/dem/meta/Shop.cpp 2010-07-17 17:37:54 +0000
+++ pkg/dem/meta/Shop.cpp 2010-07-18 11:59:57 +0000
@@ -503,6 +503,7 @@
assert(force && i);
if(!i) return i;
physMeta->explicitAction(b1->material,b2->material,i);
+ i->iterMadeReal=rb->currentIteration;
rb->interactions->insert(i);
return i;
}
=== modified file 'py/plot.py'
--- py/plot.py 2010-07-17 17:37:54 +0000
+++ py/plot.py 2010-07-18 11:59:57 +0000
@@ -148,6 +148,7 @@
# create the legend
pylab.legend([xlateLabel(_p[0]) for _p in plots_p_y1],loc=('upper left' if len(plots_p_y2)>0 else 'best'))
pylab.ylabel(','.join([xlateLabel(_p[0]) for _p in plots_p_y1]))
+ pylab.xlabel(xlateLabel(p))
# create y2 lines, if any
if len(plots_p_y2)>0:
# try to move in the color palette a little further (magenta is 5th): r,g,b,c,m,y,k
=== added file 'scripts/test/cpm-dem3dof-scgeom.py'
--- scripts/test/cpm-dem3dof-scgeom.py 1970-01-01 00:00:00 +0000
+++ scripts/test/cpm-dem3dof-scgeom.py 2010-07-18 11:59:57 +0000
@@ -0,0 +1,85 @@
+from yade import plot
+
+# setup 2 interactions on 2 otherwise identical couples of spheres
+# one is handled by Law2_Dem3DofGeom_CpmPhys_Cpm and the other by Law2_ScGeom_CpmPhys_Cpm
+
+# move the second sphere tangentially or rotate it, pick [0] or [1]
+mode=['mov','rot'][0]
+# number of steps to do (some influence on the incremental computation)
+nSteps=100
+
+
+r1,r2=1e-3,1e-3
+dist=r1+r2
+offset=Vector3(0,0,2*dist)
+p1,p2=Vector3(0,0,0),Vector3(dist,0,0)
+O.materials.append(CpmMat(young=30e9,poisson=.2,frictionAngle=atan(.8),sigmaT=3e6,relDuctility=5,epsCrackOnset=1e-4,G_over_E=.2,neverDamage=False,plTau=-1,plRateExp=0,dmgTau=-1,dmgRateExp=0))
+# first 2 spheres used for Dem3DofGeom
+# the other 2 used for ScGeom (#3 is dynamic, since ScGeom needs that)
+O.bodies.append([utils.sphere(p1,r1,dynamic=False),utils.sphere(p2,r2,dynamic=False)])
+O.bodies.append([utils.sphere(p1+offset,r1,dynamic=False),utils.sphere(p2+offset,r2,dynamic=True)])
+
+
+O.engines=[InteractionGeometryDispatcher([Ig2_Sphere_Sphere_Dem3DofGeom()]),InteractionPhysicsDispatcher([Ip2_CpmMat_CpmMat_CpmPhys()])]
+i1=utils.createInteraction(0,1) # caches functors, no need to specify them in the main loop
+O.engines=[InteractionGeometryDispatcher([Ig2_Sphere_Sphere_ScGeom()]),InteractionPhysicsDispatcher([Ip2_CpmMat_CpmMat_CpmPhys()])]
+i2=utils.createInteraction(2,3)
+
+O.engines=[
+ InteractionDispatchers([],[],[Law2_ScGeom_CpmPhys_Cpm(),Law2_Dem3DofGeom_CpmPhys_Cpm(yieldSurfType=0)]),
+ StepDisplacer(subscribedBodies=[1,3],setVelocities=True,label='jumper'), # displace non-dynamic #1, set velocity on #3
+ NewtonIntegrator(damping=0),
+ PeriodicPythonRunner(iterPeriod=1,initRun=True,command='plotData()'),
+]
+
+def plotData():
+ allData={}
+ # gather same data for both configurations, suffix their labels with -DD/-Sc
+ for i,key,sphere in zip([i1,i2],['-DD','-Sc'],[1,3]):
+ data=dict(
+ zRot=O.bodies[sphere].state.ori.toAxisAngle()[1],
+ zShift=O.bodies[sphere].state.pos[2],
+ epsT=i.phys.epsT.norm() if key=='-Sc' else i.geom.strainT().norm(),
+ Ft=i.phys.shearForce.norm(),
+ epsN=i.phys.epsN,
+ epsPlSum=i.phys.epsPlSum,
+ relResStr=i.phys.relResidualStrength,
+ dist=(O.bodies[i.id1].state.pos-O.bodies[i.id2].state.pos).norm(),
+ sigmaT=i.phys.sigmaT.norm()
+ )
+ for k in data: allData[k+key]=data[k]
+ plot.addData(allData)
+
+if mode=='mov':
+ jumper.mov=Vector3(0,0,1/(1e4*nSteps))
+elif mode=='rot':
+ jumper.rot=Quaternion(Vector3.UnitZ,1/(1e4*nSteps))
+
+O.run(nSteps,True)
+if mode=='mov':
+ plot.plots={'zShift-DD':(
+ #('epsT-DD','g-'),('epsT-Sc','r^'),
+ ('dist-DD','g-'),('dist-Sc','r^'),
+ None,
+ ('sigmaT-DD','b-'),('sigmaT-Sc','m^')
+ #('relResStr-DD','b-'),('relResStr-Sc','m^')
+ #('epsN-DD','b-'),('epsN-Sc','m^')
+ )}
+elif mode=='rot':
+ plot.plots={'zRot-DD':(
+ ('epsT-DD','g|'),('epsT-Sc','r-'),
+ None,
+ ('sigmaT-DD','b-'),('sigmaT-Sc','mv')
+ )}
+
+
+
+if 1:
+ f='/tmp/cpm-geom-'+mode+'.pdf'
+ plot.plot(noShow=True).savefig(f)
+ print 'Plot saved to '+f
+ quit()
+else:
+ plot.plot()
+
+
