From fa05944bf98ea4ed8557859362420cf90767e0d4 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Thu, 26 Jul 2018 16:09:44 +0200
Subject: [PATCH] [io][vtk] read from parallel Vtk file into a solution vector
with state
---
dumux/io/loadsolution.hh | 60 ++++++++++++++++++++++++++++++++++++----
1 file changed, 54 insertions(+), 6 deletions(-)
diff --git a/dumux/io/loadsolution.hh b/dumux/io/loadsolution.hh
index c047c386d3..4b06ce1835 100644
--- a/dumux/io/loadsolution.hh
+++ b/dumux/io/loadsolution.hh
@@ -146,11 +146,6 @@ auto loadSolutionFromParallelVtkFile(const std::string fileName,
{
const auto pvName = pvNamesFunc(pvIdx);
auto vec = vtu.readData<std::vector<Scalar>>(pvName, dataType);
- if (vec.size() != sol.size() && fvGridGeometry.gridView().comm().size() == 1)
- {
- DUNE_THROW(Dune::IOError, "Size mismatch between solution vector and read data ("
- << sol.size() << " != " << vec.size() << ")");
- }
std::vector<bool> visited;
if (dataType == VTKReader::DataType::pointData)
@@ -239,7 +234,60 @@ auto loadSolutionFromParallelVtkFile(const std::string fileName,
const FVGridGeometry& fvGridGeometry)
-> typename std::enable_if_t<decltype(isValid(Detail::hasState())(sol[0]))::value, void>
{
- DUNE_THROW(Dune::NotImplemented, "reading solution with state from a parallel Vtk file");
+ VTKReader vtu(fileName);
+
+ // get states at dof locations
+ const auto stateAtDof = vtu.readData<std::vector<int>>("phase presence", dataType);
+
+ // determine which states are present
+ std::unordered_set<int> states;
+ for (size_t i = 0; i < sol.size(); ++i)
+ states.insert(stateAtDof[i]);
+
+ // declare a vector to determine which vertices have been visited already
+ std::vector<bool> visited;
+ using GridView = typename FVGridGeometry::GridView;
+ static constexpr auto dim = GridView::dimension;
+ if (dataType == VTKReader::DataType::pointData)
+ visited.resize(fvGridGeometry.gridView().size(dim));
+
+ using PrimaryVariables = typename SolutionVector::block_type;
+ using Scalar = typename PrimaryVariables::field_type;
+ for (size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
+ {
+ std::unordered_map<int, std::vector<Scalar>> switchedPvsSol;
+ for (const auto& state : states)
+ switchedPvsSol[state] = vtu.readData<std::vector<Scalar>>(pvNamesFunc(pvIdx, state), dataType);
+
+ std::fill(visited.begin(), visited.end(), false);
+ std::size_t i = 0;
+ for (const auto& element : elements(fvGridGeometry.gridView(), Dune::Partitions::interior))
+ {
+ if (dataType == VTKReader::DataType::cellData)
+ {
+ auto eIdx = fvGridGeometry.elementMapper().index(element);
+ auto state = stateAtDof[i];
+ sol[eIdx].setState(state);
+ sol[eIdx][pvIdx] = switchedPvsSol[state][i];
+ ++i;
+ }
+ else
+ {
+ for (int vIdxLocal = 0; vIdxLocal < element.subEntities(dim); ++vIdxLocal)
+ {
+ auto vIdxGlobal = fvGridGeometry.vertexMapper().subIndex(element, vIdxLocal, dim);
+ if (!visited[vIdxGlobal])
+ {
+ auto state = stateAtDof[i];
+ sol[vIdxGlobal].setState(state);
+ sol[vIdxGlobal][pvIdx] = switchedPvsSol[state][i];
+ visited[vIdxGlobal] = true;
+ ++i;
+ }
+ }
+ }
+ }
+ }
}
/*!
--
GitLab