hedge.visualization

1 """Visualization for global DG functions. Supports VTK, Silo, etc.""" 2 3 __copyright__ = "Copyright (C) 2007 Andreas Kloeckner" 4 5 __license__ = """ 6 This program is free software: you can redistribute it and/or modify 7 it under the terms of the GNU General Public License as published by 8 the Free Software Foundation, either version 3 of the License, or 9 (at your option) any later version. 10 11 This program is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with this program. If not, see U{http://www.gnu.org/licenses/}. 18 """ 19 20 21 22 23 import hedge.tools 24 import numpy 25 26 27 28

29 -class Visualizer(object):

30 pass

31 32 33 34 35 # legacy vtk ------------------------------------------------------------------

36 -def _three_vector(x):

37 if len(x) == 3: 38 return x 39 elif len(x) == 2: 40 return x[0], x[1], 0. 41 elif len(x) == 1: 42 return x[0], 0, 0.

43 44 45 46

47 -class LegacyVtkFile(object):

48 - def __init__(self, pathname, structure, description="Hedge visualization"):

49 self.pathname = pathname 50 self.structure = structure 51 self.description = description 52 53 self.pointdata = [] 54 self.is_closed = False

55

56 - def __fin__(self):

57 self.close()

58

59 - def __enter__(self):

60 return self

61

62 - def __exit__(self, exc_type, exc_value, traceback):

63 self.close()

64

65 - def close(self):

66 if not self.is_closed: 67 from pyvtk import PointData, VtkData 68 vtk = VtkData(self.structure, 69 self.description, 70 PointData(*self.pointdata)) 71 vtk.tofile(self.pathname) 72 self.is_closed = True

73 74 75 76

77 -class LegacyVtkVisualizer(Visualizer):

78 - def __init__(self, discr):

79 from pyvtk import PolyData 80 81 points = [_three_vector(p) for p in discr.nodes] 82 polygons = [] 83 84 for eg in discr.element_groups: 85 ldis = eg.local_discretization 86 for el, (el_start, el_stop) in zip(eg.members, eg.ranges): 87 polygons += [[el_start+j for j in element] 88 for element in ldis.get_submesh_indices()] 89 90 self.structure = PolyData(points=points, polygons=polygons)

91

92 - def make_file(self, pathname, pcontext=None):

93 if pcontext is not None: 94 if len(pcontext.ranks) > 1: 95 raise RuntimeError, "Legacy VTK does not suport parallel visualization" 96 return LegacyVtkFile(pathname+".vtk", self.structure)

97

98 - def add_data(self, vtkfile, scalars=[], vectors=[], scale_factor=1):

99 from pyvtk import Scalars, Vectors 100 101 vtkfile.pointdata.extend( 102 Scalars(numpy.array(scale_factor*field), 103 name=name, lookup_table="default") 104 for name, field in scalars) 105 vtkfile.pointdata.extend( 106 Vectors([_three_vector(scale_factor*v) 107 for v in zip(field)], name=name) 108 for name, field in vectors)

109 110 111 112 113 # xml vtk ---------------------------------------------------------------------

114 -class VtkFile(hedge.tools.Closable):

115 - def __init__(self, pathname, grid, filenames=None, compressor=None):

116 """`compressor` may be what ever is accepted by L{hedge.vtk.make_vtkfile}.""" 117 hedge.tools.Closable.__init__(self) 118 self.pathname = pathname 119 self.grid = grid 120 self.compressor = compressor

121

122 - def get_head_pathname(self):

123 return self.pathname

124

125 - def do_close(self):

126 from pytools import assert_not_a_file 127 assert_not_a_file(self.pathname) 128 129 130 outf = file(self.pathname, "w") 131 132 from hedge.vtk import AppendedDataXMLGenerator 133 AppendedDataXMLGenerator(self.compressor)(self.grid).write(outf) 134 135 #from hedge.vtk import InlineXMLGenerator 136 #InlineXMLGenerator(self.compressor)(self.grid).write(outf) 137 138 outf.close()

139 140 141 142 143

144 -class ParallelVtkFile(VtkFile):

145 - def __init__(self, pathname, grid, index_pathname, pathnames=None, compressor=None):

146 VtkFile.__init__(self, pathname, grid) 147 self.index_pathname = index_pathname 148 self.pathnames = pathnames 149 self.compressor = compressor

150

151 - def get_head_pathname(self):

152 return self.index_pathname

153

154 - def do_close(self):

155 VtkFile.do_close(self) 156 157 from hedge.vtk import ParallelXMLGenerator 158 159 outf = file(self.index_pathname, "w") 160 ParallelXMLGenerator(self.pathnames)(self.grid).write(outf) 161 outf.close()

162 163 164 165 166 167

168 -class VtkVisualizer(Visualizer, hedge.tools.Closable):

169 - def __init__(self, discr, pcontext=None, basename=None, compressor=None):

170 hedge.tools.Closable.__init__(self) 171 172 from pytools import assert_not_a_file 173 174 if basename is not None: 175 self.pvd_name = basename+".pvd" 176 assert_not_a_file(self.pvd_name) 177 else: 178 self.pvd_name = None 179 180 self.pcontext = pcontext 181 self.compressor = compressor 182 183 if self.pcontext is None or self.pcontext.is_head_rank: 184 self.timestep_to_pathnames = {} 185 else: 186 self.timestep_to_pathnames = None 187 188 from hedge.vtk import UnstructuredGrid, DataArray, \ 189 VTK_LINE, VTK_TRIANGLE, VTK_TETRA, VF_LIST_OF_VECTORS 190 from hedge.mesh import Interval, Triangle, Tetrahedron 191 192 # For now, we use IntVector here because the Python allocator 193 # is somewhat reluctant to return allocated chunks of memory 194 # to the OS. 195 from hedge._internal import IntVector 196 cells = IntVector() 197 cell_types = IntVector() 198 199 for eg in discr.element_groups: 200 ldis = eg.local_discretization 201 smi = ldis.get_submesh_indices() 202 203 cells.reserve(len(cells)+len(smi)*len(eg.members)) 204 for el, el_slice in zip(eg.members, eg.ranges): 205 for element in smi: 206 for j in element: 207 cells.append(el_slice.start+j) 208 209 if ldis.geometry is Interval: 210 vtk_eltype = VTK_LINE 211 elif ldis.geometry is Triangle: 212 vtk_eltype = VTK_TRIANGLE 213 elif ldis.geometry is Tetrahedron: 214 vtk_eltype = VTK_TETRA 215 else: 216 raise RuntimeError, "unsupported element type: %s" % ldis.geometry 217 218 cell_types.extend([vtk_eltype] * len(smi) * len(eg.members)) 219 220 self.grid = UnstructuredGrid( 221 (len(discr), 222 DataArray("points", discr.nodes, vector_format=VF_LIST_OF_VECTORS)), 223 numpy.asarray(cells), 224 cell_types=numpy.asarray(cell_types, dtype=numpy.uint8))

225 226

227 - def update_pvd(self):

228 if self.pvd_name and self.timestep_to_pathnames: 229 from hedge.vtk import XMLRoot, XMLElement, make_vtkfile 230 231 collection = XMLElement("Collection") 232 233 vtkf = make_vtkfile(collection.tag, compressor=None) 234 xmlroot = XMLRoot(vtkf) 235 236 vtkf.add_child(collection) 237 238 tsteps = self.timestep_to_pathnames.keys() 239 tsteps.sort() 240 for i, time in enumerate(tsteps): 241 for part, pathname in enumerate(self.timestep_to_pathnames[time]): 242 collection.add_child(XMLElement( 243 "DataSet", 244 timestep=time, part=part, file=pathname)) 245 outf = open(self.pvd_name, "w") 246 xmlroot.write(outf) 247 outf.close()

248

249 - def do_close(self):

250 self.update_pvd()

251

252 - def make_file(self, pathname):

253 """FIXME 254 255 An appropriate extension (including the dot) is automatically 256 appended to `pathname'. 257 """ 258 if self.pcontext is None or len(self.pcontext.ranks) == 1: 259 return VtkFile(pathname+"."+self.grid.vtk_extension(), 260 self.grid.copy(), 261 compressor=self.compressor 262 ) 263 else: 264 filename_pattern = ( 265 pathname + "-%05d." + self.grid.vtk_extension()) 266 if self.pcontext.is_head_rank: 267 return ParallelVtkFile( 268 filename_pattern % self.pcontext.rank, 269 self.grid.copy(), 270 index_pathname="%s.p%s" % ( 271 pathname, self.grid.vtk_extension()), 272 pathnames=[ 273 filename_pattern % rank for rank in self.pcontext.ranks], 274 compressor=self.compressor 275 ) 276 else: 277 return VtkFile( 278 filename_pattern % self.pcontext.rank, 279 self.grid.copy(), 280 compressor=self.compressor 281 )

282

283 - def register_pathname(self, time, pathname):

284 if time is not None and self.timestep_to_pathnames is not None: 285 self.timestep_to_pathnames.setdefault(time, []).append(pathname) 286 287 # When we are run under MPI and cancelled by Ctrl+C, destructors 288 # do not get called. Therefore, we just spend the (hopefully negligible) 289 # time to update the PVD index every few data additions. 290 if len(self.timestep_to_pathnames) % 5 == 0: 291 self.update_pvd()

292

293 - def add_data(self, visf, variables=[], scalars=[], vectors=[], time=None, step=None, 294 scale_factor=1):

295 if scalars or vectors: 296 import warnings 297 warnings.warn("`scalars' and `vectors' arguments are deprecated", 298 DeprecationWarning) 299 variables = scalars + vectors 300 301 from hedge.vtk import DataArray, VF_LIST_OF_COMPONENTS 302 for name, field in variables: 303 visf.grid.add_pointdata(DataArray(name, scale_factor*field, 304 vector_format=VF_LIST_OF_COMPONENTS)) 305 306 self.register_pathname(time, visf.get_head_pathname())

307 308 309 310 311 312 # silo ------------------------------------------------------------------------

313 -class SiloMeshData(object):

314 - def __init__(self, dim, coords, element_groups):

315 self.coords = coords 316 317 from pylo import IntVector 318 self.ndims = dim 319 self.nodelist = IntVector() 320 self.shapetypes = IntVector() 321 self.shapesizes = IntVector() 322 self.shapecounts = IntVector() 323 self.nzones = 0 324 325 for nodelist_size_estimate, eg, ldis in element_groups: 326 poly_count = 0 327 poly_length = None 328 self.nodelist.reserve(len(self.nodelist) + nodelist_size_estimate) 329 for polygon in eg: 330 prev_nodelist_len = len(self.nodelist) 331 for i in polygon: 332 self.nodelist.append(i) 333 poly_count += 1 334 poly_length = len(self.nodelist) - prev_nodelist_len 335 336 if poly_count: 337 try: 338 from pylo import DB_ZONETYPE_TRIANGLE, DB_ZONETYPE_TET 339 except ImportError: 340 pass 341 else: 342 from hedge.mesh import Triangle, Tetrahedron 343 if ldis.geometry is Triangle: 344 self.shapetypes.append(DB_ZONETYPE_TRIANGLE) 345 elif ldis.geometry is Tetrahedron: 346 self.shapetypes.append(DB_ZONETYPE_TET) 347 else: 348 raise RuntimeError, "unsupported element type: %s" % ldis.geometry 349 350 self.shapesizes.append(poly_length) 351 self.shapecounts.append(poly_count) 352 self.nzones += poly_count

353

354 - def put_mesh(self, silo, zonelist_name, mesh_name, mesh_opts):

355 if self.shapetypes: 356 assert len(self.shapetypes) == len(self.shapesizes) 357 silo.put_zonelist_2(zonelist_name, self.nzones, self.ndims, self.nodelist, 358 0, 0, self.shapetypes, self.shapesizes, self.shapecounts) 359 else: 360 silo.put_zonelist(zonelist_name, self.nzones, self.ndims, self.nodelist, 361 self.shapesizes, self.shapecounts) 362 363 silo.put_ucdmesh(mesh_name, [], self.coords, self.nzones, 364 zonelist_name, None, mesh_opts)

365 366 367 368

369 -class SiloVisualizer(Visualizer):

370 - def __init__(self, discr, pcontext=None):

371 self.discr = discr 372 self.pcontext = pcontext 373 374 self.generated = False

375

376 - def _generate(self):

377 # only generate vis data when vis is really needed. 378 # saves startup time when debugging. 379 def generate_fine_elements(eg): 380 ldis = eg.local_discretization 381 smi = ldis.get_submesh_indices() 382 for el, el_slice in zip(eg.members, eg.ranges): 383 for element in smi: 384 yield [el_slice.start+j for j in element]

385 386 def generate_fine_element_groups(): 387 for eg in discr.element_groups: 388 ldis = eg.local_discretization 389 smi = ldis.get_submesh_indices() 390 nodelist_size_estimate = len(eg.members) * len(smi) * len(smi[0]) 391 yield nodelist_size_estimate, generate_fine_elements(eg), ldis

392 393 def generate_coarse_elements(eg): 394 for el in eg.members: 395 yield el.vertex_indices 396 397 def generate_coarse_element_groups(): 398 for eg in discr.element_groups: 399 if eg.members: 400 nodelist_size_estimate = len(eg.members) \ 401 * len(eg.members[0].vertex_indices) 402 else: 403 nodelist_size_estimate = 0 404 405 yield (nodelist_size_estimate, generate_coarse_elements(eg), 406 eg.local_discretization) 407 408 discr = self.discr 409 self.dim = discr.dimensions 410 if self.dim != 1: 411 self.fine_mesh = SiloMeshData(self.dim, 412 numpy.asarray(discr.nodes.T, order="C"), 413 generate_fine_element_groups()) 414 self.coarse_mesh = SiloMeshData(self.dim, 415 numpy.asarray(discr.mesh.points.T, order="C"), 416 generate_coarse_element_groups()) 417 else: 418 self.xvals = numpy.asarray(discr.nodes.T, order="C") 419 420 self.generated = True 421

422 - def close(self):

423 pass

424

425 - def make_file(self, pathname):

426 """This function returns either a pylo.SiloFile or a 427 pylo.ParallelSiloFile, depending on the ParallelContext 428 under which we are running 429 430 An extension of .silo is automatically appended to `pathname'. 431 """ 432 if not self.generated: 433 self._generate() 434 435 if self.pcontext is None or len(self.pcontext.ranks) == 1: 436 from pylo import SiloFile 437 return SiloFile(pathname+".silo") 438 else: 439 from pylo import ParallelSiloFile 440 return ParallelSiloFile( 441 pathname, 442 self.pcontext.rank, self.pcontext.ranks)

443

444 - def add_data(self, silo, variables=[], scalars=[], vectors=[], expressions=[], 445 time=None, step=None, scale_factor=1):

446 if scalars or vectors: 447 import warnings 448 warnings.warn("`scalars' and `vectors' arguments are deprecated", 449 DeprecationWarning) 450 variables = scalars + vectors 451 452 from pylo import DB_NODECENT, DBOPT_DTIME, DBOPT_CYCLE 453 454 # put mesh coordinates 455 mesh_opts = {} 456 if time is not None: 457 mesh_opts[DBOPT_DTIME] = float(time) 458 if step is not None: 459 mesh_opts[DBOPT_CYCLE] = int(step) 460 461 if self.dim == 1: 462 for name, field in variables: 463 from hedge.tools import is_obj_array 464 if is_obj_array(field): 465 AXES = ["x", "y", "z", "w"] 466 for i, f_i in enumerate(field): 467 silo.put_curve(name+AXES[i], self.xvals, 468 scale_factor*f_i, mesh_opts) 469 else: 470 silo.put_curve(name, self.xvals, 471 scale_factor*field, mesh_opts) 472 else: 473 self.fine_mesh.put_mesh(silo, "finezonelist", "finemesh", mesh_opts) 474 self.coarse_mesh.put_mesh(silo, "coarsezonelist", "mesh", mesh_opts) 475 476 from hedge.tools import log_shape 477 478 # put data 479 for name, field in variables: 480 ls = log_shape(field) 481 if ls != () and ls[0] > 1: 482 assert len(ls) == 1 483 silo.put_ucdvar(name, "finemesh", 484 ["%s_comp%d" % (name, i) 485 for i in range(ls[0])], 486 scale_factor*field, DB_NODECENT) 487 else: 488 if ls != (): 489 field = field[0] 490 silo.put_ucdvar1(name, "finemesh", scale_factor*field, DB_NODECENT) 491 492 if expressions: 493 silo.put_defvars("defvars", expressions)

494 495 496 497 498 # tools -----------------------------------------------------------------------

499 -def get_rank_partition(pcon, discr):

500 vec = discr.volume_zeros() 501 vec[:] = pcon.rank 502 return vec

503

Source Code for Module hedge.visualization