vtk_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // C++ includes
20 #include <fstream>
21 
22 // Local includes
23 #include "libmesh/libmesh_config.h"
24 #include "libmesh/vtk_io.h"
25 #include "libmesh/mesh_base.h"
27 #include "libmesh/edge_edge2.h"
28 #include "libmesh/edge_edge3.h"
29 #include "libmesh/face_tri3.h"
30 #include "libmesh/face_tri6.h"
31 #include "libmesh/face_quad4.h"
32 #include "libmesh/face_quad8.h"
33 #include "libmesh/face_quad9.h"
34 #include "libmesh/cell_tet4.h"
35 #include "libmesh/cell_tet10.h"
36 #include "libmesh/cell_prism6.h"
37 #include "libmesh/cell_prism15.h"
38 #include "libmesh/cell_prism18.h"
39 #include "libmesh/cell_pyramid5.h"
40 #include "libmesh/cell_hex8.h"
41 #include "libmesh/cell_hex20.h"
42 #include "libmesh/cell_hex27.h"
43 #include "libmesh/numeric_vector.h"
44 #include "libmesh/system.h"
45 #include "libmesh/mesh_data.h"
46 
47 #ifdef LIBMESH_HAVE_VTK
48 
49 // Tell VTK not to use old header files
50 #define VTK_LEGACY_REMOVE
51 
52 #include "vtkXMLUnstructuredGridReader.h"
53 #include "vtkXMLUnstructuredGridWriter.h"
54 #include "vtkXMLPUnstructuredGridWriter.h"
55 #include "vtkUnstructuredGrid.h"
56 #include "vtkGenericGeometryFilter.h"
57 #include "vtkCellArray.h"
58 #include "vtkCellData.h"
59 #include "vtkConfigure.h"
60 #include "vtkDoubleArray.h"
61 #include "vtkPointData.h"
62 #include "vtkPoints.h"
63 #include "vtkSmartPointer.h"
64 
65 // A convenient macro for comparing VTK versions. Returns 1 if the
66 // current VTK version is < major.minor.subminor and zero otherwise.
67 //
68 // It relies on the VTK version numbers detected during configure. Note that if
69 // LIBMESH_HAVE_VTK is not defined, none of the LIBMESH_DETECTED_VTK_VERSION_* variables will
70 // be defined either.
71 #define VTK_VERSION_LESS_THAN(major,minor,subminor) \
72  ((LIBMESH_DETECTED_VTK_VERSION_MAJOR < (major) || \
73  (LIBMESH_DETECTED_VTK_VERSION_MAJOR == (major) && (LIBMESH_DETECTED_VTK_VERSION_MINOR < (minor) || \
74  (LIBMESH_DETECTED_VTK_VERSION_MINOR == (minor) && \
75  LIBMESH_DETECTED_VTK_VERSION_SUBMINOR < (subminor))))) ? 1 : 0)
76 
77 
78 
79 
80 #endif //LIBMESH_HAVE_VTK
81 
82 namespace libMesh
83 {
84 
85 #ifdef LIBMESH_HAVE_VTK
87 {
88  vtkIdType celltype = VTK_EMPTY_CELL; // initialize to something to avoid compiler warning
89 
90  switch(type)
91  {
92  case EDGE2:
93  celltype = VTK_LINE;
94  break;
95  case EDGE3:
96  celltype = VTK_QUADRATIC_EDGE;
97  break;// 1
98  case TRI3:
99  celltype = VTK_TRIANGLE;
100  break;// 3
101  case TRI6:
102  celltype = VTK_QUADRATIC_TRIANGLE;
103  break;// 4
104  case QUAD4:
105  celltype = VTK_QUAD;
106  break;// 5
107  case QUAD8:
108  celltype = VTK_QUADRATIC_QUAD;
109  break;// 6
110  case TET4:
111  celltype = VTK_TETRA;
112  break;// 8
113  case TET10:
114  celltype = VTK_QUADRATIC_TETRA;
115  break;// 9
116  case HEX8:
117  celltype = VTK_HEXAHEDRON;
118  break;// 10
119  case HEX20:
120  celltype = VTK_QUADRATIC_HEXAHEDRON;
121  break;// 12
122  case HEX27:
123  celltype = VTK_TRIQUADRATIC_HEXAHEDRON;
124  break;
125  case PRISM6:
126  celltype = VTK_WEDGE;
127  break;// 13
128  case PRISM15:
129  celltype = VTK_QUADRATIC_WEDGE;
130  break;// 14
131  case PRISM18:
132  celltype = VTK_BIQUADRATIC_QUADRATIC_WEDGE;
133  break;// 15
134  case PYRAMID5:
135  celltype = VTK_PYRAMID;
136  break;// 16
137 #if VTK_MAJOR_VERSION > 5 || (VTK_MAJOR_VERSION == 5 && VTK_MINOR_VERSION > 0)
138  case QUAD9:
139  celltype = VTK_BIQUADRATIC_QUAD;
140  break;
141 #else
142  case QUAD9:
143 #endif
144  case EDGE4:
145  case INFEDGE2:
146  case INFQUAD4:
147  case INFQUAD6:
148  case INFHEX8:
149  case INFHEX16:
150  case INFHEX18:
151  case INFPRISM6:
152  case INFPRISM12:
153  case NODEELEM:
154  case INVALID_ELEM:
155  default:
156  {
157  libMesh::err<<"element type "<<type<<" not implemented"<<std::endl;
158  libmesh_error();
159  break;
160  }
161  }
162  return celltype;
163 }
164 
165 
166 
168 {
170 
171  // containers for points and coordinates of points
172  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
173  vtkSmartPointer<vtkDoubleArray> pcoords = vtkSmartPointer<vtkDoubleArray>::New();
174  pcoords->SetNumberOfComponents(LIBMESH_DIM);
175  points->SetNumberOfPoints(mesh.n_local_nodes()); // it seems that it needs this to prevent a segfault
176 
177  unsigned int local_node_counter = 0;
178 
181  for (; nd != nd_end; nd++, local_node_counter++)
182  {
183  Node* node = (*nd);
184 
185  double pnt[LIBMESH_DIM];
186  for (unsigned int i=0; i<LIBMESH_DIM; ++i)
187  pnt[i] = (*node)(i);
188 
189  // Fill mapping between global and local node numbers
190  _local_node_map[node->id()] = local_node_counter;
191 
192  // add point
193  pcoords->InsertNextTupleValue(pnt);
194  }
195 
196  // add coordinates to points
197  points->SetData(pcoords);
198 
199  // add points to grid
200  _vtk_grid->SetPoints(points);
201 }
202 
203 
204 
206 {
208 
209  vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
210  vtkSmartPointer<vtkIdList> pts = vtkSmartPointer<vtkIdList>::New();
211 
212  std::vector<int> types(mesh.n_active_local_elem());
213  unsigned active_element_counter = 0;
214 
215  vtkSmartPointer<vtkIntArray> elem_id = vtkSmartPointer<vtkIntArray>::New();
216  elem_id->SetName("libmesh_elem_id");
217  elem_id->SetNumberOfComponents(1);
218 
219  vtkSmartPointer<vtkIntArray> subdomain_id = vtkSmartPointer<vtkIntArray>::New();
220  subdomain_id->SetName("subdomain_id");
221  subdomain_id->SetNumberOfComponents(1);
222 
225  for (; it != end; ++it, ++active_element_counter)
226  {
227  Elem *elem = *it;
228 
229  pts->SetNumberOfIds(elem->n_nodes());
230 
231  // get the connectivity for this element
232  std::vector<dof_id_type> conn;
233  elem->connectivity(0, VTK, conn);
234 
235  for (unsigned int i=0; i<conn.size(); ++i)
236  {
237  // If the node ID is not found in the _local_node_map, we'll
238  // add it to the _vtk_grid. NOTE[JWP]: none of the examples
239  // I have actually enters this section of code...
240  if (_local_node_map.find(conn[i]) == _local_node_map.end())
241  {
242  dof_id_type global_node_id = elem->node(i);
243 
244  const Node* the_node = mesh.node_ptr(global_node_id);
245 
246  // Error checking...
247  if (the_node == NULL)
248  {
249  libMesh::err << "Error getting pointer to node "
250  << global_node_id
251  << "!" << std::endl;
252  libmesh_error();
253  }
254 
255  // InsertNextPoint accepts either a double or float array of length 3.
256  Real pt[3] = {0., 0., 0.};
257  for (unsigned int d=0; d<LIBMESH_DIM; ++d)
258  pt[d] = (*the_node)(d);
259 
260  // Insert the point into the _vtk_grid
261  vtkIdType local = _vtk_grid->GetPoints()->InsertNextPoint(pt);
262 
263  // Update the _local_node_map with the ID returned by VTK
264  _local_node_map[global_node_id] = local;
265  }
266 
267  // Otherwise, the node ID was found in the _local_node_map, so
268  // insert it into the vtkIdList.
269  pts->InsertId(i, _local_node_map[conn[i]]);
270  }
271 
272  vtkIdType vtkcellid = cells->InsertNextCell(pts);
273  types[active_element_counter] = this->get_elem_type(elem->type());
274  elem_id->InsertTuple1(vtkcellid, elem->id());
275  subdomain_id->InsertTuple1(vtkcellid, elem->subdomain_id());
276  } // end loop over active elements
277 
278  _vtk_grid->SetCells(&types[0], cells);
279  _vtk_grid->GetCellData()->AddArray(elem_id);
280  _vtk_grid->GetCellData()->AddArray(subdomain_id);
281 }
282 
283 
284 
285 /*
286  * FIXME: This is known to write nonsense on AMR meshes
287  * and it strips the imaginary parts of complex Numbers
288  */
289 void VTKIO::system_vectors_to_vtk(const EquationSystems& es, vtkUnstructuredGrid*& grid)
290 {
292  {
293  std::map<std::string, std::vector<Number> > vecs;
294  for (unsigned int i=0; i<es.n_systems(); ++i)
295  {
296  const System& sys = es.get_system(i);
299  for (; it!= v_end; ++it)
300  {
301  // for all vectors on this system
302  std::vector<Number> values;
303  // libMesh::out<<"it "<<it->first<<std::endl;
304 
305  it->second->localize_to_one(values, 0);
306  // libMesh::out<<"finish localize"<<std::endl;
307  vecs[it->first] = values;
308  }
309  }
310 
311  std::map<std::string, std::vector<Number> >::iterator it = vecs.begin();
312 
313  for (; it!=vecs.end(); ++it)
314  {
315  vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
316  data->SetName(it->first.c_str());
317  libmesh_assert_equal_to (it->second.size(), es.get_mesh().n_nodes());
318  data->SetNumberOfValues(it->second.size());
319 
320  for (unsigned int i=0; i<it->second.size(); ++i)
321  {
322 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
323  libmesh_do_once (libMesh::err << "Only writing the real part for complex numbers!\n"
324  << "if you need this support contact " << LIBMESH_PACKAGE_BUGREPORT
325  << std::endl);
326  data->SetValue(i, it->second[i].real());
327 #else
328  data->SetValue(i, it->second[i]);
329 #endif
330 
331  }
332  grid->GetPointData()->AddArray(data);
333  }
334  }
335 }
336 
337 
338 
339 // write out mesh data to the VTK file, this might come in handy to display
340 // boundary conditions and material data
341 // void VTKIO::meshdata_to_vtk(const MeshData& meshdata,
342 // vtkUnstructuredGrid* grid)
343 // {
344 // vtkPointData* pointdata = vtkPointData::New();
345 //
346 // const unsigned int n_vn = meshdata.n_val_per_node();
347 // const unsigned int n_dat = meshdata.n_node_data();
348 //
349 // pointdata->SetNumberOfTuples(n_dat);
350 // }
351 
352 #endif // LIBMESH_HAVE_VTK
353 
354 
355 
356 
357 // Constructor for reading
359  MeshInput<MeshBase> (mesh),
360  MeshOutput<MeshBase>(mesh),
361  _mesh_data(mesh_data),
362  _compress(false),
363  _local_node_map()
364 {
365  _vtk_grid = NULL;
366  libmesh_experimental();
367 }
368 
369 
370 
371 // Constructor for writing
372 VTKIO::VTKIO (const MeshBase& mesh, MeshData* mesh_data) :
373  MeshOutput<MeshBase>(mesh),
374  _mesh_data(mesh_data),
375  _compress(false),
376  _local_node_map()
377 {
378  _vtk_grid = NULL;
379  libmesh_experimental();
380 }
381 
382 
383 
384 vtkUnstructuredGrid* VTKIO::get_vtk_grid()
385 {
386  return _vtk_grid;
387 }
388 
389 
390 
392 {
393  this->_compress = b;
394 }
395 
396 
397 
398 void VTKIO::read (const std::string& name)
399 {
400  // This is a serial-only process for now;
401  // the Mesh should be read on processor 0 and
402  // broadcast later
403  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
404 
405  // Keep track of what kinds of elements this file contains
406  elems_of_dimension.clear();
407  elems_of_dimension.resize(4, false);
408 
409 #ifndef LIBMESH_HAVE_VTK
410  libMesh::err << "Cannot read VTK file: " << name
411  << "\nYou must have VTK installed and correctly configured to read VTK meshes."
412  << std::endl;
413  libmesh_error();
414 
415 #else
416  // Use a typedef, because these names are just crazy
417  typedef vtkSmartPointer<vtkXMLUnstructuredGridReader> MyReader;
418  MyReader reader = MyReader::New();
419 
420  // Pass the filename along to the reader
421  reader->SetFileName( name.c_str() );
422 
423  // Force reading
424  reader->Update();
425 
426  // read in the grid
427  _vtk_grid = reader->GetOutput();
428  // _vtk_grid->Update(); // FIXME: Necessary?
429 
430  // Get a reference to the mesh
432 
433  // Clear out any pre-existing data from the Mesh
434  mesh.clear();
435 
436  // Get the number of points from the _vtk_grid object
437  const unsigned int vtk_num_points = static_cast<unsigned int>(_vtk_grid->GetNumberOfPoints());
438 
439  // always numbered nicely??, so we can loop like this
440  // I'm pretty sure it is numbered nicely
441  for (unsigned int i=0; i<vtk_num_points; ++i)
442  {
443  // add to the id map
444  // and add the actual point
445  double * pnt = _vtk_grid->GetPoint(static_cast<vtkIdType>(i));
446  Point xyz(pnt[0], pnt[1], pnt[2]);
447  Node* newnode = mesh.add_point(xyz, i);
448 
449  // Add node to the nodes vector &
450  // tell the MeshData object the foreign node id.
451  if (this->_mesh_data != NULL)
452  this->_mesh_data->add_foreign_node_id (newnode, i);
453  }
454 
455  // Get the number of cells from the _vtk_grid object
456  const unsigned int vtk_num_cells = static_cast<unsigned int>(_vtk_grid->GetNumberOfCells());
457 
458  for (unsigned int i=0; i<vtk_num_cells; ++i)
459  {
460  vtkCell* cell = _vtk_grid->GetCell(i);
461  Elem* elem = NULL;
462  switch (cell->GetCellType())
463  {
464  case VTK_LINE:
465  elem = new Edge2;
466  break;
467  case VTK_QUADRATIC_EDGE:
468  elem = new Edge3;
469  break;
470  case VTK_TRIANGLE:
471  elem = new Tri3();
472  break;
473  case VTK_QUADRATIC_TRIANGLE:
474  elem = new Tri6();
475  break;
476  case VTK_QUAD:
477  elem = new Quad4();
478  break;
479  case VTK_QUADRATIC_QUAD:
480  elem = new Quad8();
481  break;
482 #if VTK_MAJOR_VERSION > 5 || (VTK_MAJOR_VERSION == 5 && VTK_MINOR_VERSION > 0)
483  case VTK_BIQUADRATIC_QUAD:
484  elem = new Quad9();
485  break;
486 #endif
487  case VTK_TETRA:
488  elem = new Tet4();
489  break;
490  case VTK_QUADRATIC_TETRA:
491  elem = new Tet10();
492  break;
493  case VTK_WEDGE:
494  elem = new Prism6();
495  break;
496  case VTK_QUADRATIC_WEDGE:
497  elem = new Prism15();
498  break;
499  case VTK_BIQUADRATIC_QUADRATIC_WEDGE:
500  elem = new Prism18();
501  break;
502  case VTK_HEXAHEDRON:
503  elem = new Hex8();
504  break;
505  case VTK_QUADRATIC_HEXAHEDRON:
506  elem = new Hex20();
507  break;
508  case VTK_TRIQUADRATIC_HEXAHEDRON:
509  elem = new Hex27();
510  break;
511  case VTK_PYRAMID:
512  elem = new Pyramid5();
513  break;
514  default:
515  libMesh::err << "element type not implemented in vtkinterface " << cell->GetCellType() << std::endl;
516  libmesh_error();
517  break;
518  }
519 
520  // get the straightforward numbering from the VTK cells
521  for (unsigned int j=0; j<elem->n_nodes(); ++j)
522  elem->set_node(j) = mesh.node_ptr(cell->GetPointId(j));
523 
524  // then get the connectivity
525  std::vector<dof_id_type> conn;
526  elem->connectivity(0, VTK, conn);
527 
528  // then reshuffle the nodes according to the connectivity, this
529  // two-time-assign would evade the definition of the vtk_mapping
530  for (unsigned int j=0; j<conn.size(); ++j)
531  elem->set_node(j) = mesh.node_ptr(conn[j]);
532 
533  elem->set_id(i);
534 
535  elems_of_dimension[elem->dim()] = true;
536 
537  mesh.add_elem(elem);
538  } // end loop over VTK cells
539 
540  // Set the mesh dimension to the largest encountered for an element
541  for (unsigned int i=0; i!=4; ++i)
542  if (elems_of_dimension[i])
543  mesh.set_mesh_dimension(i);
544 
545 #if LIBMESH_DIM < 3
546  if (mesh.mesh_dimension() > LIBMESH_DIM)
547  {
548  libMesh::err << "Cannot open dimension " <<
549  mesh.mesh_dimension() <<
550  " mesh file when configured without " <<
551  mesh.mesh_dimension() << "D support." <<
552  std::endl;
553  libmesh_error();
554  }
555 #endif
556 
557 #endif // LIBMESH_HAVE_VTK
558 }
559 
560 
561 
562 void VTKIO::write_nodal_data (const std::string& fname,
563 #ifdef LIBMESH_HAVE_VTK
564  const std::vector<Number>& soln,
565  const std::vector<std::string>& names
566 #else
567  const std::vector<Number>&,
568  const std::vector<std::string>&
569 #endif
570 )
571 {
572 #ifndef LIBMESH_HAVE_VTK
573 
574  libMesh::err << "Cannot write VTK file: " << fname
575  << "\nYou must have VTK installed and correctly configured to read VTK meshes."
576  << std::endl;
577  libmesh_error();
578 
579 #else
580 
582 
583  // Is this really important? If so, it should be more than an assert...
584  // libmesh_assert(fname.substr(fname.rfind("."), fname.size()) == ".pvtu");
585 
586  // we only use Unstructured grids
587  _vtk_grid = vtkUnstructuredGrid::New();
588  vtkSmartPointer<vtkXMLPUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLPUnstructuredGridWriter>::New();
589 
590  // add nodes to the grid and update _local_node_map
591  _local_node_map.clear();
592  this->nodes_to_vtk();
593 
594  // add cells to the grid
595  this->cells_to_vtk();
596 
597  // add nodal solutions to the grid, if solutions are given
598  if (names.size() > 0)
599  {
600  std::size_t num_vars = names.size();
601  dof_id_type num_nodes = mesh.n_nodes();
602 
603  for (std::size_t variable=0; variable<num_vars; ++variable)
604  {
605  vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
606  data->SetName(names[variable].c_str());
607 
608  // number of local and ghost nodes
609  data->SetNumberOfValues(_local_node_map.size());
610 
611  // loop over all nodes and get the solution for the current
612  // variable, if the node is in the current partition
613  for (dof_id_type k=0; k<num_nodes; ++k)
614  {
615  if (_local_node_map.find(k) == _local_node_map.end())
616  continue; // not a local node
617 
618  if (!soln.empty())
619  {
620 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
621  libmesh_do_once (libMesh::err << "Only writing the real part for complex numbers!\n"
622  << "if you need this support contact " << LIBMESH_PACKAGE_BUGREPORT
623  << std::endl);
624  data->SetValue(_local_node_map[k], soln[k*num_vars + variable].real());
625 #else
626  data->SetValue(_local_node_map[k], soln[k*num_vars + variable]);
627 #endif
628  }
629  else
630  {
631  data->SetValue(_local_node_map[k], 0);
632  }
633  }
634  _vtk_grid->GetPointData()->AddArray(data);
635  }
636  }
637 
638  // Tell the writer how many partitions exist and on which processor
639  // we are currently
640  writer->SetNumberOfPieces(MeshOutput<MeshBase>::mesh().n_processors());
641  writer->SetStartPiece(MeshOutput<MeshBase>::mesh().processor_id());
642  writer->SetEndPiece(MeshOutput<MeshBase>::mesh().processor_id());
643 
644  // partitions overlap by one node
645  // FIXME: According to this document
646  // http://paraview.org/Wiki/images/5/51/SC07_tut107_ParaView_Handouts.pdf
647  // the ghosts are cells rather than nodes.
648  writer->SetGhostLevel(1);
649 
650  writer->SetInput(_vtk_grid);
651  writer->SetFileName(fname.c_str());
652  writer->SetDataModeToAscii();
653 
654  // compress the output, if desired (switches also to binary)
655  if (this->_compress)
656  {
657 #if !VTK_VERSION_LESS_THAN(5,6,0)
658  writer->SetCompressorTypeToZLib();
659 #else
660  libmesh_do_once(libMesh::err << "Compression not implemented with old VTK libs!" << std::endl;);
661 #endif
662  }
663 
664  writer->Write();
665 
666  _vtk_grid->Delete();
667 #endif
668 }
669 
670 
671 
672 // Output the mesh without solutions to a .pvtu file
673 void VTKIO::write (const std::string& name)
674 {
675  std::vector<Number> soln;
676  std::vector<std::string> names;
677  this->write_nodal_data(name, soln, names);
678 }
679 
680 
681 
682 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:07 UTC

Hosted By:
SourceForge.net Logo