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" 30 #ifdef LIBMESH_HAVE_VTK 34 #include "libmesh/ignore_warnings.h" 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" 50 #ifdef LIBMESH_HAVE_MPI 52 #include "vtkMPICommunicator.h" 53 #include "vtkMPIController.h" 56 #include "libmesh/restore_warnings.h" 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) 74 #endif // LIBMESH_HAVE_VTK 85 #ifdef LIBMESH_HAVE_VTK
96 #ifdef LIBMESH_HAVE_VTK
107 std::vector<Number> soln;
108 std::vector<std::string> names;
116 #ifdef LIBMESH_HAVE_VTK 125 std::map<ElemMappingType, VTKIO::ElementMaps> all_maps;
145 #if VTK_MAJOR_VERSION > 5 || (VTK_MAJOR_VERSION == 5 && VTK_MINOR_VERSION > 0) 155 #if VTK_VERSION_LESS_THAN(9,0,0) 197 typedef vtkSmartPointer<vtkXMLPUnstructuredGridReader> MyReader;
198 MyReader reader = MyReader::New();
201 reader->SetFileName(
name.c_str());
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");
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);
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.");
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.");
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.");
246 const unsigned int vtk_num_points =
static_cast<unsigned int>(
_vtk_grid->GetNumberOfPoints());
249 std::vector<dof_id_type> vtk_node_to_libmesh;
251 vtk_node_to_libmesh.resize(vtk_num_points);
254 for (
unsigned int i=0; i<vtk_num_points; ++i)
259 _vtk_grid->GetPoint(static_cast<vtkIdType>(i), pnt);
260 Point xyz(pnt[0], pnt[1], pnt[2]);
264 auto id = node_id->GetValue(i);
272 vtk_node_to_libmesh[i] = id;
279 const unsigned int vtk_num_cells =
static_cast<unsigned int>(
_vtk_grid->GetNumberOfCells());
283 vtkSmartPointer<vtkGenericCell> cell = vtkSmartPointer<vtkGenericCell>::New();
284 for (
unsigned int i=0; i<vtk_num_cells; ++i)
289 ElemType libmesh_elem_type = element_map.find(cell->GetCellType());
293 for (
auto j : elem->node_index_range())
295 const auto vtk_point_id = cell->GetPointId(j);
297 vtk_node_to_libmesh[vtk_point_id] : vtk_point_id;
303 std::vector<dof_id_type> conn;
304 elem->connectivity(0,
VTK, conn);
308 for (
unsigned int j=0,
309 n_conn = cast_int<unsigned int>(conn.size());
315 auto id = elem_id->GetValue(i);
318 " found in libmesh_elem_ids");
326 auto sbdid = subdomain_id->GetValue(i);
327 elem->subdomain_id() = sbdid;
336 for (
unsigned char i=0; i!=4; ++i)
342 "Cannot open dimension " 344 <<
" mesh file when configured without " 347 #endif // LIBMESH_DIM < 3 353 const std::vector<Number> & soln,
354 const std::vector<std::string> & names)
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.");
366 libmesh_error_msg_if(!names.empty() && soln.empty(),
367 "Empty soln vector in VTKIO::write_nodal_data().");
373 _vtk_grid = vtkSmartPointer<vtkUnstructuredGrid>::New();
374 vtkSmartPointer<vtkXMLPUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLPUnstructuredGridWriter>::New();
375 #ifdef LIBMESH_HAVE_MPI 377 vtkSmartPointer<vtkMPICommunicator> vtk_comm = vtkSmartPointer<vtkMPICommunicator>::New();
379 vtkMPICommunicatorOpaqueComm vtk_opaque_comm(&mpi_comm);
380 vtk_comm->InitializeExternal(&vtk_opaque_comm);
382 vtkSmartPointer<vtkMPIController> vtk_mpi_ctrl = vtkSmartPointer<vtkMPIController>::New();
383 vtk_mpi_ctrl->SetCommunicator(vtk_comm);
385 writer->SetController(vtk_mpi_ctrl);
396 if (names.size() > 0)
398 std::size_t num_vars = names.size();
399 std::vector<Number> local_values;
401 #ifdef LIBMESH_USE_COMPLEX_NUMBERS 402 std::vector<Real> local_real_values;
405 for (std::size_t variable=0; variable<num_vars; ++variable)
409 #ifdef LIBMESH_USE_COMPLEX_NUMBERS 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(); });
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(); });
439 writer->SetGhostLevel(1);
444 #if VTK_VERSION_LESS_THAN(6,0,0) 450 writer->SetFileName(fname.c_str());
451 writer->SetDataModeToAscii();
456 #if !VTK_VERSION_LESS_THAN(5,6,0) 457 writer->SetCompressorTypeToZLib();
458 writer->SetDataModeToBinary();
460 libmesh_do_once(
libMesh::err <<
"Compression not implemented with old VTK libs!" << std::endl;);
489 vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
490 vtkSmartPointer<vtkDoubleArray> pcoords = vtkSmartPointer<vtkDoubleArray>::New();
492 pcoords->SetNumberOfComponents(LIBMESH_DIM);
497 #if !VTK_VERSION_LESS_THAN(9,0,0) 498 bool have_weights =
false;
499 int weight_index = 0;
500 vtkSmartPointer<vtkDoubleArray> rational_weights;
504 rational_weights = vtkSmartPointer<vtkDoubleArray>::New();
505 rational_weights->SetName(
"RationalWeights");
506 rational_weights->SetNumberOfComponents(1);
512 vtkSmartPointer<vtkIntArray> node_id = vtkSmartPointer<vtkIntArray>::New();
513 node_id->SetName(
"libmesh_node_id");
514 node_id->SetNumberOfComponents(1);
516 unsigned int local_node_counter = 0;
518 for (
const auto & node_ptr :
mesh.local_node_ptr_range())
520 const Node & node = *node_ptr;
522 double pnt[3] = {0, 0, 0};
523 for (
unsigned int i=0; i<LIBMESH_DIM; ++i)
524 pnt[i] =
double(node(i));
530 #if VTK_VERSION_LESS_THAN(7,1,0) 531 pcoords->InsertNextTupleValue(pnt);
533 pcoords->InsertNextTuple(pnt);
535 #if !VTK_VERSION_LESS_THAN(9,0,0) 539 rational_weights->InsertTuple1(local_node_counter,
double(
weight));
543 node_id->InsertTuple1(local_node_counter, node.
id());
545 ++local_node_counter;
549 points->SetData(pcoords);
553 _vtk_grid->GetPointData()->AddArray(node_id);
555 #if !VTK_VERSION_LESS_THAN(9,0,0) 557 _vtk_grid->GetPointData()->SetRationalWeights(rational_weights);
569 vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
570 vtkSmartPointer<vtkIdList> pts = vtkSmartPointer<vtkIdList>::New();
576 vtkAbstractArray * abstract_node_id =
577 _vtk_grid->GetPointData()->GetAbstractArray(
"libmesh_node_id");
578 vtkIntArray * node_id = vtkIntArray::SafeDownCast(abstract_node_id);
581 vtkSmartPointer<vtkIntArray> elem_id = vtkSmartPointer<vtkIntArray>::New();
582 elem_id->SetName(
"libmesh_elem_id");
583 elem_id->SetNumberOfComponents(1);
585 vtkSmartPointer<vtkIntArray> subdomain_id = vtkSmartPointer<vtkIntArray>::New();
586 subdomain_id->SetName(
"subdomain_id");
587 subdomain_id->SetNumberOfComponents(1);
589 vtkSmartPointer<vtkIntArray> elem_proc_id = vtkSmartPointer<vtkIntArray>::New();
590 elem_proc_id->SetName(
"processor_id");
591 elem_proc_id->SetNumberOfComponents(1);
593 unsigned active_element_counter = 0;
594 for (
const auto & elem :
mesh.active_local_element_ptr_range())
600 pts->SetNumberOfIds(elem->n_nodes());
603 std::vector<dof_id_type> conn;
604 elem->connectivity(0,
VTK, conn);
606 for (
unsigned int i=0,
607 n_conn = cast_int<unsigned int>(conn.size());
620 double pt[3] = {0., 0., 0.};
621 for (
unsigned int d=0; d<LIBMESH_DIM; ++d)
622 pt[d] =
double(the_node(d));
625 vtkIdType local =
_vtk_grid->GetPoints()->InsertNextPoint(pt);
629 cast_int<dof_id_type>(local);
631 node_id->InsertTuple1(local, global_node_id);
639 vtkIdType vtkcellid = cells->InsertNextCell(pts);
640 types[active_element_counter] = cast_int<int>(element_map.find(elem->type()));
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;
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);
655 const std::vector<Real> & local_values)
657 vtkSmartPointer<vtkDoubleArray> data = vtkSmartPointer<vtkDoubleArray>::New();
658 data->SetName(
name.c_str());
667 data->SetValue(i,
double(local_values[i]));
670 _vtk_grid->GetPointData()->AddArray(data);
674 std::size_t variable,
675 const std::vector<Number> & soln,
676 const std::vector<std::string> & names)
679 std::size_t num_vars = names.size();
682 local_values.clear();
690 local_values[local_node_it->second] = soln[k*num_vars + variable];
753 #else // !LIBMESH_HAVE_VTK 757 libmesh_error_msg(
"Cannot read VTK file: " <<
name \
758 <<
"\nYou must have VTK installed and correctly configured to read VTK meshes.");
764 const std::vector<Number> &,
765 const std::vector<std::string> &)
767 libmesh_error_msg(
"Cannot write VTK file: " << fname \
768 <<
"\nYou must have VTK installed and correctly configured to read VTK meshes.");
772 #endif // LIBMESH_HAVE_VTK std::string name(const ElemQuality q)
This function returns a string containing some name for q.
void set_compression(bool b)
Setter for compression flag.
void node_values_to_vtk(const std::string &name, const std::vector< Real > &local_values)
write the nodal values of soln to a vtkUnstructuredGrid
static std::map< ElemMappingType, ElementMaps > _element_maps
ElementMaps objects that are built statically and used by all instances of this class.
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
A Node is like a Point, but with more information.
vtkSmartPointer< vtkUnstructuredGrid > _vtk_grid
Write the system vectors to vtk.
bool _compress
Flag to indicate whether the output should be compressed.
virtual void write(const std::string &) override
Output the mesh without solutions to a .pvtu file.
dof_id_type n_local_nodes() const
const Parallel::Communicator & comm() const
This class defines an abstract interface for Mesh output.
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.
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.
ElemMappingType default_mapping_type() const
Returns the default master space to physical space mapping basis functions to be used on newly added ...
virtual void read(const std::string &) override
This method implements reading a mesh from a specified file in VTK format.
static std::map< ElemMappingType, ElementMaps > build_element_maps()
Static function used to construct _element_maps.
void nodes_to_vtk()
write the nodes from the mesh into a vtkUnstructuredGrid and update the local_node_map.
dof_id_type n_active_local_elem() const
processor_id_type n_processors() const
unsigned char default_mapping_data() const
Returns any default data value used by the master space to physical space mapping.
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)
virtual const Node * query_node_ptr(const dof_id_type i) const =0
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...
vtkUnstructuredGrid * get_vtk_grid()
Get a pointer to the VTK unstructured grid data structure.
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
std::map< ElemType, vtkIdType > writing_map
virtual void clear()
Deletes all the element and node data that is currently stored.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
VTKIO(MeshBase &mesh)
Constructor.
void associate(ElemType libmesh_type, vtkIdType vtk_type)
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.
void cells_to_vtk()
write the cells from the mesh into a vtkUnstructuredGrid
unsigned int mesh_dimension() const
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...
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.
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
virtual dof_id_type n_nodes() const =0