[Paraview] Newbie question on python filter and integrating over a dimension

[Craig Michoski]

Hi all,

I have a super newbie question, I am trying to properly learn how to
use the programmable python filter to integrate along dimensions of my
arrays.  The below seems to work, but it is extremely slow, and I'm
sure the dumbest thing I can be doing.

What I have is a scalar field in 2d called alpha_0, and a vector field
in 2d called sigma_0_0__sigma_0_1, and I want to integrate the term:

f(alpha)*g(sigma)


along a single dimension of my mesh.  In latex form, if my domain is
(x,y) \in Omega\subset R^2, then I want to compute for every x: int_y
f*g dy.  The terms f and g are defined in more detail below in the
script.

The best solution would include the one-diemsnional output graph as
well, whereas now, I just compute the values at each point and output
to print.  Any advice would be very appreciated.

Cheers,
Craig


--------------------

import numpy

alpha0 = inputs[0].PointData['alpha_0']
alpha1 = inputs[0].PointData['alpha_0']
sigma0 = inputs[0].PointData['sigma_0_0__sigma_0_1']

pdi = self.GetInputDataObject(0,0)

bg_field = 1
pdi = self.GetInputDataObject(0,0)
pdo = self.GetOutputDataObject(0)
numPoints = pdi.GetNumberOfPoints()
import time

ind=self.GetInput()
iin=ind.GetInformation()
timestep = time.clock()

dx = 100.0
dy = 100.0

# Find the dx and dy for integration on the uniform mesh (there must
be a better way?)
for (j,alpha) in enumerate(alpha0):
    x1,y1,z1 = pdi.GetPoint(j)[:3]
    x2,y2,z2 = pdi.GetPoint(j+1)[:3]
    if abs(x2-x1) <= dx and abs(x2-x1)>0:
        dx = x2 - x1
    if abs(y2-y1) <= dy and abs(y2-y1)>1e-6:
        dy = y2 - y1

counter = 0
last_step = -1.0
prev_x = 0.0
efold = numpy.zeros(numPoints)
y_first = -1.0


#integrate the term using a simple Reimann sum:
int_at_each_x_coordinate_along_y
# (1/exp(2*alpha))*sigma_x_component dy (this is very slow!)

for (j,alpha) in enumerate(alpha0):
    x,y,z = pdi.GetPoint(j)
    if x<380.0 and x>70.0 and y==0.0:
        y_first = []
        for (i,alpha_1) in enumerate(alpha1):
            x1,y1,z1 = pdi.GetPoint(i)
            if x1 == x and y1 == y:
                for (k,sigmax) in enumerate(sigma0):
                    x2,y2,z2 = pdi.GetPoint(i)
                    if x2 == x1: #and y2 == y1:
                        if y2 not in y_first:
                            efold[j] +=  ( 1.0/exp(2.0*alpha_1) ) *
sigmax[0] * dy  #exp(alpha_1)*sigmax[0]*dy
                            y_first.append(y2)

                if x!=last_step:
                    counter = counter + 1
                    print "y_integrated_chi_at1(",floor(timestep),",",
counter,") = ", efold[j]/377.0, ";"
                    print "x(",counter,") = ", x, ";"
                    last_step = x