add VTK output, use MMesh instead of YaspGrid

dune/mmdg/nonconformingp1vtkfunction.hh

+// -*- 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

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);