diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 4ef53cef1371c38ab595dfa88f5d10ed9c1a8722..7db35c580d0f28fde23a5355a6035d4a7a7ebaf1 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -46,6 +46,8 @@ public:
 
     (*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);
@@ -471,8 +473,6 @@ private:
             (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
               centerI * (interfaceUpdate1 - interfaceUpdate2);
 
-            std::cout << iElemIdxSLE + i + 1 << "\t" << interfaceUpdate3 << "\t" << centerI * (interfaceUpdate1 - interfaceUpdate2);
-
             for (int j = 0; j < i; i++)
             { //NOTE: only relevant for n=3, requires quadrature rule for n=3
               Coordinate KparDotTau_j;
@@ -593,8 +593,6 @@ private:
             (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
               centerI * (interfaceUpdate1 - 2.0 * interfaceUpdate2);
 
-            std::cout << iElemIdxSLE + i + 1 << "\t" << iFrame[i] << "\t" << center * iFrame[i] << "\t" << interfaceUpdate3 << "\t" << centerI * (interfaceUpdate1 - 2.0 * interfaceUpdate2) << "\n\n";
-
             for (int j = 0; j < i; i++)
             { //NOTE: only relevant for n=3, requires quadrature rule for n=3
               Coordinate KparDotTau_j;
@@ -670,7 +668,8 @@ private:
 
     for (const auto& iElem : elements(iGridView_))
     {
-      const int iElemIdxSLE = bulkDOF + dim * iMapper_.index(iElem);
+      const int iElemIdxOutput = dim * iMapper_.index(iElem);
+      const int iElemIdxSLE = bulkDOF + iElemIdxOutput;
       const auto& iGeo = iElem.geometry();
 
       //get basis of the interface coordinate system for evaluation of
@@ -682,21 +681,21 @@ private:
       {
         if (Base::problem_.hasExactSolution())
         {
-          exactIPressure[iElemIdxSLE + k] =
+          exactIPressure[iElemIdxOutput + k] =
             Base::problem_.exactInterfaceSolution(iGeo.corner(k));
         }
 
         //contribution of the basis function
         // phi_iElem,0 (x) = indicator(iElem);
         //at the kth corner of iElem
-        iPressure[iElemIdxSLE + k] = Base::d[iElemIdxSLE];
+        iPressure[iElemIdxOutput + k] = Base::d[iElemIdxSLE];
 
         for (int i = 0; i < dim - 1; i++)
         {
           //contribution of the basis function
           // phi_iElem,i (x) = (x * tau[i]) * indicator(iElem);
           //at the kth corner of iElem
-          iPressure[iElemIdxSLE + k] +=
+          iPressure[iElemIdxOutput + k] +=
             Base::d[iElemIdxSLE + i + 1] * (iGeo.corner(k) * iFrame[i]);
         }
       }