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

Last change on this file since 4007 was 3702, checked in by ldelgass, 11 years ago

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