yade-users team mailing list archive

[Jan Stránský <question295521@xxxxxxxxxxxxxxxxxxxxx>]

Question #295521 on Yade changed:
https://answers.launchpad.net/yade/+question/295521

    Status: Open => Answered

Jan Stránský proposed the following answer:
Hello,

there really must be something wrong in your calculations.

first of all, as I understand it, plotting normal is meaningful only for
values (0,pi), the other half may be mirrored if needed. From contact
direction point of view, vectors (1,0,0) and (-1,0,0) are identical..

with your script, I got perfect results.. either with
quadrant_X_num, quadrant_1_num+quadrant_3_num is almost exactly equal to
quadrant_2_num+quadrant_4_num (see above why there is the sum :-)

Also putting the directions into more "bins", I got e.g.
##############################################
def contactDirs(nBins=12):
normals = []
for i in O.interactions: # use only sphere-sphere contacts
if not isinstance(O.bodies[i.id1].shape,Sphere):
continue
if not isinstance(O.bodies[i.id2].shape,Sphere):
continue
normals.append(i.geom.normal)
# until now could be done in one line, but would look ugly in email
bins = [0 for i in range(nBins)]
for n in normals:
angle = atan2(n[1],n[0]) # atan2 returns value in range (-pi,pi)
if angle<0:
angle += pi
i = int(angle/pi * nBins)
bins[i] += 1
return bins
##############################################

Yade [1]: contactDirs(6)
 ->  [1]: [475, 483, 484, 533, 486, 518]

Yade [2]: contactDirs(12)
 ->  [2]: [268, 207, 264, 219, 249, 235, 262, 271, 255, 231, 235, 283]


which seems to be reasonable

cheers
Jan


2016-06-22 10:13 GMT+02:00 Fu zuoguang <question295521@xxxxxxxxxxxxxxxxxxxxx
>:

> Question #295521 on Yade changed:
> https://answers.launchpad.net/yade/+question/295521
>
> Fu zuoguang gave more information on the question:
> # encoding: utf-8
> # the script demonstrates a simple case of triaxial simulation using
> TriaxialCompressionEngine. More elaborated examples can be found in the
> triax-tutorial folder
> from yade import pack
>
> sp=pack.SpherePack()
> ## corners of the initial packing
> mn,mx=Vector3(0,0,0.5),Vector3(1,1,0.5)
>
> ## box between mn and mx, avg radius ± ½(20%), 2k spheres
> sp.makeCloud(minCorner=mn,maxCorner=mx,rRelFuzz=.2,num=2000)
>
> ## create material #0, which will be used as default
>
> O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=radians(30),density=2600,label='spheres'))
>
> O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=0,density=0,label='frictionless'))
>
> ## copy spheres from the packing into the scene
> ## use default material, don't care about that for now
> O.bodies.append([sphere(center,rad,material='spheres') for center,rad in
> sp])
>
> ## create walls around the packing
> walls=aabbWalls(thickness=1e-10,material='frictionless')
> wallIds=O.bodies.append(walls)
>
> # block the dof in z-direction.
> for b in O.bodies:
>         if isinstance(b.shape, Sphere):b.state.blockedDOFs = 'zXY'
>
> # define the engine.
> triax=TriaxialStressController(
>
> wall_bottom_id=wallIds[2],wall_top_id=wallIds[3],wall_left_id=wallIds[0],
>
> wall_right_id=wallIds[1],wall_back_id=wallIds[4],wall_front_id=wallIds[5],
>
>     # keep all the walls from moving except the back and front walls.
>     wall_back_activated = 0, wall_front_activated = 0,
>     thickness = 1e-9, internalCompaction = 0,
>
>     # stress controlling condition.
>     stressMask = 7,
>         goal1 = -100000,goal2 = -100000,goal3 = -0.0,
>
>         label="triax", )
>
> O.engines=[
>         ForceResetter(),
>         InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
>         InteractionLoop(
>                 [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
>                 [Ip2_FrictMat_FrictMat_FrictPhys()],
>                 [Law2_ScGeom_FrictPhys_CundallStrack()]
>         ),
>
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
>         triax,
>         # you can add TriaxialStateRecorder and such here…
>         NewtonIntegrator(damping=.4)
> ]
>
> # run this simulation.
> O.run(200000, True)
> ##################################  running simulation process is over!
>
> # record all the contact normals in the final state.
> contact_normal_list=[[i.geom.normal[0],i.geom.normal[1]] for i in
> O.interactions]
>
> # divide the whole Cartesian system into 4 quadrants as: (1)x>0 and y>0;
> (2)x<0 and y>0 ...
> # then summarize respectively the number of contact normals located in
> each area with
> # ignoring cursoryly the sphere-wall contacts.
>
> quadrant_1_num, quadrant_2_num, quadrant_3_num, quadrant_4_num, =0, 0,
> 0, 0
>
> for i in xrange(len(contact_normal_list)):
>     if (contact_normal_list[i][0]>0 and contact_normal_list[i][1]>0):
> quadrant_1_num=quadrant_1_num+1
>     if (contact_normal_list[i][0]<0 and contact_normal_list[i][1]>0):
> quadrant_2_num=quadrant_2_num+1
>     if (contact_normal_list[i][0]<0 and contact_normal_list[i][1]<0):
> quadrant_3_num=quadrant_3_num+1
>     if (contact_normal_list[i][0]>0 and contact_normal_list[i][1]<0):
> quadrant_4_num=quadrant_4_num+1
>     else:pass
>
> print 'the total number of contacts is: ', len(contact_normal_list)
> print 'the number of sphere-sphere contacts is:
> ',(quadrant_1_num+quadrant_2_num+quadrant_3_num+quadrant_4_num)
> print 'the number of contacts in 1 quadrant is: ', quadrant_1_num
> print 'the number of contacts in 2 quadrant is: ', quadrant_2_num
> print 'the number of contacts in 3 quadrant is: ', quadrant_3_num
> print 'the number of contacts in 4 quadrant is: ', quadrant_4_num
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> 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
>

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.