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/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 GridView, class Mapper, class Problem>
class DG
{
public:
  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>>;

  //constructor
  DG (const GridView& gridView, const Mapper& mapper,
    const Problem& problem) :
    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 K, const Scalar mu)
  {
    //assemble stiffness matrix A and load vector b
    assembleSLE(K, mu);

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

    //storage for pressure data, for each element we store the
    //pressure at the corners of the element, the VTKFunction will
    //create the output using a linear interpolation for each element
    Dune::DynamicVector<Scalar> pressure(dof, 0.0);
    Dune::DynamicVector<double> exactPressure(dof, 0.0);

    for (const auto& elem : elements(gridView_))
    {
      const int elemIdxSLE = (dim + 1)*mapper_.index(elem);
      const auto& geo = elem.geometry();

      for (int k = 0; k < geo.corners(); k++)
      {
        exactPressure[elemIdxSLE + k] = problem_.exactSolution(geo.corner(k));

        //contribution of the basis function
        // phi_elem,0 (x) = indicator(elem);
        //at the kth corner of elem
        pressure[elemIdxSLE + k] = d[elemIdxSLE];

        for (int i = 0; i < dim; i++)
        {