mmdg.hh

#ifndef DUNE_MMDG_HH
#define DUNE_MMDG_HH

#include <dune/mmdg/dg.hh>

template<class Grid, class GridView, class Mapper, class InterfaceGridView,
  class InterfaceMapper, class Problem>
class MMDG : public DG<Grid, GridView, Mapper, Problem>
{
public:
  static constexpr int dim = GridView::dimension;

  using Base = DG<Grid, GridView, Mapper, Problem>;
  using Scalar = typename Base::Scalar;
  using Matrix = typename Base::Matrix;
  using Vector = typename Base::Vector;
  using Coordinate = Dune::FieldVector<Scalar, dim>;
  using QRuleInterface = Base::QRuleEdge;
  using QRulesInterface = Base::QRuleEdges;
  static constexpr auto iSimplex = Dune::GeometryTypes::simplex(dim-1);

  //constructor
  MMDG (const Grid& grid, const GridView& gridView, const Mapper& mapper,
    const InterfaceGridView& iGridView, const InterfaceMapper& iMapper,
    const Problem& problem) :
    Base(grid, gridView, mapper, problem,
      (1 + dim) * gridView.size(0) + dim * iGridView.size(0)),
    iGridView_(iGridView), iMapper_(iMapper) {}

  const void operator() (const Scalar mu)
  {
    assembleSLE(mu);
    Base::A->compress();
  }

  const InterfaceGridView& iGridView_;
  const InterfaceMapper& iMapper_;

private:
  //assemble stiffness matrix A and load vector b
  void assembleSLE (const Scalar mu)
  {
    //quadrature rules
    const QRuleInterface& rule_qInterface =
      QRulesInterface::rule(iSimplex, problem_.quadratureOrder_qInterface());

    //we use the bulk basis
    // phi_elem,0 (x) = indicator(elem);
    // phi_elem,i (x) = x[i]*indicator(elem);
    //for elem in elements(gridView) and i = 1,...,dim

    //we use the interface basis
    // phi_iElem,0 (x) = indicator(elem);
    // phi_iElem,i (x) = (x * tau_i)*indicator(elem);
    //for iElem in elements(iGridView) and i = 1,...,dim-1
    //where {tau_i | i=1,...,dim-1} is the interface coordinate system given
    //by the method interfaceFrame(normal)

    //in total we use the basis given by the union of
    // { (phi_elem,i , 0) | elem in elements(gridView), i = 0,...,dim }
    //and
    // { (0 , phi_iElem,i ) | iElem in elements(iGridView), i = 0,...,dim-1 }

    //thus the stiffness matrix A and the load vector b exhibit the following
    //block structure
    // A = [ bulk,     coupling                b = [ bulk,
    //       coupling, interface],                   interface ]

    //call base class method for the assignment to fill stiffness matrix A
    //and load vector b with bulk entries