Commit aa18a629 authored by Václav Šmilauer's avatar Václav Šmilauer
Browse files

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.
parent 5f7326b0
......@@ -655,7 +655,7 @@ and the XML looks like this:
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 @@ Expected parameters are indicated by macro name components separated with unders
};
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
^^^^^^^^^^^^^^^^^^^
......
#!/usr/bin/python
# -*- coding: utf-8 -*-
import matplotlib
matplotlib.use('TkAgg')
O.initializers=[
BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()])
......
......@@ -271,26 +271,13 @@ Erite plain text. Paragraphs are separated by empty lines (\n\n in c strings). U
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 @@ public :
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);
......
......@@ -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
......
......@@ -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 @@ class InteractionDispatchers: public GlobalEngine {
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 @@ class InteractionDispatchers: public GlobalEngine {
#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.");
......
......@@ -27,15 +27,25 @@ CREATE_LOGGER(InteractionGeometryDispatcher);
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;
}
}
......
......@@ -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);
}
}
......
......@@ -6,18 +6,24 @@ CREATE_LOGGER(StepDisplacer);
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;
}
}
}
......
......@@ -7,8 +7,10 @@ class StepDisplacer: public PartialEngine {
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;
};
......
......@@ -42,14 +42,20 @@ class Dem3DofGeom: public GenericSpheresContact{
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);
};
......
......@@ -11,8 +11,6 @@ YADE_PLUGIN((CFpmMat)(CFpmState)(CFpmPhys)(Ip2_CFpmMat_CFpmMat_CFpmPhys)(Law2_Sc
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();
......
......@@ -136,6 +136,8 @@ Real CpmPhys::computeViscoplScalingFactor(Real sigmaTNorm, Real sigmaTYield,Real
#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 @@ void Law2_Dem3DofGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom, sh
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 @@ void Law2_ScGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom, shared_
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 @@ void CpmStateUpdater::update(Scene* _scene){
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;
......
......@@ -503,6 +503,7 @@ shared_ptr<Interaction> Shop::createExplicitInteraction(body_id_t id1, body_id_t
assert(force && i);
if(!i) return i;
physMeta->explicitAction(b1->material,b2->material,i);
i->iterMadeReal=rb->currentIteration;
rb->interactions->insert(i);
return i;
}
......
......@@ -148,6 +148,7 @@ def createPlots():
# 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
......
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()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment