problems with facet and calculation

Asked by Christian Jakob on 2011-07-20

Hi,

I created a model with 1000 spheres and one facet. There I cannot fix two problems:

1.) The facet is in xy-plane, but I want to have it normal to (.02,.02,1), not (0,0,1)!
How can I do this?
I tried different things with facet, facetBox, wall and facetPolygon, but I could not fix it.

2.) Spheres are not moving and I dont know why. timing.stats() gives:

Name Count Time Rel. time
-------------------------------------------------------------------------------------------------------
ForceResetter 0 0us
InsertionSortCollider 0 0us
InteractionLoop 0 0us
GravityEngine 0 0us
NewtonIntegrator 0 0us
TOTAL 0us

Now the script:

# script START ------------------------------------------------------------------------------------

#!/usr/bin/python
# -*- coding: utf-8 -*-

# speed test file for YADE

num_balls1D = 10

num_balls = num_balls1D*num_balls1D*num_balls1D

logfilename = 'speed-test-output%i-balls.log' % (int(num_balls))

origin_wall = num_balls1D/2

x_pos = 0

y_pos = 0

z_pos = 0
friction=0.5
angle=atan(friction)

id_FacetMat=O.materials.append(ViscElMat(kn=1e8,ks=1e8,cn=0.0,cs=0.0,frictionAngle=angle))
id_SphereMat=O.materials.append(ViscElMat(kn=1e6,ks=1e6,cn=0.0,cs=0.0,density=1000,frictionAngle=angle))

FacetMat=O.materials[id_FacetMat]
SphereMat=O.materials[id_SphereMat]

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.7),
]

id_facet=O.bodies.append(geom.facetBox((origin_wall,origin_wall,0),(100,100,0),material=FacetMat))

for ii in range(1,num_balls1D):

 x_pos = x_pos + 1

 y_pos = 0

 z_pos = 0

 for jj in range(1,num_balls1D):

  y_pos = y_pos + 1

  z_pos = 0

  for kk in range(1,num_balls1D):

   z_pos = z_pos + 1
   O.bodies.append(utils.sphere([x_pos,y_pos,z_pos], material=SphereMat, radius=0.5, wire=False, highlight=False))#, color=(50,50,0)

from yade import timing

O.dt=1e-3

start_time=O.realtime
print start_time

O.run(100,True)

stop_time=O.realtime
time_needed=stop_time-start_time
print time_needed

O.timingEnabled==True
timing.stats()

# script END ------------------------------------------------------------------------------------

Christian.

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Christian Jakob
Solved:
2011-07-22
Last query:
2011-07-22
Last reply:
2011-07-20

1/ I'm not used to facets creation, you'll need help from somebody else.
2/ Spheres are moving in your script, why do you say they don't? I see
them going down until they touch the facets (and guess what happens when
a vertical column of spheres is under gravity: they are at equilibrium.
Try e.g. gravity (1,-9.81,0) and you will see).
3/ You need to enable timing before calculation cycles if you want to
measure something. Or, at the end of your script, run a few additional
steps and type timing.stats() again.
Also, you have a typo: O.timingEnabled=True, instead of
O.timingEnabled==True.
4/ Looking at your damping value (0.7), I bet you are a PFC user? If you
are doing speed comparisons, would you mind posting some results here?
I'm all interested.

Bruno

Name
Count Time Rel. time
-------------------------------------------------------------------------------------------------------
ForceResetter 100
318us 0.32%
InsertionSortCollider 3
9891us 10.05%
InteractionLoop 100
73487us 74.67%
GravityEngine 100
1891us 1.92%
NewtonIntegrator 100
12821us 13.03%
TOTAL
98409us 100.00%

On 20/07/11 15:36, Christian Jakob wrote:
> New question #165458 on Yade:
> https://answers.launchpad.net/yade/+question/165458
>
> Hi,
>
> I created a model with 1000 spheres and one facet. There I cannot fix two problems:
>
> 1.) The facet is in xy-plane, but I want to have it normal to (.02,.02,1), not (0,0,1)!
> How can I do this?
> I tried different things with facet, facetBox, wall and facetPolygon, but I could not fix it.
>
> 2.) Spheres are not moving and I dont know why. timing.stats() gives:
>
> Name Count Time Rel. time
> -------------------------------------------------------------------------------------------------------
> ForceResetter 0 0us
> InsertionSortCollider 0 0us
> InteractionLoop 0 0us
> GravityEngine 0 0us
> NewtonIntegrator 0 0us
> TOTAL 0us
>
>
>
> Now the script:
>
> # script START ------------------------------------------------------------------------------------
>
> #!/usr/bin/python
> # -*- coding: utf-8 -*-
>
> # speed test file for YADE
>
>
> num_balls1D = 10
>
> num_balls = num_balls1D*num_balls1D*num_balls1D
>
> logfilename = 'speed-test-output%i-balls.log' % (int(num_balls))
>
> origin_wall = num_balls1D/2
>
> x_pos = 0
>
> y_pos = 0
>
> z_pos = 0
> friction=0.5
> angle=atan(friction)
>
> id_FacetMat=O.materials.append(ViscElMat(kn=1e8,ks=1e8,cn=0.0,cs=0.0,frictionAngle=angle))
> id_SphereMat=O.materials.append(ViscElMat(kn=1e6,ks=1e6,cn=0.0,cs=0.0,density=1000,frictionAngle=angle))
>
> FacetMat=O.materials[id_FacetMat]
> SphereMat=O.materials[id_SphereMat]
>
> 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.7),
> ]
>
> id_facet=O.bodies.append(geom.facetBox((origin_wall,origin_wall,0),(100,100,0),material=FacetMat))
>
> for ii in range(1,num_balls1D):
>
> x_pos = x_pos + 1
>
> y_pos = 0
>
> z_pos = 0
>
> for jj in range(1,num_balls1D):
>
> y_pos = y_pos + 1
>
> z_pos = 0
>
> for kk in range(1,num_balls1D):
>
> z_pos = z_pos + 1
> O.bodies.append(utils.sphere([x_pos,y_pos,z_pos], material=SphereMat, radius=0.5, wire=False, highlight=False))#, color=(50,50,0)
>
>
> from yade import timing
>
> O.dt=1e-3
>
> start_time=O.realtime
> print start_time
>
> O.run(100,True)
>
> stop_time=O.realtime
> time_needed=stop_time-start_time
> print time_needed
>
> O.timingEnabled==True
> timing.stats()
>
> # script END ------------------------------------------------------------------------------------
>
> Christian.
>

--
_______________
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
________________

Anton Gladky (gladky-anton) said : #2

Hi,

> 1.) The facet is in xy-plane, but I want to have it normal to (.02,.02,1), not (0,0,1)!
> How can I do this?
> I tried different things with facet, facetBox, wall and facetPolygon, but I could not fix it.

You want to rotate your plane initially? If so, you should use
"orientation" parameter in facetBox (see examples/pack/pack.py)

2.) Spheres are not moving and I dont know why. timing.stats() gives:

Add O.wait() directly after O.run()

HTH,

Anton

> 2.) Spheres are not moving and I dont know why. timing.stats() gives:
>
They are moving. See my previous post. :)
No need to O.wait() since there is already "True" in O.run(...,True), it
has the same effect.

Bruno

Christian Jakob (jakob-ifgt) said : #4

Thank you for your comments.

Using O.timingEnabled=True instead of O.timingEnabled==True fixed the problem with timing.stats().
If I had a deeper look at the packing, I see them moving now, too ;)

I tried to get the facet rotated, but I never used quaternions before. If I want to have my facet normal to (.02,.02,1), is this the right implementation then?:

#get the angle in radians:
tmp=acos(sqrt(2*.02*.02))
#rotation quaternion:
orientationFacet = Quaternion(Vector3(0.02,0.02,1),tmp)
#create facet:
id_facet=O.bodies.append(geom.facetBox((origin_wall,origin_wall,0),(100,100,0),orientationFacet,material=FacetMat))

@Bruno:

Yes I will try to compare calculation speed of PFC, YADE and another program called PASIMODO.
If the test runs are ready, I will send it all developers, if they want.

Christian

Christian Jakob (jakob-ifgt) said : #5

Ok I found out by my own how to rotate my plane facet to get a facet normal to (0.02,0.02,1):
(with a little help from google)

#rotation quaternion:
orientationFacet = Quaternion(Vector3(.01,.01,1),math.pi)
#create facet:
id_facet=O.bodies.append(geom.facetBox((origin_wall,origin_wall,0),(100,100,0),orientationFacet,material=FacetMat))

Isnt it possible to program something like a facetPlane, where the user can give either 4 vertices or a normal and a origin point?!

> Isnt it possible to program something like a facetPlane, where the user
> can give either 4 vertices or a normal and a origin point?!

I guess it's a few line in python to go from point+normal to quaternion.
Feel free to write such function, then we could add it to yade.utils.

Anton Gladky (gladky-anton) said : #7

On Fri, Jul 22, 2011 at 9:46 AM, Christian Jakob
<email address hidden> wrote:
>
> Isnt it possible to program something like a facetPlane, where the user
> can give either 4 vertices or a normal and a origin point?!
>

Sumply nobody needed that before.
Bruno is right, there should be just a couple of lines in python.

Anton