From 8f9bd0d9b08e79fcd329b38f8b47edcce65ddccd Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Wed, 8 Feb 2017 16:31:15 +0100
Subject: [PATCH] [staggeredGrid][vtk] Revise staggered vtk writer
* Now provides void addPointData() to make usage more similar to existing writer
* Will get an interface which allows use with dune sequence writer
---
dumux/freeflow/staggered/model.hh | 110 -------
dumux/freeflow/staggered/vtkoutputmodule.hh | 6 +-
dumux/io/staggeredvtkoutputmodule.hh | 166 +++++-----
dumux/io/staggeredvtkwriter.hh | 301 ++++++++++--------
.../staggered/kovasznaytestproblem.hh | 42 +--
5 files changed, 279 insertions(+), 346 deletions(-)
diff --git a/dumux/freeflow/staggered/model.hh b/dumux/freeflow/staggered/model.hh
index 2856c75b03..16dbdaef44 100644
--- a/dumux/freeflow/staggered/model.hh
+++ b/dumux/freeflow/staggered/model.hh
@@ -93,116 +93,6 @@ public:
// NonIsothermalModel::maybeAddTemperature(vtkOutputModule);
}
-
- /*!
- * \brief \copybrief Dumux::ImplicitModel::addOutputVtkFields
- *
- * Specialization for the NavierStokesModel, adding the pressure and
- * the process rank to the VTK writer.
- */
- template<class MultiWriter>
- void addOutputVtkFields(const SolutionVector &sol,
- MultiWriter &writer)
- {
- // TODO: implement vtk output properly, account for 3d
- using VectorField = Dune::BlockVector<Dune::FieldVector<double, dimWorld> >;
-
- // create the required scalar fields
- const auto numElements = this->gridView_().size(0);
- auto *p = writer.allocateManagedBuffer(numElements);
- auto *delP = writer.allocateManagedBuffer(numElements);
-
- auto *v_x_pos = writer.allocateManagedBuffer(numElements);
- auto *v_x_neg = writer.allocateManagedBuffer(numElements);
- auto *v_y_pos = writer.allocateManagedBuffer(numElements);
- auto *v_y_neg = writer.allocateManagedBuffer(numElements);
-
- VectorField *velocity = writer.template allocateManagedBuffer<double, dimWorld>(numElements);
-
- // initialize velocity field
- for (unsigned int i = 0; i < numElements; ++i)
- {
- (*velocity)[i] = Scalar(0);
- }
-
- auto *rank = writer.allocateManagedBuffer(numElements);
-
- for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
- {
- auto eIdx = this->elementMapper().index(element);
- (*rank)[eIdx] = this->gridView_().comm().rank();
-
- // get the local fv geometry
- auto fvGeometry = localView(this->globalFvGeometry());
- fvGeometry.bindElement(element);
-
- auto elemVolVars = localView(this->curGlobalVolVars());
- elemVolVars.bindElement(element, fvGeometry, this->curSol());
-
- for (auto&& scv : scvs(fvGeometry))
- {
- const auto& volVars = elemVolVars[scv];
- auto dofIdxGlobal = scv.dofIndex();
-
- (*p)[dofIdxGlobal] = volVars.pressure();
- (*delP)[dofIdxGlobal] = volVars.pressure() - 1.1e5;
-
- GlobalPosition velocityVector(0.0);
- for (auto&& scvf : scvfs(fvGeometry))
- {
- auto& origFaceVars = this->curGlobalFaceVars().faceVars(scvf.dofIndex());
- auto dirIdx = scvf.directionIndex();
-
- velocityVector[dirIdx] += 0.5*origFaceVars.velocity();
-
- if(scvf.unitOuterNormal()[dirIdx] > 0.0)
- {
- if(dirIdx == 0)
- (*v_x_pos)[dofIdxGlobal] = origFaceVars.velocity();
- if(dirIdx == 1)
- (*v_y_pos)[dofIdxGlobal] = origFaceVars.velocity();
- }
- else
- {
- if(dirIdx == 0)
- (*v_x_neg)[dofIdxGlobal] = origFaceVars.velocity();
- if(dirIdx == 1)
- (*v_y_neg)[dofIdxGlobal] = origFaceVars.velocity();
- }
- }
- (*velocity)[dofIdxGlobal] = velocityVector;
- }
- }
- writer.attachDofData(*p, "p", isBox);
- writer.attachDofData(*delP, "delP", isBox);
- writer.attachDofData(*v_x_pos, "v_x_pos", isBox);
- writer.attachDofData(*v_x_neg, "v_x_neg", isBox);
- writer.attachDofData(*v_y_pos, "v_y_pos", isBox);
- writer.attachDofData(*v_y_neg, "v_y_neg", isBox);
- writer.attachCellData(*rank, "process rank");
- writer.attachDofData(*velocity, "velocity", isBox, dim);
- }
-
- auto velocity(const Element& element) const
- {
- GlobalPosition velocityVector(0.0);
-
- // get the local fv geometry
- auto fvGeometry = localView(this->globalFvGeometry());
- fvGeometry.bindElement(element);
- for (auto&& scv : scvs(fvGeometry))
- {
- for (auto&& scvf : scvfs(fvGeometry))
- {
- auto& origFaceVars = this->curGlobalFaceVars().faceVars(scvf.dofIndex());
- auto dirIdx = scvf.directionIndex();
- velocityVector[dirIdx] += 0.5*origFaceVars.velocity();
- }
- }
- return velocityVector;
- }
-
-
};
}
diff --git a/dumux/freeflow/staggered/vtkoutputmodule.hh b/dumux/freeflow/staggered/vtkoutputmodule.hh
index f3fe5f22d3..3a55f7c828 100644
--- a/dumux/freeflow/staggered/vtkoutputmodule.hh
+++ b/dumux/freeflow/staggered/vtkoutputmodule.hh
@@ -71,13 +71,13 @@ private:
* \param face The face
*/
template<class Face>
- void getVectorData_(Data& priVarVectorData, const Face& face)
+ void getPrivarVectorData_(Data& priVarVectorData, const Face& face)
{
const int dofIdxGlobal = face.dofIndex();
const int dirIdx = directionIndex(face.unitOuterNormal());
const Scalar velocity = this->problem().model().curSol()[faceIdx][dofIdxGlobal][0];
- for (int i = 0; i < this->faceData().priVarVectorDataInfo.size(); ++i)
- priVarVectorData[i][dofIdxGlobal * this->faceData().priVarVectorDataInfo[i].pvIdx.size() + dirIdx] = velocity;
+ for (int i = 0; i < this->priVarVectorDataInfo_.size(); ++i)
+ priVarVectorData[i][dofIdxGlobal * this->priVarVectorDataInfo_[i].pvIdx.size() + dirIdx] = velocity;
}
};
diff --git a/dumux/io/staggeredvtkoutputmodule.hh b/dumux/io/staggeredvtkoutputmodule.hh
index 4a6ab7c82c..acbeb5afe3 100644
--- a/dumux/io/staggeredvtkoutputmodule.hh
+++ b/dumux/io/staggeredvtkoutputmodule.hh
@@ -55,6 +55,8 @@ class StaggeredVtkOutputModule : public VtkOutputModuleBase<TypeTag>
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;
@@ -65,31 +67,14 @@ class StaggeredVtkOutputModule : public VtkOutputModuleBase<TypeTag>
using Positions = std::vector<Scalar>;
using Data = std::vector<std::vector<Scalar>>;
- // a collection of types and variables to be passed to the staggered vtk writer
- struct WriterData
- {
- using Positions = StaggeredVtkOutputModule::Positions;
- using PriVarScalarDataInfo = StaggeredVtkOutputModule::PriVarScalarDataInfo;
- using PriVarVectorDataInfo = StaggeredVtkOutputModule::PriVarVectorDataInfo;
- using Data = StaggeredVtkOutputModule::Data;
- Data priVarScalarData;
- Data priVarVectorData;
- Data secondVarScalarData;
- Data secondVarVectorData;
- Positions positions;
- std::vector<PriVarScalarDataInfo> priVarScalarDataInfo;
- std::vector<PriVarVectorDataInfo> priVarVectorDataInfo;
- std::vector<PriVarScalarDataInfo> secondVarScalarDataInfo;
- std::vector<PriVarVectorDataInfo> secondVarVectorDataInfo;
- };
-
public:
StaggeredVtkOutputModule(const Problem& problem,
- Dune::VTK::DataMode dm = Dune::VTK::conforming) : ParentType(problem, dm), faceWriter_(problem)
+ Dune::VTK::DataMode dm = Dune::VTK::conforming) : ParentType(problem, dm), faceWriter_(coordinates_)
{
writeFaceVars_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, WriteFaceData);
+ coordinatesInitialized_ = false;
}
//////////////////////////////////////////////////////////////////////////////////////////////
@@ -97,12 +82,31 @@ public:
//! Do not call these methods after initialization
//////////////////////////////////////////////////////////////////////////////////////////////
+
+ //! Output a scalar field
+ //! \param name The name of the vtk field
+ //! \returns A reference to the resized scalar field to be filled with the actual data
+ std::vector<Scalar>& createFaceScalarField(const std::string& name)
+ {
+ faceScalarFields_.emplace_back(std::make_pair(std::vector<Scalar>(this->problem().model().numFaceDofs()), name));
+ return faceScalarFields_.back().first;
+ }
+
+ //! Output a vector field
+ //! \param name The name of the vtk field
+ //! \returns A reference to the resized vector field to be filled with the actual data
+ std::vector<GlobalPosition>& createFaceVectorField(const std::string& name)
+ {
+ faceVectorFields_.emplace_back(std::make_pair(std::vector<GlobalPosition>(this->problem().model().numFaceDofs()), name));
+ return faceVectorFields_.back().first;
+ }
+
//! Output a scalar primary variable
//! \param name The name of the vtk field
//! \param pvIdx The index in the primary variables vector
void addFacePrimaryVariable(const std::string& name, unsigned int pvIdx)
{
- faceData_.priVarScalarDataInfo.push_back(PriVarScalarDataInfo{pvIdx, name});
+ priVarScalarDataInfo_.push_back(PriVarScalarDataInfo{pvIdx, name});
}
//! Output a vector primary variable
@@ -111,7 +115,7 @@ public:
void addFacePrimaryVariable(const std::string& name, std::vector<unsigned int> pvIndices)
{
assert(pvIndices.size() < 4 && "Vtk doesn't support vector dimensions greater than 3!");
- faceData_.priVarVectorDataInfo.push_back(PriVarVectorDataInfo{pvIndices, name});
+ priVarVectorDataInfo_.push_back(PriVarVectorDataInfo{pvIndices, name});
}
void write(double time, Dune::VTK::OutputType type = Dune::VTK::ascii)
@@ -142,16 +146,25 @@ protected:
return this->problem().model().curSol()[cellCenterIdx][dofIdxGlobal][pvIdx];
}
- /*!
- * \brief Returns a reference to the face data
- */
- const WriterData& faceData() const
- {
- return faceData_;
- }
+ std::vector<PriVarScalarDataInfo> priVarScalarDataInfo_;
+ std::vector<PriVarVectorDataInfo> priVarVectorDataInfo_;
+ std::vector<PriVarScalarDataInfo> secondVarScalarDataInfo_;
+ std::vector<PriVarVectorDataInfo> secondVarVectorDataInfo_;
private:
+ void updateCoordinates_()
+ {
+ std::cout << "updating coordinates" << std::endl;
+ coordinates_.resize(this->problem().model().numFaceDofs());
+ for(auto&& facet : facets(this->problem().gridView()))
+ {
+ const int dofIdxGlobal = this->problem().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.
*/
@@ -163,13 +176,14 @@ private:
std::vector<bool> dofVisited(numPoints, false);
// get fields for all primary coordinates and variables
- Positions positions(numPoints*3);
+ if(!coordinatesInitialized_)
+ updateCoordinates_();
- Data priVarScalarData(faceData_.priVarScalarDataInfo.size(), std::vector<Scalar>(numPoints));
+ Data priVarScalarData(priVarScalarDataInfo_.size(), std::vector<Scalar>(numPoints));
- Data priVarVectorData(faceData_.priVarVectorDataInfo.size());
- for (std::size_t i = 0; i < faceData_.priVarVectorDataInfo.size(); ++i)
- priVarVectorData[i].resize(numPoints*faceData_.priVarVectorDataInfo[i].pvIdx.size());
+ Data priVarVectorData(priVarVectorDataInfo_.size());
+ for (std::size_t i = 0; i < priVarVectorDataInfo_.size(); ++i)
+ priVarVectorData[i].resize(numPoints*priVarVectorDataInfo_[i].pvIdx.size());
for(auto&& element : elements(this->problem().gridView()))
{
@@ -180,54 +194,35 @@ private:
if(dofVisited[scvf.dofIndex()])
continue;
- asImp_().getPositions_(positions, scvf);
- asImp_().getScalarData_(priVarScalarData, scvf);
- asImp_().getVectorData_(priVarVectorData, scvf);
+ asImp_().getPrivarScalarData_(priVarScalarData, scvf);
+ asImp_().getPrivarVectorData_(priVarVectorData, scvf);
+
dofVisited[scvf.dofIndex()] = true;
}
}
- faceData_.priVarScalarData = std::move(priVarScalarData);
- faceData_.priVarVectorData = std::move(priVarVectorData);
- faceData_.positions = std::move(positions);
-// results.secondVarScalarData = xxx; //TODO: implemented secondVarData
-// results.secondVarScalarData = xxx;
+ // transfer priVar scalar data to writer
+ for(int i = 0; i < priVarScalarDataInfo_.size(); ++i)
+ faceWriter_.addPointData(priVarScalarData[i], priVarScalarDataInfo_[i].name);
- faceWriter_.write(faceData_);
- }
+ // transfer priVar vector data to writer
+ for(int i = 0; i < priVarVectorDataInfo_.size(); ++i)
+ faceWriter_.addPointData(priVarVectorData[i], priVarVectorDataInfo_[i].name, priVarVectorDataInfo_[i].pvIdx.size());
- /*!
- * \brief Retrives the position (center) of the face.
- * ParaView expects three-dimensional coordinates, therefore we fill empty entries with zero for dim < 3.
- *
- * \param positions Container to store the positions
- * \param face The face
- */
- template<class Face>
- void getPositions_(Positions& positions, const Face& face)
- {
- const int dofIdxGlobal = face.dofIndex();
- auto pos = face.center();
- if(dim == 1)
- {
- positions[dofIdxGlobal*3] = pos[0];
- positions[dofIdxGlobal*3 + 1] = 0.0;
- positions[dofIdxGlobal*3 + 2] = 0.0;
- }
- else if(dim == 2)
- {
- positions[dofIdxGlobal*3] = pos[0];
- positions[dofIdxGlobal*3 + 1] = pos[1];
- positions[dofIdxGlobal*3 + 2] = 0.0;
- }
- else
- {
- positions[dofIdxGlobal*3] = pos[0];
- positions[dofIdxGlobal*3 + 1] = pos[1] ;
- positions[dofIdxGlobal*3 + 2] = pos[2] ;
- }
+ // transfer custom scalar data to writer
+ for(auto&& scalarField : faceScalarFields_)
+ faceWriter_.addPointData(scalarField.first, scalarField.second);
+
+ // transfer custom vector data to writer
+ for(auto&& vectorField : faceVectorFields_)
+ faceWriter_.addPointData(vectorField.first, vectorField.second, 3);
+
+ faceWriter_.write();
+ faceScalarFields_.clear();
+ faceVectorFields_.clear();
}
+
/*!
* \brief Retrives scalar-valued data from the face.
*
@@ -235,13 +230,11 @@ private:
* \param face The face
*/
template<class Face>
- void getScalarData_(Data& priVarScalarData, const Face& face)
+ void getPrivarScalarData_(Data& priVarScalarData, const Face& face)
{
const int dofIdxGlobal = face.dofIndex();
- for(int pvIdx = 0; pvIdx < faceData_.priVarScalarDataInfo.size(); ++pvIdx)
- {
+ for(int pvIdx = 0; pvIdx < priVarScalarDataInfo_.size(); ++pvIdx)
priVarScalarData[pvIdx][dofIdxGlobal] = this->problem().model().curSol()[faceIdx][dofIdxGlobal][pvIdx];
- }
}
/*!
@@ -251,19 +244,24 @@ private:
* \param face The face
*/
template<class Face>
- void getVectorData_(Data& priVarVectorData, const Face& face)
+ void getPrivarVectorData_(Data& priVarVectorData, const Face& face)
{
const int dofIdxGlobal = face.dofIndex();
- for (int i = 0; i < faceData_.priVarVectorDataInfo.size(); ++i)
- for (int j = 0; j < faceData_.priVarVectorDataInfo[i].pvIdx.size(); ++j)
- priVarVectorData[i][dofIdxGlobal*faceData_.priVarVectorDataInfo[i].pvIdx.size() + j]
- = this->problem().model().curSol()[faceIdx][dofIdxGlobal][faceData_.priVarVectorDataInfo[i].pvIdx[j]];
+ for (int i = 0; i < priVarVectorDataInfo_.size(); ++i)
+ for (int j = 0; j < priVarVectorDataInfo_[i].pvIdx.size(); ++j)
+ priVarVectorData[i][dofIdxGlobal*priVarVectorDataInfo_[i].pvIdx.size() + j]
+ = this->problem().model().curSol()[faceIdx][dofIdxGlobal][priVarVectorDataInfo_[i].pvIdx[j]];
}
- StaggeredVtkWriter<TypeTag, WriterData> faceWriter_;
- WriterData faceData_;
+ StaggeredVtkWriter<dimWorld> faceWriter_;
bool writeFaceVars_;
+ std::vector<GlobalPosition> coordinates_;
+ bool coordinatesInitialized_;
+
+ std::list<std::pair<std::vector<Scalar>, std::string>> faceScalarFields_;
+ std::list<std::pair<std::vector<GlobalPosition>, std::string>> faceVectorFields_;
+
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
{ return *static_cast<Implementation *>(this); }
diff --git a/dumux/io/staggeredvtkwriter.hh b/dumux/io/staggeredvtkwriter.hh
index 8db5bbbc59..f262bcd59a 100644
--- a/dumux/io/staggeredvtkwriter.hh
+++ b/dumux/io/staggeredvtkwriter.hh
@@ -27,6 +27,7 @@
#include <dumux/io/vtkoutputmodulebase.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dune/grid/io/file/vtk/common.hh>
namespace Properties
{
@@ -46,108 +47,174 @@ namespace Dumux
* initialization and/or be turned on/off using the designated properties. Additionally
* non-standardized scalar and vector fields can be added to the writer manually.
*/
-template<class TypeTag, class WriterData>
+template<int dim>
class StaggeredVtkWriter
{
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
- typename DofTypeIndices::CellCenterIdx cellCenterIdx;
- typename DofTypeIndices::FaceIdx faceIdx;
-
- enum { dim = GridView::dimension };
- enum { dimWorld = GridView::dimensionworld };
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ using Scalar = double;
+ using GlobalPosition = Dune::FieldVector<Scalar, dim>;
static constexpr unsigned int precision = 6;
static constexpr unsigned int numBeforeLineBreak = 15;
- using ScalarInfo = std::vector<typename WriterData::PriVarScalarDataInfo>;
- using VectorInfo = std::vector<typename WriterData::PriVarVectorDataInfo>;
- using Positions = typename WriterData::Positions;
- using Data = typename WriterData::Data;
+ template<class ContainerType>
+ class VTKLocalFunction
+ {
+ public:
+ VTKLocalFunction(const ContainerType& data, const std::string& name, const int numComponents) : data_(data), name_(name), numComponents_(numComponents)
+ {}
+
+ const std::string& name() const
+ {
+ return name_;
+ }
+
+ int numComponents() const
+ {
+ return numComponents_;
+ }
+
+ auto& operator() (const int dofIdx) const { return data_[dofIdx]; }
+
+ decltype(auto) begin() const
+ {
+ return data_.begin();
+ }
+
+ decltype(auto) end() const
+ {
+ return data_.end();
+ }
+
+ int size() const
+ {
+ return data_.size();
+ }
+
+ private:
+ const ContainerType& data_;
+ const std::string name_;
+ const int numComponents_;
+ };
+
public:
- StaggeredVtkWriter(const Problem& problem) : problem_(problem)
+ using ScalarLocalFunction = VTKLocalFunction<std::vector<Scalar>>;
+ using VectorLocalFunction = VTKLocalFunction<std::vector<GlobalPosition>>;
+
+
+ StaggeredVtkWriter(const std::vector<GlobalPosition>& coordinates) : coordinates_(coordinates)
{}
- void write(const WriterData& data)
+ void write(/*const WriterData& data*/)
{
- if(data.priVarScalarDataInfo.empty() &&
- data.priVarVectorDataInfo.empty() &&
- data.secondVarScalarDataInfo.empty() &&
- data.secondVarVectorDataInfo.empty())
- return;
file_.open(fileName_());
- write_(data);
+ writeHeader_();
+ writeCoordinates_(coordinates_);
+ writeDataInfo_();
+
+ for(auto&& data : scalarPointData_)
+ {
+ writeData_(data);
+ }
+ for(auto&& data :vectorPointData_)
+ {
+ writeData_(data);
+ }
+
+ file_ << "</PointData>\n";
+ file_ << "</Piece>\n";
+ file_ << "</PolyData>\n";
+ file_ << "</VTKFile>";
+
+ clear();
+
file_.close();
++curWriterNum_;
}
+
+ std::string write ( const std::string &name,
+ Dune::VTK::OutputType type = Dune::VTK::ascii )
+ {
+ return "dummy";
+ }
+
+
+ void addPointData(const std::vector<Scalar>& v, const std::string &name, int ncomps = 1)
+ {
+ assert(v.size() == ncomps * coordinates_.size());
+ scalarPointData_.push_back(ScalarLocalFunction(v, name, ncomps));
+ }
+
+ void addPointData(const std::vector<GlobalPosition>& v, const std::string &name, int ncomps = 1)
+ {
+ assert(v.size() == coordinates_.size());
+ vectorPointData_.push_back(VectorLocalFunction(v, name, ncomps));
+ }
+
+
+ void clear()
+ {
+ scalarPointData_.clear();
+ vectorPointData_.clear();
+ }
+
+
private:
- void writeHeader_(const int numPoints)
+ void writeHeader_()
{
std::string header = "<?xml version=\"1.0\"?>\n";
header += "<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
header += "<PolyData>\n";
- header += "<Piece NumberOfLines=\"0\" NumberOfPoints=\"" + std::to_string(numPoints) + "\">\n";
+ header += "<Piece NumberOfLines=\"0\" NumberOfPoints=\"" + std::to_string(coordinates_.size()) + "\">\n";
file_ << header;
}
/*!
- * \brief Writes the coordinates and the simulation data to the file
- *
- * \param priVarScalarData Container to store the scalar-valued data
- * \param priVarVectorData Container to store the vector-valued data
- * \param scalarInfo Information concerning the data (e.g. names)
- * \param vectorInfo Information concerning the data (e.g. names)
- * \param positions Container to store the positions
+ * \brief Writes information about the data to the file
*/
- void write_(const WriterData& data)
+ void writeDataInfo_()
{
- const int numPoints = problem_.model().numFaceDofs();
- writeHeader_(numPoints);
- writeCoordinates_(data.positions);
-
std::string scalarName;
std::string vectorName;
+ bool foundScalar = false;
+ bool foundVector = false;
- bool scalarValuesPresent = false;
- bool vectorValuesPresent = false;
-
-
- if(!(data.priVarScalarDataInfo.empty() && data.secondVarScalarDataInfo.empty()))
+ for(auto&& data : scalarPointData_)
{
- scalarValuesPresent = true;
- scalarName = data.priVarScalarDataInfo.empty() ? data.secondVarScalarDataInfo[0].name :
- data.priVarScalarDataInfo[0].name;
+ if(data.numComponents() == 1 && !foundScalar)
+ {
+ scalarName = data.name();
+ foundScalar = true;
+ continue;
+ }
+
+ if(data.numComponents() > 1 && !foundVector)
+ {
+ vectorName = data.name();
+ foundVector = true;
+ }
}
- if(!(data.priVarVectorDataInfo.empty() && data.secondVarVectorDataInfo.empty()))
+
+ for(auto&& data : vectorPointData_)
{
- vectorValuesPresent = true;
- vectorName = data.priVarVectorDataInfo.empty() ? data.secondVarVectorDataInfo[0].name :
- data.priVarVectorDataInfo[0].name;
+ if(data.numComponents() > 1 && !foundVector)
+ {
+ vectorName = data.name();
+ foundVector = true;
+ }
}
- if(scalarValuesPresent)
- if(vectorValuesPresent)
- file_ << "<PointData Scalars=\"" << scalarName << "\" Vectors=\"" << vectorName <<"\">";
+ if(foundScalar)
+ if(foundVector)
+ file_ << "<PointData Scalars=\"" << scalarName << "\" Vectors=\"" << vectorName <<"\">\n";
else
- file_ << "<PointData Scalars=\"" << scalarName << "\">";
- else if(vectorValuesPresent)
- file_ << "<PointData Vectors=\"" << vectorName << "\">";
+ file_ << "<PointData Scalars=\"" << scalarName << "\">\n";
+ else if(foundVector)
+ file_ << "<PointData Vectors=\"" << vectorName << "\">\n";
else
return;
-
- file_ << std::endl;
- writeAllData_(data);
-
- file_ << "</PointData>\n";
- file_ << "</Piece>\n";
- file_ << "</PolyData>\n";
- file_ << "</VTKFile>";
}
/*!
@@ -155,7 +222,7 @@ private:
*
* \param positions Container to store the positions
*/
- void writeCoordinates_(const Positions& positions)
+ void writeCoordinates_(const std::vector<GlobalPosition>& positions)
{
// write the positions to the file
file_ << "<Points>\n";
@@ -163,7 +230,13 @@ private:
int counter = 0;
for(auto&& x : positions)
{
- file_ << x << " ";
+ file_ << x ;
+
+ if(x.size() == 1)
+ file_ << " 0 0 ";
+ if(x.size() == 2)
+ file_ << " 0 ";
+
// introduce a line break after a certain time
if((++counter) > numBeforeLineBreak)
{
@@ -175,92 +248,72 @@ private:
file_ << "</Points>\n";
}
- void writeAllData_(const WriterData& data)
- {
- // write scalar-valued primary variable data
- if(!data.priVarScalarDataInfo.empty())
- writeScalarData_(data.priVarScalarData, data.priVarScalarDataInfo);
-
- // write vector-valued primary variable data
- if(!data.priVarVectorDataInfo.empty())
- writeVectorData_(data.priVarVectorData, data.priVarVectorDataInfo);
-
- // write scalar-valued secondary variable data
- if(!data.secondVarScalarDataInfo.empty())
- writeScalarData_(data.secondVarScalarData, data.secondVarScalarDataInfo);
-
- // write vector-valued primary variable data
- if(!data.secondVarVectorDataInfo.empty())
- writeVectorData_(data.secondVarVectorData, data.secondVarVectorDataInfo);
- }
-
/*!
- * \brief Writes scalar-valued data to the file
+ * \brief Writes data to the file
*
- * \param priVarScalarData Container to store the data
- * \param info Information concerning the data (e.g. names)
+ * \param data The data container which hold the data itself, as well as the name of the data set and the number of its components
*/
- void writeScalarData_(const Data& priVarScalarData, const ScalarInfo& info)
+ template<class T>
+ void writeData_(const T& data)
{
- // write the priVars to the file
- for(int pvIdx = 0; pvIdx < info.size(); ++pvIdx)
+ file_ << "<DataArray type=\"Float32\" Name=\"" << data.name() << "\" NumberOfComponents=\"" << data.numComponents() << "\" format=\"ascii\">\n";
+ int counter = 0;
+ for(auto&& value : data)
{
- file_ << "<DataArray type=\"Float32\" Name=\"" << info[pvIdx].name << "\" NumberOfComponents=\"1\" format=\"ascii\">\n";
- int counter = 0;
- for(auto&& x : priVarScalarData[pvIdx])
+ // forward to specialized function
+ writeToFile_(value);
+
+ // introduce a line break after a certain time
+ if((++counter) > numBeforeLineBreak)
{
- file_ << x << " " ;
- // introduce a line break after a certain time
- if((++counter) > numBeforeLineBreak)
- {
- file_ << std::endl;
- counter = 0;
- }
+ file_ << std::endl;
+ counter = 0;
}
- file_ << "\n</DataArray>\n";
}
+ file_ << "\n</DataArray>\n";
}
/*!
- * \brief Writes vector-valued data to the file
+ * \brief Writes a scalar to the file
*
- * \param priVarVectorData Container to store the data
- * \param info Information concerning the data (e.g. names)
+ * \param s The scalar
*/
- void writeVectorData_(const Data& priVarVectorData, const VectorInfo& info)
+ void writeToFile_(const Scalar& s)
{
- // write the priVars to the file
- for(int i = 0; i < info.size(); ++i)
- {
- file_ << "<DataArray type=\"Float32\" Name=\"" << info[i].name << "\" NumberOfComponents=\"" << info[i].pvIdx.size() << "\" format=\"ascii\">\n";
- int counter = 0;
- for(auto&& x : priVarVectorData[i])
- {
- file_ << x << " " ;
- // introduce a line break after a certain time
- if((++counter) > numBeforeLineBreak)
- {
- file_ << std::endl;
- counter = 0;
- }
- }
- file_ << "\n</DataArray>\n";
- }
+ file_ << s << " ";
+ }
+
+ /*!
+ * \brief Writes a vector to the file
+ *
+ * \param g The vector
+ */
+ void writeToFile_(const GlobalPosition& g)
+ {
+ assert(g.size() > 1 && g.size() < 4);
+ if(g.size() < 3)
+ file_ << g << " 0 ";
+ else
+ file_ << g;
}
+
/*!
* \brief Returns the file name for each timestep
*/
std::string fileName_() const
{
std::ostringstream oss;
- oss << problem_.name() << "-face-" << std::setw(5) << std::setfill('0') << curWriterNum_ << ".vtp";
+ oss << /*problem_.name() <<*/ "face-" << std::setw(5) << std::setfill('0') << curWriterNum_ << ".vtp";
return oss.str();
}
- const Problem& problem_;
+ const std::vector<GlobalPosition>& coordinates_;
std::ofstream file_;
int curWriterNum_{0};
+
+ std::list<ScalarLocalFunction> scalarPointData_;
+ std::list<VectorLocalFunction> vectorPointData_;
};
} // end namespace Dumux
diff --git a/test/freeflow/staggered/kovasznaytestproblem.hh b/test/freeflow/staggered/kovasznaytestproblem.hh
index c274b3a1c2..75c07c4961 100644
--- a/test/freeflow/staggered/kovasznaytestproblem.hh
+++ b/test/freeflow/staggered/kovasznaytestproblem.hh
@@ -299,52 +299,44 @@ public:
return values;
}
+
/*!
- * \brief Append all quantities of interest which can be derived
- * from the solution of the current time step to the VTK
- * writer.
+ * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
*/
- void addOutputVtkFields()
+ template<class VtkOutputModule>
+ void addVtkOutputFields(VtkOutputModule& outputModule) const
{
- //Here we calculate the analytical solution
- const auto numElements = this->gridView().size(0);
- auto& pExact = *(this->resultWriter().allocateManagedBuffer(numElements));
- auto& velocityExact = *(this->resultWriter()).template allocateManagedBuffer<double, dimWorld>(numElements);
-
- using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
- typename DofTypeIndices::FaceIdx faceIdx;
-
- const auto someElement = *(elements(this->gridView()).begin());
-
- auto someFvGeometry = localView(this->model().globalFvGeometry());
- someFvGeometry.bindElement(someElement);
- const auto someScv = *(scvs(someFvGeometry).begin());
+ auto& pressureExact = outputModule.createScalarField("pressureExact", 0);
+ auto& velocityExact = outputModule.createVectorField("velocityExact", dim);
- Scalar time = std::max(this->timeManager().time() + this->timeManager().timeStepSize(), 1e-10);
+ auto& scalarFaceVelocityExact = outputModule.createFaceScalarField("scalarFaceVelocityExact");
+ auto& vectorFaceVelocityExact = outputModule.createFaceVectorField("vectorFaceVelocityExact");
for (const auto& element : elements(this->gridView()))
{
auto fvGeometry = localView(this->model().globalFvGeometry());
fvGeometry.bindElement(element);
-
for (auto&& scv : scvs(fvGeometry))
{
auto dofIdxGlobal = scv.dofIndex();
- const auto& globalPos = scv.dofPosition();
-
- pExact[dofIdxGlobal] = dirichletAtPos(globalPos)[pressureIdx];
+ auto dofPosition = scv.dofPosition();
+ pressureExact[dofIdxGlobal] = dirichletAtPos(dofPosition)[pressureIdx];
GlobalPosition velocityVector(0.0);
for (auto&& scvf : scvfs(fvGeometry))
{
auto dirIdx = scvf.directionIndex();
- velocityVector[dirIdx] += 0.5*dirichletAtPos(globalPos)[faceIdx][dirIdx];
+ auto analyticalSolution = dirichletAtPos(dofPosition)[faceIdx][dirIdx];
+ velocityVector[dirIdx] += 0.5*analyticalSolution;
+ scalarFaceVelocityExact[dofIdxGlobal] = analyticalSolution;
+
+ GlobalPosition tmp(0.0);
+ tmp[dirIdx] = analyticalSolution;
+ vectorFaceVelocityExact[dofIdxGlobal] = std::move(tmp);
}
velocityExact[dofIdxGlobal] = velocityVector;
}
}
- this->resultWriter().attachDofData(pExact, "p_exact", false);
- this->resultWriter().attachDofData(velocityExact, "velocity_exact", false, dim);
}
/*!
--
GitLab