source: trunk/gui/scripts/field.tcl @ 3405

Last change on this file since 3405 was 3405, checked in by gah, 10 years ago

fix to streamlines to (again) display mulitple fields from a single VTK file

File size: 43.6 KB
Line 
1# -*- mode: tcl; indent-tabs-mode: nil -*-
2
3# ----------------------------------------------------------------------
4#  COMPONENT: field - extracts data from an XML description of a field
5#
6#  This object represents one field in an XML description of a device.
7#  It simplifies the process of extracting data vectors that represent
8#  the field.
9# ======================================================================
10#  AUTHOR:  Michael McLennan, Purdue University
11#  Copyright (c) 2004-2012  HUBzero Foundation, LLC
12#
13#  See the file "license.terms" for information on usage and
14#  redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
15# ======================================================================
16package require Itcl
17package require BLT
18
19namespace eval Rappture {
20    # forward declaration
21}
22
23#
24# Possible field dataset types:
25#
26# 2D Datasets
27#       vtk             (range of z-axis is zero).
28#       unirect2d       (deprecated except where extents > 1)
29#       cloud           (x,y point coordinates) (deprecated)
30#       mesh
31# 3D Datasets
32#       vtk
33#       unirect3d
34#       cloud           (x,y,z coordinates) (deprecated)
35#       mesh
36#       dx              (FIXME: make dx-to-vtk converter work)
37#       ucd avs
38#
39# Viewers:
40#       Format     Dim  Description                     Viewer          Server
41#       vtk         2   vtk file data.                  contour         vtkvis
42#       vtk         3   vtk file data.                  isosurface      vtkvis
43#       mesh        2   points-on-mesh                  heightmap       vtkvis
44#       mesh        3   points-on-mesh                  isosurface      vtkvis
45#       dx          3   DX                              volume          nanovis
46#       unirect2d   2   unirect3d + extents > 1 flow    flow            nanovis
47#       unirect3d   3   unirect2d + extents > 1 flow    flow            nanovis
48#       
49# The goal should be any 3D dataset can view a isosurface, volume,
50# streamlines, or flow (if extents > 1).  Any 2D dataset can view a
51# contour, heightmap, streamlines, or flow (if extents > 1). 
52#
53#
54itcl::class Rappture::Field {
55    private variable _dim       0;      # Dimension of the mesh
56    private variable _xmlobj ""  ;      # ref to XML obj with field data
57    private variable _limits     ;      # maps box name => {z0 z1} limits
58    private variable _fldObj ""
59    private variable _comp2fields ;     # cname => field names.
60    private variable _field2Components; # field name => number of components
61    private variable _field2Label;      # field name => label
62    private variable _field2Units;      # field name => units
63    private variable _hints
64    private variable _viewer "";        # Hints which viewer to use
65    constructor {xmlobj path} {
66        # defined below
67    }
68    destructor {
69        # defined below
70    }
71    public method components {args}
72    public method mesh {{cname -overall}}
73    public method values {{cname -overall}}
74    public method blob {{cname -overall}}
75    public method limits {axis}
76    public method controls {option args}
77    public method hints {{key ""}}
78    public method style { cname }
79    public method isunirect2d {}
80    public method isunirect3d {}
81    public method extents {{cname -overall}}
82    public method flowhints { cname }
83    public method type {}
84    public method viewer {} {
85        return $_viewer
86    }
87    public method vtkdata {cname}
88    public method fieldnames { cname } {
89        if { ![info exists _comp2fields($cname)] } {
90            return ""
91        }
92        return $_comp2fields($cname)
93    }
94    public method fieldinfo { fname } {
95        lappend out $_field2Label($fname)
96        lappend out $_field2Units($fname)
97        lappend out $_field2Components($fname)
98        return $out
99    }
100    protected method Build {}
101    protected method _getValue {expr}
102
103    private variable _path "";          # Path of this object in the XML
104    private variable _units ""   ;      # system of units for this field
105    private variable _zmax 0     ;# length of the device
106
107    private variable _comp2dims  ;# maps component name => dimensionality
108    private variable _comp2xy    ;# maps component name => x,y vectors
109    private variable _comp2vtk   ;# maps component name => vtk file data
110    private variable _comp2dx    ;# maps component name => OpenDX data
111    private variable _comp2unirect2d ;# maps component name => unirect2d obj
112    private variable _comp2unirect3d ;# maps component name => unirect3d obj
113    private variable _comp2style ;# maps component name => style settings
114    private variable _comp2cntls ;# maps component name => x,y control points
115    private variable _comp2extents
116    private variable _comp2limits;      #  Array of limits per component
117    private variable _type ""
118    private variable _comp2flowhints
119    private variable _comp2mesh
120    private common _counter 0    ;# counter for unique vector names
121
122    private method BuildPointsOnMesh { cname }
123    private method ConvertToVtkData { cname }
124    private method ReadVtkDataSet { cname contents }
125    private method AvsToVtk { cname contents }
126    private variable _values ""
127}
128
129# ----------------------------------------------------------------------
130# CONSTRUCTOR
131# ----------------------------------------------------------------------
132itcl::body Rappture::Field::constructor {xmlobj path} {
133    if {![Rappture::library isvalid $xmlobj]} {
134        error "bad value \"$xmlobj\": should be Rappture::library"
135    }
136    set _xmlobj $xmlobj
137    set _path $path
138    set _fldObj [$xmlobj element -as object $path]
139    set _units [$_fldObj get units]
140
141    set xunits [$xmlobj get units]
142    if {"" == $xunits || "arbitrary" == $xunits} {
143        set xunits "um"
144    }
145
146    # determine the overall size of the device
147    set z0 [set z1 0]
148    foreach elem [$_xmlobj children components] {
149        switch -glob -- $elem {
150            box* {
151                if {![regexp {[0-9]$} $elem]} {
152                    set elem "${elem}0"
153                }
154                set z0 [$_xmlobj get components.$elem.corner0]
155                set z0 [Rappture::Units::convert $z0 \
156                    -context $xunits -to $xunits -units off]
157
158                set z1 [$_xmlobj get components.$elem.corner1]
159                set z1 [Rappture::Units::convert $z1 \
160                    -context $xunits -to $xunits -units off]
161
162                set _limits($elem) [list $z0 $z1]
163            }
164        }
165    }
166    set _zmax $z1
167
168    # build up vectors for various components of the field
169    Build
170}
171
172# ----------------------------------------------------------------------
173# DESTRUCTOR
174# ----------------------------------------------------------------------
175itcl::body Rappture::Field::destructor {} {
176    itcl::delete object $_fldObj
177    # don't destroy the _xmlobj! we don't own it!
178
179    foreach name [array names _comp2xy] {
180        eval blt::vector destroy $_comp2xy($name)
181    }
182    foreach name [array names _comp2unirect2d] {
183        itcl::delete object $_comp2unirect2d($name)
184    }
185    foreach name [array names _comp2unirect3d] {
186        itcl::delete object $_comp2unirect3d($name)
187    }
188    foreach name [array names _comp2flowhints] {
189        itcl::delete object $_comp2flowhints($name)
190    }
191    foreach name [array names _comp2mesh] {
192        # Data is in the form of a mesh and a vector.
193        foreach { mesh vector } $_comp2mesh($name) break
194        # Release the mesh (may be shared)
195        set class [$mesh info class]
196        ${class}::release $mesh
197        # Destroy the vector
198        blt::vector destroy $vector
199    }
200}
201
202# ----------------------------------------------------------------------
203# USAGE: components ?-name|-dimensions|-style? ?<pattern>?
204#
205# Returns a list of names or types for the various components of
206# this field.  If the optional glob-style <pattern> is specified,
207# then it returns only the components with names matching the pattern.
208# ----------------------------------------------------------------------
209itcl::body Rappture::Field::components {args} {
210    Rappture::getopts args params {
211        flag what -name default
212        flag what -dimensions
213        flag what -style
214        flag what -particles
215        flag what -flow
216        flag what -box
217    }
218
219    set pattern *
220    if {[llength $args] > 0} {
221        set pattern [lindex $args 0]
222        set args [lrange $args 1 end]
223    }
224    if {[llength $args] > 0} {
225        error "wrong # args: should be \"components ?switches? ?pattern?\""
226    }
227
228    set rlist ""
229
230    # BE CAREFUL: return component names in proper order
231    foreach cname [$_fldObj children -type component] {
232        if {[info exists _comp2dims($cname)]
233              && [string match $pattern $cname]} {
234
235            switch -- $params(what) {
236                -name { lappend rlist $cname }
237                -dimensions { lappend rlist $_comp2dims($cname) }
238                -style { lappend rlist $_comp2style($cname) }
239            }
240        }
241    }
242    return $rlist
243}
244
245# ----------------------------------------------------------------------
246# USAGE: mesh ?<name>?
247#
248# Returns a list {xvec yvec} for the specified field component <name>.
249# If the name is not specified, then it returns the vectors for the
250# overall field (sum of all components).
251# ----------------------------------------------------------------------
252itcl::body Rappture::Field::mesh {{cname -overall}} {
253    if {$cname == "-overall" || $cname == "component0"} {
254        set cname [lindex [components -name] 0]
255    }
256    if {[info exists _comp2xy($cname)]} {
257        return [lindex $_comp2xy($cname) 0]  ;# return xv
258    }
259    if { [info exists _comp2vtk($cname)] } {
260        # FIXME: extract mesh from VTK file data.
261        error "method \"mesh\" is not implemented for VTK file data"
262    }
263    if {[info exists _comp2dx($cname)]} {
264        return ""  ;# no mesh -- it's embedded in the value data
265    }
266    if {[info exists _comp2mesh($cname)]} {
267        return ""  ;# no mesh -- it's embedded in the value data
268    }
269    if {[info exists _comp2unirect2d($cname)]} {
270        set mobj [lindex $_comp2unirect2d($cname) 0]
271        return [$mobj mesh]
272    }
273    if {[info exists _comp2unirect3d($cname)]} {
274        set mobj [lindex $_comp2unirect3d($cname) 0]
275        return [$mobj mesh]
276    }
277    error "bad option \"$cname\": should be [join [lsort [array names _comp2dims]] {, }]"
278}
279
280# ----------------------------------------------------------------------
281# USAGE: values ?<name>?
282#
283# Returns a list {xvec yvec} for the specified field component <name>.
284# If the name is not specified, then it returns the vectors for the
285# overall field (sum of all components).
286# ----------------------------------------------------------------------
287itcl::body Rappture::Field::values {{cname -overall}} {
288    if {$cname == "component0"} {
289        set cname "component"
290    }
291    if {[info exists _comp2xy($cname)]} {
292        return [lindex $_comp2xy($cname) 1]  ;# return yv
293    }
294    # VTK file data
295    if { [info exists _comp2vtk($cname)] } {
296        # FIXME: extract the values from the VTK file data
297        error "method \"values\" is not implemented for vtk file data"
298    }
299    # Points-on-mesh
300    if { [info exists _comp2mesh($cname)] } {
301        set vector [lindex $_comp2mesh($cname) 1]
302        return [$vector range 0 end]
303    }
304    if {[info exists _comp2dx($cname)]} {
305        return $_comp2dx($cname)  ;# return gzipped, base64-encoded DX data
306    }
307    if {[info exists _comp2unirect2d($cname)]} {
308        return [$_comp2unirect2d($cname) values]
309    }
310    if {[info exists _comp2unirect3d($cname)]} {
311        return [$_comp2unirect3d($cname) blob]
312    }
313    error "bad option \"$cname\": should be [join [lsort [array names _comp2dims]] {, }]"
314}
315
316# ----------------------------------------------------------------------
317# USAGE: blob ?<name>?
318#
319# Returns a string representing the blob of data for the mesh and values.
320# ----------------------------------------------------------------------
321itcl::body Rappture::Field::blob {{cname -overall}} {
322    if {$cname == "component0"} {
323        set cname "component"
324    }
325    if {[info exists _comp2xy($cname)]} {
326        return ""
327    }
328    if { [info exists _comp2vtk($cname)] } {
329        error "blob not implemented for VTK file data"
330    }
331    if {[info exists _comp2dx($cname)]} {
332        return $_comp2dx($cname)  ;# return gzipped, base64-encoded DX data
333    }
334    if {[info exists _comp2unirect2d($cname)]} {
335        set blob [$_comp2unirect2d($cname) blob]
336        lappend blob "values" $_values
337        return $blob
338    }
339    if {[info exists _comp2unirect3d($cname)]} {
340        return [$_comp2unirect3d($cname) blob]
341    }
342    error "bad option \"$cname\": should be [join [lsort [array names _comp2dims]] {, }]"
343}
344
345# ----------------------------------------------------------------------
346# USAGE: limits <axis>
347#
348# Returns a list {min max} representing the limits for the specified
349# axis.
350# ----------------------------------------------------------------------
351itcl::body Rappture::Field::limits {which} {
352    set min ""
353    set max ""
354    blt::vector tmp zero
355
356    foreach cname [array names _comp2dims] {
357        switch -- $_comp2dims($cname) {
358            1D {
359                switch -- $which {
360                    x - xlin {
361                        set pos 0; set log 0; set axis x
362                    }
363                    xlog {
364                        set pos 0; set log 1; set axis x
365                    }
366                    y - ylin - v - vlin {
367                        set pos 1; set log 0; set axis y
368                    }
369                    ylog - vlog {
370                        set pos 1; set log 1; set axis y
371                    }
372                    default {
373                        error "bad option \"$which\": should be x, xlin, xlog, y, ylin, ylog, v, vlin, vlog"
374                    }
375                }
376
377                set vname [lindex $_comp2xy($cname) $pos]
378                $vname variable vec
379
380                if {$log} {
381                    # on a log scale, use abs value and ignore 0's
382                    $vname dup tmp
383                    $vname dup zero
384                    zero expr {tmp == 0}            ;# find the 0's
385                    tmp expr {abs(tmp)}             ;# get the abs value
386                    tmp expr {tmp + zero*max(tmp)}  ;# replace 0's with abs max
387                    set axisMin [blt::vector expr min(tmp)]
388                    set axisMax [blt::vector expr max(tmp)]
389                } else {
390                    set axisMin $vec(min)
391                    set axisMax $vec(max)
392                }
393
394                if {"" == $min} {
395                    set min $axisMin
396                } elseif {$axisMin < $min} {
397                    set min $axisMin
398                }
399                if {"" == $max} {
400                    set max $axisMax
401                } elseif {$axisMax > $max} {
402                    set max $axisMax
403                }
404            }
405            2D - 3D {
406                if {[info exists _comp2unirect3d($cname)]} {
407                    set limits [$_comp2unirect3d($cname) limits $which]
408                    foreach {axisMin axisMax} $limits break
409                    set axis v
410                } elseif {[info exists _comp2limits($cname)]} {
411                    array set limits $_comp2limits($cname)
412                    switch -- $which {
413                        x - xlin - xlog {
414                            set axis x
415                            foreach {axisMin axisMax} $limits(x) break
416                        }
417                        y - ylin - ylog {
418                            set axis y
419                            foreach {axisMin axisMax} $limits(y) break
420                        }
421                        z - zlin - zlog {
422                            set axis y
423                            foreach {axisMin axisMax} $limits(z) break
424                        }
425                        v - vlin - vlog {
426                            set axis v
427                            foreach {axisMin axisMax} $limits(v) break
428                        }
429                        default {
430                            if { ![info exists limits($which)] } {
431                                error "limits: unknown axis \"$which\""
432                            }
433                            set axis v
434                            foreach {axisMin axisMax} $limits($which) break
435                        }
436                    }
437                } else {
438                    set axisMin 0  ;# HACK ALERT! must be OpenDX data
439                    set axisMax 1
440                    set axis v
441                }
442            }
443        }
444        if { "" == $min || $axisMin < $min } {
445            set min $axisMin
446        }
447        if { "" == $max || $axisMax > $max } {
448            set max $axisMax
449        }
450    }
451    blt::vector destroy tmp zero
452    set val [$_fldObj get "${axis}axis.min"]
453    if {"" != $val && "" != $min} {
454        if {$val > $min} {
455            # tool specified this min -- don't go any lower
456            set min $val
457        }
458    }
459    set val [$_fldObj get "${axis}axis.max"]
460    if {"" != $val && "" != $max} {
461        if {$val < $max} {
462            # tool specified this max -- don't go any higher
463            set max $val
464        }
465    }
466    return [list $min $max]
467}
468
469# ----------------------------------------------------------------------
470# USAGE: controls get ?<name>?
471# USAGE: controls validate <path> <value>
472# USAGE: controls put <path> <value>
473#
474# Returns a list {path1 x1 y1 val1  path2 x2 y2 val2 ...} representing
475# control points for the specified field component <name>.
476# ----------------------------------------------------------------------
477itcl::body Rappture::Field::controls {option args} {
478    switch -- $option {
479        get {
480            set cname [lindex $args 0]
481            if {[info exists _comp2cntls($cname)]} {
482                return $_comp2cntls($cname)
483            }
484            return ""
485        }
486        validate {
487            set path [lindex $args 0]
488            set value [lindex $args 1]
489            set units [$_xmlobj get $path.units]
490
491            if {"" != $units} {
492                set nv [Rappture::Units::convert \
493                    $value -context $units -to $units -units off]
494            } else {
495                set nv $value
496            }
497            if {![string is double $nv]
498                  || [regexp -nocase {^(inf|nan)$} $nv]} {
499                error "Value out of range"
500            }
501
502            set rawmin [$_xmlobj get $path.min]
503            if {"" != $rawmin} {
504                set minv $rawmin
505                if {"" != $units} {
506                    set minv [Rappture::Units::convert \
507                        $minv -context $units -to $units -units off]
508                    set nv [Rappture::Units::convert \
509                        $value -context $units -to $units -units off]
510                }
511                # fix for the case when the user tries to
512                # compare values like minv=-500 nv=-0600
513                set nv [format "%g" $nv]
514                set minv [format "%g" $minv]
515
516                if {$nv < $minv} {
517                    error "Minimum value allowed here is $rawmin"
518                }
519            }
520
521            set rawmax [$_xmlobj get $path.max]
522            if {"" != $rawmax} {
523                set maxv $rawmax
524                if {"" != $units} {
525                    set maxv [Rappture::Units::convert \
526                        $maxv -context $units -to $units -units off]
527                    set nv [Rappture::Units::convert \
528                        $value -context $units -to $units -units off]
529                }
530                # fix for the case when the user tries to
531                # compare values like maxv=-500 nv=-0600
532                set nv [format "%g" $nv]
533                set maxv [format "%g" $maxv]
534
535                if {$nv > $maxv} {
536                    error "Maximum value allowed here is $rawmax"
537                }
538            }
539
540            return "ok"
541        }
542        put {
543            set path [lindex $args 0]
544            set value [lindex $args 1]
545            $_xmlobj put $path.current $value
546            Build
547        }
548        default {
549            error "bad option \"$option\": should be get or put"
550        }
551    }
552}
553
554# ----------------------------------------------------------------------
555# USAGE: hints ?<keyword>?
556#
557# Returns a list of key/value pairs for various hints about plotting
558# this field.  If a particular <keyword> is specified, then it returns
559# the hint for that <keyword>, if it exists.
560# ----------------------------------------------------------------------
561itcl::body Rappture::Field::hints {{keyword ""}} {
562    if { ![info exists _hints] } {
563        foreach {key path} {
564            camera          camera.position
565            color           about.color
566            default         about.default
567            group           about.group
568            label           about.label
569            scale           about.scale
570            seeds           about.seeds
571            style           about.style
572            type            about.type
573            xlabel          about.xaxis.label
574            ylabel          about.yaxis.label
575            zlabel          about.zaxis.label
576            xunits          about.xaxis.units
577            yunits          about.yaxis.units
578            zunits          about.zaxis.units
579            units           units
580            updir           updir
581            vectors         about.vectors
582        } {
583            set str [$_fldObj get $path]
584            if { "" != $str } {
585                set _hints($key) $str
586            }
587        }
588        foreach {key path} {
589            toolid          tool.id
590            toolname        tool.name
591            toolcommand     tool.execute
592            tooltitle       tool.title
593            toolrevision    tool.version.application.revision
594        } {
595            set str [$_xmlobj get $path]
596            if { "" != $str } {
597                set _hints($key) $str
598            }
599        }
600        # Set toolip and path hints
601        set _hints(path) $_path
602        if { [info exists _hints(group)] && [info exists _hints(label)] } {
603            # pop-up help for each curve
604            set _hints(tooltip) $_hints(label)
605        }
606    }
607    if { $keyword != "" } {
608        if {[info exists _hints($keyword)]} {
609            return $_hints($keyword)
610        }
611        return ""
612    }
613    return [array get _hints]
614}
615
616# ----------------------------------------------------------------------
617# USAGE: Build
618#
619# Used internally to build up the vector representation for the
620# field when the object is first constructed, or whenever the field
621# data changes.  Discards any existing vectors and builds everything
622# from scratch.
623# ----------------------------------------------------------------------
624itcl::body Rappture::Field::Build {} {
625
626    # Discard any existing data
627    foreach name [array names _comp2xy] {
628        eval blt::vector destroy $_comp2xy($name)
629    }
630    array unset _comp2vtk
631    foreach name [array names _comp2unirect2d] {
632        eval itcl::delete object $_comp2unirect2d($name)
633    }
634    foreach name [array names _comp2unirect3d] {
635        eval itcl::delete object $_comp2unirect3d($name)
636    }
637    catch {unset _comp2xy}
638    catch {unset _comp2dx}
639    catch {unset _comp2dims}
640    catch {unset _comp2style}
641    array unset _comp2unirect2d
642    array unset _comp2unirect3d
643    array unset _comp2extents
644    array unset _dataobj2type
645    #
646    # Scan through the components of the field and create
647    # vectors for each part.
648    #
649    foreach cname [$_fldObj children -type component] {
650        set type ""
651        if { ([$_fldObj element $cname.constant] != "" &&
652              [$_fldObj element $cname.domain] != "") ||
653              [$_fldObj element $cname.xy] != "" } {
654            set type "1D"
655        } elseif { [$_fldObj element $cname.mesh] != "" &&
656                   [$_fldObj element $cname.values] != ""} {
657            set type "points-on-mesh"
658        } elseif { [$_fldObj element $cname.vtk] != ""} {
659            set viewer [$_fldObj get "about.view"]
660            set type "vtk"
661            if { $viewer != "" } {
662                set _viewer $viewer
663            }
664        } elseif {[$_fldObj element $cname.opendx] != ""} {
665            global env
666            if { [info exists env(VTKVOLUME)] } {
667                set type "vtkvolume"
668            } else {
669                set type "dx"
670            }
671        } elseif {[$_fldObj element $cname.dx] != ""} {
672            global env
673            if { [info exists env(VTKVOLUME)] } {
674                set type "vtkvolume"
675            } else {
676                set type "dx"
677            }
678        }
679        set _comp2style($cname) ""
680       
681        # Save the extents of the component
682        if { [$_fldObj element $cname.extents] != "" } {
683            set extents [$_fldObj get $cname.extents]
684        } else {
685            set extents 1
686        }
687        set _comp2extents($cname) $extents
688        set _type $type
689        if {$type == "1D"} {
690            #
691            # 1D data can be represented as 2 BLT vectors,
692            # one for x and the other for y.
693            #
694            set xv ""
695            set yv ""
696
697            set val [$_fldObj get $cname.constant]
698            if {$val != ""} {
699                set domain [$_fldObj get $cname.domain]
700                if {$domain == "" || ![info exists _limits($domain)]} {
701                    set z0 0
702                    set z1 $_zmax
703                } else {
704                    foreach {z0 z1} $_limits($domain) { break }
705                }
706                set xv [blt::vector create x$_counter]
707                $xv append $z0 $z1
708
709                foreach {val pcomp} [_getValue $val] break
710                set yv [blt::vector create y$_counter]
711                $yv append $val $val
712
713                if {$pcomp != ""} {
714                    set zm [expr {0.5*($z0+$z1)}]
715                    set _comp2cntls($cname) \
716                        [list $pcomp $zm $val "$val$_units"]
717                }
718            } else {
719                set xydata [$_fldObj get $cname.xy]
720                if {"" != $xydata} {
721                    set xv [blt::vector create x$_counter]
722                    set yv [blt::vector create y$_counter]
723                    set tmp [blt::vector create \#auto]
724                    $tmp set $xydata
725                    $tmp split $xv $yv
726                    blt::vector destroy $tmp
727                }
728            }
729
730            if {$xv != "" && $yv != ""} {
731                # sort x-coords in increasing order
732                $xv sort $yv
733                set _comp2dims($cname) "1D"
734                set _comp2xy($cname) [list $xv $yv]
735                incr _counter
736            }
737        } elseif {$type == "points-on-mesh"} {
738            BuildPointsOnMesh $cname
739        } elseif {$type == "vtk"} {
740            set vtkdata [$_fldObj get $cname.vtk]
741            ReadVtkDataSet $cname $vtkdata
742            set _comp2vtk($cname) $vtkdata
743            set _comp2style($cname) [$_fldObj get $cname.style]
744            incr _counter
745        } elseif {$type == "dx" } {
746            #
747            # HACK ALERT!  Extract gzipped, base64-encoded OpenDX
748            # data.  Assume that it's 3D.  Pass it straight
749            # off to the NanoVis visualizer.
750            #
751            set _viewer "nanovis"
752            set _comp2dims($cname) "3D"
753            set _comp2dx($cname)  [$_fldObj get -decode no $cname.dx]
754            if 1 {
755            set data  [$_fldObj get -decode yes $cname.dx]
756            set file "/tmp/junk.dx"
757            set f [open $file "w"]
758            puts $f $data
759            close $f
760            if { [string match "<ODX>*" $data] } {
761                set data [string range $data 5 end]
762                set _comp2dx($cname) \
763                        [Rappture::encoding::encode -as zb64 $data]
764            }
765            }
766            set _comp2style($cname) [$_fldObj get $cname.style]
767            if {[$_fldObj element $cname.flow] != ""} {
768                set _comp2flowhints($cname) \
769                    [Rappture::FlowHints ::\#auto $_fldObj $cname $_units]
770            }
771            incr _counter
772        } elseif {$type == "opendx"} {
773            #
774            # HACK ALERT!  Extract gzipped, base64-encoded OpenDX
775            # data.  Assume that it's 3D.  Pass it straight
776            # off to the NanoVis visualizer.
777            #
778            set _viewer "nanovis"
779            set _comp2dims($cname) "3D"
780            set data [$_fldObj get -decode yes $cname.opendx]
781            set data "<ODX>$data"
782            set data [Rappture::encoding::encode -as zb64 $data]
783            set _comp2dx($cname) $data
784            set _comp2style($cname) [$_fldObj get $cname.style]
785            if {[$_fldObj element $cname.flow] != ""} {
786                set _comp2flowhints($cname) \
787                    [Rapture::FlowHints ::\#auto $_fldObj $cname $_units]
788            }
789            incr _counter
790        } elseif {[$_fldObj element $cname.ucd] != ""} {
791            set _viewer "isosurface"
792            set _comp2dims($cname) "3D"
793            set contents [$_fldObj get $cname.ucd]
794            set vtkdata [AvsToVtk $cname $contents]
795            ReadVtkDataSet $cname $vtkdata
796            set _comp2vtk($cname) $vtkdata
797            set _comp2style($cname) [$_fldObj get $cname.style]
798            incr _counter
799        }
800    }
801    # Sanity check.  Verify that all components of the field have the same
802    # dimension.
803    set dim ""
804    foreach cname [array names _comp2dims] {
805        if { $dim == "" } {
806            set dim $_comp2dims($cname)
807            continue
808        }
809        if { $dim != $_comp2dims($cname) } {
810            error "field can't have components of different dimensions: [join [array get _comp2dims] ,]"
811        }
812    }
813    # FIXME: about.scalars and about.vectors are temporary.  With views
814    #        the user will specify the label, units or each field there.
815    #
816    # Override what we found in the VTK file with names that the user
817    # selected.  We override the field label and units.
818    foreach { fname label units } [$_fldObj get about.scalars] {
819        if { ![info exists _field2Name($fname)] } {
820            set _field2Name($fname) $fname
821            set _field2Components($fname) 1
822        }
823        set _field2Label($fname) $label
824        set _field2Units($fname) $units
825    }
826    foreach { fname label units } [$_fldObj get about.vectors] {
827        if { ![info exists _field2Name($fname)] } {
828            set _field2Name($fname) $fname
829            set _field2Components($fname) 3
830        }
831        set _field2Label($fname) $label
832        set _field2Units($fname) $units
833    }
834}
835
836# ----------------------------------------------------------------------
837# USAGE: _getValue <expr>
838#
839# Used internally to get the value for an expression <expr>.  Returns
840# a list of the form {val parameterPath}, where val is the numeric
841# value of the expression, and parameterPath is the XML path to the
842# parameter representing the value, or "" if the <expr> does not
843# depend on any parameters.
844# ----------------------------------------------------------------------
845itcl::body Rappture::Field::_getValue {expr} {
846    #
847    # First, look for the expression among the <parameter>'s
848    # associated with the device.
849    #
850    set found 0
851    foreach pcomp [$_xmlobj children parameters] {
852        set id [$_xmlobj element -as id parameters.$pcomp]
853        if {[string equal $id $expr]} {
854            set val [$_xmlobj get parameters.$pcomp.current]
855            if {"" == $val} {
856                set val [$_xmlobj get parameters.$pcomp.default]
857            }
858            if {"" != $val} {
859                set expr $val
860                set found 1
861                break
862            }
863        }
864    }
865    if {$found} {
866        set pcomp "parameters.$pcomp"
867    } else {
868        set pcomp ""
869    }
870
871    if {$_units != ""} {
872        set expr [Rappture::Units::convert $expr \
873            -context $_units -to $_units -units off]
874    }
875
876    return [list $expr $pcomp]
877}
878
879#
880# isunirect2d  --
881#
882# Returns if the field is a unirect2d object. 
883#
884itcl::body Rappture::Field::isunirect2d { } {
885    return [expr [array size _comp2unirect2d] > 0]
886}
887
888#
889# isunirect3d  --
890#
891# Returns if the field is a unirect3d object. 
892#
893itcl::body Rappture::Field::isunirect3d { } {
894    return [expr [array size _comp2unirect3d] > 0]
895}
896
897#
898# flowhints  --
899#
900# Returns the hints associated with a flow vector field. 
901#
902itcl::body Rappture::Field::flowhints { cname } {
903    if { [info exists _comp2flowhints($cname)] } {
904        return $_comp2flowhints($cname)
905    }
906    return ""
907}
908
909#
910# style  --
911#
912# Returns the style associated with a component of the field. 
913#
914itcl::body Rappture::Field::style { cname } {
915    if { [info exists _comp2style($cname)] } {
916        return $_comp2style($cname)
917    }
918    return ""
919}
920
921#
922# type  --
923#
924# Returns the style associated with a component of the field. 
925#
926itcl::body Rappture::Field::type {} {
927    return $_type
928}
929
930#
931# extents --
932#
933# Returns if the field is a unirect2d object. 
934#
935itcl::body Rappture::Field::extents {{cname -overall}} {
936    if {$cname == "-overall" } {
937        set max 0
938        foreach cname [$_fldObj children -type component] {
939            if { ![info exists _comp2unirect3d($cname)] &&
940                 ![info exists _comp2extents($cname)] } {
941                continue
942            }
943            set value $_comp2extents($cname)
944            if { $max < $value } {
945                set max $value
946            }
947        }
948        return $max
949    }
950    if { $cname == "component0"} {
951        set cname [lindex [components -name] 0]
952    }
953    return $_comp2extents($cname)
954}
955
956itcl::body Rappture::Field::ConvertToVtkData { cname } {
957    set ds ""
958    switch -- [typeof $cname] {
959        "unirect2d" {
960            foreach { x1 x2 xN y1 y2 yN } [$dataobj mesh $cname] break
961            set spacingX [expr {double($x2 - $x1)/double($xN - 1)}]
962            set spacingY [expr {double($y2 - $y1)/double($yN - 1)}]
963           
964            set ds [vtkImageData $this-grdataTemp]
965            $ds SetDimensions $xN $yN 1
966            $ds SetOrigin $x1 $y1 0
967            $ds SetSpacing $spacingX $spacingY 0
968            set arr [vtkDoubleArray $this-arrTemp]
969            foreach {val} [$dataobj values $cname] {
970                $arr InsertNextValue $val
971            }
972            [$ds GetPointData] SetScalars $arr
973        }
974        "unirect3d" {
975            foreach { x1 x2 xN y1 y2 yN z1 z2 zN } [$dataobj mesh $cname] break
976            set spacingX [expr {double($x2 - $x1)/double($xN - 1)}]
977            set spacingY [expr {double($y2 - $y1)/double($yN - 1)}]
978            set spacingZ [expr {double($z2 - $z1)/double($zN - 1)}]
979           
980            set ds [vtkImageData $this-grdataTemp]
981            $ds SetDimensions $xN $yN $zN
982            $ds SetOrigin $x1 $y1 $z1
983            $ds SetSpacing $spacingX $spacingY $spacingZ
984            set arr [vtkDoubleArray $this-arrTemp]
985            foreach {val} [$dataobj values $cname] {
986                $arr InsertNextValue $val
987            }
988            [$ds GetPointData] SetScalars $val
989        }
990        "contour" {
991            return [$dataobj blob $cname]
992        }
993        "dx" {
994            return [Rappture::ConvertDxToVtk $_comp2dx($cname)]
995        }
996        default {
997            set mesh [$dataobj mesh $cname]
998            switch -- [$mesh GetClassName] {
999                vtkPoints {
1000                    # handle cloud of points
1001                    set ds [vtkPolyData $this-polydataTemp]
1002                    $ds SetPoints $mesh
1003                    [$ds GetPointData] SetScalars [$dataobj values $cname]
1004                }
1005                vtkPolyData {
1006                    set ds [vtkPolyData $this-polydataTemp]
1007                    $ds ShallowCopy $mesh
1008                    [$ds GetPointData] SetScalars [$dataobj values $cname]
1009                }
1010                vtkUnstructuredGrid {
1011                    # handle 3D grid with connectivity
1012                    set ds [vtkUnstructuredGrid $this-grdataTemp]
1013                    $ds ShallowCopy $mesh
1014                    [$ds GetPointData] SetScalars [$dataobj values $cname]
1015                }
1016                vtkRectilinearGrid {
1017                    # handle 3D grid with connectivity
1018                    set ds [vtkRectilinearGrid $this-grdataTemp]
1019                    $ds ShallowCopy $mesh
1020                    [$ds GetPointData] SetScalars [$dataobj values $cname]
1021                }
1022                default {
1023                    error "don't know how to handle [$mesh GetClassName] data"
1024                }
1025            }
1026        }
1027    }
1028
1029    if {"" != $ds} {
1030        set writer [vtkDataSetWriter $this-dsWriterTmp]
1031        $writer SetInput $ds
1032        $writer SetFileTypeToASCII
1033        $writer WriteToOutputStringOn
1034        $writer Write
1035        set out [$writer GetOutputString]
1036        $ds Delete
1037        $writer Delete
1038    } else {
1039        set out ""
1040        error "No DataSet to write"
1041    }
1042
1043    append out "\n"
1044    return $out
1045}
1046
1047itcl::body Rappture::Field::ReadVtkDataSet { cname contents } {
1048    package require vtk
1049
1050    set reader $this-datasetreader
1051    vtkDataSetReader $reader
1052
1053    # Write the contents to a file just in case it's binary.
1054    set tmpfile file[pid].vtk
1055    set f [open "$tmpfile" "w"]
1056    fconfigure $f -translation binary -encoding binary
1057    puts $f $contents
1058    close $f
1059    $reader SetFileName $tmpfile
1060    $reader ReadAllScalarsOn
1061    $reader ReadAllVectorsOn
1062    $reader ReadAllFieldsOn
1063    $reader Update
1064    set dataset [$reader GetOutput]
1065    set limits {}
1066    foreach {xmin xmax ymin ymax zmin zmax} [$dataset GetBounds] break
1067    # Figure out the dimension of the mesh from the bounds.
1068    set _dim 0
1069    if { $xmax > $xmin } {
1070        incr _dim
1071    }
1072    if { $ymax > $ymin } {
1073        incr _dim
1074    }
1075    if { $zmax > $zmin } {
1076        incr _dim
1077    }
1078    if { $_viewer == "" } {
1079        if { $_dim == 2 } {
1080            set _viewer contour
1081        } else {
1082            set _viewer isosurface
1083        }
1084    }
1085    set _comp2dims($cname) ${_dim}D
1086    lappend limits x [list $xmin $xmax]
1087    lappend limits y [list $ymin $ymax]
1088    lappend limits z [list $zmin $zmax]
1089    set dataAttrs [$dataset GetPointData]
1090    if { $dataAttrs == ""} {
1091        puts stderr "No point data"
1092    }
1093    set vmin 0
1094    set vmax 1
1095    set numArrays [$dataAttrs GetNumberOfArrays]
1096    if { $numArrays > 0 } {
1097        set array [$dataAttrs GetArray 0]
1098        foreach {vmin vmax} [$array GetRange] break
1099
1100        for {set i 0} {$i < [$dataAttrs GetNumberOfArrays] } {incr i} {
1101            set array [$dataAttrs GetArray $i]
1102            set fname  [$dataAttrs GetArrayName $i]
1103            foreach {min max} [$array GetRange] break
1104            lappend limits $fname [list $min $max]
1105            set _field2Units($fname) ""
1106            set _field2Label($fname) $fname
1107            set _field2Components($fname) [$array GetNumberOfComponents]
1108            lappend _comp2fields($cname) $fname
1109        }
1110    }
1111    lappend limits v [list $vmin $vmax]
1112    set _comp2limits($cname) $limits
1113    $reader Delete
1114    file delete $tmpfile
1115}
1116
1117#
1118# vtkdata --
1119#
1120#       Returns a string representing the mesh and field data for a specific
1121#       component in the legacy VTK file format.
1122#
1123itcl::body Rappture::Field::vtkdata {cname} {
1124    if {$cname == "component0"} {
1125        set cname "component"
1126    }
1127    # DX: Convert DX to VTK
1128    if {[info exists _comp2dx($cname)]} {
1129        return [Rappture::ConvertDxToVtk $_comp2dx($cname)]
1130    }
1131    # Unirect3d: isosurface
1132    if {[info exists _comp2unirect3d($cname)]} {
1133        return [$_comp2unirect3d($cname) vtkdata]
1134    }
1135    # VTK file data:
1136    if { [info exists _comp2vtk($cname)] } {
1137        return $_comp2vtk($cname)
1138    }
1139    # Points on mesh:  Construct VTK file output.
1140    if { [info exists _comp2mesh($cname)] } {
1141        # Data is in the form mesh and vector
1142        foreach {mesh vector} $_comp2mesh($cname) break
1143        set label [hints zlabel]
1144        if { $label == "" } {
1145            set label $cname
1146        } else {
1147            regsub -all { } $label {_} label
1148        }
1149        append out "# vtk DataFile Version 3.0\n"
1150        append out "[hints label]\n"
1151        append out "ASCII\n"
1152        append out [$mesh vtkdata]
1153        append out "POINT_DATA [$vector length]\n"
1154        append out "SCALARS $label float\n"
1155        append out "LOOKUP_TABLE default\n"
1156        append out "[$vector range 0 end]\n"
1157        return $out
1158    }
1159    error "can't find vtkdata for $cname. This method should only be called by the vtkheightmap widget"
1160}
1161
1162#
1163# BuildPointsOnMesh --
1164#
1165#       Parses the field XML description to build a mesh and values vector
1166#       representing the field.  Right now we handle the deprecated types
1167#       of "cloud", "unirect2d", and "unirect3d" (mostly for flows).
1168#
1169itcl::body Rappture::Field::BuildPointsOnMesh {cname} {
1170    #
1171    # More complex 2D/3D data is represented by a mesh
1172    # object and an associated vector for field values.
1173    #
1174    set path [$_fldObj get $cname.mesh]
1175    if {[$_xmlobj element $path] == ""} {
1176        # Unknown mesh designated.
1177        return
1178    }
1179    set element [$_xmlobj element -as type $path]
1180    set name $cname
1181    regsub -all { } $name {_} name
1182    set _field2Label($name) $cname
1183    set label [hints zlabel]
1184    if { $label != "" } {
1185        set _field2Label($name) $label
1186    }
1187    set _field2Units($name) [hints zlabel]
1188    set _field2Components($name) 1
1189    lappend _comp2fields($cname) $name
1190
1191    # Handle bizarre cases that hopefully will be deprecated.
1192    if { $element == "unirect3d" } {
1193        # Special case: unirect3d (should be deprecated) + flow.
1194        if { [$_fldObj element $cname.extents] != "" } {
1195            set extents [$_fldObj get $cname.extents]
1196        } else {
1197            set extents 1
1198        }
1199        set _dim 3
1200        set _viewer flowvis
1201        set _comp2dims($cname) "3D"
1202        set _comp2unirect3d($cname) \
1203            [Rappture::Unirect3d \#auto $_xmlobj $_fldObj $cname $extents]
1204        set _comp2style($cname) [$_fldObj get $cname.style]
1205        if {[$_fldObj element $cname.flow] != ""} {
1206            set _comp2flowhints($cname) \
1207                [Rappture::FlowHints ::\#auto $_fldObj $cname $_units]
1208        }
1209        incr _counter
1210        return
1211    }
1212    if { $element == "unirect2d" && [$_fldObj element $cname.flow] != "" } {
1213        # Special case: unirect2d (normally deprecated) + flow.
1214        if { [$_fldObj element $cname.extents] != "" } {
1215            set extents [$_fldObj get $cname.extents]
1216        } else {
1217            set extents 1
1218        }
1219        set _dim 2
1220        set _viewer "flowvis"
1221        set _comp2dims($cname) "2D"
1222        set _comp2unirect2d($cname) \
1223            [Rappture::Unirect2d \#auto $_xmlobj $path]
1224        set _comp2style($cname) [$_fldObj get $cname.style]
1225        set _comp2flowhints($cname) \
1226            [Rappture::FlowHints ::\#auto $_fldObj $cname $_units]
1227        set _values [$_fldObj get $cname.values]
1228        set limits {}
1229        foreach axis { x y } {
1230            lappend limits $axis [$_comp2unirect2d($cname) limits $axis]
1231        }
1232        set xv [blt::vector create \#auto]
1233        $xv set $_values
1234        lappend limits "v" [$xv limits]
1235        blt::vector destroy $xv
1236        set _comp2limits($cname) $limits
1237        incr _counter
1238        return
1239    }
1240    set _viewer "contour"
1241    switch -- $element {
1242        "cloud" {
1243            set mesh [Rappture::Cloud::fetch $_xmlobj $path]
1244        }
1245        "mesh" {
1246            set mesh [Rappture::Mesh::fetch $_xmlobj $path]
1247        }           
1248        "unirect2d" {
1249            set mesh [Rappture::Unirect2d::fetch $_xmlobj $path]
1250            set _viewer "heightmap"
1251        }
1252    }
1253    set _dim [$mesh dimensions]
1254    if {$_dim == 1} {
1255        # Is this used anywhere?
1256        #
1257        # OOPS!  This is 1D data
1258        # Forget the cloud/field -- store BLT vectors
1259        #
1260        # Is there a natural growth path in generating output from 1D to
1261        # higher dimensions?  If there isn't, let's kill this in favor
1262        # or explicitly using a <curve> instead.  Otherwise, the features
1263        # (methods such as xmarkers) or the <curve> need to be added
1264        # to the <field>.
1265        #
1266        set xv [blt::vector create x$_counter]
1267        set yv [blt::vector create y$_counter]
1268       
1269        $yv set [$mesh points]
1270        $xv seq 0 1 [$yv length]
1271        # sort x-coords in increasing order
1272        $xv sort $yv
1273       
1274        set _comp2dims($cname) "1D"
1275        set _comp2xy($cname) [list $xv $yv]
1276        incr _counter
1277        return
1278    } elseif {$_dim == 2} {
1279        set _type "heightmap"
1280        set vector [blt::vector create \#auto]
1281        $vector set [$_fldObj get $cname.values]
1282        set _comp2dims($cname) "[$mesh dimensions]D"
1283        set _comp2mesh($cname) [list $mesh $vector]
1284        set _comp2style($cname) [$_fldObj get $cname.style]
1285        incr _counter
1286        array unset _comp2limits $cname
1287        lappend _comp2limits($cname) x [$mesh limits x]
1288        lappend _comp2limits($cname) y [$mesh limits y]
1289        lappend _comp2limits($cname) $label [$vector limits]
1290        lappend _comp2limits($cname) v [$vector limits]
1291        return
1292    } elseif {$_dim == 3} {
1293        #
1294        # 3D data: Store cloud/field as components
1295        #
1296        set values [$_fldObj get $cname.values]
1297        set farray [vtkFloatArray ::vals$_counter]
1298        foreach v $values {
1299            if {"" != $_units} {
1300                set v [Rappture::Units::convert $v \
1301                           -context $_units -to $_units -units off]
1302            }
1303            $farray InsertNextValue $v
1304        }
1305        set _viewer "isosurface"
1306        set _type "isosurface"
1307        set vector [blt::vector create \#auto]
1308        $vector set [$_fldObj get $cname.values]
1309        set _comp2dims($cname) "[$mesh dimensions]D"
1310        set _comp2mesh($cname) [list $mesh $vector]
1311        set _comp2style($cname) [$_fldObj get $cname.style]
1312        incr _counter
1313        lappend _comp2limits($cname) x [$mesh limits x]
1314        lappend _comp2limits($cname) y [$mesh limits y]
1315        lappend _comp2limits($cname) z [$mesh limits z]
1316        lappend _comp2limits($cname) $label [$vector limits]
1317        lappend _comp2limits($cname) v [$vector limits]
1318        return
1319    }
1320    error "unhandled case in field dim=$_dim element=$element"
1321}
1322
1323itcl::body Rappture::Field::AvsToVtk { comp contents } {
1324    package require vtk
1325
1326    set reader $this-datasetreader
1327    vtkAVSucdReader $reader
1328
1329    # Write the contents to a file just in case it's binary.
1330    set tmpfile file[pid].vtk
1331    set f [open "$tmpfile" "w"]
1332    fconfigure $f -translation binary -encoding binary
1333    puts $f $contents
1334    close $f
1335    $reader SetFileName $tmpfile
1336    $reader Update
1337    file delete $tmpfile
1338
1339    set output [$reader GetOutput]
1340    set pointData [$output GetPointData]
1341    set _scalars {}
1342    for { set i 0 } { $i < [$pointData GetNumberOfArrays] } { incr i } {
1343        set name [$pointData GetArrayName $i]
1344        lappend _scalars $name $name "???"
1345    }
1346    set writer $this-datasetwriter
1347    vtkDataSetWriter $writer
1348    $writer SetInputConnection [$reader GetOutputPort]
1349    $writer SetFileName $tmpfile
1350    $writer Write
1351    rename $reader ""
1352    rename $writer ""
1353    set f [open "$tmpfile" "r"]
1354    fconfigure $f -translation binary -encoding binary
1355    set vtkdata [read $f]
1356    close $f
1357    file delete $tmpfile
1358    return $vtkdata
1359}
Note: See TracBrowser for help on using the repository browser.