yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06730
[Branch ~yade-dev/yade/trunk] Rev 2659: - Add relAngVel function to ScGeom (a similar version can be called from python as done for incid...
------------------------------------------------------------
revno: 2659
committer: Chiara Modenese <c.modenese@xxxxxxxxx>
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);
};
Follow ups
