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

Last change on this file since 4259 was 4259, checked in by ldelgass, 7 years ago

Add support for 1D meshes

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