Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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/cell_c0polyhedron.h"
29 : #include "libmesh/face_c0polygon.h"
30 : #include "libmesh/face_polygon.h"
31 : #include "libmesh/enum_io_package.h"
32 : #include "libmesh/utility.h"
33 :
34 : #ifdef LIBMESH_HAVE_VTK
35 :
36 : // I get a lot of "warning: extra ';' inside a class [-Wextra-semi]" from clang
37 : // on VTK header files.
38 : #include "libmesh/ignore_warnings.h"
39 :
40 : #include "vtkXMLUnstructuredGridReader.h"
41 : #include "vtkXMLPUnstructuredGridReader.h"
42 : #include "vtkXMLUnstructuredGridWriter.h"
43 : #include "vtkXMLPUnstructuredGridWriter.h"
44 : #include "vtkXMLFileReadTester.h"
45 : #include "vtkUnstructuredGrid.h"
46 : #include "vtkIntArray.h"
47 : #include "vtkCellArray.h"
48 : #include "vtkCellData.h"
49 : #include "vtkDoubleArray.h"
50 : #include "vtkGenericCell.h"
51 : #include "vtkIdList.h"
52 : #include "vtkPointData.h"
53 : #include "vtkPoints.h"
54 : #include "vtkSmartPointer.h"
55 :
56 : #ifdef LIBMESH_HAVE_MPI
57 : #include "vtkMPI.h"
58 : #include "vtkMPICommunicator.h"
59 : #include "vtkMPIController.h"
60 : #endif
61 :
62 : #include "libmesh/restore_warnings.h"
63 :
64 : // C++ includes
65 : #include <cstring>
66 : #include <fstream>
67 :
68 :
69 : // A convenient macro for comparing VTK versions. Returns 1 if the
70 : // current VTK version is < major.minor.subminor and zero otherwise.
71 : //
72 : // It relies on the VTK version numbers detected during configure. Note that if
73 : // LIBMESH_HAVE_VTK is not defined, none of the LIBMESH_DETECTED_VTK_VERSION_* variables will
74 : // be defined either.
75 : #define VTK_VERSION_LESS_THAN(major,minor,subminor) \
76 : ((LIBMESH_DETECTED_VTK_VERSION_MAJOR < (major) || \
77 : (LIBMESH_DETECTED_VTK_VERSION_MAJOR == (major) && (LIBMESH_DETECTED_VTK_VERSION_MINOR < (minor) || \
78 : (LIBMESH_DETECTED_VTK_VERSION_MINOR == (minor) && \
79 : LIBMESH_DETECTED_VTK_VERSION_SUBMINOR < (subminor))))) ? 1 : 0)
80 :
81 : #endif // LIBMESH_HAVE_VTK
82 :
83 :
84 :
85 : namespace libMesh
86 : {
87 :
88 : // Constructor for reading
89 0 : VTKIO::VTKIO (MeshBase & mesh) :
90 : MeshInput<MeshBase> (mesh, /*is_parallel_format=*/true),
91 0 : MeshOutput<MeshBase>(mesh, /*is_parallel_format=*/true)
92 : #ifdef LIBMESH_HAVE_VTK
93 : ,_compress(false)
94 : #endif
95 : {
96 0 : }
97 :
98 :
99 :
100 : // Constructor for writing
101 0 : VTKIO::VTKIO (const MeshBase & mesh) :
102 0 : MeshOutput<MeshBase>(mesh, /*is_parallel_format=*/true)
103 : #ifdef LIBMESH_HAVE_VTK
104 : ,_compress(false)
105 : #endif
106 : {
107 0 : }
108 :
109 :
110 :
111 : // Output the mesh without solutions to a .pvtu file
112 0 : void VTKIO::write (const std::string & name)
113 : {
114 0 : std::vector<Number> soln;
115 0 : std::vector<std::string> names;
116 0 : this->write_nodal_data(name, soln, names);
117 0 : }
118 :
119 :
120 :
121 : // The rest of the file is wrapped in ifdef LIBMESH_HAVE_VTK except for
122 : // a couple of "stub" functions at the bottom.
123 : #ifdef LIBMESH_HAVE_VTK
124 :
125 : // Initialize the static _element_maps map.
126 : std::map<ElemMappingType, VTKIO::ElementMaps> VTKIO::_element_maps = VTKIO::build_element_maps();
127 :
128 : // Static function which constructs the ElementMaps object.
129 : std::map<ElemMappingType, VTKIO::ElementMaps> VTKIO::build_element_maps()
130 : {
131 : // Object to be filled up
132 : std::map<ElemMappingType, VTKIO::ElementMaps> all_maps;
133 : ElementMaps em; // Lagrange element maps
134 :
135 : em.associate(EDGE2, VTK_LINE);
136 : em.associate(EDGE3, VTK_QUADRATIC_EDGE);
137 : em.associate(TRI3, VTK_TRIANGLE);
138 : em.associate(TRI6, VTK_QUADRATIC_TRIANGLE);
139 : em.associate(QUAD4, VTK_QUAD);
140 : em.associate(QUAD8, VTK_QUADRATIC_QUAD);
141 : em.associate(TET4, VTK_TETRA);
142 : em.associate(TET10, VTK_QUADRATIC_TETRA);
143 : em.associate(HEX8, VTK_HEXAHEDRON);
144 : em.associate(HEX20, VTK_QUADRATIC_HEXAHEDRON);
145 : em.associate(HEX27, VTK_TRIQUADRATIC_HEXAHEDRON);
146 : em.associate(PRISM6, VTK_WEDGE);
147 : em.associate(PRISM15, VTK_QUADRATIC_WEDGE);
148 : em.associate(PRISM18, VTK_BIQUADRATIC_QUADRATIC_WEDGE);
149 : em.associate(PYRAMID5, VTK_PYRAMID);
150 :
151 : // VTK_BIQUADRATIC_QUAD has been around since VTK 5.0
152 : #if VTK_MAJOR_VERSION > 5 || (VTK_MAJOR_VERSION == 5 && VTK_MINOR_VERSION > 0)
153 : em.associate(QUAD9, VTK_BIQUADRATIC_QUAD);
154 : #endif
155 :
156 : // TRI3SUBDIVISION is for writing only
157 : em.writing_map[TRI3SUBDIVISION] = VTK_TRIANGLE;
158 : all_maps[ElemMappingType::LAGRANGE_MAP] = em;
159 :
160 :
161 : // VTK_BEZIER_* types were introduced in VTK 9.0
162 : #if VTK_VERSION_LESS_THAN(9,0,0)
163 : // Revert back to previous behavior when using an older version of VTK
164 : all_maps[ElemMappingType::RATIONAL_BERNSTEIN_MAP] = em;
165 : #else
166 : ElementMaps bem; // Rational Bernstein element maps
167 : bem.associate(EDGE2, VTK_LINE);
168 : bem.associate(EDGE3, VTK_BEZIER_CURVE);
169 : bem.associate(TRI3, VTK_TRIANGLE);
170 : bem.associate(TRI6, VTK_BEZIER_TRIANGLE);
171 : bem.associate(QUAD4, VTK_QUAD);
172 : bem.associate(QUAD8, VTK_QUADRATIC_QUAD);
173 : bem.associate(QUAD9, VTK_BEZIER_QUADRILATERAL);
174 : bem.associate(TET4, VTK_TETRA);
175 : bem.associate(TET10, VTK_QUADRATIC_TETRA);
176 : bem.associate(HEX8, VTK_HEXAHEDRON);
177 : bem.associate(HEX20, VTK_QUADRATIC_HEXAHEDRON);
178 : bem.associate(HEX27, VTK_BEZIER_HEXAHEDRON);
179 : bem.associate(PRISM6, VTK_WEDGE);
180 : bem.associate(PRISM15, VTK_QUADRATIC_WEDGE);
181 : bem.associate(PRISM18, VTK_BEZIER_WEDGE);
182 : bem.associate(PYRAMID5, VTK_PYRAMID);
183 : bem.writing_map[TRI3SUBDIVISION] = VTK_TRIANGLE;
184 : all_maps[ElemMappingType::RATIONAL_BERNSTEIN_MAP] = bem;
185 : #endif
186 :
187 : return all_maps;
188 : }
189 :
190 :
191 :
192 : namespace {
193 :
194 : // VTK describes a polyhedron by a "face stream" (a list of polygonal
195 : // faces, each a list of global point ids) rather than by the flat,
196 : // ordered node list used for every other supported element type. This
197 : // helper builds an equivalent libMesh C0Polyhedron from that face
198 : // stream, assigning its node pointers and adding to \p mesh any interior
199 : // "mid-element" node the C0Polyhedron construction requires.
200 : std::unique_ptr<Elem>
201 : add_vtk_polyhedron(vtkUnstructuredGrid & vtk_grid,
202 : MeshBase & mesh,
203 : const vtkIdType cell_id,
204 : vtkIntArray * node_id,
205 : const std::vector<dof_id_type> & vtk_node_to_libmesh)
206 : {
207 : // GetFaceStream() fills face_stream with
208 : // [ n_faces,
209 : // n_face0_pts, id0, id1, ...,
210 : // n_face1_pts, id0, id1, ..., ]
211 : // where all the point ids are global VTK point ids.
212 : vtkSmartPointer<vtkIdList> face_stream = vtkSmartPointer<vtkIdList>::New();
213 : vtk_grid.GetFaceStream(cell_id, face_stream);
214 :
215 : vtkIdType pos = 0;
216 : const vtkIdType n_faces = face_stream->GetId(pos++);
217 :
218 : libmesh_error_msg_if
219 : (n_faces < 4,
220 : "Error: VTK polyhedron cell " << cell_id << " has only " << n_faces <<
221 : " faces, but a polyhedron requires at least 4.");
222 :
223 : std::vector<std::shared_ptr<Polygon>> sides(cast_int<std::size_t>(n_faces));
224 : for (std::size_t s = 0; s != sides.size(); ++s)
225 : {
226 : const vtkIdType n_face_pts = face_stream->GetId(pos++);
227 : auto side = std::make_shared<C0Polygon>(cast_int<unsigned int>(n_face_pts));
228 :
229 : for (vtkIdType n = 0; n != n_face_pts; ++n)
230 : {
231 : const vtkIdType vtk_point_id = face_stream->GetId(pos++);
232 : const dof_id_type libmesh_node_id = node_id ?
233 : vtk_node_to_libmesh[vtk_point_id] :
234 : cast_int<dof_id_type>(vtk_point_id);
235 : side->set_node(cast_int<unsigned int>(n),
236 : mesh.node_ptr(libmesh_node_id));
237 : }
238 :
239 : sides[s] = std::move(side);
240 : }
241 :
242 : // Constructing the C0Polyhedron may create an interior "mid-element"
243 : // node for its default tetrahedralization; if so, it is our job to
244 : // add it to the mesh. The polyhedron's vertex node pointers were
245 : // already assigned above via its polygonal sides.
246 : std::unique_ptr<Node> mid_elem_node;
247 : auto elem = std::make_unique<C0Polyhedron>(sides, mid_elem_node);
248 : if (mid_elem_node)
249 : mesh.add_node(std::move(mid_elem_node));
250 :
251 : return elem;
252 : }
253 :
254 : } // anonymous namespace
255 :
256 :
257 :
258 : void VTKIO::read (const std::string & name)
259 : {
260 : // This is a serial-only process for now;
261 : // the Mesh should be read on processor 0 and
262 : // broadcast later
263 : libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().processor_id(), 0);
264 :
265 : // Keep track of what kinds of elements this file contains
266 : elems_of_dimension.clear();
267 : elems_of_dimension.resize(4, false);
268 :
269 : // VTK XML unstructured grid files come in two root-element forms: a
270 : // plain serial "UnstructuredGrid" file (conventionally .vtu), and a
271 : // "PUnstructuredGrid" descriptor (conventionally .pvtu) referencing one
272 : // or more serial piece files for parallel use. Each needs a different
273 : // reader class. Sniff the file's actual root <VTKFile type="..."/>
274 : // attribute rather than assuming one from the extension (or, worse,
275 : // always assuming one -- vtkXMLPUnstructuredGridReader silently yields
276 : // an empty grid, rather than any diagnostic, when handed the other
277 : // format).
278 : std::string file_type;
279 : {
280 : vtkSmartPointer<vtkXMLFileReadTester> tester =
281 : vtkSmartPointer<vtkXMLFileReadTester>::New();
282 : tester->SetFileName(name.c_str());
283 :
284 : libmesh_error_msg_if
285 : (!tester->TestReadFile(),
286 : "Error: " << name << " does not appear to be a VTK XML file.");
287 :
288 : const char * file_type_cstr = tester->GetFileDataType();
289 : libmesh_error_msg_if
290 : (!file_type_cstr,
291 : "Error: " << name << " does not specify a VTK XML data type.");
292 : file_type = file_type_cstr;
293 : }
294 :
295 : if (file_type == "PUnstructuredGrid")
296 : {
297 : vtkSmartPointer<vtkXMLPUnstructuredGridReader> reader =
298 : vtkSmartPointer<vtkXMLPUnstructuredGridReader>::New();
299 : reader->SetFileName(name.c_str());
300 : reader->Update();
301 : _vtk_grid = reader->GetOutput();
302 : }
303 : else if (file_type == "UnstructuredGrid")
304 : {
305 : vtkSmartPointer<vtkXMLUnstructuredGridReader> reader =
306 : vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
307 : reader->SetFileName(name.c_str());
308 : reader->Update();
309 : _vtk_grid = reader->GetOutput();
310 : }
311 : else
312 : libmesh_error_msg
313 : ("Error: " << name << " is a VTK XML file of type \"" << file_type
314 : << "\", but libMesh's VTKIO reader only supports \"UnstructuredGrid\" "
315 : "(.vtu) and \"PUnstructuredGrid\" (.pvtu) files.");
316 :
317 : // Get a reference to the mesh
318 : MeshBase & mesh = MeshInput<MeshBase>::mesh();
319 :
320 : // Clear out any pre-existing data from the Mesh
321 : mesh.clear();
322 :
323 : // Try to preserve any libMesh ids and subdomain ids we find in the
324 : // file. This will be null if there are none, e.g. if a non-libMesh
325 : // program wrote this file.
326 :
327 : vtkAbstractArray * abstract_elem_id =
328 : _vtk_grid->GetCellData()->GetAbstractArray("libmesh_elem_id");
329 : vtkAbstractArray * abstract_node_id =
330 : _vtk_grid->GetPointData()->GetAbstractArray("libmesh_node_id");
331 : vtkAbstractArray * abstract_subdomain_id =
332 : _vtk_grid->GetCellData()->GetAbstractArray("subdomain_id");
333 :
334 : // Get ids as integers. This will be null if they are another data
335 : // type, e.g. if a non-libMesh program used the names we thought
336 : // were unique for different data.
337 : vtkIntArray * elem_id = vtkIntArray::SafeDownCast(abstract_elem_id);
338 : vtkIntArray * node_id = vtkIntArray::SafeDownCast(abstract_node_id);
339 : vtkIntArray * subdomain_id = vtkIntArray::SafeDownCast(abstract_subdomain_id);
340 :
341 : if (abstract_elem_id && !elem_id)
342 : libmesh_warning("Found non-integral libmesh_elem_id array; forced to ignore it.\n"
343 : "This is technically valid but probably broken.");
344 :
345 : if (abstract_node_id && !node_id)
346 : libmesh_warning("Found non-integral libmesh_node_id array; forced to ignore it.\n"
347 : "This is technically valid but probably broken.");
348 :
349 : if (abstract_subdomain_id && !subdomain_id)
350 : libmesh_warning("Found non-integral subdomain_id array; forced to ignore it.\n"
351 : "This is technically valid but probably broken.");
352 :
353 : // Get the number of points from the _vtk_grid object
354 : const unsigned int vtk_num_points = static_cast<unsigned int>(_vtk_grid->GetNumberOfPoints());
355 :
356 : // Map from VTK indexing to libMesh id if necessary
357 : std::vector<dof_id_type> vtk_node_to_libmesh;
358 : if (node_id)
359 : vtk_node_to_libmesh.resize(vtk_num_points);
360 :
361 : // always numbered nicely so we can loop like this
362 : for (unsigned int i=0; i<vtk_num_points; ++i)
363 : {
364 : // add to the id map
365 : // and add the actual point
366 : double pnt[3];
367 : _vtk_grid->GetPoint(static_cast<vtkIdType>(i), pnt);
368 : Point xyz(pnt[0], pnt[1], pnt[2]);
369 :
370 : if (node_id)
371 : {
372 : auto id = node_id->GetValue(i);
373 :
374 : // It would nice to distinguish between "duplicate nodes
375 : // because one was ghosted in a parallel file segment" and
376 : // "duplicate nodes because there was a bug", but I'm not
377 : // sure how to do that with vtkXMLPUnstructuredGridReader
378 : if (!mesh.query_node_ptr(id))
379 : mesh.add_point(xyz, id);
380 : vtk_node_to_libmesh[i] = id;
381 : }
382 : else
383 : mesh.add_point(xyz, i);
384 : }
385 :
386 : // Get the number of cells from the _vtk_grid object
387 : const unsigned int vtk_num_cells = static_cast<unsigned int>(_vtk_grid->GetNumberOfCells());
388 :
389 : auto& element_map = libmesh_map_find(_element_maps, mesh.default_mapping_type());
390 :
391 : vtkSmartPointer<vtkGenericCell> cell = vtkSmartPointer<vtkGenericCell>::New();
392 : for (unsigned int i=0; i<vtk_num_cells; ++i)
393 : {
394 : _vtk_grid->GetCell(i, cell);
395 :
396 : std::unique_ptr<Elem> elem;
397 :
398 : // VTK polyhedra are described by a stream of polygonal faces
399 : // rather than by an ordered node list, so they are built
400 : // separately as C0Polyhedron elements with their nodes already
401 : // assigned via their sides.
402 : if (cell->GetCellType() == VTK_POLYHEDRON)
403 : elem = add_vtk_polyhedron(*_vtk_grid, mesh,
404 : static_cast<vtkIdType>(i),
405 : node_id, vtk_node_to_libmesh);
406 : else
407 : {
408 : // Get the libMesh element type corresponding to this VTK element type.
409 : ElemType libmesh_elem_type = element_map.find(cell->GetCellType());
410 : elem = Elem::build(libmesh_elem_type);
411 :
412 : // get the straightforward numbering from the VTK cells
413 : for (auto j : elem->node_index_range())
414 : {
415 : const auto vtk_point_id = cell->GetPointId(j);
416 : const dof_id_type libmesh_node_id = node_id ?
417 : vtk_node_to_libmesh[vtk_point_id] : vtk_point_id;
418 :
419 : elem->set_node(j, mesh.node_ptr(libmesh_node_id));
420 : }
421 :
422 : // then get the connectivity
423 : std::vector<dof_id_type> conn;
424 : elem->connectivity(0, VTK, conn);
425 :
426 : // then reshuffle the nodes according to the connectivity, this
427 : // two-time-assign would evade the definition of the vtk_mapping
428 : for (unsigned int j=0,
429 : n_conn = cast_int<unsigned int>(conn.size());
430 : j != n_conn; ++j)
431 : elem->set_node(j, mesh.node_ptr(conn[j]));
432 : }
433 :
434 : if (elem_id)
435 : {
436 : auto id = elem_id->GetValue(i);
437 : libmesh_error_msg_if
438 : (mesh.query_elem_ptr(id), "Duplicate element id " << id <<
439 : " found in libmesh_elem_ids");
440 : elem->set_id(id);
441 : }
442 : else
443 : elem->set_id(i);
444 :
445 : if (subdomain_id)
446 : {
447 : auto sbdid = subdomain_id->GetValue(i);
448 : elem->subdomain_id() = sbdid;
449 : }
450 :
451 : elems_of_dimension[elem->dim()] = true;
452 :
453 : mesh.add_elem(std::move(elem));
454 : } // end loop over VTK cells
455 :
456 : // Set the mesh dimension to the largest encountered for an element
457 : for (unsigned char i=0; i!=4; ++i)
458 : if (elems_of_dimension[i])
459 : mesh.set_mesh_dimension(i);
460 :
461 : #if LIBMESH_DIM < 3
462 : libmesh_error_msg_if(mesh.mesh_dimension() > LIBMESH_DIM,
463 : "Cannot open dimension "
464 : << mesh.mesh_dimension()
465 : << " mesh file when configured without "
466 : << mesh.mesh_dimension()
467 : << "D support.");
468 : #endif // LIBMESH_DIM < 3
469 : }
470 :
471 :
472 :
473 : void VTKIO::write_nodal_data (const std::string & fname,
474 : const std::vector<Number> & soln,
475 : const std::vector<std::string> & names)
476 : {
477 : // Warn that the .pvtu file extension should be used. Paraview
478 : // recognizes this, and it works in both serial and parallel. Only
479 : // warn about this once.
480 : if (!Utility::ends_with(fname, ".pvtu"))
481 : libmesh_do_once(libMesh::err << "The .pvtu extension should be used when writing VTK files in libMesh.");
482 :
483 : // If there are variable names being written, the solution vector
484 : // should not be empty, it should have been broadcast to all
485 : // processors by the MeshOutput base class, since VTK is a parallel
486 : // format. Verify this before going further.
487 : libmesh_error_msg_if(!names.empty() && soln.empty(),
488 : "Empty soln vector in VTKIO::write_nodal_data().");
489 :
490 : // Get a reference to the mesh
491 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
492 :
493 : // we only use Unstructured grids
494 : _vtk_grid = vtkSmartPointer<vtkUnstructuredGrid>::New();
495 : vtkSmartPointer<vtkXMLPUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLPUnstructuredGridWriter>::New();
496 : #ifdef LIBMESH_HAVE_MPI
497 : // Set VTK to the same communicator as libMesh
498 : vtkSmartPointer<vtkMPICommunicator> vtk_comm = vtkSmartPointer<vtkMPICommunicator>::New();
499 : MPI_Comm mpi_comm = mesh.comm().get();
500 : vtkMPICommunicatorOpaqueComm vtk_opaque_comm(&mpi_comm);
501 : vtk_comm->InitializeExternal(&vtk_opaque_comm);
502 :
503 : vtkSmartPointer<vtkMPIController> vtk_mpi_ctrl = vtkSmartPointer<vtkMPIController>::New();
504 : vtk_mpi_ctrl->SetCommunicator(vtk_comm);
505 :
506 : writer->SetController(vtk_mpi_ctrl);
507 : #endif
508 :
509 : // add nodes to the grid and update _local_node_map
510 : _local_node_map.clear();
511 : this->nodes_to_vtk();
512 :
513 : // add cells to the grid
514 : this->cells_to_vtk();
515 :
516 : // add nodal solutions to the grid, if solutions are given
517 : if (names.size() > 0)
518 : {
519 : std::size_t num_vars = names.size();
520 : std::vector<Number> local_values;
521 :
522 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
523 : std::vector<Real> local_real_values;
524 : #endif
525 :
526 : for (std::size_t variable=0; variable<num_vars; ++variable)
527 : {
528 : get_local_node_values(local_values, variable, soln, names);
529 :
530 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
531 : // write real part
532 : local_real_values.resize(local_values.size());
533 : std::transform(local_values.begin(), local_values.end(),
534 : local_real_values.begin(),
535 : [](Number x) { return x.real(); });
536 : node_values_to_vtk(names[variable] + "_real", local_real_values);
537 :
538 : // write imaginary part
539 : local_real_values.resize(local_values.size());
540 : std::transform(local_values.begin(), local_values.end(),
541 : local_real_values.begin(),
542 : [](Number x) { return x.imag(); });
543 : node_values_to_vtk(names[variable] + "_imag", local_real_values);
544 : #else
545 : node_values_to_vtk(names[variable], local_values);
546 : #endif
547 : }
548 : }
549 :
550 : // Tell the writer how many partitions exist and on which processor
551 : // we are currently
552 : writer->SetNumberOfPieces(mesh.n_processors());
553 : writer->SetStartPiece(mesh.processor_id());
554 : writer->SetEndPiece(mesh.processor_id());
555 :
556 : // partitions overlap by one node
557 : // FIXME: According to this document
558 : // http://paraview.org/Wiki/images/5/51/SC07_tut107_ParaView_Handouts.pdf
559 : // the ghosts are cells rather than nodes.
560 : writer->SetGhostLevel(1);
561 :
562 : // VTK 6 replaces SetInput() with SetInputData(). See
563 : // http://www.vtk.org/Wiki/VTK/VTK_6_Migration/Replacement_of_SetInput
564 : // for the full explanation.
565 : #if VTK_VERSION_LESS_THAN(6,0,0)
566 : writer->SetInput(_vtk_grid);
567 : #else
568 : writer->SetInputData(_vtk_grid);
569 : #endif
570 :
571 : writer->SetFileName(fname.c_str());
572 : writer->SetDataModeToAscii();
573 :
574 : // compress the output, if desired (switches also to binary)
575 : if (this->_compress)
576 : {
577 : #if !VTK_VERSION_LESS_THAN(5,6,0)
578 : writer->SetCompressorTypeToZLib();
579 : writer->SetDataModeToBinary();
580 : #else
581 : libmesh_do_once(libMesh::err << "Compression not implemented with old VTK libs!" << std::endl;);
582 : #endif
583 : }
584 :
585 : writer->Write();
586 :
587 : }
588 :
589 :
590 :
591 : vtkUnstructuredGrid * VTKIO::get_vtk_grid()
592 : {
593 : return _vtk_grid;
594 : }
595 :
596 :
597 :
598 : void VTKIO::set_compression(bool b)
599 : {
600 : this->_compress = b;
601 : }
602 :
603 :
604 :
605 : void VTKIO::nodes_to_vtk()
606 : {
607 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
608 :
609 : // containers for points and coordinates of points
610 : vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
611 : vtkSmartPointer<vtkDoubleArray> pcoords = vtkSmartPointer<vtkDoubleArray>::New();
612 : // if this grid is to be used in VTK then the dimension of the points should be 3
613 : pcoords->SetNumberOfComponents(LIBMESH_DIM);
614 : pcoords->Allocate(3*mesh.n_local_nodes());
615 : points->SetNumberOfPoints(mesh.n_local_nodes()); // it seems that it needs this to prevent a segfault
616 :
617 : // SetRationalWeights() was introduced in VTK 9.0
618 : #if !VTK_VERSION_LESS_THAN(9,0,0)
619 : bool have_weights = false;
620 : int weight_index = 0;
621 : vtkSmartPointer<vtkDoubleArray> rational_weights;
622 :
623 : if (mesh.default_mapping_type() == ElemMappingType::RATIONAL_BERNSTEIN_MAP)
624 : {
625 : rational_weights = vtkSmartPointer<vtkDoubleArray>::New();
626 : rational_weights->SetName("RationalWeights");
627 : rational_weights->SetNumberOfComponents(1);
628 : weight_index = static_cast<int>(mesh.default_mapping_data());
629 : have_weights = true;
630 : }
631 : #endif
632 :
633 : vtkSmartPointer<vtkIntArray> node_id = vtkSmartPointer<vtkIntArray>::New();
634 : node_id->SetName("libmesh_node_id");
635 : node_id->SetNumberOfComponents(1);
636 :
637 : unsigned int local_node_counter = 0;
638 :
639 : for (const auto & node_ptr : mesh.local_node_ptr_range())
640 : {
641 : const Node & node = *node_ptr;
642 :
643 : double pnt[3] = {0, 0, 0};
644 : for (unsigned int i=0; i<LIBMESH_DIM; ++i)
645 : pnt[i] = double(node(i));
646 :
647 : // Fill mapping between global and local node numbers
648 : _local_node_map[node.id()] = local_node_counter;
649 :
650 : // add point
651 : #if VTK_VERSION_LESS_THAN(7,1,0)
652 : pcoords->InsertNextTupleValue(pnt);
653 : #else
654 : pcoords->InsertNextTuple(pnt);
655 : #endif
656 : #if !VTK_VERSION_LESS_THAN(9,0,0)
657 : if (have_weights)
658 : {
659 : Real weight = node.get_extra_datum<Real>(weight_index);
660 : rational_weights->InsertTuple1(local_node_counter, double(weight));
661 : }
662 : #endif
663 :
664 : node_id->InsertTuple1(local_node_counter, node.id());
665 :
666 : ++local_node_counter;
667 : }
668 :
669 : // add coordinates to points
670 : points->SetData(pcoords);
671 :
672 : // add points to grid
673 : _vtk_grid->SetPoints(points);
674 : _vtk_grid->GetPointData()->AddArray(node_id);
675 :
676 : #if !VTK_VERSION_LESS_THAN(9,0,0)
677 : if (have_weights)
678 : _vtk_grid->GetPointData()->SetRationalWeights(rational_weights);
679 :
680 : #endif
681 : }
682 :
683 :
684 :
685 : void VTKIO::cells_to_vtk()
686 : {
687 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
688 :
689 : auto& element_map = libmesh_map_find(_element_maps, mesh.default_mapping_type());
690 : vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
691 : vtkSmartPointer<vtkIdList> pts = vtkSmartPointer<vtkIdList>::New();
692 :
693 : std::vector<int> types(mesh.n_active_local_elem());
694 :
695 : // We already created this but we need to add more if we have any
696 : // ghost nodes
697 : vtkAbstractArray * abstract_node_id =
698 : _vtk_grid->GetPointData()->GetAbstractArray("libmesh_node_id");
699 : vtkIntArray * node_id = vtkIntArray::SafeDownCast(abstract_node_id);
700 : libmesh_assert(node_id);
701 :
702 : vtkSmartPointer<vtkIntArray> elem_id = vtkSmartPointer<vtkIntArray>::New();
703 : elem_id->SetName("libmesh_elem_id");
704 : elem_id->SetNumberOfComponents(1);
705 :
706 : vtkSmartPointer<vtkIntArray> subdomain_id = vtkSmartPointer<vtkIntArray>::New();
707 : subdomain_id->SetName("subdomain_id");
708 : subdomain_id->SetNumberOfComponents(1);
709 :
710 : vtkSmartPointer<vtkIntArray> elem_proc_id = vtkSmartPointer<vtkIntArray>::New();
711 : elem_proc_id->SetName("processor_id");
712 : elem_proc_id->SetNumberOfComponents(1);
713 :
714 : unsigned active_element_counter = 0;
715 : for (const auto & elem : mesh.active_local_element_ptr_range())
716 : {
717 : // When using rational bernstein these hold the weights
718 : if ( elem->type() == NODEELEM )
719 : continue;
720 :
721 : pts->SetNumberOfIds(elem->n_nodes());
722 :
723 : // get the connectivity for this element
724 : std::vector<dof_id_type> conn;
725 : elem->connectivity(0, VTK, conn);
726 :
727 : for (unsigned int i=0,
728 : n_conn = cast_int<unsigned int>(conn.size());
729 : i != n_conn; ++i)
730 : {
731 : // If the node ID is not found in the _local_node_map, we'll
732 : // add it to the _vtk_grid. NOTE[JWP]: none of the examples
733 : // I have actually enters this section of code...
734 : if (!_local_node_map.count(conn[i]))
735 : {
736 : dof_id_type global_node_id = elem->node_id(i);
737 :
738 : const Point & the_node = mesh.point(global_node_id);
739 :
740 : // InsertNextPoint accepts either a double or float array of length 3.
741 : double pt[3] = {0., 0., 0.};
742 : for (unsigned int d=0; d<LIBMESH_DIM; ++d)
743 : pt[d] = double(the_node(d));
744 :
745 : // Insert the point into the _vtk_grid
746 : vtkIdType local = _vtk_grid->GetPoints()->InsertNextPoint(pt);
747 :
748 : // Update the _local_node_map with the ID returned by VTK
749 : _local_node_map[global_node_id] =
750 : cast_int<dof_id_type>(local);
751 :
752 : node_id->InsertTuple1(local, global_node_id);
753 : }
754 :
755 : // Otherwise, the node ID was found in the _local_node_map, so
756 : // insert it into the vtkIdList.
757 : pts->InsertId(i, _local_node_map[conn[i]]);
758 : }
759 :
760 : vtkIdType vtkcellid = cells->InsertNextCell(pts);
761 : types[active_element_counter] = cast_int<int>(element_map.find(elem->type()));
762 :
763 : elem_id->InsertTuple1(vtkcellid, elem->id());
764 : subdomain_id->InsertTuple1(vtkcellid, elem->subdomain_id());
765 : elem_proc_id->InsertTuple1(vtkcellid, elem->processor_id());
766 : ++active_element_counter;
767 : } // end loop over active elements
768 :
769 : _vtk_grid->SetCells(types.data(), cells);
770 : _vtk_grid->GetCellData()->AddArray(elem_id);
771 : _vtk_grid->GetCellData()->AddArray(subdomain_id);
772 : _vtk_grid->GetCellData()->AddArray(elem_proc_id);
773 : }
774 :
775 : void VTKIO::node_values_to_vtk(const std::string & name,
776 : const std::vector<Real> & local_values)
777 : {
778 : vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
779 : data->SetName(name.c_str());
780 :
781 : libmesh_assert_equal_to(_local_node_map.size(), local_values.size());
782 :
783 : // number of local and ghost nodes
784 : data->SetNumberOfValues(_local_node_map.size());
785 :
786 : // copy values into vtk
787 : for (auto i : index_range(local_values)) {
788 : data->SetValue(i, double(local_values[i]));
789 : }
790 :
791 : _vtk_grid->GetPointData()->AddArray(data);
792 : }
793 :
794 : void VTKIO::get_local_node_values(std::vector<Number> & local_values,
795 : std::size_t variable,
796 : const std::vector<Number> & soln,
797 : const std::vector<std::string> & names)
798 : {
799 : const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
800 : std::size_t num_vars = names.size();
801 : dof_id_type num_nodes = mesh.n_nodes();
802 :
803 : local_values.clear();
804 : local_values.resize(_local_node_map.size(), 0.0);
805 :
806 : // loop over all nodes and get the solution for the current
807 : // variable, if the node is in the current partition
808 : for (dof_id_type k=0; k<num_nodes; ++k)
809 : if (const auto local_node_it = _local_node_map.find(k);
810 : local_node_it != _local_node_map.end())
811 : local_values[local_node_it->second] = soln[k*num_vars + variable];
812 : }
813 :
814 :
815 :
816 : /**
817 : * FIXME: This is known to write nonsense on AMR meshes
818 : * and it strips the imaginary parts of complex Numbers
819 : *
820 : * This function is not currently used by anything, so it is commented
821 : * out, and may eventually be removed entirely.
822 : */
823 : // void VTKIO::system_vectors_to_vtk(const EquationSystems & es,
824 : // vtkUnstructuredGrid *& grid)
825 : // {
826 : // if (MeshOutput<MeshBase>::mesh().processor_id() == 0)
827 : // {
828 : // std::map<std::string, std::vector<Number>> vecs;
829 : // for (unsigned int i=0; i<es.n_systems(); ++i)
830 : // {
831 : // const System & sys = es.get_system(i);
832 : // System::const_vectors_iterator v_end = sys.vectors_end();
833 : // System::const_vectors_iterator it = sys.vectors_begin();
834 : // for (; it!= v_end; ++it)
835 : // {
836 : // // for all vectors on this system
837 : // std::vector<Number> values;
838 : // // libMesh::out<<"it "<<it->first<<std::endl;
839 : //
840 : // it->second->localize_to_one(values, 0);
841 : // // libMesh::out<<"finish localize"<<std::endl;
842 : // vecs[it->first] = values;
843 : // }
844 : // }
845 : //
846 : // std::map<std::string, std::vector<Number>>::iterator it = vecs.begin();
847 : //
848 : // for (; it!=vecs.end(); ++it)
849 : // {
850 : // vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
851 : // data->SetName(it->first.c_str());
852 : // libmesh_assert_equal_to (it->second.size(), es.get_mesh().n_nodes());
853 : // data->SetNumberOfValues(it->second.size());
854 : //
855 : // for (auto i : index_range(it->second))
856 : // {
857 : // #ifdef LIBMESH_USE_COMPLEX_NUMBERS
858 : // libmesh_do_once (libMesh::err << "Only writing the real part for complex numbers!\n"
859 : // << "if you need this support contact " << LIBMESH_PACKAGE_BUGREPORT
860 : // << std::endl);
861 : // data->SetValue(i, it->second[i].real());
862 : // #else
863 : // data->SetValue(i, it->second[i]);
864 : // #endif
865 : //
866 : // }
867 : // grid->GetPointData()->AddArray(data);
868 : // }
869 : // }
870 : // }
871 :
872 :
873 :
874 : #else // !LIBMESH_HAVE_VTK
875 :
876 0 : void VTKIO::read (const std::string & name)
877 : {
878 0 : libmesh_error_msg("Cannot read VTK file: " << name \
879 : << "\nYou must have VTK installed and correctly configured to read VTK meshes.");
880 : }
881 :
882 :
883 :
884 0 : void VTKIO::write_nodal_data (const std::string & fname,
885 : const std::vector<Number> &,
886 : const std::vector<std::string> &)
887 : {
888 0 : libmesh_error_msg("Cannot write VTK file: " << fname \
889 : << "\nYou must have VTK installed and correctly configured to read VTK meshes.");
890 : }
891 :
892 :
893 : #endif // LIBMESH_HAVE_VTK
894 :
895 :
896 :
897 : } // namespace libMesh
|
