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

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

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