24#ifndef DUMUX_IO_LOADSOLUTION_HH
25#define DUMUX_IO_LOADSOLUTION_HH
30#include <unordered_set>
31#include <unordered_map>
35#include <dune/common/exceptions.hh>
36#include <dune/common/indices.hh>
37#include <dune/grid/common/partitionset.hh>
53template <
class Container,
class EntityMapper,
int codim>
55:
public Dune::CommDataHandleIF< LoadSolutionDataHandle<Container, EntityMapper, codim>,
56 std::decay_t<decltype(std::declval<Container>()[0])> >
58 using FieldType = std::decay_t<decltype(std::declval<Container>()[0])>;
61 const EntityMapper& mapper)
63 , container_(container)
67 {
return cd == codim; }
72 template<
class EntityType>
73 size_t size (
const EntityType &e)
const
76 template<
class MessageBufferImp,
class EntityType>
77 void gather(MessageBufferImp& buff,
const EntityType& e)
const
79 const auto vIdx = mapper_.index(e);
80 buff.write(container_[vIdx]);
83 template<
class MessageBufferImp,
class EntityType>
84 void scatter(MessageBufferImp& buff,
const EntityType& e,
size_t n)
86 const auto vIdx = mapper_.index(e);
89 container_[vIdx] = tmp;
94 Container& container_;
101template <
class SolutionVector,
class PvNameFunc,
class Gr
idGeometry>
103 const std::string fileName,
104 PvNameFunc&& pvNameFunc,
105 const GridGeometry& gridGeometry,
110 using PrimaryVariables =
typename SolutionVector::block_type;
111 using Scalar =
typename PrimaryVariables::field_type;
112 constexpr auto dim = GridGeometry::GridView::dimension;
114 for (
size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
116 const auto pvName = pvNameFunc(pvIdx, 0);
117 auto vec = vtu.
readData<std::vector<Scalar>>(pvName, dataType);
122 for (
const auto& element : elements(gridGeometry.gridView(), Dune::Partitions::interior))
124 const auto eIdx = gridGeometry.elementMapper().index(element);
125 sol[eIdx][pvIdx] = vec[i++];
131 if (sol.size() != vec.size())
132 DUNE_THROW(Dune::InvalidStateException,
"Solution size (" << sol.size() <<
") does not match input size (" << vec.size() <<
")!");
134 for (std::size_t i = 0; i < sol.size(); ++i)
135 sol[i][pvIdx] = vec[i];
140 std::vector<bool> visited(gridGeometry.gridView().size(dim),
false);
141 for (
const auto& element : elements(gridGeometry.gridView(), Dune::Partitions::interior))
143 for (
int vIdxLocal = 0; vIdxLocal < element.subEntities(dim); ++vIdxLocal)
145 const auto vIdxGlobal = gridGeometry.vertexMapper().subIndex(element, vIdxLocal, dim);
146 if (!visited[vIdxGlobal])
148 sol[vIdxGlobal][pvIdx] = vec[i++];
149 visited[vIdxGlobal] =
true;
161template <
class SolutionVector,
class PvNameFunc,
class Gr
idGeometry>
163 const std::string fileName,
164 PvNameFunc&& pvNameFunc,
165 const GridGeometry& gridGeometry,
167->
typename std::enable_if_t<
decltype(
isValid(Detail::hasState())(sol[0]))::value, void>
169 VTKReader vtu(fileName);
172 const auto stateAtDof = vtu.readData<std::vector<int>>(
"phase presence", dataType);
175 std::unordered_set<int> states;
176 for (
size_t i = 0; i < stateAtDof.size(); ++i)
177 states.insert(stateAtDof[i]);
179 using PrimaryVariables =
typename SolutionVector::block_type;
180 using Scalar =
typename PrimaryVariables::field_type;
181 for (
size_t pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
183 std::unordered_map<int, std::vector<Scalar>> data;
184 for (
const auto& state : states)
185 data[state] = vtu.readData<std::vector<Scalar>>(pvNameFunc(pvIdx, state), dataType);
190 for (
const auto& element : elements(gridGeometry.gridView(), Dune::Partitions::interior))
192 const auto eIdx = gridGeometry.elementMapper().index(element);
193 const auto state = stateAtDof[i];
194 sol[eIdx][pvIdx] = data[state][i++];
195 sol[eIdx].setState(state);
201 constexpr int dim = GridGeometry::GridView::dimension;
202 std::vector<bool> visited(gridGeometry.gridView().size(dim),
false);
203 for (
const auto& element : elements(gridGeometry.gridView(), Dune::Partitions::interior))
205 for (
int vIdxLocal = 0; vIdxLocal < element.subEntities(dim); ++vIdxLocal)
207 const auto vIdxGlobal = gridGeometry.vertexMapper().subIndex(element, vIdxLocal, dim);
208 if (!visited[vIdxGlobal])
210 const auto state = stateAtDof[i];
211 sol[vIdxGlobal][pvIdx] = data[state][i++];
212 sol[vIdxGlobal].setState(state);
213 visited[vIdxGlobal] =
true;
226template<
class IOFields,
class PrimaryVariables,
class ModelTraits =
void,
class Flu
idSystem =
void,
class Sol
idSystem =
void>
228->
typename std::enable_if_t<
decltype(
isValid(
Detail::hasState())(PrimaryVariables(0)))::value, std::function<std::string(int,int)>>
230 return [paramGroup](
int pvIdx,
int state = 0)
232 static auto numStates = (1 << ModelTraits::numFluidPhases()) - 1;
233 const auto paramNameWithState =
"LoadSolution.PriVarNamesState" + std::to_string(state);
236 DUNE_THROW(Dune::NotImplemented,
"please provide LoadSolution.PriVarNamesState1..." << numStates
237 <<
" or remove LoadSolution.PriVarNames to use the model's default primary variable names");
241 const auto pvName = getParamFromGroup<std::vector<std::string>>(paramGroup, paramNameWithState);
242 return pvName[pvIdx];
245 return IOFields::template primaryVariableName<ModelTraits, FluidSystem, SolidSystem>(pvIdx, state);
254template<
class IOFields,
class PrimaryVariables,
class ModelTraits =
void,
class Flu
idSystem =
void,
class Sol
idSystem =
void>
256->
typename std::enable_if_t<!
decltype(
isValid(Detail::hasState())(PrimaryVariables(0)))::value, std::function<std::string(
int,
int)>>
260 const auto pvName = getParamFromGroup<std::vector<std::string>>(paramGroup,
"LoadSolution.PriVarNames");
261 return [n = std::move(pvName)](
int pvIdx,
int state = 0){
return n[pvIdx]; };
264 return [](
int pvIdx,
int state = 0){
return IOFields::template primaryVariableName<ModelTraits, FluidSystem, SolidSystem>(pvIdx, state); };
277template <
class SolutionVector,
class PvNameFunc,
class Gr
idGeometry>
279 const std::string& fileName,
280 PvNameFunc&& pvNameFunc,
281 const GridGeometry& gridGeometry)
283 const auto extension = fileName.substr(fileName.find_last_of(
".") + 1);
284 auto dataType = GridGeometry::discMethod == DiscretizationMethod::box ?
285 VTKReader::DataType::pointData : VTKReader::DataType::cellData;
287 if (extension ==
"vtu" || extension ==
"vtp")
289 if (GridGeometry::discMethod == DiscretizationMethod::staggered && extension ==
"vtp")
290 dataType = VTKReader::DataType::pointData;
294 else if (extension ==
"pvtu" || extension ==
"pvtp")
296 if (GridGeometry::discMethod == DiscretizationMethod::staggered)
297 DUNE_THROW(Dune::NotImplemented,
"reading staggered solution from a parallel vtk file");
302 DUNE_THROW(Dune::NotImplemented,
"loadSolution for file with extension " << extension);
305 if (gridGeometry.gridView().comm().size() > 1)
307 using GridView =
typename GridGeometry::GridView;
308 if (dataType == VTKReader::DataType::cellData)
311 dataHandle(sol, gridGeometry.elementMapper());
312 gridGeometry.gridView().communicate(dataHandle,
313 Dune::InteriorBorder_All_Interface,
314 Dune::ForwardCommunication);
319 dataHandle(sol, gridGeometry.vertexMapper());
320 gridGeometry.gridView().communicate(dataHandle,
321 Dune::InteriorBorder_All_Interface,
322 Dune::ForwardCommunication);
A helper function for class member function introspection.
Type traits to be used with matrix types.
Type traits to be used with vector types.
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
The available discretization methods in Dumux.
A vtk file reader using tinyxml2 as xml backend.
bool hasParamInGroup(const std::string ¶mGroup, const std::string ¶m)
Check whether a key exists in the parameter tree with a model group prefix.
Definition: parameters.hh:454
auto loadSolutionFromVtkFile(SolutionVector &sol, const std::string fileName, PvNameFunc &&pvNameFunc, const GridGeometry &gridGeometry, const VTKReader::DataType &dataType) -> typename std::enable_if_t<!decltype(isValid(Detail::hasState())(sol[0]))::value, void >
read from a vtk file into a solution vector with primary variables without state
Definition: loadsolution.hh:102
void loadSolution(SolutionVector &sol, const std::string &fileName, PvNameFunc &&pvNameFunc, const GridGeometry &gridGeometry)
load a solution vector from file
Definition: loadsolution.hh:278
auto createPVNameFunction(const std::string ¶mGroup="") -> typename std::enable_if_t< decltype(isValid(Detail::hasState())(PrimaryVariables(0)))::value, std::function< std::string(int, int)> >
helper function to determine the primary variable names of a model with privar state
Definition: loadsolution.hh:227
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
constexpr auto isValid(const Expression &t)
A function that creates a test functor to do class member introspection at compile time.
Definition: isvalid.hh:93
helper struct detecting if a PrimaryVariables object has a state() function
Definition: state.hh:32
a data handle to communicate the solution on ghosts and overlaps when reading from vtk file in parall...
Definition: loadsolution.hh:57
LoadSolutionDataHandle(Container &container, const EntityMapper &mapper)
Definition: loadsolution.hh:60
bool contains(int dim, int cd) const
Definition: loadsolution.hh:66
size_t size(const EntityType &e) const
Definition: loadsolution.hh:73
void scatter(MessageBufferImp &buff, const EntityType &e, size_t n)
Definition: loadsolution.hh:84
void gather(MessageBufferImp &buff, const EntityType &e) const
Definition: loadsolution.hh:77
bool fixedsize(int dim, int cd) const
Definition: loadsolution.hh:69
A vtk file reader using tinyxml2 as xml backend.
Definition: vtkreader.hh:48
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
DataType
The data array types.
Definition: vtkreader.hh:53