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