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/mmdg/nonconformingp1vtkfunction.hh>

template<class GridView, class Mapper, class Problem>
class DG
{
public:
  using Scalar = typename GridView::ctype;
  static constexpr int dim = GridView::dimension;

  //constructor
  DG (const GridView& gridView, const Mapper& mapper,
    const Problem& problem) :
    gridView_(gridView), mapper_(mapper), problem_(problem) {}

  const void operator() (const Scalar K, const Scalar mu) const
  {
    //NOTE: what is an appropriate sparse matrix type? -> BCRS
    const int dof = (1 + dim)*gridView_.size(0); //degrees of freedom
    using Matrix = Dune::DynamicMatrix<Scalar>;
    using Vector = Dune::DynamicVector<Scalar>;
    Matrix A(dof, dof, 0.0); //stiffness matrix
    Vector b(dof, 0.0); //load vector
    Vector d(dof, 0.0);

    //NOTE:
//    const int order = 3; //order of the quadrature rule

    //we use the basis
    // phi_elem,0 (x) = indicator(elem);
    // phi_elem,i (x) = x[i]*indicator(elem);
    //for elem in elements(gridView) and i = 1,...,dim
    for (const auto& elem : elements(gridView_))
    {
      const int elemIdx = mapper_.index(elem);
      const auto& geo = elem.geometry();
      const double elemVol = geo.volume();
      const auto& center = geo.center();

      //in the system of linear equations (SLE) Ad = b,
      //the index elemIdxSLE refers to the basis function phi_elem,0
      //and the indices elemIdxSLE + i + 1, i = 1,...,dim, refer to the
      //basis function phi_elem,i
      const int elemIdxSLE = (dim + 1)*elemIdx;

  /*    //TODO: can be done outside of the loop?
      const Dune::QuadratureRule<double,dim>& rule =
        Dune::QuadratureRules<double,dim>::rule(geo.type(),order);

      //NOTE: how are quadrature rules in Dune applied correctly?
      for (const auto& ip : rule)
      {
        const auto qp = ip.position();
        //NOTE: is the volume of the element taken into account
        //automatically?
        const auto weight = ip.weight();

        //quadrature for int_elem q*phi_elem,0 dV
        b[elemIdxSLE] += weight * q(qp);