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

Last change on this file since 3544 was 3530, checked in by gah, 11 years ago

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