diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 24a4c4b48c4501ded175ee59a29103037833869c..044da9845ae966c3801c3f9f6d71545ccd1ae126 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -105,6 +105,8 @@ private:
       //basis functions (0, phi_iElem,i)
       const int iElemIdxSLE = (dim + 1)*Base::gridView_.size(0) + dim*iElemIdx;
 
+      // === interface entries ===
+
       //fill load vector b with interface entries using a quadrature rule
       for (const auto& ip : rule_qInterface)
       {
@@ -112,17 +114,17 @@ private:
         const auto& qpLocal = ip.position();
         const auto& qpGlobal = iGeo.global(qpLocal);
 
-        const Scalar weight(ip.weight());
-        const Scalar integrationElem(iGeo.integrationElement(qpLocal));
+        const Scalar quadratureFactor =
+          ip.weight() * 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;
+        Base::b[iElemIdxSLE] += quadratureFactor * 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 *
+          Base::b[iElemIdxSLE + i + 1] += quadratureFactor *
             (qpGlobal * iFrame[i]) * qInterfaceEvaluation;
         }
       }
@@ -138,8 +140,8 @@ private:
         const auto& qpLocal = ip.position();
         const auto& qpGlobal = iGeo.global(qpLocal);
 
-        const Scalar weight(ip.weight());
-        const Scalar integrationElem(iGeo.integrationElement(qpLocal));
+        const Scalar quadratureFactor =
+          ip.weight() * iGeo.integrationElement(qpLocal);
 
         for (int i = 0; i < dim-1 ; i++)
         {
@@ -149,7 +151,7 @@ private:
           for (int j = 0; j <= i; j++) //TODO: symmetrize
           {
             Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
-              weight * integrationElem * Base::problem_.aperture(qpGlobal)
+              quadratureFactor * Base::problem_.aperture(qpGlobal)
               * ( K_parDotTau_i * iFrame[j] );
           }
         }
@@ -157,52 +159,220 @@ private:
 
       //TODO: check interface bilinear form
 
+      // === coupling entries ===
+
       //evaluate the coupling term
       // int_iElem beta * phi_iElem,i * phi_iElem,j ds
-      //using a quadrature rule for i = 0 and j = 0,...,dim-1 //TODO or vice versa (symmetrize)
+      //using a quadrature rule for i = 0 and j = 0,...,dim-1 or vice versa //TODO symmetrize more efficiently
       for (const auto& ip : rule_Kperp)
       {
         const auto& qpLocal = ip.position();
         const auto& qpGlobal = iGeo.global(qpLocal);
 
-        const Scalar weight(ip.weight());
-        const Scalar integrationElem(iGeo.integrationElement(qpLocal));
+        const Scalar quadratureFactor =
+          ip.weight() * iGeo.integrationElement(qpLocal);
+        const Scalar betaEvaluation = beta(qpGlobal);
 
         Base::A->entry(iElemIdxSLE, iElemIdxSLE) +=
-          weight * integrationElem * beta(qpGloabal);
+          quadratureFactor * betaEvaluation;
 
         for (int i = 0; i < dim - 1; i++)
         {
-          Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
-            weight * integrationElem * beta(qpGlobal) * (iFrame[i] * qpGlobal)
-            * (iFrame[j] * qpGlobal);
+          Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
+            quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal);
+          Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
+            quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal);
         }
       }
 
       //evaluate the coupling term
       // int_iElem beta * phi_iElem,i * phi_iElem,j ds
-      //using a quadrature rule for i,j = 1,...,dim-1 //TODO symmetrize
+      //using a quadrature rule for i,j = 1,...,dim-1 //TODO symmetrize more efficiently
       for (const auto& ip : rule_KperpPlus1)
       {
         const auto& qpLocal = ip.position();
         const auto& qpGlobal = iGeo.global(qpLocal);
 
-        const Scalar weight(ip.weight());
-        const Scalar integrationElem(iGeo.integrationElement(qpLocal));
+        const Scalar quadratureFactor =
+          ip.weight() * iGeo.integrationElement(qpLocal);
+        const Scalar betaEvaluation = beta(qpGlobal);
 
         for (int i = 0; i < dim - 1; i++)
         {
           for (int j = 0; j <= i; j++)
           {
-            Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
-              weight * integrationElem * beta(qpGlobal) * (iFrame[i] * qpGlobal) ;
+            Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+              quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal)
+              * (iFrame[j] * qpGlobal);
+
+            if (i != j)
+            {
+              Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+                quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal)
+                * (iFrame[j] * qpGlobal);
+            }
+          }
+        }
+      }
+
+      //get neighboring bulk grid elements of iELem
+      const auto elemIn = intersct.inside();
+      const auto elemOut = intersct.outside();
+
+      const int elemInIdxSLE = (dim + 1) * Base::mapper_.index(elemIn);
+      const int elemOutIdxSLE = (dim + 1) * Base::mapper_.index(elemOut);
+
+      //evalute the coupling terms
+      // int_iElem alpha * jump(phi_elem1,0) * jump(phi_elem2,i) ds
+      //and
+      // int_iElem beta * avg(phi_elem1,0) * avg(phi_elem2,i) ds
+      //for elem1, elem2 in {elemIn, elemOut} and i = 0,...,dim
+      //using a quadrature rule
+      for (const auto& ip : rule_Kperp)
+      {
+        const auto& qpLocal = ip.position();
+        const auto& qpGlobal = iGeo.global(qpLocal);
+
+        const Scalar quadratureFactor =
+          ip.weight() * iGeo.integrationElement(qpLocal);
+        const Scalar alphaEvaluation = alpha(qpGlobal);
+        const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
+
+        const Scalar bulkUpdate1 =
+          quadratureFactor * (alphaEvaluation + betaEvaluation4);
+        const Scalar bulkUpdate2 =
+          quadratureFactor * (-alphaEvaluation + betaEvaluation4);
+
+        Base::A->entry(elemInIdxSLE, elemInIdxSLE) += bulkUpdate1;
+        Base::A->entry(elemOutIdxSLE, elemOutIdxSLE) += bulkUpdate1;
+
+        Base::A->entry(elemInIdxSLE, elemOutIdxSLE) += bulkUpdate2;
+        Base::A->entry(elemOutIdxSLE, elemInIdxSLE) += bulkUpdate2;
+
+        for (int i = 0; i < dim; i++)
+        {
+          const Scalar bulkUpdate3 = qpGlobal[i] * bulkUpdate1;
+          const Scalar bulkUpdate4 = qpGlobal[i] * bulkUpdate2;
+
+          Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE) += bulkUpdate3;
+          Base::A->entry(elemInIdxSLE, elemInIdxSLE + i + 1) += bulkUpdate3;
+          Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE) += bulkUpdate3;
+          Base::A->entry(elemOutIdxSLE, elemOutIdxSLE + i + 1) += bulkUpdate3;
+
+          Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE) += bulkUpdate4;
+          Base::A->entry(elemOutIdxSLE, elemInIdxSLE + i + 1) += bulkUpdate4;
+          Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE) += bulkUpdate4;
+          Base::A->entry(elemInIdxSLE, elemOutIdxSLE + i + 1) += bulkUpdate4;
+        }
+      }
+
+      //evalute the coupling terms
+      // int_iElem alpha * jump(phi_elem1,i) * jump(phi_elem2,j) ds
+      //and
+      // int_iElem beta * avg(phi_elem1,i) * avg(phi_elem2,j) ds
+      //for elem1, elem2 in {elemIn, elemOut} and i,j = 1,...,dim
+      //using a quadrature rule
+      for (const auto& ip : rule_KperpPlus1)
+      {
+        const auto& qpLocal = ip.position();
+        const auto& qpGlobal = iGeo.global(qpLocal);
+
+        const Scalar quadratureFactor =
+          ip.weight() * iGeo.integrationElement(qpLocal);
+        const Scalar alphaEvaluation = alpha(qpGlobal);
+        const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
+
+        const Scalar bulkUpdate1 =
+          quadratureFactor * (alphaEvaluation + betaEvaluation4);
+        const Scalar bulkUpdate2 =
+          quadratureFactor * (-alphaEvaluation + betaEvaluation4);
+
+        for (int i = 0; i < dim; i++)
+        {
+          const Scalar bulkUpdate3 = qpGlobal[i] * qpGlobal[i] * bulkUpdate1;
+          const Scalar bulkUpdate4 = qpGlobal[i] * qpGlobal[i] * bulkUpdate2;
+
+          Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE + i + 1) +=
+            bulkUpdate3;
+          Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE + i + 1) +=
+            bulkUpdate3;
+
+          Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE + i + 1) +=
+            bulkUpdate4;
+          Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE + i + 1) +=
+            bulkUpdate4;
+
+          for (int j = 0; j < i; j++)
+          {
+            const Scalar bulkUpdate5 = qpGlobal[i] * qpGlobal[j] * bulkUpdate1;
+            const Scalar bulkUpdate6 = qpGlobal[i] * qpGlobal[j] * bulkUpdate2;
+
+            Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE + j + 1) +=
+              bulkUpdate5;
+            Base::A->entry(elemInIdxSLE + j + 1, elemInIdxSLE + i + 1) +=
+              bulkUpdate5;
+            Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE + j + 1) +=
+              bulkUpdate5;
+            Base::A->entry(elemOutIdxSLE + j + 1, elemOutIdxSLE + i + 1) +=
+              bulkUpdate5;
+
+            Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE + j + 1) +=
+              bulkUpdate6;
+            Base::A->entry(elemOutIdxSLE + j + 1, elemInIdxSLE + i + 1) +=
+              bulkUpdate6;
+            Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE + j + 1) +=
+              bulkUpdate6;
+            Base::A->entry(elemInIdxSLE + j + 1, elemOutIdxSLE + i + 1) +=
+              bulkUpdate6;
           }
         }
       }
 
-      const auto elem1 = intersct.inside();
-      const auto elem2 = intersct.ouside();
-      //TODO coupling entries
+      //evalute the coupling terms
+      // -int_iElem beta * phi_iElem,i * avg(phi_elem,j) ds
+      //for (i = 0 and j = 0,...,dim-1) or (i = 0,...,dim and j = 0)
+      //and for elem in {elemIn, elemOut} using a quadrature rule
+      for (const auto& ip : rule_Kperp) //TODO: symmetrize more efficiently
+      //TODO: apply quadrature rule only once
+      {
+        const auto& qpLocal = ip.position();
+        const auto& qpGlobal = iGeo.global(qpLocal);
+
+        const Scalar quadratureFactor =
+          ip.weight() * iGeo.integrationElement(qpLocal);
+        const Scalar betaEvaluation2 = 0.5 * beta(qpGlobal);
+
+        const Scalar couplingUpdate1 = quadratureFactor * betaEvaluation2;
+
+        Base::A->entry(elemInIdxSLE, iElemIdx) -= couplingUpdate1;
+        Base::A->entry(iElemIdx, elemInIdxSLE) -= couplingUpdate1;
+        Base::A->entry(elemOutIdxSLE, iElemIdx) -= couplingUpdate1;
+        Base::A->entry(iElemIdx, elemOutIdxSLE) -= couplingUpdate1;
+
+        for (int i = 0; i < dim - 1; i++)
+        {
+          const Scalar couplingUpdate2 = qpGlobal[i] * couplingUpdate1;
+          const Scalar couplingUpdate3 =
+            (iFrame[i] * qpGlobal) * couplingUpdate1;
+
+          Base::A->entry(elemInIdxSLE + i + 1, iElemIdx) -= couplingUpdate2;
+          Base::A->entry(iElemIdx, elemInIdxSLE + i + 1) -= couplingUpdate2;
+          Base::A->entry(elemOutIdxSLE + i + 1, iElemIdx) -= couplingUpdate2;
+          Base::A->entry(iElemIdx, elemOutIdxSLE + i + 1) -= couplingUpdate2;
+
+          Base::A->entry(elemInIdxSLE, iElemIdx + i + 1) -= couplingUpdate3;
+          Base::A->entry(iElemIdx + i + 1, elemInIdxSLE) -= couplingUpdate3;
+          Base::A->entry(elemOutIdxSLE, iElemIdx + i + 1) -= couplingUpdate3;
+          Base::A->entry(iElemIdx + i + 1, elemOutIdxSLE) -= couplingUpdate3;
+        }
+
+        //i = dim - 1
+        const Scalar couplingUpdate4 = qpGlobal[dim] * couplingUpdate1;
+        Base::A->entry(elemInIdxSLE + dim + 1, iElemIdx) -= couplingUpdate4;
+        Base::A->entry(iElemIdx, elemInIdxSLE + dim + 1) -= couplingUpdate4;
+        Base::A->entry(elemOutIdxSLE + dim + 1, iElemIdx) -= couplingUpdate4;
+        Base::A->entry(iElemIdx, elemOutIdxSLE + dim + 1) -= couplingUpdate4;
+      }
     }
   }