------------------------------------------------------------ revno: 2170 committer: Bruno Chareyre <bchare...@r1arduina> branch nick: trunk timestamp: Wed 2010-04-21 19:17:54 +0200 message: - Make the radius of fictious sphere equal to the one of real sphere in box-sphere geometry. - Implement energy tracing in ElasticContactLaw - Remove some "#ifdef SCG_SHEAR" modified: pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp pkg/dem/Engine/GlobalEngine/ElasticContactLaw.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/InteractionGeometry/ScGeom.cpp' --- pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2009-12-13 20:11:31 +0000 +++ pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2010-04-21 17:17:54 +0000 @@ -65,7 +65,7 @@ } #endif -void ScGeom::updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){ +Vector3r ScGeom::updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){ Vector3r axis; Real angle; @@ -113,6 +113,7 @@ Vector3r shearVelocity = relativeVelocity-normal.Dot(relativeVelocity)*normal; Vector3r shearDisplacement = shearVelocity*dt; shearForce -= ks*shearDisplacement; + return shearDisplacement; } === modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp' --- pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp 2010-04-20 14:57:37 +0000 +++ pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp 2010-04-21 17:17:54 +0000 @@ -27,7 +27,7 @@ virtual ~ScGeom(); - void updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true); + Vector3r updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true); YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom,GenericSpheresContact,"Class representing :yref:`geometry<InteractionGeometry>` of two :yref:`spheres<Sphere>` in contact. The contact has 3 DOFs (normal and 2Ãshear) and uses incremental algorithm for updating shear. (For shear formulated in total displacements and rotations, see :yref:`Dem3DofGeom` and related classes).\n\nWe use symbols $\\vec{x}$, $\\vec{v}$, $\\vec{\\omega}$ respectively for position, linear and angular velocities (all in global coordinates) and $r$ for particles radii; subscripted with 1 or 2 to distinguish 2 spheres in contact. Then we compute unit contact normal\n\n.. math::\n\n\t\\vec{n}=\\frac{\\vec{x}_2-\\vec{x}_1}{||\\vec{x}_2-\\vec{x}_1||}\n\nRelative velocity of spheres is then\n\n.. math::\n\n\t\\vec{v}_{12}=(\\vec{v}_2+\\vec{\\omega}_2\\times(-r_2\\vec{n}))-(\\vec{v}_1+\\vec{\\omega}_1\\times(r_1\\vec{n}))\n\nand its shear component\n\n.. math::\n\n\t\\Delta\\vec{v}_{12}^s=\\vec{v}_{12}-(\\vec{n}\\cdot\\vec{v}_{12})\\vec{n}.\n\nTangential displacement increment over last step then reads\n\n.. math::\n\n\t\\vec{x}_{12}^s=\\Delta t \\vec{v}_{12}^s.", ((Vector3r,contactPoint,Vector3r::ZERO,"Reference point of the contact. |ycomp|")) === modified file 'pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp' --- pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2010-02-09 20:22:04 +0000 +++ pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2010-04-21 17:17:54 +0000 @@ -93,17 +93,15 @@ if (isNew) scm = shared_ptr<ScGeom>(new ScGeom()); else scm = YADE_PTR_CAST<ScGeom>(c->interactionGeometry); - #ifdef SCG_SHEAR - if(isNew) { /* same as below */ scm->prevNormal=pt1-pt2; scm->prevNormal.Normalize(); } - else {scm->prevNormal=scm->normal;} - #endif + if(isNew) { /* same as below */ scm->prevNormal=pt1-pt2; scm->prevNormal.Normalize(); } + else {scm->prevNormal=scm->normal;} // contact point is in the middle of overlapping volumes //(in the direction of penetration, which is normal to the box surface closest to sphere center) of overlapping volumes scm->contactPoint = 0.5*(pt1+pt2); scm->normal = pt1-pt2; scm->normal.Normalize(); scm->penetrationDepth = (pt1-pt2).Length(); - scm->radius1 = s->radius*2; + scm->radius1 = s->radius; scm->radius2 = s->radius; c->interactionGeometry = scm; } else { // outside @@ -140,10 +138,9 @@ bool isNew=!c->interactionGeometry; if (isNew) scm = shared_ptr<ScGeom>(new ScGeom()); else scm = YADE_PTR_CAST<ScGeom>(c->interactionGeometry); - #ifdef SCG_SHEAR - if(isNew) { /* same as below */ scm->prevNormal=-cOnBox_sphere; } - else {scm->prevNormal=scm->normal;} - #endif + if(isNew) { /* same as below */ scm->prevNormal=-cOnBox_sphere; } + else {scm->prevNormal=scm->normal;} + scm->contactPoint = 0.5*(pt1+pt2); //scm->normal = pt1-pt2; scm->normal.Normalize(); //scm->penetrationDepth = (pt1-pt2).Length(); @@ -151,7 +148,7 @@ scm->penetrationDepth = depth; - scm->radius1 = s->radius*2; + scm->radius1 = s->radius; scm->radius2 = s->radius; c->interactionGeometry = scm; } === modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp' --- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-04-19 13:33:45 +0000 +++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-04-21 17:17:54 +0000 @@ -33,13 +33,23 @@ } } +Real Law2_ScGeom_FrictPhys_Basic::elasticEnergy() +{ + Real energy=0; + FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){ + if(!I->isReal()) continue; + FrictPhys* phys = dynamic_cast<FrictPhys*>(I->interactionPhysics.get()); + if(phys) { + energy += 0.5*(phys->normalForce.SquaredLength()/phys->kn + phys->shearForce.SquaredLength()/phys->ks);} + } + return energy; +} CREATE_LOGGER(Law2_ScGeom_FrictPhys_Basic); void Law2_ScGeom_FrictPhys_Basic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* ncb){ Real dt = Omega::instance().getTimeStep(); int id1 = contact->getId1(), id2 = contact->getId2(); -// // FIXME: mask handling should move to LawFunctor itself, outside the functors -// if( !(Body::byId(id1,ncb)->getGroupMask() & Body::byId(id2,ncb)->getGroupMask() & sdecGroupMask) ) return; + ScGeom* currentContactGeometry= static_cast<ScGeom*>(ig.get()); FrictPhys* currentContactPhysics = static_cast<FrictPhys*>(ip.get()); if(currentContactGeometry->penetrationDepth <0){ @@ -54,18 +64,42 @@ Real un=currentContactGeometry->penetrationDepth; TRVAR3(currentContactGeometry->penetrationDepth,de1->se3.position,de2->se3.position); currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal; - if(useShear){ - currentContactGeometry->updateShear(de1,de2,dt); - shearForce=currentContactPhysics->ks*currentContactGeometry->shear; - } else { - currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);} - // PFC3d SlipModel, is using friction angle. CoulombCriterion - Real maxFs = currentContactPhysics->normalForce.SquaredLength()* - std::pow(currentContactPhysics->tangensOfFrictionAngle,2); - if( shearForce.SquaredLength() > maxFs ){ - Real ratio = Mathr::Sqrt(maxFs) / shearForce.Length(); - shearForce *= ratio; - if(useShear) currentContactGeometry->shear*=ratio;} + + if (!traceEnergy){//Update force but don't compute energy terms + if(useShear){ + currentContactGeometry->updateShear(de1,de2,dt); + shearForce=currentContactPhysics->ks*currentContactGeometry->shear; + } else { + currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);} + // PFC3d SlipModel, is using friction angle. CoulombCriterion + Real maxFs = currentContactPhysics->normalForce.SquaredLength()* + std::pow(currentContactPhysics->tangensOfFrictionAngle,2); + if( shearForce.SquaredLength() > maxFs ){ + Real ratio = Mathr::Sqrt(maxFs) / shearForce.Length(); + shearForce *= ratio; + if(useShear) currentContactGeometry->shear*=ratio;} + } else {//almost the same with 2 additional Vector3r instanciated for energy tracing, duplicated block to make sure there is no cost for the instanciation of the vectors when traceEnergy==false + Vector3r prevForce=shearForce;//store prev force for definition of plastic slip + if(useShear) throw ("energy tracing not defined with useShear==true"); + /*{ + currentContactGeometry->updateShear(de1,de2,dt); + shearForce=currentContactPhysics->ks*currentContactGeometry->shear;//FIXME : energy terms if useShear? + } else {*/ + Vector3r shearDisp = currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt); + //} + // PFC3d SlipModel, is using friction angle. CoulombCriterion + Real maxFs = currentContactPhysics->normalForce.SquaredLength()* + std::pow(currentContactPhysics->tangensOfFrictionAngle,2); + if( shearForce.SquaredLength() > maxFs ){ + Real ratio = Mathr::Sqrt(maxFs) / shearForce.Length(); + //define the plastic work input and increment the total plastic energy dissipated + plasticDissipation += + (shearDisp+(1/currentContactPhysics->ks)*(shearForce-prevForce))//plastic disp. + .Dot(shearForce);//active force + shearForce *= ratio; + //if(useShear) currentContactGeometry->shear*=ratio; + } + } applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position, ncb); currentContactPhysics->prevNormal = currentContactGeometry->normal; } === modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp' --- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp 2010-04-19 13:33:45 +0000 +++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp 2010-04-21 17:17:54 +0000 @@ -21,10 +21,15 @@ class Law2_ScGeom_FrictPhys_Basic: public LawFunctor{ public: virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene*); - YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_FrictPhys_Basic,LawFunctor,"Law for linear compression, without cohesion and Mohr-Coulomb plasticity surface.\n\n.. note::\n This law uses :yref:`ScGeom`; there is also functionally equivalent :yref:`Law2_Dem3DofGeom_FrictPhys_Basic`, which uses :yref:`Dem3DofGeom`.", + Real elasticEnergy (); + YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_FrictPhys_Basic,LawFunctor,"Law for linear compression, without cohesion and Mohr-Coulomb plasticity surface.\n\n.. note::\n This law uses :yref:`ScGeom`; there is also functionally equivalent :yref:`Law2_Dem3DofGeom_FrictPhys_Basic`, which uses :yref:`Dem3DofGeom`.", ((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,useShear,false,"Use ScGeom::updateShear rather than ScGeom::updateShearForce for shear force computation.")) + ((bool,traceEnergy,false,"Define the total energy dissipated in plastic slips at all contacts.")) + ((Real,plasticDissipation,0,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:`Law2_ScGeom_FrictPhys_Basic::traceEnergy` is true. |yupdate|")) + ,, + .def("elasticEnergy",&Law2_ScGeom_FrictPhys_Basic::elasticEnergy,"Compute and return the total elastic energy in all \"FrictPhys\" contacts") ); FUNCTOR2D(ScGeom,FrictPhys); DECLARE_LOGGER;
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp