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