libMesh
vtk_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 // Local includes
20 #include "libmesh/libmesh_config.h"
21 #include "libmesh/vtk_io.h"
22 #include "libmesh/mesh_base.h"
23 #include "libmesh/equation_systems.h"
24 #include "libmesh/numeric_vector.h"
25 #include "libmesh/system.h"
26 #include "libmesh/node.h"
27 #include "libmesh/elem.h"
28 #include "libmesh/enum_io_package.h"
29 
30 #ifdef LIBMESH_HAVE_VTK
31 
32 // I get a lot of "warning: extra ';' inside a class [-Wextra-semi]" from clang
33 // on VTK header files.
34 #include "libmesh/ignore_warnings.h"
35 
36 #include "vtkXMLUnstructuredGridReader.h"
37 #include "vtkXMLPUnstructuredGridReader.h"
38 #include "vtkXMLUnstructuredGridWriter.h"
39 #include "vtkXMLPUnstructuredGridWriter.h"
40 #include "vtkUnstructuredGrid.h"
41 #include "vtkIntArray.h"
42 #include "vtkCellArray.h"
43 #include "vtkCellData.h"
44 #include "vtkDoubleArray.h"
45 #include "vtkGenericCell.h"
46 #include "vtkPointData.h"
47 #include "vtkPoints.h"
48 #include "vtkSmartPointer.h"
49 
50 #ifdef LIBMESH_HAVE_MPI
51 #include "vtkMPI.h"
52 #include "vtkMPICommunicator.h"
53 #include "vtkMPIController.h"
54 #endif
55 
56 #include "libmesh/restore_warnings.h"
57 
58 // C++ includes
59 #include <fstream>
60 
61 
62 // A convenient macro for comparing VTK versions. Returns 1 if the
63 // current VTK version is < major.minor.subminor and zero otherwise.
64 //
65 // It relies on the VTK version numbers detected during configure. Note that if
66 // LIBMESH_HAVE_VTK is not defined, none of the LIBMESH_DETECTED_VTK_VERSION_* variables will
67 // be defined either.
68 #define VTK_VERSION_LESS_THAN(major,minor,subminor) \
69  ((LIBMESH_DETECTED_VTK_VERSION_MAJOR < (major) || \
70  (LIBMESH_DETECTED_VTK_VERSION_MAJOR == (major) && (LIBMESH_DETECTED_VTK_VERSION_MINOR < (minor) || \
71  (LIBMESH_DETECTED_VTK_VERSION_MINOR == (minor) && \
72  LIBMESH_DETECTED_VTK_VERSION_SUBMINOR < (subminor))))) ? 1 : 0)
73 
74 #endif // LIBMESH_HAVE_VTK
75 
76 
77 
78 namespace libMesh
79 {
80 
81 // Constructor for reading
83  MeshInput<MeshBase> (mesh, /*is_parallel_format=*/true),
84  MeshOutput<MeshBase>(mesh, /*is_parallel_format=*/true)
85 #ifdef LIBMESH_HAVE_VTK
86  ,_compress(false)
87 #endif
88 {
89 }
90 
91 
92 
93 // Constructor for writing
94 VTKIO::VTKIO (const MeshBase & mesh) :
95  MeshOutput<MeshBase>(mesh, /*is_parallel_format=*/true)
96 #ifdef LIBMESH_HAVE_VTK
97  ,_compress(false)
98 #endif
99 {
100 }
101 
102 
103 
104 // Output the mesh without solutions to a .pvtu file
105 void VTKIO::write (const std::string & name)
106 {
107  std::vector<Number> soln;
108  std::vector<std::string> names;
109  this->write_nodal_data(name, soln, names);
110 }
111 
112 
113 
114 // The rest of the file is wrapped in ifdef LIBMESH_HAVE_VTK except for
115 // a couple of "stub" functions at the bottom.
116 #ifdef LIBMESH_HAVE_VTK
117 
118 // Initialize the static _element_maps map.
119 std::map<ElemMappingType, VTKIO::ElementMaps> VTKIO::_element_maps = VTKIO::build_element_maps();
120 
121 // Static function which constructs the ElementMaps object.
122 std::map<ElemMappingType, VTKIO::ElementMaps> VTKIO::build_element_maps()
123 {
124  // Object to be filled up
125  std::map<ElemMappingType, VTKIO::ElementMaps> all_maps;
126  ElementMaps em; // Lagrange element maps
127 
128  em.associate(EDGE2, VTK_LINE);
129  em.associate(EDGE3, VTK_QUADRATIC_EDGE);
130  em.associate(TRI3, VTK_TRIANGLE);
131  em.associate(TRI6, VTK_QUADRATIC_TRIANGLE);
132  em.associate(QUAD4, VTK_QUAD);
133  em.associate(QUAD8, VTK_QUADRATIC_QUAD);
134  em.associate(TET4, VTK_TETRA);
135  em.associate(TET10, VTK_QUADRATIC_TETRA);
136  em.associate(HEX8, VTK_HEXAHEDRON);
137  em.associate(HEX20, VTK_QUADRATIC_HEXAHEDRON);
138  em.associate(HEX27, VTK_TRIQUADRATIC_HEXAHEDRON);
139  em.associate(PRISM6, VTK_WEDGE);
140  em.associate(PRISM15, VTK_QUADRATIC_WEDGE);
141  em.associate(PRISM18, VTK_BIQUADRATIC_QUADRATIC_WEDGE);
142  em.associate(PYRAMID5, VTK_PYRAMID);
143 
144  // VTK_BIQUADRATIC_QUAD has been around since VTK 5.0
145 #if VTK_MAJOR_VERSION > 5 || (VTK_MAJOR_VERSION == 5 && VTK_MINOR_VERSION > 0)
146  em.associate(QUAD9, VTK_BIQUADRATIC_QUAD);
147 #endif
148 
149  // TRI3SUBDIVISION is for writing only
150  em.writing_map[TRI3SUBDIVISION] = VTK_TRIANGLE;
151  all_maps[ElemMappingType::LAGRANGE_MAP] = em;
152 
153 
154  // VTK_BEZIER_* types were introduced in VTK 9.0
155 #if VTK_VERSION_LESS_THAN(9,0,0)
156  // Revert back to previous behavior when using an older version of VTK
158 #else
159  ElementMaps bem; // Rational Bernstein element maps
160  bem.associate(EDGE2, VTK_LINE);
161  bem.associate(EDGE3, VTK_BEZIER_CURVE);
162  bem.associate(TRI3, VTK_TRIANGLE);
163  bem.associate(TRI6, VTK_BEZIER_TRIANGLE);
164  bem.associate(QUAD4, VTK_QUAD);
165  bem.associate(QUAD8, VTK_QUADRATIC_QUAD);
166  bem.associate(QUAD9, VTK_BEZIER_QUADRILATERAL);
167  bem.associate(TET4, VTK_TETRA);
168  bem.associate(TET10, VTK_QUADRATIC_TETRA);
169  bem.associate(HEX8, VTK_HEXAHEDRON);
170  bem.associate(HEX20, VTK_QUADRATIC_HEXAHEDRON);
171  bem.associate(HEX27, VTK_BEZIER_HEXAHEDRON);
172  bem.associate(PRISM6, VTK_WEDGE);
173  bem.associate(PRISM15, VTK_QUADRATIC_WEDGE);
174  bem.associate(PRISM18, VTK_BEZIER_WEDGE);
175  bem.associate(PYRAMID5, VTK_PYRAMID);
176  bem.writing_map[TRI3SUBDIVISION] = VTK_TRIANGLE;
178 #endif
179 
180  return all_maps;
181 }
182 
183 
184 
185 void VTKIO::read (const std::string & name)
186 {
187  // This is a serial-only process for now;
188  // the Mesh should be read on processor 0 and
189  // broadcast later
190  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
191 
192  // Keep track of what kinds of elements this file contains
193  elems_of_dimension.clear();
194  elems_of_dimension.resize(4, false);
195 
196  // Use a typedef, because these names are just crazy
197  typedef vtkSmartPointer<vtkXMLPUnstructuredGridReader> MyReader;
198  MyReader reader = MyReader::New();
199 
200  // Pass the filename along to the reader
201  reader->SetFileName(name.c_str());
202 
203  // Force reading
204  reader->Update();
205 
206  // read in the grid
207  _vtk_grid = reader->GetOutput();
208 
209  // Get a reference to the mesh
211 
212  // Clear out any pre-existing data from the Mesh
213  mesh.clear();
214 
215  // Try to preserve any libMesh ids and subdomain ids we find in the
216  // file. This will be null if there are none, e.g. if a non-libMesh
217  // program wrote this file.
218 
219  vtkAbstractArray * abstract_elem_id =
220  _vtk_grid->GetCellData()->GetAbstractArray("libmesh_elem_id");
221  vtkAbstractArray * abstract_node_id =
222  _vtk_grid->GetPointData()->GetAbstractArray("libmesh_node_id");
223  vtkAbstractArray * abstract_subdomain_id =
224  _vtk_grid->GetCellData()->GetAbstractArray("subdomain_id");
225 
226  // Get ids as integers. This will be null if they are another data
227  // type, e.g. if a non-libMesh program used the names we thought
228  // were unique for different data.
229  vtkIntArray * elem_id = vtkIntArray::SafeDownCast(abstract_elem_id);
230  vtkIntArray * node_id = vtkIntArray::SafeDownCast(abstract_node_id);
231  vtkIntArray * subdomain_id = vtkIntArray::SafeDownCast(abstract_subdomain_id);
232 
233  if (abstract_elem_id && !elem_id)
234  libmesh_warning("Found non-integral libmesh_elem_id array; forced to ignore it.\n"
235  "This is technically valid but probably broken.");
236 
237  if (abstract_node_id && !node_id)
238  libmesh_warning("Found non-integral libmesh_node_id array; forced to ignore it.\n"
239  "This is technically valid but probably broken.");
240 
241  if (abstract_subdomain_id && !subdomain_id)
242  libmesh_warning("Found non-integral subdomain_id array; forced to ignore it.\n"
243  "This is technically valid but probably broken.");
244 
245  // Get the number of points from the _vtk_grid object
246  const unsigned int vtk_num_points = static_cast<unsigned int>(_vtk_grid->GetNumberOfPoints());
247 
248  // Map from VTK indexing to libMesh id if necessary
249  std::vector<dof_id_type> vtk_node_to_libmesh;
250  if (node_id)
251  vtk_node_to_libmesh.resize(vtk_num_points);
252 
253  // always numbered nicely so we can loop like this
254  for (unsigned int i=0; i<vtk_num_points; ++i)
255  {
256  // add to the id map
257  // and add the actual point
258  double pnt[3];
259  _vtk_grid->GetPoint(static_cast<vtkIdType>(i), pnt);
260  Point xyz(pnt[0], pnt[1], pnt[2]);
261 
262  if (node_id)
263  {
264  auto id = node_id->GetValue(i);
265 
266  // It would nice to distinguish between "duplicate nodes
267  // because one was ghosted in a parallel file segment" and
268  // "duplicate nodes because there was a bug", but I'm not
269  // sure how to do that with vtkXMLPUnstructuredGridReader
270  if (!mesh.query_node_ptr(id))
271  mesh.add_point(xyz, id);
272  vtk_node_to_libmesh[i] = id;
273  }
274  else
275  mesh.add_point(xyz, i);
276  }
277 
278  // Get the number of cells from the _vtk_grid object
279  const unsigned int vtk_num_cells = static_cast<unsigned int>(_vtk_grid->GetNumberOfCells());
280 
281  auto& element_map = libmesh_map_find(_element_maps, mesh.default_mapping_type());
282 
283  vtkSmartPointer<vtkGenericCell> cell = vtkSmartPointer<vtkGenericCell>::New();
284  for (unsigned int i=0; i<vtk_num_cells; ++i)
285  {
286  _vtk_grid->GetCell(i, cell);
287 
288  // Get the libMesh element type corresponding to this VTK element type.
289  ElemType libmesh_elem_type = element_map.find(cell->GetCellType());
290  auto elem = Elem::build(libmesh_elem_type);
291 
292  // get the straightforward numbering from the VTK cells
293  for (auto j : elem->node_index_range())
294  {
295  const auto vtk_point_id = cell->GetPointId(j);
296  const dof_id_type libmesh_node_id = node_id ?
297  vtk_node_to_libmesh[vtk_point_id] : vtk_point_id;
298 
299  elem->set_node(j, mesh.node_ptr(libmesh_node_id));
300  }
301 
302  // then get the connectivity
303  std::vector<dof_id_type> conn;
304  elem->connectivity(0, VTK, conn);
305 
306  // then reshuffle the nodes according to the connectivity, this
307  // two-time-assign would evade the definition of the vtk_mapping
308  for (unsigned int j=0,
309  n_conn = cast_int<unsigned int>(conn.size());
310  j != n_conn; ++j)
311  elem->set_node(j, mesh.node_ptr(conn[j]));
312 
313  if (elem_id)
314  {
315  auto id = elem_id->GetValue(i);
316  libmesh_error_msg_if
317  (mesh.query_elem_ptr(id), "Duplicate element id " << id <<
318  " found in libmesh_elem_ids");
319  elem->set_id(id);
320  }
321  else
322  elem->set_id(i);
323 
324  if (subdomain_id)
325  {
326  auto sbdid = subdomain_id->GetValue(i);
327  elem->subdomain_id() = sbdid;
328  }
329 
330  elems_of_dimension[elem->dim()] = true;
331 
332  mesh.add_elem(std::move(elem));
333  } // end loop over VTK cells
334 
335  // Set the mesh dimension to the largest encountered for an element
336  for (unsigned char i=0; i!=4; ++i)
337  if (elems_of_dimension[i])
339 
340 #if LIBMESH_DIM < 3
341  libmesh_error_msg_if(mesh.mesh_dimension() > LIBMESH_DIM,
342  "Cannot open dimension "
343  << mesh.mesh_dimension()
344  << " mesh file when configured without "
345  << mesh.mesh_dimension()
346  << "D support.");
347 #endif // LIBMESH_DIM < 3
348 }
349 
350 
351 
352 void VTKIO::write_nodal_data (const std::string & fname,
353  const std::vector<Number> & soln,
354  const std::vector<std::string> & names)
355 {
356  // Warn that the .pvtu file extension should be used. Paraview
357  // recognizes this, and it works in both serial and parallel. Only
358  // warn about this once.
359  if (fname.substr(fname.rfind("."), fname.size()) != ".pvtu")
360  libmesh_do_once(libMesh::err << "The .pvtu extension should be used when writing VTK files in libMesh.");
361 
362  // If there are variable names being written, the solution vector
363  // should not be empty, it should have been broadcast to all
364  // processors by the MeshOutput base class, since VTK is a parallel
365  // format. Verify this before going further.
366  libmesh_error_msg_if(!names.empty() && soln.empty(),
367  "Empty soln vector in VTKIO::write_nodal_data().");
368 
369  // Get a reference to the mesh
371 
372  // we only use Unstructured grids
373  _vtk_grid = vtkSmartPointer<vtkUnstructuredGrid>::New();
374  vtkSmartPointer<vtkXMLPUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLPUnstructuredGridWriter>::New();
375 #ifdef LIBMESH_HAVE_MPI
376  // Set VTK to the same communicator as libMesh
377  vtkSmartPointer<vtkMPICommunicator> vtk_comm = vtkSmartPointer<vtkMPICommunicator>::New();
378  MPI_Comm mpi_comm = mesh.comm().get();
379  vtkMPICommunicatorOpaqueComm vtk_opaque_comm(&mpi_comm);
380  vtk_comm->InitializeExternal(&vtk_opaque_comm);
381 
382  vtkSmartPointer<vtkMPIController> vtk_mpi_ctrl = vtkSmartPointer<vtkMPIController>::New();
383  vtk_mpi_ctrl->SetCommunicator(vtk_comm);
384 
385  writer->SetController(vtk_mpi_ctrl);
386 #endif
387 
388  // add nodes to the grid and update _local_node_map
389  _local_node_map.clear();
390  this->nodes_to_vtk();
391 
392  // add cells to the grid
393  this->cells_to_vtk();
394 
395  // add nodal solutions to the grid, if solutions are given
396  if (names.size() > 0)
397  {
398  std::size_t num_vars = names.size();
399  std::vector<Number> local_values;
400 
401 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
402  std::vector<Real> local_real_values;
403 #endif
404 
405  for (std::size_t variable=0; variable<num_vars; ++variable)
406  {
407  get_local_node_values(local_values, variable, soln, names);
408 
409 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
410  // write real part
411  local_real_values.resize(local_values.size());
412  std::transform(local_values.begin(), local_values.end(),
413  local_real_values.begin(),
414  [](Number x) { return x.real(); });
415  node_values_to_vtk(names[variable] + "_real", local_real_values);
416 
417  // write imaginary part
418  local_real_values.resize(local_values.size());
419  std::transform(local_values.begin(), local_values.end(),
420  local_real_values.begin(),
421  [](Number x) { return x.imag(); });
422  node_values_to_vtk(names[variable] + "_imag", local_real_values);
423 #else
424  node_values_to_vtk(names[variable], local_values);
425 #endif
426  }
427  }
428 
429  // Tell the writer how many partitions exist and on which processor
430  // we are currently
431  writer->SetNumberOfPieces(mesh.n_processors());
432  writer->SetStartPiece(mesh.processor_id());
433  writer->SetEndPiece(mesh.processor_id());
434 
435  // partitions overlap by one node
436  // FIXME: According to this document
437  // http://paraview.org/Wiki/images/5/51/SC07_tut107_ParaView_Handouts.pdf
438  // the ghosts are cells rather than nodes.
439  writer->SetGhostLevel(1);
440 
441  // VTK 6 replaces SetInput() with SetInputData(). See
442  // http://www.vtk.org/Wiki/VTK/VTK_6_Migration/Replacement_of_SetInput
443  // for the full explanation.
444 #if VTK_VERSION_LESS_THAN(6,0,0)
445  writer->SetInput(_vtk_grid);
446 #else
447  writer->SetInputData(_vtk_grid);
448 #endif
449 
450  writer->SetFileName(fname.c_str());
451  writer->SetDataModeToAscii();
452 
453  // compress the output, if desired (switches also to binary)
454  if (this->_compress)
455  {
456 #if !VTK_VERSION_LESS_THAN(5,6,0)
457  writer->SetCompressorTypeToZLib();
458  writer->SetDataModeToBinary();
459 #else
460  libmesh_do_once(libMesh::err << "Compression not implemented with old VTK libs!" << std::endl;);
461 #endif
462  }
463 
464  writer->Write();
465 
466 }
467 
468 
469 
470 vtkUnstructuredGrid * VTKIO::get_vtk_grid()
471 {
472  return _vtk_grid;
473 }
474 
475 
476 
478 {
479  this->_compress = b;
480 }
481 
482 
483 
485 {
487 
488  // containers for points and coordinates of points
489  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
490  vtkSmartPointer<vtkDoubleArray> pcoords = vtkSmartPointer<vtkDoubleArray>::New();
491  // if this grid is to be used in VTK then the dimension of the points should be 3
492  pcoords->SetNumberOfComponents(LIBMESH_DIM);
493  pcoords->Allocate(3*mesh.n_local_nodes());
494  points->SetNumberOfPoints(mesh.n_local_nodes()); // it seems that it needs this to prevent a segfault
495 
496  // SetRationalWeights() was introduced in VTK 9.0
497 #if !VTK_VERSION_LESS_THAN(9,0,0)
498  bool have_weights = false;
499  int weight_index = 0;
500  vtkSmartPointer<vtkDoubleArray> rational_weights;
501 
503  {
504  rational_weights = vtkSmartPointer<vtkDoubleArray>::New();
505  rational_weights->SetName("RationalWeights");
506  rational_weights->SetNumberOfComponents(1);
507  weight_index = static_cast<int>(mesh.default_mapping_data());
508  have_weights = true;
509  }
510 #endif
511 
512  vtkSmartPointer<vtkIntArray> node_id = vtkSmartPointer<vtkIntArray>::New();
513  node_id->SetName("libmesh_node_id");
514  node_id->SetNumberOfComponents(1);
515 
516  unsigned int local_node_counter = 0;
517 
518  for (const auto & node_ptr : mesh.local_node_ptr_range())
519  {
520  const Node & node = *node_ptr;
521 
522  double pnt[3] = {0, 0, 0};
523  for (unsigned int i=0; i<LIBMESH_DIM; ++i)
524  pnt[i] = double(node(i));
525 
526  // Fill mapping between global and local node numbers
527  _local_node_map[node.id()] = local_node_counter;
528 
529  // add point
530 #if VTK_VERSION_LESS_THAN(7,1,0)
531  pcoords->InsertNextTupleValue(pnt);
532 #else
533  pcoords->InsertNextTuple(pnt);
534 #endif
535 #if !VTK_VERSION_LESS_THAN(9,0,0)
536  if (have_weights)
537  {
538  Real weight = node.get_extra_datum<Real>(weight_index);
539  rational_weights->InsertTuple1(local_node_counter, double(weight));
540  }
541 #endif
542 
543  node_id->InsertTuple1(local_node_counter, node.id());
544 
545  ++local_node_counter;
546  }
547 
548  // add coordinates to points
549  points->SetData(pcoords);
550 
551  // add points to grid
552  _vtk_grid->SetPoints(points);
553  _vtk_grid->GetPointData()->AddArray(node_id);
554 
555 #if !VTK_VERSION_LESS_THAN(9,0,0)
556  if (have_weights)
557  _vtk_grid->GetPointData()->SetRationalWeights(rational_weights);
558 
559 #endif
560 }
561 
562 
563 
565 {
567 
568  auto& element_map = libmesh_map_find(_element_maps, mesh.default_mapping_type());
569  vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
570  vtkSmartPointer<vtkIdList> pts = vtkSmartPointer<vtkIdList>::New();
571 
572  std::vector<int> types(mesh.n_active_local_elem());
573 
574  // We already created this but we need to add more if we have any
575  // ghost nodes
576  vtkAbstractArray * abstract_node_id =
577  _vtk_grid->GetPointData()->GetAbstractArray("libmesh_node_id");
578  vtkIntArray * node_id = vtkIntArray::SafeDownCast(abstract_node_id);
579  libmesh_assert(node_id);
580 
581  vtkSmartPointer<vtkIntArray> elem_id = vtkSmartPointer<vtkIntArray>::New();
582  elem_id->SetName("libmesh_elem_id");
583  elem_id->SetNumberOfComponents(1);
584 
585  vtkSmartPointer<vtkIntArray> subdomain_id = vtkSmartPointer<vtkIntArray>::New();
586  subdomain_id->SetName("subdomain_id");
587  subdomain_id->SetNumberOfComponents(1);
588 
589  vtkSmartPointer<vtkIntArray> elem_proc_id = vtkSmartPointer<vtkIntArray>::New();
590  elem_proc_id->SetName("processor_id");
591  elem_proc_id->SetNumberOfComponents(1);
592 
593  unsigned active_element_counter = 0;
594  for (const auto & elem : mesh.active_local_element_ptr_range())
595  {
596  // When using rational bernstein these hold the weights
597  if ( elem->type() == NODEELEM )
598  continue;
599 
600  pts->SetNumberOfIds(elem->n_nodes());
601 
602  // get the connectivity for this element
603  std::vector<dof_id_type> conn;
604  elem->connectivity(0, VTK, conn);
605 
606  for (unsigned int i=0,
607  n_conn = cast_int<unsigned int>(conn.size());
608  i != n_conn; ++i)
609  {
610  // If the node ID is not found in the _local_node_map, we'll
611  // add it to the _vtk_grid. NOTE[JWP]: none of the examples
612  // I have actually enters this section of code...
613  if (!_local_node_map.count(conn[i]))
614  {
615  dof_id_type global_node_id = elem->node_id(i);
616 
617  const Point & the_node = mesh.point(global_node_id);
618 
619  // InsertNextPoint accepts either a double or float array of length 3.
620  double pt[3] = {0., 0., 0.};
621  for (unsigned int d=0; d<LIBMESH_DIM; ++d)
622  pt[d] = double(the_node(d));
623 
624  // Insert the point into the _vtk_grid
625  vtkIdType local = _vtk_grid->GetPoints()->InsertNextPoint(pt);
626 
627  // Update the _local_node_map with the ID returned by VTK
628  _local_node_map[global_node_id] =
629  cast_int<dof_id_type>(local);
630 
631  node_id->InsertTuple1(local, global_node_id);
632  }
633 
634  // Otherwise, the node ID was found in the _local_node_map, so
635  // insert it into the vtkIdList.
636  pts->InsertId(i, _local_node_map[conn[i]]);
637  }
638 
639  vtkIdType vtkcellid = cells->InsertNextCell(pts);
640  types[active_element_counter] = cast_int<int>(element_map.find(elem->type()));
641 
642  elem_id->InsertTuple1(vtkcellid, elem->id());
643  subdomain_id->InsertTuple1(vtkcellid, elem->subdomain_id());
644  elem_proc_id->InsertTuple1(vtkcellid, elem->processor_id());
645  ++active_element_counter;
646  } // end loop over active elements
647 
648  _vtk_grid->SetCells(types.data(), cells);
649  _vtk_grid->GetCellData()->AddArray(elem_id);
650  _vtk_grid->GetCellData()->AddArray(subdomain_id);
651  _vtk_grid->GetCellData()->AddArray(elem_proc_id);
652 }
653 
654 void VTKIO::node_values_to_vtk(const std::string & name,
655  const std::vector<Real> & local_values)
656 {
657  vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
658  data->SetName(name.c_str());
659 
660  libmesh_assert_equal_to(_local_node_map.size(), local_values.size());
661 
662  // number of local and ghost nodes
663  data->SetNumberOfValues(_local_node_map.size());
664 
665  // copy values into vtk
666  for (auto i : index_range(local_values)) {
667  data->SetValue(i, double(local_values[i]));
668  }
669 
670  _vtk_grid->GetPointData()->AddArray(data);
671 }
672 
673 void VTKIO::get_local_node_values(std::vector<Number> & local_values,
674  std::size_t variable,
675  const std::vector<Number> & soln,
676  const std::vector<std::string> & names)
677 {
679  std::size_t num_vars = names.size();
680  dof_id_type num_nodes = mesh.n_nodes();
681 
682  local_values.clear();
683  local_values.resize(_local_node_map.size(), 0.0);
684 
685  // loop over all nodes and get the solution for the current
686  // variable, if the node is in the current partition
687  for (dof_id_type k=0; k<num_nodes; ++k)
688  if (const auto local_node_it = _local_node_map.find(k);
689  local_node_it != _local_node_map.end())
690  local_values[local_node_it->second] = soln[k*num_vars + variable];
691 }
692 
693 
694 
702 // void VTKIO::system_vectors_to_vtk(const EquationSystems & es,
703 // vtkUnstructuredGrid *& grid)
704 // {
705 // if (MeshOutput<MeshBase>::mesh().processor_id() == 0)
706 // {
707 // std::map<std::string, std::vector<Number>> vecs;
708 // for (unsigned int i=0; i<es.n_systems(); ++i)
709 // {
710 // const System & sys = es.get_system(i);
711 // System::const_vectors_iterator v_end = sys.vectors_end();
712 // System::const_vectors_iterator it = sys.vectors_begin();
713 // for (; it!= v_end; ++it)
714 // {
715 // // for all vectors on this system
716 // std::vector<Number> values;
717 // // libMesh::out<<"it "<<it->first<<std::endl;
718 //
719 // it->second->localize_to_one(values, 0);
720 // // libMesh::out<<"finish localize"<<std::endl;
721 // vecs[it->first] = values;
722 // }
723 // }
724 //
725 // std::map<std::string, std::vector<Number>>::iterator it = vecs.begin();
726 //
727 // for (; it!=vecs.end(); ++it)
728 // {
729 // vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
730 // data->SetName(it->first.c_str());
731 // libmesh_assert_equal_to (it->second.size(), es.get_mesh().n_nodes());
732 // data->SetNumberOfValues(it->second.size());
733 //
734 // for (auto i : index_range(it->second))
735 // {
736 // #ifdef LIBMESH_USE_COMPLEX_NUMBERS
737 // libmesh_do_once (libMesh::err << "Only writing the real part for complex numbers!\n"
738 // << "if you need this support contact " << LIBMESH_PACKAGE_BUGREPORT
739 // << std::endl);
740 // data->SetValue(i, it->second[i].real());
741 // #else
742 // data->SetValue(i, it->second[i]);
743 // #endif
744 //
745 // }
746 // grid->GetPointData()->AddArray(data);
747 // }
748 // }
749 // }
750 
751 
752 
753 #else // !LIBMESH_HAVE_VTK
754 
755 void VTKIO::read (const std::string & name)
756 {
757  libmesh_error_msg("Cannot read VTK file: " << name \
758  << "\nYou must have VTK installed and correctly configured to read VTK meshes.");
759 }
760 
761 
762 
763 void VTKIO::write_nodal_data (const std::string & fname,
764  const std::vector<Number> &,
765  const std::vector<std::string> &)
766 {
767  libmesh_error_msg("Cannot write VTK file: " << fname \
768  << "\nYou must have VTK installed and correctly configured to read VTK meshes.");
769 }
770 
771 
772 #endif // LIBMESH_HAVE_VTK
773 
774 
775 
776 } // namespace libMesh
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
void set_compression(bool b)
Setter for compression flag.
Definition: vtk_io.C:477
OStreamProxy err
const MT & mesh() const
Definition: mesh_output.h:259
void node_values_to_vtk(const std::string &name, const std::vector< Real > &local_values)
write the nodal values of soln to a vtkUnstructuredGrid
Definition: vtk_io.C:654
static std::map< ElemMappingType, ElementMaps > _element_maps
ElementMaps objects that are built statically and used by all instances of this class.
Definition: vtk_io.h:215
ElemType
Defines an enum for geometric element types.
std::map< dof_id_type, dof_id_type > _local_node_map
maps global node id to node id of partition
Definition: vtk_io.h:179
A Node is like a Point, but with more information.
Definition: node.h:52
std::vector< bool > elems_of_dimension
A vector of bools describing what dimension elements have been encountered when reading a mesh...
Definition: mesh_input.h:107
vtkSmartPointer< vtkUnstructuredGrid > _vtk_grid
Write the system vectors to vtk.
Definition: vtk_io.h:169
MeshBase & mesh
bool _compress
Flag to indicate whether the output should be compressed.
Definition: vtk_io.h:174
virtual void write(const std::string &) override
Output the mesh without solutions to a .pvtu file.
dof_id_type n_local_nodes() const
Definition: mesh_base.h:442
const Parallel::Communicator & comm() const
This class defines an abstract interface for Mesh output.
Definition: mesh_output.h:53
The libMesh namespace provides an interface to certain functionality in the library.
Helper object that holds a map from VTK to libMesh element types and vice-versa.
Definition: vtk_io.h:186
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
This is the MeshBase class.
Definition: mesh_base.h:75
ElemMappingType default_mapping_type() const
Returns the default master space to physical space mapping basis functions to be used on newly added ...
Definition: mesh_base.h:812
virtual void read(const std::string &) override
This method implements reading a mesh from a specified file in VTK format.
Definition: vtk_io.C:185
static std::map< ElemMappingType, ElementMaps > build_element_maps()
Static function used to construct _element_maps.
Definition: vtk_io.C:122
void nodes_to_vtk()
write the nodes from the mesh into a vtkUnstructuredGrid and update the local_node_map.
Definition: vtk_io.C:484
dof_id_type n_active_local_elem() const
Definition: mesh_base.h:565
processor_id_type n_processors() const
dof_id_type id() const
Definition: dof_object.h:828
This class defines an abstract interface for Mesh input.
Definition: mesh_base.h:57
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:437
unsigned char default_mapping_data() const
Returns any default data value used by the master space to physical space mapping.
Definition: mesh_base.h:830
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:444
virtual const Node * query_node_ptr(const dof_id_type i) const =0
libmesh_assert(ctx)
virtual void write_nodal_data(const std::string &, const std::vector< Number > &, const std::vector< std::string > &) override
This method implements writing a mesh with nodal data to a specified file where the nodal data and va...
Definition: vtk_io.C:352
vtkUnstructuredGrid * get_vtk_grid()
Get a pointer to the VTK unstructured grid data structure.
Definition: vtk_io.C:470
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
Definition: mesh_base.h:275
std::map< ElemType, vtkIdType > writing_map
Definition: vtk_io.h:207
virtual void clear()
Deletes all the element and node data that is currently stored.
Definition: mesh_base.C:920
communicator & get()
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
VTKIO(MeshBase &mesh)
Constructor.
Definition: vtk_io.C:82
void associate(ElemType libmesh_type, vtkIdType vtk_type)
Definition: vtk_io.h:189
virtual const Elem * query_elem_ptr(const dof_id_type i) const =0
void get_local_node_values(std::vector< Number > &local_values, std::size_t variable, const std::vector< Number > &soln, const std::vector< std::string > &names)
Extract the values of soln that correspond to the nodes.
Definition: vtk_io.C:673
void cells_to_vtk()
write the cells from the mesh into a vtkUnstructuredGrid
Definition: vtk_io.C:564
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
virtual const Point & point(const dof_id_type i) const =0
T get_extra_datum(const unsigned int index) const
Gets the value on this object of the extra datum associated with index, which should have been obtain...
Definition: dof_object.h:1146
virtual const Node * node_ptr(const dof_id_type i) const =0
processor_id_type processor_id() const
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
virtual dof_id_type n_nodes() const =0
uint8_t dof_id_type
Definition: id_types.h:67