diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index be1b24ac60acd88c0f047352fb0cc2b701049a16..dda8217a47f3c2a6ccf67adac688e63cbafd125b 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -15,8 +15,9 @@ public:
   using Matrix = typename Base::Matrix;
   using Vector = typename Base::Vector;
   using Coordinate = Dune::FieldVector<Scalar, dim>;
-  using QRuleInterface = Base::QRuleEdge;
-  using QRulesInterface = Base::QRuleEdges;
+  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
@@ -41,8 +42,8 @@ private:
   void assembleSLE (const Scalar mu)
   {
     //quadrature rules
-    const QRuleInterface& rule_qInterface =
-      QRulesInterface::rule(iSimplex, problem_.quadratureOrder_qInterface());
+    const QRuleInterface& rule_qInterface = QRulesInterface::rule(iSimplex,
+      Base::problem_.quadratureOrder_qInterface());
 
     //we use the bulk basis
     // phi_elem,0 (x) = indicator(elem);
@@ -95,16 +96,16 @@ private:
 
         const Scalar weight(ip.weight());
         const Scalar integrationElem(iGeo.integrationElement(qpLocal));
-        const Scalar qInterfaceEvaluation(problem_.qInterface(qpGlobal));
+        const Scalar qInterfaceEvaluation(Base::problem_.qInterface(qpGlobal));
 
         //quadrature for int_elem q*phi_elem,0 dV
-        b[iElemIdxSLE] += weight * integrationElem * qEvaluation;
+        Base::b[iElemIdxSLE] += weight * integrationElem * qInterfaceEvaluation;
 
         //quadrature for int_elem q*phi_elem,i dV
         for (int i = 0; i < dim - 1; i++)
         {
-          b[iElemIdxSLE + i + 1] +=
-            weight * integrationElem * (qpGlobal * iFrame[i]) * qEvaluation;
+          Base::b[iElemIdxSLE + i + 1] += weight * integrationElem *
+            (qpGlobal * iFrame[i]) * qInterfaceEvaluation;
         }
       }
 
@@ -124,8 +125,10 @@ private:
   //3d: tau_1 ~ [ normal[1] + normal[2],     tau_2 ~ [ -normal[1],
   //             -normal[0],                            normal[0] + normal[2]
   //             -normal[0]             ],             -normal[1]            ]
-  std::array<Coordinate, dim-1> interfaceFrame(const Coordinate& normal)
+  InterfaceBasis interfaceFrame(const Coordinate& normal)
   {
+    InterfaceBasis iBasis;
+
     for (int i = 0; i < dim - 1; i++)
     {
       Coordinate tau(0.0);
@@ -141,8 +144,11 @@ private:
       }
       tau /= tau.two_norm();
 
-      interfaceFrame[i] = tau;
+      iBasis[i] = tau;
     }
+
+    return iBasis;
+
 /* //NOTE: which version is better?
     if constexpr (dim == 2)
     {
@@ -171,7 +177,6 @@ private:
     }*/
   }
 
-
 };
 
 #endif
diff --git a/dune/mmdg/problems/coupleddgproblem.hh b/dune/mmdg/problems/coupleddgproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e1abfffe97541843449c9014154875354d757930
--- /dev/null
+++ b/dune/mmdg/problems/coupleddgproblem.hh
@@ -0,0 +1,58 @@
+#ifndef __DUNE_MMDG_COUPLEDDGPROBLEM_HH__
+#define __DUNE_MMDG_COUPLEDDGPROBLEM_HH__
+
+#include <dune/mmdg/problems/dgproblem.hh>
+
+//abstract class providing an interface for a coupled flow problem in a
+//fractured porous medium that is to be solved with a coupled discontinuous
+//Galerkin scheme
+template<class Vector, class Scalar = double>
+class CoupledDGProblem : public DGProblem<Vector, Scalar>
+{
+  public:
+    static constexpr int dimension = Vector::dimension;
+    using Base = DGProblem<Vector, Scalar>;
+    using Matrix = typename Base::Matrix;
+
+    //interface source term at position pos
+    virtual Scalar qInterface (const Vector& pos) const
+    {
+      return Scalar(0.0);
+    }
+
+    //aperture of the fracture at position pos
+    virtual Scalar aperture (const Vector& pos) const
+    {
+      return Scalar(1.0);
+    }
+
+    //tangential permeability tensor of the interface at position pos
+    virtual Matrix Kparallel (const Vector& pos) const
+    {
+      Matrix permeability(0.0);
+
+      for (int i = 0; i < dimension; i++)
+      {
+        permeability[i][i] = 1.0;
+      }
+
+      //return identity matrix
+      return permeability;
+    }
+
+    //permeability of the fracture in normal direction at position pos
+    virtual Scalar Kperp (const Vector& pos) const
+    {
+      return Scalar(1.0);
+    }
+
+
+    //returns the recommended quadrature order to compute an integral
+    //over x * qInterface(x)
+    virtual int quadratureOrder_qInterface () const
+    {
+      return 1;
+    }
+};
+
+#endif
diff --git a/dune/mmdg/problems/poisson.hh b/dune/mmdg/problems/poisson.hh
index 516a7c014a352181aa4ca96020b2eebc7449d3ec..8781cb38180aae3584f1f48f517da1420d7a9c5c 100644
--- a/dune/mmdg/problems/poisson.hh
+++ b/dune/mmdg/problems/poisson.hh
@@ -8,11 +8,10 @@
 
 
 //Poisson problem with right side q = -1
-template<int dim, class ctype = double>
-class Poisson : public DGProblem<Dune::FieldVector<ctype, dim>, ctype>
+template<int dim, class Scalar = double>
+class Poisson : public DGProblem<Dune::FieldVector<Scalar, dim>, Scalar>
 {
   public:
-    using Scalar = ctype;
     using Vector = Dune::FieldVector<Scalar, dim>;
 
     //the exact solution at position pos
diff --git a/src/mmdg.cc b/src/mmdg.cc
index aea7f9a8660c145c718a1f1a56addde93fccb578..ed37ba1263c0c2d16bb5327c4fc9c7dde4ff481c 100644
--- a/src/mmdg.cc
+++ b/src/mmdg.cc
@@ -18,8 +18,7 @@
 
 #include <dune/mmdg/mmdg.hh>
 #include <dune/mmdg/clockhelper.hh>
-#include <dune/mmdg/problems/dgproblem.hh>
-#include <dune/mmdg/problems/poisson.hh>
+#include <dune/mmdg/problems/coupleddgproblem.hh>
 
 
 int main(int argc, char** argv)
@@ -42,7 +41,7 @@ int main(int argc, char** argv)
     using IMapper =
       typename Dune::MultipleCodimMultipleGeomTypeMapper<IGridView>;
     using Coordinate = Dune::FieldVector<double, dim>;
-    using Problem = DGProblem<Coordinate>;
+    using Problem = CoupledDGProblem<Coordinate>;
     using MMDG = MMDG<Grid, GridView, Mapper, IGridView, IMapper, Problem>;
 
     //get parameters from file "parameterDG.ini"
@@ -67,8 +66,8 @@ int main(int argc, char** argv)
 
     //determine problem type
     if (pt["problem"] == "poisson")
-    {
-      problem = new Poisson<dim>();
+    { //TODO
+      problem = new CoupledDGProblem<Coordinate>();
     }
     else
     {