diff --git a/dune/mmdg/nonconformingp1vtkfunction.hh b/dune/mmdg/nonconformingp1vtkfunction.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9ae9ec80968a149730263320846d4e4715cbfb0e
--- /dev/null
+++ b/dune/mmdg/nonconformingp1vtkfunction.hh
@@ -0,0 +1,146 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifndef DUNE_MMFD_FUNCTION_HH
+#define DUNE_MMFD_FUNCTION_HH
+
+#include <string>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/geometry/type.hh>
+#include <dune/geometry/referenceelements.hh>
+#include <dune/geometry/multilineargeometry.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/grid/io/file/vtk/common.hh>
+#include <dune/grid/io/file/vtk/function.hh>
+
+namespace Dune
+{
+ //////////////////////////////////////////////////////////////////////
+ //
+ // Non-conforming P1VTKFunction
+ //
+
+ //! Take a vector and interpret it as point data for the VTKWriter
+ /**
+ * This class turns a generic vector containing point data into a
+ * VTKFunction. The vector must allow read access to the data via
+ * operator[]() and store the data in the order given by
+ * MultipleCodimMultipleGeomTypeMapper with a layout class that allows only
+ * vertices. Also, it must support the method size().
+ *
+ * While the number of components of the function is always 1, the vector
+ * may represent a field with multiple components of which one may be
+ * selected.
+ *
+ * \tparam GV Type of GridView the vector applies to.
+ * \tparam V Type of vector.
+ */
+ template<typename GV, typename V>
+ class NonconformingP1VTKFunction
+ : public VTKFunction< GV >
+ {
+ //! Base class
+ typedef VTKFunction< GV > Base;
+ //! Mapper for vertices
+ typedef MultipleCodimMultipleGeomTypeMapper<GV> Mapper;
+
+ //! store a reference to the vector
+ const V& v;
+ //! name of this function
+ std::string s;
+ //! number of components of the field stored in the vector
+ int ncomps_;
+ //! index of the component of the field in the vector this function is
+ //! responsible for
+ int mycomp_;
+ //! number of corners per element
+ int corners_;
+ //! precision with which to output the field
+ VTK::Precision prec_;
+ //! mapper used to map elements to indices
+ Mapper mapper;
+
+ public:
+ typedef typename Base::Entity Entity;
+ typedef typename Base::ctype ctype;
+ using Base::dim;
+
+ //! return number of components
+ int ncomps () const override
+ {
+ return 1;
+ }
+
+ //! evaluate
+ double evaluate (int comp, const Entity& e,
+ const Dune::FieldVector<ctype,dim>& xi) const override
+ {
+ const unsigned int dim = Entity::mydimension;
+ const unsigned int nVertices = e.subEntities(dim);
+
+ std::vector<FieldVector<ctype,1> > cornerValues(nVertices);
+ for (unsigned i=0; i<nVertices; ++i)
+ cornerValues[i] = v[(corners_*ncomps_)*mapper.index(e)+ncomps_*i+mycomp_];
+
+ // (Ab)use the MultiLinearGeometry class to do multi-linear interpolation between scalars
+ const MultiLinearGeometry<ctype,dim,1> interpolation(e.type(), cornerValues);
+ return interpolation.global( xi );
+ }
+
+ //! get name
+ std::string name () const override
+ {
+ return s;
+ }
+
+ //! get output precision for the field
+ VTK::Precision precision() const override
+ {
+ return prec_;
+ }
+
+ //! construct from a vector and a name
+ /**
+ * \param gv GridView to operate on (used to instantiate a
+ * MultipleCodimMultipleGeomTypeMapper, otherwise no
+ * reference or copy is stored). Note that this must be the
+ * GridView the vector applies to as well as the GridView
+ * later used by the VTKWriter -- i.e. we do not implicitly
+ * restrict or prolongate the data.
+ * \param v_ Reference to the vector holding the data. The reference
+ * is stored internally and must be valid for as long as
+ * this functions evaluate method is used.
+ * \param s_ Name of this function in the VTK file.
+ * \param ncomps Number of components of the field represented by the
+ * vector.
+ * \param mycomp Number of the field component this function is
+ * responsible for.
+ * \param prec the precision with which to output the field
+ */
+ NonconformingP1VTKFunction(const GV& gv, const V &v_, const std::string &s_,
+ int ncomps=1, int mycomp=0, VTK::Precision prec = VTK::Precision::float32)
+ : v( v_ ),
+ s( s_ ),
+ ncomps_(ncomps),
+ mycomp_(mycomp),
+ corners_( GV::dimension+1 ),
+ prec_(prec),
+ mapper( gv, mcmgElementLayout() )
+ {
+ if (v.size()!=(unsigned int)(corners_*mapper.size()*ncomps_))
+ DUNE_THROW(IOError,"NonconformingP1VTKFunction: size mismatch");
+ }
+
+ //! destructor
+ virtual ~NonconformingP1VTKFunction() {}
+ };
+
+ //! \} group VTK
+
+} // namespace Dune
+
+#endif // DUNE_MMFD_FUNCTION_HH
diff --git a/src/dg.cc b/src/dg.cc
index 23d2a9b00cfb03197f8069c021aa19079090e9b1..b51f281fc12aa94c172438f3c9e2abc21ec208a4 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -2,6 +2,9 @@
# include "config.h"
#endif
+#include <sys/stat.h>
+#include <sys/types.h>
+
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
@@ -9,14 +12,17 @@
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
-#include <dune/grid/yaspgrid.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
#include <dune/grid/io/file/dgfparser/dgfug.hh>
+#include <dune/grid/io/file/vtk.hh>
#include <dune/geometry/quadraturerules.hh>
+#include <dune/mmesh/mmesh.hh>
+
#include <dune/mmdg/clockhelper.hh>
+#include <dune/mmdg/nonconformingp1vtkfunction.hh>
int main(int argc, char** argv)
{
@@ -29,7 +35,7 @@ int main(int argc, char** argv)
static constexpr int dim = GRIDDIM;
- using Grid = Dune::YaspGrid<dim>;
+ using Grid = Dune::MovingMesh<dim>;
using GridFactory = Dune::DGFGridFactory<Grid>;
using GridView = typename Grid::LeafGridView;
using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
@@ -113,6 +119,7 @@ int main(int argc, char** argv)
const double intersctVol = intersctGeo.volume();
const auto& intersctCenter = intersctGeo.center();
+ //TODO: quadrature rule cannot be used for dim = 1!
const Dune::QuadratureRule<double,dim-1>& secondOrderRule =
Dune::QuadratureRules<double,dim-1>::rule(intersct.type(), 2);
@@ -298,7 +305,57 @@ int main(int argc, char** argv)
A.solve(d, b);
- //NOTE: what is a suitable way to create the output?
+ //storage for pressure data, for each element we store the pressure at the
+ //corners of the element, the VTKOperatorP1 will create the output using a
+ //linear interpolation for each element
+ Dune::DynamicVector<double> pressure(3*gridView.size(0), 0.0);
+
+ for (const auto& elem : elements(gridView))
+ {
+ const int elemIdx = mapper.index(elem);
+ const auto& geo = elem.geometry();
+
+ for (int k = 0; k < geo.corners(); k++)
+ {
+ //contribution of the basis function
+ // phi_elem,0 (x) = indicator(elem);
+ //at the kth corner of elem
+ pressure[3*elemIdx + k] = d[(1 + dim)*elemIdx];
+
+ for (int i = 0; i < dim; i++)
+ {
+ //contribution of the basis function
+ // phi_elem,i (x) = x[i]*indicator(elem);
+ //at the kth corner of elem
+ pressure[3*elemIdx + k] +=
+ d[(1 + dim)*elemIdx + i + 1] * geo.corner(k)[i];
+ }
+ }
+ }
+
+/*
+ typedef Dune::VTKFunction< GridView > VTKFunction;
+ typedef Dune::NonconformingP1VTKFunction<GridView,
+ Dune::DynamicVector<double>> P1Function;
+ VTKFunction* p1 = new P1Function(gridView, lp, "linearpressure");
+ vtkwriter.addVertexData( std::shared_ptr< const VTKFunction >(p1) );
+ */
+
+ //create output directory if necessary
+ mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+
+ auto vtkWriter = std::make_shared<Dune::VTKWriter<GridView>>(gridView);
+ Dune::VTKSequenceWriter<GridView>
+ vtkSqWriter(vtkWriter, "pressure", "data", "data");
+
+ using VTKOperatorP1 = Dune::VTKFunction<GridView>;
+ using P1Function = Dune::NonconformingP1VTKFunction<GridView,
+ Dune::DynamicVector<double>>;
+ VTKOperatorP1* vtkOperator = new P1Function(gridView, pressure, "pressure");
+
+
+ vtkSqWriter.addVertexData( std::shared_ptr<const VTKOperatorP1>(vtkOperator) );
+ vtkSqWriter.write(0.0);
//print total run time
const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);