import k3d
import vtk
from vtk.util import numpy_support
import numpy as np
plot = k3d.plot(camera_auto_fit=False)
plot.display()
Output()
reader = vtk.vtkSTLReader()
reader.SetFileName('./cfd/c0006.stl')
reader.Update()
geometry_wss = reader.GetOutput()
reader = vtk.vtkXMLImageDataReader()
reader.SetFileName('./cfd/000005595.vti')
reader.Update()
vti = reader.GetOutput()
bounds = vti.GetBounds()
shape = vti.GetDimensions()
spacing = vti.GetSpacing()
seeds_points = vtk.vtkPoints()
velocity = numpy_support.vtk_to_numpy(vti.GetPointData().GetArray('v [m/s]'))
velocity = velocity.flatten().reshape(shape[2], shape[1], shape[0], 3)
speed = np.linalg.norm(velocity, axis=3)
speed[np.isnan(speed)] = 0.0
points = vtk.vtkPoints()
points.SetNumberOfPoints(vti.GetNumberOfPoints())
for i in range(points.GetNumberOfPoints()):
points.SetPoint(i, vti.GetPoint(i))
grid = vtk.vtkStructuredGrid()
grid.SetDimensions(shape)
grid.SetPoints(points)
grid.GetPointData().SetVectors(
vti.GetPointData().GetVectors('v [m/s]')
)
seeds = vtk.vtkPolyData()
points = vtk.vtkPoints()
for (z,x) in zip(*np.where(speed[:,0,:] > 0.0)):
y = 0
if x % 4 == 0 and z % 4 == 0:
points.InsertNextPoint((x * spacing[0] + bounds[0],
y * spacing[1] + bounds[2],
z * spacing[2] + bounds[4]))
seeds.SetPoints(points)
streamer = vtk.vtkStreamTracer()
streamer.SetInputData(grid)
streamer.SetSourceData(seeds)
streamer.SetIntegrationDirectionToForward()
streamer.SetComputeVorticity(0)
streamer.SetIntegrator(vtk.vtkRungeKutta4())
streamer.Update()
streamline = streamer.GetOutput()
streamlines_points = numpy_support.vtk_to_numpy(streamline.GetPoints().GetData())
streamlines_velocity = numpy_support.vtk_to_numpy(streamline.GetPointData().GetArray('v [m/s]'))
streamlines_speed = np.linalg.norm(streamlines_velocity, axis=1)
C:\Tools\Anaconda\lib\site-packages\numpy\linalg\linalg.py:2286: RuntimeWarning: overflow encountered in multiply s = (x.conj() * x).real
vtkLines = streamline.GetLines()
vtkLines.InitTraversal()
point_list = vtk.vtkIdList()
lines = []
lines_attributes = []
while vtkLines.GetNextCell(point_list):
start_id = point_list.GetId(0)
end_id = point_list.GetId(point_list.GetNumberOfIds() - 1)
l= []
v= []
for i in range(start_id, end_id):
l.append(streamlines_points[i])
v.append(streamlines_speed[i])
lines.append(np.array(l))
lines_attributes.append(np.array(v))
cm = k3d.matplotlib_color_maps.Inferno
for l, a in zip(lines, lines_attributes):
plot += k3d.line(l, attribute=a, width=0.0001,
color_map=cm, color_range=[0,0.5],
shader='mesh',
compression_level=9)
plot += k3d.vtk_poly_data(geometry_wss, opacity=0.1, wireframe=True, color=0xaaaaaa)