------------------------------------------------------------ revno: 2659 committer: Chiara Modenese <c.moden...@gmail.com> branch nick: yade timestamp: Mon 2011-01-17 18:21:10 +0000 message: - Add relAngVel function to ScGeom (a similar version can be called from python as done for incidentVel) @Bruno: I would like to add the incremental formulation to CFLaw, how do you suggest to do that? Maybe I can use a bool something like (useIncrementalForm)? modified: pkg/dem/HertzMindlin.cpp pkg/dem/ScGeom.cpp pkg/dem/ScGeom.hpp
-- 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/dem/HertzMindlin.cpp' --- pkg/dem/HertzMindlin.cpp 2011-01-13 16:46:20 +0000 +++ pkg/dem/HertzMindlin.cpp 2011-01-17 18:21:10 +0000 @@ -415,7 +415,9 @@ /********************************************/ if (includeMoment){ // new code to compute relative particle rotation (similar to the way the shear is computed) - Vector3r relAngVel = (b2->state->angVel-b1->state->angVel); + // use scg function to compute relAngVel + Vector3r relAngVel = scg->getRelAngVel(de1,de2,dt); + //Vector3r relAngVel = (b2->state->angVel-b1->state->angVel); relAngVel = relAngVel - scg->normal.dot(relAngVel)*scg->normal; // keep only the bending part Vector3r relRot = relAngVel*dt; // relative rotation due to rolling behaviour // incremental formulation for the bending moment (as for the shear part) === modified file 'pkg/dem/ScGeom.cpp' --- pkg/dem/ScGeom.cpp 2010-12-31 14:35:21 +0000 +++ pkg/dem/ScGeom.cpp 2011-01-17 18:21:10 +0000 @@ -79,6 +79,17 @@ avoidGranularRatcheting); } +Vector3r ScGeom::getRelAngVel(const State* rbp1, const State* rbp2, Real dt){ + Vector3r relAngVel = (rbp2->angVel-rbp1->angVel); + return relAngVel; +} + +Vector3r ScGeom::getRelAngVel_py(shared_ptr<Interaction> i){ + if(i->geom.get()!=this) throw invalid_argument("ScGeom object is not the same as Interaction.geom."); + Scene* scene=Omega::instance().getScene().get(); + return getRelAngVel(Body::byId(i->getId1(),scene)->state.get(),Body::byId(i->getId2(),scene)->state.get(),scene->dt); +} + //!Precompute relative rotations (and precompute ScGeom3D) void ScGeom6D::precomputeRotations(const State& rbp1, const State& rbp2, bool isNew, bool creep){ if (isNew) { === modified file 'pkg/dem/ScGeom.hpp' --- pkg/dem/ScGeom.hpp 2010-12-31 14:35:21 +0000 +++ pkg/dem/ScGeom.hpp 2011-01-17 18:21:10 +0000 @@ -43,9 +43,12 @@ Vector3r getIncidentVel(const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, const Vector3r& shift2, bool avoidGranularRatcheting=true); // Implement another version of getIncidentVel which does not handle periodicity. Vector3r getIncidentVel(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true); - + // Add function to get the relative angular velocity (useful to determine bending moment at the contact level) + Vector3r getRelAngVel(const State* rbp1, const State* rbp2, Real dt); + // convenience version to be called from python Vector3r getIncidentVel_py(shared_ptr<Interaction> i, bool avoidGranularRatcheting); + Vector3r getRelAngVel_py(shared_ptr<Interaction> i); YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom,GenericSpheresContact,"Class representing :yref:`geometry<IGeom>` of a contact point between two :yref:`bodies<Body>`s. It is more general than sphere-sphere contact : even though it is primarily focused on spheres interactions (reason for the 'Sc' namming), it is also used for representing contacts of a :yref:`Sphere` with a non-spherical bodies (:yref:`Facet`, :yref:`Plane`, :yref:`Box`, :yref:`ChainedCylinder`), or between non-spherical bodies. The contact has 3 DOFs (normal and 2Ãshear) and uses incremental algorithm for updating shear.\n\nWe use symbols $\\vec{x}$, $\\vec{v}$, $\\vec{\\omega}$ respectively for position, linear and angular velocities (all in global coordinates) and $r$ for particles radii; subscripted with 1 or 2 to distinguish 2 spheres in contact. Then we compute unit contact normal\n\n.. math::\n\n\t\\vec{n}=\\frac{\\vec{x}_2-\\vec{x}_1}{||\\vec{x}_2-\\vec{x}_1||}\n\nRelative velocity of spheres is then\n\n.. math::\n\n\t\\vec{v}_{12}=(\\vec{v}_2+\\vec{\\omega}_2\\times(-r_2\\vec{n}))-(\\vec{v}_1+\\vec{\\omega}_1\\times(r_1\\vec{n}))\n\nand its shear component\n\n.. math::\n\n\t\\Delta\\vec{v}_{12}^s=\\vec{v}_{12}-(\\vec{n}\\cdot\\vec{v}_{12})\\vec{n}.\n\nTangential displacement increment over last step then reads\n\n.. math::\n\n\t\\vec{x}_{12}^s=\\Delta t \\vec{v}_{12}^s.", ((Real,penetrationDepth,NaN,(Attr::noSave|Attr::readonly),"Penetration distance of spheres (positive if overlapping)")) @@ -54,6 +57,7 @@ /* extra initializers */ ((radius1,GenericSpheresContact::refR1)) ((radius2,GenericSpheresContact::refR2)), /* ctor */ createIndex(); twist_axis=orthonormal_axis=Vector3r::Zero();, /* py */ .def("incidentVel",&ScGeom::getIncidentVel_py,(python::arg("i"),python::arg("avoidGranularRatcheting")=true),"Return incident velocity of the interaction.") + .def("relAngVel",&ScGeom::getRelAngVel_py,(python::arg("i")),"Return relative angular velocity of the interaction.") ); REGISTER_CLASS_INDEX(ScGeom,GenericSpheresContact); };
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp