Export Local Porosity to vtk file

Hello,

> VoxalPorosity

please use actual names. Both lower/upper case and correct "spelling".
Here it is clear that voxelPorosity is probably what you meant, but it might not be always like that.

> I have tried to use the TesselationWrapper and the VoxalPorosity to calculate it, with no satisfactory results.

please provide the code or at least description what you did and why it was not satisfactory, and which part was not satisfactory (calculation, export, visualization, ...).
voxelPorosity should suit perfectly for this task.

Cheers
Jan

Hello Jan and thanks for coming back to me.

Here is the code I am using right now:

# gravity deposition in box, showing how to plot and save history of data,
# and how to control the simulation while it is running by calling
# python functions from within the simulation loop

# import yade modules that we will use below
from yade import pack, plot
from yade import utils, export
from minieigen import Vector3

# create rectangular box from facets
O.bodies.append(geom.facetBox((.5, .5, .5), (.5, .5, .5), wallMask=31))

# create empty sphere packing
# sphere packing is not equivalent to particles in simulation, it contains only the pure geometry
sp = pack.SpherePack()
# generate randomly spheres with uniform radius distribution
sp.makeCloud((0, 0, 0), (1, 1, 1), rMean=.05, rRelFuzz=.5)
# add the sphere pack to the simulation
sp.toSimulation()

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
        InteractionLoop(
                # handle sphere+sphere and facet+sphere collisions
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.4),
        # call the checkUnbalanced function (defined below) every 2 seconds
        PyRunner(command='checkUnbalanced()', realPeriod=2),
        # call the addPlotData function every 200 steps
        PyRunner(command='addPlotData()', iterPeriod=100),
        VTKRecorder(iterPeriod=2000, recorders=['all'], fileName='/home/antonsie/DEM/Test/M01/-')
]
O.dt = .5 * PWaveTimeStep()

# enable energy tracking; any simulation parts supporting it
# can create and update arbitrary energy types, which can be
# accessed as O.energy['energyName'] subsequently
O.trackEnergy = True

# if the unbalanced forces goes below .05, the packing
# is considered stabilized, therefore we stop collected
# data history and stop
def checkUnbalanced():
 if unbalancedForce() < .05:
  O.pause()
  plot.saveDataTxt('bbb.txt.bz2')
  # plot.saveGnuplot('bbb') is also possible

# collect history of data which will be plotted
def addPlotData():
 # each item is given a names, by which it can be the unsed in plot.plots
 # the **O.energy converts dictionary-like O.energy to plot.addData arguments
 plot.addData(i=O.iter, unbalanced=unbalancedForce(), **O.energy)

# define how to plot data: 'i' (step number) on the x-axis, unbalanced force
# on the left y-axis, all energies on the right y-axis
# (O.energy.keys is function which will be called to get all defined energies)
# None separates left and right y-axis
plot.plots = {'i': ('unbalanced', None, O.energy.keys)}

# show the plot on the screen, and update while the simulation runs
plot.plot()

O.saveTmp()

# overall porosity
pa = utils.voxelPorosity(200,(0, 0, 0),(1, 1, 1))
pb = utils.porosity()
print('porosity of whole sample computed with voxel porosity', pa)
print('porosity of whole sample computed with porosity', pb)

# porosity from sub volumes (8 total)
p1 = utils.voxelPorosity(200,(0, 0, 0),(0.5, 0.5, 0.5))
p2 = utils.voxelPorosity(200,(0, 0.5, 0),(0.5, 1, 0.5))
p3 = utils.voxelPorosity(200,(0.5, 0, 0),(1, 0.5, 0.5))
p4 = utils.voxelPorosity(200,(0.5, 0.5, 0),(1, 1, 0.5))
p5 = utils.voxelPorosity(200,(0, 0, 0.5),(0.5, 0.5, 1))
p6 = utils.voxelPorosity(200,(0.5, 0, 0.5),(1, 0.5, 1))
p7 = utils.voxelPorosity(200,(0, 0.5, 0.5),(0.5, 1, 1))
p8 = utils.voxelPorosity(200,(0.5, 0.5, 0.5),(1, 1, 1))

# print porosity obtained by subvolumina
print(p1, p2, p3, p4, p5, p6, p7, p8)

#vtk = yade.export.VTKExporter('/home/antonsie/DEM/Test/M01/pa-')
#vtk.export('pa')

The code is working fine except for the last line. It returns the porosity calculated with utils.voxelPorosity() for the whole sample. From the documentation I understood, that it is also calculating the porosity for various subvolumina (in this case 200), but it is not clear to me, how I can access this data. Can you help me out here?

So I turned to defining subvolumina in the code which I find to be a tedious process. My resort would be to save them as a txt document and convert it to vtk.

The VTKRecorder can't return the values I would need to calculate the porosity, which are the volumes of the subvolumina defined through utils.voxelPorosity() so I tried to use the VTKExporter.
The yade documentation says, with VTKExporter, the user is able to export any user defined variable. However, it is not clear to me, how to do that. I searched in the documentation and in examples on GitLab. Could you provide more infomation?

Thanks in advance!

Hello,

thanks for the code, now it is much more clear

> except for the last line

see below

> From the documentation I understood, that it is also calculating the porosity for various subvolumina (in this case 200), but it is not clear to me, how I can access this data.

No.
voxelPorosity calculates porosity in given cubic region (determined by coordinates of "start" and "end" corners). It returns one number.
200 here is resolution. In this case, it divides the region into 200x200x200 voxels and for each voxel it determines if it is particle or void. From these values, porosity is computed and returned.
The higher resolution, the more accurate is the result, but also the more RAM and more time is used.

> So I turned to defining subvolumina in the code

This is how I meant it to be used :-)

> which I find to be a tedious process.

In what sense is it tedious?
Code writing of each command separately? If so, use a for loop.

> My resort would be to save them as a txt document and convert it to vtk.

good idea, or I think Paraview can just load the data from .txt file (?)

> so I tried to use the VTKExporter.
> The yade documentation says, with VTKExporter, the user is able to export any user defined variable.
> any variable, but related to existing spheres (taking spheres as the example).

VTKExporter can export any value, but related to existing spheres (if spheres are exported).
So e.g. you can export sphere color, material properties, volume etc.
It is not meant to export literally any variables (e.g. here porosity defined in arbitrary points in space).

I would stick to the txt file first. If it does not work, you can convert the txt file to vtk, or change export directly to vtk (either "manually" as VTKExporter is implemented, or using some external package)

Cheers
Jan