Commit 8f9bd0d9 authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Timo Koch
Browse files

[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
parent 079f8329
......@@ -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;
}
};
}
......
......@@ -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;
}
};
......
......@@ -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); }
......
......@@ -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)
{