yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #04141
[Branch ~yade-dev/yade/trunk] Rev 2182: - Fix cohesiveFrictional crasher
------------------------------------------------------------
revno: 2182
committer: Bruno Chareyre <bchareyre@r1arduina>
branch nick: trunk
timestamp: Mon 2010-04-26 11:02:51 +0200
message:
- Fix cohesiveFrictional crasher
- Initialise saveSImulation correctly in TCE
- Remove initialKn/Ks assignment Ip2_2xCohFrictMat_CohFrictPhys
modified:
pkg/dem/DataClass/InteractionPhysics/FrictPhys.hpp
pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_CohFrictPhys.cpp
pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.hpp
pkg/dem/Engine/PartialEngine/TriaxialCompressionEngine.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/FrictPhys.hpp'
--- pkg/dem/DataClass/InteractionPhysics/FrictPhys.hpp 2010-04-25 13:18:11 +0000
+++ pkg/dem/DataClass/InteractionPhysics/FrictPhys.hpp 2010-04-26 09:02:51 +0000
@@ -1,4 +1,10 @@
-// © 2004 Olivier Galizzi <olivier.galizzi@xxxxxxx>
+/*************************************************************************
+* Copyright (C) 2007 by Bruno CHAREYRE *
+* bruno.chareyre@xxxxxxxxxxx *
+* *
+* 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/NormShearPhys.hpp>
=== modified file 'pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_CohFrictPhys.cpp'
--- pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_CohFrictPhys.cpp 2010-04-22 16:40:45 +0000
+++ pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_CohFrictPhys.cpp 2010-04-26 09:02:51 +0000
@@ -50,9 +50,6 @@
// Victor Donze, "Calibration procedure for spherical
// discrete elements using a local moemnt law".
Real Kr = Da*Db*Ks*2.0; // just like "2.0" above - it's an arbitrary parameter
-
- contactPhysics->initialKn = Kn;
- contactPhysics->initialKs = Ks;
contactPhysics->frictionAngle = std::min(fa,fb);
contactPhysics->tangensOfFrictionAngle = std::tan(contactPhysics->frictionAngle);
@@ -73,8 +70,8 @@
contactPhysics->orientationToContact2 = contactPhysics->initialOrientation2.Conjugate() * contactPhysics->initialContactOrientation;
}
contactPhysics->prevNormal = interactionGeometry->normal;
- contactPhysics->kn = contactPhysics->initialKn;
- contactPhysics->ks = contactPhysics->initialKs;
+ 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;
@@ -89,8 +86,6 @@
else // !isNew, but if setCohesionNow, all contacts are initialized like if they were newly created
{
CohFrictPhys* contactPhysics = YADE_CAST<CohFrictPhys*>(interaction->interactionPhysics.get());
- contactPhysics->kn = contactPhysics->initialKn;
- contactPhysics->ks = contactPhysics->initialKs;
if (setCohesionNow && sdec1->isCohesive && sdec2->isCohesive)
{
contactPhysics->cohesionBroken = false;
@@ -100,10 +95,6 @@
contactPhysics->initialOrientation2 = Body::byId(interaction->getId2())->state->ori;
contactPhysics->initialPosition1 = Body::byId(interaction->getId1())->state->pos;
contactPhysics->initialPosition2 = Body::byId(interaction->getId2())->state->pos;
- Real Da = interactionGeometry->radius1;
- Real Db = interactionGeometry->radius2;
- Real Kr = Da*Db*contactPhysics->ks*2.0; // just like "2.0" above - it's an arbitrary parameter
- contactPhysics->kr = Kr;
contactPhysics->initialContactOrientation.Align(Vector3r(1.0,0.0,0.0),interactionGeometry->normal);
contactPhysics->currentContactOrientation = contactPhysics->initialContactOrientation;
contactPhysics->orientationToContact1 = contactPhysics->initialOrientation1.Conjugate() * contactPhysics->initialContactOrientation;
=== modified file 'pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp'
--- pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp 2010-04-19 13:23:53 +0000
+++ pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp 2010-04-26 09:02:51 +0000
@@ -24,22 +24,14 @@
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());
+ 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");
-
+ Real Da,Db; Vector3r normal;
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;}
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-04-25 13:18:11 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-04-26 09:02:51 +0000
@@ -19,21 +19,6 @@
Vector3r translation_vect_ ( 0.10,0,0 );
-/*
-CohesiveFrictionalContactLaw::CohesiveFrictionalContactLaw() : GlobalEngine()
-{
- sdecGroupMask=1;
- erosionActivated = false;
- detectBrokenBodies = false;
- always_use_moment_law = false;
-
-//CREEP
- shear_creep=false;
- twist_creep=false;
- creep_viscosity = 1.0;
-}*/
-
-
void out ( Quaternionr q )
{
Vector3r axis;
@@ -69,123 +54,106 @@
void Law2_ScGeom_CohFrictPhys_ElasticPlastic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* ncb)
{
-
- shared_ptr<BodyContainer>& bodies = scene->bodies;
-
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;
- 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;
- }
-
- 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 = 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;
- }
+// 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,ncb).get();
+ Body* b2 = Body::byId(id2,ncb).get();
+ ScGeom* currentContactGeometry = YADE_CAST<ScGeom*> (ig.get());
+ CohFrictPhys* currentContactPhysics = YADE_CAST<CohFrictPhys*> (ip.get());
+
+ Vector3r& shearForce = currentContactPhysics->shearForce;
+
+ if (contact->isFresh(ncb)) 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
+ ncb->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;
+ if (maxFs>1) cerr << "maxFs>1!!" << endl;
+ shearForce *= maxFs;
+ if (Fn<0) currentContactPhysics->normalForce = Vector3r::ZERO;//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;
+ }
+
+ 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
+
+ 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 = 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;
}
}
-
-
-
-
-
-
-
-
-
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.hpp 2010-04-19 13:23:53 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.hpp 2010-04-26 09:02:51 +0000
@@ -23,7 +23,7 @@
((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..."))
+ ((Real,creep_viscosity,1,"creep viscosity [Pa.s/m]. probably should be moved to Ip2_2xCohFrictMat_CohFrictPhys..."))
((bool,erosionActivated,false,""))
((bool,detectBrokenBodies,false,""))
=== modified file 'pkg/dem/Engine/PartialEngine/TriaxialCompressionEngine.hpp'
--- pkg/dem/Engine/PartialEngine/TriaxialCompressionEngine.hpp 2010-04-22 16:40:45 +0000
+++ pkg/dem/Engine/PartialEngine/TriaxialCompressionEngine.hpp 2010-04-26 09:02:51 +0000
@@ -120,6 +120,7 @@
firstRun=true;
previousSigmaIso=sigma_iso;
boxVolume=0;
+ saveSimulation=false;
,
.def("setContactProperties",&TriaxialCompressionEngine::setContactProperties,"Assign a new friction angle (degrees) to dynamic bodies and relative interactions")
)
