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

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

These was are all related to the omenwire example.

o Added validity test for fields, meshes, clouds, and unirect2ds. There is

now a "isvalid" method that viewers should use to verify that the data object
can be plotted.

In some cases with fields this means that the widget won't even be created.
The resultviewer tests for the dimensionality which is by default 0.

o Thanks to Leif for pointing this out, it's not enough to check if the field

is valid. Individual components of the field may be invalid. Added check so
that viewers are never passed the names of invalid field components.

o Changed many "error" commands to just print to stderr and tolerantly deal

with the error.

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