diff --git a/dumux/freeflow/staggered/vtkoutputfields.hh b/dumux/freeflow/staggered/vtkoutputfields.hh
index e4273b45348d1fe2dd518cc37adb95daea289410..8e848827ad2ca35d4441da2ddb3e489268ac58b3 100644
--- a/dumux/freeflow/staggered/vtkoutputfields.hh
+++ b/dumux/freeflow/staggered/vtkoutputfields.hh
@@ -23,7 +23,8 @@
#ifndef DUMUX_NAVIER_STOKES_VTK_OUTPUT_FIELDS_HH
#define DUMUX_NAVIER_STOKES_VTK_OUTPUT_FIELDS_HH
-#include <dumux/common/basicproperties.hh>
+#include <dumux/common/properties.hh>
+#include <dune/common/fvector.hh>
namespace Dumux
{
@@ -36,13 +37,35 @@ template<class TypeTag>
class NavierStokesVtkOutputFields
{
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
+
+ using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
+
public:
template <class VtkOutputModule>
static void init(VtkOutputModule& vtk)
{
vtk.addVolumeVariable([](const VolumeVariables& v){ return v.pressure(); }, "p");
+ const bool writeFaceVars_ = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Vtk.WriteFaceData", false);
+ if(writeFaceVars_)
+ {
+ vtk.addFaceVariable([](const FaceVariables& f){ return f.velocitySelf(); }, "scalarFaceVelocity");
+
+ auto faceVelocityVector = [](const SubControlVolumeFace& scvf, const FaceVariables& f)
+ {
+ GlobalPosition velocity(0.0);
+ velocity[scvf.directionIndex()] = f.velocitySelf();
+ return velocity;
+ };
+
+ vtk.addFaceVariable(faceVelocityVector, "vectorFaceVelocity");
+ }
+
if(GET_PROP_VALUE(TypeTag, EnableEnergyBalance))
vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.temperature(); },"temperature");
}
diff --git a/dumux/io/staggeredvtkoutputmodule.hh b/dumux/io/staggeredvtkoutputmodule.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a531c138196dcaa90f640b279dfa6799b5315372
--- /dev/null
+++ b/dumux/io/staggeredvtkoutputmodule.hh
@@ -0,0 +1,232 @@
+// -*- 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
+ * \brief A VTK output module to simplify writing dumux simulation data to VTK format. Specialization for staggered grids with dofs on faces.
+ */
+#ifndef STAGGERED_VTK_OUTPUT_MODULE_HH
+#define STAGGERED_VTK_OUTPUT_MODULE_HH
+
+#include <dune/common/fvector.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/pointcloudvtkwriter.hh>
+#include <dumux/io/vtksequencewriter.hh>
+
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup InputOutput
+ * \brief A VTK output module to simplify writing dumux simulation data to VTK format
+ * Specialization for staggered grids with dofs on faces.
+ */
+template<typename TypeTag>
+class StaggeredVtkOutputModule : public VtkOutputModule<TypeTag>
+{
+ friend class VtkOutputModule<TypeTag>;
+ using ParentType = VtkOutputModule<TypeTag>;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+
+
+ enum { dim = GridView::dimension };
+ enum { dimWorld = GridView::dimensionworld };
+
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+ struct PriVarScalarDataInfo { unsigned int pvIdx; std::string name; };
+ struct PriVarVectorDataInfo { std::vector<unsigned int> pvIdx; std::string name; };
+
+ struct FaceVarScalarDataInfo { std::function<Scalar(const FaceVariables&)> get; std::string name; };
+ struct FaceVarVectorDataInfo { std::function<GlobalPosition(const SubControlVolumeFace& scvf, const FaceVariables&)> get; std::string name; };
+
+ using Positions = std::vector<Scalar>;
+ using Data = std::vector<std::vector<Scalar>>;
+
+public:
+
+ StaggeredVtkOutputModule(const Problem& problem,
+ const FVGridGeometry& fvGridGeometry,
+ const GridVariables& gridVariables,
+ const SolutionVector& sol,
+ const std::string& name,
+ bool verbose = true,
+ Dune::VTK::DataMode dm = Dune::VTK::conforming)
+ : ParentType(problem, fvGridGeometry, gridVariables, sol, name, verbose, dm)
+ , problem_(problem)
+ , gridGeom_(fvGridGeometry)
+ , gridVariables_(gridVariables)
+ , sol_(sol)
+ , faceWriter_(std::make_shared<PointCloudVtkWriter<Scalar, dim>>(coordinates_))
+ , sequenceWriter_(faceWriter_, problem.name() + "-face", "","",
+ fvGridGeometry.gridView().comm().rank(),
+ fvGridGeometry.gridView().comm().size() )
+
+ {
+ writeFaceVars_ = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Vtk.WriteFaceData", false);
+ coordinatesInitialized_ = false;
+ }
+
+ //////////////////////////////////////////////////////////////////////////////////////////////
+ //! Methods to conveniently add primary and secondary variables upon problem initialization
+ //! Do not call these methods after initialization
+ //////////////////////////////////////////////////////////////////////////////////////////////
+
+ template<typename Vector>
+ void addFaceField(const Vector& v, const std::string& name, int nComp = 1)
+ {
+ if (v.size() == this->gridGeom_.gridView().size(1))
+ int x;
+ // fields_.emplace_back(gridGeom_.gridView(), gridGeom_.elementMapper(), v, name, nComp, 0);
+ else
+ DUNE_THROW(Dune::RangeError, "Size mismatch of added field!");
+ }
+
+ void addFaceVariable(std::function<Scalar(const FaceVariables&)>&& f, const std::string& name)
+ {
+ faceVarScalarDataInfo_.push_back(FaceVarScalarDataInfo{f, name});
+ }
+
+ void addFaceVariable(std::function<GlobalPosition(const SubControlVolumeFace& scvf, const FaceVariables&)>&& f, const std::string& name)
+ {
+ faceVarVectorDataInfo_.push_back(FaceVarVectorDataInfo{f, name});
+ }
+
+
+ void write(double time, Dune::VTK::OutputType type = Dune::VTK::ascii)
+ {
+ ParentType::write(time, type);
+ if(writeFaceVars_)
+ getFaceDataAndWrite_(time);
+ }
+
+
+private:
+
+ void updateCoordinates_()
+ {
+ coordinates_.resize(gridGeom_.numFaceDofs());
+ for(auto&& facet : facets(gridGeom_.gridView()))
+ {
+ const int dofIdxGlobal = gridGeom_.gridView().indexSet().index(facet);
+ coordinates_[dofIdxGlobal] = facet.geometry().center();
+ }
+ coordinatesInitialized_ = true;
+ }
+
+ /*!
+ * \brief Gathers all face-related data and invokes the face vtk-writer using these data.
+ */
+ void getFaceDataAndWrite_(const Scalar time)
+ {
+ const auto numPoints = gridGeom_.numFaceDofs();
+
+ // make sure not to iterate over the same dofs twice
+ std::vector<bool> dofVisited(numPoints, false);
+
+ // get fields for all primary coordinates and variables
+ if(!coordinatesInitialized_)
+ updateCoordinates_();
+
+ std::vector<std::vector<Scalar>> faceVarScalarData;
+ std::vector<std::vector<GlobalPosition>> faceVarVectorData;
+
+ if(!faceVarScalarDataInfo_.empty())
+ faceVarScalarData.resize(faceVarScalarDataInfo_.size(), std::vector<Scalar>(numPoints));
+
+ if(!faceVarVectorDataInfo_.empty())
+ faceVarVectorData.resize(faceVarVectorDataInfo_.size(), std::vector<GlobalPosition>(numPoints));
+
+ for (const auto& element : elements(gridGeom_.gridView(), Dune::Partitions::interior))
+ {
+ auto fvGeometry = localView(gridGeom_);
+ auto elemFaceVars = localView(gridVariables_.curGridFaceVars());
+
+ if (!faceVarScalarDataInfo_.empty() || !faceVarVectorDataInfo_.empty())
+ {
+ fvGeometry.bind(element);
+ elemFaceVars.bindElement(element, fvGeometry, sol_);
+
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ const auto dofIdxGlobal = scvf.dofIndex();
+ if(dofVisited[dofIdxGlobal])
+ continue;
+
+ dofVisited[dofIdxGlobal] = true;
+
+ const auto& faceVars = elemFaceVars[scvf];
+
+ for (std::size_t i = 0; i < faceVarScalarDataInfo_.size(); ++i)
+ faceVarScalarData[i][dofIdxGlobal] = faceVarScalarDataInfo_[i].get(faceVars);
+
+ for (std::size_t i = 0; i < faceVarVectorDataInfo_.size(); ++i)
+ faceVarVectorData[i][dofIdxGlobal] = faceVarVectorDataInfo_[i].get(scvf, faceVars);
+ }
+ }
+ }
+
+ if(!faceVarScalarDataInfo_.empty())
+ for (std::size_t i = 0; i < faceVarScalarDataInfo_.size(); ++i)
+ faceWriter_->addPointData(faceVarScalarData[i], faceVarScalarDataInfo_[i].name);
+
+ if(!faceVarVectorDataInfo_.empty())
+ for (std::size_t i = 0; i < faceVarVectorDataInfo_.size(); ++i)
+ faceWriter_->addPointData(faceVarVectorData[i], faceVarVectorDataInfo_[i].name, 3);
+
+ sequenceWriter_.write(time);
+ // faceScalarFields_.clear();
+ // faceVectorFields_.clear();
+ }
+
+
+ const Problem& problem_;
+ const FVGridGeometry& gridGeom_;
+ const GridVariables& gridVariables_;
+ const SolutionVector& sol_;
+
+ std::shared_ptr<PointCloudVtkWriter<Scalar, dim>> faceWriter_;
+
+ VTKSequenceWriter<PointCloudVtkWriter<Scalar, dim>> sequenceWriter_;
+
+ bool writeFaceVars_;
+
+ std::vector<GlobalPosition> coordinates_;
+ bool coordinatesInitialized_;
+
+ std::vector<FaceVarScalarDataInfo> faceVarScalarDataInfo_;
+ std::vector<FaceVarVectorDataInfo> faceVarVectorDataInfo_;
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/freeflow/staggered/test_donea.cc b/test/freeflow/staggered/test_donea.cc
index b3afe34f9627861b2836c03668fafb7a99b78a57..70ba090206f6f734e5d50991475336c0119c38fb 100644
--- a/test/freeflow/staggered/test_donea.cc
+++ b/test/freeflow/staggered/test_donea.cc
@@ -53,6 +53,7 @@
#include <dumux/discretization/methods.hh>
#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
@@ -157,7 +158,7 @@ int main(int argc, char** argv) try
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
- VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ StaggeredVtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //! Add model specific output fields
vtkWriter.addField(problem->getAnalyticalPressureSolution(), "pressureExact", 1);
vtkWriter.addField(problem->getAnalyticalVelocitySolution(), "velocityExact", GridView::dimensionworld);