yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11054
[Branch ~yade-pkg/yade/git-trunk] Rev 4078: Join all InelastCohFrict files into InelastCohFrictPM.
------------------------------------------------------------
revno: 4078
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Mon 2014-07-14 22:08:34 +0200
message:
Join all InelastCohFrict files into InelastCohFrictPM.
removed:
pkg/dem/InelastCohFrictMat.cpp
pkg/dem/InelastCohFrictMat.hpp
pkg/dem/InelastCohFrictPhys.cpp
pkg/dem/InelastCohFrictPhys.hpp
pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.cpp
pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp
pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.cpp
pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp
added:
pkg/dem/InelastCohFrictPM.cpp
pkg/dem/InelastCohFrictPM.hpp
modified:
pkg/dem/Ip2_ElastMat.cpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== removed file 'pkg/dem/InelastCohFrictMat.cpp'
--- pkg/dem/InelastCohFrictMat.cpp 2012-12-14 14:38:44 +0000
+++ pkg/dem/InelastCohFrictMat.cpp 1970-01-01 00:00:00 +0000
@@ -1,7 +0,0 @@
-#include "InelastCohFrictMat.hpp"
-
-InelastCohFrictMat::~InelastCohFrictMat()
-{
-}
-
-YADE_PLUGIN((InelastCohFrictMat));
\ No newline at end of file
=== removed file 'pkg/dem/InelastCohFrictMat.hpp'
--- pkg/dem/InelastCohFrictMat.hpp 2013-06-25 08:02:13 +0000
+++ pkg/dem/InelastCohFrictMat.hpp 1970-01-01 00:00:00 +0000
@@ -1,50 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* 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. *
-*************************************************************************/
-
-
-#pragma once
-
-#include "../common/ElastMat.hpp"
-
-class InelastCohFrictMat : public FrictMat
-{
- public :
- virtual ~InelastCohFrictMat ();
-
-/// Serialization
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(InelastCohFrictMat,FrictMat,"",
- ((Real,tensionModulus,0.0,,"Tension elasticity modulus"))
- ((Real,compressionModulus,0.0,,"Compresion elasticity modulus"))
- ((Real,shearModulus,0.0,,"shear elasticity modulus"))
- ((Real,alphaKr,2.0,,"Dimensionless coefficient used for the rolling stiffness."))
- ((Real,alphaKtw,2.0,,"Dimensionless coefficient used for the twist stiffness."))
-
- ((Real,nuBending,0.0,,"Bending elastic stress limit"))
- ((Real,nuTwist,0.0,,"Twist elastic stress limit"))
- ((Real,sigmaTension,0.0,,"Tension elastic stress limit"))
- ((Real,sigmaCompression,0.0,,"Compression elastic stress limit"))
- ((Real,shearCohesion,0.0,,"Shear elastic stress limit"))
-
- ((Real,creepTension,0.0,,"Tension/compression creeping coefficient. Usual values between 0 and 1."))
- ((Real,creepBending,0.0,,"Bending creeping coefficient. Usual values between 0 and 1."))
- ((Real,creepTwist,0.0,,"Twist creeping coefficient. Usual values between 0 and 1."))
-
- ((Real,unloadTension,0.0,,"Tension/compression plastic unload coefficient. Usual values between 0 and +infinity."))
- ((Real,unloadBending,0.0,,"Bending plastic unload coefficient. Usual values between 0 and +infinity."))
- ((Real,unloadTwist,0.0,,"Twist plastic unload coefficient. Usual values between 0 and +infinity."))
-
- ((Real,epsilonMaxTension,0.0,,"Maximal plastic strain tension"))
- ((Real,epsilonMaxCompression,0.0,,"Maximal plastic strain compression"))
- ((Real,etaMaxBending,0.0,,"Maximal plastic bending strain"))
- ((Real,etaMaxTwist,0.0,,"Maximal plastic twist strain")),
- createIndex();
- );
-/// Indexable
- REGISTER_CLASS_INDEX(InelastCohFrictMat,FrictMat);
-};
-
-REGISTER_SERIALIZABLE(InelastCohFrictMat);
=== added file 'pkg/dem/InelastCohFrictPM.cpp'
--- pkg/dem/InelastCohFrictPM.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/InelastCohFrictPM.cpp 2014-07-14 20:08:34 +0000
@@ -0,0 +1,271 @@
+#include "InelastCohFrictPM.hpp"
+
+YADE_PLUGIN((InelastCohFrictMat)(InelastCohFrictPhys)(Ip2_2xInelastCohFrictMat_InelastCohFrictPhys)(Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment));
+
+
+void Ip2_2xInelastCohFrictMat_InelastCohFrictPhys::go(const shared_ptr<Material>& b1 // InelastCohFrictMat
+ , const shared_ptr<Material>& b2 // InelastCohFrictMat
+ , const shared_ptr<Interaction>& interaction)
+{
+
+ InelastCohFrictMat* sdec1 = static_cast<InelastCohFrictMat*>(b1.get());
+ InelastCohFrictMat* sdec2 = static_cast<InelastCohFrictMat*>(b2.get());
+ ScGeom6D* geom = YADE_CAST<ScGeom6D*>(interaction->geom.get());
+
+ //FIXME : non cohesive contact are not implemented, it would be useful to use setCohesionNow, setCohesionOnNewContacts etc ...
+
+ if (geom) {
+ if (!interaction->phys) {
+ interaction->phys = shared_ptr<InelastCohFrictPhys>(new InelastCohFrictPhys());
+ InelastCohFrictPhys* contactPhysics = YADE_CAST<InelastCohFrictPhys*>(interaction->phys.get());
+ Real pi = 3.14159265;
+ Real r1 = geom->radius1;
+ Real r2 = geom->radius2;
+ Real f1 = sdec1->frictionAngle;
+ Real f2 = sdec2->frictionAngle;
+
+ contactPhysics->tangensOfFrictionAngle = tan(min(f1,f2));
+
+ // harmonic average of modulus
+ contactPhysics->knC = 2.0*sdec1->compressionModulus*r1*sdec2->compressionModulus*r2/(sdec1->compressionModulus*r1+sdec2->compressionModulus*r2);
+ contactPhysics->knT = 2.0*sdec1->tensionModulus*r1*sdec2->tensionModulus*r2/(sdec1->tensionModulus*r1+sdec2->tensionModulus*r2);
+ contactPhysics->ks = 2.0*sdec1->shearModulus*r1*sdec2->shearModulus*r2/(sdec1->shearModulus*r1+sdec2->shearModulus*r2);
+
+ // harmonic average of coeficients for bending and twist coeficients
+ Real AlphaKr = 2.0*sdec1->alphaKr*sdec2->alphaKr/(sdec1->alphaKr+sdec2->alphaKr);
+ Real AlphaKtw = 2.0*sdec1->alphaKtw*sdec2->alphaKtw/(sdec1->alphaKtw+sdec2->alphaKtw);
+
+ contactPhysics->kr = r1*r2*contactPhysics->ks*AlphaKr;
+ contactPhysics->ktw = r1*r2*contactPhysics->ks*AlphaKtw;
+
+ contactPhysics->kTCrp = contactPhysics->knT*min(sdec1->creepTension,sdec2->creepTension);
+ contactPhysics->kRCrp = contactPhysics->kr*min(sdec1->creepBending,sdec2->creepBending);
+ contactPhysics->kTwCrp = contactPhysics->ktw*min(sdec1->creepTwist,sdec2->creepTwist);
+
+ contactPhysics->kRUnld = contactPhysics->kr*min(sdec1->unloadBending,sdec2->unloadBending);
+ contactPhysics->kTUnld = contactPhysics->knT*min(sdec1->unloadTension,sdec2->unloadTension);
+ contactPhysics->kTwUnld = contactPhysics->ktw*min(sdec1->unloadTwist,sdec2->unloadTwist);
+
+ contactPhysics->maxElC = min(sdec1->sigmaCompression,sdec2->sigmaCompression)*pow(min(r2, r1),2);
+ contactPhysics->maxElT = min(sdec1->sigmaTension,sdec2->sigmaTension)*pow(min(r2, r1),2);
+ contactPhysics->maxElB = min(sdec1->nuBending,sdec2->nuBending)*pow(min(r2, r1),3);
+ contactPhysics->maxElTw = min(sdec1->nuTwist,sdec2->nuTwist)*pow(min(r2, r1),3);
+
+ contactPhysics->shearAdhesion = min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(min(r1, r2),2);
+
+ contactPhysics->maxExten = min(sdec1->epsilonMaxTension*r1,sdec2->epsilonMaxTension*r2);
+ contactPhysics->maxContract = min(sdec1->epsilonMaxCompression*r1,sdec2->epsilonMaxCompression*r2);
+
+ contactPhysics->maxBendMom = min(sdec1->etaMaxBending,sdec2->etaMaxBending)*pow(min(r2, r1),3);
+ contactPhysics->maxTwist = 2*pi*min(sdec1->etaMaxTwist,sdec2->etaMaxTwist);
+ }
+ }
+};
+
+
+
+Real Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::normElastEnergy()
+{ //FIXME : this have to be checked and adapted
+ Real normEnergy=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ InelastCohFrictPhys* phys = YADE_CAST<InelastCohFrictPhys*>(I->phys.get());
+ if (phys) {
+ normEnergy += 0.5*(phys->normalForce.squaredNorm()/phys->kn);
+ }
+ }
+ return normEnergy;
+}
+
+Real Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::shearElastEnergy()
+{ //FIXME : this have to be checked and adapted
+ Real shearEnergy=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ InelastCohFrictPhys* phys = YADE_CAST<InelastCohFrictPhys*>(I->phys.get());
+ if (phys) {
+ shearEnergy += 0.5*(phys->shearForce.squaredNorm()/phys->ks);
+ }
+ }
+ return shearEnergy;
+}
+
+
+void Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact)
+{
+//FIXME : non cohesive contact are not implemented, it would be useful to use setCohesionNow, setCohesionOnNewContacts etc ...
+ const int &id1 = contact->getId1();
+ const int &id2 = contact->getId2();
+ const Real& dt = scene->dt;
+ ScGeom6D* geom = YADE_CAST<ScGeom6D*> (ig.get());
+ InelastCohFrictPhys* phys = YADE_CAST<InelastCohFrictPhys*> (ip.get());
+ if (contact->isFresh(scene)) phys->shearForce = Vector3r::Zero();
+
+ Real un = geom->penetrationDepth-phys->unp;
+ Real Fn;
+
+ State* de1 = Body::byId(id1,scene)->state.get();
+ State* de2 = Body::byId(id2,scene)->state.get();
+
+
+ if(un<=0){/// tension ///
+ if(-un>phys->maxExten || phys->isBroken){//plastic failure.
+ phys->isBroken=1;
+ phys->normalForce=phys->shearForce=phys->moment_twist=phys->moment_bending=Vector3r(0,0,0);
+ scene->interactions->requestErase(contact);
+ return;
+ }
+ Fn=phys->knT*un; //elasticity
+ if(-Fn>phys->maxElT || phys->onPlastT){ //so we are on plastic deformation.
+ phys->onPlastT=1;
+ phys->onPlastC=1; //if plasticity is reached on tension, set it to compression too.
+ if(phys->maxCrpRchdT[0]<un){ //unloading/reloading on plastic deformation.
+ Fn = phys->kTUnld*(un-phys->maxCrpRchdT[0])+phys->maxCrpRchdT[1];
+ }
+ else{//loading on plastic deformation : creep.
+ Fn = -phys->maxElT+phys->kTCrp*(un+phys->maxElT/phys->knT);
+ phys->maxCrpRchdT[0]=un; //new maximum is reached.
+ phys->maxCrpRchdT[1]=Fn;
+ }
+ if (Fn>0){ //so the contact just passed the equilibrium state, set new "unp" who stores the plastic equilibrium state.
+ phys->unp=geom->penetrationDepth;
+ phys->maxCrpRchdT[0]=1e20;
+ phys->maxElT=0;
+ }
+ }
+ else{ //elasticity
+ phys->maxCrpRchdT[0]=un;
+ phys->maxCrpRchdT[1]=Fn;
+ }
+ }
+
+ else{/// compression /// similar to tension.
+ if(un>phys->maxContract || phys->isBroken){
+ phys->isBroken=1;
+ phys->normalForce=phys->shearForce=phys->moment_twist=phys->moment_bending=Vector3r(0,0,0);
+ if(geom->penetrationDepth<=0){ //do not erase the contact while penetrationDepth<0 because it would be recreated at next timestep.
+ scene->interactions->requestErase(contact);
+ }
+ return;
+ }
+ Fn=phys->knC*un;
+ if(Fn>phys->maxElC || phys->onPlastC){
+ phys->onPlastC=1;
+ if(phys->maxCrpRchdC[0]>un){
+ Fn = phys->kTUnld*(un-phys->maxCrpRchdC[0])+phys->maxCrpRchdC[1];
+ }
+ else{
+ Fn = phys->maxElC+phys->kTCrp*(un-phys->maxElC/phys->knC);
+ phys->maxCrpRchdC[0]=un;
+ phys->maxCrpRchdC[1]=Fn;
+ }
+ if (Fn<0){
+ phys->unp=geom->penetrationDepth;
+ phys->maxCrpRchdC[0]=-1e20;
+ phys->maxElC=0;
+ }
+ }
+ else{
+ phys->maxCrpRchdC[0]=un;
+ phys->maxCrpRchdC[1]=Fn;
+ }
+ }
+
+ /// Shear ///
+ Vector3r shearForce = geom->rotate(phys->shearForce);
+ const Vector3r& dus = geom->shearIncrement();
+
+ //Linear elasticity giving "trial" shear force
+ shearForce += phys->ks*dus;
+ Real Fs = shearForce.norm();
+ Real maxFs = phys->shearAdhesion;
+ if (maxFs==0)maxFs = Fn*phys->tangensOfFrictionAngle;
+ maxFs = std::max((Real) 0, maxFs);
+ if (Fs > maxFs) {//Plasticity condition on shear force
+ if (!phys->cohesionBroken) {
+ phys->cohesionBroken=1;
+ phys->shearAdhesion=0;
+ maxFs = max((Real) 0, Fn*phys->tangensOfFrictionAngle);
+ }
+ maxFs = maxFs / Fs;
+ shearForce *= maxFs;
+ }
+
+ //rotational moment are only applied if the cohesion is not broken.
+ /// Twist /// the twist law is driven by twist displacement ("getTwist()").
+ if(!phys->cohesionBroken){
+ Real twist = geom->getTwist() - phys->twp;
+ Real twistM=twist*phys->ktw; //elastic twist moment.
+ bool sgnChanged=0; //whether the twist moment just passed the equilibrium state.
+ if(!contact->isFresh(scene) && phys->moment_twist.dot(twistM*geom->normal)<0)sgnChanged=1;
+ if(std::abs(twist)>phys->maxTwist){
+ phys->cohesionBroken=1;
+ twistM=0;
+ }
+ else{
+ if(std::abs(twistM)>phys->maxElTw || phys->onPlastTw){ //plastic deformation.
+ phys->onPlastTw=1;
+ if(std::abs(phys->maxCrpRchdTw[0])>std::abs(twist)){ //unloading/reloading
+ twistM = phys->kTwUnld*(twist-phys->maxCrpRchdTw[0])+phys->maxCrpRchdTw[1];
+ }
+ else{//creep loading.
+ int sign = twist<0?-1:1;
+ twistM = sign*phys->maxElTw+phys->kTwCrp*(twist-sign*phys->maxElTw/phys->ktw); //creep
+ phys->maxCrpRchdTw[0]=twist; //new maximum reached
+ phys->maxCrpRchdTw[1]=twistM;
+ }
+ if(sgnChanged){
+ phys->maxElTw=0;
+ phys->twp=geom->getTwist();
+ phys->maxCrpRchdTw[0]=0;
+ }
+ }
+ else{ //elasticity
+ phys->maxCrpRchdTw[0]=twist;
+ phys->maxCrpRchdTw[1]=twistM;
+ }
+ }
+ phys->moment_twist = twistM * geom->normal;
+ }
+ else phys->moment_twist=Vector3r(0,0,0);
+
+ /// Bending /// incremental form.
+ if(!phys->cohesionBroken){
+ Vector3r bendM = phys->moment_bending;
+ Vector3r relAngVel = geom->getRelAngVel(de1,de2,dt);
+ Vector3r relRotBend = (relAngVel - geom->normal.dot(relAngVel)*geom->normal)*dt; // relative rotation due to rolling behaviour
+ bendM = geom->rotate(phys->moment_bending); // rotate moment vector (updated)
+ phys->pureCreep=geom->rotate(phys->pureCreep); // pure creep is updated to compute the damage.
+ Vector3r bendM_elast = bendM-phys->kr*relRotBend;
+ if(bendM_elast.norm()>phys->maxElB || phys->onPlastB){ // plastic behavior
+ phys->onPlastB=1;
+ bendM=bendM-phys->kDam*relRotBend; //trial bending
+ if(bendM.norm()<phys->moment_bending.norm()){ // if bending decreased, we are unloading ...
+ bendM = bendM+phys->kDam*relRotBend-phys->kRUnld*relRotBend; // ... so undo bendM and apply unload coefficient.
+ Vector3r newPureCreep = phys->pureCreep-phys->kRCrp*relRotBend; // trial pure creep.
+ phys->pureCreep = newPureCreep.norm()<phys->pureCreep.norm()?newPureCreep:phys->pureCreep+phys->kRCrp*relRotBend; // while unloading, pure creep must decrease.
+ phys->kDam=phys->kr+(phys->kRCrp-phys->kr)*(phys->maxCrpRchdB.norm()-phys->maxElB)/(phys->maxBendMom-phys->maxElB); // compute the damage coefficient.
+ }
+ else{ // bending increased, so we are loading (bendM has to be unchanged).
+ Vector3r newPureCreep = phys->pureCreep-phys->kRCrp*relRotBend;
+ phys->pureCreep = newPureCreep.norm()>phys->pureCreep.norm()?newPureCreep:phys->pureCreep+phys->kRCrp*relRotBend; // while loading, pure creep must increase.
+ if(phys->pureCreep.norm()<bendM.norm()) bendM=phys->pureCreep; // bending moment can't be greather than pure creep.
+ if(phys->pureCreep.norm()>phys->maxCrpRchdB.norm()) phys->maxCrpRchdB=phys->pureCreep; // maxCrpRchdB must follow the maximum of pure creep.
+ if(phys->pureCreep.norm()>phys->maxBendMom){
+ phys->cohesionBroken=1;
+ bendM=bendM_elast=Vector3r(0,0,0);
+ }
+ }
+ phys->moment_bending=bendM;
+ }
+ else{//elasticity
+ phys->pureCreep=phys->moment_bending=phys->maxCrpRchdB=bendM_elast;
+ phys->kDam=phys->kRCrp;
+ }
+ }
+ phys->shearForce=shearForce;
+ phys->normalForce=-Fn*geom->normal;
+ applyForceAtContactPoint(phys->normalForce+phys->shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position + (scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): Vector3r::Zero()));
+ scene->forces.addTorque(id1,-phys->moment_bending-phys->moment_twist);
+ scene->forces.addTorque(id2,phys->moment_bending+phys->moment_twist);
+}
=== added file 'pkg/dem/InelastCohFrictPM.hpp'
--- pkg/dem/InelastCohFrictPM.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/InelastCohFrictPM.hpp 2014-07-14 20:08:34 +0000
@@ -0,0 +1,149 @@
+/*************************************************************************
+* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
+* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
+* 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. *
+*************************************************************************/
+
+
+#pragma once
+
+#include <yade/pkg/common/ElastMat.hpp>
+#include <yade/pkg/dem/ScGeom.hpp>
+#include <yade/pkg/dem/FrictPhys.hpp>
+#include "CohesiveFrictionalContactLaw.hpp"
+
+class InelastCohFrictMat : public FrictMat
+{
+ public :
+ virtual ~InelastCohFrictMat () {};
+
+/// Serialization
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(InelastCohFrictMat,FrictMat,"",
+ ((Real,tensionModulus,0.0,,"Tension elasticity modulus"))
+ ((Real,compressionModulus,0.0,,"Compresion elasticity modulus"))
+ ((Real,shearModulus,0.0,,"shear elasticity modulus"))
+ ((Real,alphaKr,2.0,,"Dimensionless coefficient used for the rolling stiffness."))
+ ((Real,alphaKtw,2.0,,"Dimensionless coefficient used for the twist stiffness."))
+
+ ((Real,nuBending,0.0,,"Bending elastic stress limit"))
+ ((Real,nuTwist,0.0,,"Twist elastic stress limit"))
+ ((Real,sigmaTension,0.0,,"Tension elastic stress limit"))
+ ((Real,sigmaCompression,0.0,,"Compression elastic stress limit"))
+ ((Real,shearCohesion,0.0,,"Shear elastic stress limit"))
+
+ ((Real,creepTension,0.0,,"Tension/compression creeping coefficient. Usual values between 0 and 1."))
+ ((Real,creepBending,0.0,,"Bending creeping coefficient. Usual values between 0 and 1."))
+ ((Real,creepTwist,0.0,,"Twist creeping coefficient. Usual values between 0 and 1."))
+
+ ((Real,unloadTension,0.0,,"Tension/compression plastic unload coefficient. Usual values between 0 and +infinity."))
+ ((Real,unloadBending,0.0,,"Bending plastic unload coefficient. Usual values between 0 and +infinity."))
+ ((Real,unloadTwist,0.0,,"Twist plastic unload coefficient. Usual values between 0 and +infinity."))
+
+ ((Real,epsilonMaxTension,0.0,,"Maximal plastic strain tension"))
+ ((Real,epsilonMaxCompression,0.0,,"Maximal plastic strain compression"))
+ ((Real,etaMaxBending,0.0,,"Maximal plastic bending strain"))
+ ((Real,etaMaxTwist,0.0,,"Maximal plastic twist strain")),
+ createIndex();
+ );
+/// Indexable
+ REGISTER_CLASS_INDEX(InelastCohFrictMat,FrictMat);
+};
+
+REGISTER_SERIALIZABLE(InelastCohFrictMat);
+
+class InelastCohFrictPhys : public FrictPhys
+{
+ public :
+ virtual ~InelastCohFrictPhys() {};
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(InelastCohFrictPhys,FrictPhys,"",
+ ((bool,cohesionBroken,false,,"is cohesion active? will be set false when a fragile contact is broken"))
+
+ ((Real,knT,0,,"tension stiffness"))
+ ((Real,knC,0,,"compression stiffness"))
+ ((Real,ktw,0,,"twist shear stiffness"))
+ ((Real,ks,0,,"shear stiffness"))
+ ((Real,kr,0,,"bending stiffness"))
+
+ ((Real,maxElB,0.0,,"Maximum bending elastic moment."))
+ ((Real,maxElTw,0.0,,"Maximum twist elastic moment."))
+ ((Real,maxElT,0.0,,"Maximum tension elastic force."))
+ ((Real,maxElC,0.0,,"Maximum compression elastic force."))
+ ((Real,shearAdhesion,0,,"Maximum elastic shear force (cohesion)."))
+
+ ((Real,kTCrp,0.0,,"Tension/compression creep stiffness"))
+ ((Real,kRCrp,0.0,,"Bending creep stiffness"))
+ ((Real,kTwCrp,0.0,,"Twist creep stiffness"))
+
+ ((Real,kTUnld,0.0,,"Tension/compression plastic unload stiffness"))
+ ((Real,kRUnld,0.0,,"Bending plastic unload stiffness"))
+ ((Real,kTwUnld,0.0,,"Twist plastic unload stiffness"))
+
+ ((Real,maxExten,0.0,,"Plastic failure extension (stretching)."))
+ ((Real,maxContract,0.0,,"Plastic failure contraction (shrinkage)."))
+ ((Real,maxBendMom,0.0,,"Plastic failure bending moment."))
+ ((Real,maxTwist,0.0,,"Plastic failure twist angle"))
+
+ ((bool,isBroken,false,,"true if compression plastic fracture achieved"))
+
+ ((Real,unp,0,,"plastic normal penetration depth describing the equilibrium state."))
+ ((Real,twp,0,,"plastic twist penetration depth describing the equilibrium state."))
+
+ ((bool,onPlastB,false,Attr::readonly,"true if plasticity achieved on bending"))
+ ((bool,onPlastTw,false,Attr::readonly,"true if plasticity achieved on twisting"))
+ ((bool,onPlastT,false,Attr::readonly,"true if plasticity achieved on traction"))
+ ((bool,onPlastC,false,Attr::readonly,"true if plasticity achieved on compression"))
+
+ ((Vector2r,maxCrpRchdT,Vector2r(0,0),Attr::readonly,"maximal extension reached on plastic deformation. maxCrpRchdT[0] stores un and maxCrpRchdT[1] stores Fn."))
+ ((Vector2r,maxCrpRchdC,Vector2r(0,0),Attr::readonly,"maximal compression reached on plastic deformation. maxCrpRchdC[0] stores un and maxCrpRchdC[1] stores Fn."))
+ ((Vector2r,maxCrpRchdTw,Vector2r(0,0),Attr::readonly,"maximal twist reached on plastic deformation. maxCrpRchdTw[0] stores twist angle and maxCrpRchdTw[1] stores twist moment."))
+ ((Vector3r,maxCrpRchdB,Vector3r(0,0,0),Attr::readonly,"maximal bending moment reached on plastic deformation."))
+
+ ((Vector3r,moment_twist,Vector3r(0,0,0),(Attr::readonly),"Twist moment"))
+ ((Vector3r,moment_bending,Vector3r(0,0,0),(Attr::readonly),"Bending moment"))
+ ((Vector3r,pureCreep,Vector3r(0,0,0),(Attr::readonly),"Pure creep curve, used for comparison in calculation."))
+ ((Real,kDam,0,(Attr::readonly),"Damage coefficient on bending, computed from maximum bending moment reached and pure creep behaviour. Its values will vary between :yref:`InelastCohFrictPhys::kr` and :yref:`InelastCohFrictPhys::kRCrp` ."))
+ // internal attributes
+ ,
+ createIndex();
+ );
+/// Indexable
+ REGISTER_CLASS_INDEX(InelastCohFrictPhys,FrictPhys);
+
+};
+
+REGISTER_SERIALIZABLE(InelastCohFrictPhys);
+
+class Ip2_2xInelastCohFrictMat_InelastCohFrictPhys : public IPhysFunctor
+{
+ public :
+ virtual void go(const shared_ptr<Material>& b1,
+ const shared_ptr<Material>& b2,
+ const shared_ptr<Interaction>& interaction);
+ int cohesionDefinitionIteration;
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_2xInelastCohFrictMat_InelastCohFrictPhys,IPhysFunctor,
+ "Generates cohesive-frictional interactions with moments. Used in the contact law :yref:`Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment`.",
+ ,
+ cohesionDefinitionIteration = -1;
+ );
+ FUNCTOR2D(InelastCohFrictMat,InelastCohFrictMat);
+};
+
+REGISTER_SERIALIZABLE(Ip2_2xInelastCohFrictMat_InelastCohFrictPhys);
+
+
+class Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment: public LawFunctor{
+ public:
+ Real normElastEnergy();
+ Real shearElastEnergy();
+ virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment,LawFunctor,"This law is currently under developpement. Final version and documentation will come before the end of 2014.",
+ ,,
+ .def("normElastEnergy",&Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::normElastEnergy,"Compute normal elastic energy.")
+ .def("shearElastEnergy",&Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::shearElastEnergy,"Compute shear elastic energy.")
+ );
+ FUNCTOR2D(ScGeom6D,InelastCohFrictPhys);
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment);
=== removed file 'pkg/dem/InelastCohFrictPhys.cpp'
--- pkg/dem/InelastCohFrictPhys.cpp 2013-06-25 08:02:13 +0000
+++ pkg/dem/InelastCohFrictPhys.cpp 1970-01-01 00:00:00 +0000
@@ -1,13 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* 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. *
-*************************************************************************/
-
-
-#include "InelastCohFrictPhys.hpp"
-InelastCohFrictPhys::~InelastCohFrictPhys()
-{
-}
-YADE_PLUGIN((InelastCohFrictPhys));
\ No newline at end of file
=== removed file 'pkg/dem/InelastCohFrictPhys.hpp'
--- pkg/dem/InelastCohFrictPhys.hpp 2013-06-25 08:02:13 +0000
+++ pkg/dem/InelastCohFrictPhys.hpp 1970-01-01 00:00:00 +0000
@@ -1,74 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* 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. *
-*************************************************************************/
-
-
-#pragma once
-
-#include "FrictPhys.hpp"
-
-
-class InelastCohFrictPhys : public FrictPhys
-{
- public :
- virtual ~InelastCohFrictPhys();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(InelastCohFrictPhys,FrictPhys,"",
- ((bool,cohesionBroken,false,,"is cohesion active? will be set false when a fragile contact is broken"))
-
- ((Real,knT,0,,"tension stiffness"))
- ((Real,knC,0,,"compression stiffness"))
- ((Real,ktw,0,,"twist shear stiffness"))
- ((Real,ks,0,,"shear stiffness"))
- ((Real,kr,0,,"bending stiffness"))
-
- ((Real,maxElB,0.0,,"Maximum bending elastic moment."))
- ((Real,maxElTw,0.0,,"Maximum twist elastic moment."))
- ((Real,maxElT,0.0,,"Maximum tension elastic force."))
- ((Real,maxElC,0.0,,"Maximum compression elastic force."))
- ((Real,shearAdhesion,0,,"Maximum elastic shear force (cohesion)."))
-
- ((Real,kTCrp,0.0,,"Tension/compression creep stiffness"))
- ((Real,kRCrp,0.0,,"Bending creep stiffness"))
- ((Real,kTwCrp,0.0,,"Twist creep stiffness"))
-
- ((Real,kTUnld,0.0,,"Tension/compression plastic unload stiffness"))
- ((Real,kRUnld,0.0,,"Bending plastic unload stiffness"))
- ((Real,kTwUnld,0.0,,"Twist plastic unload stiffness"))
-
- ((Real,maxExten,0.0,,"Plastic failure extension (stretching)."))
- ((Real,maxContract,0.0,,"Plastic failure contraction (shrinkage)."))
- ((Real,maxBendMom,0.0,,"Plastic failure bending moment."))
- ((Real,maxTwist,0.0,,"Plastic failure twist angle"))
-
- ((bool,isBroken,false,,"true if compression plastic fracture achieved"))
-
- ((Real,unp,0,,"plastic normal penetration depth describing the equilibrium state."))
- ((Real,twp,0,,"plastic twist penetration depth describing the equilibrium state."))
-
- ((bool,onPlastB,false,Attr::readonly,"true if plasticity achieved on bending"))
- ((bool,onPlastTw,false,Attr::readonly,"true if plasticity achieved on twisting"))
- ((bool,onPlastT,false,Attr::readonly,"true if plasticity achieved on traction"))
- ((bool,onPlastC,false,Attr::readonly,"true if plasticity achieved on compression"))
-
- ((Vector2r,maxCrpRchdT,Vector2r(0,0),Attr::readonly,"maximal extension reached on plastic deformation. maxCrpRchdT[0] stores un and maxCrpRchdT[1] stores Fn."))
- ((Vector2r,maxCrpRchdC,Vector2r(0,0),Attr::readonly,"maximal compression reached on plastic deformation. maxCrpRchdC[0] stores un and maxCrpRchdC[1] stores Fn."))
- ((Vector2r,maxCrpRchdTw,Vector2r(0,0),Attr::readonly,"maximal twist reached on plastic deformation. maxCrpRchdTw[0] stores twist angle and maxCrpRchdTw[1] stores twist moment."))
- ((Vector3r,maxCrpRchdB,Vector3r(0,0,0),Attr::readonly,"maximal bending moment reached on plastic deformation."))
-
- ((Vector3r,moment_twist,Vector3r(0,0,0),(Attr::readonly),"Twist moment"))
- ((Vector3r,moment_bending,Vector3r(0,0,0),(Attr::readonly),"Bending moment"))
- ((Vector3r,pureCreep,Vector3r(0,0,0),(Attr::readonly),"Pure creep curve, used for comparison in calculation."))
- ((Real,kDam,0,(Attr::readonly),"Damage coefficient on bending, computed from maximum bending moment reached and pure creep behaviour. Its values will vary between :yref:`InelastCohFrictPhys::kr` and :yref:`InelastCohFrictPhys::kRCrp` ."))
- // internal attributes
- ,
- createIndex();
- );
-/// Indexable
- REGISTER_CLASS_INDEX(InelastCohFrictPhys,FrictPhys);
-
-};
-
-REGISTER_SERIALIZABLE(InelastCohFrictPhys);
=== removed file 'pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.cpp'
--- pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.cpp 2013-06-25 08:02:13 +0000
+++ pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.cpp 1970-01-01 00:00:00 +0000
@@ -1,72 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* 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. *
-*************************************************************************/
-
-
-#include "Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp"
-#include<yade/pkg/dem/ScGeom.hpp>
-
-
-
-void Ip2_2xInelastCohFrictMat_InelastCohFrictPhys::go(const shared_ptr<Material>& b1 // InelastCohFrictMat
- , const shared_ptr<Material>& b2 // InelastCohFrictMat
- , const shared_ptr<Interaction>& interaction)
-{
-
- InelastCohFrictMat* sdec1 = static_cast<InelastCohFrictMat*>(b1.get());
- InelastCohFrictMat* sdec2 = static_cast<InelastCohFrictMat*>(b2.get());
- ScGeom6D* geom = YADE_CAST<ScGeom6D*>(interaction->geom.get());
-
- //FIXME : non cohesive contact are not implemented, it would be useful to use setCohesionNow, setCohesionOnNewContacts etc ...
-
- if (geom) {
- if (!interaction->phys) {
- interaction->phys = shared_ptr<InelastCohFrictPhys>(new InelastCohFrictPhys());
- InelastCohFrictPhys* contactPhysics = YADE_CAST<InelastCohFrictPhys*>(interaction->phys.get());
- Real pi = 3.14159265;
- Real r1 = geom->radius1;
- Real r2 = geom->radius2;
- Real f1 = sdec1->frictionAngle;
- Real f2 = sdec2->frictionAngle;
-
- contactPhysics->tangensOfFrictionAngle = tan(min(f1,f2));
-
- // harmonic average of modulus
- contactPhysics->knC = 2.0*sdec1->compressionModulus*r1*sdec2->compressionModulus*r2/(sdec1->compressionModulus*r1+sdec2->compressionModulus*r2);
- contactPhysics->knT = 2.0*sdec1->tensionModulus*r1*sdec2->tensionModulus*r2/(sdec1->tensionModulus*r1+sdec2->tensionModulus*r2);
- contactPhysics->ks = 2.0*sdec1->shearModulus*r1*sdec2->shearModulus*r2/(sdec1->shearModulus*r1+sdec2->shearModulus*r2);
-
- // harmonic average of coeficients for bending and twist coeficients
- Real AlphaKr = 2.0*sdec1->alphaKr*sdec2->alphaKr/(sdec1->alphaKr+sdec2->alphaKr);
- Real AlphaKtw = 2.0*sdec1->alphaKtw*sdec2->alphaKtw/(sdec1->alphaKtw+sdec2->alphaKtw);
-
- contactPhysics->kr = r1*r2*contactPhysics->ks*AlphaKr;
- contactPhysics->ktw = r1*r2*contactPhysics->ks*AlphaKtw;
-
- contactPhysics->kTCrp = contactPhysics->knT*min(sdec1->creepTension,sdec2->creepTension);
- contactPhysics->kRCrp = contactPhysics->kr*min(sdec1->creepBending,sdec2->creepBending);
- contactPhysics->kTwCrp = contactPhysics->ktw*min(sdec1->creepTwist,sdec2->creepTwist);
-
- contactPhysics->kRUnld = contactPhysics->kr*min(sdec1->unloadBending,sdec2->unloadBending);
- contactPhysics->kTUnld = contactPhysics->knT*min(sdec1->unloadTension,sdec2->unloadTension);
- contactPhysics->kTwUnld = contactPhysics->ktw*min(sdec1->unloadTwist,sdec2->unloadTwist);
-
- contactPhysics->maxElC = min(sdec1->sigmaCompression,sdec2->sigmaCompression)*pow(min(r2, r1),2);
- contactPhysics->maxElT = min(sdec1->sigmaTension,sdec2->sigmaTension)*pow(min(r2, r1),2);
- contactPhysics->maxElB = min(sdec1->nuBending,sdec2->nuBending)*pow(min(r2, r1),3);
- contactPhysics->maxElTw = min(sdec1->nuTwist,sdec2->nuTwist)*pow(min(r2, r1),3);
-
- contactPhysics->shearAdhesion = min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(min(r1, r2),2);
-
- contactPhysics->maxExten = min(sdec1->epsilonMaxTension*r1,sdec2->epsilonMaxTension*r2);
- contactPhysics->maxContract = min(sdec1->epsilonMaxCompression*r1,sdec2->epsilonMaxCompression*r2);
-
- contactPhysics->maxBendMom = min(sdec1->etaMaxBending,sdec2->etaMaxBending)*pow(min(r2, r1),3);
- contactPhysics->maxTwist = 2*pi*min(sdec1->etaMaxTwist,sdec2->etaMaxTwist);
- }
- }
-};
-YADE_PLUGIN((Ip2_2xInelastCohFrictMat_InelastCohFrictPhys));
\ No newline at end of file
=== removed file 'pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp'
--- pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp 2014-07-14 20:08:33 +0000
+++ pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp 1970-01-01 00:00:00 +0000
@@ -1,32 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* 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. *
-*************************************************************************/
-
-
-#pragma once
-#include "CohesiveFrictionalContactLaw.hpp"
-#include "InelastCohFrictMat.hpp"
-#include "InelastCohFrictPhys.hpp"
-#include<yade/pkg/dem/ScGeom.hpp>
-
-
-class Ip2_2xInelastCohFrictMat_InelastCohFrictPhys : public IPhysFunctor
-{
- public :
- virtual void go( const shared_ptr<Material>& b1,
- const shared_ptr<Material>& b2,
- const shared_ptr<Interaction>& interaction);
- int cohesionDefinitionIteration;
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_2xInelastCohFrictMat_InelastCohFrictPhys,IPhysFunctor,
- "Generates cohesive-frictional interactions with moments. Used in the contact law :yref:`Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment`.",
- ,
- cohesionDefinitionIteration = -1;
- );
- FUNCTOR2D(InelastCohFrictMat,InelastCohFrictMat);
-};
-
-REGISTER_SERIALIZABLE(Ip2_2xInelastCohFrictMat_InelastCohFrictPhys);
=== modified file 'pkg/dem/Ip2_ElastMat.cpp'
--- pkg/dem/Ip2_ElastMat.cpp 2013-09-23 17:39:59 +0000
+++ pkg/dem/Ip2_ElastMat.cpp 2014-07-14 20:08:34 +0000
@@ -1,7 +1,6 @@
#include "Ip2_ElastMat.hpp"
#include <yade/pkg/common/NormShearPhys.hpp>
-#include <yade/pkg/common/NormShearPhys.hpp>
#include <yade/pkg/dem/DemXDofGeom.hpp>
YADE_PLUGIN((Ip2_ElastMat_ElastMat_NormPhys)(Ip2_ElastMat_ElastMat_NormShearPhys));
=== removed file 'pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.cpp'
--- pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.cpp 1970-01-01 00:00:00 +0000
@@ -1,226 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* 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. *
-*************************************************************************/
-
-
-
-
-#include "Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp"
-
-Real Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::normElastEnergy()
-{ //FIXME : this have to be checked and adapted
- Real normEnergy=0;
- FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
- if(!I->isReal()) continue;
- InelastCohFrictPhys* phys = YADE_CAST<InelastCohFrictPhys*>(I->phys.get());
- if (phys) {
- normEnergy += 0.5*(phys->normalForce.squaredNorm()/phys->kn);
- }
- }
- return normEnergy;
-}
-Real Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::shearElastEnergy()
-{ //FIXME : this have to be checked and adapted
- Real shearEnergy=0;
- FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
- if(!I->isReal()) continue;
- InelastCohFrictPhys* phys = YADE_CAST<InelastCohFrictPhys*>(I->phys.get());
- if (phys) {
- shearEnergy += 0.5*(phys->shearForce.squaredNorm()/phys->ks);
- }
- }
- return shearEnergy;
-}
-
-
-
-void Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact)
-{
-//FIXME : non cohesive contact are not implemented, it would be useful to use setCohesionNow, setCohesionOnNewContacts etc ...
- const int &id1 = contact->getId1();
- const int &id2 = contact->getId2();
- const Real& dt = scene->dt;
- ScGeom6D* geom = YADE_CAST<ScGeom6D*> (ig.get());
- InelastCohFrictPhys* phys = YADE_CAST<InelastCohFrictPhys*> (ip.get());
- if (contact->isFresh(scene)) phys->shearForce = Vector3r::Zero();
-
- Real un = geom->penetrationDepth-phys->unp;
- Real Fn;
-
- State* de1 = Body::byId(id1,scene)->state.get();
- State* de2 = Body::byId(id2,scene)->state.get();
-
-
- if(un<=0){/// tension ///
- if(-un>phys->maxExten || phys->isBroken){//plastic failure.
- phys->isBroken=1;
- phys->normalForce=phys->shearForce=phys->moment_twist=phys->moment_bending=Vector3r(0,0,0);
- scene->interactions->requestErase(contact);
- return;
- }
- Fn=phys->knT*un; //elasticity
- if(-Fn>phys->maxElT || phys->onPlastT){ //so we are on plastic deformation.
- phys->onPlastT=1;
- phys->onPlastC=1; //if plasticity is reached on tension, set it to compression too.
- if(phys->maxCrpRchdT[0]<un){ //unloading/reloading on plastic deformation.
- Fn = phys->kTUnld*(un-phys->maxCrpRchdT[0])+phys->maxCrpRchdT[1];
- }
- else{//loading on plastic deformation : creep.
- Fn = -phys->maxElT+phys->kTCrp*(un+phys->maxElT/phys->knT);
- phys->maxCrpRchdT[0]=un; //new maximum is reached.
- phys->maxCrpRchdT[1]=Fn;
- }
- if (Fn>0){ //so the contact just passed the equilibrium state, set new "unp" who stores the plastic equilibrium state.
- phys->unp=geom->penetrationDepth;
- phys->maxCrpRchdT[0]=1e20;
- phys->maxElT=0;
- }
- }
- else{ //elasticity
- phys->maxCrpRchdT[0]=un;
- phys->maxCrpRchdT[1]=Fn;
- }
- }
-
- else{/// compression /// similar to tension.
- if(un>phys->maxContract || phys->isBroken){
- phys->isBroken=1;
- phys->normalForce=phys->shearForce=phys->moment_twist=phys->moment_bending=Vector3r(0,0,0);
- if(geom->penetrationDepth<=0){ //do not erase the contact while penetrationDepth<0 because it would be recreated at next timestep.
- scene->interactions->requestErase(contact);
- }
- return;
- }
- Fn=phys->knC*un;
- if(Fn>phys->maxElC || phys->onPlastC){
- phys->onPlastC=1;
- if(phys->maxCrpRchdC[0]>un){
- Fn = phys->kTUnld*(un-phys->maxCrpRchdC[0])+phys->maxCrpRchdC[1];
- }
- else{
- Fn = phys->maxElC+phys->kTCrp*(un-phys->maxElC/phys->knC);
- phys->maxCrpRchdC[0]=un;
- phys->maxCrpRchdC[1]=Fn;
- }
- if (Fn<0){
- phys->unp=geom->penetrationDepth;
- phys->maxCrpRchdC[0]=-1e20;
- phys->maxElC=0;
- }
- }
- else{
- phys->maxCrpRchdC[0]=un;
- phys->maxCrpRchdC[1]=Fn;
- }
- }
-
- /// Shear ///
- Vector3r shearForce = geom->rotate(phys->shearForce);
- const Vector3r& dus = geom->shearIncrement();
-
- //Linear elasticity giving "trial" shear force
- shearForce += phys->ks*dus;
- Real Fs = shearForce.norm();
- Real maxFs = phys->shearAdhesion;
- if (maxFs==0)maxFs = Fn*phys->tangensOfFrictionAngle;
- maxFs = std::max((Real) 0, maxFs);
- if (Fs > maxFs) {//Plasticity condition on shear force
- if (!phys->cohesionBroken) {
- phys->cohesionBroken=1;
- phys->shearAdhesion=0;
- maxFs = max((Real) 0, Fn*phys->tangensOfFrictionAngle);
- }
- maxFs = maxFs / Fs;
- shearForce *= maxFs;
- }
-
- //rotational moment are only applied if the cohesion is not broken.
- /// Twist /// the twist law is driven by twist displacement ("getTwist()").
- if(!phys->cohesionBroken){
- Real twist = geom->getTwist() - phys->twp;
- Real twistM=twist*phys->ktw; //elastic twist moment.
- bool sgnChanged=0; //whether the twist moment just passed the equilibrium state.
- if(!contact->isFresh(scene) && phys->moment_twist.dot(twistM*geom->normal)<0)sgnChanged=1;
- if(std::abs(twist)>phys->maxTwist){
- phys->cohesionBroken=1;
- twistM=0;
- }
- else{
- if(std::abs(twistM)>phys->maxElTw || phys->onPlastTw){ //plastic deformation.
- phys->onPlastTw=1;
- if(std::abs(phys->maxCrpRchdTw[0])>std::abs(twist)){ //unloading/reloading
- twistM = phys->kTwUnld*(twist-phys->maxCrpRchdTw[0])+phys->maxCrpRchdTw[1];
- }
- else{//creep loading.
- int sign = twist<0?-1:1;
- twistM = sign*phys->maxElTw+phys->kTwCrp*(twist-sign*phys->maxElTw/phys->ktw); //creep
- phys->maxCrpRchdTw[0]=twist; //new maximum reached
- phys->maxCrpRchdTw[1]=twistM;
- }
- if(sgnChanged){
- phys->maxElTw=0;
- phys->twp=geom->getTwist();
- phys->maxCrpRchdTw[0]=0;
- }
- }
- else{ //elasticity
- phys->maxCrpRchdTw[0]=twist;
- phys->maxCrpRchdTw[1]=twistM;
- }
- }
- phys->moment_twist = twistM * geom->normal;
- }
- else phys->moment_twist=Vector3r(0,0,0);
-
- /// Bending /// incremental form.
- if(!phys->cohesionBroken){
- Vector3r bendM = phys->moment_bending;
- Vector3r relAngVel = geom->getRelAngVel(de1,de2,dt);
- Vector3r relRotBend = (relAngVel - geom->normal.dot(relAngVel)*geom->normal)*dt; // relative rotation due to rolling behaviour
- bendM = geom->rotate(phys->moment_bending); // rotate moment vector (updated)
- phys->pureCreep=geom->rotate(phys->pureCreep); // pure creep is updated to compute the damage.
- Vector3r bendM_elast = bendM-phys->kr*relRotBend;
- if(bendM_elast.norm()>phys->maxElB || phys->onPlastB){ // plastic behavior
- phys->onPlastB=1;
- bendM=bendM-phys->kDam*relRotBend; //trial bending
- if(bendM.norm()<phys->moment_bending.norm()){ // if bending decreased, we are unloading ...
- bendM = bendM+phys->kDam*relRotBend-phys->kRUnld*relRotBend; // ... so undo bendM and apply unload coefficient.
- Vector3r newPureCreep = phys->pureCreep-phys->kRCrp*relRotBend; // trial pure creep.
- phys->pureCreep = newPureCreep.norm()<phys->pureCreep.norm()?newPureCreep:phys->pureCreep+phys->kRCrp*relRotBend; // while unloading, pure creep must decrease.
- phys->kDam=phys->kr+(phys->kRCrp-phys->kr)*(phys->maxCrpRchdB.norm()-phys->maxElB)/(phys->maxBendMom-phys->maxElB); // compute the damage coefficient.
- }
- else{ // bending increased, so we are loading (bendM has to be unchanged).
- Vector3r newPureCreep = phys->pureCreep-phys->kRCrp*relRotBend;
- phys->pureCreep = newPureCreep.norm()>phys->pureCreep.norm()?newPureCreep:phys->pureCreep+phys->kRCrp*relRotBend; // while loading, pure creep must increase.
- if(phys->pureCreep.norm()<bendM.norm()) bendM=phys->pureCreep; // bending moment can't be greather than pure creep.
- if(phys->pureCreep.norm()>phys->maxCrpRchdB.norm()) phys->maxCrpRchdB=phys->pureCreep; // maxCrpRchdB must follow the maximum of pure creep.
- if(phys->pureCreep.norm()>phys->maxBendMom){
- phys->cohesionBroken=1;
- bendM=bendM_elast=Vector3r(0,0,0);
- }
- }
- phys->moment_bending=bendM;
- }
- else{//elasticity
- phys->pureCreep=phys->moment_bending=phys->maxCrpRchdB=bendM_elast;
- phys->kDam=phys->kRCrp;
- }
- }
- phys->shearForce=shearForce;
- phys->normalForce=-Fn*geom->normal;
- applyForceAtContactPoint(phys->normalForce+phys->shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position + (scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): Vector3r::Zero()));
- scene->forces.addTorque(id1,-phys->moment_bending-phys->moment_twist);
- scene->forces.addTorque(id2,phys->moment_bending+phys->moment_twist);
-}
-
-YADE_PLUGIN((Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment));
-
-
-
-
-
-
-
=== removed file 'pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp'
--- pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp 2014-05-26 09:09:48 +0000
+++ pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp 1970-01-01 00:00:00 +0000
@@ -1,26 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* 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. *
-*************************************************************************/
-
-
-#pragma once
-#include "CohesiveFrictionalContactLaw.hpp"
-#include "InelastCohFrictPhys.hpp"
-
-class Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment: public LawFunctor{
- public:
- Real normElastEnergy();
- Real shearElastEnergy();
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
- YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment,LawFunctor,"This law is currently under developpement. Final version and documentation will come before the end of 2014.",
- ,,
- .def("normElastEnergy",&Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::normElastEnergy,"Compute normal elastic energy.")
- .def("shearElastEnergy",&Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::shearElastEnergy,"Compute shear elastic energy.")
- );
- FUNCTOR2D(ScGeom6D,InelastCohFrictPhys);
- DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment);
