dg.hh

#ifndef DUNE_MMDG_DG_HH
#define DUNE_MMDG_DG_HH

#include <cmath>

#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>

#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/io/file/vtk.hh>

#include <dune/geometry/quadraturerules.hh>

#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/solver.hh>
#include <dune/istl/umfpack.hh>

#include <dune/mmdg/nonconformingp1vtkfunction.hh>

template<class Grid, class GridView, class Mapper, class Problem>
class DG
{ //TODO: Ordnung (private / public)
public:
  //using GridView = typename Grid::LeafGridView;
  //using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
  using Scalar = typename GridView::ctype;
  static constexpr int dim = GridView::dimension;
  using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
  using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
  using VTKFunction = Dune::VTKFunction<GridView>;
  using P1Function = Dune::NonconformingP1VTKFunction<GridView,
    Dune::DynamicVector<Scalar>>;
  using QRuleElem = Dune::QuadratureRule<Scalar, dim>;
  using QRulesElem = Dune::QuadratureRules<Scalar, dim>;
  using QRuleEdge = Dune::QuadratureRule<Scalar, dim-1>;
  using QRulesEdge = Dune::QuadratureRules<Scalar, dim-1>;
  static constexpr auto simplex = Dune::GeometryTypes::simplex(dim);
  static constexpr auto edge = Dune::GeometryTypes::simplex(dim-1);

  //constructor
  DG (const Grid& grid, const GridView& gridView, const Mapper mapper, const Problem& problem) :
  //NOTE: why does this not work?
    //grid_(grid), gridView_(grid.leafGridView()), mapper_(*(new Mapper(gridView_, Dune::mcmgElementLayout()))), problem_(problem), dof_((1 + dim) * gridView_.size(0))
    grid_(grid), gridView_(gridView), mapper_(mapper), problem_(problem),
    dof_((1 + dim) * gridView.size(0))
    {
      //initialize stiffness matrix A, load vector b and solution vector d
      A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
      b = Vector(dof_);
      d = Vector(dof_);
    }

  const void operator() (const Scalar mu)
  {
    //assemble stiffness matrix A and load vector b
    assembleSLE(mu);
    A->compress();

    Dune::InverseOperatorResult result;
    Dune::UMFPack<Matrix> solver(*A);
    solver.apply(d, b, result);

    //write solution to a vtk file
    writeVTKoutput();
  }

  const Grid& grid_;
  const GridView& gridView_;
  const Mapper& mapper_; //element mapper for gridView_
  const Problem& problem_; //the DG problem to be solved