source: branches/Rappture 1.2/gui/scripts/mesh.tcl @ 3288

Last change on this file since 3288 was 3288, checked in by gah, 12 years ago

make placeholder for isosurface viewer
File size: 34.1 KB
Line 
1
2# ----------------------------------------------------------------------
3#  COMPONENT: mesh - represents a structured mesh for a device
4#
5#  This object represents a mesh in an XML description of simulator
6#  output.  A mesh is a structured arrangement of points, as elements
7#  composed of nodes representing coordinates.  This is a little
8#  different from a cloud, which is an unstructured arrangement
9#  (shotgun blast) of points.
10# ======================================================================
11#  AUTHOR:  Michael McLennan, Purdue University
12#  Copyright (c) 2004-2012  HUBzero Foundation, LLC
13#
14#  See the file "license.terms" for information on usage and
15#  redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
16# ======================================================================
17package require Itcl
18
19namespace eval Rappture {
20    # forward declaration
21}
22
23itcl::class Rappture::Mesh {
24    private variable _xmlobj ""  ;      # Ref to XML obj with device data
25    private variable _mesh ""    ;      # Lib obj representing this mesh
26    private variable _dim       0;      # Dimension of mesh (1, 2, or 3)
27    private variable _type "";          # Indicates the type of mesh.
28    private variable _units "m m m" ;   # System of units for x, y, z
29    private variable _limits        ;   # Array of mesh limits. Keys are
30                                        # xmin, xmax, ymin, ymax, ...
31    private variable _numPoints 0   ;   # # of points in mesh
32    private variable _vtkdata "";       # Mesh in vtk file format.
33    constructor {xmlobj path} {
34        # defined below
35    }
36    destructor {
37        # defined below
38    }
39    public method points {}
40    public method elements {}
41    public method mesh {{-type "vtk"}}
42    public method size {{what -points}}
43    public method dimensions {}
44    public method limits {which}
45    public method hints {{key ""}}
46
47    public proc fetch {xmlobj path}
48    public proc release {obj}
49    public method vtkdata {}
50    public method type {} {
51        return $_type
52    }
53    public method numpoints {} {
54        return $_numPoints
55    }
56
57    private method ReadNodesElements {xmlobj path}
58    private method _buildRectMesh {xmlobj path}
59    private method _getVtkElement {npts}
60
61    private variable _pts2elem   ;      # maps # points => vtk element
62
63    private variable _points ""      ;  # name of vtkPointData object
64
65    private common _xp2obj       ;      # used for fetch/release ref counting
66    private common _obj2ref      ;      # used for fetch/release ref counting
67    private variable _isValid 0
68    private variable _vtkoutput ""
69    private variable _vtkreader ""
70    private variable _vtkpoints ""
71    private variable _vtkgrid   ""
72    private variable _xv        ""
73    private variable _yv        ""
74    private variable _zv        ""
75    private variable _bltmesh   ""
76    private variable _xValues   "";     # For the blt contour only
77    private variable _yValues   "";     # For the blt contour only
78   
79    private method GetDimension {}
80    private method GetDouble { path }
81    private method GetInt { path }
82    private method ReadBltCloudMesh { path }
83    private method ReadBltGridMesh { path }
84    private method ReadBltMesh { path }
85    private method ReadBltTriangleMesh { path }
86    private method ReadCells { xmlobj path }
87    private method ReadCloud { xmlobj path }
88    private method ReadGrid { xmlobj path }
89    private method ReadTriangles { xmlobj path }
90    private method ReadVtk { xmlobj path }
91    private method _buildNodesElements {xmlobj path}
92}
93
94# ----------------------------------------------------------------------
95# USAGE: Rappture::Mesh::fetch <xmlobj> <path>
96#
97# Clients use this instead of a constructor to fetch the Mesh for
98# a particular <path> in the <xmlobj>.  When the client is done with
99# the mesh, he calls "release" to decrement the reference count.
100# When the mesh is no longer needed, it is cleaned up automatically.
101# ----------------------------------------------------------------------
102itcl::body Rappture::Mesh::fetch {xmlobj path} {
103    set handle "$xmlobj|$path"
104    if {[info exists _xp2obj($handle)]} {
105        set obj $_xp2obj($handle)
106        incr _obj2ref($obj)
107        return $obj
108    }
109    set obj [Rappture::Mesh ::#auto $xmlobj $path]
110    set _xp2obj($handle) $obj
111    set _obj2ref($obj) 1
112    return $obj
113}
114
115# ----------------------------------------------------------------------
116# USAGE: Rappture::Mesh::release <obj>
117#
118# Clients call this when they're no longer using a Mesh fetched
119# previously by the "fetch" proc.  This decrements the reference
120# count for the mesh and destroys the object when it is no longer
121# in use.
122# ----------------------------------------------------------------------
123itcl::body Rappture::Mesh::release {obj} {
124    if {[info exists _obj2ref($obj)]} {
125        incr _obj2ref($obj) -1
126        if {$_obj2ref($obj) <= 0} {
127            unset _obj2ref($obj)
128            foreach handle [array names _xp2obj] {
129                if {$_xp2obj($handle) == $obj} {
130                    unset _xp2obj($handle)
131                }
132            }
133            itcl::delete object $obj
134        }
135    } else {
136        error "can't find reference count for $obj"
137    }
138}
139
140# ----------------------------------------------------------------------
141# CONSTRUCTOR
142# ----------------------------------------------------------------------
143itcl::body Rappture::Mesh::constructor {xmlobj path} {
144    package require vtk
145    if {![Rappture::library isvalid $xmlobj]} {
146        error "bad value \"$xmlobj\": should be Rappture::library"
147    }
148    set _xmlobj $xmlobj
149    set _mesh [$xmlobj element -as object $path]
150
151    # Initialize mesh bounds to empty
152    foreach lim {xmin xmax ymin ymax zmin zmax} {
153        set _limits($lim) ""
154    }
155    GetDimension
156    set u [$_mesh get units]
157    if {"" != $u} {
158        while {[llength $u] < 3} {
159            lappend u [lindex $u end]
160        }
161        set _units $u
162    }
163
164    foreach lim {xmin xmax ymin ymax zmin zmax} {
165        set _limits($lim) ""
166    }
167    # Meshes comes in a variety of flavors
168    #
169    # Dimensionality is determined from the <dimension> tag. 
170    #
171    # <vtk> described mesh
172    # <element> +  <node> definitions
173    # <cloud>           x,y coordinates or x,y,z coordinates
174    # <grid>            rectangular mesh
175    # <triangles>       triangular mesh
176    # <cells>           homogeneous cell type mesh.
177
178    # Check that only one mesh type was defined.
179    set subcount 0
180    foreach cname [$_mesh children] {
181        foreach type { vtk cloud grid triangles cells } {
182            if { [$_mesh element $type] != "" } {
183                incr subcount
184                break
185            }
186        }
187    }
188    set elemcount 0
189    foreach cname [$_mesh children] {
190        foreach type { node element } {
191            if { [$_mesh element $type] != "" } {
192                incr elemcount
193                break
194            }
195        }
196    }
197    if { $elemcount > 0 } {
198        incr $subcount
199    }
200    if { $subcount > 1 } {
201        puts stderr "too many mesh types specified: picking first found."
202    }
203    if { [$_mesh element vtk] != ""} {
204        ReadVtk $xmlobj $path
205    } elseif { [$_mesh element cloud] != "" } {
206        ReadCloud $xmlobj $path
207    } elseif {[$_mesh element grid] != "" } {
208        ReadGrid $xmlobj $path
209    } elseif {[$_mesh element triangles] != "" } {
210        ReadTriangles $xmlobj $path
211    } elseif {[$_mesh element cells] != "" } {
212        ReadCells $xmlobj $path
213    } elseif {[$_mesh element node] != "" && [$_mesh element element] != ""} {
214        ReadNodesElements $xmlobj $path
215    }
216}
217
218# ----------------------------------------------------------------------
219# DESTRUCTOR
220# ----------------------------------------------------------------------
221itcl::body Rappture::Mesh::destructor {} {
222    # don't destroy the _xmlobj! we don't own it!
223    itcl::delete object $_mesh
224
225    if { $_xValues != "" } {
226        blt::vector destroy $_xValues
227    }
228    if { $_yValues != "" } {
229        blt::vector destroy $_yValues
230    }
231    if { $_vtkoutput != "" } {
232        rename $_vtkoutput ""
233    }
234    if { $_vtkreader != "" } {
235        rename $_vtkreader ""
236    }
237    if { $_points != "" } {
238        rename $_points ""
239    }
240    if { $_bltmesh != "" } {
241        blt::mesh destroy $_bltmesh
242    }
243    catch {rename $this-points ""}
244    catch {rename $this-grid ""}
245    catch {rename $this-dset ""}
246    foreach type [array names _pts2elem] {
247        rename $_pts2elem($type) ""
248    }
249}
250
251#
252# vtkdata --
253#
254#       This is called by the field object to generate a VTK file to send to
255#       the remote render server.  Returns the vtkDataSet object containing
256#       (at this point) just the mesh.  The field object doesn't know (or
257#       care) what type of mesh is used.  The field object will add field
258#       arrays before generating output to send to the remote render server.
259#
260itcl::body Rappture::Mesh::vtkdata {} {
261    return $_vtkdata
262}
263
264# ----------------------------------------------------------------------
265# USAGE: points
266#
267# Returns the vtk object containing the points for this mesh.
268# ----------------------------------------------------------------------
269itcl::body Rappture::Mesh::points {} {
270    return $_points
271}
272
273# ----------------------------------------------------------------------
274# USAGE: elements
275#
276# Returns a list of the form {plist r plist r ...} for each element
277# in this mesh.  Each plist is a list of {x y x y ...} coordinates
278# for the mesh.
279# ----------------------------------------------------------------------
280itcl::body Rappture::Mesh::elements {} {
281    # build a map for region number => region type
282    foreach comp [$_mesh children -type region] {
283        set id [$_mesh element -as id $comp]
284        set regions($id) [$_mesh get $comp]
285    }
286    set regions() "unknown"
287
288    set rlist ""
289    foreach comp [$_mesh children -type element] {
290        set rid [$_mesh get $comp.region]
291
292        #
293        # HACK ALERT!
294        #
295        # Prophet puts out nodes in a funny "Z" shaped order,
296        # not in proper clockwise fashion.  Switch the last
297        # two nodes for now to make them correct.
298        #
299        set nlist [$_mesh get $comp.nodes]
300        set n2 [lindex $nlist 2]
301        set n3 [lindex $nlist 3]
302        set nlist [lreplace $nlist 2 3 $n3 $n2]
303        lappend nlist [lindex $nlist 0]
304
305        set plist ""
306        foreach nid $nlist {
307            eval lappend plist [$_mesh get node($nid)]
308        }
309        lappend rlist $plist $regions($rid)
310    }
311    return $rlist
312}
313
314# ----------------------------------------------------------------------
315# USAGE: mesh
316#
317# Returns the vtk object representing the mesh.
318# ----------------------------------------------------------------------
319itcl::body Rappture::Mesh::mesh { {type "vtk"} } {
320    switch $type {
321        "vtk" {
322            return $_vtkgrid
323        }
324        "blt" {
325            if { $_bltmesh == "" } {
326                ReadBltMesh
327                return $_bltmesh
328            }
329        }
330        default {
331            error "Requested mesh type \"$type\" is unknown."
332        }
333    }
334}
335
336# ----------------------------------------------------------------------
337# USAGE: size ?-points|-elements?
338#
339# Returns the number of points in this mesh.
340# ----------------------------------------------------------------------
341itcl::body Rappture::Mesh::size {{what -points}} {
342    switch -- $what {
343        -points {
344            return [$_vtkpoints GetNumberOfPoints]
345        }
346        -elements {
347            return [$_vtkgrid GetNumberOfCells]
348        }
349        default {
350            error "bad option \"$what\": should be -points or -elements"
351        }
352    }
353}
354
355# ----------------------------------------------------------------------
356# USAGE: dimensions
357#
358# Returns the number of dimensions for this object: 1, 2, or 3.
359# ----------------------------------------------------------------------
360itcl::body Rappture::Mesh::dimensions {} {
361    return $_dim
362}
363
364# ----------------------------------------------------------------------
365# USAGE: limits x|y|z
366#
367# Returns the {min max} values for the limits of the specified axis.
368# ----------------------------------------------------------------------
369itcl::body Rappture::Mesh::limits {which} {
370    if {![info exists _limits(${which}min)]} {
371        error "bad axis \"$which\": should be x, y, z"
372    }
373    return [list $_limits(${which}min) $_limits(${which}max)]
374}
375
376# ----------------------------------------------------------------------
377# USAGE: hints ?<keyword>?
378#
379# Returns a list of key/value pairs for various hints about plotting
380# this field.  If a particular <keyword> is specified, then it returns
381# the hint for that <keyword>, if it exists.
382# ----------------------------------------------------------------------
383itcl::body Rappture::Mesh::hints {{keyword ""}} {
384    foreach key {label color units} {
385        set str [$_mesh get $key]
386        if {"" != $str} {
387            set hints($key) $str
388        }
389    }
390
391    if {$keyword != ""} {
392        if {[info exists hints($keyword)]} {
393            return $hints($keyword)
394        }
395        return ""
396    }
397    return [array get hints]
398}
399
400itcl::body Rappture::Mesh::GetDimension {} {
401    set string [$_xmlobj get dimension]
402    if { [scan $string "%d" _dim] != 1 } {
403        set _dim -1
404        return
405    }
406    if { $_dim < 1 || $_dim > 4 } {
407        error "bad dimension in mesh"
408    } 
409}
410
411itcl::body Rappture::Mesh::GetDouble { path } {
412    set string [$_xmlobj get $path]
413    if { [scan $string "%g" value] != 1 } {
414        puts stderr "can't get double value \"$string\" of \"$path\""
415        return 0.0
416    }
417    return $value
418}
419
420itcl::body Rappture::Mesh::GetInt { path } {
421    set string [$_xmlobj get $path]
422    if { [scan $string "%d" value] != 1 } {
423        puts stderr "can't get integer value \"$string\" of \"$path\""
424        return 0.0
425    }
426    return $value
427}
428
429# ----------------------------------------------------------------------
430# USAGE: _buildNodesElements <xmlobj> <path>
431#
432# Used internally to build a mesh representation based on nodes and
433# elements stored in the XML.
434# ----------------------------------------------------------------------
435itcl::body Rappture::Mesh::_buildNodesElements {xmlobj path} {
436    # create the vtk objects containing points and connectivity
437
438    set _vtkpoints $this-points
439    vtkPoints $_vtkpoints
440
441    #
442    # Extract each node and add it to the points list.
443    #
444    foreach comp [$xmlobj children -type node $path] {
445        set xyz [$xmlobj get $path.$comp]
446
447        # make sure we have x,y,z
448        while {[llength $xyz] < 3} {
449            lappend xyz "0"
450        }
451
452        # extract each point and add it to the points list
453        foreach {x y z} $xyz break
454        foreach dim {x y z} units $_units {
455            set v [Rappture::Units::convert [set $dim] \
456                -context $units -to $units -units off]
457
458            set $dim $v  ;# save back to real x/y/z variable
459
460            if {"" == $_limits(${dim}min)} {
461                set _limits(${dim}min) $v
462                set _limits(${dim}max) $v
463            } else {
464                if {$v < $_limits(${dim}min)} { set _limits(${dim}min) $v }
465                if {$v > $_limits(${dim}max)} { set _limits(${dim}max) $v }
466            }
467        }
468        $_vtkpoints InsertNextPoint $x $y $z
469    }
470
471    set _vtkgrid $this-grid
472    vtkUnstructuredGrid $_vtkgrid
473    $_vtkgrid SetPoints $_vtkpoints
474
475    #
476    # Extract each element and add it to the mesh.
477    #
478    foreach comp [$xmlobj children -type element $path] {
479        set nlist [$_mesh get $comp.nodes]
480        set elem [_getVtkElement [llength $nlist]]
481
482        set i 0
483        foreach n $nlist {
484            [$elem GetPointIds] SetId $i $n
485            incr i
486        }
487        $_vtkgrid InsertNextCell [$elem GetCellType] [$elem GetPointIds]
488    }
489}
490
491# ----------------------------------------------------------------------
492# USAGE: _buildRectMesh <xmlobj> <path>
493#
494# Used internally to build a mesh representation based on a native
495# VTK file for a rectangular mesh.
496# ----------------------------------------------------------------------
497itcl::body Rappture::Mesh::_buildRectMesh {xmlobj path} {
498    vtkRectilinearGridReader $this-gr
499    $this-gr ReadFromInputStringOn
500    $this-gr SetInputString [$xmlobj get $path.vtk]
501
502    set _vtkgrid [$this-gr GetOutput]
503    set _points [$_vtkgrid GetPointData]
504
505    $_vtkgrid Update
506    foreach name {xmin xmax ymin ymax zmin zmax} val [$_vtkgrid GetBounds] {
507        set _limits($name) $val
508    }
509}
510
511# ----------------------------------------------------------------------
512# USAGE: _getVtkElement <npts>
513#
514# Used internally to find (or allocate, if necessary) a vtk element
515# that can be used to build up a mesh.  The element depends on the
516# number of points passed in.  4 points is a tetrahedron, 5 points
517# is a pyramid, etc.
518# ----------------------------------------------------------------------
519itcl::body Rappture::Mesh::_getVtkElement {npts} {
520    if {![info exists _pts2elem($npts)]} {
521        switch -- $npts {
522            3 {
523                set _pts2elem($npts) $this-elem3
524                vtkTriangle $_pts2elem($npts)
525            }
526            4 {
527                set _pts2elem($npts) $this-elem4
528                vtkTetra $_pts2elem($npts)
529            }
530            5 {
531                set _pts2elem($npts) $this-elem5
532                vtkPyramid $_pts2elem($npts)
533            }
534            6 {
535                set _pts2elem($npts) $this-elem6
536                vtkWedge $_pts2elem($npts)
537            }
538            8 {
539                set _pts2elem($npts) $this-elem8
540                vtkVoxel $_pts2elem($npts)
541            }
542        }
543    }
544    return $_pts2elem($npts)
545}
546
547itcl::body Rappture::Mesh::ReadCloud { xmlobj path } {
548    set _type "cloud"
549    set data [$xmlobj get $path.cloud.points]
550    Rappture::ReadPoints $data _dim values
551       
552    set all [blt::vector create \#auto]
553    $all set $values
554    set xv [blt::vector create \#auto]
555    set yv [blt::vector create \#auto]
556    set zv [blt::vector create \#auto]
557    if { $_dim == 2 } {
558        $all split $xv $yv
559        $zv seq 0 0 [$xv length]
560    } elseif { $_dim == 3 } {
561        $all split $xv $yv $zv
562    } else {
563        error "invalid dimension \"$_dim\" of cloud mesh"
564    }
565    set _numPoints [$xv length]
566    if { $_dim == 2 } {
567        append out "DATASET POLYDATA\n"
568        append out "POINTS $_numPoints float\n"
569        foreach x [$xv range 0 end] y [$yv range 0 end]  {
570            append out "$x $y 0\n"
571        }
572        append out "POINT_DATA $_numPoints\n"
573        set _vtkdata $out
574        foreach axis {x y} {
575            set vector [set ${axis}v]
576            set _limits(${axis}min) [${vector} min]
577            set _limits(${axis}max) [${vector} max]
578        }
579    } elseif { $_dim == 3 } {
580        append out "DATASET POLYDATA\n"
581        append out "POINTS $_numPoints float\n"
582        foreach x [$xv range 0 end] y [$yv range 0 end] z [$zv range 0 end] {
583            append out "$x $y $z\n"
584        }
585        append out "POINT_DATA $_numPoints\n"
586        set _vtkdata $out
587        foreach axis {x y z} {
588            set vector [set ${axis}v]
589            set _limits(${axis}min) [${vector} min]
590            set _limits(${axis}max) [${vector} max]
591        }
592    }
593    blt::vector destroy $xv $yv $zv $all
594}
595
596itcl::body Rappture::Mesh::ReadVtk { xmlobj path } {
597    set _type "vtk"
598    if { $_dim < 1 } {
599        error "can't read vtk mesh: no dimension specified"
600    }
601
602    # Write the contents to a file just in case it's binary.
603    set tmpfile file[pid].vtk
604    set f [open "$tmpfile" "w"]
605    fconfigure $f -translation binary -encoding binary
606    puts $f [$xmlobj get vtk]
607    close $f
608
609    # Read the data into a polydata output
610    set _vtkreader $this-datasetreader
611    vtkDataSetReader $_vtkreader
612    $_vtkreader SetFileName $tmpfile
613    $_vtkreader Update
614    set _vtkoutput [$_vtkreader GetOutput]
615    set _vtkpoints [$_vtkoutput GetPointData]
616    file delete $tmpfile
617    set _isValid 1
618}
619
620itcl::body Rappture::Mesh::ReadGrid { xmlobj path } {
621    set _type "grid"
622    set numUniform 0
623    set numNonUniform 0
624    foreach axis { x y z } {
625        set min    [$xmlobj get "$path.grid.${axis}axis.min"]
626        set max    [$xmlobj get "$path.grid.${axis}axis.max"]
627        set num    [$xmlobj get "$path.grid.${axis}axis.numpoints"]
628        set values [$xmlobj get "$path.grid.${axis}values"]
629        if { $min != "" && $max != "" && $num != "" && $num > 0 } {
630            set ${axis}Min $min
631            set ${axis}Max $max
632            set ${axis}Num $num
633            incr numUniform
634        } elseif { $values != "" } {
635            incr numNonUniform
636            set ${axis}Values $values
637        }
638    }
639    set _dim [expr $numNonUniform + $numUniform]
640    if { $_dim == 0 } {
641        # No data found.
642        return 0
643    }
644    if { $numNonUniform == 0 && $numUniform > 0 } {
645        # Specially handle the case where all axes are uniform.  This means we
646        # can use a vtkImageData to represent the data instead of a
647        # vtkRectlinearGrid.  Vtk points are always 3D
648        if { $_dim == 2 } {
649            append out "DATASET STRUCTURED_POINTS\n"
650            append out "DIMENSIONS $xNum $yNum 1\n"
651            set xSpace [expr ($xMax - $xMin) / double($xNum - 1)]
652            set ySpace [expr ($yMax - $yMin) / double($yNum - 1)]
653            append out "SPACING $xSpace $ySpace 0\n"
654            append out "ORIGIN 0 0 0\n"
655            set _numPoints [expr $xNum * $yNum]
656            append out "POINT_DATA $_numPoints\n"
657            set _vtkdata $out
658            foreach axis {x y} {
659                set _limits(${axis}min) [set ${axis}Min]
660                set _limits(${axis}max) [set ${axis}Max]
661            }
662        } elseif { $_dim == 3 } {
663            append out "DATASET STRUCTURED_POINTS\n"
664            append out "DIMENSIONS $xNum $yNum $zNum\n"
665            set xSpace [expr ($xMax - $xMin) / double($xNum - 1)]
666            set ySpace [expr ($yMax - $yMin) / double($yNum - 1)]
667            set zSpace [expr ($zMax - $zMin) / double($zNum - 1)]
668            append out "SPACING $xSpace $ySpace $zSpace\n"
669            append out "ORIGIN 0 0 0\n"
670            set _numPoints [expr $xNum * $yNum * zNum]
671            append out "POINT_DATA $_numPoints\n"
672            set _vtkdata $out
673            foreach axis {x y z} {
674                set _limits(${axis}min) [set ${axis}Min]
675                set _limits(${axis}max) [set ${axis}Max]
676            }
677        } else {
678            error "bad dimension of mesh \"$_dim\""
679        }
680        set _isValid 1
681        return 1
682    }
683
684    set xv [blt::vector create \#auto]
685    if { [info exists xMin] } {
686        $xv seq $xMin $xMax $xNum
687    } else {
688        $xv set [$xmlobj get $path.grid.xvalues]
689        set xMin [$xv min]
690        set xMax [$xv max]
691        set xNum [$xv length]
692    }
693    set yv [blt::vector create \#auto]
694    if { [info exists yMin] } {
695        $yv seq $yMin $yMax $yNum
696    } else {
697        $yv set [$xmlobj get $path.grid.yvalues]
698        set yMin [$yv min]
699        set yMax [$yv max]
700        set yNum [$yv length]
701    }
702    set zv [blt::vector create \#auto]
703    if { $_dim == 3 } {
704        if { [info exists zMin] } {
705            $zv seq $zMin $zMax $zNum
706        }  else {
707            $zv set [$xmlobj get $path.grid.zvalues]
708            set zMin [$zv min]
709            set zMax [$zv max]
710            set zNum [$zv length]
711        }
712    } else {
713        set zNum 1
714    }
715    if { $_dim == 3 } {
716        set _numPoints [expr $xNum * $yNum * $zNum]
717        append out "DATASET RECTILINEAR_GRID\n"
718        append out "DIMENSIONS $xNum $yNum $zNum\n"
719        append out "X_COORDINATES $xNum double\n"
720        append out [$xv range 0 end]
721        append out "\n"
722        append out "Y_COORDINATES $yNum double\n"
723        append out [$yv range 0 end]
724        append out "\n"
725        append out "Z_COORDINATES $zNum double\n"
726        append out [$zv range 0 end]
727        append out "\n"
728        append out "POINT_DATA $_numPoints\n"
729        set _vtkdata $out
730        foreach axis {x y z} {
731            if { [info exists ${axis}Min] } {
732                set _limits(${axis}min) [set ${axis}Min]
733                set _limits(${axis}max) [set ${axis}Max]
734            }
735        }
736    } elseif { $_dim == 2 } {
737        set _numPoints [expr $xNum * $yNum]
738        append out "DATASET RECTILINEAR_GRID\n"
739        append out "DIMENSIONS $xNum $yNum 1\n"
740        append out "X_COORDINATES $xNum double\n"
741        append out [$xv range 0 end]
742        append out "\n"
743        append out "Y_COORDINATES $yNum double\n"
744        append out [$yv range 0 end]
745        append out "\n"
746        append out "Z_COORDINATES 1 double\n"
747        append out "0\n"
748        append out "POINT_DATA $_numPoints\n"
749        foreach axis {x y} {
750            if { [info exists ${axis}Min] } {
751                set _limits(${axis}min) [set ${axis}Min]
752                set _limits(${axis}max) [set ${axis}Max]
753            }
754        }
755        set _vtkdata $out
756    } else {
757        error "invalid dimension of grid \"$_dim\""
758    }
759    blt::vector destroy $xv $yv $zv
760    set _isValid 1
761    return 1
762}
763
764itcl::body Rappture::Mesh::ReadCells { xmlobj path } {
765    set _type "cells"
766    set celltype [$xmlobj get $path.cells.celltype]
767    array set type2vertices {
768        triangle 3
769        quad     4
770        tetraheadron 4
771        voxel 8
772        hexaheadron 8
773        wedge 6
774        pyramid 5
775    }
776    array set type2dim {
777        triangle 2
778        quad     2
779        tetraheadron 3
780        voxel 3
781        hexaheadron 3
782        wedge 3
783        pyramid 3
784    }
785    array set type2vtkid {
786        triangle 5
787        quad     9
788        tetraheadron 10
789        voxel 11
790        hexaheadron 12
791        wedge 13
792        pyramid 14
793    }
794    if { ![info exists type2vertices($celltype)] } {
795        puts stderr "unknown celltype \"$celltype\""
796        return 0
797    }
798    set numVertices $type2vertices($celltype)
799    set _dim $type2dim($celltype)
800    set vtkid $type2vtkid($celltype)
801   
802    set points [$xmlobj get $path.cells.points]
803    set xvalues [$xmlobj get $path.cells.xvalues]
804    set yvalues [$xmlobj get $path.cells.yvalues]
805    if { $_dim == 2 } {
806        if { $xvalues != "" && $yvalues != "" } {
807            set xv [blt::vector create \#auto]
808            set yv [blt::vector create \#auto]
809            $xv set $xvalues
810            $yx set $yvalues
811        } elseif { $points != "" } {
812            set all [blt::vector create \#auto]
813            $all set $points
814            set xv [blt::vector create \#auto]
815            set yv [blt::vector create \#auto]
816            $all split $xv $yv
817        } else {
818            puts stderr "missing either xvalues, yvalues, or points for $celltype cell mesh"
819            return 0
820        }
821        set _numPoints [$xv length]
822        set numIndices 0
823        set data {}
824        set indices [$xmlobj get $path.cells.indices]
825        foreach { a b c } $indices {
826            append data "    $numVertices $a $b $c\n"
827            incr numIndices
828        }
829        set celltypes {}
830        for { set i 0 } { $i < $numIndices } { incr i } {
831            append celltypes "$vtkid\n"
832        }
833        append out "DATASET UNSTRUCTURED_GRID\n"
834        append out "POINTS $_numPoints float\n"
835        foreach x [$xv range 0 end] y [$yv range 0 end] {
836            append out "    $x $y 0\n"
837        }
838        append out "CELLS $numIndices [expr $numIndices * 4]\n"
839        append out $data
840        append out "CELL_TYPES $numIndices\n"
841        append out $celltypes
842        append out "POINT_DATA $_numPoints\n"
843        set _limits(xmin) [$xv min]
844        set _limits(xmax) [$xv max]
845        set _limits(ymin) [$yv min]
846        set _limits(ymax) [$yv max]
847    } else {
848        set zvalues [$xmlobj get $path.cells.zvalues]
849        if { $xvalues != "" && $yvalues != ""  && $zvalues != "" } {
850            set xv [blt::vector create \#auto]
851            set yv [blt::vector create \#auto]
852            set zv [blt::vector create \#auto]
853            $xv set $xvalues
854            $yx set $yvalues
855            $zx set $yvalues
856        } elseif { $points != "" } {
857            set all [blt::vector create \#auto]
858            $all set $points
859            set xv [blt::vector create \#auto]
860            set yv [blt::vector create \#auto]
861            set zv [blt::vector create \#auto]
862            $all split $xv $yv $zv
863        } else {
864            puts stderr "missing either xvalues, yvalues, zvalues, or points for $celltype cell mesh"
865            return 0
866        }
867        set _numPoints [$xv length]
868        set data {}
869        set indices [$xmlobj get $path.cells.indices]
870        set numIndices 0
871        foreach { a b c } $indices {
872            append data "    $numVertices $a $b $c\n"
873            incr numIndices
874        }
875        set celltypes {}
876        for { set i 0 } { $i < $numIndices } { incr i } {
877            append celltypes "$vtkid\n"
878        }
879        append out "DATASET UNSTRUCTURED_GRID\n"
880        append out "POINTS $_numPoints float\n"
881        foreach x [$xv range 0 end] y [$yv range 0 end] z [$zv range 0 end] {
882            append out "    $x $y $z\n"
883        }
884        append out "CELLS $numIndices [expr $numIndices * 4]\n"
885        append out $data
886        append out "CELL_TYPES $numIndices\n"
887        append out $celltypes
888        append out "POINT_DATA $_numPoints\n"
889        set _limits(xmin) [$xv min]
890        set _limits(xmax) [$xv max]
891        set _limits(ymin) [$yv min]
892        set _limits(ymax) [$yv max]
893        set _limits(zmin) [$zv min]
894        set _limits(zmax) [$zv max]
895    }
896    set _vtkdata $out
897    set _isValid 1
898}
899
900itcl::body Rappture::Mesh::ReadTriangles { xmlobj path } {
901    set _type "triangles"
902    set _dim 2
903    set triangles [$xmlobj get $path.triangles.indices]
904    if { $triangles == "" } {
905        puts stderr "no triangle indices specified in mesh"
906        return 0
907    }
908    set xvalues [$xmlobj get $path.triangles.xvalues]
909    set yvalues [$xmlobj get $path.triangles.yvalues]
910    set points [$xmlobj get $path.triangles.points]
911    if { $xvalues != "" && $yvalues != "" } {
912        set xv [blt::vector create \#auto]
913        set yv [blt::vector create \#auto]
914        $xv set $xvalues
915        $yx set $yvalues
916    } elseif { $points != "" } {
917        set all [blt::vector create \#auto]
918        $all set $points
919        set xv [blt::vector create \#auto]
920        set yv [blt::vector create \#auto]
921        $all split $xv $yv
922    } else {
923        puts stderr "missing either xvalues, yvalues, or points for triangular mesh"
924        return 0
925    }
926    set _numPoints [$xv length]
927    set count 0
928    set data {}
929    foreach { a b c } $triangles {
930        append data "    3 $a $b $c\n"
931        incr count
932    }
933    set celltypes {}
934    for { set i 0 } { $i < $count } { incr i } {
935        append celltypes "5\n"
936    }
937    append out "DATASET UNSTRUCTURED_GRID\n"
938    append out "POINTS $_numPoints float\n"
939    foreach x [$xv range 0 end] y [$yv range 0 end] {
940        append out "    $x $y 0\n"
941    }
942    append out "CELLS $count [expr $count * 4]\n"
943    append out $data
944    append out "CELL_TYPES $count\n"
945    append out $celltypes
946    append out "POINT_DATA $_numPoints\n"
947    set _limits(xmin) [$xv min]
948    set _limits(xmax) [$xv max]
949    set _limits(ymin) [$yv min]
950    set _limits(ymax) [$yv max]
951    set _vtkdata $out
952    set _isValid 1
953}
954
955itcl::body Rappture::Mesh::ReadBltGridMesh { path } {
956    if { $_dim != 2 } {
957        error "must be 2D mesh for BLT"
958    }
959    # Vtk points are always 3D
960    set xNum 0
961    if { $_xValues != "" } {
962        blt::vector destroy $_xValues
963    }
964    set yNum 0
965    if { $_yValues != "" } {
966        blt::vector destroy $_yValues
967    }
968    set string [$xmlobj get grid.xvalues]
969    set _xValues [blt::vector create \#auto]
970    if { [llength $string] == 0 } {
971        set xMin [GetDouble grid.xaxis.min]
972        set xMax [GetDouble grid.xaxis.max]
973        set xNum [GetInt grid.xaxis.numpoints]
974        $xValues seq $xMin $xMax $xNum
975    } else {
976        $xValues set $string
977        set xNum [$xValues length]
978        set xMax [$xValues max]
979        set xMin [$xValues min]
980    }
981    set _yValues [blt::vector create \#auto]
982    set string [$xmlobj get grid.yvalues]
983    if { [llength $string] == 0 } {
984        set string [$xmlobj get grid.yaxis.min]
985        if { $string != "" } {
986            set yMin [GetDouble grid.yaxis.min]
987            set yMax [GetDouble grid.yaxis.max]
988            set yNum [GetInt grid.yaxis.numpoints]
989            $yValues seq $yMin $yMax $yNum
990        }
991    } else {
992        set yValues [blt::vector create \#auto]
993        $yValues set $string
994        set yNum [$yValues length]
995        set yMax [$yValues max]
996        set yMin [$yValues min]
997    }
998    if { [$yValues length] == 0 } {
999        # Not 2D
1000    }
1001    foreach axis {x y} {
1002        set _limits(${axis}min) [set ${axis}Min]
1003        set _limits(${axis}max) [set ${axis}Max]
1004    }
1005    set mesh [blt::mesh create irregular]
1006    blt::mesh configure $mesh -x $_xValues -y $_yValues
1007    return $mesh
1008}
1009
1010itcl::body Rappture::Mesh::ReadBltTriangleMesh { path } {
1011    set _xv [blt::vector create \#auto]
1012    set _yv [blt::vector create \#auto]
1013    $xv set [$_xmlobj get $path.triangle.xaxis.values]
1014    $yv set [$_xmlobj get $path.triangle.yaxis.values]
1015    set triangles [$_xmlobj get $path.triangle.triangles]
1016    set _limits(xmin) [$_xv min]
1017    set _limits(xmax) [$_xv max]
1018    set _limits(ymin) [$_yv min]
1019    set _limits(ymax) [$_yv max]
1020    set mesh [blt::mesh create triangle]
1021    blt::mesh configure $mesh -x $xv -y $yv -triangles $triangles
1022    return $mesh
1023}
1024
1025itcl::body Rappture::Mesh::ReadBltCloudMesh { path } {
1026    # Read the points into a vector and split them into individual coordinate
1027    # arrays.
1028    set all [blt::vector create \#auto]
1029    $all set [$_xmlobj get $path.cloud.points]
1030    if { [$all length] & 0x1 } {
1031        error "odd number of points in 2D cloud mesh"
1032    }
1033    if { $_xv != "" } {
1034        blt::vector destroy $_xv
1035    }
1036    if { $_yv != "" } {
1037        blt::vector destroy $_yv
1038    }
1039    set _xv [blt::vector create \#auto]
1040    set _yv [blt::vector create \#auto]
1041    $all split $_xv $_yv
1042    set _limits(xmin) [$_xv min]
1043    set _limits(xmax) [$_xv max]
1044    set _limits(ymin) [$_yv min]
1045    set _limits(ymax) [$_yv max]
1046    set mesh [blt::mesh create cloud]
1047    blt::mesh configure $mesh -x $_xv -y $_yv
1048    blt::vector destroy $all
1049    return $mesh
1050}
1051
1052itcl::body Rappture::Mesh::ReadBltMesh { path } {
1053    if { $_dim != 2 } {
1054        error "dimension of mesh $_dim is not 2"
1055    }
1056    switch -- $_type {
1057        "cloud" {
1058            set _bltmesh [ReadBltCloudMesh $path ]
1059        }
1060        "unirect" {
1061            set _bltmesh [ReadBltRegularMesh $path]
1062        }
1063        "rectangle" {
1064            set _bltmesh [ReadBltIrregularMesh $path]
1065        }
1066        "triangle" {
1067            set _bltmesh [ReadBltTriangleMesh $path]
1068        }
1069        default {
1070            error "bad mesh type \"$_type\": should be cloud, unirect, rectangle, or triangle"
1071        }
1072    }
1073}
1074
1075# ----------------------------------------------------------------------
1076# USAGE: ReadNodesElements <xmlobj> <path>
1077#
1078# Used internally to build a mesh representation based on nodes and
1079# elements stored in the XML.
1080# ----------------------------------------------------------------------
1081itcl::body Rappture::Mesh::ReadNodesElements {xmlobj path} {
1082    set type "nodeselements"
1083    # create the vtk objects containing points and connectivity
1084    set _dim 2
1085    #
1086    # Extract each node and add it to the points list.
1087    #
1088    set points {}
1089    set count 0
1090    foreach comp [$xmlobj children -type node $path] {
1091        set xyz [$xmlobj get $path.$comp]
1092
1093        # make sure we have x,y,z
1094        while {[llength $xyz] < 3} {
1095            lappend xyz "0"
1096        }
1097
1098        # extract each point and add it to the points list
1099        foreach {x y z} $xyz break
1100        foreach dim {x y z} units $_units {
1101            set v [Rappture::Units::convert [set $dim] \
1102                -context $units -to $units -units off]
1103
1104            set $dim $v  ;# save back to real x/y/z variable
1105
1106            if {"" == $_limits(${dim}min)} {
1107                set _limits(${dim}min) $v
1108                set _limits(${dim}max) $v
1109            } else {
1110                if {$v < $_limits(${dim}min)} { set _limits(${dim}min) $v }
1111                if {$v > $_limits(${dim}max)} { set _limits(${dim}max) $v }
1112            }
1113        }
1114        append points "   $x $y $z\n"
1115        incr count
1116    }
1117    set _numPoints $count
1118    #
1119    # Extract each element and add it to the mesh.
1120    #
1121    array set node2celltype {
1122        3 5
1123        4 10
1124        8 11
1125        6 13
1126        5 14
1127    }
1128    set count 0
1129    set celltypes {}
1130    set data {}
1131    foreach comp [$xmlobj children -type element $path] {
1132        set nlist [$_mesh get $comp.nodes]
1133        set numNodes [llength $nlist]
1134        if { ![info exists node2celltype($numNodes)] } {
1135            return 0
1136        }
1137        set celltype $node2celltype($numNodes)
1138        append celltypes "  $celltype\n"
1139        append data "  $numNodes $nlist\n"
1140        incr count
1141    }
1142    append out "DATASET UNSTRUCTURED_GRID\n"
1143    append out "POINTS $_numPoints float\n"
1144    append out $points
1145    append out "CELLS $count [expr $count * 4]\n"
1146    append out $data
1147    append out "CELL_TYPES $count\n"
1148    append out $celltypes
1149    append out "\n"
1150    append out "POINT_DATA $_numPoints\n"
1151    foreach axis {x y} {
1152        if { [info exists ${axis}Min] } {
1153            set _limits(${axis}min) [set ${axis}Min]
1154            set _limits(${axis}max) [set ${axis}Max]
1155        }
1156    }
1157    set _vtkdata $out
1158    set _isValid 1
1159}
Note: See TracBrowser for help on using the repository browser.