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