yade-dev team mailing list archive

[Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>]

(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/yade-users@xxxxxxxxxxxxxxxxxxx/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     : yade-users@xxxxxxxxxxxxxxxxxxx
> 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