Hi
We're written a new Filter for mayavi (1.x) which allows one to plot
arbitrary fields (scalars and vectors) on iso-surfaces. It also gives
the volume fluxes and scalar fluxes through the iso-surface. The
result...some interesting diagnostics and cool pictures :-) Note: flux
fields are calculated using point wise quantities (using the normal at a
vertex) while integrals are performed over triangular elements in the
iso-surface.
For those interested in the hacking aspect, it's interesting because you
have the change the active scalar in a vtk dataset, so there cannot be
one dataset, ie. the active scalar required for the iso-surface is
different from the scalar you want to plot on the iso-surface. The only
way around this appears to be to make a deep copy of the output of the
contour filter, modify the active scalar of this copy and make this
object the output of the filter.
Prabhu, after people have had a change to give feedback, are you happy
with this being committed? We have started playing with mayavi2 (very
nice) so at a later stage we'll port it across.
Enjoy!
Gerard Gorman
Colin Cotter
"""
This module uses an iso-surface as a probe allowing arbitary scalars,
vectors and tensors to be plotted on that iso-surface. Fluxes are also
calculated. This will work for any dataset.
This code is distributed under the conditions of the BSD license. See
LICENSE.txt for details.
Copyright (c) 2001-2002, Prabhu Ramachandran.
"""
__author__ = "Gerard Gorman, Colin Cottor annd Prabhu Ramachandran"
__version__ = "$Revision: 1.0 $"
__date__ = "$Date: 2007/03/08 18:30:13 $"
import Base.Objects, Common
import Tkinter, tkColorChooser, math
import vtk
import vtkPipeline.vtkMethodParser
debug = Common.debug
class IsoSurfaceProbe (Base.Objects.Filter):
""" This module finds an iso-surface of scalar data, and
calculates the normal flux across it due to a vector field. This
will work for any dataset."""
def initialize (self):
debug ("In IsoSurfaceProbe::initialize ()")
Common.state.busy ()
self.cont_fil = vtk.vtkContourFilter ()
self.fil = self.cont_fil
if not self.mod_m.get_scalar_data_name ():
msg = "Warning: No scalar data present to contour!"
Common.print_err (msg)
self.iso = vtk.vtkPolyData ()
self.data_range = self.prev_fil.GetOutput ().GetScalarRange ()
self.SetInput (self.prev_fil)
self.cont_fil.SetValue (0, self.midrange)
attributes = self._get_attribute_list (self.prev_fil.GetOutput ().GetPointData ())
self.scalar_lst = attributes['scalars']
self.visible_scalar = ""
if self.scalar_lst:
self.visible_scalar = self.scalar_lst[0]
else:
self.visible_scalar = "Volume flux"
self._gui_init ()
self.scalar_var = Tkinter.StringVar ()
self.scalar_var.set (self.visible_scalar)
# used for the pipeline browser
self.pipe_objs = self.cont_fil
self.Update ()
Common.state.idle ()
def _gui_init (self):
debug ("In IsoSurfaceProbe::_gui_init ()")
self.slider = None
self.slider_pos = self.midrange
self.contour_var = Tkinter.DoubleVar ()
self.contour_var.set (self.slider_pos)
self.c_entry_var = Tkinter.DoubleVar ()
self.c_entry_var.set (self.slider_pos)
dr = self.mod_m.get_scalar_data_range ()
self.min_cnt = Tkinter.DoubleVar ()
self.min_cnt.set (dr[0])
self.max_cnt = Tkinter.DoubleVar ()
self.max_cnt.set (dr[1])
self.volume_flux_value = Tkinter.DoubleVar ()
self.volume_flux_value.set (0.0)
self.flux_values = {}
for sc in self.scalar_lst:
self.flux_values[sc+"Flux"] = Tkinter.DoubleVar ()
self.flux_values[sc+"Flux"].set (0.0)
self._auto_sweep_init ()
self.sweep_step.set (15)
def SetInput (self, source):
debug ("In IsoSurfaceProbe::SetInput ()")
Common.state.busy ()
self.cont_fil.SetInput (source.GetOutput ())
self.midrange = (self.data_range[0] + self.data_range[1])*0.5
dr = source.GetOutput ().GetScalarRange ()
if (dr[0] != self.data_range[0]) or (dr[1] != self.data_range[1]):
self.data_range = dr
self.min_cnt.set (dr[0])
self.max_cnt.set (dr[1])
self.c_entry_var.set ((dr[0] + dr[1])*0.5)
self.change_contour_limits ()
self.change_contour_entry ()
Common.state.idle ()
def save_config (self, file):
debug ("In IsoSurfaceProbe::save_config ()")
file.write ("%f, %f, %f, '%s'\n"%(self.c_entry_var.get (),
self.min_cnt.get (),
self.max_cnt.get (),
self.scalar_var.get()))
p = vtkPipeline.vtkMethodParser.VtkPickler ()
p.dump (self.cont_fil, file)
def load_config (self, file):
debug ("In IsoSurfaceProbe::load_config ()")
c_val, min_cnt, max_cnt, self.visible_scalar = eval (file.readline ())
self.slider_pos = c_val
self.contour_var.set (c_val)
self.c_entry_var.set (c_val)
self.min_cnt.set (min_cnt)
self.max_cnt.set (max_cnt)
self.scalar_var.set (self.visible_scalar)
p = vtkPipeline.vtkMethodParser.VtkPickler ()
p.load (self.cont_fil, file)
self.change_contour_slider ()
def config_changed (self):
debug ("In IsoSurfaceProbe::config_changed ()")
pass
def make_main_gui (self):
debug ("In IsoSurfaceProbe::make_main_gui ()")
frame = Tkinter.Frame (self.root, relief='ridge', bd=2)
frame.pack (side='top', fill='both', expand=1)
name = "Scalar Variable: " + self.mod_m.get_scalar_data_name ()
dr = (self.min_cnt.get (), self.max_cnt.get ())
resolution = (dr[1] - dr[0])/1000
self.slider = Tkinter.Scale (frame, label=name,
resolution=resolution,
variable=self.contour_var,
from_=dr[0], to=dr[1], length="8c",
orient='horizontal')
self.slider.grid (row=0, column=0, columnspan=2, sticky='ew')
self.slider.bind ("<ButtonRelease>", self.change_contour_slider )
lab = Tkinter.Label (frame, text="Iso-surface value :")
lab.grid (row=1, column=0, sticky='w')
entry = Tkinter.Entry (frame, width=10, relief='sunken',
textvariable=self.c_entry_var)
entry.grid (row=1, column=1, sticky='ew')
entry.bind ("<Return>", self.change_contour_entry)
lab = Tkinter.Label (frame, text="Minimum contour:")
lab.grid (row=2, column=0, sticky='w')
entry = Tkinter.Entry (frame, width=10, relief='sunken',
textvariable=self.min_cnt)
entry.grid (row=2, column=1, sticky='we')
entry.bind ("<Return>", self.change_contour_limits)
lab = Tkinter.Label (frame, text="Maximum contour:")
lab.grid (row=3, column=0, sticky='w')
entry = Tkinter.Entry (frame, width=10, relief='sunken',
textvariable=self.max_cnt)
entry.grid (row=3, column=1, sticky='we')
entry.bind ("<Return>", self.change_contour_limits)
self.scalar_gui (self.root)
def scalar_gui (self, master):
debug ("In IsoSurfaceProbe::scalar_gui ()")
if not self.scalar_lst:
return
rw = 0
frame = Tkinter.Frame (master, relief='ridge', bd=2)
frame.pack (side='top')
Tkinter.Label (frame, text="Select Scalar").grid (row=rw,
column=0,
sticky='ew')
rw = rw + 1
for sc in self.scalar_lst:
rb0 = Tkinter.Radiobutton (frame, text=sc,
variable=self.scalar_var, value=sc,
command=self.set_scalar_gui)
rb0.grid (row=rw, column=0, sticky='w')
rw = rw + 1
rb1 = Tkinter.Radiobutton (frame, text=sc+" advective flux",
variable=self.scalar_var, value=sc+"Flux",
command=self.set_scalar_gui)
rb1.grid (row=rw, column=0, sticky='w')
entry = Tkinter.Entry (frame, width=10, relief='sunken',
textvariable=self.flux_values[sc+"Flux"])
entry.grid (row=rw, column=1, sticky='we')
entry.bind ("<Return>", self.calc_flux)
rw = rw + 1
rbv = Tkinter.Radiobutton (frame, text="Volume flux",
variable=self.scalar_var, value="VolumeFlux",
command=self.set_scalar_gui)
rbv.grid (row=rw, column=0, sticky='w')
entryv = Tkinter.Entry (frame, width=10, relief='sunken',
textvariable=self.volume_flux_value)
entryv.grid (row=rw, column=1, sticky='we')
entryv.bind ("<Return>", self.calc_flux)
rw = rw + 1
def set_scalar_gui (self):
debug ("In IsoSurfaceProbe::set_scalar_gui ()")
Common.state.busy ()
self.visible_scalar = self.scalar_var.get()
self.Update ()
self.renwin.Render ()
Common.state.idle ()
def GetOutput (self):
debug ("In Projection::GetOutput ()")
return self.iso
def get_array_type(self,arr):
"""Returns if the array is a scalar ('scalars'), vector
('vectors') or tensor ('tensors'). It looks at the number of
components to decide. If it has a wierd number of components it
returns the empty string."""
n = arr.GetNumberOfComponents()
if n == 1:
return 'scalars'
elif n == 3:
return 'vectors'
elif n == 9:
return 'tensors'
else:
return ''
def _get_attribute_list(self, data):
""" Gets scalar, vector and tensor information from the given data
(either cell or point data). """
debug ("In _get_attribute_list ()")
attr = {'scalars':[], 'vectors':[], 'tensors':[]}
if data:
n = data.GetNumberOfArrays()
for i in range(n):
name = data.GetArrayName(i)
type = self.get_array_type(data.GetArray(i))
if type:
attr[type].extend([name])
return attr
def change_contour_limits (self, event=None):
debug ("In IsoSurfaceProbe::change_contour_limits ()")
min_cnt = self.min_cnt.get ()
max_cnt = self.max_cnt.get ()
val = self.contour_var.get ()
dr = self.data_range
if max_cnt < val:
msg = "Error: max. contour value less than current "\
"contour value."
debug (msg)
max_cnt = dr[1]
self.max_cnt.set (max_cnt)
if min_cnt > val:
msg = "Error: min. contour value larger than current "\
"contour value."
debug (msg)
min_cnt = dr[0]
self.min_cnt.set (min_cnt)
resolution = (max_cnt - min_cnt)/1000
if self.slider:
self.slider.config (from_=min_cnt, to=max_cnt,
resolution=resolution)
def change_contour_slider (self, event=None):
debug ("In IsoSurfaceProbe::change_contour_slider ()")
Common.state.busy ()
self.slider_pos = self.contour_var.get ()
self.c_entry_var.set (self.slider_pos)
self.cont_fil.SetValue (0, self.slider_pos)
self.Update ()
Common.state.idle ()
def change_contour_entry (self, event=None):
debug ("In IsoSurfaceProbe::change_contour_entry ()")
Common.state.busy ()
self.slider_pos = self.c_entry_var.get ()
self.contour_var.set (self.slider_pos)
self.cont_fil.SetValue (0, self.slider_pos)
self.Update ()
Common.state.idle ()
def do_sweep (self, event=None):
debug ("In IsoSurfaceProbe::do_sweep ()")
if self.sweep_var.get ():
val = int (1000*self.sweep_delay.get ())
self.root.after (val, self.update_sweep)
def Update (self):
debug ("In IsoSurfaceProbe::Update ()")
Common.state.busy ()
self.cont_fil.Update ()
self.mod_m.Update ()
self.iso.DeepCopy (self.cont_fil.GetOutput ())
self.calc_flux ()
self.iso.GetPointData ().SetActiveScalars (self.visible_scalar)
self.renwin.Render ()
Common.state.idle ()
def update_sweep (self):
debug ("In IsoSurfaceProbe::update_sweep ()")
if self.sweep_var.get ():
min_cnt = self.min_cnt.get ()
max_cnt = self.max_cnt.get ()
d_pos = (max_cnt - min_cnt)/self.sweep_step.get ()
pos = self.slider_pos + d_pos
if (d_pos > 0) and (pos >= max_cnt):
pos = min_cnt
elif (d_pos < 0) and (pos <= min_cnt):
pos = max_cnt
self.contour_var.set (pos)
self.change_contour_slider ()
val = int (1000*self.sweep_delay.get ())
self.root.after (val, self.update_sweep)
def close_gui (self, event=None):
debug ("In IsoSurfaceProbe::close_gui ()")
Base.Objects.Filter.close_gui (self, event)
self.slider = None
def calc_flux (self, event=None):
debug ("In IsoSurfaceProbe::calc_flux ()")
self.cont_fil.Update ()
n_points = self.GetOutput ().GetNumberOfPoints ()
point_normal_calculator = vtk.vtkPolyDataNormals ()
point_normal_calculator.SetInput (self.GetOutput ())
point_normal_calculator.Update ()
normals = point_normal_calculator.GetOutput ().GetPointData ().GetNormals ()
volume_flux_array = vtk.vtkDoubleArray()
volume_flux_array.SetNumberOfTuples(n_points)
volume_flux_array.SetName('VolumeFlux')
vecs = self.GetOutput().GetPointData().GetVectors()
for i in range(n_points):
(u,v,w) = vecs.GetTuple3(i)
(n1,n2,n3) = normals.GetTuple3 (i)
volume_flux = u*n1 + v*n2 + w*n3
volume_flux_array.SetTuple1 (i, volume_flux)
self.GetOutput ().GetPointData ().AddArray(volume_flux_array)
self.GetOutput ().GetPointData ().SetActiveAttribute(volume_flux_array.GetName(),0)
self.GetOutput ().Update ()
flux_array = vtk.vtkDoubleArray()
flux_array.SetNumberOfTuples(n_points)
for sc in self.scalar_lst:
flux_name = sc+"Flux"
flux_array.SetName (flux_name)
scal = self.cont_fil.GetOutput ().GetPointData ().GetScalars (sc)
fluxscal = self.GetOutput ().GetPointData ().GetScalars ("VolumeFlux")
if self.visible_scalar==flux_name:
for i in range(n_points):
s = scal.GetTuple1 (i)
f = fluxscal.GetTuple1 (i)
flux = s*f
flux_array.SetTuple1 (i, flux)
self.GetOutput ().GetPointData ().AddArray(flux_array)
self.GetOutput ().GetPointData ().SetActiveAttribute(flux_array.GetName(),0)
self.GetOutput ().Update ()
# Integration loop
n_cells = self.GetOutput ().GetNumberOfCells ()
Integral_volume_flux = 0.0
scalar_advected_flux = 0.0
flux_name = None
flux_field = None
for sc in self.scalar_lst:
if self.visible_scalar==sc+"Flux":
flux_name = sc+"Flux"
flux_field = sc
break
vecs = self.GetOutput ().GetPointData ().GetVectors ()
if flux_name != None:
fluxf = self.GetOutput().GetPointData().GetScalars(sc)
for cell_no in range(n_cells):
Cell = self.GetOutput().GetCell(cell_no)
if(Cell.GetNumberOfPoints() != 3):
msg = "We have non-triangular elements in surface! You are too interesting. Go away."
Common.print_err (msg)
Cell_points = Cell.GetPoints ()
Area = Cell.TriangleArea (Cell_points.GetPoint(0), \
Cell_points.GetPoint(1), \
Cell_points.GetPoint(2))
n = [0.0,0.0,0.0]
Cell.ComputeNormal (Cell_points.GetPoint(0), \
Cell_points.GetPoint(1), \
Cell_points.GetPoint(2), n)
(u0,v0,w0) = vecs.GetTuple3 (Cell.GetPointId (0))
(u1,v1,w1) = vecs.GetTuple3 (Cell.GetPointId (1))
(u2,v2,w2) = vecs.GetTuple3 (Cell.GetPointId (2))
u = (u0+u1+u2)/3.0
v = (v0+v1+v2)/3.0
w = (w0+w1+w2)/3.0
f = Area*(u*n[0] + v*n[1] + w*n[2])
Integral_volume_flux = Integral_volume_flux + f
if flux_name != None:
s0 = fluxf.GetTuple1 (Cell.GetPointId (0))
s1 = fluxf.GetTuple1 (Cell.GetPointId (1))
s2 = fluxf.GetTuple1 (Cell.GetPointId (2))
s = (s0+s1+s2)/3.0
scalar_advected_flux = scalar_advected_flux + s*f
self.volume_flux_value.set (Integral_volume_flux)
if flux_name != None:
self.flux_values[flux_name].set (scalar_advected_flux)
-------------------------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems? Stop.
Now Search log events and configuration files using AJAX and a browser.
Download your FREE copy of Splunk now >> http://get.splunk.com/
_______________________________________________
MayaVi-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/mayavi-users
