deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
utilities.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2025 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
13
14#include <algorithm>
15#include <fstream>
16#include <string>
17#include <vector>
18
19#ifdef DEAL_II_WITH_VTK
21
23
24# include <deal.II/fe/fe.h>
25# include <deal.II/fe/fe_dgq.h>
27# include <deal.II/fe/fe_q.h>
29# include <deal.II/fe/fe_system.h>
30
31# include <deal.II/grid/grid_in.h>
37
38# include <deal.II/lac/vector.h>
39
40// Make sure that the VTK version macros are available.
41# include <vtkVersion.h> // Do not convert for module purposes
42
43// VTK_VERSION_CHECK is defined by the header above for 9.3.0 and above, but
44// we provide a fallback older versions.
45# ifndef VTK_VERSION_CHECK
46# define VTK_VERSION_CHECK(major, minor, build) \
47 (10000000000ULL * (major) + 100000000ULL * (minor) + (build))
48# endif
49
50// Normalize the version macro name to cover both the quick (>=9.3) and
51// legacy (<9.3) headers.
52# if defined(VTK_VERSION_NUMBER_QUICK)
53# define DEAL_II_VTK_VERSION_NUMBER VTK_VERSION_NUMBER_QUICK
54# elif defined(VTK_VERSION_NUMBER)
55# define DEAL_II_VTK_VERSION_NUMBER VTK_VERSION_NUMBER
56# else
57# define DEAL_II_VTK_VERSION_NUMBER 0
58# endif
59
60# include <vtkAppendDataSets.h>
61# include <vtkCellData.h>
62# if DEAL_II_VTK_VERSION_NUMBER >= VTK_VERSION_CHECK(9, 3, 0)
63# include <vtkCleanUnstructuredGrid.h>
64# include <vtkGenericDataObjectReader.h>
65# include <vtkXMLGenericDataObjectReader.h>
66# endif
67# include <vtkCellType.h>
68# include <vtkDataArray.h>
69# include <vtkDataSet.h>
70# include <vtkDataSetReader.h>
71# include <vtkFieldData.h>
72# include <vtkGenericCell.h>
73# include <vtkHexahedron.h>
74# include <vtkIntArray.h>
75# include <vtkLine.h>
76# include <vtkPointData.h>
77# include <vtkPoints.h>
78# include <vtkPolyData.h>
79# include <vtkPolygon.h>
80# include <vtkRectilinearGrid.h>
81# include <vtkSmartPointer.h>
82# include <vtkStructuredGrid.h>
83# include <vtkStructuredPoints.h>
84# include <vtkUnstructuredGrid.h>
85# include <vtkUnstructuredGridWriter.h>
86# include <vtkXMLUnstructuredGridWriter.h>
87
88# include <fstream>
89# include <stdexcept>
90
91#endif
92
94
95#ifdef DEAL_II_WITH_VTK
96namespace VTKWrappers
97{
98 namespace internal
99 {
100 vtkSmartPointer<vtkUnstructuredGrid>
102 const vtkSmartPointer<vtkDataObject> &data_object)
103 {
104 auto *dataset = dynamic_cast<vtkDataSet *>(data_object.Get());
105 AssertThrow(dataset,
107 std::string(
108 "VTK file does not contain a supported dataset type: ") +
109 data_object->GetClassName()));
110
111 auto append_filter = vtkSmartPointer<vtkAppendDataSets>::New();
112 append_filter->SetOutputDataSetType(VTK_UNSTRUCTURED_GRID);
113 append_filter->AddInputData(dataset);
114 append_filter->Update();
115
116 auto out = vtkSmartPointer<vtkUnstructuredGrid>::New();
117 out->ShallowCopy(append_filter->GetOutput());
118
119 AssertThrow(out->GetNumberOfPoints() == dataset->GetNumberOfPoints() ||
120 dataset->GetNumberOfPoints() == 0,
121 ExcMessage("VTK dataset could not be converted to an "
122 "unstructured grid without losing points."));
123
124 AssertThrow(out->GetNumberOfCells() == dataset->GetNumberOfCells() ||
125 dataset->GetNumberOfCells() == 0,
126 ExcMessage("VTK dataset could not be converted to an "
127 "unstructured grid without losing cells."));
128
129 out->GetPointData()->PassData(dataset->GetPointData());
130 out->GetCellData()->PassData(dataset->GetCellData());
131 out->GetFieldData()->PassData(dataset->GetFieldData());
132
133 return out;
134 }
135
136
137
138 vtkSmartPointer<vtkUnstructuredGrid>
139 load_vtk_file(const std::string &vtk_filename,
140 const bool cleanup,
141 const double relative_tolerance)
142 {
143 // check that the file exists
144 std::ifstream file(vtk_filename);
145 AssertThrow(file.good(),
146 ExcMessage("VTK file not found: " + vtk_filename));
147
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));
152
153# if DEAL_II_VTK_VERSION_NUMBER >= VTK_VERSION_CHECK(9, 3, 0)
154 if (ext == "vtu")
155 {
156 auto reader = vtkSmartPointer<vtkXMLGenericDataObjectReader>::New();
157 reader->SetFileName(vtk_filename.c_str());
158 reader->Update();
159 data_object = reader->GetOutput();
160 }
161 else if (ext == "vtk")
162 {
163 // Read the legacy VTK format without assuming a specific dataset
164 // type; the result is normalized to an unstructured grid below.
165 auto reader = vtkSmartPointer<vtkGenericDataObjectReader>::New();
166 reader->SetFileName(vtk_filename.c_str());
167 reader->Update();
168 data_object = reader->GetOutput();
169 }
170 else
171 AssertThrow(false,
172 ExcMessage("Unsupported file extension '" + ext +
173 "'. Use '.vtu' for VTK XML format or '.vtk' "
174 "for legacy VTK format."));
175# else
176 AssertThrow(ext == "vtk",
177 ExcMessage("Unsupported file extension '" + ext +
178 "'. With VTK < 9.3 only '.vtk' legacy files can "
179 "be read."));
180 {
181 auto reader = vtkSmartPointer<vtkDataSetReader>::New();
182 reader->SetFileName(vtk_filename.c_str());
183 reader->Update();
184 data_object = reader->GetOutput();
185 }
186# endif
187
188 auto out = convert_to_unstructured_grid(data_object);
189
190# if DEAL_II_VTK_VERSION_NUMBER >= VTK_VERSION_CHECK(9, 3, 0)
191 if (cleanup)
192 {
193 auto cleaner = vtkSmartPointer<vtkCleanUnstructuredGrid>::New();
194 cleaner->SetToleranceIsAbsolute(false);
195 cleaner->SetTolerance(relative_tolerance);
196 cleaner->SetInputData(out);
197 cleaner->Update();
198 out->ShallowCopy(cleaner->GetOutput());
199 }
200# else
201 (void)cleanup; // avoid unused variable warning
202 (void)relative_tolerance; // avoid unused variable warning
203 deallog << "VTK version < 9.3: skipping cleanup step." << std::endl;
204# endif
205
206 AssertThrow(out, ExcMessage("Failed to read VTK file: " + vtk_filename));
207 return out;
208 }
209
210
211
212 void
213 write_vtk(const std::string &vtk_filename,
214 const vtkSmartPointer<vtkDataObject> &data_object)
215 {
216 // Convert any supported dataset to unstructured grid
217 auto grid = convert_to_unstructured_grid(data_object);
218
219 // Determine output format from file extension
220 const auto dot_pos = vtk_filename.find_last_of('.');
221 Assert(dot_pos != std::string::npos,
222 ExcMessage("Cannot determine output format: no file extension "
223 "found in filename '" +
224 vtk_filename + "'."));
225 const std::string ext = vtk_filename.substr(dot_pos + 1);
226
227 if (ext == "vtu")
228 {
229 auto writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
230 writer->SetFileName(vtk_filename.c_str());
231 writer->SetInputData(grid);
232 writer->Write();
233 }
234 else if (ext == "vtk")
235 {
236 auto writer = vtkSmartPointer<vtkUnstructuredGridWriter>::New();
237 writer->SetFileName(vtk_filename.c_str());
238 writer->SetInputData(grid);
239 writer->Write();
240 }
241 else
242 {
243 AssertThrow(false,
244 ExcMessage("Unsupported file extension '" + ext +
245 "'. Use '.vtu' for VTK XML format or '.vtk' "
246 "for legacy VTK format."));
247 }
248 }
249
250 } // namespace internal
251
252
253
254 template <int dim, int spacedim>
255 void
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)
262 {
263 auto &grid = const_cast<vtkUnstructuredGrid &>(unstructured_grid);
264
265 auto get_cell_data_array =
266 [&grid](const std::string &name) -> vtkDataArray * {
267 if (name.empty())
268 return nullptr;
269
270 vtkCellData *cell_data = grid.GetCellData();
271 if (cell_data == nullptr)
272 return nullptr;
273
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."));
279
280 return data_array;
281 };
282
283 // Get observing pointers to the arrays that store
284 // material/boundary/manifold ids.
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);
288
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)
293 return types::material_id(0);
294
295 return static_cast<types::material_id>(material_value);
296 };
297
298 // Read points
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)
303 {
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)
310 coords[d] == 0.0,
312 "VTK grid has non-zero coordinate in unused dimension."));
313 }
314
315 // Read cells
316 std::vector<CellData<dim>> cells;
317 SubCellData subcell_data;
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(
324 n_points, numbers::flat_manifold_id);
325 for (vtkIdType vtk_id = 0; vtk_id < n_cells; ++vtk_id)
326 {
327 vtkCell *cell = grid.GetCell(vtk_id);
328 if (cell->GetCellDimension() < dim)
329 {
330 if (boundary_ids == nullptr && manifold_ids == nullptr)
331 continue;
332
333 auto set_subcell_ids = [vtk_id, boundary_ids, manifold_ids](
334 auto &cell_data) {
335 if (boundary_ids != nullptr)
336 {
337 double bval = boundary_ids->GetComponent(vtk_id, 0);
338 if (bval < 0)
339 cell_data.boundary_id = numbers::internal_face_boundary_id;
340 else
341 cell_data.boundary_id =
342 static_cast<types::boundary_id>(bval);
343 }
344 else if (manifold_ids != nullptr)
345 cell_data.boundary_id = numbers::internal_face_boundary_id;
346
347 if (manifold_ids != nullptr)
348 {
349 double mval = manifold_ids->GetComponent(vtk_id, 0);
350 if (mval < 0)
351 cell_data.manifold_id = numbers::flat_manifold_id;
352 else
353 cell_data.manifold_id =
354 static_cast<types::manifold_id>(mval);
355 }
356 };
357
358 if constexpr (dim == 1)
359 {
360 if (cell->GetCellType() == VTK_VERTEX)
361 {
363 cell->GetNumberOfPoints() == 1,
365 "Only vertex subcells with 1 point are supported."));
366
367 const vtkIdType vertex_index = cell->GetPointId(0);
368 AssertIndexRange(vertex_index, n_points);
369
370 if (boundary_ids != nullptr)
371 {
372 const double boundary_value =
373 boundary_ids->GetComponent(vtk_id, 0);
374 if (boundary_value >= 0)
375 {
376 boundary_vertex_ids_present[vertex_index] = true;
377 boundary_vertex_ids[vertex_index] =
378 static_cast<types::boundary_id>(boundary_value);
379 }
380 }
381
382 if (manifold_ids != nullptr)
383 {
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 ?
390 static_cast<types::manifold_id>(manifold_value));
391 }
392 }
393 }
394 else if constexpr (dim == 2)
395 {
396 if (cell->GetCellType() == VTK_LINE)
397 {
399 cell->GetNumberOfPoints() == 2,
401 "Only line subcells with 2 points are supported."));
402
403 CellData<1> cell_data(2);
404 for (unsigned int j = 0; j < 2; ++j)
405 cell_data.vertices[j] = cell->GetPointId(j);
406
407 set_subcell_ids(cell_data);
408 subcell_data.boundary_lines.push_back(cell_data);
409 }
410 }
411 else if constexpr (dim == 3)
412 {
413 if (cell->GetCellType() == VTK_LINE)
414 {
416 cell->GetNumberOfPoints() == 2,
418 "Only line subcells with 2 points are supported."));
419
420 CellData<1> cell_data(2);
421 for (unsigned int j = 0; j < 2; ++j)
422 cell_data.vertices[j] = cell->GetPointId(j);
423
424 set_subcell_ids(cell_data);
425 subcell_data.boundary_lines.push_back(cell_data);
426 }
427 else if (cell->GetCellType() == VTK_TRIANGLE ||
428 cell->GetCellType() == VTK_QUAD)
429 {
431 cell->GetNumberOfPoints() == 3 ||
432 cell->GetNumberOfPoints() == 4,
434 "Only triangle and quad face subcells are supported."));
435
436 CellData<2> cell_data(cell->GetNumberOfPoints());
437 for (unsigned int j = 0; j < cell->GetNumberOfPoints(); ++j)
438 cell_data.vertices[j] = cell->GetPointId(j);
439
440 // VTK and deal.II use different vertex ordering for
441 // quad faces; match deal.II ordering by swapping the
442 // last two vertices (consistent with quad cell handling
443 // above).
444 if (cell->GetCellType() == VTK_QUAD)
445 std::swap(cell_data.vertices[2], cell_data.vertices[3]);
446
447 set_subcell_ids(cell_data);
448 subcell_data.boundary_quads.push_back(cell_data);
449 }
450 }
451
452 continue;
453 }
454
455 AssertThrow(cell->GetCellDimension() == dim,
456 ExcMessage("Unsupported VTK cell dimension."));
457
458 if constexpr (dim == 1)
459 {
460 if (cell->GetCellType() != VTKCellType::VTK_LINE)
461 AssertThrow(false,
463 "Unsupported cell type in 1D VTK file: only "
464 "VTK_LINE is supported."));
465 AssertThrow(cell->GetNumberOfPoints() == 2,
467 "Only line cells with 2 points are supported."));
468 CellData<1> cell_data(2);
469 for (unsigned int j = 0; j < 2; ++j)
470 cell_data.vertices[j] = cell->GetPointId(j);
471 cell_data.material_id = 0;
472 if (material_ids != nullptr)
473 cell_data.material_id =
474 material_id_from_vtk(*material_ids, vtk_id);
475 if (boundary_ids != nullptr)
476 {
477 const double boundary_value =
478 boundary_ids->GetComponent(vtk_id, 0);
479 if (boundary_value >= 0)
480 cell_data.boundary_id =
481 static_cast<types::boundary_id>(boundary_value);
482 }
483 if (manifold_ids != nullptr)
484 {
485 double mval = manifold_ids->GetComponent(vtk_id, 0);
486 if (mval < 0)
488 else
489 cell_data.manifold_id = static_cast<types::manifold_id>(mval);
490 }
491 cells.push_back(cell_data);
492 }
493 else if constexpr (dim == 2)
494 {
495 if (cell->GetCellType() == VTKCellType::VTK_QUAD)
496 {
497 AssertThrow(cell->GetNumberOfPoints() == 4,
499 "Only quad cells with 4 points are supported."));
500 CellData<2> cell_data(4);
501 for (unsigned int j = 0; j < 4; ++j)
502 cell_data.vertices[j] = cell->GetPointId(j);
503 std::swap(cell_data.vertices[2], cell_data.vertices[3]);
504 cell_data.material_id = 0;
505 if (material_ids != nullptr)
506 cell_data.material_id =
507 material_id_from_vtk(*material_ids, vtk_id);
508 if (manifold_ids != nullptr)
509 {
510 double mval = manifold_ids->GetComponent(vtk_id, 0);
511 if (mval < 0)
513 else
514 cell_data.manifold_id =
515 static_cast<types::manifold_id>(mval);
516 }
517 cells.push_back(cell_data);
518 }
519 else if (cell->GetCellType() == VTKCellType::VTK_TRIANGLE)
520 {
522 cell->GetNumberOfPoints() == 3,
524 "Only triangle cells with 3 points are supported."));
525 CellData<2> cell_data(3);
526 for (unsigned int j = 0; j < 3; ++j)
527 cell_data.vertices[j] = cell->GetPointId(j);
528 cell_data.material_id = 0;
529 if (material_ids != nullptr)
530 cell_data.material_id =
531 material_id_from_vtk(*material_ids, vtk_id);
532 if (manifold_ids != nullptr)
533 {
534 double mval = manifold_ids->GetComponent(vtk_id, 0);
535 if (mval < 0)
537 else
538 cell_data.manifold_id =
539 static_cast<types::manifold_id>(mval);
540 }
541 cells.push_back(cell_data);
542 }
543 else
544 AssertThrow(false,
546 "Unsupported cell type in 2D VTK file: only "
547 "VTK_QUAD and VTK_TRIANGLE are supported."));
548 }
549 else if constexpr (dim == 3)
550 {
551 if (cell->GetCellType() == VTKCellType::VTK_HEXAHEDRON)
552 {
553 AssertThrow(cell->GetNumberOfPoints() == 8,
555 "Only hex cells with 8 points are supported."));
556 CellData<3> cell_data(8);
557 for (unsigned int j = 0; j < 8; ++j)
558 cell_data.vertices[j] = cell->GetPointId(j);
559 cell_data.material_id = 0;
560 // Numbering of vertices in VTK files is different from
561 // deal.II
562 std::swap(cell_data.vertices[2], cell_data.vertices[3]);
563 std::swap(cell_data.vertices[6], cell_data.vertices[7]);
564 if (material_ids != nullptr)
565 cell_data.material_id =
566 material_id_from_vtk(*material_ids, vtk_id);
567 if (manifold_ids != nullptr)
568 {
569 double mval = manifold_ids->GetComponent(vtk_id, 0);
570 if (mval < 0)
572 else
573 cell_data.manifold_id =
574 static_cast<types::manifold_id>(mval);
575 }
576 cells.push_back(cell_data);
577 }
578 else if (cell->GetCellType() == VTKCellType::VTK_TETRA)
579 {
581 cell->GetNumberOfPoints() == 4,
583 "Only tetrahedron cells with 4 points are supported."));
584 CellData<3> cell_data(4);
585 for (unsigned int j = 0; j < 4; ++j)
586 cell_data.vertices[j] = cell->GetPointId(j);
587 cell_data.material_id = 0;
588 if (material_ids != nullptr)
589 cell_data.material_id =
590 material_id_from_vtk(*material_ids, vtk_id);
591 if (manifold_ids != nullptr)
592 {
593 double mval = manifold_ids->GetComponent(vtk_id, 0);
594 if (mval < 0)
596 else
597 cell_data.manifold_id =
598 static_cast<types::manifold_id>(mval);
599 }
600 cells.push_back(cell_data);
601 }
602 else if (cell->GetCellType() == VTKCellType::VTK_WEDGE)
603 {
604 AssertThrow(cell->GetNumberOfPoints() == 6,
606 "Only prism cells with 6 points are supported."));
607 CellData<3> cell_data(6);
608 for (unsigned int j = 0; j < 6; ++j)
609 cell_data.vertices[j] = cell->GetPointId(j);
610 cell_data.material_id = 0;
611 if (material_ids != nullptr)
612 cell_data.material_id =
613 material_id_from_vtk(*material_ids, vtk_id);
614 if (manifold_ids != nullptr)
615 {
616 double mval = manifold_ids->GetComponent(vtk_id, 0);
617 if (mval < 0)
619 else
620 cell_data.manifold_id =
621 static_cast<types::manifold_id>(mval);
622 }
623 cells.push_back(cell_data);
624 }
625 else if (cell->GetCellType() == VTKCellType::VTK_PYRAMID)
626 {
628 cell->GetNumberOfPoints() == 5,
630 "Only pyramid cells with 5 points are supported."));
631 CellData<3> cell_data(5);
632 for (unsigned int j = 0; j < 5; ++j)
633 cell_data.vertices[j] = cell->GetPointId(j);
634 cell_data.material_id = 0;
635 if (material_ids != nullptr)
636 cell_data.material_id =
637 material_id_from_vtk(*material_ids, vtk_id);
638 if (manifold_ids != nullptr)
639 {
640 double mval = manifold_ids->GetComponent(vtk_id, 0);
641 if (mval < 0)
643 else
644 cell_data.manifold_id =
645 static_cast<types::manifold_id>(mval);
646 }
647 cells.push_back(cell_data);
648 }
649 else
651 false,
653 "Unsupported cell type in 3D VTK file: only "
654 "VTK_HEXAHEDRON, VTK_TETRA, VTK_WEDGE, and VTK_PYRAMID are supported."));
655 }
656 else
657 {
658 AssertThrow(false, ExcMessage("Unsupported dimension."));
659 }
660 }
661
662 // Create triangulation
663 tria.create_triangulation(points, cells, subcell_data);
664
665 if constexpr (dim == 1)
666 for (const auto &cell : tria.active_cell_iterators())
667 for (unsigned int f = 0; f < cell->n_faces(); ++f)
668 if (cell->face(f)->at_boundary())
669 {
670 const unsigned int vertex_index = cell->face(f)->vertex_index(0);
671
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]);
678 }
679 }
680
681
682
683 template <int dim, int spacedim>
684 vtkSmartPointer<vtkUnstructuredGrid>
687 const std::string &material_id_field,
688 const std::string &boundary_id_field,
689 const std::string &manifold_id_field)
690 {
691 auto grid = vtkSmartPointer<vtkUnstructuredGrid>::New();
692 auto points = vtkSmartPointer<vtkPoints>::New();
693
694 points->SetNumberOfPoints(tria.n_vertices());
695 for (unsigned int i = 0; i < tria.n_vertices(); ++i)
696 {
697 std::array<double, 3> coords = {{0.0, 0.0, 0.0}};
698 for (unsigned int d = 0; d < spacedim; ++d)
699 coords[d] = tria.get_vertices()[i][d];
700 points->SetPoint(i, coords.data());
701 }
702 grid->SetPoints(points);
703
704 // Prepare optional VTK arrays for cell data. We will append values in
705 // the same order that we insert cells: first full (codim 0) cells, then
706 // any subcells (boundary faces) we append afterwards.
707 vtkSmartPointer<vtkIntArray> material_array;
708 vtkSmartPointer<vtkIntArray> boundary_array;
709 vtkSmartPointer<vtkIntArray> manifold_array;
710
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();
714
715 if (output_material)
716 {
717 material_array = vtkSmartPointer<vtkIntArray>::New();
718 material_array->SetName(material_id_field.c_str());
719 material_array->SetNumberOfComponents(1);
720 }
721 if (output_boundary)
722 {
723 boundary_array = vtkSmartPointer<vtkIntArray>::New();
724 boundary_array->SetName(boundary_id_field.c_str());
725 boundary_array->SetNumberOfComponents(1);
726 }
727 if (output_manifold)
728 {
729 manifold_array = vtkSmartPointer<vtkIntArray>::New();
730 manifold_array->SetName(manifold_id_field.c_str());
731 manifold_array->SetNumberOfComponents(1);
732 }
733
734 for (const auto &cell : tria.active_cell_iterators())
735 {
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));
740
741 int vtk_cell_type = -1;
742 if constexpr (dim == 1)
743 {
744 AssertThrow(n_vertices == 2,
745 ExcMessage("Unsupported 1D cell with != 2 vertices."));
746 vtk_cell_type = VTKCellType::VTK_LINE;
747 }
748 else if constexpr (dim == 2)
749 {
750 if (n_vertices == 4)
751 {
752 vtk_cell_type = VTKCellType::VTK_QUAD;
753 std::swap(point_ids[2], point_ids[3]);
754 }
755 else if (n_vertices == 3)
756 vtk_cell_type = VTKCellType::VTK_TRIANGLE;
757 else
758 AssertThrow(false,
759 ExcMessage("Unsupported 2D cell type: only "
760 "quads and triangles are supported."));
761 }
762 else if constexpr (dim == 3)
763 {
764 if (n_vertices == 8)
765 {
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]);
769 }
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;
776 else
777 AssertThrow(false,
778 ExcMessage("Unsupported 3D cell type: only "
779 "hexes, tets, wedges and pyramids are "
780 "supported."));
781 }
782 else
783 AssertThrow(false, ExcMessage("Unsupported dimension."));
784
785 grid->InsertNextCell(vtk_cell_type, n_vertices, point_ids.data());
786
787 // For codim-0 cells write material/manifold values and zero for
788 // boundary id, since boundary ids live on codim-1 subcells (vertices
789 // in 1d, edges in 2d, and faces in 3d).
790 if (output_material)
791 material_array->InsertNextValue(
792 static_cast<int>(cell->material_id()));
793 if (output_boundary)
794 boundary_array->InsertNextValue(0);
795 if (output_manifold)
796 manifold_array->InsertNextValue(
797 static_cast<int>(numbers::flat_manifold_id));
798 }
799
800 // Append boundary faces as separate VTK cells when requested. We only
801 // append codim-1 subcells (faces / edges) for which either a non-zero
802 // boundary id or a non-flat manifold id exists, depending on the fields
803 // requested by the caller.
804 for (const auto &cell : tria.active_cell_iterators())
805 for (unsigned int f = 0; f < cell->n_faces(); ++f)
806 if (cell->face(f)->at_boundary())
807 {
808 const types::boundary_id face_bid = cell->face(f)->boundary_id();
809 const types::manifold_id face_mid = cell->face(f)->manifold_id();
810
811 const bool include_by_boundary = output_boundary && (face_bid != 0);
812 const bool include_by_manifold =
813 output_manifold && (face_mid != numbers::flat_manifold_id);
814
815 if (!(include_by_boundary || include_by_manifold))
816 continue;
817
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)
821 face_point_ids[v] =
822 static_cast<vtkIdType>(cell->face(f)->vertex_index(v));
823
824 int vtk_face_type = -1;
825 if constexpr (dim == 1)
826 {
827 Assert(nfv == 1, ExcInternalError());
828 vtk_face_type = VTK_VERTEX;
829 }
830 else if constexpr (dim == 2)
831 vtk_face_type = VTK_LINE;
832 else if constexpr (dim == 3)
833 {
834 if (nfv == 3)
835 vtk_face_type = VTK_TRIANGLE;
836 else if (nfv == 4)
837 {
838 vtk_face_type = VTK_QUAD;
839 std::swap(face_point_ids[2], face_point_ids[3]);
840 }
841 else
843 }
844
845 grid->InsertNextCell(vtk_face_type, nfv, face_point_ids.data());
846
847 if (output_material)
848 material_array->InsertNextValue(
849 static_cast<int>(cell->material_id()));
850 if (output_boundary)
851 boundary_array->InsertNextValue(static_cast<int>(face_bid));
852 if (output_manifold)
853 manifold_array->InsertNextValue(static_cast<int>(face_mid));
854 }
855
856 // Attach arrays to cell data if created.
857 if (output_material)
858 grid->GetCellData()->AddArray(material_array);
859 if (output_boundary)
860 grid->GetCellData()->AddArray(boundary_array);
861 if (output_manifold)
862 grid->GetCellData()->AddArray(manifold_array);
863
864 return grid;
865 }
866
867
868
869 template <int dim, int spacedim>
870 void
871 write_vtk(const std::string &vtk_filename,
873 const std::string &material_id_field,
874 const std::string &boundary_id_field,
875 const std::string &manifold_id_field)
876 {
878 tria, material_id_field, boundary_id_field, manifold_id_field);
879 internal::write_vtk(vtk_filename, grid);
880 }
881
882
883
884 template <int dim, int spacedim>
885 void
886 read_tria(const std::string &vtk_filename,
888 const bool cleanup,
889 const double relative_tolerance,
890 const std::string &material_id_field,
891 const std::string &boundary_id_field,
892 const std::string &manifold_id_field)
893 {
894 vtkSmartPointer<vtkUnstructuredGrid> grid =
895 internal::load_vtk_file(vtk_filename, cleanup, relative_tolerance);
897 *grid, tria, material_id_field, boundary_id_field, manifold_id_field);
898 }
899
900
901
902 void
903 read_cell_data(const std::string &vtk_filename,
904 const std::string &cell_data_name,
905 Vector<double> &output_vector,
906 const bool cleanup,
907 const double relative_tolerance)
908 {
909 vtkSmartPointer<vtkUnstructuredGrid> grid =
910 internal::load_vtk_file(vtk_filename, cleanup, relative_tolerance);
911 vtkDataArray *data_array =
912 grid->GetCellData()->GetArray(cell_data_name.c_str());
913 AssertThrow(data_array,
914 ExcMessage("Cell data array '" + cell_data_name +
915 "' not found in VTK file: " + vtk_filename));
916 vtkIdType n_tuples = data_array->GetNumberOfTuples();
917 int n_components = data_array->GetNumberOfComponents();
918 output_vector.reinit(n_tuples * n_components);
919 for (vtkIdType i = 0; i < n_tuples; ++i)
920 for (int j = 0; j < n_components; ++j)
921 output_vector[i * n_components + j] = data_array->GetComponent(i, j);
922 }
923
924
925
926 void
927 read_vertex_data(const std::string &vtk_filename,
928 const std::string &vertex_data_name,
929 Vector<double> &output_vector,
930 const bool cleanup,
931 const double relative_tolerance)
932 {
933 vtkSmartPointer<vtkUnstructuredGrid> grid =
934 internal::load_vtk_file(vtk_filename, cleanup, relative_tolerance);
935 vtkDataArray *data_array =
936 grid->GetPointData()->GetArray(vertex_data_name.c_str());
937 AssertThrow(data_array,
938 ExcMessage("Point data array '" + vertex_data_name +
939 "' not found in VTK file: " + vtk_filename));
940 vtkIdType n_tuples = data_array->GetNumberOfTuples();
941 int n_components = data_array->GetNumberOfComponents();
942 output_vector.reinit(n_tuples * n_components);
943 for (vtkIdType i = 0; i < n_tuples; ++i)
944 for (int j = 0; j < n_components; ++j)
945 output_vector[i * n_components + j] = data_array->GetComponent(i, j);
946 }
947
948
949
950 void
951 read_all_data(const std::string &vtk_filename,
952 Vector<double> &output_vector,
953 const bool cleanup,
954 const double relative_tolerance)
955 {
956 vtkSmartPointer<vtkUnstructuredGrid> grid =
957 internal::load_vtk_file(vtk_filename, cleanup, relative_tolerance);
958
959 std::vector<double> data;
960
961 vtkPointData *point_data = grid->GetPointData();
962 if (point_data)
963 {
964 for (int i = 0; i < point_data->GetNumberOfArrays(); ++i)
965 {
966 vtkDataArray *data_array = point_data->GetArray(i);
967 if (!data_array)
968 continue;
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);
977 }
978 }
979
980 vtkCellData *cell_data = grid->GetCellData();
981 if (cell_data)
982 {
983 for (int i = 0; i < cell_data->GetNumberOfArrays(); ++i)
984 {
985 vtkDataArray *data_array = cell_data->GetArray(i);
986 if (!data_array)
987 continue;
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);
996 }
997 }
998 output_vector.reinit(data.size());
999 std::copy(data.begin(), data.end(), output_vector.begin());
1000 }
1001
1002
1003
1004 template <int dim, int spacedim>
1005 std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>,
1006 std::vector<std::string>>
1007 vtk_to_finite_element(const std::string &vtk_filename)
1008 {
1009 std::vector<std::string> data_names;
1010
1011 vtkSmartPointer<vtkUnstructuredGrid> grid =
1012 internal::load_vtk_file(vtk_filename, false, 0.0);
1013
1014 vtkCellData *cell_data = grid->GetCellData();
1015 vtkPointData *point_data = grid->GetPointData();
1016
1017 std::vector<std::shared_ptr<FiniteElement<dim, spacedim>>> fe_collection;
1018 std::vector<unsigned int> n_components_collection;
1019
1020 bool is_simplex = false;
1021 const vtkIdType n_cells_check = grid->GetNumberOfCells();
1022 for (vtkIdType i = 0; i < n_cells_check; ++i)
1023 {
1024 vtkCell *cell = grid->GetCell(i);
1025 if (!cell)
1026 continue;
1027 const int cell_type = cell->GetCellType();
1028 if constexpr (dim == 2)
1029 {
1030 if (cell_type == VTKCellType::VTK_TRIANGLE)
1031 {
1032 is_simplex = true;
1033 break;
1034 }
1035 }
1036 else if constexpr (dim == 3)
1037 {
1038 if (cell_type == VTKCellType::VTK_TETRA)
1039 {
1040 is_simplex = true;
1041 break;
1042 }
1043 }
1044 else
1045 {
1046 is_simplex = false;
1047 break;
1048 }
1049 }
1050
1051 // Query point data fields
1052 for (int i = 0; i < point_data->GetNumberOfArrays(); ++i)
1053 {
1054 vtkDataArray *arr = point_data->GetArray(i);
1055 if (!arr)
1056 continue;
1057 std::string name = arr->GetName();
1058 int n_comp = arr->GetNumberOfComponents();
1059
1060 if (is_simplex)
1061 {
1062 if (n_comp == 1)
1063 fe_collection.push_back(
1064 std::make_shared<FE_SimplexP<dim, spacedim>>(1));
1065 else
1066 // Use FESystem for vector fields
1067 fe_collection.push_back(std::make_shared<FESystem<dim, spacedim>>(
1068 FE_SimplexP<dim, spacedim>(1), n_comp));
1069 }
1070 else
1071 {
1072 if (n_comp == 1)
1073 fe_collection.push_back(std::make_shared<FE_Q<dim, spacedim>>(1));
1074 else
1075 // Use FESystem for vector fields
1076 fe_collection.push_back(std::make_shared<FESystem<dim, spacedim>>(
1077 FE_Q<dim, spacedim>(1), n_comp));
1078 }
1079 n_components_collection.push_back(n_comp);
1080 data_names.push_back(name);
1081 }
1082
1083 // Query cell data fields
1084 for (int i = 0; i < cell_data->GetNumberOfArrays(); ++i)
1085 {
1086 vtkDataArray *arr = cell_data->GetArray(i);
1087 if (!arr)
1088 continue;
1089 std::string name = arr->GetName();
1090 int n_comp = arr->GetNumberOfComponents();
1091 if (is_simplex)
1092 {
1093 if (n_comp == 1)
1094 fe_collection.push_back(
1095 std::make_shared<FE_SimplexDGP<dim, spacedim>>(0));
1096 else
1097 // Use FESystem for vector fields
1098 fe_collection.push_back(std::make_shared<FESystem<dim, spacedim>>(
1099 FE_SimplexDGP<dim, spacedim>(0), n_comp));
1100 }
1101 else
1102 {
1103 if (n_comp == 1)
1104 fe_collection.push_back(
1105 std::make_shared<FE_DGQ<dim, spacedim>>(0));
1106 else
1107 // Use FESystem for vector fields
1108 fe_collection.push_back(std::make_shared<FESystem<dim, spacedim>>(
1109 FE_DGQ<dim, spacedim>(0), n_comp));
1110 }
1111 n_components_collection.push_back(n_comp);
1112 data_names.push_back(name);
1113 }
1114
1115
1116 // Build a FESystem with all fields
1117 std::vector<const FiniteElement<dim, spacedim> *> fe_ptrs;
1118 std::vector<unsigned int> multiplicities;
1119 for (const auto &fe : fe_collection)
1120 {
1121 fe_ptrs.push_back(fe.get());
1122 multiplicities.push_back(1);
1123 }
1124 if (fe_ptrs.empty())
1125 return std::make_pair(std::make_unique<FE_Nothing<dim, spacedim>>(),
1126 std::vector<std::string>());
1127 else
1128 {
1129 return std::make_pair(
1130 std::make_unique<FESystem<dim, spacedim>>(fe_ptrs, multiplicities),
1131 data_names);
1132 }
1133 }
1134
1135
1136 template <int dim, int spacedim>
1137 void
1138 read_vtk(const std::string &vtk_filename,
1139 DoFHandler<dim, spacedim> &dof_handler,
1140 Vector<double> &output_vector,
1141 std::vector<std::string> &data_names,
1142 const bool cleanup,
1143 const double relative_tolerance,
1144 const std::string &material_id_field,
1145 const std::string &boundary_id_field,
1146 const std::string &manifold_id_field)
1147 {
1148 // Get a non-const reference to the triangulation
1149 auto &tria = const_cast<Triangulation<dim, spacedim> &>(
1150 dof_handler.get_triangulation());
1151
1152 // Make sure the triangulation is actually a serial triangulation
1153 auto parallel_tria =
1154 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria);
1155 AssertThrow(parallel_tria == nullptr,
1156 ExcMessage(
1157 "The input triangulation must be a serial triangulation."));
1158
1159 // Clear the triangulation to ensure it is empty before reading
1160 tria.clear();
1161 // Read the mesh from the VTK file
1162 read_tria(vtk_filename,
1163 tria,
1164 cleanup,
1165 relative_tolerance,
1166 material_id_field,
1167 boundary_id_field,
1168 manifold_id_field);
1169
1170 Vector<double> raw_data_vector;
1171 read_all_data(vtk_filename, raw_data_vector, cleanup, relative_tolerance);
1172
1173 auto [fe, data_names_from_fe] =
1174 vtk_to_finite_element<dim, spacedim>(vtk_filename);
1175
1176 dof_handler.distribute_dofs(*fe);
1177 output_vector.reinit(dof_handler.n_dofs());
1178 data_to_dealii_vector(tria, raw_data_vector, dof_handler, output_vector);
1179
1180 AssertDimension(dof_handler.n_dofs(), output_vector.size());
1181 AssertDimension(dof_handler.get_fe().n_blocks(), data_names_from_fe.size());
1182 data_names = data_names_from_fe;
1183 }
1184
1185# include "vtk/utilities.inst"
1186
1187} // namespace VTKWrappers
1188
1189#endif
1190
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
Definition fe_q.h:552
unsigned int n_blocks() const
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
const std::vector< Point< spacedim > > & get_vertices() const
unsigned int n_vertices() const
virtual size_type size() const override
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
#define DEAL_II_ASSERT_UNREACHABLE()
IteratorRange< active_cell_iterator > active_cell_iterators() const
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
LogStream deallog
Definition logstream.cc:36
std::vector< index_type > data
Definition mpi.cc:734
vtkSmartPointer< vtkUnstructuredGrid > convert_to_unstructured_grid(const vtkSmartPointer< vtkDataObject > &data_object)
Definition utilities.cc:101
void write_vtk(const std::string &vtk_filename, const vtkSmartPointer< vtkDataObject > &data_object)
Definition utilities.cc:213
vtkSmartPointer< vtkUnstructuredGrid > load_vtk_file(const std::string &vtk_filename, const bool cleanup=true, const double relative_tolerance=0.0)
Definition utilities.cc:139
void read_tria(const std::string &vtk_filename, Triangulation< dim, spacedim > &tria, const bool cleanup=true, const double relative_tolerance=0.0, const std::string &material_id_field="", const std::string &boundary_id_field="", const std::string &manifold_id_field="")
Definition utilities.cc:886
void read_vtk(const std::string &vtk_filename, DoFHandler< dim, spacedim > &dof_handler, Vector< double > &output_vector, std::vector< std::string > &data_names, const bool cleanup=true, const double relative_tolerance=0.0, const std::string &material_id_field="", const std::string &boundary_id_field="", const std::string &manifold_id_field="")
void read_cell_data(const std::string &vtk_filename, const std::string &cell_data_name, Vector< double > &output_vector, const bool cleanup=true, const double relative_tolerance=0.0)
Definition utilities.cc:903
void read_vertex_data(const std::string &vtk_filename, const std::string &vertex_data_name, Vector< double > &output_vector, const bool cleanup=true, const double relative_tolerance=0.0)
Definition utilities.cc:927
vtkSmartPointer< vtkUnstructuredGrid > dealii_triangulation_to_unstructured_grid(const Triangulation< dim, spacedim > &tria, const std::string &material_id_field="", const std::string &boundary_id_field="", const std::string &manifold_id_field="")
Definition utilities.cc:685
void write_vtk(const std::string &vtk_filename, const Triangulation< dim, spacedim > &tria, const std::string &material_id_field="", const std::string &boundary_id_field="", const std::string &manifold_id_field="")
Definition utilities.cc:871
void unstructured_grid_to_dealii_triangulation(const vtkUnstructuredGrid &unstructured_grid, Triangulation< dim, spacedim > &tria, const std::string &material_id_field="", const std::string &boundary_id_field="", const std::string &manifold_id_field="")
Definition utilities.cc:256
void data_to_dealii_vector(const Triangulation< dim, spacedim > &serial_tria, const Vector< double > &data, const DoFHandler< dim, spacedim > &dh, VectorType &output_vector)
void read_all_data(const std::string &vtk_filename, Vector< double > &output_vector, const bool cleanup=true, const double relative_tolerance=0.0)
Definition utilities.cc:951
std::pair< std::unique_ptr< FiniteElement< dim, spacedim > >, std::vector< std::string > > vtk_to_finite_element(const std::string &vtk_filename)
constexpr types::boundary_id internal_face_boundary_id
Definition types.h:319
constexpr types::manifold_id flat_manifold_id
Definition types.h:332
unsigned int material_id
Definition types.h:182
types::manifold_id manifold_id
Definition cell_data.h:125
std_cxx26::inplace_vector< unsigned int, ReferenceCells::max_n_vertices< structdim >()> vertices
Definition cell_data.h:84
types::material_id material_id
Definition cell_data.h:103
types::boundary_id boundary_id
Definition cell_data.h:114
std::vector< CellData< 2 > > boundary_quads
Definition cell_data.h:247
std::vector< CellData< 1 > > boundary_lines
Definition cell_data.h:231