vtkProbeFilter getting stuck in Python script

Hello,

I have a 3-D unstructured mesh of a sphere (a model of the earth) in vtu format (a .pvtu file with 16 associated .vtu files). I would like to interpolate these data onto a spherical shell of a given radius, to extract a depth slice of the data. I’m trying to use a vtkProbeFilterfilter in Python to do this, but it’s giving me problems. See below code. In particular, it’s getting stuck on the probe.Update() command. I can get it to work on small models, but on the bigger ones I’m actually interested in it doesn’t move past that line, even if I wait for hours. I’ve also tried using an vtkPointInterpolator filter to achieve the same goals, which produces output in a reasonable timeframe, but it gives a weirdly discontinuous field that I think has to do with the unstructured mesh. Any suggestions for how to do this interpolation in as efficient a manner as possible?

Thanks,
Sam

import vtk
from vtk.util.numpy_support import vtk_to_numpy
import numpy as np

dx = 0.1 # Resolution to interpolate onto (degrees)

# Read in pvtu file
reader = vtk.vtkXMLPUnstructuredGridReader()
reader.SetFileName(folder + '/solution/solution-00000.pvtu')
reader.Update()

# Define sphere at radius R to interpolate onto
R = 6371.e3-200.e3
sphere = vtk.vtkSphereSource()
sphere.SetCenter(0.0,0.0,0.0)
sphere.SetRadius(R)
sphere.SetPhiResolution(int(180.0/dx)+1)
sphere.SetThetaResolution(int(360.0/dx)+1)
sphere.Update()

# Interpolate data onto sphere
probe = vtk.vtkProbeFilter()
probe.SetInputConnection(sphere.GetOutputPort())
probe.SetSourceData(reader.GetOutput())
probe.Update()

# Convert to numpy arrays
V = vtk_to_numpy(probe.GetOutput().GetPointData().GetArray('velocity'))
p = vtk_to_numpy(probe.GetOutput().GetPointData().GetArray('p'))
points = vtk_to_numpy(probe.GetOutput().GetPoints().GetData())

# Split arrays into vectors
x = points[:,0]
y = points[:,1]
z = points[:,2]
del points

# Convert to spherical coordinates
pi = 3.14159265358979323846
r = (x**2.0+y**2.0+z**2.0)**0.5
lon = np.arctan2(y,x)*180.0/pi
del x
del y
lat = np.arcsin(z/r)*180.0/pi
del z