From 910a16ef4489ef54ec12a69e0f0c16f1d0d7c978 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Wed, 11 Jul 2018 08:47:07 +0200
Subject: [PATCH] [io][vtk][xml] Implement vtu reader to separate xml interface
from loadsolution algorithm
---
dumux/io/loadsolution.hh | 197 +++++++++---------------
dumux/io/vtk/vtureader.hh | 114 ++++++++++++++
dumux/io/{tinyxml2 => xml}/tinyxml2.cpp | 0
dumux/io/{tinyxml2 => xml}/tinyxml2.h | 0
4 files changed, 190 insertions(+), 121 deletions(-)
create mode 100644 dumux/io/vtk/vtureader.hh
rename dumux/io/{tinyxml2 => xml}/tinyxml2.cpp (100%)
rename dumux/io/{tinyxml2 => xml}/tinyxml2.h (100%)
diff --git a/dumux/io/loadsolution.hh b/dumux/io/loadsolution.hh
index 25132ba57d..ac74e63f4f 100644
--- a/dumux/io/loadsolution.hh
+++ b/dumux/io/loadsolution.hh
@@ -21,17 +21,19 @@
* \ingroup InputOutput
* \brief read from a file into a solution vector
*/
-#ifndef DUMUX_LOADSOLUTION_HH
-#define DUMUX_LOADSOLUTION_HH
+#ifndef DUMUX_IO_LOADSOLUTION_HH
+#define DUMUX_IO_LOADSOLUTION_HH
#include <string>
#include <iostream>
+#include <vector>
+#include <unordered_set>
#include <dune/common/exceptions.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/typetraits/isvalid.hh>
-#include <dumux/io/tinyxml2/tinyxml2.h>
+#include <dumux/io/vtk/vtureader.hh>
namespace Dumux {
@@ -46,139 +48,78 @@ struct hasState
/*!
* \ingroup InputOutput
- * \brief helper class to read from a file into a solution vector
+ * \brief helper function to read from a file into a solution vector
*/
-class LoadSolutionHelper {
-public:
- template <class FVGridGeometry, class SolutionVector, class PvNamesFunc>
- static void loadFromVtkFile(const std::string fileName,
- const FVGridGeometry& fvGridGeometry,
- PvNamesFunc pvNamesFunc,
- SolutionVector& sol)
- {
- using namespace tinyxml2;
-
- XMLDocument xmlDoc;
- auto eResult = xmlDoc.LoadFile(fileName.c_str());
- if (eResult != XML_SUCCESS)
- DUNE_THROW(Dune::IOError, "Couldn't open XML file " << fileName << ".");
-
- XMLElement *pieceNode = xmlDoc.FirstChildElement("VTKFile")->FirstChildElement("UnstructuredGrid")->FirstChildElement("Piece");
- if (pieceNode == nullptr)
- DUNE_THROW(Dune::IOError, "Couldn't get Piece node in " << fileName << ".");
-
- XMLElement *dataNode = nullptr;
- if (FVGridGeometry::discMethod == DiscretizationMethod::box)
- dataNode = pieceNode->FirstChildElement("PointData");
- else
- dataNode = pieceNode->FirstChildElement("CellData");
+template <class SolutionVector, class PvNamesFunc>
+auto loadFromVtuFile(const std::string fileName,
+ const VTUReader::DataType& dataType,
+ PvNamesFunc&& pvNamesFunc,
+ SolutionVector& sol)
+-> typename std::enable_if_t<!decltype(isValid(hasState())(sol[0]))::value, void>
+{
+ VTUReader vtu(fileName);
+ using PrimaryVariables = typename SolutionVector::block_type;
+ using Scalar = typename PrimaryVariables::field_type;
- if (dataNode == nullptr)
- DUNE_THROW(Dune::IOError, "Couldn't get data node in " << fileName << ".");
+ for (size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
+ {
+ const auto pvName = pvNamesFunc(pvIdx);
+ const auto vec = vtu.readData<std::vector<Scalar>>(pvName, dataType);
+ if (vec.size() != sol.size())
+ DUNE_THROW(Dune::IOError, "Size mismatch between solution vector and read data (" << sol.size() << " != " << vec.size() << ")");
- setPrimaryVariables_(dataNode, fvGridGeometry, pvNamesFunc, sol);
+ for (std::size_t i = 0; i < sol.size(); ++i)
+ sol[i][pvIdx] = vec[i];
}
+}
-private:
-
- template <class Scalar, class FVGridGeometry>
- static std::vector<Scalar> extractDataToVector_(tinyxml2::XMLElement *xmlNode,
- const std::string& name,
- const FVGridGeometry& fvGridGeometry)
+/*!
+ * \ingroup InputOutput
+ * \brief helper function to read from a file into a solution vector
+ */
+template <class SolutionVector, class PvNamesFunc>
+auto loadFromVtuFile(const std::string fileName,
+ const VTUReader::DataType& dataType,
+ PvNamesFunc&& pvNamesFunc,
+ SolutionVector& sol)
+-> typename std::enable_if_t<decltype(isValid(hasState())(sol[0]))::value, void>
+{
+ VTUReader vtu(fileName);
+ const auto vec = vtu.readData<std::vector<int>>("phase presence", dataType);
+ std::unordered_set<int> states;
+ for (size_t i = 0; i < sol.size(); ++i)
{
- // loop over XML node siblings to find the correct data array
- tinyxml2::XMLElement *dataArray = xmlNode->FirstChildElement("DataArray");
- for (; dataArray != nullptr; dataArray = dataArray->NextSiblingElement("DataArray"))
- {
- const char *attributeText = dataArray->Attribute("Name");
-
- if (attributeText == nullptr)
- DUNE_THROW(Dune::IOError, "Couldn't get Name attribute.");
-
- if (std::string(attributeText) == name)
- break;
- }
- if (dataArray == nullptr)
- DUNE_THROW(Dune::IOError, "Couldn't find the data array " << name << ".");
-
- std::stringstream dataStream(dataArray->GetText());
- std::vector<Scalar> vec(fvGridGeometry.numDofs());
- for (auto& val : vec)
- {
- dataStream >> val;
- }
-
- return vec;
+ const int state = vec[i];
+ sol[i].setState(state);
+ states.insert(state);
}
- template<class SolutionVector, class FVGridGeometry, class PvNamesFunc>
- static auto setPrimaryVariables_(tinyxml2::XMLElement *node,
- const FVGridGeometry& fvGridGeometry,
- PvNamesFunc pvNamesFunc,
- SolutionVector& sol)
- -> typename std::enable_if_t<!decltype(isValid(hasState())(sol[0]))::value, void>
+ using PrimaryVariables = typename SolutionVector::block_type;
+ using Scalar = typename PrimaryVariables::field_type;
+ for (size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
{
- using PrimaryVariables = typename SolutionVector::block_type;
- using Scalar = typename PrimaryVariables::field_type;
-
- for (size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
+ if (pvNamesFunc(pvIdx, 1) == pvNamesFunc(pvIdx, 2))
{
- auto pvName = pvNamesFunc(pvIdx);
-
- auto vec = extractDataToVector_<Scalar>(node, pvName, fvGridGeometry);
-
+ const auto vec = vtu.readData<std::vector<Scalar>>(pvNamesFunc(pvIdx, 1), dataType);
for (size_t i = 0; i < sol.size(); ++i)
sol[i][pvIdx] = vec[i];
}
- }
-
- template<class SolutionVector, class FVGridGeometry, class PvNamesFunc>
- static auto setPrimaryVariables_(tinyxml2::XMLElement *node,
- const FVGridGeometry& fvGridGeometry,
- PvNamesFunc pvNamesFunc,
- SolutionVector& sol)
- -> typename std::enable_if_t<decltype(isValid(hasState())(sol[0]))::value, void>
- {
- auto vec = extractDataToVector_<int>(node, "phase presence", fvGridGeometry);
-
- int minState = std::numeric_limits<int>::max();
- int maxState = std::numeric_limits<int>::min();
- for (size_t i = 0; i < sol.size(); ++i)
+ else
{
- int state = vec[i];
- sol[i].setState(state);
- using std::min;
- minState = min(minState, state);
- using std::max;
- maxState = max(maxState, state);
- }
+ std::unordered_map<int, std::vector<Scalar>> switchedPvsSol;
+ for (const auto& state : states)
+ switchedPvsSol[state] = vtu.readData<std::vector<Scalar>>(pvNamesFunc(pvIdx, state), dataType);
- using PrimaryVariables = typename SolutionVector::block_type;
- using Scalar = typename PrimaryVariables::field_type;
- for (size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
- {
- if (pvNamesFunc(pvIdx, 1) == pvNamesFunc(pvIdx, 2))
- {
- auto vec = extractDataToVector_<Scalar>(node, pvNamesFunc(pvIdx, 1), fvGridGeometry);
-
- for (size_t i = 0; i < sol.size(); ++i)
- sol[i][pvIdx] = vec[i];
- }
- else
- {
- std::vector<std::vector<Scalar>> switchedPvsSol;
- for (int state = minState; state <= maxState; ++state)
- switchedPvsSol.push_back(extractDataToVector_<Scalar>(node,
- pvNamesFunc(pvIdx, state),
- fvGridGeometry));
-
- for (size_t i = 0; i < sol.size(); ++i)
- sol[i][pvIdx] = switchedPvsSol[sol[i].state() - minState][i];
- }
+ for (size_t i = 0; i < sol.size(); ++i)
+ sol[i][pvIdx] = switchedPvsSol[sol[i].state()][i];
}
}
-};
+}
+/*!
+ * \ingroup InputOutput
+ * \brief helper function to determine the primray variable names of a model with privar state
+ */
template<class ModelTraits, class FluidSystem>
std::string primaryVariableName(int pvIdx, int state)
{
@@ -188,6 +129,10 @@ std::string primaryVariableName(int pvIdx, int state)
return ModelTraits::template primaryVariableName<FluidSystem>(pvIdx, state);
}
+/*!
+ * \ingroup InputOutput
+ * \brief helper function to determine the primray variable names of a model
+ */
template<class ModelTraits>
std::string primaryVariableName(int pvIdx)
{
@@ -197,16 +142,26 @@ std::string primaryVariableName(int pvIdx)
return ModelTraits::primaryVariableName(pvIdx);
}
+
+/*!
+ * \ingroup InputOutput
+ * \brief load a solution vector from file
+ * \note Supports the following file extensions: *.vtu
+ */
template <class FVGridGeometry, class SolutionVector, class PvNamesFunc>
void loadSolution(const std::string& fileName,
const FVGridGeometry& fvGridGeometry,
- PvNamesFunc pvNamesFunc,
+ PvNamesFunc&& pvNamesFunc,
SolutionVector& sol)
{
- auto extension = fileName.substr(fileName.find_last_of(".") + 1);
+ const auto extension = fileName.substr(fileName.find_last_of(".") + 1);
if (extension == "vtu")
- LoadSolutionHelper::loadFromVtkFile(fileName, fvGridGeometry, pvNamesFunc, sol);
+ {
+ const auto dataType = FVGridGeometry::discMethod == DiscretizationMethod::box
+ ? VTUReader::DataType::pointData : VTUReader::DataType::cellData;
+ loadFromVtuFile(fileName, dataType, pvNamesFunc, sol);
+ }
else
DUNE_THROW(Dune::NotImplemented, "loadSolution for extension " << extension);
}
diff --git a/dumux/io/vtk/vtureader.hh b/dumux/io/vtk/vtureader.hh
new file mode 100644
index 0000000000..5ce9f7e6a6
--- /dev/null
+++ b/dumux/io/vtk/vtureader.hh
@@ -0,0 +1,114 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup InputOutput
+ * \brief A vtu reader using tinyxml2 as xml backend
+ */
+#ifndef DUMUX_IO_VTK_VTUREADER_HH
+#define DUMUX_IO_VTK_VTUREADER_HH
+
+#include <iostream>
+#include <iterator>
+#include <algorithm>
+
+#include <dune/common/exceptions.hh>
+#include <dumux/io/xml/tinyxml2.h>
+
+namespace Dumux {
+
+/*!
+ * \ingroup InputOutput
+ * \brief A vtu reader using tinyxml2 as xml backend
+ */
+class VTUReader
+{
+public:
+ /*!
+ * \brief The data array types
+ */
+ enum class DataType { cellData, pointData };
+
+ /*!
+ * \brief The contructor creates a tinyxml2::XMLDocument from file
+ */
+ VTUReader(const std::string& fileName)
+ : fileName_(fileName)
+ {
+ using namespace tinyxml2;
+ const auto eResult = doc_.LoadFile(fileName.c_str());
+ if (eResult != XML_SUCCESS)
+ DUNE_THROW(Dune::IOError, "Couldn't open XML file " << fileName << ".");
+ }
+
+ /*!
+ * \brief The contructor creates a tinyxml2::XMLDocument from file
+ * \tparam Container a container type that has begin(), end(), push_back(), e.g. std::vector<>
+ * \param name the name attribute of the data array to read
+ * \param type the data array type
+ */
+ template<class Container>
+ Container readData(const std::string& name, const DataType& type)
+ {
+ using namespace tinyxml2;
+ XMLElement *pieceNode = doc_.FirstChildElement("VTKFile")->FirstChildElement("UnstructuredGrid")->FirstChildElement("Piece");
+ if (pieceNode == nullptr)
+ DUNE_THROW(Dune::IOError, "Couldn't get Piece node in " << fileName_ << ".");
+
+ XMLElement *dataNode = nullptr;
+ if (type == DataType::pointData)
+ dataNode = pieceNode->FirstChildElement("PointData");
+ else if (type == DataType::cellData)
+ dataNode = pieceNode->FirstChildElement("CellData");
+ else
+ DUNE_THROW(Dune::IOError, "Only cell and point data are supported.");
+
+ if (dataNode == nullptr)
+ DUNE_THROW(Dune::IOError, "Couldn't get data node in " << fileName_ << ".");
+
+ // loop over XML node siblings to find the correct data array
+ XMLElement *dataArray = dataNode->FirstChildElement("DataArray");
+ for (; dataArray != nullptr; dataArray = dataArray->NextSiblingElement("DataArray"))
+ {
+ const char *attributeText = dataArray->Attribute("Name");
+
+ if (attributeText == nullptr)
+ DUNE_THROW(Dune::IOError, "Couldn't get Name attribute of a data array.");
+
+ if (attributeText == name)
+ break;
+ }
+ if (dataArray == nullptr)
+ DUNE_THROW(Dune::IOError, "Couldn't find the data array " << name << ".");
+
+ Container data;
+ std::stringstream dataStream(dataArray->GetText());
+ std::istream_iterator<typename Container::value_type> it(dataStream);
+ std::copy(it, std::istream_iterator<typename Container::value_type>(), std::back_inserter(data));
+ return data;
+ }
+
+private:
+ const std::string fileName_; //!< the vtu file name
+ tinyxml2::XMLDocument doc_; //!< the xml document created from file with name fileName_
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/io/tinyxml2/tinyxml2.cpp b/dumux/io/xml/tinyxml2.cpp
similarity index 100%
rename from dumux/io/tinyxml2/tinyxml2.cpp
rename to dumux/io/xml/tinyxml2.cpp
diff --git a/dumux/io/tinyxml2/tinyxml2.h b/dumux/io/xml/tinyxml2.h
similarity index 100%
rename from dumux/io/tinyxml2/tinyxml2.h
rename to dumux/io/xml/tinyxml2.h
--
GitLab