mmdg.hh

#ifndef DUNE_MMDG_HH
#define DUNE_MMDG_HH

#include <dune/mmdg/dg.hh>

template<class GridView, class Mapper, class InterfaceGridView,
  class InterfaceMapper, class Problem>
class MMDG : public DG<GridView, Mapper, Problem>
{
public:
  static constexpr int dim = GridView::dimension;
  static constexpr auto iSimplex = Dune::GeometryTypes::simplex(dim-1);

  using Base = DG<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 InterfaceBasis = std::array<Coordinate, dim-1>;
  using QRuleInterface = typename Base::QRuleEdge;
  using QRulesInterface = typename Base::QRulesEdge;
  using VTKFunction = Dune::VTKFunction<InterfaceGridView>;
  using P1Function = Dune::NonconformingP1VTKFunction<InterfaceGridView,
    Dune::DynamicVector<Scalar>>;

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

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

    for (int i = 0; i < Base::dof_; i++)
      for (int j = 0; j < i; j++)
        /*assert*/if(std::abs((*Base::A)[i][j] - (*Base::A)[j][i]) >
          std::numeric_limits<Scalar>::epsilon())
            std::cout << "NON-symmetric!: " << i << ", " << j << std::endl;

    std::cout << *Base::A << std::endl;

    (*Base::A).solve(Base::d, Base::b);

    std::cout << "\n\n" << Base::b << "\n\n" << Base::d << "\n\n";

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

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

  //compute L2 error using quadrature rules
  // norm = sqrt( ||.||^2_L2(bulk) + ||.||^2_L2(interface) )
  Scalar computeL2error(const int quadratureOrder = 5) const
  {
    const Scalar bulkError = Base::computeL2error(quadratureOrder);
    Scalar interfaceError(0.0);

    for (const auto& iElem : elements(iGridView_))
    {
      const auto& iGeo = iElem.geometry();
      const int iElemIdxSLE =
        (dim + 1) * Base::gridView_.size(0) + dim * iMapper_.index(iElem);