(moving to yade-dev) Hi Vaclav,
I like reading your scripts because I learn something at each line (this video-generating code is awesome!). The attached script is comparing L6 and ScGeom6 in the L-shaped beam example. The results are : free end position in L6GeomInc : Vector3(3.3886340639637584,7.8792525473155308,-6.7166071545283756) free end position in ScGeom6D : Vector3(3.3823175232092986,7.8859989218194588,-6.7148657033006511) Now, more annoying : type this once the beam is stabilized : O.bodies[0].state.angVel=Vector3(0,0,0.01) In both cases you will see growing errors in the form of a weird deformed shape after enough 360° turns (earlier with ScGeom, I confess). Not sure if this is due to approximations in shearDisp, rotations, both, or still something else. It would be interesting to investigate this, though I doubt it really invalidate simulations of dense packings. Side questions on this. 1. It seems that L3Geom::applyLocalForceTorque will always use radius and normal to define x1c. Isn't it therefore restricting to spherical bodies, as in Law2_ScGeom_FrictPhys_CundallStrack::sphericalBodies=true? (see http://www.mail-archive.com/[email protected]/msg02715.html) 2. It is not clear for me how local coordinates make things simpler. For linear elasticity, it looks exactly the same (3 lines in both cases). 3. I thought that L3/L6 was for research purpose. Don't you think that, if users start implementing more complete laws (not just linear elastic) on this basis, it will result in duplicating existing code (Ig2_Sphere_Sphere_L3Geom_Inc itself is in a large extent a copy of ScGeom code - not saying you pasted, but it is really the same)? Bruno On 07/12/10 14:53, Václav Šmilauer wrote: > Hi there, > > I commited functional version of L6Geom recently > (http://www.youtube.com/watch?v=YYGHjiZRUpM, shown in > scripts/test/beam-l6geom.py), contact geometry with 6 degrees of freedom > with local coordinates (L). You (or Giulia) could consider using it as > the base for writing a new contact law; it is ridiculously easy due to > local coordinates > (http://bazaar.launchpad.net/~yade-dev/yade/trunk/annotate/head%3A/pkg/dem/L3Geom.cpp#L256 > is the one used in the movie). > > Cheers, Vaclav >> quick question. A colleague of mine needs to use the linear contact >> law with moment transfer rotation. Which law in Yade is the most >> updated and satisfy these requirements? Thanks, so we avoid >> duplicates. If there is something to add/modify let me know, I can do it. > > > _______________________________________________ > Mailing list: https://launchpad.net/~yade-users > Post to : [email protected] > Unsubscribe : https://launchpad.net/~yade-users > More help : https://help.launchpad.net/ListHelp > -- _______________ Bruno Chareyre Associate Professor ENSE³ - Grenoble INP Lab. 3SR BP 53 - 38041, Grenoble cedex 9 - France Tél : +33 4 56 52 86 21 Fax : +33 4 76 82 70 43 ________________
# # L6Geom / ScGeom6D comparison # import numpy # radius, number and distance of spheres rad,num=1,6; dist=1.9999*rad young=1.0e7 poisson=1 density=2.60e3 frictionAngle=radians(30) O.materials.append(CohFrictMat(alphaKr=1.0, alphaKtw=1.0,young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,normalCohesion=1e13,shearCohesion=1e13,momentRotationLaw=True,label='mat')) # one arm O.bodies.append([utils.sphere((0,y,0),rad,wire=True) for y in numpy.arange(0,2*num-1,dist)]) # the lateral arm O.bodies.append([utils.sphere((x,(num-1)*2*rad,0),rad,wire=True) for x in numpy.arange(dist,1+num/2,dist)]) # support sphere O.bodies[0].state.blockedDOFs=['x','y','z','rx','ry','rz'] # small dt to see in realtime how it swings; real critical is higher, but much less than p-wave O.dt=.1*utils.PWaveTimeStep() O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb()]), InteractionLoop([Ig2_Sphere_Sphere_L6Geom_Inc()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_L6Geom_FrictPhys_Linear(charLen=1)]), GravityEngine(gravity=(0,0,-9.81)), NewtonIntegrator(damping=0.4), ] O.saveTmp() O.run(50000,True) print "free end position in L6GeomInc : ",O.bodies[7].state.pos O.reload() O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Sphere_ChainedCylinder_CylScGeom()], [Ip2_2xCohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True)], [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(label='law')] ), GravityEngine(gravity=(0,0,-9.81)), NewtonIntegrator(damping=0.4), ] O.run(50000,True) print "free end position in ScGeom6D : ",O.bodies[7].state.pos
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
