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 InterfaceBasis = std::array<Coordinate, dim-1>;
  using QRuleInterface = typename Base::QRuleEdge;
  using QRulesInterface = typename Base::QRulesEdge;
  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,
      Base::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
    Base::assembleSLE(mu);

    for (const auto& iElem : elements(iGridView_))
    {
      const int iElemIdx = iMapper_.index(iElem);
      const auto& iGeo = iElem.geometry();

      //get basis of the interface coordinate system for evaluation of
      //basis function phi_iElem,i with i=1,...,dim-1
      const std::array<Coordinate, dim-1> iFrame = interfaceFrame(
        Base::grid_.asIntersection(iElem).centerUnitOuterNormal());

      //in the total system of linear equations (SLE) Ad = b,
      //the index iElemIdxSLE refers to the basis function (0, phi_iElem,0)
      //and the indices iElemIdxSLE + i + 1, i = 1,...,dim-1, refer to the
      //basis functions (0, phi_iElem,i)
      const int iElemIdxSLE = (dim + 1)*Base::gridView_.size(0) + dim*iElemIdx;

      //fill load vector b with interface entries using a quadrature rule
      for (const auto& ip : rule_qInterface)
      {
        //quadrature point in local and global coordinates
        const auto& qpLocal = ip.position();
        const auto& qpGlobal = iGeo.global(qpLocal);

        const Scalar weight(ip.weight());
        const Scalar integrationElem(iGeo.integrationElement(qpLocal));
        const Scalar qInterfaceEvaluation(Base::problem_.qInterface(qpGlobal));

        //quadrature for int_elem q*phi_elem,0 dV
        Base::b[iElemIdxSLE] += weight * integrationElem * qInterfaceEvaluation;

        //quadrature for int_elem q*phi_elem,i dV
        for (int i = 0; i < dim - 1; i++)
        {
          Base::b[iElemIdxSLE + i + 1] += weight * integrationElem *
            (qpGlobal * iFrame[i]) * qInterfaceEvaluation;
        }
      }

      //TODO: fill stiffness matrix A with interface entries

      //fill stiffness matrix A with coupling entries
      for (const auto& elem : elements(Base::gridView_))
      {
        //TODO
      }
    }
  }

  //returns a basis of the interface coordinate system
  //2d: tau = [ normal[1],
  //           -normal[0] ]
  //3d: tau_1 ~ [ normal[1] + normal[2],     tau_2 ~ [ -normal[1],
  //             -normal[0],                            normal[0] + normal[2]
  //             -normal[0]             ],             -normal[1]            ]
  InterfaceBasis interfaceFrame(const Coordinate& normal)
  {
    InterfaceBasis iBasis;

    for (int i = 0; i < dim - 1; i++)
    {
      Coordinate tau(0.0);
      for (int j = 0; j < dim - 1; j++)
      {
        if (j == i)
        {
          continue;
        }

        tau[i] += normal[j];
        tau[j] = -normal[i];
      }
      tau /= tau.two_norm();

      iBasis[i] = tau;
    }

    return iBasis;

/* //NOTE: which version is better?
    if constexpr (dim == 2)
    {
      Coordinate tau;
      tau[0] = normal[1];
      tau[1] = -normal[0];

      interfaceFrame = {tau};
    }
    else if constexpr (dim == 3)
    {
      Coordinate tau1;
      tau1[0] = normal[1] + normal[2];
      tau1[1] = -normal[0];
      tau1[2] = -normal[0];
      Coordinate tau2;
      tau2[0] = -normal[1];
      tau2[1] = normal[0] + normal[2];
      tau2[2] = -normal[1];
      interfaceFrame = {tau1, tau2};
    }
    else
    {
      DUNE_THROW(Dune::Exception,
        "Invalid dimension: dim = " + std::to_string(dim));
    }*/
  }

};

#endif