Context Navigation

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

Last change on this file since 3534 was 3530, checked in by gah, 12 years ago
leif's fix for voxels to hexahedrons
File size: 36.5 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: