[Paraview] Programmable filter issue with modifying grid data
Burlen Loring
bloring at lbl.gov
Fri Oct 28 19:23:14 EDT 2011
> I actually am increasing my extents (extent[1]+1 for example) to
> compensate for this,
yes you added one. that is needed to compute the number of cells from
the extents. you still need to add one to the result (the number of
cells) to get the number of points. As is, you are not looking up the
values you think you are in your input array. You won't crash since
you'll stay within the point array bounds.
for cells compute the number of them like this:
ncx = ext[1]-ext[0]+1
ncy = ext[3]-ext[2]+1
ncxy = ncx*ncy
ncz = ext[3]-ext[2]+1
and in your loop index a cell at i,j,k like this:
cidx = k*ncxy + j*ncx+i
For points compute the number of them like this:
npx = ncx+1
npy = ncy+1
npxy = npx*npy
npz = ncz+1
and in your loop index a point at i,j,k like this:
pidx = k*npxy + j*npx+i
> I agree that i should be the fastest changing index in some sense, but
> because I have to step over all values of i,j,k to read the entire
> array into my python array, I don't see why any other looping
> direction would be quicker?
This has to do with the cache inside your cpu. When you use a memory
location the cpu automatically loads the adjacent locations(a cache
line) into the cache. In memory i direction changes fastest, meaning
location (i+1,j,k) is adjacent following location (i,j,k). when you
access the adjacent memory locations in your loops sequentially your
data are more likely to be in the cache and hence a faster. On the other
hand memory location (i,j,k+1) is npx*npy elements away from location
(i,j,k). That's usually quite a number of cache lines away and hence
less likely to already be in the cache resulting in a cache miss.
Execution stalls while the cpu loads the cache line. the cache is fairly
small so a cache miss likely results in some data being evicted to make
room for the new data. When this happens on every array access it's
called thrashing.This is what happens when you use the wrong order in
your loops.
> I think you're right, it does seem to construct a list explicitly, but
> I don't really know what the alternative is?
I am not a great python programmer, but one way would be to use a while
loop:
k=ext[4]
while(k<=ext[5]):
#
j=ext[2]
while(j<=ext[3]):
#
i=ext[0]
while(i<=ext[1]):
#
pidx = k*npxy + j*npx+i
# do something here
i+=1
#
j+=1
#
k+=1
#
On 10/28/2011 03:27 PM, stofem wrote:
> Hi,
>
> Thanks for your help both of you.
>
> Cheers for pointing out my Cell-Point mistake Burlen, I actually am
> increasing my extents (extent[1]+1 for example) to compensate for
> this, but I did make a mistake in my final call
>
> *out.GetCellData().AddArray(newArray)*
> *
> *
> which should be getting PointData, not CellData. My filter now works
> as planned :)
>
> However I am interested in some of your comments. Being a new
> Python/Paraview/VTK programmer I would appreciate a bit of further advice:
>
> In the two first loops you are treating 'i' as the slowest
> changing direction when it is the fastest. That will kill your
> performance. Last loop has correct order.
>
>
> I agree that i should be the fastest changing index in some sense, but
> because I have to step over all values of i,j,k to read the entire
> array into my python array, I don't see why any other looping
> direction would be quicker?
>
> You may want to avoid the range function in your loop since if my
> recollection is correct that explicitly constructs a list of values
>
>
> I think you're right, it does seem to construct a list explicitly, but
> I don't really know what the alternative is?
>
> Thanks again.
>
> Looking forward to seeing that new filter built in!
>
> On Sat, Oct 29, 2011 at 3:08 AM, Burlen Loring <bloring at lbl.gov
> <mailto:bloring at lbl.gov>> wrote:
>
> Hi,
>
> You haven't configured your output, assuming rectilinear grid you
> need to set the extent and provide coordinate arrays.
>
> You are indexing a point based array with cell based index, so the
> lookup you make into vtkVarray is incorrect. 'Extent' tells you
> the number of cells, not the number of points, for that use
> 'Dimensions' or add one in each direction. Your script would
> certainly be more readable if you used some variables such as:
>
> ncx = ext[1]-ext[0]+1
> ncy = ext[3] -ext[2] +1
> ncxy = nx*ny
> ncz = ext[5]-ext[4] +1
> cidx = k*nxy + j*nx+i
>
> for number of cells and a cell array index.
>
> In the two first loops you are treating 'i' as the slowest
> changing direction when it is the fastest. That will kill your
> performance. Last loop has correct order.
>
> You may want to avoid the range function in your loop since if my
> recollection is correct that explicitly constructs a list of values.
>
> I hope this is helpful.
> Burlen
>
>
> On 10/26/2011 10:41 PM, Mr FancyPants wrote:
>> Hi there,
>>
>> I am trying to write a programmable filter which will take my
>> RectilinearGrid data and modify one of the data
>> arrays. Specifically I want to sum over one dimension. I have
>> written something to do this, but my data comes out garbled.
>>
>> Below is my code. Basically all I am doing is extracting all the
>> data, doing my sum over the z direction and then putting this
>> data into another rectilinear grid.
>>
>>
>> *data=self.GetInput()*
>> *out=self.GetOutput()*
>> *extent=data.GetExtent()*
>> *vtkVarray=data.GetPointData().GetArray('velocity')*
>> *
>> *
>> *import numpy*
>> *pyVarray_x = numpy.zeros([extent[1]+1,extent[3]+1,extent[5]+1])*
>> *pyVarray_y = numpy.zeros([extent[1]+1,extent[3]+1,extent[5]+1])*
>> *pyVarray_z = numpy.zeros([extent[1]+1,extent[3]+1,extent[5]+1])*
>> *
>> *
>> *for i in range(0,extent[1]+1):*
>> *
>> for j in range(0,extent[3]+1):*
>> *
>> for k in range(0,extent[5]+1):*
>> *
>> tup=vtkVarray.GetTuple(i+j*(extent[1]+1)+k*(extent[1]+1)*(extent[3]+1))*
>> *
>> pyVarray_x[i,j,k]=tup[0]*
>> *
>> pyVarray_y[i,j,k]=tup[1]*
>> *
>> pyVarray_z[i,j,k]=tup[2]*
>> *
>> *
>> *
>> *
>> *
>> *
>> *pyVarray_x_noz = numpy.zeros([extent[1]+1,extent[3]+1,1])*
>> *pyVarray_y_noz = numpy.zeros([extent[1]+1,extent[3]+1,1])*
>> *pyVarray_z_noz = numpy.zeros([extent[1]+1,extent[3]+1,1])*
>> *
>> *
>> *for i in range(0,extent[1]+1):*
>> *
>> for j in range(0,extent[3]+1):*
>> *
>> pyVarray_x_noz[i,j]=pyVarray_x[i,j,:].sum()*
>> *
>> pyVarray_y_noz[i,j]=pyVarray_y[i,j,:].sum()*
>> *
>> pyVarray_z_noz[i,j]=pyVarray_z[i,j,:].sum()*
>> *
>> *
>> *newArray=vtk.vtkDoubleArray()*
>> *newArray.SetName("test")*
>> *newArray.SetNumberOfComponents(3)*
>> *
>> *
>> *print newArray*
>> *
>> *
>> *for k in range(0,extent[5]+1):*
>> *
>> for j in range(0,extent[3]+1):*
>> *
>> for i in range(0,extent[1]+1):*
>> *
>> tup0=pyVarray_x_noz[i,j]*
>> *tup1=pyVarray_y_noz[i,j]*
>> *
>> tup2=pyVarray_z_noz[i,j]*
>> *
>> tup=(tup0,tup1,tup2)*
>> *
>> newArray.InsertNextTuple((tup0,tup1,tup2))*
>> *
>> *
>> *print newArray*
>> *print vtkVarray*
>> *
>> *
>> *out.GetCellData().AddArray(newArray)*
>> *
>> *
>> I have no idea whats going wrong with this. Any help is much
>> appreciated, new to using paraview.
>>
>> Thanks, James
>>
>>
>> _______________________________________________
>> Powered bywww.kitware.com <http://www.kitware.com>
>> Visit other Kitware open-source projects at
>> http://www.kitware.com/opensource/opensource.html Please keep
>> messages on-topic and check the ParaView Wiki at:
>> http://paraview.org/Wiki/ParaView Follow this link to
>> subscribe/unsubscribe:
>> http://www.paraview.org/mailman/listinfo/paraview
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.paraview.org/pipermail/paraview/attachments/20111028/377c2be6/attachment-0001.htm>
More information about the ParaView
mailing list