1#!/usr/bin/env python 2import vtk 3from vtk.util.misc import vtkGetDataRoot 4VTK_DATA_ROOT = vtkGetDataRoot() 5 6# 7# The dataset read by this exercise ("combVectors.vtk") has field data 8# associated with the pointdata, namely two vector fields. In this exercise, 9# you will convert both sets of field data into attribute data. Mappers only 10# process attribute data, not field data. So we must convert the field data to 11# attribute data in order to display it. (You'll need to determine the "names" 12# of the two vector fields in the field data.) 13# 14# If there is time remaining, you might consider adding a programmable filter 15# to convert the two sets of vectors into a single scalar field, representing 16# the angle between the two vector fields. 17# 18# You will most likely use vtkFieldDataToAttributeDataFilter, vtkHedgeHog, 19# and vtkProgrammableAttributeDataFilter. 20# 21# Create the RenderWindow, Renderer and interactor 22# 23ren1 = vtk.vtkRenderer() 24renWin = vtk.vtkRenderWindow() 25renWin.SetMultiSamples(0) 26renWin.AddRenderer(ren1) 27iren = vtk.vtkRenderWindowInteractor() 28iren.SetRenderWindow(renWin) 29# create pipeline 30# 31# get the pressure gradient vector field 32pl3d_gradient = vtk.vtkMultiBlockPLOT3DReader() 33pl3d_gradient.SetXYZFileName("" + str(VTK_DATA_ROOT) + "/Data/combxyz.bin") 34pl3d_gradient.SetQFileName("" + str(VTK_DATA_ROOT) + "/Data/combq.bin") 35pl3d_gradient.SetScalarFunctionNumber(100) 36pl3d_gradient.SetVectorFunctionNumber(210) 37pl3d_gradient.Update() 38pl3d_g_output = pl3d_gradient.GetOutput().GetBlock(0) 39# get the velocity vector field 40pl3d_velocity = vtk.vtkMultiBlockPLOT3DReader() 41pl3d_velocity.SetXYZFileName("" + str(VTK_DATA_ROOT) + "/Data/combxyz.bin") 42pl3d_velocity.SetQFileName("" + str(VTK_DATA_ROOT) + "/Data/combq.bin") 43pl3d_velocity.SetScalarFunctionNumber(100) 44pl3d_velocity.SetVectorFunctionNumber(200) 45pl3d_velocity.Update() 46pl3d_v_output = pl3d_velocity.GetOutput().GetBlock(0) 47# contour the scalar fields 48contour = vtk.vtkContourFilter() 49contour.SetInputData(pl3d_g_output) 50contour.SetValue(0,0.225) 51# probe the vector fields to get data at the contour surface 52probe_gradient = vtk.vtkProbeFilter() 53probe_gradient.SetInputConnection(contour.GetOutputPort()) 54probe_gradient.SetSourceData(pl3d_g_output) 55probe_velocity = vtk.vtkProbeFilter() 56probe_velocity.SetInputConnection(contour.GetOutputPort()) 57probe_velocity.SetSourceData(pl3d_v_output) 58# 59# To display the vector fields, we use vtkHedgeHog to create lines. 60# 61velocity = vtk.vtkHedgeHog() 62velocity.SetInputConnection(probe_velocity.GetOutputPort()) 63velocity.SetScaleFactor(0.0015) 64pressureGradient = vtk.vtkHedgeHog() 65pressureGradient.SetInputConnection(probe_gradient.GetOutputPort()) 66pressureGradient.SetScaleFactor(0.00002) 67def ExecuteDot (__vtk__temp0=0,__vtk__temp1=0): 68 # proc for ProgrammableAttributeDataFilter. Note the use of "double()" 69 # in the calculations. This protects us from Python using ints and 70 # overflowing. 71 inputs = dotProduct.GetInputList() 72 input0 = inputs.GetDataSet(0) 73 input1 = inputs.GetDataSet(1) 74 numPts = input0.GetNumberOfPoints() 75 vectors0 = input0.GetPointData().GetVectors() 76 vectors1 = input1.GetPointData().GetVectors() 77 scalars = vtk.vtkFloatArray() 78 i = 0 79 while i < numPts: 80 v0 = vectors0.GetTuple3(i) 81 v1 = vectors1.GetTuple3(i) 82 v0x = lindex(v0,0) 83 v0y = lindex(v0,1) 84 v0z = lindex(v0,2) 85 v1x = lindex(v1,0) 86 v1y = lindex(v1,1) 87 v1z = lindex(v1,2) 88 l0 = expr.expr(globals(), locals(),["double","(","v0x",")*","double","(","v0x",")","+","double","(","v0y",")*","double","(","v0y",")","+","double","(","v0z",")*","double","(","v0z",")"]) 89 l1 = expr.expr(globals(), locals(),["double","(","v1x",")*","double","(","v1x",")","+","double","(","v1y",")*","double","(","v1y",")","+","double","(","v1z",")*","double","(","v1z",")"]) 90 l0 = expr.expr(globals(), locals(),["sqrt","(","double","(","l0","))"]) 91 l1 = expr.expr(globals(), locals(),["sqrt","(","double","(","l1","))"]) 92 if (l0 > 0.0 and l1 > 0.0): 93 d = expr.expr(globals(), locals(),["(","double","(","v0x",")*","double","(","v1x",")","+","double","(","v0y",")*","double","(","v1y",")","+","double","(","v0z",")*","double","(","v1z","))/(","l0","*","l1",")"]) 94 pass 95 else: 96 d = 0.0 97 pass 98 scalars.InsertValue(i,d) 99 i = i + 1 100 101 dotProduct.GetOutput().GetPointData().SetScalars(scalars) 102 del scalars 103 104# 105# We use the ProgrammableAttributeDataFilter to compute the cosine 106# of the angle between the two vector fields (i.e. the dot product 107# normalized by the product of the vector lengths). 108# 109# 110dotProduct = vtk.vtkProgrammableAttributeDataFilter() 111dotProduct.SetInputConnection(probe_velocity.GetOutputPort()) 112dotProduct.AddInput(probe_velocity.GetOutput()) 113dotProduct.AddInput(probe_gradient.GetOutput()) 114dotProduct.SetExecuteMethod(ExecuteDot) 115# 116# Create the mappers and actors. Note the call to GetPolyDataOutput when 117# setting up the mapper for the ProgrammableAttributeDataFilter 118# 119velocityMapper = vtk.vtkPolyDataMapper() 120velocityMapper.SetInputConnection(velocity.GetOutputPort()) 121velocityMapper.ScalarVisibilityOff() 122velocityActor = vtk.vtkLODActor() 123velocityActor.SetMapper(velocityMapper) 124velocityActor.SetNumberOfCloudPoints(1000) 125velocityActor.GetProperty().SetColor(1,0,0) 126pressureGradientMapper = vtk.vtkPolyDataMapper() 127pressureGradientMapper.SetInputConnection(pressureGradient.GetOutputPort()) 128pressureGradientMapper.ScalarVisibilityOff() 129pressureGradientActor = vtk.vtkLODActor() 130pressureGradientActor.SetMapper(pressureGradientMapper) 131pressureGradientActor.SetNumberOfCloudPoints(1000) 132pressureGradientActor.GetProperty().SetColor(0,1,0) 133dotMapper = vtk.vtkPolyDataMapper() 134dotMapper.SetInputConnection(dotProduct.GetOutputPort()) 135dotMapper.SetScalarRange(-1,1) 136dotActor = vtk.vtkLODActor() 137dotActor.SetMapper(dotMapper) 138dotActor.SetNumberOfCloudPoints(1000) 139# 140# The PLOT3DReader is used to draw the outline of the original dataset. 141# 142pl3d = vtk.vtkMultiBlockPLOT3DReader() 143pl3d.SetXYZFileName("" + str(VTK_DATA_ROOT) + "/Data/combxyz.bin") 144pl3d.Update() 145pl3d_output = pl3d.GetOutput().GetBlock(0) 146outline = vtk.vtkStructuredGridOutlineFilter() 147outline.SetInputData(pl3d_output) 148outlineMapper = vtk.vtkPolyDataMapper() 149outlineMapper.SetInputConnection(outline.GetOutputPort()) 150outlineActor = vtk.vtkActor() 151outlineActor.SetMapper(outlineMapper) 152outlineActor.GetProperty().SetColor(0,0,0) 153# 154# Add the actors to the renderer, set the background and size 155# 156ren1.AddActor(outlineActor) 157ren1.AddActor(velocityActor) 158ren1.AddActor(pressureGradientActor) 159ren1.AddActor(dotActor) 160ren1.SetBackground(1,1,1) 161renWin.SetSize(500,500) 162#ren1 SetBackground 0.1 0.2 0.4 163cam1 = ren1.GetActiveCamera() 164cam1.SetClippingRange(3.95297,50) 165cam1.SetFocalPoint(9.71821,0.458166,29.3999) 166cam1.SetPosition(-21.6807,-22.6387,35.9759) 167cam1.SetViewUp(-0.0158865,0.293715,0.955761) 168# render the image 169# 170renWin.Render() 171renWin.SetWindowName("Multidimensional Visualization Exercise") 172# prevent the tk window from showing up then start the event loop 173# --- end of script -- 174