yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #04038
[Branch ~yade-dev/yade/trunk] Rev 2158: - Add a "cohesiveFrictional" functor for usage in interaction dispatching (functionally replace t...
------------------------------------------------------------
revno: 2158
committer: Bruno Chareyre <bchareyre@r1arduina>
branch nick: trunk
timestamp: Mon 2010-04-19 15:23:53 +0200
message:
- Add a "cohesiveFrictional" functor for usage in interaction dispatching (functionally replace the global engine CohesiveFrictionalContactLaw, which is still here and
used in cohesive preprocessor).
- Register "cohesiveFrictional" classes and rename them.
- Some cleaning and a fix in the law for brittle failure.
modified:
pkg/dem/DataClass/InteractionPhysics/NormalInelasticityPhys.hpp
pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.cpp
pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.hpp
pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.hpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'pkg/dem/DataClass/InteractionPhysics/NormalInelasticityPhys.hpp'
--- pkg/dem/DataClass/InteractionPhysics/NormalInelasticityPhys.hpp 2010-04-04 22:05:37 +0000
+++ pkg/dem/DataClass/InteractionPhysics/NormalInelasticityPhys.hpp 2010-04-19 13:23:53 +0000
@@ -12,7 +12,7 @@
/*! \brief Interaction for using Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity
-This interaction is similar to CohesiveFrictionalContactInteraction. Among the differences are the unMax, previousun and previousFn (allowing to describe the inelastic unloadings in compression), no more shear and tension Adhesion, no more "fragile", "cohesionBroken" and "cohesionDisablesFriction"
+This interaction is similar to CohFrictPhys. Among the differences are the unMax, previousun and previousFn (allowing to describe the inelastic unloadings in compression), no more shear and tension Adhesion, no more "fragile", "cohesionBroken" and "cohesionDisablesFriction"
*/
class NormalInelasticityPhys : public FrictPhys
=== modified file 'pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.cpp'
--- pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.cpp 2010-02-26 23:29:53 +0000
+++ pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.cpp 2010-04-19 13:23:53 +0000
@@ -9,7 +9,7 @@
#include"Ip2_2xCohFrictMat_NormalInelasticityPhys.hpp"
#include<yade/pkg-dem/ScGeom.hpp>
#include<yade/pkg-dem/NormalInelasticityPhys.hpp>
-#include<yade/pkg-dem/CohesiveFrictionalMat.hpp>
+#include<yade/pkg-dem/CohFrictMat.hpp>
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
@@ -23,12 +23,12 @@
//
//
-void Ip2_2xCohFrictMat_NormalInelasticityPhys::go( const shared_ptr<Material>& b1 // CohesiveFrictionalMat
- , const shared_ptr<Material>& b2 // CohesiveFrictionalMat
+void Ip2_2xCohFrictMat_NormalInelasticityPhys::go( const shared_ptr<Material>& b1 // CohFrictMat
+ , const shared_ptr<Material>& b2 // CohFrictMat
, const shared_ptr<Interaction>& interaction)
{
- CohesiveFrictionalMat* sdec1 = static_cast<CohesiveFrictionalMat*>(b1.get());
- CohesiveFrictionalMat* sdec2 = static_cast<CohesiveFrictionalMat*>(b2.get());
+ CohFrictMat* sdec1 = static_cast<CohFrictMat*>(b1.get());
+ CohFrictMat* sdec2 = static_cast<CohFrictMat*>(b2.get());
ScGeom* interactionGeometry = YADE_CAST<ScGeom*>(interaction->interactionGeometry.get());
//Create cohesive interractions only once
=== modified file 'pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.hpp'
--- pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.hpp 2010-03-23 14:19:56 +0000
+++ pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.hpp 2010-04-19 13:23:53 +0000
@@ -27,7 +27,7 @@
int cohesionDefinitionIteration; //useful is you want to use setCohesionNow
- FUNCTOR2D(CohesiveFrictionalMat,CohesiveFrictionalMat);
+ FUNCTOR2D(CohFrictMat,CohFrictMat);
YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_2xCohFrictMat_NormalInelasticityPhys,
InteractionPhysicsFunctor,
"The RelationShips for using Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity\n\n In these RelationShips all the attributes of the interactions (which are of NormalInelasticityPhys type) are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<InteractionPhysicsFunctor>`, most of the attributes are computed only once, when the interaction is new.",
=== modified file 'pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp'
--- pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp 2010-04-08 09:14:11 +0000
+++ pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp 2010-04-19 13:23:53 +0000
@@ -13,6 +13,7 @@
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
#include<yade/pkg-common/ElastMat.hpp>
+#include <cassert>
@@ -21,36 +22,38 @@
, const shared_ptr<Interaction>& interaction)
{
if(interaction->interactionPhysics) return;
-
- const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(b1);
- const shared_ptr<FrictMat>& mat2 = YADE_PTR_CAST<FrictMat>(b2);
- if (!interaction->interactionPhysics)
- interaction->interactionPhysics = shared_ptr<FrictPhys>(new FrictPhys());
- const shared_ptr<FrictPhys>& contactPhysics = YADE_PTR_CAST<FrictPhys>(interaction->interactionPhysics);
-
- Real Ea = mat1->young;
- Real Eb = mat2->young;
- Real Va = mat1->poisson;
- Real Vb = mat2->poisson;
-
- Real Da,Db; Vector3r normal;
- //FIXME : dynamic casts here???!!!
- ScGeom* scg=dynamic_cast<ScGeom*>(interaction->interactionGeometry.get());
- Dem3DofGeom* d3dg=dynamic_cast<Dem3DofGeom*>(interaction->interactionGeometry.get());
- if(scg){Da=scg->radius1; Db=scg->radius2; normal=scg->normal;}
- else if(d3dg){Da=d3dg->refR1>0?d3dg->refR1:2*d3dg->refR2; Db=d3dg->refR2>0?d3dg->refR2:d3dg->refR1; normal=d3dg->normal;}
- else throw runtime_error("Ip2_FrictMat_FrictMat_FrictPhys: geometry is neither ScGeom nor Dem3DofGeom");
- //harmonic average of the two stiffnesses when (Di.Ei/2) is the stiffness of sphere "i"
- Real Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);
- //same for shear stiffness
- Real Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);
-
- contactPhysics->frictionAngle = std::min(mat1->frictionAngle,mat2->frictionAngle);
- contactPhysics->tangensOfFrictionAngle = std::tan(contactPhysics->frictionAngle);
- contactPhysics->prevNormal = normal;
- contactPhysics->kn = Kn;
- contactPhysics->ks = Ks;
- return;
+ const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(b1);
+ const shared_ptr<FrictMat>& mat2 = YADE_PTR_CAST<FrictMat>(b2);
+ if (!interaction->interactionPhysics)
+ interaction->interactionPhysics = shared_ptr<FrictPhys>(new FrictPhys());
+ const shared_ptr<FrictPhys>& contactPhysics = YADE_PTR_CAST<FrictPhys>(interaction->interactionPhysics);
+ Real Ea = mat1->young;
+ Real Eb = mat2->young;
+ Real Va = mat1->poisson;
+ Real Vb = mat2->poisson;
+
+ Real Da,Db; Vector3r normal;
+ //FIXME : dynamic casts here???!!!
+// ScGeom* scg=dynamic_cast<ScGeom*>(interaction->interactionGeometry.get());
+// Dem3DofGeom* d3dg=dynamic_cast<Dem3DofGeom*>(interaction->interactionGeometry.get());
+// if(scg){Da=scg->radius1; Db=scg->radius2; normal=scg->normal;}
+// else if(d3dg){Da=d3dg->refR1>0?d3dg->refR1:2*d3dg->refR2; Db=d3dg->refR2>0?d3dg->refR2:d3dg->refR1; normal=d3dg->normal;}
+// else throw runtime_error("Ip2_FrictMat_FrictMat_FrictPhys: geometry is neither ScGeom nor Dem3DofGeom");
+
+ assert(dynamic_cast<GenericSpheresContact*>(interaction->interactionGeometry.get()));//only in debug mode
+ GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->interactionGeometry.get());
+ {Da=sphCont->refR1>0?sphCont->refR1:sphCont->refR2; Db=sphCont->refR2>0?sphCont->refR2:sphCont->refR1; normal=sphCont->normal;}
+
+ //harmonic average of the two stiffnesses when (Di.Ei/2) is the stiffness of a contact point on sphere "i"
+ Real Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);
+ //same for shear stiffness
+ Real Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);
+
+ contactPhysics->frictionAngle = std::min(mat1->frictionAngle,mat2->frictionAngle);
+ contactPhysics->tangensOfFrictionAngle = std::tan(contactPhysics->frictionAngle);
+ contactPhysics->prevNormal = normal;
+ contactPhysics->kn = Kn;
+ contactPhysics->ks = Ks;
};
YADE_PLUGIN((Ip2_FrictMat_FrictMat_FrictPhys));
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-03-20 12:40:44 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-04-19 13:23:53 +0000
@@ -7,19 +7,22 @@
*************************************************************************/
#include "CohesiveFrictionalContactLaw.hpp"
-#include<yade/pkg-dem/CohesiveFrictionalMat.hpp>
+#include<yade/pkg-dem/CohFrictMat.hpp>
#include<yade/pkg-dem/ScGeom.hpp>
-#include<yade/pkg-dem/CohesiveFrictionalContactInteraction.hpp>
+#include<yade/pkg-dem/CohFrictPhys.hpp>
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
+
+YADE_PLUGIN((CohesiveFrictionalContactLaw)(Law2_ScGeom_CohFrictPhys_ElasticPlastic));
+CREATE_LOGGER(Law2_ScGeom_CohFrictPhys_ElasticPlastic);
+
Vector3r translation_vect_ ( 0.10,0,0 );
-
+/*
CohesiveFrictionalContactLaw::CohesiveFrictionalContactLaw() : GlobalEngine()
{
sdecGroupMask=1;
- momentRotationLaw = true;
erosionActivated = false;
detectBrokenBodies = false;
always_use_moment_law = false;
@@ -28,7 +31,7 @@
shear_creep=false;
twist_creep=false;
creep_viscosity = 1.0;
-}
+}*/
void out ( Quaternionr q )
@@ -44,196 +47,145 @@
std::cout << " axis: " << axis[0] << " " << axis[1] << " " << axis[2] << ", length: " << axis.Length() << " | ";
}
-void CohesiveFrictionalContactLaw::action ()
-{
+
+
+void CohesiveFrictionalContactLaw::action()
+{
+ if(!functor) functor=shared_ptr<Law2_ScGeom_CohFrictPhys_ElasticPlastic>(new Law2_ScGeom_CohFrictPhys_ElasticPlastic);
+ functor->sdecGroupMask=sdecGroupMask;
+ functor->erosionActivated = erosionActivated;
+ functor->detectBrokenBodies = detectBrokenBodies;
+ functor->always_use_moment_law = always_use_moment_law;
+ functor->shear_creep=shear_creep;
+ functor->twist_creep=twist_creep;
+ functor->creep_viscosity = creep_viscosity;
+
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ functor->go(I->interactionGeometry, I->interactionPhysics, I.get(), scene);
+ }
+}
+
+
+void Law2_ScGeom_CohFrictPhys_ElasticPlastic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* ncb)
+{
+
shared_ptr<BodyContainer>& bodies = scene->bodies;
- Real dt = Omega::instance().getTimeStep();
- InteractionContainer::iterator ii = scene->interactions->begin();
- InteractionContainer::iterator iiEnd = scene->interactions->end();
- for ( ; ii!=iiEnd ; ++ii )
- {
- //if ((*ii)->interactionGeometry && (*ii)->interactionPhysics)
- if ( ( *ii )->isReal() )
- {
- if ( detectBrokenBodies
- && (*bodies)[(*ii)->getId1()]->shape->getClassName() != "box"
- && (*bodies)[(*ii)->getId2()]->shape->getClassName() != "box"
- )
- {
- YADE_CAST<CohesiveFrictionalMat*> ( ( *bodies ) [ ( *ii )->getId1() ]->material.get() )->isBroken = false;
- YADE_CAST<CohesiveFrictionalMat*> ( ( *bodies ) [ ( *ii )->getId2() ]->material.get() )->isBroken = false;
- }
-
- const shared_ptr<Interaction>& contact = *ii;
- int id1 = contact->getId1();
- int id2 = contact->getId2();
-
- if (!((*bodies)[id1]->getGroupMask() & (*bodies)[id2]->getGroupMask() & sdecGroupMask))
- continue; // skip other groups,
- Body* b1 = ( *bodies ) [id1].get();
- Body* b2 = ( *bodies ) [id2].get();
- ScGeom* currentContactGeometry = YADE_CAST<ScGeom*> ( contact->interactionGeometry.get() );
- CohesiveFrictionalContactInteraction* currentContactPhysics = YADE_CAST<CohesiveFrictionalContactInteraction*> ( contact->interactionPhysics.get() );
- Vector3r& shearForce = currentContactPhysics->shearForce;
- if ( contact->isFresh ( scene ) ) shearForce = Vector3r::ZERO;
-
- Real un = currentContactGeometry->penetrationDepth;
- Real Fn = currentContactPhysics->kn*un;
- currentContactPhysics->normalForce = Fn*currentContactGeometry->normal;
- if ( un < 0
- && ( currentContactPhysics->normalForce.SquaredLength() > pow ( currentContactPhysics->normalAdhesion,2 )
- || currentContactPhysics->normalAdhesion==0
- )
- )
- {
- // BREAK due to tension
- scene->interactions->requestErase ( contact->getId1(),contact->getId2() );
- // contact->interactionPhysics was reset now; currentContactPhysics still hold the object, but is not associated with the interaction anymore
- currentContactPhysics->cohesionBroken = true;
- currentContactPhysics->normalForce = Vector3r::ZERO;
- currentContactPhysics->shearForce = Vector3r::ZERO;
- }
- else
- {
- Vector3r axis;
- Real angle;
- /// Here is the code with approximated rotations ///
- axis = currentContactPhysics->prevNormal.Cross ( currentContactGeometry->normal );
- shearForce -= shearForce.Cross ( axis );
- angle = dt*0.5*currentContactGeometry->normal.Dot ( Body::byId ( id1 )->state->angVel+Body::byId ( id2 )->state->angVel );
- axis = angle*currentContactGeometry->normal;
- shearForce -= shearForce.Cross ( axis );
- Vector3r x = currentContactGeometry->contactPoint;
- Vector3r c1x = ( x - b1->state->pos );
- Vector3r c2x = ( x - b2->state->pos );
- /// The following definition of c1x and c2x is to avoid "granular ratcheting"
- /// (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN,
- /// "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena",
- /// ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors)
- Vector3r _c1x_ = currentContactGeometry->radius1*currentContactGeometry->normal;
- Vector3r _c2x_ = -currentContactGeometry->radius2*currentContactGeometry->normal;
- Vector3r relativeVelocity = ( b2->state->vel+b2->state->angVel.Cross ( _c2x_ ) ) - ( b1->state->vel+b1->state->angVel.Cross ( _c1x_ ) );
- Vector3r shearVelocity = relativeVelocity-currentContactGeometry->normal.Dot ( relativeVelocity ) *currentContactGeometry->normal;
- Vector3r shearDisplacement = shearVelocity*dt;
-
-
-///////////////////////// CREEP START (commented out) ///////////
- if ( shear_creep )
- {
- shearForce -= currentContactPhysics->ks* ( shearDisplacement + shearForce*dt/creep_viscosity );
- }
- else
- {
- shearForce -= currentContactPhysics->ks*shearDisplacement;
- }
-///////////////////////// CREEP END /////////////////////////////
-
- // cerr << "shearForce = " << shearForce << endl;
- Real maxFs = 0;
- Real Fs = currentContactPhysics->shearForce.Length();
- maxFs = std::max ( ( Real ) 0, currentContactPhysics->shearAdhesion + Fn*currentContactPhysics->tangensOfFrictionAngle );
- // if (!currentContactPhysics->cohesionDisablesFriction)
- // maxFs += Fn * currentContactPhysics->tangensOfFrictionAngle;
- //} else
- // maxFs = Fn * currentContactPhysics->tangensOfFrictionAngle;
- // cerr << "maxFs = " << maxFs << " Fs = " << Fs<< endl;
- if ( Fs > maxFs )
- {
- currentContactPhysics->cohesionBroken = true;
-
- // brokenStatus[
- //if (currentContactPhysics->fragile && !currentContactPhysics->cohesionBroken)
-
+ const Real dt = Omega::instance().getTimeStep();
+
+ if (contact->isReal()) {
+ 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 = (*bodies)[id1].get();
+ Body* b2 = (*bodies)[id2].get();
+ ScGeom* currentContactGeometry = YADE_CAST<ScGeom*> (ig.get());
+ CohFrictPhys* currentContactPhysics = YADE_CAST<CohFrictPhys*> (ip.get());
+ Vector3r& shearForce = currentContactPhysics->shearForce;
+ if (contact->isFresh(scene)) shearForce = Vector3r::ZERO;
+
+ Real un = currentContactGeometry->penetrationDepth;
+ Real Fn = currentContactPhysics->kn*un;
+ currentContactPhysics->normalForce = Fn*currentContactGeometry->normal;
+ if (un < 0 && (currentContactPhysics->normalForce.SquaredLength() > pow(currentContactPhysics->normalAdhesion,2)
+ || currentContactPhysics->normalAdhesion==0) ) {
+ // BREAK due to tension
+ scene->interactions->requestErase(contact->getId1(),contact->getId2());
+ // contact->interactionPhysics was reset now; currentContactPhysics still hold the object, but is not associated with the interaction anymore
+ currentContactPhysics->cohesionBroken = true;
+ currentContactPhysics->normalForce = Vector3r::ZERO;
+ currentContactPhysics->shearForce = Vector3r::ZERO;
+ } else {
+ State* de1 = Body::byId(id1,ncb)->state.get();
+ State* de2 = Body::byId(id2,ncb)->state.get();
+ ///////////////////////// CREEP START ///////////
+ if (shear_creep) shearForce -= currentContactPhysics->ks*(shearForce*dt/creep_viscosity);
+ ///////////////////////// CREEP END ////////////
+ currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);
+
+ Real Fs = currentContactPhysics->shearForce.Length();
+ Real maxFs = currentContactPhysics->shearAdhesion;
+ if (!currentContactPhysics->cohesionDisablesFriction || maxFs==0)
+ maxFs += Fn*currentContactPhysics->tangensOfFrictionAngle;
+ maxFs = std::max((Real) 0, maxFs);
+
+ if (Fs > maxFs) {//Plasticity condition on shear force
+ if (currentContactPhysics->fragile && !currentContactPhysics->cohesionBroken) {
currentContactPhysics->SetBreakingState();
- maxFs = max ( ( Real ) 0, Fn * currentContactPhysics->tangensOfFrictionAngle );
- maxFs = maxFs / Fs;
- // cerr << "maxFs = " << maxFs << endl;
- if ( maxFs>1 )
- cerr << "maxFs>1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
- shearForce *= maxFs;
- if ( Fn<0 ) currentContactPhysics->normalForce = Vector3r::ZERO;
+ maxFs = max((Real) 0, Fn*currentContactPhysics->tangensOfFrictionAngle);}
+ maxFs = maxFs / Fs;
+ if (maxFs>1) cerr << "maxFs>1!!" << endl;
+ shearForce *= maxFs;
+ if (Fn<0) currentContactPhysics->normalForce = Vector3r::ZERO;}
+
+ applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position, ncb);
+
+ /// 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;
}
- ////////// PFC3d SlipModel
- Vector3r f = currentContactPhysics->normalForce + shearForce;
- // cerr << "currentContactPhysics->normalForce= " << currentContactPhysics->normalForce << endl;
- // cerr << "shearForce " << shearForce << endl;
- // cerr << "f= " << f << endl;
- // it will be some macro( body->physicalActions, ActionType , bodyId )
- scene->forces.addForce ( id1,-f );
- scene->forces.addForce ( id2, f );
- scene->forces.addTorque ( id1,-c1x.Cross ( f ) );
- scene->forces.addTorque ( id2, c2x.Cross ( f ) );
-
-
-
- /// Moment law ///
- if ( momentRotationLaw && ( currentContactPhysics->cohesionBroken == false || 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;
- }
-
- Vector3r axis; // axis of rotation - this is the Moment direction UNIT vector.
- Real angle; // angle represents the power of resistant ELASTIC moment
- delta.ToAxisAngle ( axis,angle );
- if ( angle > Mathr::PI ) angle -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi
-
- //This indentation is a rewrite of original equations (the two commented lines), should work exactly the same.
+
+ Vector3r axis; // axis of rotation - this is the Moment direction UNIT vector.
+ Real angle; // angle represents the power of resistant ELASTIC moment
+ delta.ToAxisAngle(axis,angle);
+ if (angle > Mathr::PI) angle -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi
+
+ //This indentation is a rewrite of original equations (the two commented lines), should work exactly the same.
//Real elasticMoment = currentContactPhysics->kr * std::abs(angle); // positive value (*)
- Real angle_twist ( angle * 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 ( currentContactGeometry->normal , angle_twist );
- //Quaternionr q_twist_creeped(currentContactGeometry->normal , angle_twist*0.996);
- Quaternionr q_twist_creeped ( currentContactGeometry->normal , angle_twist_creeped );
- 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();
- }
-
-
- Vector3r moment_twist ( axis_twist * currentContactPhysics->kr );
-
- Vector3r axis_bending ( angle*axis - axis_twist );
- Vector3r moment_bending ( axis_bending * currentContactPhysics->kr );
-
-//Vector3r moment = axis * elasticMoment * (angle<0.0?-1.0:1.0); // restore sign. (*)
-
- Vector3r moment = moment_twist + moment_bending;
- currentContactPhysics->moment_twist = moment_twist;
- currentContactPhysics->moment_bending = moment_bending;
- scene->forces.addTorque ( id1,-moment );
- scene->forces.addTorque ( id2, moment );
+ Real angle_twist(angle * 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(currentContactGeometry->normal , angle_twist);
+ //Quaternionr q_twist_creeped(currentContactGeometry->normal , angle_twist*0.996);
+ Quaternionr q_twist_creeped(currentContactGeometry->normal , angle_twist_creeped);
+ 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();
}
- /// Moment law END ///
-
- currentContactPhysics->prevNormal = currentContactGeometry->normal;
+ Vector3r moment_twist(axis_twist * currentContactPhysics->kr);
+ Vector3r axis_bending(angle*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;
+ scene->forces.addTorque(id1,-moment);
+ scene->forces.addTorque(id2, moment);
}
+ /// Moment law END ///
+ currentContactPhysics->prevNormal = currentContactGeometry->normal;
}
}
}
-YADE_PLUGIN ((CohesiveFrictionalContactLaw));
+
+
+
+
+
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.hpp 2010-03-20 12:40:44 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.hpp 2010-04-19 13:23:53 +0000
@@ -9,27 +9,51 @@
#pragma once
#include<yade/core/GlobalEngine.hpp>
-
+#include<yade/pkg-common/LawFunctor.hpp>
#include <set>
#include <boost/tuple/tuple.hpp>
+class Law2_ScGeom_CohFrictPhys_ElasticPlastic: public LawFunctor{
+ public:
+ virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene*);
+ YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_CohFrictPhys_ElasticPlastic,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. Can be elastic-fragile or perfectly elastic-plastic. Creep at contact can be enabled.\n\n.. note::\n This law uses :yref:`ScGeom`.",
+ ((int,sdecGroupMask,1,"Bitmask for allowing collision between particles :yref:`Body::groupMask`"))
+ ((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:`CohesiveFrictionalContactLaw::always_use_moment_law` 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."))
+ ((bool,shear_creep,false,"activate creep on the shear force, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
+ ((bool,twist_creep,false,"activate creep on the twisting moment, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
+ ((Real,creep_viscosity,false,"creep viscosity [Pa.s/m]. probably should be moved to Ip2_2xCohFrictMat_CohFrictPhys..."))
+ ((bool,erosionActivated,false,""))
+ ((bool,detectBrokenBodies,false,""))
+
+ );
+ FUNCTOR2D(ScGeom,CohFrictPhys);
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Law2_ScGeom_CohFrictPhys_ElasticPlastic);
+
class CohesiveFrictionalContactLaw : public GlobalEngine
{
-/// Attributes
- public :
- int sdecGroupMask;
- bool momentRotationLaw, erosionActivated, detectBrokenBodies,always_use_moment_law;
- bool shear_creep,twist_creep;
- Real creep_viscosity; /// probably should be moved to CohesiveFrictionalRelationships...
+ shared_ptr<Law2_ScGeom_CohFrictPhys_ElasticPlastic> functor;
+
+ public :
long iter;/// used for checking if new iteration
-
- CohesiveFrictionalContactLaw();
void action();
-
- REGISTER_ATTRIBUTES(GlobalEngine,(sdecGroupMask)(momentRotationLaw)(erosionActivated)(detectBrokenBodies)(always_use_moment_law)(shear_creep)(twist_creep)(creep_viscosity));
- REGISTER_CLASS_NAME(CohesiveFrictionalContactLaw);
- REGISTER_BASE_CLASS_NAME(GlobalEngine);
+
+ YADE_CLASS_BASE_DOC_ATTRS(CohesiveFrictionalContactLaw,GlobalEngine,"[DEPRECATED] Loop over interactions applying :yref:`Law2_ScGeom_CohFrictPhys_ElasticPlastic` on all interactions.\n\n.. note::\n Use :yref:`InteractionDispatchers` and :yref:`Law2_ScGeom_CohFrictPhys_ElasticPlastic` instead of this class for performance reasons.",
+ ((int,sdecGroupMask,1,"Bitmask for allowing collision between particles :yref:`Body::groupMask`"))
+ ((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:`CohesiveFrictionalContactLaw::always_use_moment_law` 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."))
+ ((bool,shear_creep,false,"activate creep on the shear force, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
+ ((bool,twist_creep,false,"activate creep on the twisting moment, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
+ ((Real,creep_viscosity,false,"creep viscosity [Pa.s/m]. probably should be moved to Ip2_2xCohFrictMat_CohFrictPhys..."))
+ ((bool,erosionActivated,false,""))
+ ((bool,detectBrokenBodies,false,""))
+
+ );
};
REGISTER_SERIALIZABLE(CohesiveFrictionalContactLaw);
