------------------------------------------------------------ revno: 2182 committer: Bruno Chareyre <bchare...@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 <[email protected]> +/************************************************************************* +* Copyright (C) 2007 by Bruno CHAREYRE * +* [email protected] * +* * +* 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") )
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
