yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05943
Re: [Branch ~yade-dev/yade/trunk] Rev 2504: - Add new IG class ScGeom6D precomputing relative rotations, and associated Ig functors.
Merci beaucoup, Bruno! This new class is exactly what I needed and I think
also a good introduction to the software development.
Thanks again. Chiara
On 20 October 2010 12:38, <noreply@xxxxxxxxxxxxx> wrote:
> ------------------------------------------------------------
> 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<https://code.launchpad.net/%7Eyade-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<https://code.launchpad.net/%7Eyade-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;
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev<https://launchpad.net/%7Eyade-dev>
> Post to : yade-dev@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-dev<https://launchpad.net/%7Eyade-dev>
> More help : https://help.launchpad.net/ListHelp
>
>
References
