yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05387
[Branch ~yade-dev/yade/trunk] Rev 2369: - Revert r2365 as a whole (and r2367).
------------------------------------------------------------
revno: 2369
committer: bchareyre <bchareyre@dt-rv020>
branch nick: trunk
timestamp: Thu 2010-07-15 16:26:48 +0200
message:
- Revert r2365 as a whole (and r2367).
- All changes in ScGeom and related functors will be commited again when I have time.
modified:
pkg/common/DataClass/Shape/Cylinder.cpp
pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.cpp
pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp
pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp
pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp
pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp
pkg/dem/Engine/Functor/Ig2_Facet_Sphere_ScGeom.cpp
pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp
pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp
pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp
pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.cpp
pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp
pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp
pkg/dem/meta/Shop.cpp
py/_eudoxos.cpp
py/_utils.cpp
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/DataClass/Shape/Cylinder.cpp'
--- pkg/common/DataClass/Shape/Cylinder.cpp 2010-07-14 19:02:58 +0000
+++ pkg/common/DataClass/Shape/Cylinder.cpp 2010-07-15 14:26:48 +0000
@@ -52,13 +52,18 @@
Real norm=normal.norm(); normal/=norm; // normal is unit vector now
Real penetrationDepth=s1->radius+s2->radius-norm;
scm->contactPoint=se31.position+(s1->radius-0.5*penetrationDepth)*normal;//0.5*(pt1+pt2);
-// if(isNew) scm->prevNormal=normal;
-// else scm->prevNormal=scm->normal;
-// scm->normal=normal;
+ if(isNew) scm->prevNormal=normal;
+ else scm->prevNormal=scm->normal;
+ scm->normal=normal;
scm->penetrationDepth=penetrationDepth;
scm->radius1=s1->radius;
scm->radius2=s2->radius;
- scm->precompute(state1,state2,scene,c,normal,isNew,true);
+#ifdef IGCACHE
+// if (scene->isPeriodic) {
+// Vector3r shiftVel = scene->cell->velGrad*scene->cell->Hsize*c->cellDist.cast<Real>();
+// scm->precompute(state1,state2,scene->dt,shiftVel,true);}
+ /*else */scm->precompute(state1,state2,scene,c,true);
+#endif
return true;
}
return false;
@@ -109,13 +114,15 @@
else { scm=shared_ptr<ScGeom>(new ScGeom()); c->interactionGeometry=scm; }
Real length=(bchain2.pos-bchain1.pos).norm();
Vector3r segt =pChain2->pos-pChain1->pos;
- if(isNew) {/*scm->normal=scm->prevNormal=segt/length;*/bs1->initLength=length;}
+ if(isNew) {scm->normal=scm->prevNormal=segt/length;bs1->initLength=length;}
+ else {scm->prevNormal=scm->normal; scm->normal=segt/length;}
scm->radius1=revert ? 0:bs1->initLength;
scm->radius2=revert ? bs1->initLength:0;
scm->penetrationDepth=bs1->initLength-length;
scm->contactPoint=bchain2.pos;
//bs1->segment used for fast BBs and projections + display
- bs1->segment= bchain2.pos-bchain1.pos;
+ if (revert) segt = -segt;
+ bs1->segment=segt;
#ifdef YADE_OPENGL
//bs1->length and s1->chainedOrientation used for display only,
bs1->length=length;
@@ -123,7 +130,7 @@
#endif
#ifdef IGCACHE
- scm->precompute(state1,state2,scene,c,segt/length,isNew,true);
+ scm->precompute(state1,state2,scene,c,true);
#else
cerr<<"this is not supposed to work without #define IGCACHE"<<endl;
#endif
=== modified file 'pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.cpp'
--- pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.cpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.cpp 2010-07-15 14:26:48 +0000
@@ -1,5 +1,5 @@
#include"DemXDofGeom.hpp"
-YADE_PLUGIN((Dem3DofGeom)/*(GenericSpheresContact)*/);
+YADE_PLUGIN((GenericSpheresContact)(Dem3DofGeom));
Real Dem3DofGeom::displacementN(){throw;}
Dem3DofGeom::~Dem3DofGeom(){}
-// GenericSpheresContact::~GenericSpheresContact(){}
+GenericSpheresContact::~GenericSpheresContact(){}
=== modified file 'pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp'
--- pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp 2010-07-15 14:26:48 +0000
@@ -1,36 +1,31 @@
// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
#pragma once
#include<yade/core/InteractionGeometry.hpp>
-#include<yade/pkg-dem/ScGeom.hpp>
/*! Abstract class that unites ScGeom and Dem3DofGeom,
created for the purposes of GlobalStiffnessTimeStepper.
It might be removed in the future. */
-// class GenericSpheresContact: public InteractionGeometry{
-// YADE_CLASS_BASE_DOC_ATTRS_CTOR(GenericSpheresContact,InteractionGeometry,
-// "Class uniting :yref:`ScGeom` and :yref:`Dem3DofGeom`, for the purposes of :yref:`GlobalStiffnessTimeStepper`. (It might be removed inthe future). Do not use this class directly.",
-// ((Vector3r,normal,,"Unit vector oriented along the interaction. |yupdate|"))
-// // ((Real,refR1,,"Reference radius of particle #1. |ycomp|"))
-// // ((Real,refR2,,"Reference radius of particle #2. |ycomp|")),
-// ,createIndex();
-// );
-// REGISTER_CLASS_INDEX(GenericSpheresContact,InteractionGeometry);
-//
-// virtual ~GenericSpheresContact(); // vtable
-// };
-// REGISTER_SERIALIZABLE(GenericSpheresContact);
+class GenericSpheresContact: public InteractionGeometry{
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(GenericSpheresContact,InteractionGeometry,
+ "Class uniting :yref:`ScGeom` and :yref:`Dem3DofGeom`, for the purposes of :yref:`GlobalStiffnessTimeStepper`. (It might be removed inthe future). Do not use this class directly.",
+ ((Vector3r,normal,,"Unit vector oriented along the interaction. |yupdate|"))
+ ((Real,refR1,,"Reference radius of particle #1. |ycomp|"))
+ ((Real,refR2,,"Reference radius of particle #2. |ycomp|")),
+ createIndex();
+ );
+ REGISTER_CLASS_INDEX(GenericSpheresContact,InteractionGeometry);
+
+ virtual ~GenericSpheresContact(); // vtable
+};
+REGISTER_SERIALIZABLE(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)
*/
-class ScGeom;
-
-class Dem3DofGeom: public ScGeom {
+class Dem3DofGeom: public GenericSpheresContact{
public:
virtual ~Dem3DofGeom();
- Real &refR1, &refR2;
-
// API that needs to be implemented in derived classes
virtual Real displacementN();
virtual Vector3r displacementT(){throw;}
@@ -47,17 +42,16 @@
Vector3r strainT(){return displacementT()/refLength;}
Real slipToStrainTMax(Real strainTMax){return slipToDisplacementTMax(strainTMax*refLength)/refLength;}
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(Dem3DofGeom,ScGeom,"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).",
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(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|"))
+ ((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|")),
- ((refR1,ScGeom::radius1))
- ((refR2,ScGeom::radius2)),
- createIndex();,
+ createIndex()
);
- REGISTER_CLASS_INDEX(Dem3DofGeom,ScGeom);
+ REGISTER_CLASS_INDEX(Dem3DofGeom,InteractionGeometry);
};
REGISTER_SERIALIZABLE(Dem3DofGeom);
=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2010-07-15 14:26:48 +0000
@@ -7,7 +7,10 @@
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
+
YADE_PLUGIN((ScGeom));
+
+// At least one virtual method must be in the .cpp file (!!!)
ScGeom::~ScGeom(){};
Vector3r& ScGeom::rotate(Vector3r& shearForce) const {
@@ -20,14 +23,12 @@
}
//!Precompute data needed for rotating tangent vectors attached to the interaction
-void ScGeom::precompute(const State& rbp1, const State& rbp2, const Scene* scene, const shared_ptr<Interaction>& c, const Vector3r& currentNormal, bool isNew, bool avoidGranularRatcheting){
- if(!isNew) {
- orthonormal_axis = normal.cross(currentNormal);
- Real angle = scene->dt*0.5*normal.dot(rbp1.angVel + rbp2.angVel);
- twist_axis = angle*normal;}
- else twist_axis=orthonormal_axis=Vector3r::Zero();
- //Update contact normal
- normal=currentNormal;
+void ScGeom::precompute(const State& rbp1, const State& rbp2, const Scene* scene, const shared_ptr<Interaction>& c, bool avoidGranularRatcheting){
+ orthonormal_axis = prevNormal.cross(normal);
+ Real angle = scene->dt*0.5*normal.dot(rbp1.angVel + rbp2.angVel);
+ twist_axis = angle*normal;
+ //The previous normal could be updated here after being used, them remove it from Ig2's?
+ //prevNormal=normal;
//Precompute shear increment
Vector3r& x = contactPoint;
Vector3r c1x, c2x;
@@ -63,7 +64,73 @@
shearInc = relativeVelocity*scene->dt;
}
-void ScGeom::updateShear(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){rotate(shear); shear -= shearInc;}
+
+//Do we need a function for such simple operation? BTW, the original was BAD duplicated code. ;-)
+Vector3r ScGeom::updateShear(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){rotate(shear); shear -= shearInc;}
+
+#ifndef IGCACHE
+Vector3r ScGeom::rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){
+ //Just pass null shift to the periodic version
+ return rotateAndGetShear(shearForce,prevNormal,rbp1,rbp2,dt,Vector3r::Zero(),avoidGranularRatcheting);
+}
+#endif
+
+//Removing this function will need to adapt all laws
+const Vector3r& ScGeom::rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting){
+#ifdef IGCACHE
+ rotate(shearForce);
+ return shearInc;
+ #else
+ rotate(shearForce,prevNormal,rbp1,rbp2,dt);
+ Vector3r& x = contactPoint;
+ Vector3r c1x, c2x;
+ if(avoidGranularRatcheting){
+ // FIXME: For sphere-facet contact this will give an erroneous value of relative velocity...
+ c1x = radius1*normal;
+ c2x = -radius2*normal;
+ } else {
+ // FIXME: It is correct for sphere-sphere and sphere-facet contact
+ c1x = (x - rbp1->pos);
+ c2x = (x - rbp2->pos);
+ }
+ Vector3r relativeVelocity = (rbp2->vel+rbp2->angVel.cross(c2x)) - (rbp1->vel+rbp1->angVel.cross(c1x));
+
+ //update relative velocity for interactions across periods
+ relativeVelocity+=shiftVel;
+ //keep the shear part only
+ relativeVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
+ Vector3r shearDisplacement = relativeVelocity*dt;
+ return shearDisplacement;
+#endif //IGCACHE
+}
+
+Vector3r& ScGeom::rotate(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt){
+#ifdef IGCACHE
+ return rotate(shearForce);
+#else
+ Vector3r axis;
+ Real angle;
+ // approximated rotations
+ axis = prevNormal.cross(normal);
+ shearForce -= shearForce.cross(axis);
+ angle = dt*0.5*normal.dot(rbp1->angVel + rbp2->angVel);
+ axis = angle*normal;
+ shearForce -= shearForce.cross(axis);
+ // exact rotations
+
+// Quaternionr q;
+// axis = prevNormal.Cross(normal);
+// angle = acos(normal.Dot(prevNormal));
+// q.FromAngleAxis(angle,axis);
+// shearForce = shearForce*q;
+// angle = dt*0.5*normal.dot(rbp1->angVel+rbp2->angVel);
+// axis = normal;
+// q.FromAngleAxis(angle,axis);
+// shearForce = q*shearForce;
+
+ return shearForce;
+#endif
+}
Vector3r ScGeom::getIncidentVel(const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting){
Vector3r& x = contactPoint;
@@ -87,3 +154,22 @@
//Just pass null shift to the periodic version
return getIncidentVel(rbp1,rbp2,dt,Vector3r::Zero(),avoidGranularRatcheting);
}
+
+/* keep this for reference; declarations in the header */
+#if 0
+ Vector3r ScGeom::relRotVector() const{
+ Quaternionr relOri12=ori1.Conjugate()*ori2;
+ Quaternionr oriDiff=initRelOri12.Conjugate()*relOri12;
+ Vector3r axis; Real angle;
+ oriDiff.ToAxisAngle(axis,angle);
+ if(angle>Mathr::PI)angle-=Mathr::TWO_PI;
+ return angle*axis;
+ }
+
+ void ScGeom::bendingTorsionAbs(Vector3r& bend, Real& tors){
+ Vector3r relRot=relRotVector();
+ tors=relRot.Dot(normal);
+ bend=relRot-tors*normal;
+ }
+#endif
+
=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp 2010-07-15 14:26:48 +0000
@@ -7,8 +7,8 @@
#include<yade/core/InteractionGeometry.hpp>
#include<yade/core/State.hpp>
#include<yade/lib-base/Math.hpp>
-
-/*! Class representing geometry of two bodies in contact and defining values needed for the incremental algorithm.
+#include<yade/pkg-dem/DemXDofGeom.hpp>
+/*! Class representing geometry of two bodies in contact.
*
* The code under SCG_SHEAR is experimental and is used only if ElasticContactLaw::useShear is explicitly true
*/
@@ -16,26 +16,35 @@
#define SCG_SHEAR
#define IGCACHE
-class ScGeom: public InteractionGeometry {
+class ScGeom: public GenericSpheresContact {
private:
//cached values
Vector3r twist_axis;//rotation vector arounf normal
Vector3r orthonormal_axis;//rotation vector in contact plane
Vector3r shearInc;//shear disp. increment
public:
+ // inherited from GenericSpheresContact: Vector3r& normal;
+ Real &radius1, &radius2;
Real penetrationDepth;
virtual ~ScGeom();
- //!precompute values of shear increment and interaction rotation data. Update contact normal to the vurrentNormal value. Precondition : the value of normal is not updated outside (and before) this function.
- void precompute(const State& rbp1, const State& rbp2, const Scene* scene, const shared_ptr<Interaction>& c, const Vector3r& currentNormal, bool isNew, bool avoidGranularRatcheting=true);
+ //!precompute values of shear increment and interaction rotation data
+ void precompute(const State& rbp1, const State& rbp2, const Scene* scene, const shared_ptr<Interaction>& c, bool avoidGranularRatcheting=true);
+ //! Periodic variant. Needs the velocity shift across periods for periodic BCs (else it is safe to pass Vector3r::Zero()). Typically obtained as scene->cell->velGrad*scene->cell->Hsize*cellDist.
+// void precompute(const State& rbp1, const State& rbp2, const Real& dt, const Vector3r& shiftVel, bool avoidGranularRatcheting=true);
//! Rotates a "shear" vector to keep track of contact orientation. Returns reference of the updated vector.
Vector3r& rotate(Vector3r& tangentVector) const;
const Vector3r& shearIncrement() const {return shearInc;}
//!________________REMOVED IN NEXT STEP___________________________________
+ Vector3r& rotate(Vector3r& tangentVector, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt);
//! update shear on this contact given dynamic parameters of bodies. Should be called from constitutive law, exactly once per iteration. Returns shear increment in this update, which is already added to shear.
- void updateShear(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
+ Vector3r updateShear(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
+ //! Update the shear force orientation (or any other vector, think "shear") on this contact given dynamic parameters of bodies. Compute and return shear displacement increment.
+ const Vector3r& rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt,bool avoidGranularRatcheting=true);
+ //! Periodic variant. Needs the velocity shift across periods for periodic BCs (else it is safe to pass Vector3r::Zero()). Typically obtained as scene->cell->velGrad*scene->cell->Hsize*cellDist. It would be better to define the shift transparently inside the function, but it needs scene and interaction pointers, which we don't have here.
+ const Vector3r& rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting=true);
//!________________END REMOVED IN NEXT STEP___________________________________
// Add method which returns the impact velocity (then, inside the contact law, this can be split into shear and normal component). Handle periodicity.
@@ -43,19 +52,17 @@
// Implement another version of getIncidentVel which does not handle periodicity.
Vector3r getIncidentVel(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
- YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(ScGeom,InteractionGeometry,"Class representing :yref:`geometry<InteractionGeometry>` of two :yref:`bodies<Body>` in contact. The contact has 3 DOFs (normal and 2×shear) and uses incremental algorithm for updating shear. (For shear formulated in total displacements and rotations, see :yref:`Dem3DofGeom` and related classes).\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.",
- ((Vector3r,normal,,"Unit vector oriented along the interaction. Updated by :yref:`ScGeom::precompute` |yupdate|"))
- ((Real,radius1,0,"Reference radius of particle #1. |ycomp|"))
- ((Real,radius2,0,"Reference radius of particle #2. |ycomp|"))
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom,GenericSpheresContact,"Class representing :yref:`geometry<InteractionGeometry>` of two :yref:`bodies<Body>` in contact. The contact has 3 DOFs (normal and 2×shear) and uses incremental algorithm for updating shear. (For shear formulated in total displacements and rotations, see :yref:`Dem3DofGeom` and related classes).\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.",
((Vector3r,contactPoint,Vector3r::Zero(),"Reference point of the contact. |ycomp|"))
- ((Vector3r,shear,Vector3r::Zero(),"Total value of the current shear. Update the value using ScGeom::updateShear."))
- /* ctor */
- ,createIndex(); shearInc=Vector3r::Zero(); twist_axis=orthonormal_axis=Vector3r::Zero();,
+ ((Vector3r,shear,Vector3r::Zero(),"Total value of the current shear. Update the value using ScGeom::updateShear. This is experimental. You have to update it yourself in constitutive laws using e.g. :yref:`ScGeom::rotate` and :yref:`ScGeom::shearIncrement`"))
+ ((Vector3r,prevNormal,Vector3r::Zero(),"Normal of the contact in the previous step. |ycomp|")) ,
+ /* extra initializers */ ((radius1,GenericSpheresContact::refR1)) ((radius2,GenericSpheresContact::refR2)),
+ /* ctor */ createIndex(); shearInc=Vector3r::Zero(); twist_axis=orthonormal_axis=Vector3r::Zero();,
/* py */
.def_readwrite("penetrationDepth",&ScGeom::penetrationDepth,"documentation")
.def_readonly("shearIncrement",&ScGeom::shearInc,"documentation")
);
- REGISTER_CLASS_INDEX(ScGeom,InteractionGeometry);
+ REGISTER_CLASS_INDEX(ScGeom,GenericSpheresContact);
};
REGISTER_SERIALIZABLE(ScGeom);
=== modified file 'pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2010-07-15 14:26:48 +0000
@@ -99,16 +99,18 @@
bool isNew=!c->interactionGeometry;
if (isNew) scm = shared_ptr<ScGeom>(new ScGeom());
else scm = YADE_PTR_CAST<ScGeom>(c->interactionGeometry);
+
+ if(isNew) { /* same as below */ scm->prevNormal=normal;}
+ else {scm->prevNormal=scm->normal;}
// 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);
-// scm->normal = normal;
+ scm->normal = normal;
scm->penetrationDepth = (pt1-pt2).norm();
scm->radius1 = s->radius;
scm->radius2 = s->radius;
c->interactionGeometry = scm;
- scm->precompute(state1,state2,scene,c,normal,isNew,true);
} else { // outside
Vector3r cOnBox_box = boxAxisT*cOnBox_boxLocal; // projection of sphere's center on closest box surface - relative to box's origin, but GLOBAL orientation!
Vector3r cOnBox_sphere = cOnBox_box-relPos21; // same, but origin in sphere's center
@@ -142,13 +144,23 @@
bool isNew=!c->interactionGeometry;
if (isNew) scm = shared_ptr<ScGeom>(new ScGeom());
else scm = YADE_PTR_CAST<ScGeom>(c->interactionGeometry);
+ if(isNew) { /* same as below */ scm->prevNormal=-cOnBox_sphere; }
+ else scm->prevNormal=scm->normal;
scm->contactPoint = 0.5*(pt1+pt2);
+ //scm->normal = pt1-pt2; scm->normal.normalize();
+ //scm->penetrationDepth = (pt1-pt2).norm();
+ scm->normal = -cOnBox_sphere;
scm->penetrationDepth = depth;
scm->radius1 = s->radius;
scm->radius2 = s->radius;
c->interactionGeometry = scm;
- scm->precompute(state1,state2,scene,c,-cOnBox_sphere,isNew,true);
}
+#ifdef IGCACHE
+// if (scene->isPeriodic){
+// Vector3r shiftVel = scene->cell->velGrad*scene->cell->Hsize*c->cellDist.cast<Real>();
+// scm->precompute(state1,state2,scene->dt,shiftVel,true);}
+ /*else */scm->precompute(state1,state2,scene,c,true);
+#endif
return true;
}
=== modified file 'pkg/dem/Engine/Functor/Ig2_Facet_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Facet_Sphere_ScGeom.cpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/Engine/Functor/Ig2_Facet_Sphere_ScGeom.cpp 2010-07-15 14:26:48 +0000
@@ -103,11 +103,19 @@
normal = facetAxisT*normal; // in global orientation
scm->contactPoint = se32.position - (sphereRadius-0.5*penetrationDepth)*normal;
+ //Update prevNormal (mandatory in ScGeom algoritms)
+ if(isNew) { scm->prevNormal=normal;}
+ else scm->prevNormal=scm->normal;
+ scm->normal = normal;
scm->penetrationDepth = penetrationDepth;
scm->radius1 = 2*sphereRadius;
scm->radius2 = sphereRadius;
- if (isNew) c->interactionGeometry = scm;
- scm->precompute(state1,state2,scene,c,normal,isNew,true);
+ if (!c->interactionGeometry)
+ c->interactionGeometry = scm;
+// if (scene->isPeriodic){
+// Vector3r shiftVel = scene->cell->velGrad*scene->cell->Hsize*c->cellDist.cast<Real>();
+// scm->precompute(state1,state2,scene->dt,shiftVel,true);}
+ /*else */scm->precompute(state1,state2,scene,c,true);
return true;
}
return false;
=== modified file 'pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp 2010-07-15 14:26:48 +0000
@@ -36,10 +36,23 @@
Real norm=normal.norm(); normal/=norm; // normal is unit vector now
Real penetrationDepth=s1->radius+s2->radius-norm;
scm->contactPoint=se31.position+(s1->radius-0.5*penetrationDepth)*normal;//0.5*(pt1+pt2);
+ if(isNew) scm->prevNormal=normal;
+ else scm->prevNormal=scm->normal;
+ scm->normal=normal;
scm->penetrationDepth=penetrationDepth;
scm->radius1=s1->radius;
scm->radius2=s2->radius;
- scm->precompute(state1,state2,scene,c,normal,isNew,true);
+#ifdef IGCACHE
+// if (scene->isPeriodic) {
+// Vector3r shiftVel = scene->cell->velGrad*scene->cell->Hsize*c->cellDist.cast<Real>();
+// scm->precompute(state1,state2,scene->dt,shiftVel,true);}
+ /*else */scm->precompute(state1,state2,scene,c,true);
+#endif
+
+ // keep this for reference on how to compute bending and torsion from relative orientation; parts in ScGeom header
+ #if 0
+ scm->initRelOri12=se31.orientation.Conjugate()*se32.orientation;
+ #endif
return true;
}
return false;
=== modified file 'pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp'
--- pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp 2010-07-15 14:26:48 +0000
@@ -32,9 +32,9 @@
Real Vb = mat2->poisson;
Real Ra,Rb; Vector3r normal;
- assert(dynamic_cast<ScGeom*>(interaction->interactionGeometry.get()));//only in debug mode
- ScGeom* sphCont=YADE_CAST<ScGeom*>(interaction->interactionGeometry.get());
- {Ra=sphCont->radius1>0?sphCont->radius1:sphCont->radius2; Rb=sphCont->radius2>0?sphCont->radius2:sphCont->radius1; normal=sphCont->normal;}
+ assert(dynamic_cast<GenericSpheresContact*>(interaction->interactionGeometry.get()));//only in debug mode
+ GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->interactionGeometry.get());
+ {Ra=sphCont->refR1>0?sphCont->refR1:sphCont->refR2; Rb=sphCont->refR2>0?sphCont->refR2:sphCont->refR1; normal=sphCont->normal;}
//harmonic average of the two stiffnesses when (Ri.Ei/2) is the stiffness of a contact point on sphere "i"
Real Kn = 2*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-07-15 14:26:48 +0000
@@ -74,10 +74,12 @@
///////////////////////// CREEP START ///////////
if (shear_creep) shearForce -= currentContactPhysics->ks*(shearForce*dt/creep_viscosity);
///////////////////////// CREEP END ////////////
-
+#ifdef IGCACHE
Vector3r& shearForce = currentContactGeometry->rotate(currentContactPhysics->shearForce);
const Vector3r& dus = currentContactGeometry->shearIncrement();
-
+#else
+ Vector3r dus = currentContactGeometry->rotateAndGetShear(shearForce,currentContactPhysics->prevNormal,de1,de2,dt);
+#endif
//Linear elasticity giving "trial" shear force
shearForce -= currentContactPhysics->ks*dus;
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-07-15 14:26:48 +0000
@@ -11,7 +11,7 @@
CREATE_LOGGER(Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM);
void Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* scene){
-// const Real& dt = scene->dt;
+ const Real& dt = scene->dt;
ScGeom* geom = static_cast<ScGeom*>(ig.get());
CFpmPhys* phys = static_cast<CFpmPhys*>(ip.get());
@@ -67,15 +67,15 @@
State* st1 = Body::byId(id1,scene)->state.get();
State* st2 = Body::byId(id2,scene)->state.get();
-// #ifdef IGCACHE
+#ifdef IGCACHE
geom->rotate(phys->shearForce);
const Vector3r& dus = geom->shearIncrement();
-// #else
-// // define shift to handle periodicity
-// Vector3r shiftVel = scene->isPeriodic ? (Vector3r)((scene->cell->velGrad*scene->cell->Hsize)*Vector3r((Real) contact->cellDist[0],(Real) contact->cellDist[1],(Real) contact->cellDist[2])) : Vector3r::Zero();
-// Vector3r dus = geom->rotateAndGetShear(shearForce, phys->prevNormal, st1, st2, dt, shiftVel, preventGranularRatcheting);
-// #endif
-
+#else
+ // define shift to handle periodicity
+ Vector3r shiftVel = scene->isPeriodic ? (Vector3r)((scene->cell->velGrad*scene->cell->Hsize)*Vector3r((Real) contact->cellDist[0],(Real) contact->cellDist[1],(Real) contact->cellDist[2])) : Vector3r::Zero();
+ Vector3r dus = geom->rotateAndGetShear(shearForce, phys->prevNormal, st1, st2, dt, shiftVel, preventGranularRatcheting);
+#endif
+ /// before changes to adapt periodic scene it was like that: Vector3r dus = geom->rotateAndGetShear(shearForce, phys->prevNormal, st1, st2, dt, preventGranularRatcheting);
//Linear elasticity giving "trial" shear force
shearForce -= phys->ks*dus;
// needed for the next timestep
=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-07-15 14:26:48 +0000
@@ -64,10 +64,15 @@
Real& un=currentContactGeometry->penetrationDepth;
TRVAR3(currentContactGeometry->penetrationDepth,de1->se3.position,de2->se3.position);
currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal;
-
+#ifdef IGCACHE
Vector3r& shearForce = currentContactGeometry->rotate(currentContactPhysics->shearForce);
const Vector3r& shearDisp = currentContactGeometry->shearIncrement();
-
+#else
+ Vector3r& shearForce = currentContactPhysics->shearForce;
+ Vector3r shiftVel = scene->isPeriodic ? (Vector3r)((scene->cell->velGrad*scene->cell->Hsize)*Vector3r((Real) contact->cellDist[0],(Real) contact->cellDist[1],(Real) contact->cellDist[2])) : Vector3r::Zero();
+ Vector3r shearDisp = currentContactGeometry->rotateAndGetShear(shearForce,currentContactPhysics->prevNormal,de1,de2,dt,shiftVel,true);
+// cerr << "shearForce "<<shearForce<<" shearIncrement "<<currentContactGeometry->shearIncrement<<endl;
+#endif
if (!traceEnergy){//Update force but don't compute energy terms (see below))
if(useShear){
currentContactGeometry->updateShear(de1,de2,dt);
@@ -81,7 +86,7 @@
shearForce *= ratio;
if(useShear) currentContactGeometry->shear*=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
+ //almost the same with 2 additional Vector3r instanciated for energy tracing, duplicated block to make sure there is no cost for the instanciation of the vectors when traceEnergy==false
if(useShear) throw ("energy tracing not defined with useShear==true");
shearForce -= currentContactPhysics->ks*shearDisp;
Real maxFs = currentContactPhysics->normalForce.squaredNorm()*
=== modified file 'pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.cpp'
--- pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.cpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.cpp 2010-07-15 14:26:48 +0000
@@ -102,10 +102,10 @@
FOREACH(const shared_ptr<Interaction>& contact, *rb->interactions){
if(!contact->isReal()) continue;
- ScGeom* geom=YADE_CAST<ScGeom*>(contact->interactionGeometry.get()); assert(geom);
+ GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(contact->interactionGeometry.get()); assert(geom);
NormShearPhys* phys=YADE_CAST<NormShearPhys*>(contact->interactionPhysics.get()); assert(phys);
// all we need for getting stiffness
- Vector3r& normal=geom->normal; Real& kn=phys->kn; Real& ks=phys->ks; Real& radius1=geom->radius1; Real& radius2=geom->radius1;
+ Vector3r& normal=geom->normal; Real& kn=phys->kn; Real& ks=phys->ks; Real& radius1=geom->refR1; Real& radius2=geom->refR2;
Real fn = (static_cast<NormShearPhys *> (contact->interactionPhysics.get()))->normalForce.squaredNorm();
if (fn==0) continue;//Is it a problem with some laws? I still don't see why.
=== modified file 'pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp'
--- pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp 2010-07-15 14:26:48 +0000
@@ -112,8 +112,7 @@
// Define shift to handle periodicity
Vector3r shiftVel = scene->isPeriodic ? (Vector3r)((scene->cell->velGrad*scene->cell->Hsize)*Vector3r((Real) contact->cellDist[0],(Real) contact->cellDist[1],(Real) contact->cellDist[2])) : Vector3r::Zero();
// 1. Rotate shear force
-// shearElastic = scg->rotate(shearElastic, phys->prevNormal, de1, de2, dt);
- shearElastic = scg->rotate(shearElastic);
+ shearElastic = scg->rotate(shearElastic, phys->prevNormal, de1, de2, dt);
// 2. Get incident velocity, get shear and normal components
// FIXME: Concerning the possibility to avoid granular ratcheting, it is not clear how this works in the case of HM. See the thread http://www.mail-archive.com/yade-users@xxxxxxxxxxxxxxxxxxx/msg01947.html
Vector3r incidentV = scg->getIncidentVel(de1, de2, dt, shiftVel, preventGranularRatcheting);
=== modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2010-07-15 14:26:48 +0000
@@ -118,12 +118,12 @@
FOREACH(const shared_ptr<Interaction>&I, *scene->interactions){
if ( !I->isReal() ) continue;
NormShearPhys* nsi=YADE_CAST<NormShearPhys*> ( I->interactionPhysics.get() );
- ScGeom* gsc=YADE_CAST<ScGeom*> ( I->interactionGeometry.get() );
+ GenericSpheresContact* gsc=YADE_CAST<GenericSpheresContact*> ( I->interactionGeometry.get() );
//Contact force
Vector3r f= ( reversedForces?-1.:1. ) * ( nsi->normalForce+nsi->shearForce );
//branch vector, FIXME : the first definition generalizes to non-spherical bodies but needs wrapped coords.
// Vector3r branch=(Body::byId(I->getId1())->state->pos-Body::byId(I->getId2())->state->pos);
- Vector3r branch= gsc->normal* ( gsc->radius1+gsc->radius2 );
+ Vector3r branch= gsc->normal* ( gsc->refR1+gsc->refR2 );
#if 0
// remove this block later
// tensorial product f*branch (hand-write the tensor product to prevent matrix instanciation inside the loop by makeTensorProduct)
=== modified file 'pkg/dem/meta/Shop.cpp'
--- pkg/dem/meta/Shop.cpp 2010-07-14 19:02:58 +0000
+++ pkg/dem/meta/Shop.cpp 2010-07-15 14:26:48 +0000
@@ -40,9 +40,8 @@
#include<yade/pkg-dem/GlobalStiffnessTimeStepper.hpp>
#include<yade/pkg-dem/ElasticContactLaw.hpp>
-#include<yade/pkg-dem/DemXDofGeom.hpp>
+
#include<yade/pkg-dem/ScGeom.hpp>
-
#include<yade/pkg-dem/FrictPhys.hpp>
=== modified file 'py/_eudoxos.cpp'
--- py/_eudoxos.cpp 2010-07-14 19:02:58 +0000
+++ py/_eudoxos.cpp 2010-07-15 14:26:48 +0000
@@ -4,7 +4,6 @@
#include<boost/python.hpp>
#include<yade/extra/boost_python_len.hpp>
#include<yade/pkg-dem/Shop.hpp>
-#include <yade/pkg-dem/DemXDofGeom.hpp>
#ifdef YADE_VTK
#include<vtkPointLocator.h>
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2010-07-14 19:18:03 +0000
+++ py/_utils.cpp 2010-07-15 14:26:48 +0000
@@ -127,7 +127,7 @@
const shared_ptr<Body>& b1=Body::byId(i->getId1(),rb), b2=Body::byId(i->getId2(),rb);
if(!b1->maskOk(mask) || !b2->maskOk(mask)) continue;
if(useBB && !isInBB(b1->state->pos,bbMin,bbMax) && !isInBB(b2->state->pos,bbMin,bbMax)) continue;
- ScGeom* geom=dynamic_cast<ScGeom*>(i->interactionGeometry.get());
+ GenericSpheresContact* geom=dynamic_cast<GenericSpheresContact*>(i->interactionGeometry.get());
if(!geom) continue;
Vector3r n(geom->normal); n[axis]=0.; Real nLen=n.norm();
if(nLen<minProjLen) continue; // this interaction is (almost) exactly parallel to our axis; skip that one
=== modified file 'py/utils.py'
--- py/utils.py 2010-07-14 19:02:58 +0000
+++ py/utils.py 2010-07-15 14:26:48 +0000
@@ -231,14 +231,14 @@
: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).
+ 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
b=Body()
b.shape=ChainedCylinder(radius=radius,length=segment.norm(),color=color if color else randomColor(),wire=wire,highlight=highlight)
- b.shape.segment=segment
+ b.shape.segment=b.shape.length*Vector3(0.,0.,1.)
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()
=== modified file 'scripts/test/chained-cylinder-spring.py'
--- scripts/test/chained-cylinder-spring.py 2010-07-14 19:02:58 +0000
+++ scripts/test/chained-cylinder-spring.py 2010-07-15 14:26:48 +0000
@@ -38,7 +38,7 @@
]
#Generate a spiral
-Ne=200
+Ne=400
for i in range(0, Ne):
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;
