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