get porosity of a packing

Asked by ceguo on 2012-03-07

Hi,

The problem seems trivial but when I try using `utils.porosity` I can only get something like <Boost.Python.function object at ...>. The documentation says it should return a float.

Besides, as I am using the trunk version with a 2D packing. All the quantities like mass, porosity, stress (averaged by area rather than volume) etc. should be modified accordingly. Will the code take care of these right now? If not, maybe we could do scaling (3D to 2D) outside the code. So this is not a big problem.

Another problem is that the documentation says we'd better use ScGeom in stead of Dem3DofGeom. But as I try using `Ig2_Sphere_Sphere_ScGeom()`, `Ip2_FrictMat_FrictMat_FrictPhys()` and `Law2_ScGeom_FrictPhys_CundallStrack()`, the program runs without error but there is no stress generated.

Ning

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
ceguo
Solved:
2012-03-07
Last query:
2012-03-07
Last reply:
2012-03-07
Luc Scholtès (luc) said : #1

Hi Ceguo,

Try utils.porosity(X), with X the volume of your assembly and it should
work.

2012/3/7 ceguo <email address hidden>

> New question #189913 on Yade:
> https://answers.launchpad.net/yade/+question/189913
>
> Hi,
>
> The problem seems trivial but when I try using `utils.porosity` I can only
> get something like <Boost.Python.function object at ...>. The documentation
> says it should return a float.
>
> Besides, as I am using the trunk version with a 2D packing. All the
> quantities like mass, porosity, stress (averaged by area rather than
> volume) etc. should be modified accordingly. Will the code take care of
> these right now? If not, maybe we could do scaling (3D to 2D) outside the
> code. So this is not a big problem.
>
> Another problem is that the documentation says we'd better use ScGeom in
> stead of Dem3DofGeom. But as I try using `Ig2_Sphere_Sphere_ScGeom()`,
> `Ip2_FrictMat_FrictMat_FrictPhys()` and
> `Law2_ScGeom_FrictPhys_CundallStrack()`, the program runs without error but
> there is no stress generated.
>
> Ning
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

If you are in periodic bondaries, X is not needed. But since porosity is a function, it needs at least the bracket: "utils.porosity()"
Also note that "utils." is useless, "porosity()" will do.

>the program runs without error but there is no stress generated.

We can't guess anything without an example script Ning.

All quantities are 3D, you need to convert to 2D.

ceguo (hhh-guo) said : #3

Thanks,

I notice the missing of brackets for a function. Regarding to the last problem, I use the example script just with modification from Dem3DofGeom to ScGeom (Dem3DofGeom works, but as I read the doc which suggest us using ScGeom instead.)
#============================================================
"Test and demonstrate use of PeriTriaxController."
from yade import *
from yade import pack,log,qt
log.setLevel('PeriTriaxController',log.TRACE)
O.periodic=True
O.cell.refSize=Vector3(.1,.1,.1)
#O.cell.Hsize=Matrix3(0.1,0,0, 0,0.1,0, 0,0,0.1)
sp=pack.SpherePack()
radius=5e-3
num=sp.makeCloud(Vector3().Zero,O.cell.refSize,radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic
O.bodies.append([utils.sphere(s[0],s[1]) for s in sp])

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()],nBins=5,sweepLength=.05*radius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 PeriTriaxController(dynCell=True,mass=0.2,maxUnbalanced=0.01,relStressTol=0.02,goal=(-1e4,-1e4,0),stressMask=3,globUpdate=5,maxStrainRate=(1.,1.,1.),doneHook='triaxDone()',label='triax'),
 NewtonIntegrator(damping=.2),
]
O.dt=utils.PWaveTimeStep()
O.run();
qt.View()

phase=0
def triaxDone():
 global phase
 if phase==0:
  print 'Here we are: stress',triax.stress,'strain',triax.strain,'stiffness',triax.stiff
  print 'Now 蔚z will go from 0 to .2 while 蟽x and 蟽y will be kept the same.'
  print 'Porosity',utils.porosity
  triax.goal=(-1e4,-1e4,-0.2)
  phase+=1
 elif phase==1:
  print 'Here we are: stress',triax.stress,'strain',triax.strain,'stiffness',triax.stiff
  print 'Porosity',utils.porosity
  print 'Done, pausing now.'
  O.pause()

The problem is the compression will go on endlessly as there is no stress generated to reach the goal.

Ning

Jan Stránský (honzik) said : #4

try

utils.porosity()

instead of

utils.porosity

as already suggested Bruno

Jan

ceguo píše v St 07. 03. 2012 v 12:30 +0000:
> Question #189913 on Yade changed:
> https://answers.launchpad.net/yade/+question/189913
>
> Status: Answered => Open
>
> ceguo is still having a problem:
> Thanks,
>
> I notice the missing of brackets for a function. Regarding to the last problem, I use the example script just with modification from Dem3DofGeom to ScGeom (Dem3DofGeom works, but as I read the doc which suggest us using ScGeom instead.)
> #============================================================
> "Test and demonstrate use of PeriTriaxController."
> from yade import *
> from yade import pack,log,qt
> log.setLevel('PeriTriaxController',log.TRACE)
> O.periodic=True
> O.cell.refSize=Vector3(.1,.1,.1)
> #O.cell.Hsize=Matrix3(0.1,0,0, 0,0.1,0, 0,0,0.1)
> sp=pack.SpherePack()
> radius=5e-3
> num=sp.makeCloud(Vector3().Zero,O.cell.refSize,radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic
> O.bodies.append([utils.sphere(s[0],s[1]) for s in sp])
>
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb()],nBins=5,sweepLength=.05*radius),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom()],
> [Ip2_FrictMat_FrictMat_FrictPhys()],
> [Law2_ScGeom_FrictPhys_CundallStrack()]
> ),
> PeriTriaxController(dynCell=True,mass=0.2,maxUnbalanced=0.01,relStressTol=0.02,goal=(-1e4,-1e4,0),stressMask=3,globUpdate=5,maxStrainRate=(1.,1.,1.),doneHook='triaxDone()',label='triax'),
> NewtonIntegrator(damping=.2),
> ]
> O.dt=utils.PWaveTimeStep()
> O.run();
> qt.View()
>
> phase=0
> def triaxDone():
> global phase
> if phase==0:
> print 'Here we are: stress',triax.stress,'strain',triax.strain,'stiffness',triax.stiff
> print 'Now 蔚z will go from 0 to .2 while 蟽x and 蟽y will be kept the same.'
> print 'Porosity',utils.porosity
> triax.goal=(-1e4,-1e4,-0.2)
> phase+=1
> elif phase==1:
> print 'Here we are: stress',triax.stress,'strain',triax.strain,'stiffness',triax.stiff
> print 'Porosity',utils.porosity
> print 'Done, pausing now.'
> O.pause()
>
> The problem is the compression will go on endlessly as there is no stress generated to reach the goal.
>
> Ning
>

Jérôme Duriez (jduriez) said : #5

Hello,

I guess you probably broke something in the contact law computations with your change between De3DofGeom and ScGeom, so that contact forces at the boundaries are not computed in your example => they remain nill..

The point is that I do not know PeriTriaxController (how are "boundaries" taken into account ?), so that I can not help you more...

Jérôme Duriez (jduriez) said : #6

Or, give also a look at
https://www.yade-dem.org/doc/yade.wrapper.html#yade.wrapper.PeriTriaxController.reversedForces

Maybe were you actually performing an isotropic traction ??

ceguo (hhh-guo) said : #7

Thank you all.

The first problem is due to the missing of brackets when I call a function as suggested by Luc, Bruno and Jan.
The last problem is solved after adding `reversedForces=True` in PeriTriaxController because ScGeom uses a reversed sign for forces, as pointed out by Jduriez.

Ning

Try (see the additional "reversedForces" attribute):
 PeriTriaxController(dynCell=True,mass=0.2,maxUnbalanced=0.01,relStressTol=0.02,goal=(-1e4,-1e4,0),stressMask=3,globUpdate=5,maxStrainRate=(1.,1.,1.),doneHook='triaxDone()',label='triax',reversedForces=True)

A sign convention difference between different functors. We already concluded that reversedForces should be true by default a while ago. This annoying "feature" will be fixed, eventually...

Oooops, I didn't update the page, Jerome has been faster.