source: trunk/gui/scripts/mesh.tcl @ 4045

Last change on this file since 4045 was 3996, checked in by gah, 11 years ago

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