GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/io/vtk/vtkreader.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 194 217 89.4%
Functions: 26 34 76.5%
Branches: 170 860 19.8%

Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup InputOutput
10 * \brief A vtk file reader using tinyxml2 as xml backend
11 */
12 #ifndef DUMUX_IO_VTK_VTKREADER_HH
13 #define DUMUX_IO_VTK_VTKREADER_HH
14
15 #include <iostream>
16 #include <iterator>
17 #include <algorithm>
18 #include <memory>
19 #include <type_traits>
20 #include <unordered_map>
21 #include <utility>
22 #include <filesystem>
23
24 #include <dune/common/parallel/mpihelper.hh>
25 #include <dune/common/exceptions.hh>
26 #include <dune/grid/common/capabilities.hh>
27 #include <dune/grid/io/file/vtk/common.hh>
28 #include <dune/grid/common/gridfactory.hh>
29
30 #include <dumux/io/container.hh>
31
32 #if DUMUX_HAVE_GRIDFORMAT
33
34 #include <gridformat/gridformat.hpp>
35 #include <gridformat/traits/dune.hpp>
36 #include <gridformat/reader.hpp>
37 #include <gridformat/decorators/reader_polylines_subdivider.hpp>
38
39 namespace Dumux::Detail::VTKReader {
40
41 template<GridFormat::Concepts::Communicator C>
42 auto makeGridformatReaderFactory(const C& c)
43 {
44 return [fac = GridFormat::AnyReaderFactory<C>{c}] (const std::string& f)
45 -> std::unique_ptr<GridFormat::GridReader> {
46 // use adapter for poly data that splits polylines into segments
47 if (f.ends_with("vtp"))
48 return std::make_unique<GridFormat::ReaderDecorators::PolylinesSubdivider>(fac.make_for(f));
49 return fac.make_for(f);
50 };
51 }
52
53 } // end namespace Dumux::Detail::VTKReader
54
55 namespace Dumux {
56
57 /*!
58 * \ingroup InputOutput
59 * \brief A vtk file reader using gridformat
60 */
61 class VTKReader
62 {
63 public:
64 enum class DataType { cellData, pointData, fieldData };
65
66 //! the cell / point data type for point data read from a grid file
67 using Data = std::unordered_map<std::string, std::vector<double>>;
68
69 explicit VTKReader(const std::string& fileName)
70 : reader_(GridFormat::Reader::from(fileName,
71 Detail::VTKReader::makeGridformatReaderFactory(
72 Dune::MPIHelper::instance().getCommunicator()
73 )
74 ))
75 {}
76
77 /*!
78 * \brief Reviews data from the vtk file to check if there is a data array with a specified name
79 * \param name the name attribute of the data array to read
80 * \param type the data array type
81 */
82 bool hasData(const std::string& name, const DataType& type) const
83 {
84 if (type == DataType::cellData)
85 return std::ranges::any_of(
86 cell_field_names(reader_),
87 [&] (const auto& n) { return name == n; }
88 );
89 else if (type == DataType::pointData)
90 return std::ranges::any_of(
91 point_field_names(reader_),
92 [&] (const auto& n) { return name == n; }
93 );
94 else if (type == DataType::fieldData)
95 return std::ranges::any_of(
96 meta_data_field_names(reader_),
97 [&] (const auto& n) { return name == n; }
98 );
99 else
100 DUNE_THROW(Dune::IOError, "Unknown data type for field '" << name << "'");
101 }
102
103 /*!
104 * \brief read data from the vtk file to a container, e.g. std::vector<double>
105 * \tparam Container a container type that has begin(), end(), push_back(), e.g. std::vector<>
106 * \param name the name attribute of the data array to read
107 * \param type the data array type
108 */
109 template<class Container>
110 Container readData(const std::string& name, const DataType& type) const
111 {
112 if (type == DataType::cellData)
113 {
114 Container values(reader_.number_of_cells());
115 reader_.cell_field(name)->export_to(values);
116 return values;
117 }
118 else if (type == DataType::pointData)
119 {
120 Container values(reader_.number_of_points());
121 reader_.point_field(name)->export_to(values);
122 return values;
123 }
124 else if (type == DataType::fieldData)
125 {
126 DUNE_THROW(Dune::NotImplemented, "Reading meta data not yet implemented");
127 }
128 else
129 DUNE_THROW(Dune::IOError, "Unknown data type for field '" << name << "'");
130 }
131
132 /*!
133 * \brief Read a grid from a vtk/vtu/vtp file, ignoring cell and point data
134 * \param verbose if the output should be verbose
135 */
136 template<class Grid>
137 std::unique_ptr<Grid> readGrid(bool verbose = false) const
138 {
139 Dune::GridFactory<Grid> factory;
140 return readGrid(factory, verbose);
141 }
142
143 /*!
144 * \brief Read a grid from a vtk/vtu/vtp file, ignoring cell and point data
145 * \note use this signature if the factory might be needed outside to interpret the data via the factory's insertion indices
146 * \param verbose if the output should be verbose
147 * \param factory the (empty) grid factory
148 */
149 template<class Grid>
150 std::unique_ptr<Grid> readGrid(Dune::GridFactory<Grid>& factory, bool verbose = false) const
151 {
152 if (Dune::MPIHelper::instance().rank() == 0)
153 {
154 if (verbose)
155 std::cout << "Reading " << Grid::dimension << "d grid from vtk file " << reader_.filename() << "." << std::endl;
156
157 {
158 GridFormat::Dune::GridFactoryAdapter<Grid> adapter{ factory };
159 reader_.export_grid(adapter);
160 }
161 }
162
163 return std::unique_ptr<Grid>(factory.createGrid());
164 }
165
166 /*!
167 * \brief Read a grid from a vtk/vtu/vtp file, reading all cell and point data
168 * \note the factory will be needed outside to interpret the data via the factory's insertion indices
169 * \param factory the (empty) grid factory
170 * \param cellData the cell data arrays to be filled
171 * \param pointData the point data arrays to be filled
172 * \param verbose if the output should be verbose
173 */
174 template<class Grid>
175 std::unique_ptr<Grid> readGrid(Dune::GridFactory<Grid>& factory, Data& cellData, Data& pointData, bool verbose = false) const
176 {
177 auto grid = readGrid(factory, verbose);
178 if (Dune::MPIHelper::instance().rank() == 0)
179 {
180 for (const auto& [name, field_ptr] : cell_fields(reader_))
181 field_ptr->export_to(cellData[name]);
182
183 for (const auto& [name, field_ptr] : point_fields(reader_))
184 field_ptr->export_to(pointData[name]);
185 }
186 return std::unique_ptr<Grid>(std::move(grid));
187 }
188
189 private:
190 GridFormat::Reader reader_;
191 };
192
193 } // namespace Dumux
194
195 #else // DUMUX_HAVE_GRIDFORMAT
196
197 #include <dumux/io/xml/tinyxml2.h>
198 // fallback to simple vtk reader using tinyxml2
199 // this reader only support ascii files and limited VTK file formats
200
201 namespace Dumux::Detail::VTKReader {
202
203 /*!
204 * \brief Get the piece node an xml document
205 * \note Returns nullptr if the piece node wasn't found
206 * \param doc an xml document
207 * \param fileName a file name the doc was created from
208 */
209 491 const tinyxml2::XMLElement* getPieceNode(const tinyxml2::XMLDocument& doc, const std::string& fileName)
210 {
211 491 using namespace tinyxml2;
212
213 491 const XMLElement* pieceNode = doc.FirstChildElement("VTKFile");
214
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 491 times.
491 if (pieceNode == nullptr)
215 DUNE_THROW(Dune::IOError, "Couldn't get 'VTKFile' node in " << fileName << ".");
216
217 491 pieceNode = pieceNode->FirstChildElement("UnstructuredGrid");
218
2/2
✓ Branch 0 taken 124 times.
✓ Branch 1 taken 367 times.
491 if (pieceNode == nullptr)
219 248 pieceNode = doc.FirstChildElement("VTKFile")->FirstChildElement("PolyData");
220
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 122 times.
124 if (pieceNode == nullptr)
221 4 pieceNode = doc.FirstChildElement("VTKFile")->FirstChildElement("PUnstructuredGrid");
222
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 369 times.
369 if (pieceNode == nullptr)
223 pieceNode = doc.FirstChildElement("VTKFile")->FirstChildElement("PPolyData");
224 if (pieceNode == nullptr)
225 DUNE_THROW(Dune::IOError, "Couldn't get 'UnstructuredGrid', 'PUnstructuredGrid', 'PolyData', or 'PPolyData' node in " << fileName << ".");
226
227 982 return pieceNode->FirstChildElement("Piece");
228 }
229
230 /*!
231 * \brief get the vtk filename for the current processor
232 */
233 2 std::string getProcessPieceFileName(const std::string& pvtkFileName)
234 {
235 2 using namespace tinyxml2;
236
237 2 XMLDocument pDoc;
238
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
2 const auto eResult = pDoc.LoadFile(pvtkFileName.c_str());
239
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 if (eResult != XML_SUCCESS)
240 DUNE_THROW(Dune::IOError, "Couldn't open XML file " << pvtkFileName << ".");
241
242 // get the first piece node
243
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
2 const XMLElement* pieceNode = getPieceNode(pDoc, pvtkFileName);
244
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
2 const auto myrank = Dune::MPIHelper::instance().rank();
245
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
3 for (int rank = 0; rank < myrank; ++rank)
246 {
247 1 pieceNode = pieceNode->NextSiblingElement("Piece");
248
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if (pieceNode == nullptr)
249
0/34
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
1 DUNE_THROW(Dune::IOError, "Couldn't find 'Piece' node for rank "
250 << rank << " in " << pvtkFileName << ".");
251 }
252
253 2 const char *vtkFileName = pieceNode->Attribute("Source");
254
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 if (vtkFileName == nullptr)
255 DUNE_THROW(Dune::IOError, "Couldn't get 'Source' attribute of 'Piece' node no. " << myrank << " in " << pvtkFileName);
256
257
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
4 return vtkFileName;
258 2 }
259
260 } // end namespace Dumux::Detail::VTKReader
261
262 namespace Dumux {
263
264 /*!
265 * \ingroup InputOutput
266 * \brief A vtk file reader using tinyxml2 as xml backend
267 */
268 69 class VTKReader
269 {
270 public:
271 /*!
272 * \brief The data array types
273 */
274 enum class DataType { cellData, pointData };
275
276 //! the cell / point data type for point data read from a grid file
277 using Data = std::unordered_map<std::string, std::vector<double>>;
278
279 /*!
280 * \brief The constructor creates a tinyxml2::XMLDocument from file
281 */
282 69 explicit VTKReader(const std::string& fileName)
283 69 {
284 69 using namespace tinyxml2;
285
2/4
✓ Branch 1 taken 69 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 69 times.
✗ Branch 5 not taken.
138 const auto ext = std::filesystem::path(fileName).extension().string();
286 // If in parallel and the file to read is a parallel piece collection (pvtu/pvtp)
287 // read only the piece belonging to the own process. For this to work, the files
288 // have to have exactly the same amount of pieces than processes.
289
6/8
✓ Branch 1 taken 69 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 63 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 4 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
69 fileName_ = Dune::MPIHelper::instance().size() > 1 && ext[1] == 'p' ?
290
1/2
✓ Branch 1 taken 69 times.
✗ Branch 2 not taken.
138 Detail::VTKReader::getProcessPieceFileName(fileName) : fileName;
291
292
1/2
✓ Branch 1 taken 69 times.
✗ Branch 2 not taken.
69 const auto eResult = doc_.LoadFile(fileName_.c_str());
293
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 69 times.
69 if (eResult != tinyxml2::XML_SUCCESS)
294 DUNE_THROW(Dune::IOError, "Couldn't open XML file " << fileName_ << ".");
295
296 138 const XMLElement* pieceNode = getPieceNode_();
297
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 69 times.
69 if (pieceNode == nullptr)
298
0/30
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
69 DUNE_THROW(Dune::IOError, "Couldn't get 'Piece' node in " << fileName_ << ".");
299 69 }
300
301 /*!
302 * \brief Reviews data from the vtk file to check if there is a data array with a specified name
303 * \param name the name attribute of the data array to read
304 * \param type the data array type
305 */
306 273 bool hasData(const std::string& name, const DataType& type) const
307 {
308 273 using namespace tinyxml2;
309
310 546 const XMLElement* pieceNode = getPieceNode_();
311 273 const XMLElement* dataNode = getDataNode_(pieceNode, type);
312
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 273 times.
273 if (dataNode == nullptr)
313 return false;
314
315 273 const XMLElement* dataArray = findDataArray_(dataNode, name);
316
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 271 times.
273 if (dataArray == nullptr)
317 return false;
318
319 return true;
320 }
321
322 /*!
323 * \brief read data from the vtk file to a container, e.g. std::vector<double>
324 * \tparam Container a container type that has begin(), end(), push_back(), e.g. std::vector<>
325 * \param name the name attribute of the data array to read
326 * \param type the data array type
327 */
328 template<class Container>
329 139 Container readData(const std::string& name, const DataType& type) const
330 {
331 using namespace tinyxml2;
332
333 139 const XMLElement* pieceNode = getPieceNode_();
334 139 const XMLElement* dataNode = getDataNode_(pieceNode, type);
335
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 125 times.
139 if (dataNode == nullptr)
336 DUNE_THROW(Dune::IOError, "Couldn't get 'PointData' or 'CellData' node in " << fileName_ << ".");
337
338 139 const XMLElement* dataArray = findDataArray_(dataNode, name);
339
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 125 times.
139 if (dataArray == nullptr)
340 DUNE_THROW(Dune::IOError, "Couldn't find the data array " << name << ".");
341
342 139 return parseDataArray_<Container>(dataArray);
343 }
344
345 /*!
346 * \brief Read a grid from a vtk/vtu/vtp file, ignoring cell and point data
347 * \param verbose if the output should be verbose
348 */
349 template<class Grid>
350 std::unique_ptr<Grid> readGrid(bool verbose = false) const
351 {
352 static_assert(!Dune::Capabilities::isCartesian<Grid>::v, "Grid reader only supports unstructured grid implementations");
353
354 // make a grid factory
355 Dune::GridFactory<Grid> factory;
356
357 // only read on rank 0
358 if (Dune::MPIHelper::instance().rank() == 0)
359 {
360 if (verbose)
361 std::cout << "Reading " << Grid::dimension << "d grid from vtk file " << fileName_ << "." << std::endl;
362 readGrid_(factory, verbose);
363 }
364
365 return std::unique_ptr<Grid>(factory.createGrid());
366 }
367
368 /*!
369 * \brief Read a grid from a vtk/vtu/vtp file, ignoring cell and point data
370 * \note use this signature if the factory might be needed outside to interpret the data via the factory's insertion indices
371 * \param verbose if the output should be verbose
372 * \param factory the (empty) grid factory
373 */
374 template<class Grid>
375 std::unique_ptr<Grid> readGrid(Dune::GridFactory<Grid>& factory, bool verbose = false) const
376 {
377 static_assert(!Dune::Capabilities::isCartesian<Grid>::v, "Grid reader only supports unstructured grid implementations");
378
379 if (Dune::MPIHelper::instance().rank() == 0)
380 {
381 if (verbose)
382 std::cout << "Reading " << Grid::dimension << "d grid from vtk file " << fileName_ << "." << std::endl;
383 readGrid_(factory, verbose);
384 }
385
386 return std::unique_ptr<Grid>(factory.createGrid());
387 }
388
389 /*!
390 * \brief Read a grid from a vtk/vtu/vtp file, reading all cell and point data
391 * \note the factory will be needed outside to interpret the data via the factory's insertion indices
392 * \param factory the (empty) grid factory
393 * \param cellData the cell data arrays to be filled
394 * \param pointData the point data arrays to be filled
395 * \param verbose if the output should be verbose
396 */
397 template<class Grid>
398 13 std::unique_ptr<Grid> readGrid(Dune::GridFactory<Grid>& factory, Data& cellData, Data& pointData, bool verbose = false) const
399 {
400 static_assert(!Dune::Capabilities::isCartesian<Grid>::v, "Grid reader only supports unstructured grid implementations");
401
402
2/2
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 2 times.
13 if (Dune::MPIHelper::instance().rank() == 0)
403 {
404
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
11 if (verbose)
405 11 std::cout << "Reading " << Grid::dimension << "d grid from vtk file " << fileName_ << "." << std::endl;
406 11 readGrid_(factory, verbose);
407 11 readGridData_(cellData, pointData, verbose);
408 }
409
410 13 return std::unique_ptr<Grid>(factory.createGrid());
411 }
412
413 private:
414 /*!
415 * \brief Read a grid from a vtk/vtu/vtp file
416 * \param factory the (empty) grid factory
417 * \param verbose if the output should be verbose
418 */
419 template<class Grid>
420 11 void readGrid_(Dune::GridFactory<Grid>& factory, bool verbose = false) const
421 {
422 using namespace tinyxml2;
423
424 11 const XMLElement* pieceNode = getPieceNode_();
425 22 const XMLElement* pointsNode = pieceNode->FirstChildElement("Points")->FirstChildElement("DataArray");
426
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
11 if (pointsNode == nullptr)
427 DUNE_THROW(Dune::IOError, "Couldn't get data array of points in " << fileName_ << ".");
428
429 using Point3D = Dune::FieldVector<double, 3>;
430
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
11 std::vector<Point3D> points3D;
431
4/10
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 11 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 11 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
22 std::stringstream dataStream(pointsNode->GetText());
432
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
11 std::istream_iterator<Point3D> it(dataStream);
433
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
11 std::copy(it, std::istream_iterator<Point3D>(), std::back_inserter(points3D));
434
435 // adapt point dimensions if grid dimension is smaller than 3
436
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
11 auto points = adaptPointDimension_<Grid::dimensionworld>(std::move(points3D));
437
438
7/11
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 5 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 5 times.
✓ Branch 6 taken 6 times.
✓ Branch 7 taken 5 times.
✓ Branch 9 taken 6 times.
✗ Branch 10 not taken.
✗ Branch 2 not taken.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
17 if (verbose) std::cout << "Found " << points.size() << " vertices." << std::endl;
439
440 // insert vertices to the grid factory
441
2/2
✓ Branch 0 taken 18389 times.
✓ Branch 1 taken 11 times.
18400 for (auto&& point : points)
442
1/2
✓ Branch 1 taken 18389 times.
✗ Branch 2 not taken.
18389 factory.insertVertex(std::move(point));
443
444 11 const XMLElement* cellsNode = pieceNode->FirstChildElement("Cells");
445 11 const XMLElement* linesNode = pieceNode->FirstChildElement("Lines");
446
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 3 times.
11 if (cellsNode)
447 {
448
3/6
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 8 times.
✗ Branch 9 not taken.
8 const XMLElement* connectivityNode = findDataArray_(cellsNode, "connectivity");
449
3/6
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 8 times.
✗ Branch 9 not taken.
8 const XMLElement* offsetsNode = findDataArray_(cellsNode, "offsets");
450
2/6
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
8 const XMLElement* typesNode = findDataArray_(cellsNode, "types");
451
452
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
8 const auto connectivity = parseDataArray_<std::vector<unsigned int>>(connectivityNode);
453
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
8 const auto offsets = parseDataArray_<std::vector<unsigned int>>(offsetsNode);
454
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
8 const auto types = parseDataArray_<std::vector<unsigned int>>(typesNode);
455
456
4/8
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 8 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
16 if (verbose) std::cout << "Found " << offsets.size() << " elements." << std::endl;
457
458 unsigned int lastOffset = 0;
459
2/2
✓ Branch 0 taken 18827 times.
✓ Branch 1 taken 8 times.
18835 for (unsigned int i = 0; i < offsets.size(); ++i)
460 {
461
1/2
✓ Branch 1 taken 18827 times.
✗ Branch 2 not taken.
18827 const auto geomType = vtkToDuneGeomType_(types[i]);
462
1/2
✓ Branch 1 taken 18827 times.
✗ Branch 2 not taken.
18827 unsigned int offset = offsets[i];
463
1/2
✓ Branch 1 taken 18827 times.
✗ Branch 2 not taken.
18827 std::vector<unsigned int> corners; corners.resize(offset-lastOffset);
464
2/2
✓ Branch 0 taken 77461 times.
✓ Branch 1 taken 18827 times.
96288 for (unsigned int j = 0; j < offset-lastOffset; ++j)
465 77461 corners[Dune::VTK::renumber(geomType, j)] = connectivity[lastOffset+j];
466
1/2
✓ Branch 1 taken 18827 times.
✗ Branch 2 not taken.
18827 factory.insertElement(geomType, std::move(corners));
467
1/2
✓ Branch 0 taken 15231 times.
✗ Branch 1 not taken.
18827 lastOffset = offset;
468 }
469 13 }
470 // for poly data
471
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
3 else if (linesNode)
472 {
473 // sanity check
474 if (Grid::dimension != 1)
475
1/28
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
3 DUNE_THROW(Dune::IOError, "Grid expects dimension " << Grid::dimension
476 << " but " << fileName_ << " contains a 1D grid.");
477
478
3/6
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 3 times.
✗ Branch 9 not taken.
3 const XMLElement* connectivityNode = findDataArray_(linesNode, "connectivity");
479
2/6
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
3 const XMLElement* offsetsNode = findDataArray_(linesNode, "offsets");
480
481
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
3 const auto connectivity = parseDataArray_<std::vector<unsigned int>>(connectivityNode);
482
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
3 const auto offsets = parseDataArray_<std::vector<unsigned int>>(offsetsNode);
483
484
4/8
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
6 if (verbose) std::cout << "Found " << offsets.size() << " polylines." << std::endl;
485
486 unsigned int lastOffset = 0;
487
2/2
✓ Branch 0 taken 3479 times.
✓ Branch 1 taken 3 times.
3482 for (unsigned int i = 0; i < offsets.size(); ++i)
488 {
489 // a polyline can have many points in the VTK format
490 // split the line in segments with two points
491 3479 unsigned int offset = offsets[i];
492
2/2
✓ Branch 0 taken 3496 times.
✓ Branch 1 taken 3479 times.
6975 for (unsigned int j = 0; j < offset-lastOffset-1; ++j)
493
2/4
✓ Branch 1 taken 3496 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3496 times.
✗ Branch 5 not taken.
3496 factory.insertElement(Dune::GeometryTypes::line,
494
1/4
✓ Branch 1 taken 3496 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
6992 std::vector<unsigned int>({ connectivity[lastOffset+j], connectivity[lastOffset+j+1] }));
495 3479 lastOffset = offset;
496 }
497 4 }
498 else
499
2/40
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✓ Branch 39 taken 1 times.
✗ Branch 40 not taken.
✓ Branch 0 taken 6 times.
✗ Branch 3 not taken.
✗ Branch 6 not taken.
✗ Branch 9 not taken.
✗ Branch 12 not taken.
✗ Branch 15 not taken.
✗ Branch 18 not taken.
✗ Branch 21 not taken.
✗ Branch 24 not taken.
✗ Branch 27 not taken.
✗ Branch 30 not taken.
✗ Branch 33 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
7 DUNE_THROW(Dune::IOError, "No Cells or Lines element found in " << fileName_);
500 27 }
501
502 /*!
503 * \brief Read a grid data from a vtk/vtu/vtp file
504 * \param cellData the cell data arrays to be filled
505 * \param pointData the point data arrays to be filled
506 * \param verbose if the output should be verbose
507 */
508 11 void readGridData_(Data& cellData, Data& pointData, bool verbose = false) const
509 {
510 11 using namespace tinyxml2;
511
512 22 const XMLElement* pieceNode = getPieceNode_();
513 11 const XMLElement* cellDataNode = getDataNode_(pieceNode, DataType::cellData);
514
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
11 if (cellDataNode != nullptr)
515 {
516 11 const XMLElement* cellsNode = pieceNode->FirstChildElement("Cells");
517 11 const XMLElement* linesNode = pieceNode->FirstChildElement("Lines");
518
519
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 3 times.
11 if (cellsNode)
520 {
521 8 const XMLElement* dataArray = cellDataNode->FirstChildElement("DataArray");
522
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 8 times.
54 for (; dataArray != nullptr; dataArray = dataArray->NextSiblingElement("DataArray"))
523 {
524 23 const char *attributeText = dataArray->Attribute("Name");
525
526
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 23 times.
23 if (attributeText == nullptr)
527 DUNE_THROW(Dune::IOError, "Couldn't get Name attribute of a cell data array.");
528
529
2/6
✓ Branch 2 taken 23 times.
✗ Branch 3 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 23 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
46 cellData[std::string(attributeText)] = parseDataArray_<std::vector<double>>(dataArray);
530
531
1/2
✓ Branch 0 taken 23 times.
✗ Branch 1 not taken.
23 if (verbose)
532 23 std::cout << "Read cell data field " << attributeText << std::endl;
533 }
534 }
535 // for poly data
536
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
3 else if (linesNode)
537 {
538 // first parse all the cell data (each cell in this sense can be a polyline)
539 3 const XMLElement* dataArray = cellDataNode->FirstChildElement("DataArray");
540
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
3 if (dataArray)
541 {
542 3 Data polyLineCellData;
543
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 3 times.
13 for (; dataArray != nullptr; dataArray = dataArray->NextSiblingElement("DataArray"))
544 {
545 7 const char *attributeText = dataArray->Attribute("Name");
546
547
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
7 if (attributeText == nullptr)
548 DUNE_THROW(Dune::IOError, "Couldn't get Name attribute of a cell data array.");
549
550
3/8
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 7 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
14 polyLineCellData[std::string(attributeText)] = parseDataArray_<std::vector<double>>(dataArray);
551
552
1/2
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
7 if (verbose)
553
3/6
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.
7 std::cout << "Read cell data field " << attributeText << std::endl;
554 }
555
556 // a polyline can have many points in the VTK format
557 // we split the line in segments with two points
558 // so we also need to duplicate the cell data to fit the increased line number
559
2/4
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
3 const XMLElement* offsetsNode = findDataArray_(linesNode, "offsets");
560
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
3 const auto offsets = parseDataArray_<std::vector<unsigned int>>(offsetsNode);
561
562
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
3 if (offsets.size() != polyLineCellData.begin()->second.size())
563 DUNE_THROW(Dune::IOError, "Expected the same number of cell data entries (is "
564 << polyLineCellData.begin()->second.size()
565 << ") as polylines (" << offsets.size() << ")!");
566
567 // count the number of Dune cells to be able to resize the data vectors
568 unsigned int lastOffset = 0;
569 std::size_t numCells = 0;
570
2/2
✓ Branch 0 taken 3479 times.
✓ Branch 1 taken 3 times.
3482 for (unsigned int i = 0; i < offsets.size(); ++i)
571 {
572 3479 unsigned int offset = offsets[i];
573
2/2
✓ Branch 0 taken 3496 times.
✓ Branch 1 taken 3479 times.
6975 for (unsigned int j = 0; j < offset-lastOffset-1; ++j)
574 3496 ++numCells;
575 3479 lastOffset = offset;
576 }
577
578 // create the data arrays
579
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 3 times.
10 for (const auto& dArray : polyLineCellData)
580 {
581
3/8
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 7 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
7 cellData[dArray.first] = std::vector<double>(numCells);
582
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
7 auto& cd = cellData[dArray.first];
583 const auto& pd = dArray.second;
584
585 lastOffset = 0;
586 std::size_t cellIdx = 0;
587
2/2
✓ Branch 0 taken 8696 times.
✓ Branch 1 taken 7 times.
8703 for (unsigned int i = 0; i < offsets.size(); ++i)
588 {
589 8696 unsigned int offset = offsets[i];
590
2/2
✓ Branch 0 taken 8730 times.
✓ Branch 1 taken 8696 times.
17426 for (unsigned int j = 0; j < offset-lastOffset-1; ++j)
591 8730 cd[cellIdx++] = pd[i];
592 8696 lastOffset = offset;
593 }
594 }
595 3 }
596 }
597 else
598 DUNE_THROW(Dune::IOError, "No Cells or Lines element found in " << fileName_);
599 }
600
601 11 const XMLElement* pointDataNode = getDataNode_(pieceNode, DataType::pointData);
602
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 4 times.
11 if (pointDataNode != nullptr)
603 {
604 7 const XMLElement* dataArray = pointDataNode->FirstChildElement("DataArray");
605
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 7 times.
39 for (; dataArray != nullptr; dataArray = dataArray->NextSiblingElement("DataArray"))
606 {
607 16 const char *attributeText = dataArray->Attribute("Name");
608
609
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
16 if (attributeText == nullptr)
610 DUNE_THROW(Dune::IOError, "Couldn't get Name attribute of a point data array.");
611
612
2/6
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
32 pointData[std::string(attributeText)] = parseDataArray_<std::vector<double>>(dataArray);
613
614
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
16 if (verbose)
615 16 std::cout << "Read point data field " << attributeText << std::endl;
616 }
617 }
618 11 }
619
620 /*!
621 * \brief Get the piece node of the XMLDocument
622 * \note Returns nullptr if the piece node wasn't found
623 */
624 489 const tinyxml2::XMLElement* getPieceNode_() const
625
2/4
✓ Branch 3 taken 62 times.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
489 { return Detail::VTKReader::getPieceNode(doc_, fileName_); }
626
627 /*!
628 * \brief Get the piece node of the XMLDocument
629 * \param pieceNode the pieceNode of the vtk file
630 * \param type the vtk data type (cell data or point data)
631 * \note Returns nullptr if the data node wasn't found
632 */
633 420 const tinyxml2::XMLElement* getDataNode_(const tinyxml2::XMLElement* pieceNode, const DataType& type) const
634 {
635 420 using namespace tinyxml2;
636
637 420 const XMLElement* dataNode = nullptr;
638
2/2
✓ Branch 0 taken 132 times.
✓ Branch 1 taken 288 times.
420 if (type == DataType::pointData)
639 132 dataNode = pieceNode->FirstChildElement("PointData");
640
1/2
✓ Branch 0 taken 288 times.
✗ Branch 1 not taken.
288 else if (type == DataType::cellData)
641 288 dataNode = pieceNode->FirstChildElement("CellData");
642 else
643 DUNE_THROW(Dune::IOError, "Only cell and point data are supported.");
644
645 420 return dataNode;
646 }
647
648 /*!
649 * \brief Find a data array with a specific name
650 * \param dataNode a cell or point data node
651 * \param name the name of the data array to be found
652 * \note Returns nullptr if the data array wasn't found
653 */
654 431 const tinyxml2::XMLElement* findDataArray_(const tinyxml2::XMLElement* dataNode, const std::string& name) const
655 {
656 431 using namespace tinyxml2;
657
658 // loop over XML node siblings to find the correct data array
659 431 const XMLElement* dataArray = dataNode->FirstChildElement("DataArray");
660
2/2
✓ Branch 0 taken 2805 times.
✓ Branch 1 taken 2 times.
5183 for (; dataArray != nullptr; dataArray = dataArray->NextSiblingElement("DataArray"))
661 {
662 2805 const char *attributeText = dataArray->Attribute("Name");
663
664
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2805 times.
2805 if (attributeText == nullptr)
665 DUNE_THROW(Dune::IOError, "Couldn't get Name attribute of a data array.");
666
667
2/2
✓ Branch 0 taken 2376 times.
✓ Branch 1 taken 429 times.
2805 if (attributeText == name)
668 break;
669 }
670
671 431 return dataArray;
672 }
673
674 /*!
675 * \brief Parses the text of a data array into a container
676 * \tparam Container a container type that has begin(), end(), push_back(), e.g. std::vector<double>
677 * \param dataArray the data array node to be parsed
678 */
679 template<class Container>
680 297 Container parseDataArray_(const tinyxml2::XMLElement* dataArray) const
681 {
682
1/2
✓ Branch 3 taken 204 times.
✗ Branch 4 not taken.
297 std::stringstream dataStream(dataArray->GetText());
683
1/2
✓ Branch 1 taken 204 times.
✗ Branch 2 not taken.
594 return readStreamToContainer<Container>(dataStream);
684 297 }
685
686 /*!
687 * \brief Return the Dune::GeometryType for a given VTK geometry type
688 * \param vtkCellType the vtk cell type
689 */
690 18827 Dune::GeometryType vtkToDuneGeomType_(unsigned int vtkCellType) const
691 {
692
3/9
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 10639 times.
✓ Branch 3 taken 4990 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3198 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
18827 switch (vtkCellType)
693 {
694 case Dune::VTK::GeometryType::vertex: return Dune::GeometryTypes::vertex;
695 case Dune::VTK::GeometryType::line: return Dune::GeometryTypes::line;
696 10639 case Dune::VTK::GeometryType::triangle: return Dune::GeometryTypes::triangle;
697 4990 case Dune::VTK::GeometryType::quadrilateral: return Dune::GeometryTypes::quadrilateral;
698 case Dune::VTK::GeometryType::tetrahedron: return Dune::GeometryTypes::tetrahedron;
699 3198 case Dune::VTK::GeometryType::hexahedron: return Dune::GeometryTypes::hexahedron;
700 case Dune::VTK::GeometryType::prism: return Dune::GeometryTypes::prism;
701 case Dune::VTK::GeometryType::pyramid: return Dune::GeometryTypes::pyramid;
702 default: DUNE_THROW(Dune::NotImplemented, "VTK cell type " << vtkCellType);
703 }
704 }
705
706 template<int dim, std::enable_if_t<(dim < 3), int> = 0>
707 std::vector<Dune::FieldVector<double, dim>>
708 6 adaptPointDimension_(std::vector<Dune::FieldVector<double, 3>>&& points3D) const
709 {
710 6 std::vector<Dune::FieldVector<double, dim>> points(points3D.size());
711
2/2
✓ Branch 0 taken 10490 times.
✓ Branch 1 taken 6 times.
10496 for (std::size_t i = 0; i < points.size(); ++i)
712
2/2
✓ Branch 0 taken 20980 times.
✓ Branch 1 taken 10490 times.
31470 for (int j = 0; j < dim; ++j)
713 20980 points[i][j] = points3D[i][j];
714 6 return points;
715 }
716
717 template<int dim, std::enable_if_t<(dim == 3), int> = 0>
718 std::vector<Dune::FieldVector<double, dim>>
719
1/4
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
5 adaptPointDimension_(std::vector<Dune::FieldVector<double, 3>>&& points3D) const
720
1/4
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
5 { return std::move(points3D); }
721
722 std::string fileName_; //!< the vtk file name
723 tinyxml2::XMLDocument doc_; //!< the xml document created from file with name fileName_
724 };
725
726 } // end namespace Dumux
727
728 #endif // DUMUX_HAVE_GRIDFORMAT
729
730 #endif
731