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

Last change on this file since 3443 was 3421, checked in by gah, 12 years ago

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