yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05942
[Branch ~yade-dev/yade/trunk] Rev 2504: - Add new IG class ScGeom6D precomputing relative rotations, and associated Ig functors.
------------------------------------------------------------
revno: 2504
committer: bchareyre <bchareyre@dt-rv020>
branch nick: yade
timestamp: Wed 2010-10-20 13:31:49 +0200
message:
- Add new IG class ScGeom6D precomputing relative rotations, and associated Ig functors.
- Update related classes (CFLaw)
- add kdev4 configuration for starting debug from kdevelop
modified:
.kdev4/yade.kdev4
doc/README
doc/yade-articles.bib
pkg/common/Cylinder.cpp
pkg/dem/CohFrictPhys.hpp
pkg/dem/CohesiveFrictionalContactLaw.cpp
pkg/dem/CohesiveFrictionalContactLaw.hpp
pkg/dem/CohesiveTriaxialTest.cpp
pkg/dem/Ig2_Sphere_Sphere_ScGeom.cpp
pkg/dem/Ig2_Sphere_Sphere_ScGeom.hpp
pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.cpp
pkg/dem/ScGeom.cpp
pkg/dem/ScGeom.hpp
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 '.kdev4/yade.kdev4'
--- .kdev4/yade.kdev4 2010-10-18 19:27:07 +0000
+++ .kdev4/yade.kdev4 2010-10-20 11:31:49 +0000
@@ -1,6 +1,31 @@
[Buildset]
BuildItems=@Variant(\x00\x00\x00\t\x00\x00\x00\x00\x01\x00\x00\x00\x0b\x00\x00\x00\x00\x01\x00\x00\x00\x08\x00y\x00a\x00d\x00e)
+[Launch]
+Launch Configurations=Launch Configuration 0
+
+[Launch][Launch Configuration 0]
+Configured Launch Modes=execute
+Configured Launchers=nativeAppLauncher
+Name=yade-debug
+Type=Native Application
+
+[Launch][Launch Configuration 0][Data]
+Arguments=\s/home/3S-LAB/bchareyre/yade/yade-test-kdev/bin/yade-true-true --debug -j 1 /home/3S-LAB/bchareyre2/yade/yade-test-kdev/yade/scripts/test/chained-cylinder-spring.py
+Debugger Shell=
+Dependencies=@Variant(\x00\x00\x00\t\x00\x00\x00\x00\x00)
+Dependency Action=Nothing
+Display Demangle Names=true
+Display Static Members=false
+EnvironmentGroup=default
+Executable=python
+GDB Path=gdb
+Remote GDB Config Script=
+Remote GDB Run Script=
+Remote GDB Shell Script=
+Try Setting Breakpoints On Loading Libraries=true
+Working Directory=
+isExecutable=true
[MakeBuilder]
Make Binary=scons
Number Of Jobs=3
=== modified file 'doc/README'
--- doc/README 2010-01-09 13:13:52 +0000
+++ doc/README 2010-10-20 11:31:49 +0000
@@ -1,9 +1,6 @@
All documentation is at http://www.yade-dem.org, in particular see
https://yade-dem.org/index.php/Reference_documentation
-To generate documentation (doxygen for c++ and epydoc for python classes), run
-
- doxygen Doxyfile
- yade-* yade-epydoc.py
+To generate documentation see sphinx/README
=== modified file 'doc/yade-articles.bib'
--- doc/yade-articles.bib 2010-07-16 18:12:30 +0000
+++ doc/yade-articles.bib 2010-10-20 11:31:49 +0000
@@ -264,3 +264,10 @@
number = {B},
year = {2004}
}
+
+@Article{ Hassan2010,
+ author = "A. Hassan and B. Chareyre and F. Darve and J. Meyssonier and F. Flin",
+ title = "Microtomography-based Discrete Element Modelling of Creep in Snow",
+ journal = "Granular Matter",
+ year = "2010 (submitted)"
+}
=== modified file 'pkg/common/Cylinder.cpp'
--- pkg/common/Cylinder.cpp 2010-10-16 18:31:17 +0000
+++ pkg/common/Cylinder.cpp 2010-10-20 11:31:49 +0000
@@ -103,10 +103,10 @@
ChainedCylinder *bs1=static_cast<ChainedCylinder*>(revert? cm2.get():cm1.get());
- shared_ptr<ScGeom> scm;
+ shared_ptr<ScGeom6D> scm;
bool isNew = !c->geom;
- if(!isNew) scm=YADE_PTR_CAST<ScGeom>(c->geom);
- else { scm=shared_ptr<ScGeom>(new ScGeom()); c->geom=scm; }
+ if(!isNew) scm=YADE_PTR_CAST<ScGeom6D>(c->geom);
+ else { scm=shared_ptr<ScGeom6D>(new ScGeom6D()); c->geom=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;}
@@ -122,7 +122,7 @@
bs1->chainedOrientation.setFromTwoVectors(Vector3r::UnitZ(),bchain1.ori.conjugate()*segt);
#endif
scm->precompute(state1,state2,scene,c,segt/length,isNew,shift2,true);
-
+ scm->precomputeRotations(state1,state2,isNew,false);
//Set values that will be considered in Ip2 functor, geometry (precomputed) is really defined with values above
scm->radius1 = scm->radius2 = bs1->initLength*0.5;
return true;
=== modified file 'pkg/dem/CohFrictPhys.hpp'
--- pkg/dem/CohFrictPhys.hpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/CohFrictPhys.hpp 2010-10-20 11:31:49 +0000
@@ -27,15 +27,6 @@
// internal attributes
((Vector3r,moment_twist,Vector3r(0,0,0),,""))
((Vector3r,moment_bending,Vector3r(0,0,0),,""))
- ((Vector3r,initialPosition1,Vector3r(0,0,0),,""))
- ((Vector3r,initialPosition2,Vector3r(0,0,0),,""))
- ((Quaternionr,orientationToContact1,Quaternionr(1.0,0.0,0.0,0.0),,""))
- ((Quaternionr,orientationToContact2,Quaternionr(1.0,0.0,0.0,0.0),,""))
- ((Quaternionr,initialOrientation1,Quaternionr(1.0,0.0,0.0,0.0),,""))
- ((Quaternionr,initialOrientation2,Quaternionr(1.0,0.0,0.0,0.0),,""))
- ((Quaternionr,twistCreep,Quaternionr(1.0,0.0,0.0,0.0),,""))
- ((Quaternionr,currentContactOrientation,Quaternionr(1.0,0.0,0.0,0.0),,""))
- ((Quaternionr,initialContactOrientation,Quaternionr(1.0,0.0,0.0,0.0),,""))
,
createIndex();
);
=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/CohesiveFrictionalContactLaw.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.cpp 2010-10-20 11:31:49 +0000
@@ -37,20 +37,12 @@
}
}
-
void Law2_ScGeom_CohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact)
{
const Real& dt = scene->dt;
-// if (detectBrokenBodies //Experimental, has no effect
-// && (*bodies)[contact->getId1()]->shape->getClassName() != "box"
-// && (*bodies)[contact->getId2()]->shape->getClassName() != "box") {
-// YADE_CAST<CohFrictMat*> ((*bodies)[contact->getId1()]->material.get())->isBroken = false;
-// YADE_CAST<CohFrictMat*> ((*bodies)[contact->getId2()]->material.get())->isBroken = false;}
const int &id1 = contact->getId1();
const int &id2 = contact->getId2();
- Body* b1 = Body::byId(id1,scene).get();
- Body* b2 = Body::byId(id2,scene).get();
- ScGeom* currentContactGeometry = YADE_CAST<ScGeom*> (ig.get());
+ ScGeom6D* currentContactGeometry = YADE_CAST<ScGeom6D*> (ig.get());
CohFrictPhys* currentContactPhysics = YADE_CAST<CohFrictPhys*> (ip.get());
Vector3r& shearForce = currentContactPhysics->shearForce;
@@ -61,6 +53,7 @@
if (un < 0 && (currentContactPhysics->normalForce.squaredNorm() > pow(currentContactPhysics->normalAdhesion,2)
|| currentContactPhysics->normalAdhesion==0)) {
// BREAK due to tension
+
scene->interactions->requestErase(contact->getId1(),contact->getId2());
} else {
State* de1 = Body::byId(id1,scene)->state.get();
@@ -93,45 +86,17 @@
/// Moment law ///
if (momentRotationLaw && (!currentContactPhysics->cohesionBroken || always_use_moment_law)) {
- // Not necessary. OK.
- //{// updates only orientation of contact (local coordinate system)
- // Vector3r axis = currentContactPhysics->prevNormal.UnitCross(currentContactGeometry->normal);
- // Real angle = unitVectorsAngle(currentContactPhysics->prevNormal,currentContactGeometry->normal);
- // Quaternionr align(axis,angle);
- // currentContactPhysics->currentContactOrientation = align * currentContactPhysics->currentContactOrientation;
- //}
- Quaternionr delta((b1->state->ori * (currentContactPhysics->initialOrientation1.conjugate())) *
- (currentContactPhysics->initialOrientation2 * (b2->state->ori.conjugate())));
- if (twist_creep) {
- delta = delta * currentContactPhysics->twistCreep;
- }
-
- AngleAxisr aa(delta); // axis of rotation - this is the Moment direction UNIT vector; // angle represents the power of resistant ELASTIC moment
- //Eigen::AngleAxisr(q) returns nan's when q close to identity, next tline fixes the pb.
- if (isnan(aa.angle())) aa.angle()=0;
- if (aa.angle() > Mathr::PI) aa.angle() -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi
- Real angle_twist(aa.angle() * aa.axis().dot(currentContactGeometry->normal));
- Vector3r axis_twist(angle_twist * currentContactGeometry->normal);
-
if (twist_creep) {
Real viscosity_twist = creep_viscosity * std::pow((2 * std::min(currentContactGeometry->radius1,currentContactGeometry->radius2)),2) / 16.0;
- Real angle_twist_creeped = angle_twist * (1 - dt/viscosity_twist);
- Quaternionr q_twist(AngleAxisr(angle_twist,currentContactGeometry->normal));
- //Quaternionr q_twist_creeped(currentContactGeometry->normal , angle_twist*0.996);
+ Real angle_twist_creeped = currentContactGeometry->getTwist() * (1 - dt/viscosity_twist);
+ Quaternionr q_twist(AngleAxisr(currentContactGeometry->getTwist(),currentContactGeometry->normal));
Quaternionr q_twist_creeped(AngleAxisr(angle_twist_creeped,currentContactGeometry->normal));
Quaternionr q_twist_delta(q_twist_creeped * q_twist.conjugate());
- currentContactPhysics->twistCreep = currentContactPhysics->twistCreep * q_twist_delta;
- // modify the initialRelativeOrientation to substract some twisting
- // currentContactPhysics->initialRelativeOrientation = currentContactPhysics->initialRelativeOrientation * q_twist_delta;
- //currentContactPhysics->initialOrientation1 = currentContactPhysics->initialOrientation1 * q_twist_delta;
- //currentContactPhysics->initialOrientation2 = currentContactPhysics->initialOrientation2 * q_twist_delta.Conjugate();
+ currentContactGeometry->twistCreep = currentContactGeometry->twistCreep * q_twist_delta;
}
- Vector3r moment_twist(axis_twist * currentContactPhysics->kr);
- Vector3r axis_bending(aa.angle()*aa.axis() - axis_twist);
- Vector3r moment_bending(axis_bending * currentContactPhysics->kr);
- Vector3r moment = moment_twist + moment_bending;
- currentContactPhysics->moment_twist = moment_twist;
- currentContactPhysics->moment_bending = moment_bending;
+ currentContactPhysics->moment_twist = (currentContactGeometry->getTwist()*currentContactPhysics->kr)*currentContactGeometry->normal;
+ currentContactPhysics->moment_bending = currentContactGeometry->getBending() * currentContactPhysics->kr;
+ Vector3r moment = currentContactPhysics->moment_twist + currentContactPhysics->moment_bending;
scene->forces.addTorque(id1,-moment);
scene->forces.addTorque(id2, moment);
}
=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/CohesiveFrictionalContactLaw.hpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.hpp 2010-10-20 11:31:49 +0000
@@ -16,7 +16,7 @@
class Law2_ScGeom_CohFrictPhys_CohesionMoment: public LawFunctor{
public:
virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
- YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_CohFrictPhys_CohesionMoment,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. This law adds adhesion and moments to :yref:`Law2_ScGeom_FrictPhys_CundallStrack`.\n\nThe normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, a_n)$, with $a_n$ the normal adhesion. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force, by default $F_s^{max}=F_n*tan(\\phi)+a_s$, with $\\phi$ the friction angle and $a_n$ the shear adhesion. If :yref:`CohFrictPhys::cohesionDisableFriction` is True, friction is ignored as long as adhesion is active, and the maximum shear force is only $F_s^{max}=a_s$.\n\nIf the maximum tensile or maximum shear force is reached and :yref:`CohFrictPhys::fragile` =True (default), the cohesive link is broken, and $a_n, a_s$ are set back to zero. If a tensile force is present, the contact is lost, else the shear strength is $F_s^{max}=F_n*tan(\\phi)$. If :yref:`CohFrictPhys::fragile` =False (in course of implementation), the behaviour is perfectly plastic, and the shear strength is kept constant.\n\nIf :yref:`Law2_ScGeom_CohFrictPhys_CohesionMoment::momentRotationLaw` =True, bending and twisting moments are computed using a linear law with moduli respectively $k_t$ and $k_r$ (the two values are the same currently), so that the moments are : $M_b=k_b*\\Theta_b$ and $M_t=k_t*\\Theta_t$, with $\\Theta_{b,t}$ the relative rotations between interacting bodies. There is no maximum value of moments in the current implementation, though they could be added in the future.\n\nCreep at contact is implemented in this law, as defined in [Hassan2010]_. If activated, there is a viscous behaviour of the shear and twisting components, and the evolution of the elastic parts of shear displacement and relative twist is given by $du_{s,e}/dt=F_s/\\nu_s$ and $d\\Theta_{t,e}/dt=M_t/\\nu_t$.",
+ YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_CohFrictPhys_CohesionMoment,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. This law adds adhesion and moments to :yref:`Law2_ScGeom_FrictPhys_CundallStrack`.\n\nThe normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, a_n)$, with $a_n$ the normal adhesion. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force, by default $F_s^{max}=F_n*tan(\\phi)+a_s$, with $\\phi$ the friction angle and $a_n$ the shear adhesion. If :yref:`CohFrictPhys::cohesionDisableFriction` is True, friction is ignored as long as adhesion is active, and the maximum shear force is only $F_s^{max}=a_s$.\n\nIf the maximum tensile or maximum shear force is reached and :yref:`CohFrictPhys::fragile` =True (default), the cohesive link is broken, and $a_n, a_s$ are set back to zero. If a tensile force is present, the contact is lost, else the shear strength is $F_s^{max}=F_n*tan(\\phi)$. If :yref:`CohFrictPhys::fragile` =False (in course of implementation), the behaviour is perfectly plastic, and the shear strength is kept constant.\n\nIf :yref:`Law2_ScGeom_CohFrictPhys_CohesionMoment::momentRotationLaw` =True, bending and twisting moments are computed using a linear law with moduli respectively $k_t$ and $k_r$ (the two values are the same currently), so that the moments are : $M_b=k_b*\\Theta_b$ and $M_t=k_t*\\Theta_t$, with $\\Theta_{b,t}$ the relative rotations between interacting bodies. There is no maximum value of moments in the current implementation, though they could be added in the future.\n\nCreep at contact is implemented in this law, as defined in [Hassan2010]_. If activated, there is a viscous behaviour of the shear and twisting components, and the evolution of the elastic parts of shear displacement and relative twist is given by $du_{s,e}/dt=-F_s/\\nu_s$ and $d\\Theta_{t,e}/dt=-M_t/\\nu_t$.",
((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
((bool,momentRotationLaw,false,,"use bending/twisting moment at contacts. See :yref:`CohFrictPhys::cohesionDisablesFriction` for details."))
((bool,always_use_moment_law,false,,"If true, use bending/twisting moments at all contacts. If false, compute moments only for cohesive contacts."))
@@ -24,7 +24,7 @@
((bool,twist_creep,false,,"activate creep on the twisting moment, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
((Real,creep_viscosity,1,,"creep viscosity [Pa.s/m]. probably should be moved to Ip2_2xCohFrictMat_CohFrictPhys..."))
);
- FUNCTOR2D(ScGeom,CohFrictPhys);
+ FUNCTOR2D(ScGeom6D,CohFrictPhys);
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(Law2_ScGeom_CohFrictPhys_CohesionMoment);
=== modified file 'pkg/dem/CohesiveTriaxialTest.cpp'
--- pkg/dem/CohesiveTriaxialTest.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/CohesiveTriaxialTest.cpp 2010-10-20 11:31:49 +0000
@@ -287,7 +287,7 @@
{
shared_ptr<IGeomDispatcher> interactionGeometryDispatcher(new IGeomDispatcher);
- shared_ptr<IGeomFunctor> s1(new Ig2_Sphere_Sphere_ScGeom);
+ shared_ptr<IGeomFunctor> s1(new Ig2_Sphere_Sphere_ScGeom6D);
interactionGeometryDispatcher->add(s1);
shared_ptr<IGeomFunctor> s2(new Ig2_Box_Sphere_ScGeom);
interactionGeometryDispatcher->add(s2);
=== modified file 'pkg/dem/Ig2_Sphere_Sphere_ScGeom.cpp'
--- pkg/dem/Ig2_Sphere_Sphere_ScGeom.cpp 2010-10-16 18:31:17 +0000
+++ pkg/dem/Ig2_Sphere_Sphere_ScGeom.cpp 2010-10-20 11:31:49 +0000
@@ -1,4 +1,3 @@
-// © 2004 Olivier Galizzi <olivier.galizzi@xxxxxxx>
// © 2004 Janek Kozicki <cosurgi@xxxxxxxxxx>
// © 2007 Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
// © 2008 Václav Šmilauer <eudoxos@xxxxxxxx>
@@ -23,7 +22,6 @@
{
#endif
const Se3r& se31=state1.se3; const Se3r& se32=state2.se3;
-
Sphere *s1=static_cast<Sphere*>(cm1.get()), *s2=static_cast<Sphere*>(cm2.get());
Vector3r normal=(se32.position+shift2)-se31.position;
Real penetrationDepthSq=pow(interactionDetectionFactor*(s1->radius+s2->radius),2) - normal.squaredNorm();
@@ -45,7 +43,6 @@
return false;
}
-
bool Ig2_Sphere_Sphere_ScGeom::goReverse( const shared_ptr<Shape>& cm1,
const shared_ptr<Shape>& cm2,
const State& state1,
@@ -58,3 +55,37 @@
}
YADE_PLUGIN((Ig2_Sphere_Sphere_ScGeom));
+
+#ifdef YADE_DEVIRT_FUNCTORS
+bool Ig2_Sphere_Sphere_ScGeom6D::go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){ throw runtime_error("Do not call Ig2_Sphere_Sphere_ScGeom6D::go, use getStaticFunctorPtr and call that function instead."); }
+bool Ig2_Sphere_Sphere_ScGeom6D::goStatic(IGeomFunctor* _self, const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){
+ const Ig2_Sphere_Sphere_ScGeom6D* self=static_cast<Ig2_Sphere_Sphere_ScGeom*>(_self);
+ const Real& interactionDetectionFactor=self->interactionDetectionFactor;
+#else
+bool Ig2_Sphere_Sphere_ScGeom6D::go( const shared_ptr<Shape>& cm1,
+ const shared_ptr<Shape>& cm2,
+ const State& state1, const State& state2, const Vector3r& shift2, const bool& force,
+ const shared_ptr<Interaction>& c)
+{
+#endif
+ bool isNew = !c->geom;
+ if (Ig2_Sphere_Sphere_ScGeom::go(cm1,cm2,state1,state2,shift2,force,c)){//the 3 DOFS from ScGeom are updated here
+ shared_ptr<ScGeom6D> scm=YADE_PTR_CAST<ScGeom6D>(c->geom);
+ if (updateRotations) scm->precomputeRotations(state1,state2,isNew,creep);
+ return true;
+ }
+ else return false;
+}
+
+bool Ig2_Sphere_Sphere_ScGeom6D::goReverse( const shared_ptr<Shape>& cm1,
+ const shared_ptr<Shape>& cm2,
+ const State& state1,
+ const State& state2,
+ const Vector3r& shift2,
+ const bool& force,
+ const shared_ptr<Interaction>& c)
+{
+ return go(cm1,cm2,state2,state1,-shift2,force,c);
+}
+
+YADE_PLUGIN((Ig2_Sphere_Sphere_ScGeom6D));
=== modified file 'pkg/dem/Ig2_Sphere_Sphere_ScGeom.hpp'
--- pkg/dem/Ig2_Sphere_Sphere_ScGeom.hpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/Ig2_Sphere_Sphere_ScGeom.hpp 2010-10-20 11:31:49 +0000
@@ -1,10 +1,6 @@
-/*************************************************************************
-* Copyright (C) 2004 by Olivier Galizzi *
-* olivier.galizzi@xxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
+// © 2004 Janek Kozicki <cosurgi@xxxxxxxxxx>
+// © 2007 Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
+// © 2008 Václav Šmilauer <eudoxos@xxxxxxxx>
#pragma once
@@ -38,3 +34,22 @@
};
REGISTER_SERIALIZABLE(Ig2_Sphere_Sphere_ScGeom);
+class Ig2_Sphere_Sphere_ScGeom6D: public Ig2_Sphere_Sphere_ScGeom{
+ public:
+ virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ virtual bool goReverse( const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ #ifdef YADE_DEVIRT_FUNCTORS
+ void* getStaticFuncPtr(){ return (void*)&Ig2_Sphere_Sphere_ScGeom6D::goStatic; }
+ static bool goStatic(IGeomFunctor* self, const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& se32, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ #endif
+ YADE_CLASS_BASE_DOC_ATTRS(Ig2_Sphere_Sphere_ScGeom6D,Ig2_Sphere_Sphere_ScGeom,"Create/update a :yref:`ScGeom6D` instance representing intersection of two :yref:`Spheres<Sphere>`.",
+ ((bool,updateRotations,true,,"Precompute relative rotations. Turning this false can speed up simulations when rotations are not needed in constitutive laws (e.g. when spheres are compressed without cohesion and moment in early stage of a triaxial test), but is not foolproof. Change this value only if you know what you are doing."))
+ ((bool,creep,false,,"Substract rotational creep from relative rotation. The rotational creep :yref:`ScGeom6D::twistCreep` is a quaternion and has to be updated inside a constitutive law, see for instance :yref:`Law2_ScGeom_CohFrictPhys_CohesionMoment`."
+ ))
+ );
+ FUNCTOR2D(Sphere,Sphere);
+ // needed for the dispatcher, even if it is symmetric
+ DEFINE_FUNCTOR_ORDER_2D(Sphere,Sphere);
+};
+REGISTER_SERIALIZABLE(Ig2_Sphere_Sphere_ScGeom6D);
+
=== modified file 'pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.cpp'
--- pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.cpp 2010-10-20 11:31:49 +0000
@@ -19,7 +19,7 @@
{
CohFrictMat* sdec1 = static_cast<CohFrictMat*>(b1.get());
CohFrictMat* sdec2 = static_cast<CohFrictMat*>(b2.get());
- ScGeom* geom = YADE_CAST<ScGeom*>(interaction->geom.get());
+ ScGeom6D* geom = YADE_CAST<ScGeom6D*>(interaction->geom.get());
//Create cohesive interractions only once
if (setCohesionNow && cohesionDefinitionIteration==-1) {
@@ -29,10 +29,8 @@
cohesionDefinitionIteration = -1;
setCohesionNow = 0;}
- if (geom)
- {
- if (!interaction->phys)
- {
+ if (geom) {
+ if (!interaction->phys) {
interaction->phys = shared_ptr<CohFrictPhys>(new CohFrictPhys());
CohFrictPhys* contactPhysics = YADE_CAST<CohFrictPhys*>(interaction->phys.get());
Real Ea = sdec1->young;
@@ -61,47 +59,22 @@
contactPhysics->cohesionBroken = false;
contactPhysics->normalAdhesion = normalCohesion*pow(std::min(Db, Da),2);
contactPhysics->shearAdhesion = shearCohesion*pow(std::min(Db, Da),2);;
-
- contactPhysics->initialOrientation1 = Body::byId(interaction->getId1())->state->ori;
- contactPhysics->initialOrientation2 = Body::byId(interaction->getId2())->state->ori;
- contactPhysics->initialPosition1 = Body::byId(interaction->getId1())->state->pos;
- contactPhysics->initialPosition2 = Body::byId(interaction->getId2())->state->pos;
- contactPhysics->kr = Kr;
- contactPhysics->initialContactOrientation.setFromTwoVectors(Vector3r(1.0,0.0,0.0),geom->normal);
- contactPhysics->currentContactOrientation = contactPhysics->initialContactOrientation;
- contactPhysics->orientationToContact1 = contactPhysics->initialOrientation1.conjugate() * contactPhysics->initialContactOrientation;
- contactPhysics->orientationToContact2 = contactPhysics->initialOrientation2.conjugate() * contactPhysics->initialContactOrientation;
+ geom->initRotations(*(Body::byId(interaction->getId1(),scene)->state),*(Body::byId(interaction->getId2(),scene)->state));
}
- contactPhysics->prevNormal = geom->normal;
contactPhysics->kn = Kn;
contactPhysics->ks = Ks;
- contactPhysics->initialOrientation1 = Body::byId(interaction->getId1())->state->ori;
- contactPhysics->initialOrientation2 = Body::byId(interaction->getId2())->state->ori;
- contactPhysics->initialPosition1 = Body::byId(interaction->getId1())->state->pos;
- contactPhysics->initialPosition2 = Body::byId(interaction->getId2())->state->pos;
contactPhysics->kr = Kr;
- contactPhysics->initialContactOrientation.setFromTwoVectors(Vector3r(1.0,0.0,0.0),geom->normal);
- contactPhysics->currentContactOrientation = contactPhysics->initialContactOrientation;
- contactPhysics->orientationToContact1 = contactPhysics->initialOrientation1.conjugate() * contactPhysics->initialContactOrientation;
- contactPhysics->orientationToContact2 = contactPhysics->initialOrientation2.conjugate() * contactPhysics->initialContactOrientation;
//contactPhysics->elasticRollingLimit = elasticRollingLimit;
+
}
- else // !isNew, but if setCohesionNow, all contacts are initialized like if they were newly created
- {
+ else {// !isNew, but if setCohesionNow, all contacts are initialized like if they were newly created
CohFrictPhys* contactPhysics = YADE_CAST<CohFrictPhys*>(interaction->phys.get());
if (setCohesionNow && sdec1->isCohesive && sdec2->isCohesive)
{
contactPhysics->cohesionBroken = false;
contactPhysics->normalAdhesion = normalCohesion*pow(std::min(geom->radius2, geom->radius1),2);
contactPhysics->shearAdhesion = shearCohesion*pow(std::min(geom->radius2, geom->radius1),2);
- contactPhysics->initialOrientation1 = Body::byId(interaction->getId1())->state->ori;
- contactPhysics->initialOrientation2 = Body::byId(interaction->getId2())->state->ori;
- contactPhysics->initialPosition1 = Body::byId(interaction->getId1())->state->pos;
- contactPhysics->initialPosition2 = Body::byId(interaction->getId2())->state->pos;
- contactPhysics->initialContactOrientation.setFromTwoVectors(Vector3r(1.0,0.0,0.0),geom->normal);
- contactPhysics->currentContactOrientation = contactPhysics->initialContactOrientation;
- contactPhysics->orientationToContact1 = contactPhysics->initialOrientation1.conjugate() * contactPhysics->initialContactOrientation;
- contactPhysics->orientationToContact2 = contactPhysics->initialOrientation2.conjugate() * contactPhysics->initialContactOrientation;
+ geom->initRotations(*(Body::byId(interaction->getId1(),scene)->state),*(Body::byId(interaction->getId2(),scene)->state));
}
}
}
=== modified file 'pkg/dem/ScGeom.cpp'
--- pkg/dem/ScGeom.cpp 2010-10-16 18:31:17 +0000
+++ pkg/dem/ScGeom.cpp 2010-10-20 11:31:49 +0000
@@ -7,8 +7,9 @@
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
-YADE_PLUGIN((ScGeom));
+YADE_PLUGIN((ScGeom)(ScGeom6D));
ScGeom::~ScGeom(){};
+ScGeom6D::~ScGeom6D(){};
Vector3r& ScGeom::rotate(Vector3r& shearForce) const {
// approximated rotations
@@ -83,3 +84,32 @@
scene->isPeriodic ? Vector3r(scene->cell->velGrad*scene->cell->Hsize*i->cellDist.cast<Real>()) : Vector3r::Zero(), // shiftVel
avoidGranularRatcheting);
}
+
+//!Precompute relative rotations (and precompute ScGeom3D)
+void ScGeom6D::precomputeRotations(const State& rbp1, const State& rbp2, bool isNew, bool creep){
+ if (isNew) {
+ initRotations(rbp1,rbp2);
+ } else {
+ Quaternionr delta((rbp1.ori * (initialOrientation1.conjugate())) * (initialOrientation2 * (rbp2.ori.conjugate())));
+ if (creep) delta = delta * twistCreep;
+ AngleAxisr aa(delta); // axis of rotation - this is the Moment direction UNIT vector; // angle represents the power of resistant ELASTIC moment
+ //Eigen::AngleAxisr(q) returns nan's when q close to identity, next tline fixes the pb.
+ if (isnan(aa.angle())) aa.angle()=0;
+ if (aa.angle() > Mathr::PI) aa.angle() -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi
+ twist = (aa.angle() * aa.axis().dot(normal));
+ bending = Vector3r(aa.angle()*aa.axis() - twist*normal);
+ }
+}
+
+void ScGeom6D::initRotations(const State& state1, const State& state2)
+{
+ initialOrientation1 = state1.ori;
+ initialOrientation2 = state2.ori;
+ initialContactOrientation.setFromTwoVectors(Vector3r(1.0,0.0,0.0),normal);
+ currentContactOrientation = initialContactOrientation;
+ orientationToContact1 = initialOrientation1.conjugate() * initialContactOrientation;
+ orientationToContact2 = initialOrientation2.conjugate() * initialContactOrientation;
+ twist=0;
+ bending=Vector3r::Zero();
+ twistCreep=Quaternionr(1.0,0.0,0.0,0.0);
+}
=== modified file 'pkg/dem/ScGeom.hpp'
--- pkg/dem/ScGeom.hpp 2010-10-16 18:31:17 +0000
+++ pkg/dem/ScGeom.hpp 2010-10-20 11:31:49 +0000
@@ -53,3 +53,31 @@
};
REGISTER_SERIALIZABLE(ScGeom);
+class ScGeom6D: public ScGeom {
+ public:
+ virtual ~ScGeom6D();
+ const Real& getTwist() const {return twist;}
+ const Vector3r& getBending() const {return bending;}
+
+ void precomputeRotations(const State& rbp1, const State& rbp2, bool isNew, bool creep=false);
+ void initRotations(const State& rbp1, const State& rbp2);
+
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom6D,ScGeom,"Class representing :yref:`geometry<IGeom>` of two :yref:`bodies<Body>` in contact. The contact has 6 DOFs (normal, 2×shear, twist, 2xbending) and uses :yref:`ScGeom` incremental algorithm for updating shear.",
+ ((Quaternionr,orientationToContact1,Quaternionr(1.0,0.0,0.0,0.0),,""))
+ ((Quaternionr,orientationToContact2,Quaternionr(1.0,0.0,0.0,0.0),,""))
+ ((Quaternionr,initialOrientation1,Quaternionr(1.0,0.0,0.0,0.0),,""))
+ ((Quaternionr,initialOrientation2,Quaternionr(1.0,0.0,0.0,0.0),,""))
+ ((Quaternionr,twistCreep,Quaternionr(1.0,0.0,0.0,0.0),,""))
+ ((Quaternionr,currentContactOrientation,Quaternionr(1.0,0.0,0.0,0.0),,""))
+ ((Quaternionr,initialContactOrientation,Quaternionr(1.0,0.0,0.0,0.0),,""))
+ ((Real,twist,0,(Attr::noSave | Attr::readonly),"Elastic twist angle of the contact."))
+ ((Vector3r,bending,Vector3r::Zero(),(Attr::noSave | Attr::readonly),"Bending at contact as a vector defining axis of rotation and angle (angle=norm)."))
+ ,
+ /* extra initializers */,
+ /* ctor */ createIndex(); twist=0;bending=Vector3r::Zero();,
+ /* py */
+ );
+ REGISTER_CLASS_INDEX(ScGeom6D,ScGeom);
+};
+REGISTER_SERIALIZABLE(ScGeom6D);
+
=== modified file 'scripts/test/chained-cylinder-spring.py'
--- scripts/test/chained-cylinder-spring.py 2010-10-16 18:31:17 +0000
+++ scripts/test/chained-cylinder-spring.py 2010-10-20 11:31:49 +0000
@@ -5,7 +5,7 @@
# Experiment beam-like behaviour with chained cylinders + CohFrict connexions
from yade import utils
-young=1.0e5
+young=1.0e6
poisson=5
density=2.60e3
frictionAngle=radians(30)
@@ -28,11 +28,11 @@
## Motion equation
NewtonIntegrator(damping=0.15),
PyRunner(iterPeriod=500,command='history()'),
- PyRunner(iterPeriod=5000,command='if O.iter<21000 : yade.qt.center()')
+ #PyRunner(iterPeriod=5000,command='if O.iter<21000 : yade.qt.center()')
]
#Generate a spiral
-Ne=400
+Ne=200
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;
Follow ups
