24#ifndef DUMUX_IO_VTK_OUTPUT_MODULE_HH
25#define DUMUX_IO_VTK_OUTPUT_MODULE_HH
31#include <dune/common/timer.hh>
32#include <dune/common/fvector.hh>
33#include <dune/common/typetraits.hh>
35#include <dune/geometry/type.hh>
36#include <dune/geometry/multilineargeometry.hh>
38#include <dune/grid/common/mcmgmapper.hh>
39#include <dune/grid/common/partitionset.hh>
40#include <dune/grid/io/file/vtk/vtkwriter.hh>
41#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
42#include <dune/grid/common/partitionset.hh>
59template<
class Gr
idGeometry>
62 using GridView =
typename GridGeometry::GridView;
63 using Field = Vtk::template Field<GridView>;
64 static constexpr int dim = GridView::dimension;
74 const std::string&
name,
76 Dune::VTK::DataMode dm = Dune::VTK::conforming,
84 const auto precisionString = getParamFromGroup<std::string>(
paramGroup,
"Vtk.Precision",
"Float32");
87 writer_ = std::make_shared<Dune::VTKWriter<GridView>>(
gridGeometry.gridView(), dm, coordPrecision);
88 sequenceWriter_ = std::make_unique<Dune::VTKSequenceWriter<GridView>>(writer_,
name);
95 {
return paramGroup_; }
106 template<
typename Vector>
108 const std::string&
name,
122 template<
typename Vector>
124 const std::string&
name,
129 const auto nComp = getNumberOfComponents_(v);
131 const auto numElemDofs =
gridGeometry().elementMapper().size();
132 const auto numVertexDofs =
gridGeometry().vertexMapper().size();
137 if(numElemDofs == numVertexDofs)
138 DUNE_THROW(Dune::InvalidStateException,
"Automatic deduction of FieldType failed. Please explicitly specify FieldType::element or FieldType::vertex.");
140 if(v.size() == numElemDofs)
142 else if(v.size() == numVertexDofs)
145 DUNE_THROW(Dune::RangeError,
"Size mismatch of added field!");
151 if(v.size() != numElemDofs)
152 DUNE_THROW(Dune::RangeError,
"Size mismatch of added field!");
155 if(v.size() != numVertexDofs)
156 DUNE_THROW(Dune::RangeError,
"Size mismatch of added field!");
161 fields_.emplace_back(gridGeometry_.gridView(), gridGeometry_.elementMapper(), v,
name, nComp, 0, dm_,
precision);
163 fields_.emplace_back(gridGeometry_.gridView(), gridGeometry_.vertexMapper(), v,
name, nComp, dim, dm_,
precision);
171 void write(
double time, Dune::VTK::OutputType type = Dune::VTK::ascii)
176 if (dm_ == Dune::VTK::conforming)
177 writeConforming_(time, type);
178 else if (dm_ == Dune::VTK::nonconforming)
179 writeNonConforming_(time, type);
181 DUNE_THROW(Dune::NotImplemented,
"Output for provided VTK data mode");
186 std::cout << Fmt::format(
"Writing output for problem \"{}\". Took {:.2g} seconds.\n", name_, timer.elapsed());
193 const std::string&
name()
const {
return name_; }
194 Dune::VTK::DataMode
dataMode()
const {
return dm_; }
195 Dumux::Vtk::Precision
precision()
const {
return precision_; }
197 Dune::VTKWriter<GridView>&
writer() {
return *writer_; }
200 const std::vector<Field>&
fields()
const {
return fields_; }
204 virtual void writeConforming_(
double time, Dune::VTK::OutputType type)
211 static bool addProcessRank = getParamFromGroup<bool>(this->
paramGroup(),
"Vtk.AddProcessRank");
212 std::vector<int> rank;
215 if (!fields_.empty() || addProcessRank)
217 const auto numCells = gridGeometry_.gridView().size(0);
222 rank.resize(numCells);
224 for (
const auto& element : elements(gridGeometry_.gridView(), Dune::Partitions::interior))
226 const auto eIdxGlobal = gridGeometry_.elementMapper().index(element);
227 rank[eIdxGlobal] = gridGeometry_.gridView().comm().rank();
237 this->
sequenceWriter().addCellData(Field(gridGeometry_.gridView(), gridGeometry_.elementMapper(), rank,
"process rank", 1, 0).get());
240 for (
auto&& field : fields_)
242 if (field.codim() == 0)
244 else if (field.codim() == dim)
247 DUNE_THROW(Dune::RangeError,
"Cannot add wrongly sized vtk scalar field!");
263 virtual void writeNonConforming_(
double time, Dune::VTK::OutputType type)
265 DUNE_THROW(Dune::NotImplemented,
"Non-conforming VTK output");
269 template<
class Vector>
270 std::size_t getNumberOfComponents_(
const Vector& v)
272 if constexpr (IsIndexable<decltype(std::declval<Vector>()[0])>::value)
278 const GridGeometry& gridGeometry_;
280 const std::string paramGroup_;
281 Dune::VTK::DataMode dm_;
283 Dumux::Vtk::Precision precision_;
285 std::shared_ptr<Dune::VTKWriter<GridView>> writer_;
286 std::unique_ptr<Dune::VTKSequenceWriter<GridView>> sequenceWriter_;
288 std::vector<Field> fields_;
303template<
class Gr
idVariables,
class SolutionVector>
307 using GridGeometry =
typename GridVariables::GridGeometry;
309 using VV =
typename GridVariables::VolumeVariables;
310 using Scalar =
typename GridVariables::Scalar;
312 using GridView =
typename GridGeometry::GridView;
315 dim = GridView::dimension,
316 dimWorld = GridView::dimensionworld
319 using Element =
typename GridView::template Codim<0>::Entity;
320 using VolVarsVector = Dune::FieldVector<Scalar, dimWorld>;
323 static constexpr int dofCodim = isBox ? dim : 0;
325 struct VolVarScalarDataInfo { std::function<Scalar(
const VV&)> get; std::string
name; Dumux::Vtk::Precision precision_; };
326 struct VolVarVectorDataInfo { std::function<VolVarsVector(
const VV&)> get; std::string
name; Dumux::Vtk::Precision precision_; };
327 using Field = Vtk::template Field<GridView>;
336 const SolutionVector&
sol,
337 const std::string&
name,
339 Dune::VTK::DataMode dm = Dune::VTK::conforming,
365 const std::string&
name)
367 volVarScalarDataInfo_.push_back(VolVarScalarDataInfo{f,
name, this->
precision()});
374 template<
class VVV = VolVarsVector,
typename std::enable_if_t<(VVV::dimension > 1),
int> = 0>
376 const std::string&
name)
378 volVarVectorDataInfo_.push_back(VolVarVectorDataInfo{f,
name, this->
precision()});
383 const auto&
problem()
const {
return gridVariables_.curGridVolVars().problem(); }
385 const GridGeometry&
gridGeometry()
const {
return gridVariables_.gridGeometry(); }
386 const SolutionVector&
sol()
const {
return sol_; }
397 void writeConforming_(
double time, Dune::VTK::OutputType type)
override
399 const Dune::VTK::DataMode dm = Dune::VTK::conforming;
406 std::vector<VelocityVector> velocity(velocityOutput_->numFluidPhases());
409 static bool addProcessRank = getParamFromGroup<bool>(this->
paramGroup(),
"Vtk.AddProcessRank");
410 std::vector<double> rank;
413 std::vector<std::vector<Scalar>> volVarScalarData;
414 std::vector<std::vector<VolVarsVector>> volVarVectorData;
417 if (!volVarScalarDataInfo_.empty()
418 || !volVarVectorDataInfo_.empty()
419 || !this->fields().empty()
420 || velocityOutput_->enableOutput()
423 const auto numCells =
gridGeometry().gridView().size(0);
424 const auto numDofs = numDofs_();
427 if (!volVarScalarDataInfo_.empty())
428 volVarScalarData.resize(volVarScalarDataInfo_.size(), std::vector<Scalar>(numDofs));
429 if (!volVarVectorDataInfo_.empty())
430 volVarVectorData.resize(volVarVectorDataInfo_.size(), std::vector<VolVarsVector>(numDofs));
432 if (velocityOutput_->enableOutput())
434 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
436 if(isBox && dim == 1)
437 velocity[phaseIdx].resize(numCells);
439 velocity[phaseIdx].resize(numDofs);
444 if (addProcessRank) rank.resize(numCells);
446 for (
const auto& element : elements(
gridGeometry().gridView(), Dune::Partitions::interior))
448 const auto eIdxGlobal =
gridGeometry().elementMapper().index(element);
451 auto elemVolVars =
localView(gridVariables_.curGridVolVars());
455 if (velocityOutput_->enableOutput())
457 fvGeometry.bind(element);
458 elemVolVars.bind(element, fvGeometry, sol_);
462 fvGeometry.bindElement(element);
463 elemVolVars.bindElement(element, fvGeometry, sol_);
466 if (!volVarScalarDataInfo_.empty()
467 || !volVarVectorDataInfo_.empty())
469 for (
auto&& scv : scvs(fvGeometry))
471 const auto dofIdxGlobal = scv.dofIndex();
472 const auto& volVars = elemVolVars[scv];
475 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
476 volVarScalarData[i][dofIdxGlobal] = volVarScalarDataInfo_[i].get(volVars);
479 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
480 volVarVectorData[i][dofIdxGlobal] = volVarVectorDataInfo_[i].get(volVars);
485 if (velocityOutput_->enableOutput())
487 auto elemFluxVarsCache =
localView(gridVariables_.gridFluxVarsCache());
488 elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
490 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
491 velocityOutput_->calculateVelocity(velocity[phaseIdx], element, fvGeometry, elemVolVars, elemFluxVarsCache, phaseIdx);
496 rank[eIdxGlobal] =
static_cast<double>(
gridGeometry().gridView().comm().rank());
506 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
508 volVarScalarDataInfo_[i].
name, 1, dim, dm, this->
precision()).get() );
509 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
511 volVarVectorDataInfo_[i].
name, dimWorld, dim, dm, this->
precision()).get() );
515 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
517 volVarScalarDataInfo_[i].
name, 1, 0,dm, this->
precision()).get() );
518 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
520 volVarVectorDataInfo_[i].
name, dimWorld, 0,dm, this->
precision()).get() );
524 if (velocityOutput_->enableOutput())
526 if (isBox && dim > 1)
528 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
530 "velocity_" + velocityOutput_->phaseName(phaseIdx) +
" (m/s)",
531 dimWorld, dim, dm, this->precision()).get() );
536 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
538 "velocity_" + velocityOutput_->phaseName(phaseIdx) +
" (m/s)",
539 dimWorld, 0, dm, this->precision()).get() );
548 for (
auto&& field : this->
fields())
550 if (field.codim() == 0)
552 else if (field.codim() == dim)
555 DUNE_THROW(Dune::RangeError,
"Cannot add wrongly sized vtk scalar field!");
571 void writeNonConforming_(
double time, Dune::VTK::OutputType type)
override
573 const Dune::VTK::DataMode dm = Dune::VTK::nonconforming;
576 DUNE_THROW(Dune::InvalidStateException,
"Non-conforming output makes no sense for cell-centered schemes!");
583 bool enableVelocityOutput = getParamFromGroup<bool>(this->
paramGroup(),
"Vtk.AddVelocity");
584 if (enableVelocityOutput ==
true && !velocityOutput_->enableOutput())
585 std::cerr <<
"Warning! Velocity output was enabled in the input file"
586 <<
" but no velocity output policy was set for the VTK output module:"
587 <<
" There will be no velocity output."
588 <<
" Use the addVelocityOutput member function of the VTK output module." << std::endl;
590 std::vector<VelocityVector> velocity(velocityOutput_->numFluidPhases());
593 static bool addProcessRank = getParamFromGroup<bool>(this->
paramGroup(),
"Vtk.AddProcessRank");
594 std::vector<double> rank;
597 using ScalarDataContainer = std::vector< std::vector<Scalar> >;
598 using VectorDataContainer = std::vector< std::vector<VolVarsVector> >;
599 std::vector< ScalarDataContainer > volVarScalarData;
600 std::vector< VectorDataContainer > volVarVectorData;
603 if (!volVarScalarDataInfo_.empty()
604 || !volVarVectorDataInfo_.empty()
605 || !this->fields().empty()
606 || velocityOutput_->enableOutput()
609 const auto numCells =
gridGeometry().gridView().size(0);
610 const auto numDofs = numDofs_();
613 if (!volVarScalarDataInfo_.empty())
614 volVarScalarData.resize(volVarScalarDataInfo_.size(), ScalarDataContainer(numCells));
615 if (!volVarVectorDataInfo_.empty())
616 volVarVectorData.resize(volVarVectorDataInfo_.size(), VectorDataContainer(numCells));
618 if (velocityOutput_->enableOutput())
620 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
622 if(isBox && dim == 1)
623 velocity[phaseIdx].resize(numCells);
625 velocity[phaseIdx].resize(numDofs);
630 if (addProcessRank) rank.resize(numCells);
632 for (
const auto& element : elements(
gridGeometry().gridView(), Dune::Partitions::interior))
634 const auto eIdxGlobal =
gridGeometry().elementMapper().index(element);
635 const auto numCorners = element.subEntities(dim);
638 auto elemVolVars =
localView(gridVariables_.curGridVolVars());
641 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
642 volVarScalarData[i][eIdxGlobal].resize(numCorners);
643 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
644 volVarVectorData[i][eIdxGlobal].resize(numCorners);
648 if (velocityOutput_->enableOutput())
650 fvGeometry.bind(element);
651 elemVolVars.bind(element, fvGeometry, sol_);
655 fvGeometry.bindElement(element);
656 elemVolVars.bindElement(element, fvGeometry, sol_);
659 if (!volVarScalarDataInfo_.empty()
660 || !volVarVectorDataInfo_.empty())
662 for (
auto&& scv : scvs(fvGeometry))
664 const auto& volVars = elemVolVars[scv];
667 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
668 volVarScalarData[i][scv.elementIndex()][scv.localDofIndex()] = volVarScalarDataInfo_[i].get(volVars);
671 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
672 volVarVectorData[i][scv.elementIndex()][scv.localDofIndex()] = volVarVectorDataInfo_[i].get(volVars);
677 if (velocityOutput_->enableOutput())
679 auto elemFluxVarsCache =
localView(gridVariables_.gridFluxVarsCache());
680 elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
682 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
683 velocityOutput_->calculateVelocity(velocity[phaseIdx], element, fvGeometry, elemVolVars, elemFluxVarsCache, phaseIdx);
688 rank[eIdxGlobal] =
static_cast<double>(
gridGeometry().gridView().comm().rank());
696 for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i)
698 volVarScalarDataInfo_[i].
name, 1, dim, dm, this->
precision()).get() );
700 for (std::size_t i = 0; i < volVarVectorDataInfo_.size(); ++i)
702 volVarVectorDataInfo_[i].
name, dimWorld, dim, dm, this->
precision()).get() );
705 if (velocityOutput_->enableOutput())
709 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
711 "velocity_" + velocityOutput_->phaseName(phaseIdx) +
" (m/s)",
712 dimWorld, dim, dm, this->precision()).get() );
716 for (
int phaseIdx = 0; phaseIdx < velocityOutput_->numFluidPhases(); ++phaseIdx)
718 "velocity_" + velocityOutput_->phaseName(phaseIdx) +
" (m/s)",
719 dimWorld, 0,dm, this->precision()).get());
727 for (
auto&& field : this->
fields())
729 if (field.codim() == 0)
731 else if (field.codim() == dim)
734 DUNE_THROW(Dune::RangeError,
"Cannot add wrongly sized vtk scalar field!");
750 std::size_t numDofs_()
const {
return dofCodim == dim ?
gridGeometry().vertexMapper().size() :
gridGeometry().elementMapper().size(); }
752 const GridVariables& gridVariables_;
753 const SolutionVector& sol_;
755 std::vector<VolVarScalarDataInfo> volVarScalarDataInfo_;
756 std::vector<VolVarVectorDataInfo> volVarVectorDataInfo_;
758 std::shared_ptr<VelocityOutput> velocityOutput_;
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
The available discretization methods in Dumux.
Formatting based on the fmt-library which implements std::format of C++20.
Dune style VTK functions.
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
Precision stringToPrecision(std::string_view precisionName)
Maps a string (e.g. from input) to a Dune precision type.
Definition: vtkprecision.hh:39
Velocity output for implicit (porous media) models.
Definition: io/velocityoutput.hh:41
std::vector< Dune::FieldVector< Scalar, dimWorld > > VelocityVector
Definition: io/velocityoutput.hh:50
A VTK output module to simplify writing dumux simulation data to VTK format.
Definition: io/vtkoutputmodule.hh:61
Dumux::Vtk::Precision precision() const
Definition: io/vtkoutputmodule.hh:195
const std::string & paramGroup() const
the parameter group for getting parameter from the parameter tree
Definition: io/vtkoutputmodule.hh:94
void write(double time, Dune::VTK::OutputType type=Dune::VTK::ascii)
Definition: io/vtkoutputmodule.hh:171
VtkOutputModuleBase(const GridGeometry &gridGeometry, const std::string &name, const std::string ¶mGroup="", Dune::VTK::DataMode dm=Dune::VTK::conforming, bool verbose=true)
Definition: io/vtkoutputmodule.hh:73
FieldType
export field type
Definition: io/vtkoutputmodule.hh:69
Dune::VTK::DataMode dataMode() const
Definition: io/vtkoutputmodule.hh:194
void addField(const Vector &v, const std::string &name, Dumux::Vtk::Precision precision, FieldType fieldType=FieldType::automatic)
Add a scalar or vector valued vtk field.
Definition: io/vtkoutputmodule.hh:123
Dune::VTKWriter< GridView > & writer()
Definition: io/vtkoutputmodule.hh:197
Dune::VTKSequenceWriter< GridView > & sequenceWriter()
Definition: io/vtkoutputmodule.hh:198
const std::string & name() const
Definition: io/vtkoutputmodule.hh:193
const std::vector< Field > & fields() const
Definition: io/vtkoutputmodule.hh:200
virtual ~VtkOutputModuleBase()=default
bool verbose() const
Definition: io/vtkoutputmodule.hh:192
void addField(const Vector &v, const std::string &name, FieldType fieldType=FieldType::automatic)
Add a scalar or vector valued vtk field.
Definition: io/vtkoutputmodule.hh:107
const GridGeometry & gridGeometry() const
Definition: io/vtkoutputmodule.hh:190
A VTK output module to simplify writing dumux simulation data to VTK format.
Definition: io/vtkoutputmodule.hh:305
const auto & problem() const
Definition: io/vtkoutputmodule.hh:383
const GridVariables & gridVariables() const
Definition: io/vtkoutputmodule.hh:384
void addVolumeVariable(std::function< VolVarsVector(const VolumeVariables &)> &&f, const std::string &name)
Definition: io/vtkoutputmodule.hh:375
VV VolumeVariables
export type of the volume variables for the outputfields
Definition: io/vtkoutputmodule.hh:333
void addVolumeVariable(std::function< Scalar(const VolumeVariables &)> &&f, const std::string &name)
Definition: io/vtkoutputmodule.hh:364
void addVelocityOutput(std::shared_ptr< VelocityOutputType > velocityOutput)
Add a velocity output policy.
Definition: io/vtkoutputmodule.hh:358
const VelocityOutput & velocityOutput() const
Definition: io/vtkoutputmodule.hh:392
const SolutionVector & sol() const
Definition: io/vtkoutputmodule.hh:386
const GridGeometry & gridGeometry() const
Definition: io/vtkoutputmodule.hh:385
VtkOutputModule(const GridVariables &gridVariables, const SolutionVector &sol, const std::string &name, const std::string ¶mGroup="", Dune::VTK::DataMode dm=Dune::VTK::conforming, bool verbose=true)
Definition: io/vtkoutputmodule.hh:335
const std::vector< VolVarScalarDataInfo > & volVarScalarDataInfo() const
Definition: io/vtkoutputmodule.hh:388
const std::vector< VolVarVectorDataInfo > & volVarVectorDataInfo() const
Definition: io/vtkoutputmodule.hh:389
Velocity output for porous media models.