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

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

merge (by hand) with Rappture1.2 branch
File size: 32.1 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 float\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 float\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        return 0
518    }
519    if { $numCurvilinear > 0 } {
520        # This is the 2D/3D curilinear case. We found <xdim>, <ydim>, or <zdim>
521        if { $numRectilinear > 0 || $numUniform > 0 } {
522            error "can't mix curvilinear and rectilinear grids."
523        }
524        if { $numCurvilinear < 2 } {
525            error "curvilinear grid must be 2D or 3D."
526        }
527        set points [$xmlobj get $path.grid.points]
528        if { $points == "" } {
529            error "missing <points> from grid description."
530        }
531        if { ![info exists xNum] || ![info exists yNum] } {
532            error "invalid dimensions for curvilinear grid: missing <xdim> or <ydim> from grid description."
533        }
534        set all [blt::vector create \#auto]
535        set xv [blt::vector create \#auto]
536        set yv [blt::vector create \#auto]
537        set zv [blt::vector create \#auto]
538        $all set $points
539        set numCoords [$all length]
540        if { [info exists zNum] } {
541            set _dim 3
542            set _numPoints [expr $xNum * $yNum * $zNum]
543            $all set $points
544            if { ($_numPoints*3) != $numCoords } {
545                error "invalid grid: \# of points does not match dimensions <xdim> * <ydim>"
546            }
547            if { ($numCoords % 3) != 0 } {
548                error "wrong \# of coordinates for 3D grid"
549            }
550            $all split $xv $yv $zv
551            foreach axis {x y z} {
552                set vector [set ${axis}v]
553                set _limits($axis) [$vector limits]
554            }
555            append out "DATASET STRUCTURED_GRID\n"
556            append out "DIMENSIONS $xNum $yNum $zNum\n"
557            append out "POINTS $_numPoints double\n"
558            append out [$all range 0 end]
559            append out "\n"
560            set _vtkdata $out
561        } else {
562            set _dim 2
563            set _numPoints [expr $xNum * $yNum]
564            if { ($_numPoints*2) != $numCoords } {
565                error "invalid grid: \# of points does not match dimensions <xdim> * <ydim> * <zdim>"
566            }
567            if { ($numCoords % 2) != 0 } {
568                error "wrong \# of coordinates for 2D grid"
569            }
570            foreach axis {x y} {
571                set vector [set ${axis}v]
572                set _limits($axis) [$vector limits]
573            }
574            $zv seq 0 0 [$xv length]
575            $all merge $xv $yv $zv
576            append out "DATASET STRUCTURED_GRID\n"
577            append out "DIMENSIONS $xNum $yNum 1\n"
578            append out "POINTS $_numPoints double\n"
579            append out [$all range 0 end]
580            append out "\n"
581            set _vtkdata $out
582        }
583        blt::vector destroy $all $xv $yv $zv
584        set _isValid 1
585        return 1
586    }
587    if { $numRectilinear == 0 && $numUniform > 0} {
588        # This is the special case where all axes 2D/3D are uniform. 
589        # This results in a STRUCTURE_POINTS
590        if { $_dim == 2 } {
591            set xSpace [expr ($xMax - $xMin) / double($xNum - 1)]
592            set ySpace [expr ($yMax - $yMin) / double($yNum - 1)]
593            set _numPoints [expr $xNum * $yNum]
594            append out "DATASET STRUCTURED_POINTS\n"
595            append out "DIMENSIONS $xNum $yNum 1\n"
596            append out "SPACING $xSpace $ySpace 0\n"
597            append out "ORIGIN 0 0 0\n"
598            set _vtkdata $out
599            foreach axis {x y} {
600                set _limits($axis) [list [set ${axis}Min] [set ${axis}Max]]
601            }
602        } elseif { $_dim == 3 } {
603            set xSpace [expr ($xMax - $xMin) / double($xNum - 1)]
604            set ySpace [expr ($yMax - $yMin) / double($yNum - 1)]
605            set zSpace [expr ($zMax - $zMin) / double($zNum - 1)]
606            set _numPoints [expr $xNum * $yNum * zNum]
607            append out "DATASET STRUCTURED_POINTS\n"
608            append out "DIMENSIONS $xNum $yNum $zNum\n"
609            append out "SPACING $xSpace $ySpace $zSpace\n"
610            append out "ORIGIN 0 0 0\n"
611            set _vtkdata $out
612            foreach axis {x y z} {
613                set _limits($axis) [list [set ${axis}Min] [set ${axis}Max]]
614            }
615        } else {
616            error "bad dimension of mesh \"$_dim\""
617        }
618        set _isValid 1
619        return 1
620    }
621    # This is the hybrid case.  Some axes are uniform, others are nonuniform.
622    set xv [blt::vector create \#auto]
623    if { [info exists xMin] } {
624        $xv seq $xMin $xMax $xNum
625    } else {
626        $xv set [$xmlobj get $path.grid.xcoords]
627        set xMin [$xv min]
628        set xMax [$xv max]
629        set xNum [$xv length]
630    }
631    set yv [blt::vector create \#auto]
632    if { [info exists yMin] } {
633        $yv seq $yMin $yMax $yNum
634    } else {
635        $yv set [$xmlobj get $path.grid.ycoords]
636        set yMin [$yv min]
637        set yMax [$yv max]
638        set yNum [$yv length]
639    }
640    set zv [blt::vector create \#auto]
641    if { $_dim == 3 } {
642        if { [info exists zMin] } {
643            $zv seq $zMin $zMax $zNum
644        }  else {
645            $zv set [$xmlobj get $path.grid.zcoords]
646            set zMin [$zv min]
647            set zMax [$zv max]
648            set zNum [$zv length]
649        }
650    } else {
651        set zNum 1
652    }
653    if { $_dim == 3 } {
654        set _numPoints [expr $xNum * $yNum * $zNum]
655        append out "DATASET RECTILINEAR_GRID\n"
656        append out "DIMENSIONS $xNum $yNum $zNum\n"
657        append out "X_COORDINATES $xNum double\n"
658        append out [$xv range 0 end]
659        append out "\n"
660        append out "Y_COORDINATES $yNum double\n"
661        append out [$yv range 0 end]
662        append out "\n"
663        append out "Z_COORDINATES $zNum double\n"
664        append out [$zv range 0 end]
665        append out "\n"
666        set _vtkdata $out
667        foreach axis {x y z} {
668            if { [info exists ${axis}Min] } {
669                set _limits($axis) [list [set ${axis}Min] [set ${axis}Max]]
670            }
671        }
672    } elseif { $_dim == 2 } {
673        set _numPoints [expr $xNum * $yNum]
674        append out "DATASET RECTILINEAR_GRID\n"
675        append out "DIMENSIONS $xNum $yNum 1\n"
676        append out "X_COORDINATES $xNum double\n"
677        append out [$xv range 0 end]
678        append out "\n"
679        append out "Y_COORDINATES $yNum double\n"
680        append out [$yv range 0 end]
681        append out "\n"
682        append out "Z_COORDINATES 1 double\n"
683        append out "0\n"
684        foreach axis {x y} {
685            if { [info exists ${axis}Min] } {
686                set _limits($axis) [list [set ${axis}Min] [set ${axis}Max]]
687            }
688        }
689        set _vtkdata $out
690    } else {
691        error "invalid dimension of grid \"$_dim\""
692    }
693    blt::vector destroy $xv $yv $zv
694    set _isValid 1
695    return 1
696}
697
698itcl::body Rappture::Mesh::ReadCells { xmlobj path } {
699    set _type "cells"
700
701    set data [$xmlobj get $path.cells.points]
702    Rappture::ReadPoints $data _dim points
703    if { $points == "" } {
704        puts stderr "no <points> found for <cells> mesh"
705        return 0
706    }
707    array set celltype2vertices {
708        triangles 3
709        quads    4
710        tetraheadrons 4
711        voxels 8
712        hexaheadrons 8
713        wedges 6
714        pyramids 5
715    }
716    array set celltype2vtkid {
717        triangles 5
718        quads    9
719        tetraheadrons 10
720        voxels 11
721        hexaheadrons 12
722        wedges 13
723        pyramids 14
724    }
725    set count 0
726
727    # Try to detect the celltype used.  There should be a single XML tag
728    # containing all the cells defined.  It's an error if there's more than
729    # one of these, or if there isn't one. 
730    foreach type {
731        triangles quads tetraheadrons voxels hexaheadrons wedges pyramids }  {
732        if { [$xmlobj get $path.cells.$type] != "" } {
733            set celltype $type
734            incr count
735        }
736    }
737    if { $count == 0 } {
738        puts stderr "no cell types found"
739        return 0
740    }
741    if { $count > 1 } {
742        puts stderr "too many cell types tags found"
743        return 0
744    }
745    set cells [$xmlobj get $path.cells.$celltype]
746    if { $cells == "" } {
747        puts stderr "no data for celltype \"$celltype\""
748        return 0
749    }
750    # Get the number of vertices for each cell and the corresponding VTK id.
751    set numVertices $celltype2vertices($celltype)
752    set vtkid       $celltype2vtkid($celltype)
753   
754    if { $_dim == 2 } {
755        # Load the points into one big vector, then split out the x and y
756        # coordinates.
757        set all [blt::vector create \#auto]
758        $all set $points
759        set xv [blt::vector create \#auto]
760        set yv [blt::vector create \#auto]
761        $all split $xv $yv
762        set _numPoints [$xv length]
763
764        set data {}
765        set _numCells 0
766        set celltypes {}
767        # Build cell information.  Look line by line.  Most cell types
768        # have a specific number of vertices, but there are a few that
769        # have a variable number of vertices.
770        foreach line [split $cells \n] {
771            set numVerticies [llength $line]
772            if { $numVerticies == 0 } {
773                continue;               # Skip blank lines
774            }
775            append data "    $numVertices $line\n"
776            append celltypes "$vtkid\n"
777            incr _numCells
778        }
779        append out "DATASET UNSTRUCTURED_GRID\n"
780        append out "POINTS $_numPoints float\n"
781        foreach x [$xv range 0 end] y [$yv range 0 end] {
782            append out "    $x $y 0\n"
783        }
784        append out "CELLS $_numCells [expr $_numCells * ($numVertices + 1)]\n"
785        append out $data
786        append out "CELL_TYPES $_numCells\n"
787        append out $celltypes
788        set _limits(x) [$xv limits]
789        set _limits(y) [$yv limits]
790        blt::vector destroy $xv $yv $all
791    } else {
792        # Load the points into one big vector, then split out the x, y, and z
793        # coordinates.
794        set all [blt::vector create \#auto]
795        $all set $points
796        set xv [blt::vector create \#auto]
797        set yv [blt::vector create \#auto]
798        set zv [blt::vector create \#auto]
799        $all split $xv $yv $zv
800        set _numPoints [$xv length]
801
802        set data {}
803        set _numCells 0
804        set celltypes {}
805        # Build cell information.  Look line by line.  Most cell types
806        # have a specific number of vertices, but there are a few that
807        # have a variable number of vertices.
808        foreach line [split $cells \n] {
809            set numVerticies [llength $line]
810            if { $numVerticies == 0 } {
811                continue;               # Skip blank lines
812            }
813            append data "    $numVertices $line\n"
814            append celltypes "$vtkid\n"
815            incr _numCells
816        }
817        append out "DATASET UNSTRUCTURED_GRID\n"
818        append out "POINTS $_numPoints float\n"
819        foreach x [$xv range 0 end] y [$yv range 0 end] z [$zv range 0 end] {
820            append out "    $x $y $z\n"
821        }
822        append out "CELLS $_numCells [expr $_numCells * ($numVerticies + 1)]\n"
823        append out $data
824        append out "CELL_TYPES $_numCells\n"
825        append out $celltypes
826        set _limits(x) [$xv limits]
827        set _limits(y) [$yv limits]
828        set _limits(z) [$zv limits]
829        blt::vector destroy $xv $yv $zv $all
830    }
831    set _vtkdata $out
832    set _isValid 1
833}
834
835
836itcl::body Rappture::Mesh::ReadTriangles { xmlobj path } {
837    set _type "triangles"
838    set _dim 2
839    set triangles [$xmlobj get $path.triangles.indices]
840    if { $triangles == "" } {
841        puts stderr "no triangle indices specified in mesh"
842        return 0
843    }
844    set xcoords [$xmlobj get $path.triangles.xcoords]
845    set ycoords [$xmlobj get $path.triangles.ycoords]
846    set points [$xmlobj get $path.triangles.points]
847    set xv [blt::vector create \#auto]
848    set yv [blt::vector create \#auto]
849    if { $xcoords != "" && $ycoords != "" } {
850        $xv set $xcoords
851        $yx set $ycoords
852    } elseif { $points != "" } {
853        set all [blt::vector create \#auto]
854        $all set $points
855        $all split $xv $yv
856    } else {
857        puts stderr "missing either <xcoords>, <ycoords>, or <points> for triangular mesh"
858        return 0
859    }
860    set _numPoints [$xv length]
861    set count 0
862    set data {}
863    foreach { a b c } $triangles {
864        append data "    3 $a $b $c\n"
865        incr count
866    }
867    set celltypes {}
868    for { set i 0 } { $i < $count } { incr i } {
869        append celltypes "5\n"
870    }
871    append out "DATASET UNSTRUCTURED_GRID\n"
872    append out "POINTS $_numPoints float\n"
873    foreach x [$xv range 0 end] y [$yv range 0 end] {
874        append out "    $x $y 0\n"
875    }
876    append out "CELLS $count [expr $count * 4]\n"
877    append out $data
878    append out "CELL_TYPES $count\n"
879    append out $celltypes
880    set _limits(x) [$xv limits]
881    set _limits(y) [$yv limits]
882    set _vtkdata $out
883    set _isValid 1
884}
885
886itcl::body Rappture::Mesh::ReadUnstructuredGrid { xmlobj path } {
887    set _type "unstructured_grid"
888
889    set data [$xmlobj get $path.unstructured_grid.points]
890    Rappture::ReadPoints $data _dim points
891    if { $points == "" } {
892        puts stderr "no <points> found for <cells> mesh"
893        return 0
894    }
895    set cells [$xmlobj get $path.unstructured_grid.cells]
896    if { $cells == "" } {
897        puts stderr "no <cells> description found unstructured_grid"
898        return 0
899    }
900    set celltypes [$xmlobj get $path.unstructured_grid.celltypes]
901    if { $cells == "" } {
902        puts stderr "no <cells> description found unstructured_grid."
903        return 0
904    }
905    set lines [split $cells \n]
906   
907    set all [blt::vector create \#auto]
908    set xv [blt::vector create \#auto]
909    set yv [blt::vector create \#auto]
910    set zv [blt::vector create \#auto]
911    $all set $points
912    if { $_dim == 2 } {
913        # Load the points into one big vector, then split out the x and y
914        # coordinates.
915        $all split $xv $yv
916        set _numPoints [$xv length]
917
918        set data {}
919        set count 0
920        set _numCells 0
921        foreach line $lines {
922            set length [llength $line]
923            if { $length == 0 } {
924                continue
925            }
926            set celltype [lindex $celltypes $_numCells]
927            if { $celltype == "" } {
928                puts stderr "can't find celltype for line \"$line\""
929                return 0
930            }
931            append data " $line\n"
932            incr count $length
933            incr _numCells
934        }
935        append out "DATASET UNSTRUCTURED_GRID\n"
936        append out "POINTS $_numPoints float\n"
937        $zv seq 0 0 [$xv length];       # Make an all zeroes vector.
938        $all merge $xv $yv $zv
939        append out [$all range 0 end]
940        append out "CELLS $_numCells $count\n"
941        append out $data
942        append out "CELL_TYPES $_numCells\n"
943        append out $celltypes
944        set _limits(x) [$xv limits]
945        set _limits(y) [$yv limits]
946    } else {
947        # Load the points into one big vector, then split out the x, y, and z
948        # coordinates.
949        $all split $xv $yv $zv
950        set _numPoints [$xv length]
951
952        set data {}
953        set count 0
954        set _numCells 0
955        foreach line $lines {
956            set length [llength $line]
957            if { $length == 0 } {
958                continue
959            }
960            set celltype [lindex $celltypes $_numCells]
961            if { $celltype == "" } {
962                puts stderr "can't find celltype for line \"$line\""
963                return 0
964            }
965            append data " $line\n"
966            incr count $length
967            incr _numCells
968        }
969        append out "DATASET UNSTRUCTURED_GRID\n"
970        append out "POINTS $_numPoints float\n"
971        $all merge $xv $yv $zv
972        append out [$all range 0 end]
973        append out "\n"
974        append out "CELLS $_numCells $count\n"
975        append out $data
976        append out "CELL_TYPES $_numCells\n"
977        append out $celltypes
978        set _limits(x) [$xv limits]
979        set _limits(y) [$yv limits]
980        set _limits(z) [$zv limits]
981    }
982    blt::vector destroy $xv $yv $zv $all
983    set _vtkdata $out
984    set _isValid 1
985}
986
987# ----------------------------------------------------------------------
988# USAGE: ReadNodesElements <xmlobj> <path>
989#
990# Used internally to build a mesh representation based on nodes and
991# elements stored in the XML.
992# ----------------------------------------------------------------------
993itcl::body Rappture::Mesh::ReadNodesElements {xmlobj path} {
994    set type "nodeselements"
995    set count 0
996
997    # Scan for nodes.  Each node represents a vertex.
998    set data {}
999    foreach cname [$xmlobj children -type node $path] {
1000        append data "[$xmlobj get $path.$cname]\n"
1001    }   
1002    Rappture::ReadPoints $data _dim points
1003    if { $_dim == 2 } {
1004        set all [blt::vector create \#auto]
1005        set xv [blt::vector create \#auto]
1006        set yv [blt::vector create \#auto]
1007        set zv [blt::vector create \#auto]
1008        $all set $points
1009        $all split $xv $yv
1010        set _numPoints [$xv length]
1011        foreach axis {x y} {
1012            set vector [set ${axis}v]
1013            set _limits($axis) [${vector} limits]
1014        }
1015        # 2D Dataset. All Z coordinates are 0
1016        $zv seq 0.0 0.0 $_numPoints
1017        $all merge $xv $yv $zv
1018        set points [$all range 0 end]
1019        blt::vector destroy $all $xv $yv $zv
1020    } elseif { $_dim == 3 } {
1021        set all [blt::vector create \#auto]
1022        set xv [blt::vector create \#auto]
1023        set yv [blt::vector create \#auto]
1024        set zv [blt::vector create \#auto]
1025        $all set $points
1026        $all split $xv $yv $zv
1027        set _numPoints [$xv length]
1028        foreach axis {x y z} {
1029            set vector [set ${axis}v]
1030            set _limits($axis) [${vector} limits]
1031        }
1032        set points [$all range 0 end]
1033        blt::vector destroy $all $xv $yv $zv
1034    } else {
1035        error "bad dimension \"$_dim\" for nodes mesh"
1036    }
1037    array set node2celltype {
1038        3 5
1039        4 10
1040        8 11
1041        6 13
1042        5 14
1043    }
1044    set count 0
1045    set _numCells 0
1046    set celltypes {}
1047    set data {}
1048    # Next scan for elements.  Each element represents a cell.
1049    foreach cname [$xmlobj children -type element $path] {
1050        set nodeList [$_mesh get $cname.nodes]
1051        set numNodes [llength $nodeList]
1052        if { ![info exists node2celltype($numNodes)] } {
1053            puts stderr "unknown number of indices \"$_numNodes\": should be 3, 4, 5, 6, or 8"
1054            return 0
1055        }
1056        set celltype $node2celltype($numNodes)
1057        append celltypes "  $celltype\n"
1058        append data "  $numNodes $nodeList\n"
1059        incr _numCells
1060        incr count $numNodes
1061        incr count;                     # One extra for the VTK celltype id.
1062    }
1063
1064    append out "DATASET UNSTRUCTURED_GRID\n"
1065    append out "POINTS $_numPoints float\n"
1066    append out $points
1067    append out "CELLS $_numCells $count\n"
1068    append out $data
1069    append out "CELL_TYPES $_numCells\n"
1070    append out $celltypes
1071    append out "\n"
1072    set _vtkdata $out
1073    set _isValid 1
1074}
Note: See TracBrowser for help on using the repository browser.