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

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

new mesh changes (version 1,999,999), untested quads, tetrahedrons, hexahedrons, wedges, and pyramids
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 ReadTriangles { xv yv zv triangles }
77    private method ReadQuads { xv yv zv quads }
78    private method ReadTetrahedrons { xv yv zv tetrahedrons }
79    private method ReadHexahedrons { xv yv zv hexhedrons }
80    private method ReadWedges { xv yv zv wedges }
81    private method ReadPyramids { xv yv zv pyramids }
82    private method ReadCells { xv yv zv cells celltypes }
83    private method ReadCloud { 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::ReadCloud { 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::ReadTriangles { 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::ReadQuads { 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::ReadTetrahedrons { 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::ReadHexahedrons { 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::ReadWedges { 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::ReadPyramids { 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 " 6 $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 * 5]\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::ReadCells { 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.
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 <zcoord> with a 2 dimensional 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 <cells> mesh"
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: Read the cells and write the vtk data.
1037    if { $numCells == 0 } {
1038        puts stderr "numCells=$numCells call ReadCloud"
1039        ReadCloud $xv $yv $zv
1040    } elseif { $type == "cells" } {
1041        set cells [$_xmlobj get $path.unstructured.cells]
1042        ReadCells $xv $yv $zv $cells $celltypes
1043    } else {
1044        puts stderr "type=$type"
1045        set cmd "Read[string totitle $type]"
1046        set cells [$_xmlobj get $path.unstructured.$type]
1047        $cmd $xv $yv $zv $cells
1048    }
1049    # Clean up.
1050    blt::vector destroy $xv $yv $zv
1051}
1052
1053# ----------------------------------------------------------------------
1054# USAGE: ReadNodesElements <path>
1055#
1056# Used internally to build a mesh representation based on nodes and
1057# elements stored in the XML.
1058# ----------------------------------------------------------------------
1059itcl::body Rappture::Mesh::ReadNodesElements {path} {
1060    set type "nodeselements"
1061    set count 0
1062
1063    # Scan for nodes.  Each node represents a vertex.
1064    set data {}
1065    foreach cname [$_xmlobj children -type node $path] {
1066        append data "[$_xmlobj get $path.$cname]\n"
1067    }   
1068    Rappture::ReadPoints $data _dim points
1069    if { $_dim == 2 } {
1070        set all [blt::vector create \#auto]
1071        set xv [blt::vector create \#auto]
1072        set yv [blt::vector create \#auto]
1073        set zv [blt::vector create \#auto]
1074        $all set $points
1075        $all split $xv $yv
1076        set _numPoints [$xv length]
1077        set _limits(x) [$xv limits]
1078        set _limits(y) [$yv limits]
1079        # 2D Dataset. All Z coordinates are 0
1080        $zv seq 0.0 0.0 $_numPoints
1081        $all merge $xv $yv $zv
1082        set points [$all range 0 end]
1083        blt::vector destroy $all $xv $yv $zv
1084    } elseif { $_dim == 3 } {
1085        set all [blt::vector create \#auto]
1086        set xv [blt::vector create \#auto]
1087        set yv [blt::vector create \#auto]
1088        set zv [blt::vector create \#auto]
1089        $all set $points
1090        $all split $xv $yv $zv
1091        set _numPoints [$xv length]
1092        set _limits(x) [$xv limits]
1093        set _limits(y) [$yv limits]
1094        set _limits(z) [$zv limits]
1095        set points [$all range 0 end]
1096        blt::vector destroy $all $xv $yv $zv
1097    } else {
1098        error "bad dimension \"$_dim\" for nodes mesh"
1099    }
1100    array set node2celltype {
1101        3 5
1102        4 10
1103        8 11
1104        6 13
1105        5 14
1106    }
1107    set count 0
1108    set _numCells 0
1109    set celltypes {}
1110    set data {}
1111    # Next scan for elements.  Each element represents a cell.
1112    foreach cname [$_xmlobj children -type element $path] {
1113        set nodeList [$_mesh get $cname.nodes]
1114        set numNodes [llength $nodeList]
1115        if { ![info exists node2celltype($numNodes)] } {
1116            puts stderr "unknown number of indices \"$_numNodes\": should be 3, 4, 5, 6, or 8"
1117            return 0
1118        }
1119        set celltype $node2celltype($numNodes)
1120        append celltypes "  $celltype\n"
1121        append data "  $numNodes $nodeList\n"
1122        incr _numCells
1123        incr count $numNodes
1124        incr count;                     # One extra for the VTK celltype id.
1125    }
1126
1127    append out "DATASET UNSTRUCTURED_GRID\n"
1128    append out "POINTS $_numPoints double\n"
1129    append out $points
1130    append out "CELLS $_numCells $count\n"
1131    append out $data
1132    append out "CELL_TYPES $_numCells\n"
1133    append out $celltypes
1134    append out "\n"
1135    set _vtkdata $out
1136    set _isValid 1
1137}
1138
1139itcl::body Rappture::Mesh::GetCellType { name } {
1140    array set name2type {
1141        "triangle"     5
1142        "quad"         9
1143        "tetrahedron"  10
1144        "hexahedron"   11
1145        "wedge"        13
1146        "pyramid"      14
1147    }
1148    if { [info exists name2type($name)] } {
1149        return $name2type($name)
1150    }
1151    if { [string is int $name] == 1 && $name >= 1 && $name <= 16 } {
1152        return $name
1153    }
1154    error "invalid celltype \"$name\""
1155}
1156
1157itcl::body Rappture::Mesh::GetNumIndices { type } {
1158    array set type2indices {
1159        1       1
1160        2       0
1161        3       2
1162        4       0
1163        5       3
1164        6       0
1165        7       0
1166        8       4
1167        9       4
1168        10      4
1169        11      8
1170        12      8
1171        13      6
1172        14      5
1173        15      0
1174        16      0
1175    }
1176    if { [info exists type2indices($type)] } {
1177        return $type2indices($type)
1178    }
1179    error "invalid celltype \"$type\""
1180}
Note: See TracBrowser for help on using the repository browser.