------------------------------------------------------------ revno: 2519 committer: Sergei D. <s...@think> branch nick: trunk timestamp: Wed 2010-10-27 22:54:33 +0400 message: 1. Fix ViscElBasic law to handle periodicity; 2. Add some regression tests added: scripts/test/facet-sphere-ViscElBasic-peri.py scripts/test/facet-sphere-ViscElBasic.py scripts/test/sphere-sphere-ViscElBasic-peri.py modified: pkg/dem/Ig2_Facet_Sphere_ScGeom.cpp pkg/dem/ViscoelasticPM.cpp
-- 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/Ig2_Facet_Sphere_ScGeom.cpp' --- pkg/dem/Ig2_Facet_Sphere_ScGeom.cpp 2010-10-25 20:31:09 +0000 +++ pkg/dem/Ig2_Facet_Sphere_ScGeom.cpp 2010-10-27 18:54:33 +0000 @@ -107,7 +107,7 @@ scm->radius1 = 2*sphereRadius; scm->radius2 = sphereRadius; if (isNew) c->geom = scm; - scm->precompute(state1,state2,scene,c,normal,isNew,shift2,true); + scm->precompute(state1,state2,scene,c,normal,isNew,shift2,false/*avoidGranularRatcheting only for sphere-sphere*/); return true; } return false; === modified file 'pkg/dem/ViscoelasticPM.cpp' --- pkg/dem/ViscoelasticPM.cpp 2010-10-25 20:31:09 +0000 +++ pkg/dem/ViscoelasticPM.cpp 2010-10-27 18:54:33 +0000 @@ -67,20 +67,17 @@ axis = angle*geom.normal; shearForce -= shearForce.cross(axis); - //const Vector3r c1x = (geom.contactPoint - de1.pos); - //const Vector3r c2x = (geom.contactPoint - de2.pos); - const Vector3r c1x = geom.radius1*geom.normal; - const Vector3r c2x = -geom.radius2*geom.normal; - /// The following definition of c1x and c2x is to avoid "granular ratcheting" - /// (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN, - /// "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena", - /// ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors) - //Vector3r _c1x_ = geom->radius1*geom->normal; - //Vector3r _c2x_ = -geom->radius2*geom->normal; - //Vector3r relativeVelocity = (de2->velocity+de2->angularVelocity.Cross(_c2x_)) - (de1->velocity+de1->angularVelocity.Cross(_c1x_)); - const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) ; + // Handle periodicity. + const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(I->cellDist): Vector3r::Zero(); + const Vector3r shiftVel = scene->isPeriodic ? scene->cell->intrShiftVel(I->cellDist): Vector3r::Zero(); + + const Vector3r c1x = (geom.contactPoint - de1.pos); + const Vector3r c2x = (geom.contactPoint - de2.pos - shift2); + + const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) + shiftVel; const Real normalVelocity = geom.normal.dot(relativeVelocity); const Vector3r shearVelocity = relativeVelocity-normalVelocity*geom.normal; + // As Chiara Modenese suggest, we store the elastic part // and then add the viscous part if we pass the Mohr-Coulomb criterion. // See http://www.mail-archive.com/[email protected]/msg01391.html === added file 'scripts/test/facet-sphere-ViscElBasic-peri.py' --- scripts/test/facet-sphere-ViscElBasic-peri.py 1970-01-01 00:00:00 +0000 +++ scripts/test/facet-sphere-ViscElBasic-peri.py 2010-10-27 18:54:33 +0000 @@ -0,0 +1,53 @@ +# -*- coding: utf-8 + +# Testing facet-sphere interaction in periodic case. +# Pass, if the sphere is rolling from left to right through the period. + +from yade import utils + +sphereRadius=0.1 +tc=0.001# collision time +en=0.3 # normal restitution coefficient +es=0.3 # tangential restitution coefficient +density=2700 +frictionAngle=radians(35)# +params=utils.getViscoelasticFromSpheresInteraction(tc,en,es) +facetMat=O.materials.append(ViscElMat(frictionAngle=frictionAngle,**params)) +sphereMat=O.materials.append(ViscElMat(density=density,frictionAngle=frictionAngle,**params)) + +#floor +n=5. +s=1./n +for i in range(0,n): + for j in range(0,n): + O.bodies.append([ + utils.facet( [(i*s,j*s,0.1),(i*s,(j+1)*s,0.1),((i+1)*s,(j+1)*s,0.1)],material=facetMat), + utils.facet( [(i*s,j*s,0.1),((i+1)*s,j*s,0.1),((i+1)*s,(j+1)*s,0.1)],material=facetMat), + ]) + +# Spheres +sphId=O.bodies.append([ + utils.sphere( (0.5,0.5,0.2), 0.1, material=sphereMat), + ]) +O.bodies[sphId[-1]].state.vel=(0.5,0,0) + +## Engines +O.engines=[ + ForceResetter(), + InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]), + InteractionLoop( + [Ig2_Facet_Sphere_ScGeom()], + [Ip2_ViscElMat_ViscElMat_ViscElPhys()], + [Law2_ScGeom_ViscElPhys_Basic()], + ), + GravityEngine(gravity=[0,0,-9.81]), + NewtonIntegrator(damping=0), +] + +O.periodic=True +O.cell.refSize=Vector3(1,1,1) + +O.dt=.01*tc + +O.saveTmp() + === added file 'scripts/test/facet-sphere-ViscElBasic.py' --- scripts/test/facet-sphere-ViscElBasic.py 1970-01-01 00:00:00 +0000 +++ scripts/test/facet-sphere-ViscElBasic.py 2010-10-27 18:54:33 +0000 @@ -0,0 +1,48 @@ +# -*- coding: utf-8 -*- + +# Testing facet-sphere interaction. +# A facet is rotated around Z axis. Test pass, if a sphere at (0,0) position is not moving (because in this case no transfer moment from the facet to the sphere), but a sphere at facet's edge moves with the facet (for this sphere blocked the rotation DOFs in order to remove rolling). + +## PhysicalParameters +Density=2400 +frictionAngle=radians(35) +tc = 0.001 +en = 0.3 +es = 0.3 + +## Import wall's geometry +params=utils.getViscoelasticFromSpheresInteraction(tc,en,es) +facetMat=O.materials.append(ViscElMat(frictionAngle=frictionAngle,**params)) +sphereMat=O.materials.append(ViscElMat(density=Density,frictionAngle=frictionAngle,**params)) + +facetId=O.bodies.append(utils.facet( [ (-1,0,0), (1,1,0), (1,-1,0)], material=facetMat,color=(1,0,0))) + +sphIds=O.bodies.append([ + utils.sphere( (0,0,0.1),0.1, material=sphereMat,color=(0,1,0)), + utils.sphere( (0.9,0,0.1),0.1, material=sphereMat,color=(0,1,0)) + ]) + +O.bodies[sphIds[1]].state.blockedDOFs=['rx','ry','rz'] + +## Timestep +O.dt=.1*tc + +## Engines +O.engines=[ + ForceResetter(), + InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]), + InteractionLoop( + [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()], + [Ip2_ViscElMat_ViscElMat_ViscElPhys()], + [Law2_ScGeom_ViscElPhys_Basic()], + ), + GravityEngine(gravity=[0,0,-9.81]), + NewtonIntegrator(damping=0), + RotationEngine(subscribedBodies=[facetId],rotationAxis=[0,0,1],rotateAroundZero=True,angularVelocity=0.1) +] + +from yade import qt +qt.View() +O.saveTmp() +#O.run() + === added file 'scripts/test/sphere-sphere-ViscElBasic-peri.py' --- scripts/test/sphere-sphere-ViscElBasic-peri.py 1970-01-01 00:00:00 +0000 +++ scripts/test/sphere-sphere-ViscElBasic-peri.py 2010-10-27 18:54:33 +0000 @@ -0,0 +1,44 @@ +# -*- coding: utf-8 + +# Testing sphere-sphere interaction in periodic case. +# Pass, if the spheres moves along the X axis, interacting through the period. + +from yade import utils + +sphereRadius=0.1 +tc=0.001# collision time +en=1 # normal restitution coefficient +es=1 # tangential restitution coefficient +density=2700 +frictionAngle=radians(35)# +params=utils.getViscoelasticFromSpheresInteraction(tc,en,es) +sphereMat=O.materials.append(ViscElMat(density=density,frictionAngle=frictionAngle,**params)) + + +# Spheres +sphId=O.bodies.append([ + utils.sphere( (0.4,0.5,0.5), 0.1, material=sphereMat), + utils.sphere( (0.6,0.5,0.5), 0.1, material=sphereMat) + ]) +O.bodies[sphId[-1]].state.vel=(0.5,0,0) +O.bodies[sphId[0]].state.vel=(-0.5,0,0) + +## Engines +O.engines=[ + ForceResetter(), + InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]), + InteractionLoop( + [Ig2_Sphere_Sphere_ScGeom()], + [Ip2_ViscElMat_ViscElMat_ViscElPhys()], + [Law2_ScGeom_ViscElPhys_Basic()], + ), + NewtonIntegrator(damping=0), +] + +O.periodic=True +O.cell.refSize=Vector3(1,1,1) + +O.dt=.01*tc + +O.saveTmp() +
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
