TesselationWrapper and porosity

You can get a map of porosity with TW, but not a map of porosity gradient.
Porosity gradient at the particle scale is a very unusual quantity, I don't know how it could be defined.
Bruno

map of porosity is what I'd like to get. How can I get this by means of TW?

TW.volume(i) is the volume assigned to one sphere. From this you can
easily define the corresponding porosity.

B

I see, but I cannot get how this can be mapped or exported to vtk or whatever
Please check out the following script.
Is this what you'r talking about? Then what do I do with 'density matrix'?

Cheers,

from yade import utils, plot
from yade import pack, qt
from datetime import datetime

#================= define the materials =======================

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= False, young=6.81e5,
density=1523.6, poisson=0.3, frictionAngle= 20, fragile=True, label='wall'))

#========= creating walls ======================

walls=aabbWalls([(-0.4,-0.4,0.0),(0.4,0.4,1.0)],thickness=0.01,oversizeFactor=1.0,material='wall')
wallIds=O.bodies.append(walls)

#==================================== make the binder ==================================

pred1=pack.inAlignedBox((-0.4,-0.4,0.0),(0.4,0.4,0.32))
O.bodies.append(pack.regularOrtho(pred1,radius=3e-2,gap=1e-4, material='wall'))

pred2=pack.inAlignedBox((-0.4,-0.4,0.31),(0.4,0.4,0.65))
O.bodies.append(pack.regularOrtho(pred2,radius=3e-2,gap=-1e-4, material='wall'))

pred3=pack.inAlignedBox((-0.4,-0.4,0.64),(0.4,0.4,0.999))
O.bodies.append(pack.regularOrtho(pred3,radius=3e-2,gap=-6e-3, material='wall'))

j=0
for x in range(len(O.bodies)):
 if (O.bodies[x]):
  if isinstance(O.bodies[x].shape,Sphere):
   if (O.bodies[x].isStandalone):
    j=j+1
print ''
print ''
print 'standalone ones are *** ',j

v_s=j*1.3333333333*pi*(0.03**3)
############################
### DEFINING ENGINES ###
############################

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop([Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
 [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True, setCohesionOnNewContacts=True)],
 [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 ),
 NewtonIntegrator(damping=0.4,gravity=[0,0,-10])
]

O.dt=3e-6
#===============================

TW=TesselationWrapper()
TW.triangulate()
TW.computeVolumes()

density=[]

for i in O.bodies:
 if isinstance(i.shape,Sphere):
  density.append((1.13097e-4)/TW.volume(i.id))

I am assuming that you know how to export any value attached to particles using VTK exporter (else please see the doc, this one is well documented).
Then the only question is how to compute porosity for one particle.
I don't understand the equation (1.13097e-4)/TW.volume(i.id)
If Vi the volume of the sphere "i", then I would think the porosity is [1 - Vi/TW.volume(i.id)]

Bruno

yes the equation is fine. How to export it however, is not clear to me.
I'd like to get the triangulated scene with this parameter (porosity) assigned to them.

The following script does not work!

v_s=1.13e-4 # this is the volume of the spheres (it's a monosize assembly)

O.engines=O.engines+[PyRunner(command="vtkExporter.exportSpheres(what=['porosity','1-v_s/TW.volume(i.id)'])", iterPeriod=100)]

TW=TesselationWrapper()
TW.triangulate()
TW.computeVolumes()
vtkExporter=export.VTKExporter('test')

I understand that to visualize the deformation for example we can write:

#TW.setState(0)
#O.run(10,True)
#TW.setState(1)
#TW.defToVtk("strainvtk")

But there's not function like PorosityToVtk and documentation is not helping me to see how can I make one.

And if I need to use vtk exporter along with TesselatioinWrapper, the link between them is not clear and there's no example or something in documentation to help

I even tried the following. The problem is that vtkExporter is not accepting any function from outside. In the following script I get the error of TW is not defined.
Then, I tried to put the values of TW.volume(x) in a list named vols. Then, again, vtkExporter is not accepting the list vols and I get the error of vols is not defined.

TW=TesselationWrapper()
TW.triangulate()
TW.computeVolumes()

vtkExporter = export.VTKExporter('Test')
vtkExporter.exportSpheres(what=[('myParameter','b.state.mass*TW.volume(b.id)')])

It must be a scope issue, try something around globals()["TW"].volume(id)?

Else it should be possible to add the attribute:
b.volume=TW.volume(b.id) #then use it to export

I just find no way to get this thing working. Even with globals it gives an error!

TW=TesselationWrapper()
TW.triangulate()
TW.computeVolumes()

global xy
xy=[]

for i in O.bodies:
 if isinstance(i.shape,Sphere):
  xy.append(i.state.pos[0])

vtkExporter = export.VTKExporter('Rep')
vtkExporter.exportSpheres(what=[('xDist','globals()["xy"][b.id]')])

output:

Traceback (most recent call last):
  File "/home/bemaj3/sw/bin/yade2015", line 182, in runScript
    execfile(script,globals())
  File "51.py", line 84, in <module>
    vtkExporter.exportSpheres(what=[('xDist','globals()["xy"][b.id]')])
  File "/home/bemaj3/old-trunk/build/lib/x86_64-linux-gnu/yade-Unknown/py/yade/export.py", line 431, in exportSpheres
    test = eval(command) # ... eval one example to see what type (float, Vector3, Matrix3) the result is ...
  File "<string>", line 1, in <module>
KeyError: 'xy'

This:

global xy
xy=[]

is just not doing what you expect. There is no 'xy' in globals after that:
KeyError: 'xy'

Did you try the second method?

Bruno

I don't know what do you mean by second method!

b.volume=TW.volume(b.id) : This, as far as I understand, is the same as the other. We're assigning the TW parameters to a new parameter and then we take that one to export.

In the user manual for tesselationWrapper we've got this as well:

#"b" being a body
TW=TesselationWrapper()
TW.computeVolumes()
s=bodyStressTensors()
stress = s[b.id]**4.*pi/3.*b.shape.radius**3/TW.volume(b.id)

But, the main part has not been given. The question is how we can export a user defined parameter, stress, porosity or whatever the same way as the strain?

b.volume=TW.volume(b.id) is not the same in terms of scope, and the problems you have are mainly scope problems.

>"The question is how we can export a user defined parameter the same way as the strain?"
We cannot, since there cannot be a specific export function implemented in advance for every possible _user defined_ parameter.

Bruno

OK. Thank you Bruno

Hey behzad,

        Do you solve the problem?

Weimin,

No! I could not find a solution finally.