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

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

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