diff --git a/dumux/io/vtkoutputmodule.hh b/dumux/io/vtkoutputmodule.hh
index daf41b0b1933458f5d65ac27d752d8d95830fe03..f1f115f9d7e906541b89d57dfa36335a234c56fb 100644
--- a/dumux/io/vtkoutputmodule.hh
+++ b/dumux/io/vtkoutputmodule.hh
@@ -235,6 +235,24 @@ public:
}
}
+protected:
+ // some return functions for differing implementations to use
+ const Problem& problem() const { return problem_; }
+ const FVGridGeometry& fvGridGeometry() const { return gridGeom_; }
+ const GridVariables& gridVariables() const { return gridVariables_; }
+ const SolutionVector& sol() const { return sol_; }
+
+ const bool verbose() const { return verbose_; }
+ const std::string& name() const { return name_; }
+ const Dune::VTK::DataMode dataMode() const { return dm_; }
+
+ Dune::VTKWriter<GridView>& writer() { return *writer_; }
+ Dune::VTKSequenceWriter<GridView>& sequenceWriter() { return sequenceWriter_; }
+
+ const std::vector<VolVarScalarDataInfo>& volVarScalarDataInfo() const { return volVarScalarDataInfo_; }
+ const std::vector<VolVarVectorDataInfo>& volVarVectorDataInfo() const { return volVarVectorDataInfo_; }
+ const std::vector<Field>& fields() const { return fields_; }
+
private:
//! Assembles the fields and adds them to the writer (conforming output)
diff --git a/dumux/porousmediumflow/boxdfm/vtkoutputmodule.hh b/dumux/porousmediumflow/boxdfm/vtkoutputmodule.hh
new file mode 100644
index 0000000000000000000000000000000000000000..59b06c098804654602a9b2c772f7e5c1fe258d45
--- /dev/null
+++ b/dumux/porousmediumflow/boxdfm/vtkoutputmodule.hh
@@ -0,0 +1,416 @@
+// -*- 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
+ * \ingroup InputOutput
+ * \brief A VTK output module to simplify writing dumux simulation data to VTK format
+ */
+#ifndef POROUSMEDIUMFLOW_BOXDFM_VTK_OUTPUT_MODULE_HH
+#define POROUSMEDIUMFLOW_BOXDFM_VTK_OUTPUT_MODULE_HH
+
+#include <set>
+
+#include <dune/grid/common/gridfactory.hh>
+#include <dune/geometry/referenceelements.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup InputOutput
+ * \brief A VTK output module to simplify writing dumux simulation data to VTK format.
+ * This output module is specialized for writing out data obtained by the box-dfm
+ * scheme. It writes out separate vtk files for the solution on the fracture. For
+ * this, a grid type to be used the fracture-conforming lower-dimensional grid has
+ * to be provided.
+ *
+ * \tparam TypeTag The TypeTag of the problem implementation
+ * \tparam FractureGrid The Type used for the lower-dimensional grid
+ * \tparam phaseIdxOffset Used for single-phase problems to retrieve the right phase name
+ *
+ * Handles the output of scalar and vector fields to VTK formatted file for multiple
+ * variables and timesteps. Certain predefined fields can be registered on
+ * 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 FractureGrid, int phaseIdxOffset = 0>
+class BoxDfmVtkOutputModule : public VtkOutputModule<TypeTag, phaseIdxOffset>
+{
+ using ParentType = VtkOutputModule<TypeTag, phaseIdxOffset>;
+
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using VelocityOutput = typename GET_PROP_TYPE(TypeTag, VelocityOutput);
+ using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+ static constexpr int numPhases = ModelTraits::numPhases();
+
+ using VV = typename GridVariables::VolumeVariables;
+ using FluidSystem = typename VV::FluidSystem;
+ using Scalar = typename GridVariables::Scalar;
+
+ using GridView = typename FVGridGeometry::GridView;
+ using FractureGridView = typename FractureGrid::LeafGridView;
+ using FractureVertexMapper = Dune::MultipleCodimMultipleGeomTypeMapper<FractureGridView>;
+
+ enum
+ {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using ReferenceElements = typename Dune::ReferenceElements<typename GridView::ctype, dim>;
+
+ using Field = Vtk::template Field<GridView>;
+ using FractureField = Vtk::template Field<FractureGridView>;
+
+ static_assert(dim > 1, "Box-Dfm output only works for dim > 1");
+ static_assert(FractureGrid::dimension == int(dim-1), "Fracture grid must be of codimension one!");
+ static_assert(FractureGrid::dimensionworld == int(dimWorld), "Fracture grid has to has the same coordinate dimension!");
+ static_assert(FVGridGeometry::discMethod == DiscretizationMethod::box, "Box-Dfm output module can only be used with the box scheme!");
+public:
+
+ BoxDfmVtkOutputModule(const Problem& problem,
+ const FVGridGeometry& fvGridGeometry,
+ const GridVariables& gridVariables,
+ const SolutionVector& sol,
+ const std::string& name,
+ const std::string& paramGroup = "",
+ Dune::VTK::DataMode dm = Dune::VTK::conforming,
+ bool verbose = true)
+ : ParentType(problem, fvGridGeometry, gridVariables, sol, name, paramGroup, dm, verbose)
+ {
+ const auto& gridView = fvGridGeometry.gridView();
+ Dune::GridFactory<FractureGrid> gridFactory;
+
+ // insert fracture vertices
+ std::size_t fracVertexCount = 0;
+ vertexToFractureVertexIdx_.resize(gridView.size(dim));
+ for (const auto& v : vertices(gridView))
+ {
+ const auto vIdxGlobal = fvGridGeometry.vertexMapper().index(v);
+ if (fvGridGeometry.dofOnFracture(vIdxGlobal))
+ {
+ gridFactory.insertVertex(v.geometry().center());
+ vertexToFractureVertexIdx_[vIdxGlobal] = fracVertexCount++;
+ }
+ }
+
+ // insert fracture elements
+ std::set< std::pair<IndexType, unsigned int> > handledFacets;
+ auto vertIsOnFracture = [&fvGridGeometry] (auto idx) { return fvGridGeometry.dofOnFracture(idx); };
+ for (const auto& element : elements(gridView))
+ {
+ const auto referenceElement = ReferenceElements::general(element.geometry().type());
+
+ for (const auto& is : intersections(gridView, element))
+ {
+ // obtain all vertex indices on this intersection
+ const auto& isGeometry = is.geometry();
+ const auto numCorners = isGeometry.corners();
+
+ std::vector<IndexType> isVertexIndices(numCorners);
+ for (unsigned int i = 0; i < numCorners; ++i)
+ {
+ const auto vIdxLocal = referenceElement.subEntity(is.indexInInside(), 1, i, dim);
+ isVertexIndices[i] = fvGridGeometry.vertexMapper().subIndex(element, vIdxLocal, dim);
+ }
+
+ // determine if this is a fracture facet & if it has to be inserted
+ bool insertFacet = false;
+ if ( std::all_of(isVertexIndices.begin(), isVertexIndices.end(), vertIsOnFracture) )
+ {
+ insertFacet = true;
+ if (!is.boundary())
+ {
+ // only proceed if facet has not been handled yet
+ const auto pair = std::make_pair<IndexType, unsigned int>(fvGridGeometry.elementMapper().index(is.outside()),
+ is.indexInOutside());
+ if (handledFacets.count( pair ) != 0)
+ insertFacet = false;
+ }
+ }
+
+ if (insertFacet)
+ {
+ // transform intersection vertex indices to frac grid indices
+ std::for_each( isVertexIndices.begin(),
+ isVertexIndices.end(),
+ [&] (auto& idx) { idx = this->vertexToFractureVertexIdx_[idx]; } );
+
+ // insert the element
+ gridFactory.insertElement(isGeometry.type(), isVertexIndices);
+
+ // insert to set of handled facets
+ handledFacets.insert( std::make_pair<IndexType, unsigned int>(fvGridGeometry.elementMapper().index(element),
+ is.indexInInside()) );
+ }
+ }
+ }
+
+ auto grid = gridFactory.createGrid();
+ fractureGridView_ = std::make_unique<FractureGridView>(grid->leafGridView());
+
+ // update fracture vertex mapper
+ fractureVertexMapper_ = std::make_unique<FractureVertexMapper>(*fractureGridView_, Dune::mcmgVertexLayout());
+
+ // obtain map fracture vertex insertion index -> fracture grid vertex index
+ std::vector<IndexType> insToVertexIdx(fractureGridView_->size(FractureGridView::dimension));
+ for (const auto& v : vertices(*fractureGridView_))
+ insToVertexIdx[ gridFactory.insertionIndex(v) ] = fractureVertexMapper_->index(v);
+
+ // update the index map
+ for (IndexType dofIdx = 0; dofIdx < gridView.size(GridView::dimension); ++dofIdx)
+ if (fvGridGeometry.dofOnFracture(dofIdx))
+ vertexToFractureVertexIdx_[dofIdx] = insToVertexIdx[ vertexToFractureVertexIdx_[dofIdx] ];
+
+ // instantiate writers for the fracture
+ fractureWriter_ = std::make_shared< Dune::VTKWriter<FractureGridView> >(*fractureGridView_, dm);
+ fractureSequenceWriter_ = std::make_unique< Dune::VTKSequenceWriter<FractureGridView> >(fractureWriter_, name + "_fracture");
+ }
+
+ //////////////////////////////////////////////////////////////////////////////////////////////
+ //! Writing data
+ //////////////////////////////////////////////////////////////////////////////////////////////
+
+ //! Write the data for this timestep to file in four steps
+ //! (1) We assemble all registered variable fields
+ //! (2) We register them with the vtk writer
+ //! (3) The writer writes the output for us
+ //! (4) Clear the writer for the next time step
+ void write(double time, Dune::VTK::OutputType type = Dune::VTK::ascii)
+ {
+ Dune::Timer timer;
+
+ // write to file depending on data mode
+ const auto dm = this->dataMode();
+ if (dm == Dune::VTK::conforming)
+ writeConforming_(time, type);
+ else if (dm == Dune::VTK::nonconforming)
+ writeNonConforming_(time, type);
+ else
+ DUNE_THROW(Dune::NotImplemented, "Output for provided vtk data mode");
+
+ //! output
+ timer.stop();
+ if (this->verbose())
+ std::cout << "Writing output for problem \"" << this->name() << "\". Took " << timer.elapsed() << " seconds." << std::endl;
+ }
+
+private:
+ //! Assembles the fields and adds them to the writer (conforming output)
+ void writeConforming_(double time, Dune::VTK::OutputType type)
+ {
+ //////////////////////////////////////////////////////////////
+ //! (1) Assemble all variable fields and add to writer
+ //////////////////////////////////////////////////////////////
+
+ // instatiate the velocity output
+ VelocityOutput velocityOutput(this->problem(), this->fvGridGeometry(), this->gridVariables(), this->sol());
+ std::array<std::vector<GlobalPosition>, numPhases> velocity;
+
+ // process rank
+ static bool addProcessRank = getParamFromGroup<bool>(this->paramGroup(), "Vtk.AddProcessRank");
+ std::vector<double> rank;
+
+ // volume variable data
+ std::vector<std::vector<Scalar>> volVarScalarData;
+ std::vector<std::vector<Scalar>> volVarScalarDataFracture;
+ std::vector<std::vector<GlobalPosition>> volVarVectorData;
+ std::vector<std::vector<GlobalPosition>> volVarVectorDataFracture;
+
+ // some references for convenience
+ const auto& gridView = this->fvGridGeometry().gridView();
+ const auto& volVarScalarDataInfo = this->volVarScalarDataInfo();
+ const auto& volVarVectorDataInfo = this->volVarVectorDataInfo();
+
+ //! Abort if no data was registered
+ if (!volVarScalarDataInfo.empty()
+ || !volVarVectorDataInfo.empty()
+ || !this->fields().empty()
+ || velocityOutput.enableOutput()
+ || addProcessRank)
+ {
+ const auto numCells = gridView.size(0);
+ const auto numDofs = gridView.size(dim);
+ const auto numFractureVert = fractureGridView_->size(FractureGridView::dimension);
+
+ // get fields for all volume variables
+ if (!this->volVarScalarDataInfo().empty())
+ {
+ volVarScalarData.resize(volVarScalarDataInfo.size(), std::vector<Scalar>(numDofs));
+ volVarScalarDataFracture.resize(volVarScalarDataInfo.size(), std::vector<Scalar>(numFractureVert));
+ }
+ if (!this->volVarVectorDataInfo().empty())
+ {
+ volVarVectorData.resize(volVarVectorDataInfo.size(), std::vector<GlobalPosition>(numDofs));
+ volVarVectorDataFracture.resize(volVarVectorDataInfo.size(), std::vector<GlobalPosition>(numFractureVert));
+ }
+
+ if (velocityOutput.enableOutput())
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ velocity[phaseIdx].resize(numDofs);
+
+ // maybe allocate space for the process rank
+ if (addProcessRank) rank.resize(numCells);
+
+ for (const auto& element : elements(gridView, Dune::Partitions::interior))
+ {
+ const auto eIdxGlobal = this->fvGridGeometry().elementMapper().index(element);
+
+ auto fvGeometry = localView(this->fvGridGeometry());
+ auto elemVolVars = localView(this->gridVariables().curGridVolVars());
+
+ // If velocity output is enabled we need to bind to the whole stencil
+ // otherwise element-local data is sufficient
+ if (velocityOutput.enableOutput())
+ {
+ fvGeometry.bind(element);
+ elemVolVars.bind(element, fvGeometry, this->sol());
+ }
+ else
+ {
+ fvGeometry.bindElement(element);
+ elemVolVars.bindElement(element, fvGeometry, this->sol());
+ }
+
+ if (!volVarScalarDataInfo.empty() || !volVarVectorDataInfo.empty())
+ {
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto dofIdxGlobal = scv.dofIndex();
+ const auto& volVars = elemVolVars[scv];
+
+ // get the scalar-valued data
+ if (!scv.isOnFracture())
+ {
+ for (std::size_t i = 0; i < volVarScalarDataInfo.size(); ++i)
+ volVarScalarData[i][dofIdxGlobal] = volVarScalarDataInfo[i].get(volVars);
+ }
+ else
+ {
+ for (std::size_t i = 0; i < volVarScalarDataInfo.size(); ++i)
+ volVarScalarDataFracture[i][vertexToFractureVertexIdx_[dofIdxGlobal]] = volVarScalarDataInfo[i].get(volVars);
+ }
+
+ // get the vector-valued data
+ if (!scv.isOnFracture())
+ {
+ for (std::size_t i = 0; i < volVarVectorDataInfo.size(); ++i)
+ volVarVectorData[i][dofIdxGlobal] = volVarVectorDataInfo[i].get(volVars);
+ }
+ else
+ {
+ for (std::size_t i = 0; i < volVarVectorDataInfo.size(); ++i)
+ volVarVectorDataFracture[i][vertexToFractureVertexIdx_[dofIdxGlobal]] = volVarVectorDataInfo[i].get(volVars);
+ }
+ }
+ }
+
+ // velocity output
+ if (velocityOutput.enableOutput())
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ velocityOutput.calculateVelocity(velocity[phaseIdx], elemVolVars, fvGeometry, element, phaseIdx);
+
+ //! the rank
+ if (addProcessRank)
+ rank[eIdxGlobal] = static_cast<double>(gridView.comm().rank());
+ }
+
+ //////////////////////////////////////////////////////////////
+ //! (2) Register data fields with the vtk writer
+ //////////////////////////////////////////////////////////////
+
+ // volume variables if any
+ for (std::size_t i = 0; i < volVarScalarDataInfo.size(); ++i)
+ {
+ this->sequenceWriter().addVertexData(volVarScalarData[i], volVarScalarDataInfo[i].name);
+ fractureSequenceWriter_->addVertexData(volVarScalarDataFracture[i], volVarScalarDataInfo[i].name);
+ }
+
+ for (std::size_t i = 0; i < volVarVectorDataInfo.size(); ++i)
+ {
+ this->sequenceWriter().addVertexData( Field(gridView, this->fvGridGeometry().vertexMapper(), volVarVectorData[i],
+ volVarVectorDataInfo[i].name, /*numComp*/dimWorld, /*codim*/dim).get() );
+ fractureSequenceWriter_->addVertexData( FractureField(*fractureGridView_, *fractureVertexMapper_, volVarVectorDataFracture[i],
+ volVarVectorDataInfo[i].name, /*numComp*/dimWorld, /*codim*/dim-1).get() );
+ }
+
+ // the velocity field
+ if (velocityOutput.enableOutput())
+ {
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ this->sequenceWriter().addVertexData( Field(gridView, this->fvGridGeometry().vertexMapper(), velocity[phaseIdx],
+ "velocity_" + std::string(FluidSystem::phaseName(phaseIdx+phaseIdxOffset)) + " (m/s)",
+ /*numComp*/dimWorld, /*codim*/dim).get() );
+ }
+
+ // the process rank
+ if (addProcessRank)
+ this->sequenceWriter().addCellData( Field(gridView, this->fvGridGeometry().elementMapper(), rank,
+ "process rank", /*numComp*/1, /*codim*/0).get() );
+
+ // also register additional (non-standardized) user fields if any (only on matrix grid)
+ for (auto&& field : this->fields())
+ {
+ if (field.codim() == 0)
+ this->sequenceWriter().addCellData(field.get());
+ else if (field.codim() == dim)
+ this->sequenceWriter().addVertexData(field.get());
+ else
+ DUNE_THROW(Dune::RangeError, "Cannot add wrongly sized vtk scalar field!");
+ }
+ }
+
+ //////////////////////////////////////////////////////////////
+ //! (2) The writer writes the output for us
+ //////////////////////////////////////////////////////////////
+ this->sequenceWriter().write(time, type);
+ fractureSequenceWriter_->write(time, type);
+
+ //////////////////////////////////////////////////////////////
+ //! (3) Clear the writer
+ //////////////////////////////////////////////////////////////
+ this->writer().clear();
+ fractureWriter_->clear();
+ }
+
+ //! Assembles the fields and adds them to the writer (conforming output)
+ void writeNonConforming_(double time, Dune::VTK::OutputType type)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Non-conforming output for Box-Dfm model");
+ }
+
+ std::unique_ptr<FractureVertexMapper> fractureVertexMapper_;
+ std::unique_ptr<FractureGridView> fractureGridView_;
+ std::shared_ptr<Dune::VTKWriter<FractureGridView>> fractureWriter_;
+ std::unique_ptr< Dune::VTKSequenceWriter<FractureGridView> > fractureSequenceWriter_;
+
+ std::vector<IndexType> vertexToFractureVertexIdx_;
+};
+
+} // end namespace Dumux
+
+#endif