141 const double relative_tolerance)
144 std::ifstream file(vtk_filename);
146 ExcMessage(
"VTK file not found: " + vtk_filename));
148 vtkSmartPointer<vtkDataObject> data_object;
149 const auto dot_pos = vtk_filename.find_last_of(
'.');
150 const std::string ext =
151 (dot_pos == std::string::npos ?
"" : vtk_filename.substr(dot_pos + 1));
153# if DEAL_II_VTK_VERSION_NUMBER >= VTK_VERSION_CHECK(9, 3, 0)
156 auto reader = vtkSmartPointer<vtkXMLGenericDataObjectReader>::New();
157 reader->SetFileName(vtk_filename.c_str());
159 data_object = reader->GetOutput();
161 else if (ext ==
"vtk")
165 auto reader = vtkSmartPointer<vtkGenericDataObjectReader>::New();
166 reader->SetFileName(vtk_filename.c_str());
168 data_object = reader->GetOutput();
172 ExcMessage(
"Unsupported file extension '" + ext +
173 "'. Use '.vtu' for VTK XML format or '.vtk' "
174 "for legacy VTK format."));
177 ExcMessage(
"Unsupported file extension '" + ext +
178 "'. With VTK < 9.3 only '.vtk' legacy files can "
181 auto reader = vtkSmartPointer<vtkDataSetReader>::New();
182 reader->SetFileName(vtk_filename.c_str());
184 data_object = reader->GetOutput();
190# if DEAL_II_VTK_VERSION_NUMBER >= VTK_VERSION_CHECK(9, 3, 0)
193 auto cleaner = vtkSmartPointer<vtkCleanUnstructuredGrid>::New();
194 cleaner->SetToleranceIsAbsolute(
false);
195 cleaner->SetTolerance(relative_tolerance);
196 cleaner->SetInputData(out);
198 out->ShallowCopy(cleaner->GetOutput());
202 (void)relative_tolerance;
203 deallog <<
"VTK version < 9.3: skipping cleanup step." << std::endl;
257 const vtkUnstructuredGrid &unstructured_grid,
259 const std::string &material_id_field,
260 const std::string &boundary_id_field,
261 const std::string &manifold_id_field)
263 auto &grid =
const_cast<vtkUnstructuredGrid &
>(unstructured_grid);
265 auto get_cell_data_array =
266 [&grid](
const std::string &name) -> vtkDataArray * {
270 vtkCellData *cell_data = grid.GetCellData();
271 if (cell_data ==
nullptr)
274 vtkDataArray *data_array = cell_data->GetArray(name.c_str());
275 if (data_array !=
nullptr)
276 AssertThrow(data_array->GetNumberOfComponents() == 1,
277 ExcMessage(
"The VTK cell data array '" + name +
278 "' must be scalar."));
285 vtkDataArray *material_ids = get_cell_data_array(material_id_field);
286 vtkDataArray *boundary_ids = get_cell_data_array(boundary_id_field);
287 vtkDataArray *manifold_ids = get_cell_data_array(manifold_id_field);
289 auto material_id_from_vtk = [](vtkDataArray &material_ids,
290 const vtkIdType vtk_id) {
291 const double material_value = material_ids.GetComponent(vtk_id, 0);
292 if (material_value < 0)
299 vtkPoints *vtk_points = grid.GetPoints();
300 const vtkIdType n_points = vtk_points->GetNumberOfPoints();
301 std::vector<Point<spacedim>> points(n_points);
302 for (vtkIdType i = 0; i < n_points; ++i)
304 std::array<double, 3> coords = {{0, 0, 0}};
305 vtk_points->GetPoint(i, coords.data());
306 for (
unsigned int d = 0; d < spacedim; ++d)
307 points[i][d] = coords[d];
308 for (
unsigned int d = spacedim; d < 3; ++d)
312 "VTK grid has non-zero coordinate in unused dimension."));
316 std::vector<CellData<dim>> cells;
318 const vtkIdType n_cells = grid.GetNumberOfCells();
319 std::vector<bool> boundary_vertex_ids_present(n_points,
false);
320 std::vector<bool> manifold_vertex_ids_present(n_points,
false);
321 std::vector<types::boundary_id> boundary_vertex_ids(
323 std::vector<types::manifold_id> manifold_vertex_ids(
325 for (vtkIdType vtk_id = 0; vtk_id < n_cells; ++vtk_id)
327 vtkCell *cell = grid.GetCell(vtk_id);
328 if (cell->GetCellDimension() < dim)
330 if (boundary_ids ==
nullptr && manifold_ids ==
nullptr)
333 auto set_subcell_ids = [vtk_id, boundary_ids, manifold_ids](
335 if (boundary_ids !=
nullptr)
337 double bval = boundary_ids->GetComponent(vtk_id, 0);
341 cell_data.boundary_id =
344 else if (manifold_ids !=
nullptr)
347 if (manifold_ids !=
nullptr)
349 double mval = manifold_ids->GetComponent(vtk_id, 0);
353 cell_data.manifold_id =
358 if constexpr (dim == 1)
360 if (cell->GetCellType() == VTK_VERTEX)
363 cell->GetNumberOfPoints() == 1,
365 "Only vertex subcells with 1 point are supported."));
367 const vtkIdType vertex_index = cell->GetPointId(0);
370 if (boundary_ids !=
nullptr)
372 const double boundary_value =
373 boundary_ids->GetComponent(vtk_id, 0);
374 if (boundary_value >= 0)
376 boundary_vertex_ids_present[vertex_index] =
true;
377 boundary_vertex_ids[vertex_index] =
382 if (manifold_ids !=
nullptr)
384 const double manifold_value =
385 manifold_ids->GetComponent(vtk_id, 0);
386 manifold_vertex_ids_present[vertex_index] =
true;
387 manifold_vertex_ids[vertex_index] =
388 (manifold_value < 0 ?
394 else if constexpr (dim == 2)
396 if (cell->GetCellType() == VTK_LINE)
399 cell->GetNumberOfPoints() == 2,
401 "Only line subcells with 2 points are supported."));
404 for (
unsigned int j = 0; j < 2; ++j)
405 cell_data.
vertices[j] = cell->GetPointId(j);
407 set_subcell_ids(cell_data);
411 else if constexpr (dim == 3)
413 if (cell->GetCellType() == VTK_LINE)
416 cell->GetNumberOfPoints() == 2,
418 "Only line subcells with 2 points are supported."));
421 for (
unsigned int j = 0; j < 2; ++j)
422 cell_data.
vertices[j] = cell->GetPointId(j);
424 set_subcell_ids(cell_data);
427 else if (cell->GetCellType() == VTK_TRIANGLE ||
428 cell->GetCellType() == VTK_QUAD)
431 cell->GetNumberOfPoints() == 3 ||
432 cell->GetNumberOfPoints() == 4,
434 "Only triangle and quad face subcells are supported."));
437 for (
unsigned int j = 0; j < cell->GetNumberOfPoints(); ++j)
438 cell_data.
vertices[j] = cell->GetPointId(j);
444 if (cell->GetCellType() == VTK_QUAD)
447 set_subcell_ids(cell_data);
456 ExcMessage(
"Unsupported VTK cell dimension."));
458 if constexpr (dim == 1)
460 if (cell->GetCellType() != VTKCellType::VTK_LINE)
463 "Unsupported cell type in 1D VTK file: only "
464 "VTK_LINE is supported."));
467 "Only line cells with 2 points are supported."));
469 for (
unsigned int j = 0; j < 2; ++j)
470 cell_data.
vertices[j] = cell->GetPointId(j);
472 if (material_ids !=
nullptr)
474 material_id_from_vtk(*material_ids, vtk_id);
475 if (boundary_ids !=
nullptr)
477 const double boundary_value =
478 boundary_ids->GetComponent(vtk_id, 0);
479 if (boundary_value >= 0)
483 if (manifold_ids !=
nullptr)
485 double mval = manifold_ids->GetComponent(vtk_id, 0);
491 cells.push_back(cell_data);
493 else if constexpr (dim == 2)
495 if (cell->GetCellType() == VTKCellType::VTK_QUAD)
499 "Only quad cells with 4 points are supported."));
501 for (
unsigned int j = 0; j < 4; ++j)
502 cell_data.
vertices[j] = cell->GetPointId(j);
505 if (material_ids !=
nullptr)
507 material_id_from_vtk(*material_ids, vtk_id);
508 if (manifold_ids !=
nullptr)
510 double mval = manifold_ids->GetComponent(vtk_id, 0);
517 cells.push_back(cell_data);
519 else if (cell->GetCellType() == VTKCellType::VTK_TRIANGLE)
522 cell->GetNumberOfPoints() == 3,
524 "Only triangle cells with 3 points are supported."));
526 for (
unsigned int j = 0; j < 3; ++j)
527 cell_data.
vertices[j] = cell->GetPointId(j);
529 if (material_ids !=
nullptr)
531 material_id_from_vtk(*material_ids, vtk_id);
532 if (manifold_ids !=
nullptr)
534 double mval = manifold_ids->GetComponent(vtk_id, 0);
541 cells.push_back(cell_data);
546 "Unsupported cell type in 2D VTK file: only "
547 "VTK_QUAD and VTK_TRIANGLE are supported."));
549 else if constexpr (dim == 3)
551 if (cell->GetCellType() == VTKCellType::VTK_HEXAHEDRON)
555 "Only hex cells with 8 points are supported."));
557 for (
unsigned int j = 0; j < 8; ++j)
558 cell_data.
vertices[j] = cell->GetPointId(j);
564 if (material_ids !=
nullptr)
566 material_id_from_vtk(*material_ids, vtk_id);
567 if (manifold_ids !=
nullptr)
569 double mval = manifold_ids->GetComponent(vtk_id, 0);
576 cells.push_back(cell_data);
578 else if (cell->GetCellType() == VTKCellType::VTK_TETRA)
581 cell->GetNumberOfPoints() == 4,
583 "Only tetrahedron cells with 4 points are supported."));
585 for (
unsigned int j = 0; j < 4; ++j)
586 cell_data.
vertices[j] = cell->GetPointId(j);
588 if (material_ids !=
nullptr)
590 material_id_from_vtk(*material_ids, vtk_id);
591 if (manifold_ids !=
nullptr)
593 double mval = manifold_ids->GetComponent(vtk_id, 0);
600 cells.push_back(cell_data);
602 else if (cell->GetCellType() == VTKCellType::VTK_WEDGE)
606 "Only prism cells with 6 points are supported."));
608 for (
unsigned int j = 0; j < 6; ++j)
609 cell_data.
vertices[j] = cell->GetPointId(j);
611 if (material_ids !=
nullptr)
613 material_id_from_vtk(*material_ids, vtk_id);
614 if (manifold_ids !=
nullptr)
616 double mval = manifold_ids->GetComponent(vtk_id, 0);
623 cells.push_back(cell_data);
625 else if (cell->GetCellType() == VTKCellType::VTK_PYRAMID)
628 cell->GetNumberOfPoints() == 5,
630 "Only pyramid cells with 5 points are supported."));
632 for (
unsigned int j = 0; j < 5; ++j)
633 cell_data.
vertices[j] = cell->GetPointId(j);
635 if (material_ids !=
nullptr)
637 material_id_from_vtk(*material_ids, vtk_id);
638 if (manifold_ids !=
nullptr)
640 double mval = manifold_ids->GetComponent(vtk_id, 0);
647 cells.push_back(cell_data);
653 "Unsupported cell type in 3D VTK file: only "
654 "VTK_HEXAHEDRON, VTK_TETRA, VTK_WEDGE, and VTK_PYRAMID are supported."));
665 if constexpr (dim == 1)
667 for (
unsigned int f = 0; f < cell->n_faces(); ++f)
668 if (cell->face(f)->at_boundary())
670 const unsigned int vertex_index = cell->face(f)->vertex_index(0);
672 if (boundary_vertex_ids_present[vertex_index])
673 cell->face(f)->set_boundary_id(
674 boundary_vertex_ids[vertex_index]);
675 if (manifold_vertex_ids_present[vertex_index])
676 cell->face(f)->set_manifold_id(
677 manifold_vertex_ids[vertex_index]);
687 const std::string &material_id_field,
688 const std::string &boundary_id_field,
689 const std::string &manifold_id_field)
691 auto grid = vtkSmartPointer<vtkUnstructuredGrid>::New();
692 auto points = vtkSmartPointer<vtkPoints>::New();
695 for (
unsigned int i = 0; i < tria.
n_vertices(); ++i)
697 std::array<double, 3> coords = {{0.0, 0.0, 0.0}};
698 for (
unsigned int d = 0; d < spacedim; ++d)
700 points->SetPoint(i, coords.data());
702 grid->SetPoints(points);
707 vtkSmartPointer<vtkIntArray> material_array;
708 vtkSmartPointer<vtkIntArray> boundary_array;
709 vtkSmartPointer<vtkIntArray> manifold_array;
711 const bool output_material = !material_id_field.empty();
712 const bool output_boundary = !boundary_id_field.empty();
713 const bool output_manifold = !manifold_id_field.empty();
717 material_array = vtkSmartPointer<vtkIntArray>::New();
718 material_array->SetName(material_id_field.c_str());
719 material_array->SetNumberOfComponents(1);
723 boundary_array = vtkSmartPointer<vtkIntArray>::New();
724 boundary_array->SetName(boundary_id_field.c_str());
725 boundary_array->SetNumberOfComponents(1);
729 manifold_array = vtkSmartPointer<vtkIntArray>::New();
730 manifold_array->SetName(manifold_id_field.c_str());
731 manifold_array->SetNumberOfComponents(1);
736 const unsigned int n_vertices = cell->n_vertices();
737 std::vector<vtkIdType> point_ids(n_vertices);
738 for (
unsigned int i = 0; i < n_vertices; ++i)
739 point_ids[i] =
static_cast<vtkIdType
>(cell->vertex_index(i));
741 int vtk_cell_type = -1;
742 if constexpr (dim == 1)
745 ExcMessage(
"Unsupported 1D cell with != 2 vertices."));
746 vtk_cell_type = VTKCellType::VTK_LINE;
748 else if constexpr (dim == 2)
752 vtk_cell_type = VTKCellType::VTK_QUAD;
753 std::swap(point_ids[2], point_ids[3]);
755 else if (n_vertices == 3)
756 vtk_cell_type = VTKCellType::VTK_TRIANGLE;
760 "quads and triangles are supported."));
762 else if constexpr (dim == 3)
766 vtk_cell_type = VTKCellType::VTK_HEXAHEDRON;
767 std::swap(point_ids[2], point_ids[3]);
768 std::swap(point_ids[6], point_ids[7]);
770 else if (n_vertices == 4)
771 vtk_cell_type = VTKCellType::VTK_TETRA;
772 else if (n_vertices == 6)
773 vtk_cell_type = VTKCellType::VTK_WEDGE;
774 else if (n_vertices == 5)
775 vtk_cell_type = VTKCellType::VTK_PYRAMID;
779 "hexes, tets, wedges and pyramids are "
785 grid->InsertNextCell(vtk_cell_type, n_vertices, point_ids.data());
791 material_array->InsertNextValue(
792 static_cast<int>(cell->material_id()));
794 boundary_array->InsertNextValue(0);
796 manifold_array->InsertNextValue(
805 for (
unsigned int f = 0; f < cell->n_faces(); ++f)
806 if (cell->face(f)->at_boundary())
811 const bool include_by_boundary = output_boundary && (face_bid != 0);
812 const bool include_by_manifold =
815 if (!(include_by_boundary || include_by_manifold))
818 const unsigned int nfv = cell->face(f)->n_vertices();
819 std::vector<vtkIdType> face_point_ids(nfv);
820 for (
unsigned int v = 0; v < nfv; ++v)
822 static_cast<vtkIdType
>(cell->face(f)->vertex_index(v));
824 int vtk_face_type = -1;
825 if constexpr (dim == 1)
828 vtk_face_type = VTK_VERTEX;
830 else if constexpr (dim == 2)
831 vtk_face_type = VTK_LINE;
832 else if constexpr (dim == 3)
835 vtk_face_type = VTK_TRIANGLE;
838 vtk_face_type = VTK_QUAD;
839 std::swap(face_point_ids[2], face_point_ids[3]);
845 grid->InsertNextCell(vtk_face_type, nfv, face_point_ids.data());
848 material_array->InsertNextValue(
849 static_cast<int>(cell->material_id()));
851 boundary_array->InsertNextValue(
static_cast<int>(face_bid));
853 manifold_array->InsertNextValue(
static_cast<int>(face_mid));
858 grid->GetCellData()->AddArray(material_array);
860 grid->GetCellData()->AddArray(boundary_array);
862 grid->GetCellData()->AddArray(manifold_array);
954 const double relative_tolerance)
956 vtkSmartPointer<vtkUnstructuredGrid> grid =
959 std::vector<double>
data;
961 vtkPointData *point_data = grid->GetPointData();
964 for (
int i = 0; i < point_data->GetNumberOfArrays(); ++i)
966 vtkDataArray *data_array = point_data->GetArray(i);
969 vtkIdType n_tuples = data_array->GetNumberOfTuples();
970 int n_components = data_array->GetNumberOfComponents();
971 unsigned int current_size =
data.size();
972 data.resize(current_size + n_tuples * n_components, 0.0);
973 for (vtkIdType tuple_idx = 0; tuple_idx < n_tuples; ++tuple_idx)
974 for (
int comp_idx = 0; comp_idx < n_components; ++comp_idx)
975 data[current_size + tuple_idx * n_components + comp_idx] =
976 data_array->GetComponent(tuple_idx, comp_idx);
980 vtkCellData *cell_data = grid->GetCellData();
983 for (
int i = 0; i < cell_data->GetNumberOfArrays(); ++i)
985 vtkDataArray *data_array = cell_data->GetArray(i);
988 vtkIdType n_tuples = data_array->GetNumberOfTuples();
989 int n_components = data_array->GetNumberOfComponents();
990 unsigned int current_size =
data.size();
991 data.resize(current_size + n_tuples * n_components,
true);
992 for (vtkIdType tuple_idx = 0; tuple_idx < n_tuples; ++tuple_idx)
993 for (
int comp_idx = 0; comp_idx < n_components; ++comp_idx)
994 data[current_size + tuple_idx * n_components + comp_idx] =
995 data_array->GetComponent(tuple_idx, comp_idx);
1009 std::vector<std::string> data_names;
1011 vtkSmartPointer<vtkUnstructuredGrid> grid =
1014 vtkCellData *cell_data = grid->GetCellData();
1015 vtkPointData *point_data = grid->GetPointData();
1017 std::vector<std::shared_ptr<FiniteElement<dim, spacedim>>> fe_collection;
1018 std::vector<unsigned int> n_components_collection;
1020 bool is_simplex =
false;
1021 const vtkIdType n_cells_check = grid->GetNumberOfCells();
1022 for (vtkIdType i = 0; i < n_cells_check; ++i)
1024 vtkCell *cell = grid->GetCell(i);
1027 const int cell_type = cell->GetCellType();
1028 if constexpr (dim == 2)
1030 if (cell_type == VTKCellType::VTK_TRIANGLE)
1036 else if constexpr (dim == 3)
1038 if (cell_type == VTKCellType::VTK_TETRA)
1052 for (
int i = 0; i < point_data->GetNumberOfArrays(); ++i)
1054 vtkDataArray *arr = point_data->GetArray(i);
1057 std::string name = arr->GetName();
1058 int n_comp = arr->GetNumberOfComponents();
1063 fe_collection.push_back(
1079 n_components_collection.push_back(n_comp);
1080 data_names.push_back(name);
1084 for (
int i = 0; i < cell_data->GetNumberOfArrays(); ++i)
1086 vtkDataArray *arr = cell_data->GetArray(i);
1089 std::string name = arr->GetName();
1090 int n_comp = arr->GetNumberOfComponents();
1094 fe_collection.push_back(
1104 fe_collection.push_back(
1111 n_components_collection.push_back(n_comp);
1112 data_names.push_back(name);
1117 std::vector<const FiniteElement<dim, spacedim> *> fe_ptrs;
1118 std::vector<unsigned int> multiplicities;
1119 for (
const auto &fe : fe_collection)
1121 fe_ptrs.push_back(fe.get());
1122 multiplicities.push_back(1);
1124 if (fe_ptrs.empty())
1126 std::vector<std::string>());
1129 return std::make_pair(