ParaView Voidfraction Filter by richti83

ntry=5000 #increase for better accuracy, decrease for faster performance

from random import uniform
input1 = self.GetInputDataObject(0, 0)
output = self.GetOutputDataObject(0)
b=input1.GetBounds()
xlow,xhi,ylow,yhi,zlow,zhi=b[:6]
P=input1.GetPointData()
COUNT=input1.GetNumberOfPoints()
#step 1: estimate bounding volume
for i in range(0,COUNT):
    coord = input1.GetPoint(i)
    x, y, z = coord[:3]
    r=P.GetArray("radius").GetValue(i)
    if x<=xlow:
        xlow=x-r
    if y<=ylow:
        ylow=y-r
    if z<=zlow:
        zlow=z-r
    if x>=xhi:
        xhi=x+r
    if y>=yhi:
        yhi=y+r
    if z>=zhi:
        zhi=z+r
Vb=(xhi-xlow)*(yhi-ylow)*(zhi-zlow)

#Step 2: Monte Carlo Integration of the volume of the spheres
hit=0.0

for t in range(0,ntry):
    xtry=uniform(xlow,xhi)
    ytry=uniform(ylow,yhi)
    ztry=uniform(zlow,zhi)
    known=0
    for i in range(0,COUNT):    #for all centers of spheres check distance to X_try
     coord = input1.GetPoint(i)
     x, y, z = coord[:3]
     r=P.GetArray("radius").GetValue(i)
     d=0.0
     d+=(xtry-x)*(xtry-x)
     d+=(ytry-y)*(ytry-y)
     d+=(ztry-z)*(ztry-z)
     if (d<r*r):
        #point is inside this sphere
        hit+=1.0
        known=1
     if (known==1):
        break

vol=hit/float(ntry)*Vb
void=(Vb-vol)*100/Vb

outputarray = vtk.vtkStringArray()
outputarray.SetName("Text")
outputarray.SetNumberOfTuples(1)
outputarray.SetValue(0,"Spheres=%d\nV_S=%1.6f m3\nV_B=%1.6f m3\nVoidfraction=%1.1f%%" % (COUNT,vol,Vb,void))
output.GetRowData().AddArray(outputarray)