source: branches/1.3/gui/scripts/mesh.tcl @ 4474

Last change on this file since 4474 was 4474, checked in by ldelgass, 10 years ago

Merge mesh/field fixes from trunk, add mesh viewer
File size: 42.6 KB
Line 
1# -*- mode: tcl; indent-tabs-mode: nil -*-
2
3# ----------------------------------------------------------------------
4#  COMPONENT: mesh - represents a structured mesh for a device
5#
6#  This object represents a mesh in an XML description of simulator
7#  output.  A mesh is a structured arrangement of points, as elements
8#  composed of nodes representing coordinates.  This is a little
9#  different from a cloud, which is an unstructured arrangement
10#  (shotgun blast) of points.
11# ======================================================================
12#  AUTHOR:  Michael McLennan, Purdue University
13#  Copyright (c) 2004-2012  HUBzero Foundation, LLC
14#
15#  See the file "license.terms" for information on usage and
16#  redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
17# ======================================================================
18package require Itcl
19
20namespace eval Rappture {
21    # forward declaration
22}
23
24itcl::class Rappture::Mesh {
25    private variable _xmlobj ""  ;      # Ref to XML obj with device data
26    private variable _mesh ""    ;      # Lib obj representing this mesh
27    private variable _dim       0;      # Dimension of mesh (1, 2, or 3)
28    private variable _type "";          # Indicates the type of mesh.
29    private variable _axis2units;       # System of units for x, y, z
30    private variable _axis2labels;      #
31    private variable _hints
32    private variable _limits        ;   # Array of mesh limits. Keys are
33                                        # xmin, xmax, ymin, ymax, ...
34    private variable _numPoints 0   ;   # # of points in mesh
35    private variable _numCells 0   ;    # # of cells in mesh
36    private variable _vtkdata "";       # Mesh in vtk file format.
37    private variable _isValid 0;        # Indicates if the mesh is valid.
38    constructor {xmlobj path} {
39        # defined below
40    }
41    destructor {
42        # defined below
43    }
44    public method points {}
45    public method elements {}
46    public method mesh {{-type "vtk"}}
47    public method size {{what -points}}
48    public method dimensions {}
49    public method limits {which}
50    public method units { axis }
51    public method label { axis }
52    public method hints {{key ""}}
53    public method isvalid {} {
54        return $_isValid
55    }
56    public proc fetch {xmlobj path}
57    public proc release {obj}
58    public method vtkdata {{what -partial}}
59    public method type {} {
60        return $_type
61    }
62    public method numpoints {} {
63        return $_numPoints
64    }
65    public method numcells {} {
66        return $_numCells
67    }
68
69    private common _xp2obj       ;      # used for fetch/release ref counting
70    private common _obj2ref      ;      # used for fetch/release ref counting
71    private variable _xv        ""
72    private variable _yv        ""
73    private variable _zv        ""
74    private variable _xCoords   "";     # For the blt contour only
75    private variable _yCoords   "";     # For the blt contour only
76   
77    private method ReadNodesElements {path}
78    private method GetDimension { path }
79    private method GetDouble { path }
80    private method GetInt { path }
81    private method InitHints {}
82    private method ReadGrid { path }
83    private method ReadUnstructuredGrid { path }
84    private method ReadVtk { path }
85    private method WriteTriangles { path xv yv zv triangles }
86    private method WriteQuads { path xv yv zv quads }
87    private method WriteTetrahedrons { path xv yv zv tetrahedrons }
88    private method WriteHexahedrons { path xv yv zv hexhedrons }
89    private method WriteWedges { path xv yv zv wedges }
90    private method WritePyramids { path xv yv zv pyramids }
91    private method WriteHybridCells { path xv yv zv cells celltypes }
92    private method WritePointCloud { path xv yv zv }
93    private method GetCellType { name }
94    private method GetNumIndices { type }
95}
96
97# ----------------------------------------------------------------------
98# USAGE: Rappture::Mesh::fetch <xmlobj> <path>
99#
100# Clients use this instead of a constructor to fetch the Mesh for
101# a particular <path> in the <xmlobj>.  When the client is done with
102# the mesh, he calls "release" to decrement the reference count.
103# When the mesh is no longer needed, it is cleaned up automatically.
104# ----------------------------------------------------------------------
105itcl::body Rappture::Mesh::fetch {xmlobj path} {
106    set handle "$xmlobj|$path"
107    if {[info exists _xp2obj($handle)]} {
108        set obj $_xp2obj($handle)
109        incr _obj2ref($obj)
110        return $obj
111    }
112    set obj [Rappture::Mesh ::#auto $xmlobj $path]
113    set _xp2obj($handle) $obj
114    set _obj2ref($obj) 1
115    return $obj
116}
117
118# ----------------------------------------------------------------------
119# USAGE: Rappture::Mesh::release <obj>
120#
121# Clients call this when they're no longer using a Mesh fetched
122# previously by the "fetch" proc.  This decrements the reference
123# count for the mesh and destroys the object when it is no longer
124# in use.
125# ----------------------------------------------------------------------
126itcl::body Rappture::Mesh::release {obj} {
127    if {[info exists _obj2ref($obj)]} {
128        incr _obj2ref($obj) -1
129        if {$_obj2ref($obj) <= 0} {
130            unset _obj2ref($obj)
131            foreach handle [array names _xp2obj] {
132                if {$_xp2obj($handle) == $obj} {
133                    unset _xp2obj($handle)
134                }
135            }
136            itcl::delete object $obj
137        }
138    } else {
139        error "can't find reference count for $obj"
140    }
141}
142
143# ----------------------------------------------------------------------
144# CONSTRUCTOR
145# ----------------------------------------------------------------------
146itcl::body Rappture::Mesh::constructor {xmlobj path} {
147    package require vtk
148    if {![Rappture::library isvalid $xmlobj]} {
149        error "bad value \"$xmlobj\": should be Rappture::library"
150    }
151    set _xmlobj $xmlobj
152    set _mesh [$xmlobj element -as object $path]
153
154    # Initialize mesh bounds to empty
155    foreach axis {x y z} {
156        set _limits($axis) ""
157    }
158    set units [$_mesh get units]
159    set first [lindex $units 0]
160    foreach u $units axis { x y z } {
161        if { $u != "" } {
162            set _axis2units($axis) $u
163        } else {
164            set _axis2units($axis) $first
165        }
166    }
167    foreach label [$_mesh get labels] axis { x y z } {
168        if { $label != "" } {
169            set _axis2labels($axis) $label
170        } else {
171            set _axis2labels($axis) [string toupper $axis]
172        }
173    }
174
175    # Meshes comes in a variety of flavors
176    #
177    # Dimensionality is determined from the <dimension> tag. 
178    #
179    # <vtk> described mesh
180    # <element> +  <node> definitions
181    # <grid>            rectangular mesh
182    # <unstructured>    homogeneous cell type mesh.
183
184    # Check that only one mesh type was defined.
185    set subcount 0
186    foreach cname [$_mesh children] {
187        foreach type { vtk grid unstructured } {
188            if { $cname == $type } {
189                incr subcount
190                break
191            }
192        }
193    }
194    if {[$_mesh element "node"] != "" ||
195        [$_mesh element "element"] != ""} {
196        incr subcount
197    }
198
199    if { $subcount == 0 } {
200        puts stderr "WARNING: no mesh specified for \"$path\"."
201        return
202    }
203    if { $subcount > 1 } {
204        puts stderr "WARNING: too many mesh types specified for \"$path\"."
205        return
206    }
207    set result 0
208    if { [$_mesh element "vtk"] != ""} {
209        set result [ReadVtk $path]
210    } elseif {[$_mesh element "grid"] != "" } {
211        set result [ReadGrid $path]
212    } elseif {[$_mesh element "unstructured"] != "" } {
213        set result [ReadUnstructuredGrid $path]
214    } elseif {[$_mesh element "node"] != "" && [$_mesh element "element"] != ""} {
215        set result [ReadNodesElements $path]
216    }
217    set _isValid $result
218    InitHints
219}
220
221# ----------------------------------------------------------------------
222# DESTRUCTOR
223# ----------------------------------------------------------------------
224itcl::body Rappture::Mesh::destructor {} {
225    # don't destroy the _xmlobj! we don't own it!
226    itcl::delete object $_mesh
227
228    if { $_xCoords != "" } {
229        blt::vector destroy $_xCoords
230    }
231    if { $_yCoords != "" } {
232        blt::vector destroy $_yCoords
233    }
234}
235
236#
237# vtkdata --
238#
239#       This is called by the field object to generate a VTK file to send to
240#       the remote render server.  Returns the vtkDataSet object containing
241#       (at this point) just the mesh.  The field object doesn't know (or
242#       care) what type of mesh is used.  The field object will add field
243#       arrays before generating output to send to the remote render server.
244#
245itcl::body Rappture::Mesh::vtkdata {{what -partial}} {
246    if {$what == "-full"} {
247        append out "# vtk DataFile Version 3.0\n"
248        append out "[hints label]\n"
249        append out "ASCII\n"
250        append out $_vtkdata
251        return $out
252    } else {
253        return $_vtkdata
254    }
255}
256
257# ----------------------------------------------------------------------
258# USAGE: points
259#
260# Returns the vtk object containing the points for this mesh.
261# ----------------------------------------------------------------------
262itcl::body Rappture::Mesh::points {} {
263    return ""
264}
265
266#
267# units --
268#
269#       Returns the units of the given axis.
270#
271itcl::body Rappture::Mesh::units { axis } {
272    if { ![info exists _axis2units($axis)] } {
273        return ""
274    }
275    return $_axis2units($axis)
276}
277
278#
279# label --
280#
281#       Returns the label of the given axis.
282#
283itcl::body Rappture::Mesh::label { axis } {
284    if { ![info exists _axis2labels($axis)] } {
285        return ""
286    }
287    return $_axis2labels($axis)
288}
289
290# ----------------------------------------------------------------------
291# USAGE: elements
292#
293# Returns a list of the form {plist r plist r ...} for each element
294# in this mesh.  Each plist is a list of {x y x y ...} coordinates
295# for the mesh.
296# ----------------------------------------------------------------------
297itcl::body Rappture::Mesh::elements {} {
298    # build a map for region number => region type
299    foreach comp [$_mesh children -type region] {
300        set id [$_mesh element -as id $comp]
301        set regions($id) [$_mesh get $comp]
302    }
303    set regions() "unknown"
304
305    set rlist ""
306    foreach comp [$_mesh children -type element] {
307        set rid [$_mesh get $comp.region]
308
309        #
310        # HACK ALERT!
311        #
312        # Prophet puts out nodes in a funny "Z" shaped order,
313        # not in proper clockwise fashion.  Switch the last
314        # two nodes for now to make them correct.
315        #
316        set nlist [$_mesh get $comp.nodes]
317        set n2 [lindex $nlist 2]
318        set n3 [lindex $nlist 3]
319        set nlist [lreplace $nlist 2 3 $n3 $n2]
320        lappend nlist [lindex $nlist 0]
321
322        set plist ""
323        foreach nid $nlist {
324            eval lappend plist [$_mesh get node($nid)]
325        }
326        lappend rlist $plist $regions($rid)
327    }
328    return $rlist
329}
330
331# ----------------------------------------------------------------------
332# USAGE: mesh
333#
334# Returns the vtk object representing the mesh.
335# ----------------------------------------------------------------------
336itcl::body Rappture::Mesh::mesh { {type "vtk"} } {
337    switch $type {
338        "vtk" {
339            return ""
340        }
341        default {
342            error "Requested mesh type \"$type\" is unknown."
343        }
344    }
345}
346
347# ----------------------------------------------------------------------
348# USAGE: size ?-points|-elements?
349#
350# Returns the number of points in this mesh.
351# ----------------------------------------------------------------------
352itcl::body Rappture::Mesh::size {{what -points}} {
353    switch -- $what {
354        -points {
355            return $_numPoints
356        }
357        -elements {
358            return $_numCells
359        }
360        default {
361            error "bad option \"$what\": should be -points or -elements"
362        }
363    }
364}
365
366# ----------------------------------------------------------------------
367# USAGE: dimensions
368#
369# Returns the number of dimensions for this object: 1, 2, or 3.
370# ----------------------------------------------------------------------
371itcl::body Rappture::Mesh::dimensions {} {
372    return $_dim
373}
374
375# ----------------------------------------------------------------------
376# USAGE: limits x|y|z
377#
378# Returns the {min max} coords for the limits of the specified axis.
379# ----------------------------------------------------------------------
380itcl::body Rappture::Mesh::limits {axis} {
381    if {![info exists _limits($axis)]} {
382        error "bad axis \"$which\": should be x, y, z"
383    }
384    return $_limits($axis)
385}
386
387# ----------------------------------------------------------------------
388# USAGE: hints ?<keyword>?
389#
390# Returns a list of key/value pairs for various hints about plotting
391# this field.  If a particular <keyword> is specified, then it returns
392# the hint for that <keyword>, if it exists.
393# ----------------------------------------------------------------------
394itcl::body Rappture::Mesh::hints {{keyword ""}} {
395    if {$keyword != ""} {
396        if {[info exists _hints($keyword)]} {
397            return $_hints($keyword)
398        }
399        return ""
400    }
401    return [array get _hints]
402}
403
404# ----------------------------------------------------------------------
405# USAGE: InitHints
406#
407# Returns a list of key/value pairs for various hints about plotting
408# this mesh.  If a particular <keyword> is specified, then it returns
409# the hint for that <keyword>, if it exists.
410# ----------------------------------------------------------------------
411itcl::body Rappture::Mesh::InitHints {} {
412    foreach {key path} {
413        camera       camera.position
414        color        about.color
415        label        about.label
416        style        about.style
417        units        units
418    } {
419        set str [$_mesh get $path]
420        if {"" != $str} {
421            set _hints($key) $str
422        }
423    }
424}
425
426itcl::body Rappture::Mesh::GetDimension { path } {
427    set string [$_xmlobj get $path.dim]
428    if { $string == "" } {
429        puts stderr "WARNING: no tag <dim> found in mesh \"$path\"."
430        return 0
431    }
432    if { [scan $string "%d" _dim] == 1 } {
433        if { $_dim == 1 || $_dim == 2 || $_dim == 3 } {
434            return 1
435        }
436    }
437    puts stderr "WARNING: bad <dim> tag value \"$string\": should be 1, 2 or 3."
438    return 0
439}
440
441itcl::body Rappture::Mesh::GetDouble { path } {
442    set string [$_xmlobj get $path]
443    if { [scan $string "%g" value] != 1 } {
444        puts stderr "ERROR: can't get double value \"$string\" of \"$path\""
445        return 0.0
446    }
447    return $value
448}
449
450itcl::body Rappture::Mesh::GetInt { path } {
451    set string [$_xmlobj get $path]
452    if { [scan $string "%d" value] != 1 } {
453        puts stderr "ERROR: can't get integer value \"$string\" of \"$path\""
454        return 0.0
455    }
456    return $value
457}
458
459itcl::body Rappture::Mesh::ReadVtk { path } {
460    set _type "vtk"
461
462    if { ![GetDimension $path] } {
463        return 0
464    }
465    # Create a VTK file with the mesh in it. 
466    set _vtkdata [$_xmlobj get $path.vtk]
467    append out "# vtk DataFile Version 3.0\n"
468    append out "mesh\n"
469    append out "ASCII\n"
470    append out "$_vtkdata\n"
471
472     # Write the contents to a file just in case it's binary.
473    set tmpfile file[pid].vtk
474    set f [open "$tmpfile" "w"]
475    fconfigure $f -translation binary -encoding binary
476    puts $f $out
477    close $f
478
479    # Read the data back into a vtk dataset and query the bounds.
480    set reader $this-datasetreader
481    vtkDataSetReader $reader
482    $reader SetFileName $tmpfile
483    $reader Update
484    set output [$reader GetOutput]
485    set _numPoints [$output GetNumberOfPoints]
486    set _numCells [$output GetNumberOfCells]
487    foreach { xmin xmax ymin ymax zmin zmax } [$output GetBounds] break
488    set _limits(x) [list $xmin $xmax]
489    set _limits(y) [list $ymin $ymax]
490    set _limits(z) [list $zmin $zmax]
491    file delete $tmpfile
492    rename $output ""
493    rename $reader ""
494    return 1
495}
496
497itcl::body Rappture::Mesh::ReadGrid { path } {
498    set _type "grid"
499
500    if { ![GetDimension $path] } {
501        return 0
502    }
503    set numUniform 0
504    set numRectilinear 0
505    set numCurvilinear 0
506    foreach axis { x y z } {
507        set min    [$_xmlobj get "$path.grid.${axis}axis.min"]
508        set max    [$_xmlobj get "$path.grid.${axis}axis.max"]
509        set num    [$_xmlobj get "$path.grid.${axis}axis.numpoints"]
510        set coords [$_xmlobj get "$path.grid.${axis}coords"]
511        set dim    [$_xmlobj get "$path.grid.${axis}dim"]
512        if { $min != "" && $max != "" && $num != "" && $num > 0 } {
513            set ${axis}Min $min
514            set ${axis}Max $max
515            set ${axis}Num $num
516            incr numUniform
517        } elseif { $coords != "" } {
518            incr numRectilinear
519            set ${axis}Coords $coords
520        } elseif { $dim != "" } {
521            set ${axis}Num $dim
522            incr numCurvilinear
523        }
524    }
525    set _dim [expr $numRectilinear + $numUniform + $numCurvilinear]
526    if { $_dim == 0 } {
527        # No data found.
528        puts stderr "WARNING: bad grid \"$path\": no data found"
529        return 0
530    }
531    if { $numCurvilinear > 0 } {
532        # This is the 2D/3D curilinear case. We found <xdim>, <ydim>, or <zdim>
533        if { $numRectilinear > 0 || $numUniform > 0 } {
534            puts stderr "WARNING: bad grid \"$path\": can't mix curvilinear and rectilinear grids."
535            return 0
536        }
537        set points [$_xmlobj get $path.grid.points]
538        if { $points == "" } {
539            puts stderr "WARNING: bad grid \"$path\": no <points> found."
540            return 0
541        }
542        if { ![info exists xNum] } {
543            puts stderr "WARNING: bad grid \"$path\": invalid dimensions for curvilinear grid: missing <xdim> from grid description."
544            return 0
545        }
546        set all [blt::vector create \#auto]
547        set xv [blt::vector create \#auto]
548        set yv [blt::vector create \#auto]
549        set zv [blt::vector create \#auto]
550        $all set $points
551        set numCoords [$all length]
552        if { [info exists zNum] } {
553            set _dim 3
554            set _numPoints [expr $xNum * $yNum * $zNum]
555            set _numCells [expr ($xNum - 1) * ($yNum - 1) * ($zNum - 1)]
556            if { ($_numPoints*3) != $numCoords } {
557                puts stderr "WARNING: bad grid \"$path\": invalid grid: \# of points does not match dimensions <xdim> * <ydim> * <zdim>"
558                return 0
559            }
560            if { ($numCoords % 3) != 0 } {
561                puts stderr "WARNING: bad grid \"$path\": wrong \# of coordinates for 3D grid"
562                return 0
563            }
564            $all split $xv $yv $zv
565            foreach axis {x y z} {
566                set vector [set ${axis}v]
567                set _limits($axis) [$vector limits]
568            }
569            append out "DATASET STRUCTURED_GRID\n"
570            append out "DIMENSIONS $xNum $yNum $zNum\n"
571            append out "POINTS $_numPoints double\n"
572            append out [$all range 0 end]
573            append out "\n"
574            set _vtkdata $out
575        } elseif { [info exists yNum] } {
576            set _dim 2
577            set _numPoints [expr $xNum * $yNum]
578            set _numCells [expr ($xNum - 1) * ($yNum - 1)]
579            if { ($_numPoints*2) != $numCoords } {
580                puts stderr "WARNING: bad grid \"$path\": \# of points does not match dimensions <xdim> * <ydim>"
581                return 0
582            }
583            if { ($numCoords % 2) != 0 } {
584                puts stderr "WARNING: bad grid \"$path\": wrong \# of coordinates for 2D grid"
585                return 0
586            }
587            foreach axis {x y} {
588                set vector [set ${axis}v]
589                set _limits($axis) [$vector limits]
590            }
591            set _limits(z) [list 0 0]
592            $zv seq 0 0 [$xv length]
593            $all merge $xv $yv $zv
594            append out "DATASET STRUCTURED_GRID\n"
595            append out "DIMENSIONS $xNum $yNum 1\n"
596            append out "POINTS $_numPoints double\n"
597            append out [$all range 0 end]
598            append out "\n"
599            set _vtkdata $out
600        } else {
601            set _dim 1
602            set _numPoints $xNum
603            set _numCells [expr $xNum - 1]
604            if { $_numPoints != $numCoords } {
605                puts stderr "WARNING: bad grid \"$path\": \# of points does not match <xdim>"
606                return 0
607            }
608            set _limits(x) [$xv limits]
609            set _limits(y) [list 0 0]
610            set _limits(z) [list 0 0]
611            $yv seq 0 0 [$xv length]
612            $zv seq 0 0 [$xv length]
613            $all merge $xv $yv $zv
614            append out "DATASET STRUCTURED_GRID\n"
615            append out "DIMENSIONS $xNum 1 1\n"
616            append out "POINTS $_numPoints double\n"
617            append out [$all range 0 end]
618            append out "\n"
619            set _vtkdata $out
620        }
621        blt::vector destroy $all $xv $yv $zv
622        return 1
623    }
624    if { $numRectilinear == 0 && $numUniform > 0} {
625        # This is the special case where all axes 2D/3D are uniform. 
626        # This results in a STRUCTURED_POINTS
627        if { $_dim == 1 } {
628            set xSpace [expr ($xMax - $xMin) / double($xNum - 1)]
629            set _numPoints $xNum
630            set _numCells [expr $xNum - 1]
631            append out "DATASET STRUCTURED_POINTS\n"
632            append out "DIMENSIONS $xNum 1 1\n"
633            append out "ORIGIN $xMin 0 0\n"
634            append out "SPACING $xSpace 0 0\n"
635            set _vtkdata $out
636            set _limits(x) [list $xMin $xMax]
637            set _limits(y) [list 0 0]
638            set _limits(z) [list 0 0]
639        } elseif { $_dim == 2 } {
640            set xSpace [expr ($xMax - $xMin) / double($xNum - 1)]
641            set ySpace [expr ($yMax - $yMin) / double($yNum - 1)]
642            set _numPoints [expr $xNum * $yNum]
643            set _numCells [expr ($xNum - 1) * ($yNum - 1)]
644            append out "DATASET STRUCTURED_POINTS\n"
645            append out "DIMENSIONS $xNum $yNum 1\n"
646            append out "ORIGIN $xMin $yMin 0\n"
647            append out "SPACING $xSpace $ySpace 0\n"
648            set _vtkdata $out
649            foreach axis {x y} {
650                set _limits($axis) [list [set ${axis}Min] [set ${axis}Max]]
651            }
652            set _limits(z) [list 0 0]
653        } elseif { $_dim == 3 } {
654            set xSpace [expr ($xMax - $xMin) / double($xNum - 1)]
655            set ySpace [expr ($yMax - $yMin) / double($yNum - 1)]
656            set zSpace [expr ($zMax - $zMin) / double($zNum - 1)]
657            set _numPoints [expr $xNum * $yNum * $zNum]
658            set _numCells [expr ($xNum - 1) * ($yNum - 1) * ($zNum - 1)]
659            append out "DATASET STRUCTURED_POINTS\n"
660            append out "DIMENSIONS $xNum $yNum $zNum\n"
661            append out "ORIGIN $xMin $yMin $zMin\n"
662            append out "SPACING $xSpace $ySpace $zSpace\n"
663            set _vtkdata $out
664            foreach axis {x y z} {
665                set _limits($axis) [list [set ${axis}Min] [set ${axis}Max]]
666            }
667        } else {
668            puts stderr "WARNING: bad grid \"$path\": bad dimension \"$_dim\""
669            return 0
670        }
671        return 1
672    }
673    # This is the hybrid case.  Some axes are uniform, others are nonuniform.
674    set xv [blt::vector create \#auto]
675    if { [info exists xMin] } {
676        $xv seq $xMin $xMax $xNum
677    } else {
678        $xv set [$_xmlobj get $path.grid.xcoords]
679        set xMin [$xv min]
680        set xMax [$xv max]
681        set xNum [$xv length]
682    }
683    set yv [blt::vector create \#auto]
684    if { $_dim > 1 } {
685        if { [info exists yMin] } {
686            $yv seq $yMin $yMax $yNum
687        } else {
688            $yv set [$_xmlobj get $path.grid.ycoords]
689            set yMin [$yv min]
690            set yMax [$yv max]
691            set yNum [$yv length]
692        }
693    } else {
694        set yNum 1
695    }
696    set zv [blt::vector create \#auto]
697    if { $_dim == 3 } {
698        if { [info exists zMin] } {
699            $zv seq $zMin $zMax $zNum
700        }  else {
701            $zv set [$_xmlobj get $path.grid.zcoords]
702            set zMin [$zv min]
703            set zMax [$zv max]
704            set zNum [$zv length]
705        }
706    } else {
707        set zNum 1
708    }
709    if { $_dim == 3 } {
710        set _numPoints [expr $xNum * $yNum * $zNum]
711        set _numCells [expr ($xNum - 1) * ($yNum - 1) * ($zNum - 1)]
712        append out "DATASET RECTILINEAR_GRID\n"
713        append out "DIMENSIONS $xNum $yNum $zNum\n"
714        append out "X_COORDINATES $xNum double\n"
715        append out [$xv range 0 end]
716        append out "\n"
717        append out "Y_COORDINATES $yNum double\n"
718        append out [$yv range 0 end]
719        append out "\n"
720        append out "Z_COORDINATES $zNum double\n"
721        append out [$zv range 0 end]
722        append out "\n"
723        set _vtkdata $out
724        foreach axis {x y z} {
725            if { [info exists ${axis}Min] } {
726                set _limits($axis) [list [set ${axis}Min] [set ${axis}Max]]
727            }
728        }
729    } elseif { $_dim == 2 } {
730        set _numPoints [expr $xNum * $yNum]
731        set _numCells [expr ($xNum - 1) * ($yNum - 1)]
732        append out "DATASET RECTILINEAR_GRID\n"
733        append out "DIMENSIONS $xNum $yNum 1\n"
734        append out "X_COORDINATES $xNum double\n"
735        append out [$xv range 0 end]
736        append out "\n"
737        append out "Y_COORDINATES $yNum double\n"
738        append out [$yv range 0 end]
739        append out "\n"
740        append out "Z_COORDINATES 1 double\n"
741        append out "0\n"
742        foreach axis {x y} {
743            if { [info exists ${axis}Min] } {
744                set _limits($axis) [list [set ${axis}Min] [set ${axis}Max]]
745            }
746        }
747        set _limits(z) [list 0 0]
748        set _vtkdata $out
749    } elseif { $_dim == 1 } {
750        set _numPoints $xNum
751        set _numCells [expr $xNum - 1]
752        append out "DATASET RECTILINEAR_GRID\n"
753        append out "DIMENSIONS $xNum 1 1\n"
754        append out "X_COORDINATES $xNum double\n"
755        append out [$xv range 0 end]
756        append out "\n"
757        append out "Y_COORDINATES 1 double\n"
758        append out "0\n"
759        append out "Z_COORDINATES 1 double\n"
760        append out "0\n"
761        if { [info exists xMin] } {
762            set _limits(x) [list $xMin $xMax]
763        }
764        set _limits(y) [list 0 0]
765        set _limits(z) [list 0 0]
766        set _vtkdata $out
767    } else {
768        puts stderr "WARNING: bad grid \"$path\": invalid dimension \"$_dim\""
769        return 0
770    }
771    blt::vector destroy $xv $yv $zv
772    return 1
773}
774
775itcl::body Rappture::Mesh::WritePointCloud { path xv yv zv } {
776    set _type "cloud"
777    set _numPoints [$xv length]
778    append out "DATASET POLYDATA\n"
779    append out "POINTS $_numPoints double\n"
780    foreach x [$xv range 0 end] y [$yv range 0 end] z [$zv range 0 end] {
781        append out "$x $y $z\n"
782    }
783    set _vtkdata $out
784    set _limits(x) [$xv limits]
785    if { $_dim > 1 } {
786        set _limits(y) [$yv limits]
787    } else {
788        set _limits(y) [list 0 0]
789    }
790    if { $_dim == 3 } {
791        set _limits(z) [$zv limits]
792    } else {
793        set _limits(z) [list 0 0]
794    }
795    return 1
796}
797
798itcl::body Rappture::Mesh::WriteTriangles { path xv yv zv triangles } {
799    set _type "triangles"
800    set _numPoints [$xv length]
801    set _numCells 0
802    set data {}
803    set celltypes {}
804    foreach { a b c } $triangles {
805        append data " 3 $a $b $c\n"
806        append celltypes "5\n"
807        incr _numCells
808    }
809    append out "DATASET UNSTRUCTURED_GRID\n"
810    append out "POINTS $_numPoints double\n"
811    foreach x [$xv range 0 end] y [$yv range 0 end] z [$zv range 0 end] {
812        append out " $x $y $z\n"
813    }
814    set count [expr $_numCells * 4]
815    append out "CELLS $_numCells $count\n"
816    append out $data
817    append out "CELL_TYPES $_numCells\n"
818    append out $celltypes
819    set _limits(x) [$xv limits]
820    set _limits(y) [$yv limits]
821    if { $_dim == 3 } {
822        set _limits(z) [$zv limits]
823    } else {
824        set _limits(z) [list 0 0]
825    }
826    set _vtkdata $out
827    return 1
828}
829
830itcl::body Rappture::Mesh::WriteQuads { path xv yv zv quads } {
831    set _type "quads"
832    set _numPoints [$xv length]
833    set _numCells 0
834    set data {}
835    set celltypes {}
836    foreach { a b c d } $quads {
837        append data " 4 $a $b $c $d\n"
838        append celltypes "9\n"
839        incr _numCells
840    }
841    append out "DATASET UNSTRUCTURED_GRID\n"
842    append out "POINTS $_numPoints double\n"
843    foreach x [$xv range 0 end] y [$yv range 0 end] z [$zv range 0 end] {
844        append out " $x $y $z\n"
845    }
846    set count [expr $_numCells * 5]
847    append out "CELLS $_numCells $count\n"
848    append out $data
849    append out "CELL_TYPES $_numCells\n"
850    append out $celltypes
851    set _limits(x) [$xv limits]
852    set _limits(y) [$yv limits]
853    if { $_dim == 3 } {
854        set _limits(z) [$zv limits]
855    } else {
856        set _limits(z) [list 0 0]
857    }
858    set _vtkdata $out
859    return 1
860}
861
862itcl::body Rappture::Mesh::WriteTetrahedrons { path xv yv zv tetras } {
863    set _type "tetrahedrons"
864    set _numPoints [$xv length]
865    set _numCells 0
866    set data {}
867    set celltypes {}
868    foreach { a b c d } $tetras {
869        append data " 4 $a $b $c $d\n"
870        append celltypes "10\n"
871        incr _numCells
872    }
873    append out "DATASET UNSTRUCTURED_GRID\n"
874    append out "POINTS $_numPoints double\n"
875    foreach x [$xv range 0 end] y [$yv range 0 end] z [$zv range 0 end] {
876        append out " $x $y $z\n"
877    }
878    set count [expr $_numCells * 5]
879    append out "CELLS $_numCells $count\n"
880    append out $data
881    append out "CELL_TYPES $_numCells\n"
882    append out $celltypes
883    set _limits(x) [$xv limits]
884    set _limits(y) [$yv limits]
885    set _limits(z) [$zv limits]
886
887    set _vtkdata $out
888    return 1
889}
890
891itcl::body Rappture::Mesh::WriteHexahedrons { path xv yv zv hexas } {
892    set _type "hexahedrons"
893    set _numPoints [$xv length]
894    set _numCells 0
895    set data {}
896    set celltypes {}
897    foreach { a b c d e f g h } $hexas {
898        append data " 8 $a $b $c $d $e $f $g $h\n"
899        append celltypes "12\n"
900        incr _numCells
901    }
902    append out "DATASET UNSTRUCTURED_GRID\n"
903    append out "POINTS $_numPoints double\n"
904    foreach x [$xv range 0 end] y [$yv range 0 end] z [$zv range 0 end] {
905        append out " $x $y $z\n"
906    }
907    set count [expr $_numCells * 9]
908    append out "CELLS $_numCells $count\n"
909    append out $data
910    append out "CELL_TYPES $_numCells\n"
911    append out $celltypes
912    set _limits(x) [$xv limits]
913    set _limits(y) [$yv limits]
914    set _limits(z) [$zv limits]
915
916    set _vtkdata $out
917    return 1
918}
919
920itcl::body Rappture::Mesh::WriteWedges { path xv yv zv wedges } {
921    set _type "wedges"
922    set _numPoints [$xv length]
923    set _numCells 0
924    set data {}
925    set celltypes {}
926    foreach { a b c d e f } $wedges {
927        append data " 6 $a $b $c $d $e $f\n"
928        append celltypes "13\n"
929        incr _numCells
930    }
931    append out "DATASET UNSTRUCTURED_GRID\n"
932    append out "POINTS $_numPoints double\n"
933    foreach x [$xv range 0 end] y [$yv range 0 end] z [$zv range 0 end] {
934        append out " $x $y $z\n"
935    }
936    set count [expr $_numCells * 7]
937    append out "CELLS $_numCells $count\n"
938    append out $data
939    append out "CELL_TYPES $_numCells\n"
940    append out $celltypes
941    set _limits(x) [$xv limits]
942    set _limits(y) [$yv limits]
943    set _limits(z) [$zv limits]
944
945    set _vtkdata $out
946    return 1
947}
948
949itcl::body Rappture::Mesh::WritePyramids { path xv yv zv pyramids } {
950    set _type "pyramids"
951    set _numPoints [$xv length]
952    set _numCells 0
953    set data {}
954    set celltypes {}
955    foreach { a b c d e } $pyramids {
956        append data " 5 $a $b $c $d $e\n"
957        append celltypes "14\n"
958        incr _numCells
959    }
960    append out "DATASET UNSTRUCTURED_GRID\n"
961    append out "POINTS $_numPoints double\n"
962    foreach x [$xv range 0 end] y [$yv range 0 end] z [$zv range 0 end] {
963        append out " $x $y $z\n"
964    }
965    set count [expr $_numCells * 6]
966    append out "CELLS $_numCells $count\n"
967    append out $data
968    append out "CELL_TYPES $_numCells\n"
969    append out $celltypes
970    set _limits(x) [$xv limits]
971    set _limits(y) [$yv limits]
972    set _limits(z) [$zv limits]
973
974    set _vtkdata $out
975    return 1
976}
977
978itcl::body Rappture::Mesh::WriteHybridCells { path xv yv zv cells celltypes } {
979    set _type "unstructured"
980    set lines [split $cells \n]
981    set numCellTypes [llength $celltypes]
982    if { $numCellTypes == 1 } {
983        set celltype [GetCellType $celltypes]
984    }
985
986    set _numPoints [$xv length]
987    set data {}
988    set count 0
989    set _numCells 0
990    set celltypes {}
991    foreach line $lines {
992        set length [llength $line]
993        if { $length == 0 } {
994            continue
995        }
996        if { $numCellTypes > 1 } {
997            set cellType [GetCellType [lindex $cellTypes $_numCells]]
998        }
999        set numIndices [GetNumIndices $celltype]
1000        if { $numIndices > 0 && $numIndices != $length } {
1001            puts stderr "WARNING: bad unstructured grid \"$path\": wrong \# of indices specified for celltype $celltype on line \"$line\""
1002            return 0
1003        } else {
1004            set numIndices $length
1005        }
1006        append data " $numIndices $line\n"
1007        lappend celltypes $celltype
1008        incr count $length;         # Include the indices
1009        incr count;                 # and the number of indices
1010        incr _numCells
1011    }
1012    append out "DATASET UNSTRUCTURED_GRID\n"
1013    append out "POINTS $_numPoints double\n"
1014    set all [blt::vector create \#auto]
1015    $all merge $xv $yv $zv
1016    append out [$all range 0 end]
1017    blt::vector destroy $all
1018    append out "CELLS $_numCells $count\n"
1019    append out $data
1020    append out "CELL_TYPES $_numCells\n"
1021    append out $celltypes
1022    set _limits(x) [$xv limits]
1023    set _limits(y) [$yv limits]
1024    if { $_dim < 3 } {
1025        set _limits(z) [list 0 0]
1026    } else {
1027        set _limits(z) [$zv limits]
1028    }
1029
1030    set _vtkdata $out
1031    return 1
1032}
1033
1034itcl::body Rappture::Mesh::ReadUnstructuredGrid { path } {
1035    set _type "unstructured"
1036
1037    if { ![GetDimension $path] } {
1038        return 0
1039    }
1040    # Step 1: Verify that there's only one cell tag of any kind.
1041    set numCells 0
1042    foreach type { cells triangles quads tetrahedrons
1043        hexahedrons wedges pyramids } {
1044        set data [$_xmlobj get $path.unstructured.$type]
1045        if { $data != "" } {
1046            incr numCells
1047        }
1048    }
1049    # The generic <cells> tag requires there be a <celltypes> tag too.
1050    set celltypes [$_xmlobj get $path.unstructured.celltypes]
1051    if { $numCells == 0 && $celltypes != "" } {
1052        puts stderr "WARNING: bad unstuctured grid \"$path\": no <cells> description found."
1053        return 0
1054    }
1055    if { $numCells > 1 } {
1056        puts stderr "WARNING: bad unstructured grid \"$path\": too many <cells>, <triangles>, <quads>... descriptions found: should be only one."
1057        return 0
1058    }
1059    foreach type { cells triangles quads tetrahedrons
1060        hexahedrons wedges pyramids } {
1061        set data [$_xmlobj get $path.unstructured.$type]
1062        if { $data != "" } {
1063            break
1064        }
1065    }
1066    # Step 2: Allow points to be specified as <points> or
1067    #         <xcoords>, <ycoords>, <zcoords>.  Split and convert into
1068    #         3 vectors, one for each coordinate.
1069    if { $_dim == 1 } {
1070        set xcoords [$_xmlobj get $path.unstructured.xcoords]
1071        set ycoords [$_xmlobj get $path.unstructured.ycoords]
1072        set zcoords [$_xmlobj get $path.unstructured.zcoords]
1073        set data    [$_xmlobj get $path.unstructured.points]
1074        if { $ycoords != "" } {
1075            put stderr "can't specify <ycoords> with a 1D mesh"
1076            return 0
1077        }
1078        if { $zcoords != "" } {
1079            put stderr "can't specify <zcoords> with a 1D mesh"
1080            return 0
1081        }
1082        if { $xcoords != "" } {
1083            set xv [blt::vector create \#auto]
1084            $xv set $xcoords
1085        } elseif { $data != "" } {
1086            Rappture::ReadPoints $data dim points
1087            if { $points == "" } {
1088                puts stderr "WARNING: bad unstructured grid \"$path\": no <points> found."
1089                return 0
1090            }
1091            if { $dim != 1 } {
1092                puts stderr "WARNING: bad unstructured grid \"$path\": \# of coordinates per point is \"$dim\": does not agree with dimension specified for mesh \"$_dim\""
1093                return 0
1094            }
1095            set xv [blt::vector create \#auto]
1096            $xv set $points
1097        } else {
1098            puts stderr "WARNING: bad unstructured grid \"$path\": no points specified."
1099            return 0
1100        }
1101        set yv [blt::vector create \#auto]
1102        set zv [blt::vector create \#auto]
1103        $yv seq 0 0 [$xv length];       # Make an all zeroes vector.
1104        $zv seq 0 0 [$xv length];       # Make an all zeroes vector.
1105    } elseif { $_dim == 2 } {
1106        set xcoords [$_xmlobj get $path.unstructured.xcoords]
1107        set ycoords [$_xmlobj get $path.unstructured.ycoords]
1108        set zcoords [$_xmlobj get $path.unstructured.zcoords]
1109        set data    [$_xmlobj get $path.unstructured.points]
1110        if { $zcoords != "" } {
1111            put stderr "can't specify <zcoords> with a 2D mesh"
1112            return 0
1113        }
1114        if { $xcoords != "" && $ycoords != "" } {
1115            set xv [blt::vector create \#auto]
1116            set yv [blt::vector create \#auto]
1117            $xv set $xcoords
1118            $yv set $ycoords
1119        } elseif { $data != "" } {
1120            Rappture::ReadPoints $data dim points
1121            if { $points == "" } {
1122                puts stderr "WARNING: bad unstructured grid \"$path\": no <points> found."
1123                return 0
1124            }
1125            if { $dim != 2 } {
1126                puts stderr "WARNING: bad unstructured grid \"$path\": \# of coordinates per point is \"$dim\": does not agree with dimension specified for mesh \"$_dim\""
1127                return 0
1128            }
1129            set all [blt::vector create \#auto]
1130            set xv [blt::vector create \#auto]
1131            set yv [blt::vector create \#auto]
1132            $all set $points
1133            $all split $xv $yv
1134            blt::vector destroy $all
1135        } else {
1136            puts stderr "WARNING: bad unstructured grid \"$path\": no points specified."
1137            return 0
1138        }
1139        set zv [blt::vector create \#auto]
1140        $zv seq 0 0 [$xv length];       # Make an all zeroes vector.
1141    } elseif { $_dim == 3 } {
1142        set xcoords [$_xmlobj get $path.unstructured.xcoords]
1143        set ycoords [$_xmlobj get $path.unstructured.ycoords]
1144        set zcoords [$_xmlobj get $path.unstructured.zcoords]
1145        set data    [$_xmlobj get $path.unstructured.points]
1146        if { $xcoords != "" && $ycoords != "" && $zcoords != "" } {
1147            set xv [blt::vector create \#auto]
1148            set yv [blt::vector create \#auto]
1149            set zv [blt::vector create \#auto]
1150            $xv set $xcoords
1151            $yv set $ycoords
1152            $zv set $zcoords
1153        } elseif { $data != "" }  {
1154            Rappture::ReadPoints $data dim points
1155            if { $points == "" } {
1156                puts stderr "WARNING: bad unstructured grid \"$path\": no <points> found."
1157                return 0
1158            }
1159            if { $dim != 3 } {
1160                puts stderr "WARNING: bad unstructured grid \"$path\": \# of coordinates per point is \"$dim\": does not agree with dimension specified for mesh \"$_dim\""
1161                return 0
1162            }
1163            set all [blt::vector create \#auto]
1164            set xv [blt::vector create \#auto]
1165            set yv [blt::vector create \#auto]
1166            set zv [blt::vector create \#auto]
1167            $all set $points
1168            $all split $xv $yv $zv
1169            blt::vector destroy $all
1170        } else {
1171            puts stderr "WARNING: bad unstructured grid \"$path\": no points specified."
1172            return 0
1173        }
1174    }
1175    set _numPoints [$xv length]
1176
1177    # Step 3: Write the points and cells as vtk data.
1178    if { $numCells == 0 } {
1179        set result [WritePointCloud $path $xv $yv $zv]
1180    } elseif { $type == "cells" } {
1181        set cells [$_xmlobj get $path.unstructured.cells]
1182        if { $cells == "" } {
1183            puts stderr "WARNING: bad unstructured grid \"$path\": no cells found."
1184            return 0
1185        }
1186        set result [WriteHybridCells $path $xv $yv $zv $cells $celltypes]
1187    } else {
1188        set cmd "Write[string totitle $type]"
1189        set cells [$_xmlobj get $path.unstructured.$type]
1190        if { $cells == "" } {
1191            puts stderr "WARNING: bad unstructured grid \"$path\": no cells found."
1192            return 0
1193        }
1194        set result [$cmd $path $xv $yv $zv $cells]
1195    }
1196    # Clean up.
1197    blt::vector destroy $xv $yv $zv
1198    return $result
1199}
1200
1201# ----------------------------------------------------------------------
1202# USAGE: ReadNodesElements <path>
1203#
1204# Used internally to build a mesh representation based on nodes and
1205# elements stored in the XML.
1206# ----------------------------------------------------------------------
1207itcl::body Rappture::Mesh::ReadNodesElements {path} {
1208    set _type "nodeselements"
1209    set count 0
1210
1211    # Scan for nodes.  Each node represents a vertex.
1212    set data {}
1213    foreach cname [$_xmlobj children -type node $path] {
1214        append data "[$_xmlobj get $path.$cname]\n"
1215    }   
1216    Rappture::ReadPoints $data _dim points
1217    if { $_dim == 2 } {
1218        set all [blt::vector create \#auto]
1219        set xv [blt::vector create \#auto]
1220        set yv [blt::vector create \#auto]
1221        set zv [blt::vector create \#auto]
1222        $all set $points
1223        $all split $xv $yv
1224        set _numPoints [$xv length]
1225        set _limits(x) [$xv limits]
1226        set _limits(y) [$yv limits]
1227        set _limits(z) [list 0 0]
1228        # 2D Dataset. All Z coordinates are 0
1229        $zv seq 0.0 0.0 $_numPoints
1230        $all merge $xv $yv $zv
1231        set points [$all range 0 end]
1232        blt::vector destroy $all $xv $yv $zv
1233    } elseif { $_dim == 3 } {
1234        set all [blt::vector create \#auto]
1235        set xv [blt::vector create \#auto]
1236        set yv [blt::vector create \#auto]
1237        set zv [blt::vector create \#auto]
1238        $all set $points
1239        $all split $xv $yv $zv
1240        set _numPoints [$xv length]
1241        set _limits(x) [$xv limits]
1242        set _limits(y) [$yv limits]
1243        set _limits(z) [$zv limits]
1244        set points [$all range 0 end]
1245        blt::vector destroy $all $xv $yv $zv
1246    } else {
1247        error "bad dimension \"$_dim\" for nodes mesh"
1248    }
1249    array set node2celltype {
1250        3 5
1251        4 10
1252        8 12
1253        6 13
1254        5 14
1255    }
1256    set count 0
1257    set _numCells 0
1258    set celltypes {}
1259    set data {}
1260    # Next scan for elements.  Each element represents a cell.
1261    foreach cname [$_xmlobj children -type element $path] {
1262        set nodeList [$_mesh get $cname.nodes]
1263        set numNodes [llength $nodeList]
1264        if { ![info exists node2celltype($numNodes)] } {
1265            puts stderr "WARNING: bad nodes/elements mesh \$path\": unknown number of indices \"$_numNodes\": should be 3, 4, 5, 6, or 8"
1266            return 0
1267        }
1268        set celltype $node2celltype($numNodes)
1269        append celltypes "  $celltype\n"
1270        if { $celltype == 12 } {
1271            # Formerly used voxels instead of hexahedrons. We're converting
1272            # it here to be backward compatible and still fault-tolerant of
1273            # non-axis aligned cells.
1274            set newList {}
1275            foreach i { 0 1 3 2 4 5 7 6 } {
1276                lappend newList [lindex $nodeList $i]
1277            }
1278            set nodeList $newList
1279        }
1280        append data "  $numNodes $nodeList\n"
1281        incr _numCells
1282        incr count $numNodes
1283        incr count;                     # One extra for the VTK celltype id.
1284    }
1285
1286    append out "DATASET UNSTRUCTURED_GRID\n"
1287    append out "POINTS $_numPoints double\n"
1288    append out $points
1289    append out "\n"
1290    append out "CELLS $_numCells $count\n"
1291    append out $data
1292    append out "CELL_TYPES $_numCells\n"
1293    append out $celltypes
1294    append out "\n"
1295    set _vtkdata $out
1296    set _isValid 1
1297}
1298
1299itcl::body Rappture::Mesh::GetCellType { name } {
1300    array set name2type {
1301        "triangle"     5
1302        "quad"         9
1303        "tetrahedron"  10
1304        "hexahedron"   12
1305        "wedge"        13
1306        "pyramid"      14
1307    }
1308    if { [info exists name2type($name)] } {
1309        return $name2type($name)
1310    }
1311    if { [string is int $name] == 1 && $name >= 1 && $name <= 16 } {
1312        return $name
1313    }
1314    error "invalid celltype \"$name\""
1315}
1316
1317itcl::body Rappture::Mesh::GetNumIndices { type } {
1318    array set type2indices {
1319        1       1
1320        2       0
1321        3       2
1322        4       0
1323        5       3
1324        6       0
1325        7       0
1326        8       4
1327        9       4
1328        10      4
1329        11      8
1330        12      8
1331        13      6
1332        14      5
1333        15      0
1334        16      0
1335    }
1336    if { [info exists type2indices($type)] } {
1337        return $type2indices($type)
1338    }
1339    error "invalid celltype \"$type\""
1340}
Note: See TracBrowser for help on using the repository browser.