24#ifndef DUMUX_IO_VTK_VTKREADER_HH
25#define DUMUX_IO_VTK_VTKREADER_HH
32#include <unordered_map>
34#include <dune/common/parallel/mpihelper.hh>
35#include <dune/common/exceptions.hh>
36#include <dune/grid/common/capabilities.hh>
37#include <dune/grid/io/file/vtk/common.hh>
38#include <dumux/io/xml/tinyxml2.h>
39#include <dune/grid/common/gridfactory.hh>
56 using Data = std::unordered_map<std::string, std::vector<double>>;
63 fileName_ = Dune::MPIHelper::getCollectiveCommunication().size() > 1 ?
64 getProcessFileName_(fileName) : fileName;
66 const auto eResult = doc_.LoadFile(fileName_.c_str());
67 if (eResult != tinyxml2::XML_SUCCESS)
68 DUNE_THROW(Dune::IOError,
"Couldn't open XML file " << fileName_ <<
".");
77 template<
class Container>
80 using namespace tinyxml2;
82 const XMLElement* pieceNode = getPieceNode_();
83 if (pieceNode ==
nullptr)
84 DUNE_THROW(Dune::IOError,
"Couldn't get 'Piece' node in " << fileName_ <<
".");
86 const XMLElement* dataNode = getDataNode_(pieceNode, type);
87 if (dataNode ==
nullptr)
88 DUNE_THROW(Dune::IOError,
"Couldn't get 'PointData' or 'CellData' node in " << fileName_ <<
".");
90 const XMLElement* dataArray = findDataArray_(dataNode, name);
91 if (dataArray ==
nullptr)
92 DUNE_THROW(Dune::IOError,
"Couldn't find the data array " << name <<
".");
94 return parseDataArray_<Container>(dataArray);
102 std::unique_ptr<Grid>
readGrid(
bool verbose =
false)
const
104 static_assert(!Dune::Capabilities::isCartesian<Grid>::v,
"Grid reader only supports unstructured grid implementations");
106 if (verbose) std::cout <<
"Reading " << Grid::dimension <<
"d grid from vtk file " << fileName_ <<
"." << std::endl;
109 Dune::GridFactory<Grid> factory;
111 readGrid_(factory, verbose);
113 return std::unique_ptr<Grid>(factory.createGrid());
123 std::unique_ptr<Grid>
readGrid(Dune::GridFactory<Grid>& factory,
bool verbose =
false)
const
125 static_assert(!Dune::Capabilities::isCartesian<Grid>::v,
"Grid reader only supports unstructured grid implementations");
127 if (verbose) std::cout <<
"Reading " << Grid::dimension <<
"d grid from vtk file " << fileName_ <<
"." << std::endl;
129 readGrid_(factory, verbose);
131 return std::unique_ptr<Grid>(factory.createGrid());
143 std::unique_ptr<Grid>
readGrid(Dune::GridFactory<Grid>& factory,
Data& cellData,
Data& pointData,
bool verbose =
false)
const
145 static_assert(!Dune::Capabilities::isCartesian<Grid>::v,
"Grid reader only supports unstructured grid implementations");
147 if (verbose) std::cout <<
"Reading " << Grid::dimension <<
"d grid from vtk file " << fileName_ <<
"." << std::endl;
149 readGrid_(factory, verbose);
150 readGridData_(cellData, pointData, verbose);
152 return std::unique_ptr<Grid>(factory.createGrid());
159 std::string getProcessFileName_(
const std::string& pvtkFileName)
161 using namespace tinyxml2;
164 const auto eResult = pDoc.LoadFile(pvtkFileName.c_str());
165 if (eResult != XML_SUCCESS)
166 DUNE_THROW(Dune::IOError,
"Couldn't open XML file " << pvtkFileName <<
".");
169 const XMLElement* pieceNode = getPieceNode_(pDoc, pvtkFileName);
170 if (pieceNode ==
nullptr)
171 DUNE_THROW(Dune::IOError,
"Couldn't get 'Piece' node in " << pvtkFileName <<
".");
173 const auto myrank = Dune::MPIHelper::getCollectiveCommunication().rank();
174 for (
int rank = 0; rank < myrank; ++rank)
176 pieceNode = pieceNode->NextSiblingElement(
"Piece");
177 if (pieceNode ==
nullptr)
178 DUNE_THROW(Dune::IOError,
"Couldn't find 'Piece' node for rank "
179 << rank <<
" in " << pvtkFileName <<
".");
182 const char *vtkFileName = pieceNode->Attribute(
"Source");
183 if (vtkFileName ==
nullptr)
184 DUNE_THROW(Dune::IOError,
"Couldn't get 'Source' attribute of 'Piece' node no. " << myrank <<
" in " << pvtkFileName);
195 void readGrid_(Dune::GridFactory<Grid>& factory,
bool verbose =
false)
const
197 using namespace tinyxml2;
199 const XMLElement* pieceNode = getPieceNode_();
200 if (pieceNode ==
nullptr)
201 DUNE_THROW(Dune::IOError,
"Couldn't get 'Piece' node in " << fileName_ <<
".");
203 const XMLElement* pointsNode = pieceNode->FirstChildElement(
"Points")->FirstChildElement(
"DataArray");
204 if (pointsNode ==
nullptr)
205 DUNE_THROW(Dune::IOError,
"Couldn't get data array of points in " << fileName_ <<
".");
207 using Point3D = Dune::FieldVector<double, 3>;
208 std::vector<Point3D> points3D;
209 std::stringstream dataStream(pointsNode->GetText());
210 std::istream_iterator<Point3D> it(dataStream);
211 std::copy(it, std::istream_iterator<Point3D>(), std::back_inserter(points3D));
214 auto points = adaptPointDimension_<Grid::dimensionworld>(std::move(points3D));
216 if (verbose) std::cout <<
"Found " << points.size() <<
" vertices." << std::endl;
219 for (
auto&& point : points)
220 factory.insertVertex(std::move(point));
222 const XMLElement* cellsNode = pieceNode->FirstChildElement(
"Cells");
223 const XMLElement* linesNode = pieceNode->FirstChildElement(
"Lines");
226 const XMLElement* connectivityNode = findDataArray_(cellsNode,
"connectivity");
227 const XMLElement* offsetsNode = findDataArray_(cellsNode,
"offsets");
228 const XMLElement* typesNode = findDataArray_(cellsNode,
"types");
230 const auto connectivity = parseDataArray_<std::vector<unsigned int>>(connectivityNode);
231 const auto offsets = parseDataArray_<std::vector<unsigned int>>(offsetsNode);
232 const auto types = parseDataArray_<std::vector<unsigned int>>(typesNode);
234 if (verbose) std::cout <<
"Found " << offsets.size() <<
" element." << std::endl;
236 unsigned int lastOffset = 0;
237 for (
unsigned int i = 0; i < offsets.size(); ++i)
239 const auto geomType = vtkToDuneGeomType_(types[i]);
240 unsigned int offset = offsets[i];
241 std::vector<unsigned int> corners; corners.resize(offset-lastOffset);
242 for (
unsigned int j = 0; j < offset-lastOffset; ++j)
243 corners[Dune::VTK::renumber(geomType, j)] = connectivity[lastOffset+j];
244 factory.insertElement(geomType, std::move(corners));
252 if (Grid::dimension != 1)
253 DUNE_THROW(Dune::IOError,
"Grid expects dimension " << Grid::dimension
254 <<
" but " << fileName_ <<
" contains a 1D grid.");
256 const XMLElement* connectivityNode = findDataArray_(linesNode,
"connectivity");
257 const XMLElement* offsetsNode = findDataArray_(linesNode,
"offsets");
259 const auto connectivity = parseDataArray_<std::vector<unsigned int>>(connectivityNode);
260 const auto offsets = parseDataArray_<std::vector<unsigned int>>(offsetsNode);
262 if (verbose) std::cout <<
"Found " << offsets.size() <<
" polylines." << std::endl;
264 unsigned int lastOffset = 0;
265 for (
unsigned int i = 0; i < offsets.size(); ++i)
269 unsigned int offset = offsets[i];
270 for (
unsigned int j = 0; j < offset-lastOffset-1; ++j)
271 factory.insertElement(Dune::GeometryTypes::line,
272 std::vector<unsigned int>({ connectivity[lastOffset+j], connectivity[lastOffset+j+1] }));
277 DUNE_THROW(Dune::IOError,
"No Cells or Lines element found in " << fileName_);
286 void readGridData_(Data& cellData, Data& pointData,
bool verbose =
false)
const
288 using namespace tinyxml2;
290 const XMLElement* pieceNode = getPieceNode_();
291 if (pieceNode ==
nullptr)
292 DUNE_THROW(Dune::IOError,
"Couldn't get 'Piece' node in " << fileName_ <<
".");
294 const XMLElement* cellDataNode = getDataNode_(pieceNode, DataType::cellData);
295 if (cellDataNode !=
nullptr)
297 const XMLElement* cellsNode = pieceNode->FirstChildElement(
"Cells");
298 const XMLElement* linesNode = pieceNode->FirstChildElement(
"Lines");
302 const XMLElement* dataArray = cellDataNode->FirstChildElement(
"DataArray");
303 for (; dataArray !=
nullptr; dataArray = dataArray->NextSiblingElement(
"DataArray"))
305 const char *attributeText = dataArray->Attribute(
"Name");
307 if (attributeText ==
nullptr)
308 DUNE_THROW(Dune::IOError,
"Couldn't get Name attribute of a cell data array.");
310 cellData[std::string(attributeText)] = parseDataArray_<std::vector<double>>(dataArray);
317 const XMLElement* dataArray = cellDataNode->FirstChildElement(
"DataArray");
320 Data polyLineCellData;
321 for (; dataArray !=
nullptr; dataArray = dataArray->NextSiblingElement(
"DataArray"))
323 const char *attributeText = dataArray->Attribute(
"Name");
325 if (attributeText ==
nullptr)
326 DUNE_THROW(Dune::IOError,
"Couldn't get Name attribute of a cell data array.");
328 polyLineCellData[std::string(attributeText)] = parseDataArray_<std::vector<double>>(dataArray);
334 const XMLElement* offsetsNode = findDataArray_(linesNode,
"offsets");
335 const auto offsets = parseDataArray_<std::vector<unsigned int>>(offsetsNode);
337 if (offsets.size() != polyLineCellData.begin()->second.size())
338 DUNE_THROW(Dune::IOError,
"Expected the same number of cell data entries (is "
339 << polyLineCellData.begin()->second.size()
340 <<
") as polylines (" << offsets.size() <<
")!");
343 unsigned int lastOffset = 0;
344 std::size_t numCells = 0;
345 for (
unsigned int i = 0; i < offsets.size(); ++i)
347 unsigned int offset = offsets[i];
348 for (
unsigned int j = 0; j < offset-lastOffset-1; ++j)
354 for (
const auto& dArray : polyLineCellData)
356 cellData[dArray.first] = std::vector<double>(numCells);
357 auto& cd = cellData[dArray.first];
358 const auto& pd = dArray.second;
361 std::size_t cellIdx = 0;
362 for (
unsigned int i = 0; i < offsets.size(); ++i)
364 unsigned int offset = offsets[i];
365 for (
unsigned int j = 0; j < offset-lastOffset-1; ++j)
366 cd[cellIdx++] = pd[i];
373 DUNE_THROW(Dune::IOError,
"No Cells or Lines element found in " << fileName_);
376 const XMLElement* pointDataNode = getDataNode_(pieceNode, DataType::pointData);
377 if (pointDataNode !=
nullptr)
379 const XMLElement* dataArray = pointDataNode->FirstChildElement(
"DataArray");
380 for (; dataArray !=
nullptr; dataArray = dataArray->NextSiblingElement(
"DataArray"))
382 const char *attributeText = dataArray->Attribute(
"Name");
384 if (attributeText ==
nullptr)
385 DUNE_THROW(Dune::IOError,
"Couldn't get Name attribute of a point data array.");
387 pointData[std::string(attributeText)] = parseDataArray_<std::vector<double>>(dataArray);
396 const tinyxml2::XMLElement* getPieceNode_()
const
397 {
return getPieceNode_(doc_, fileName_); }
405 const tinyxml2::XMLElement* getPieceNode_(
const tinyxml2::XMLDocument& doc,
const std::string& fileName)
const
407 using namespace tinyxml2;
409 const XMLElement* pieceNode = doc.FirstChildElement(
"VTKFile");
410 if (pieceNode ==
nullptr)
411 DUNE_THROW(Dune::IOError,
"Couldn't get 'VTKFile' node in " << fileName <<
".");
413 pieceNode = pieceNode->FirstChildElement(
"UnstructuredGrid");
414 if (pieceNode ==
nullptr)
415 pieceNode = doc.FirstChildElement(
"VTKFile")->FirstChildElement(
"PolyData");
416 if (pieceNode ==
nullptr)
417 pieceNode = doc.FirstChildElement(
"VTKFile")->FirstChildElement(
"PUnstructuredGrid");
418 if (pieceNode ==
nullptr)
419 pieceNode = doc.FirstChildElement(
"VTKFile")->FirstChildElement(
"PPolyData");
420 if (pieceNode ==
nullptr)
421 DUNE_THROW(Dune::IOError,
"Couldn't get 'UnstructuredGrid', 'PUnstructuredGrid', 'PolyData', or 'PPolyData' node in " << fileName <<
".");
423 return pieceNode->FirstChildElement(
"Piece");
432 const tinyxml2::XMLElement* getDataNode_(
const tinyxml2::XMLElement* pieceNode,
const DataType& type)
const
434 using namespace tinyxml2;
436 const XMLElement* dataNode =
nullptr;
437 if (type == DataType::pointData)
438 dataNode = pieceNode->FirstChildElement(
"PointData");
439 else if (type == DataType::cellData)
440 dataNode = pieceNode->FirstChildElement(
"CellData");
442 DUNE_THROW(Dune::IOError,
"Only cell and point data are supported.");
453 const tinyxml2::XMLElement* findDataArray_(
const tinyxml2::XMLElement* dataNode,
const std::string& name)
const
455 using namespace tinyxml2;
458 const XMLElement* dataArray = dataNode->FirstChildElement(
"DataArray");
459 for (; dataArray !=
nullptr; dataArray = dataArray->NextSiblingElement(
"DataArray"))
461 const char *attributeText = dataArray->Attribute(
"Name");
463 if (attributeText ==
nullptr)
464 DUNE_THROW(Dune::IOError,
"Couldn't get Name attribute of a data array.");
466 if (attributeText == name)
478 template<
class Container>
479 Container parseDataArray_(
const tinyxml2::XMLElement* dataArray)
const
482 std::stringstream dataStream(dataArray->GetText());
483 std::istream_iterator<typename Container::value_type> it(dataStream);
484 std::copy(it, std::istream_iterator<typename Container::value_type>(), std::back_inserter(data));
492 Dune::GeometryType vtkToDuneGeomType_(
unsigned int vtkCellType)
const
496 case Dune::VTK::GeometryType::vertex:
return Dune::GeometryTypes::vertex;
497 case Dune::VTK::GeometryType::line:
return Dune::GeometryTypes::line;
498 case Dune::VTK::GeometryType::triangle:
return Dune::GeometryTypes::triangle;
499 case Dune::VTK::GeometryType::quadrilateral:
return Dune::GeometryTypes::quadrilateral;
500 case Dune::VTK::GeometryType::hexahedron:
return Dune::GeometryTypes::hexahedron;
501 case Dune::VTK::GeometryType::prism:
return Dune::GeometryTypes::prism;
502 case Dune::VTK::GeometryType::pyramid:
return Dune::GeometryTypes::pyramid;
503 default: DUNE_THROW(Dune::NotImplemented,
"VTK cell type " << vtkCellType);
507 template<
int dim, std::enable_if_t<(dim < 3),
int> = 0>
508 std::vector<Dune::FieldVector<
double, dim>>
509 adaptPo
intDimension_(std::vector<Dune::FieldVector<
double, 3>>&& po
ints3D) const
511 std::vector<Dune::FieldVector<
double, dim>> po
ints(po
ints3D.size());
512 for (std::
size_t i = 0; i < po
ints.size(); ++i)
513 for (
int j = 0; j < dim; ++j)
514 po
ints[i][j] = po
ints3D[i][j];
518 template<
int dim, std::enable_if_t<(dim == 3),
int> = 0>
519 std::vector<Dune::FieldVector<double, dim>>
520 adaptPointDimension_(std::vector<Dune::FieldVector<double, 3>>&& points3D)
const
523 std::string fileName_;
524 tinyxml2::XMLDocument doc_;
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
A vtk file reader using tinyxml2 as xml backend.
Definition: vtkreader.hh:48
std::unique_ptr< Grid > readGrid(Dune::GridFactory< Grid > &factory, Data &cellData, Data &pointData, bool verbose=false) const
Read a grid from a vtk/vtu/vtp file, reading all cell and point data.
Definition: vtkreader.hh:143
Container readData(const std::string &name, const DataType &type) const
read data from the vtk file to a container, e.g. std::vector<double>
Definition: vtkreader.hh:78
VTKReader(const std::string &fileName)
The contructor creates a tinyxml2::XMLDocument from file.
Definition: vtkreader.hh:61
std::unique_ptr< Grid > readGrid(Dune::GridFactory< Grid > &factory, bool verbose=false) const
Read a grid from a vtk/vtu/vtp file, ignoring cell and point data.
Definition: vtkreader.hh:123
std::unique_ptr< Grid > readGrid(bool verbose=false) const
Read a grid from a vtk/vtu/vtp file, ignoring cell and point data.
Definition: vtkreader.hh:102
DataType
The data array types.
Definition: vtkreader.hh:53
std::unordered_map< std::string, std::vector< double > > Data
the cell / point data type for point data read from a grid file
Definition: vtkreader.hh:56