diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index f0f3495f52760ba457e04954462a5f2bc8ee7d21..344c397e7575d645a5d17096a14080e5d0b9b3d2 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -26,8 +26,12 @@ public:
   static constexpr auto edge = Dune::GeometryTypes::simplex(dim-1);
 
   using Scalar = typename GridView::ctype;
-  using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
-  using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
+  // using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
+  // using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
+
+  using Matrix = Dune::DynamicMatrix<Scalar>;
+  using Vector = Dune::DynamicVector<Scalar>;
+
   using VTKFunction = Dune::VTKFunction<GridView>;
   using P1Function = Dune::NonconformingP1VTKFunction<GridView,
     Dune::DynamicVector<Scalar>>;
@@ -41,21 +45,34 @@ public:
     gridView_(gridView), mapper_(mapper), problem_(problem),
     dof_((1 + dim) * gridView.size(0))
     {
+
+      A = std::make_shared<Matrix>(dof_, dof_, 0.0); //stiffness matrix
+      b = Vector(dof_, 0.0); //load vector
+      d = Vector(dof_, 0.0); //solution vector
+
       //initialize stiffness matrix A, load vector b and solution vector d
-      A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
-      b = Vector(dof_);
-      d = Vector(dof_);
+      // A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
+      // b = Vector(dof_);
+      // d = Vector(dof_);
     }
 
   void operator() (const Scalar mu)
   {
     //assemble stiffness matrix A and load vector b
     assembleSLE(mu);
-    A->compress();
+    //A->compress();
 
-    Dune::InverseOperatorResult result;
-    Dune::UMFPack<Matrix> solver(*A);
-    solver.apply(d, b, result);
+    (*A).solve(d, b);
+
+    for (int i = 0; i < dof_; i++)
+      for (int j = 0; j < i; j++)
+        /*assert*/if(std::abs((*A)[i][j] - (*A)[j][i]) >
+          std::numeric_limits<Scalar>::epsilon())
+            std::cout << i << ", " << j << std::endl;
+
+    // Dune::InverseOperatorResult result;
+    // Dune::UMFPack<Matrix> solver(*A);
+    // solver.apply(d, b, result);
 
     //write solution to a vtk file
     writeVTKoutput();
@@ -115,10 +132,14 @@ protected:
     const int dof) :
     gridView_(gridView), mapper_(mapper), problem_(problem), dof_(dof)
     {
+      A = std::make_shared<Matrix>(dof_, dof_, 0.0); //stiffness matrix
+      b = Vector(dof_, 0.0); //load vector
+      d = Vector(dof_, 0.0); //solution vector
+
       //initialize stiffness matrix A, load vector b and solution vector d
-      A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
-      b = Vector(dof_);
-      d = Vector(dof_);
+      // A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
+      // b = Vector(dof_);
+      // d = Vector(dof_);
     }
 
   //assemble stiffness matrix A and load vector b
@@ -190,7 +211,7 @@ protected:
         {
           for (int j = 0; j <= i; j++)
           {
-            A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
+            (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
               quadratureFactor * problem_.K(qpGlobal)[j][i];
           }
         }
@@ -303,7 +324,7 @@ protected:
 
         //exact evaluation of
         // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
-        A->entry(elemIdxSLE, elemIdxSLE) += mu * intersctVol;
+        (*A)[elemIdxSLE][elemIdxSLE] += mu * intersctVol;
 
         if (intersct.neighbor()) //intersct has neighboring element
         {
@@ -318,7 +339,7 @@ protected:
             //and
             // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds
             // = 0.5 * int_intersct K[:,i] * normal ds
-            A->entry(elemIdxSLE + i + 1, elemIdxSLE) +=
+            (*A)[elemIdxSLE + i + 1][elemIdxSLE] +=
               mu * linearIntegrals[i] - 0.5 * KIntegrals[i];
 
             for (int j = 0; j <= i; j++)
@@ -329,7 +350,7 @@ protected:
               //and
               // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,j) ds
               // = 0.5 * int_intersct x[j] * ( K[:,i] * normal ) ds
-              A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1)
+              (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
                 += mu * quadraticIntregrals[i][j]
                   - 0.5 * ( KIntegralsX[i][j] + KIntegralsX[j][i] );
             }
@@ -342,11 +363,11 @@ protected:
 
           //exact evaluation of
           // int_intersct mu * jump(phi_elem,0) * jump(phi_neighbor,0) ds
-          A->entry(elemIdxSLE, neighborIdxSLE) += -mu * intersctVol;
+          (*A)[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;
 
           //stiffness matrix A is symmetric
-          A->entry(neighborIdxSLE, elemIdxSLE) +=
-            A->entry(elemIdxSLE, neighborIdxSLE);
+          (*A)[neighborIdxSLE][elemIdxSLE] +=
+            (*A)[elemIdxSLE][neighborIdxSLE];
 
           for (int i = 0; i < dim; i++)
           {
@@ -356,7 +377,7 @@ protected:
             //and
             // int_intersct avg(K*grad(phi_neighbor,i)) * jump(phi_elem,0) ds
             // = -0.5 * int_intersct K[:,i] * normal ds
-            A->entry(elemIdxSLE + i + 1, neighborIdxSLE) +=
+            (*A)[elemIdxSLE + i + 1][neighborIdxSLE] +=
               -mu * linearIntegrals[i] + 0.5 * KIntegrals[i];
 
             //we use the relations
@@ -365,14 +386,14 @@ protected:
             //and
             // int_intersct avg(K*grad(phi_neighbor,i)) * jump(phi_elem,0) ds
             // = 0.5 * int_intersct K[:,i] * normal ds
-            A->entry(elemIdxSLE, neighborIdxSLE + i + 1) +=
+            (*A)[elemIdxSLE][neighborIdxSLE + i + 1] +=
               -mu * linearIntegrals[i] - 0.5 * KIntegrals[i];
 
             //stiffness matrix A is symmetric
-            A->entry(neighborIdxSLE, elemIdxSLE + i + 1) +=
-              A->entry(elemIdxSLE + i + 1, neighborIdxSLE);
-            A->entry(neighborIdxSLE + i + 1, elemIdxSLE) +=
-              A->entry(elemIdxSLE, neighborIdxSLE + i + 1);
+            (*A)[neighborIdxSLE][elemIdxSLE + i + 1] +=
+              (*A)[elemIdxSLE + i + 1][neighborIdxSLE];
+            (*A)[neighborIdxSLE + i + 1][elemIdxSLE] +=
+              (*A)[elemIdxSLE][neighborIdxSLE + i + 1];
 
             for (int j = 0; j <= i; j++)
             {
@@ -385,13 +406,13 @@ protected:
               //as well as
               // int_intersct avg(K*grad(phi_elem,i)) * jump(phi_neighbor,j) ds
               // = -0.5 * int_intersct x[j] * ( K[:,i] * normal ) ds
-              A->entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1) +=
+              (*A)[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
                 -mu * quadraticIntregrals[i][j]
                 - 0.5 * ( KIntegralsX[j][i] - KIntegralsX[i][j] );
 
               //stiffness matrix A is symmetric
-              A->entry(neighborIdxSLE + j + 1, elemIdxSLE + i + 1) +=
-                A->entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1);
+              (*A)[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
+                (*A)[elemIdxSLE + i + 1][neighborIdxSLE + j + 1];
 
               if (i != j)
               {
@@ -404,13 +425,13 @@ protected:
                 //as well as
                 // int_intersct avg(K*grad(phi_elem,j))*jump(phi_neighbor,i) ds
                 // = -0.5 * int_intersct x[i] * ( K[:,j] * normal ) ds
-                A->entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1) +=
+                (*A)[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] +=
                   -mu * quadraticIntregrals[i][j]
                   - 0.5 * ( KIntegralsX[i][j] - KIntegralsX[j][i] );
 
                 //stiffness matrix A is symmetric
-                A->entry(neighborIdxSLE + i + 1, elemIdxSLE + j + 1) +=
-                  A->entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1);
+                (*A)[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] +=
+                  (*A)[elemIdxSLE + j + 1][neighborIdxSLE + i + 1];
               }
             }
           }
@@ -452,7 +473,7 @@ protected:
             //and
             // int_intersct phi_elem,0 * K*grad(phi_elem,i) * normal ds
             //with K*grad(phi_elem,i) = K[:,i]
-            A->entry(elemIdxSLE + i + 1, elemIdxSLE) +=
+            (*A)[elemIdxSLE + i + 1][elemIdxSLE] +=
               mu * linearIntegrals[i] - KIntegrals[i];
 
             for (int j = 0; j <= i; j++)
@@ -463,7 +484,7 @@ protected:
                 // int_intersct (phi_elem,i * K[:,j]
                 //    + phi_elem,j * K[:,i]) * normal ds
                 //with K*grad(phi_elem,i) = K[:,i]
-              A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
+              (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
                 mu * quadraticIntregrals[i][j]
                 - ( KIntegralsX[i][j] + KIntegralsX[j][i] );
             }
@@ -474,13 +495,13 @@ protected:
       //stiffness matrix A is symmetric
       for (int i = 0; i < dim; i++)
       {
-        A->entry(elemIdxSLE, elemIdxSLE + i + 1) =
-          A->entry(elemIdxSLE + i + 1, elemIdxSLE);
+        (*A)[elemIdxSLE][elemIdxSLE + i + 1] =
+          (*A)[elemIdxSLE + i + 1][elemIdxSLE];
 
         for(int j = 0; j < i; j++)
         {
-          A->entry(elemIdxSLE + j + 1, elemIdxSLE + i + 1) =
-            A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1);
+          (*A)[elemIdxSLE + j + 1][elemIdxSLE + i + 1] =
+            (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1];
         }
       }
     }
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index ccaef4d5945090dba159bf9a51f19f741a2d4813..c830287cfacdf8542b50acebc5467ac614e7a03e 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -34,11 +34,21 @@ public:
   const void operator() (const Scalar mu, const Scalar xi)
   {
     assembleSLE(mu, xi);
-    Base::A->compress();
+    //Base::A->compress();
 
-    Dune::InverseOperatorResult result;
-    Dune::UMFPack<Matrix> solver(*Base::A);
-    solver.apply(Base::d, Base::b, result);
+    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 << i << ", " << j << std::endl;
+
+    std::cout << *Base::A << std::endl;
+
+    (*Base::A).solve(Base::d, Base::b);
+
+    // Dune::InverseOperatorResult result;
+    // Dune::UMFPack<Matrix> solver(*Base::A);
+    // solver.apply(Base::d, Base::b, result);
 
     //write solution to a vtk file
     writeVTKoutput();
@@ -96,15 +106,10 @@ private:
   //assemble stiffness matrix A and load vector b
   void assembleSLE (const Scalar mu, const Scalar xi)
   {
-    //modell parameters //NOTE: what should be the capture of alpha?
-    const auto alpha = [this](const Coordinate& pos) -> Scalar
-    {
-      return Base::problem_.Kperp(pos) / Base::problem_.aperture(pos);
-    };
-
-    const auto beta = [&alpha, &xi](const Coordinate& pos) -> Scalar
+    //modell parameter
+    const auto beta = [this, &xi](const Coordinate& pos) -> Scalar
     {
-      return 4 * alpha(pos) / (2 * xi - 1);
+      return 4 * Base::problem_.Kperp(pos) / (2 * xi - 1);
     };
 
     //quadrature rules
@@ -157,7 +162,7 @@ private:
       //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)
+      //basis functions (0, phi_iElem,i)interfaceinterface
       const int iElemIdxSLE = (dim + 1)*Base::gridView_.size(0) + dim*iElemIdx;
 
       // === coupling entries ===
@@ -174,14 +179,14 @@ private:
           ip.weight() * iGeo.integrationElement(qpLocal);
         const Scalar betaEvaluation = beta(qpGlobal);
 
-        Base::A->entry(iElemIdxSLE, iElemIdxSLE) +=
+        (*Base::A)[iElemIdxSLE][iElemIdxSLE] +=
           quadratureFactor * betaEvaluation;
 
         for (int i = 0; i < dim - 1; i++)
         {
-          Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
+          (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
             quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal);
-          Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
+          (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] +=
             quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal);
         }
       }
@@ -202,13 +207,13 @@ private:
         {
           for (int j = 0; j <= i; j++)
           {
-            Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+            (*Base::A)[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) +=
+              (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
                 quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal)
                 * (iFrame[j] * qpGlobal);
             }
@@ -224,7 +229,7 @@ private:
       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
+      // int_iElem K_perp * 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
@@ -236,39 +241,39 @@ private:
 
         const Scalar quadratureFactor =
           ip.weight() * iGeo.integrationElement(qpLocal);
-        const Scalar alphaEvaluation = alpha(qpGlobal);
+        const Scalar KEvaluation = Base::problem_.Kperp(qpGlobal);
         const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
 
         const Scalar bulkUpdate1 =
-          quadratureFactor * (alphaEvaluation + betaEvaluation4);
+          quadratureFactor * (KEvaluation + betaEvaluation4);
         const Scalar bulkUpdate2 =
-          quadratureFactor * (-alphaEvaluation + betaEvaluation4);
+          quadratureFactor * (-KEvaluation + betaEvaluation4);
 
-        Base::A->entry(elemInIdxSLE, elemInIdxSLE) += bulkUpdate1;
-        Base::A->entry(elemOutIdxSLE, elemOutIdxSLE) += bulkUpdate1;
+        (*Base::A)[elemInIdxSLE][elemInIdxSLE] += bulkUpdate1;
+        (*Base::A)[elemOutIdxSLE][elemOutIdxSLE] += bulkUpdate1;
 
-        Base::A->entry(elemInIdxSLE, elemOutIdxSLE) += bulkUpdate2;
-        Base::A->entry(elemOutIdxSLE, elemInIdxSLE) += bulkUpdate2;
+        (*Base::A)[elemInIdxSLE][elemOutIdxSLE] += bulkUpdate2;
+        (*Base::A)[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)[elemInIdxSLE + i + 1][elemInIdxSLE] += bulkUpdate3;
+          (*Base::A)[elemInIdxSLE][elemInIdxSLE + i + 1] += bulkUpdate3;
+          (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE] += bulkUpdate3;
+          (*Base::A)[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;
+          (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE] += bulkUpdate4;
+          (*Base::A)[elemOutIdxSLE][elemInIdxSLE + i + 1] += bulkUpdate4;
+          (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE] += bulkUpdate4;
+          (*Base::A)[elemInIdxSLE][elemOutIdxSLE + i + 1] += bulkUpdate4;
         }
       }
 
       //evalute the coupling terms
-      // int_iElem alpha * jump(phi_elem1,i) * jump(phi_elem2,j) ds
+      // int_iElem K_perp * 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
@@ -280,27 +285,27 @@ private:
 
         const Scalar quadratureFactor =
           ip.weight() * iGeo.integrationElement(qpLocal);
-        const Scalar alphaEvaluation = alpha(qpGlobal);
+        const Scalar KEvaluation = Base::problem_.Kperp(qpGlobal);
         const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
 
         const Scalar bulkUpdate1 =
-          quadratureFactor * (alphaEvaluation + betaEvaluation4);
+          quadratureFactor * (KEvaluation + betaEvaluation4);
         const Scalar bulkUpdate2 =
-          quadratureFactor * (-alphaEvaluation + betaEvaluation4);
+          quadratureFactor * (-KEvaluation + 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) +=
+          (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE + i + 1] +=
             bulkUpdate3;
-          Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE + i + 1) +=
+          (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE + i + 1] +=
             bulkUpdate3;
 
-          Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE + i + 1) +=
+          (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + i + 1] +=
             bulkUpdate4;
-          Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE + i + 1) +=
+          (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE + i + 1] +=
             bulkUpdate4;
 
           for (int j = 0; j < i; j++)
@@ -308,22 +313,22 @@ private:
             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) +=
+            (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE + j + 1] +=
               bulkUpdate5;
-            Base::A->entry(elemInIdxSLE + j + 1, elemInIdxSLE + i + 1) +=
+            (*Base::A)[elemInIdxSLE + j + 1][elemInIdxSLE + i + 1] +=
               bulkUpdate5;
-            Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE + j + 1) +=
+            (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE + j + 1] +=
               bulkUpdate5;
-            Base::A->entry(elemOutIdxSLE + j + 1, elemOutIdxSLE + i + 1) +=
+            (*Base::A)[elemOutIdxSLE + j + 1][elemOutIdxSLE + i + 1] +=
               bulkUpdate5;
 
-            Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE + j + 1) +=
+            (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + j + 1] +=
               bulkUpdate6;
-            Base::A->entry(elemOutIdxSLE + j + 1, elemInIdxSLE + i + 1) +=
+            (*Base::A)[elemOutIdxSLE + j + 1][elemInIdxSLE + i + 1] +=
               bulkUpdate6;
-            Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE + j + 1) +=
+            (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE + j + 1] +=
               bulkUpdate6;
-            Base::A->entry(elemInIdxSLE + j + 1, elemOutIdxSLE + i + 1) +=
+            (*Base::A)[elemInIdxSLE + j + 1][elemOutIdxSLE + i + 1] +=
               bulkUpdate6;
           }
         }
@@ -345,10 +350,10 @@ private:
 
         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;
+        (*Base::A)[elemInIdxSLE][iElemIdx] -= couplingUpdate1;
+        (*Base::A)[iElemIdx][elemInIdxSLE] -= couplingUpdate1;
+        (*Base::A)[elemOutIdxSLE][iElemIdx] -= couplingUpdate1;
+        (*Base::A)[iElemIdx][elemOutIdxSLE] -= couplingUpdate1;
 
         for (int i = 0; i < dim - 1; i++)
         {
@@ -356,23 +361,23 @@ private:
           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)[elemInIdxSLE + i + 1][iElemIdx] -= couplingUpdate2;
+          (*Base::A)[iElemIdx][elemInIdxSLE + i + 1] -= couplingUpdate2;
+          (*Base::A)[elemOutIdxSLE + i + 1][iElemIdx] -= couplingUpdate2;
+          (*Base::A)[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;
+          (*Base::A)[elemInIdxSLE][iElemIdx + i + 1] -= couplingUpdate3;
+          (*Base::A)[iElemIdx + i + 1][elemInIdxSLE] -= couplingUpdate3;
+          (*Base::A)[elemOutIdxSLE][iElemIdx + i + 1] -= couplingUpdate3;
+          (*Base::A)[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;
+        (*Base::A)[elemInIdxSLE + dim + 1][iElemIdx] -= couplingUpdate4;
+        (*Base::A)[iElemIdx][elemInIdxSLE + dim + 1] -= couplingUpdate4;
+        (*Base::A)[elemOutIdxSLE + dim + 1][iElemIdx] -= couplingUpdate4;
+        (*Base::A)[iElemIdx][elemOutIdxSLE + dim + 1] -= couplingUpdate4;
       }
 
       // === interface entries ===
@@ -418,12 +423,12 @@ private:
 
           for (int j = 0; j <= i; j++)
           {
-            Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+            (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
               quadratureFactor * ( KparDotTau_i * iFrame[j] );
 
             if (i != j)
             {
-              Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+              (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
                 quadratureFactor * ( KparDotTau_i * iFrame[j] );
             }
           }
@@ -440,7 +445,7 @@ private:
 
         //evaluate the interface term
         // int_intersct mu^Gamma * jump(phi_iElem,0) * jump(phi_iElem,0) dr
-        Base::A->entry(iElemIdxSLE, iElemIdxSLE) += mu;
+        (*Base::A)[iElemIdxSLE][iElemIdxSLE] += mu;
 
         if (intersct.neighbor()) //inner interface edge
         {
@@ -459,13 +464,15 @@ private:
             const Scalar interfaceUpdate3 =
               interfaceUpdate1 - 0.5 * interfaceUpdate2;
 
-            Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
+            (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
               interfaceUpdate3;
-            Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
+            (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] +=
               interfaceUpdate3;
-            Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + i + 1) +=
+            (*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;
@@ -475,9 +482,9 @@ private:
                 - 0.5 * (interfaceUpdate2 * (center * iFrame[j])
                 + (KparDotTau_j * normal) * centerI);
 
-              Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+              (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
                 interfaceUpdate4;
-              Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+              (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
                 interfaceUpdate4;
             }
           }
@@ -502,8 +509,8 @@ private:
           // -int_intersct avg(Kparallel * grad(phi_neighbor,i))
           //    * jump(phi_iElem,j) dr
           //for i,j = 0,...,dim-1
-          Base::A->entry(iElemIdxSLE, neighborIdxSLE) += -mu;
-          Base::A->entry(neighborIdxSLE, iElemIdxSLE) += -mu;
+          (*Base::A)[iElemIdxSLE][neighborIdxSLE] += -mu;
+          (*Base::A)[neighborIdxSLE][iElemIdxSLE] += -mu;
 
           for (int i = 0; i < dim - 1; i++)
           {
@@ -516,16 +523,16 @@ private:
             const Scalar interfaceUpdate4 =
               -interfaceUpdate1 - 0.5 * interfaceUpdate2;
 
-            Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE) +=
+            (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE] +=
               interfaceUpdate3;
-            Base::A->entry(neighborIdxSLE, iElemIdxSLE + i + 1) +=
+            (*Base::A)[neighborIdxSLE][iElemIdxSLE + i + 1] +=
               interfaceUpdate3;
-            Base::A->entry(neighborIdxSLE + i + 1, iElemIdxSLE) +=
+            (*Base::A)[neighborIdxSLE + i + 1][iElemIdxSLE] +=
               -interfaceUpdate4;
-            Base::A->entry(iElemIdxSLE, neighborIdxSLE + i + 1) +=
+            (*Base::A)[iElemIdxSLE][neighborIdxSLE + i + 1] +=
               -interfaceUpdate4;
 
-            Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE + i + 1) +=
+            (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE + i + 1] +=
               -centerI * interfaceUpdate1;
 
             for (int j = 0; j < i; j++)
@@ -541,13 +548,13 @@ private:
               const Scalar interfaceUpdate8 =
                 -interfaceUpdate5 + 0.5 * interfaceUpdate6;
 
-              Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE + j + 1) +=
+              (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
                 interfaceUpdate7;
-              Base::A->entry(neighborIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+              (*Base::A)[neighborIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
                 interfaceUpdate7;
-              Base::A->entry(iElemIdxSLE + j + 1, neighborIdxSLE + i + 1) +=
+              (*Base::A)[iElemIdxSLE + j + 1][neighborIdxSLE + i + 1] +=
                 interfaceUpdate8;
-              Base::A->entry(neighborIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+              (*Base::A)[neighborIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
                 interfaceUpdate8;
             }
           }
@@ -579,13 +586,15 @@ private:
 
             Base::b[iElemIdxSLE] += dgGamma * (mu * centerI - interfaceUpdate2);
 
-            Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
+            (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
               interfaceUpdate3;
-            Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
+            (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] +=
               interfaceUpdate3;
-            Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + i + 1) +=
+            (*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;
@@ -595,9 +604,9 @@ private:
                 - (interfaceUpdate2 * (center * iFrame[j])
                   + (KparDotTau_j * normal) * centerI);
 
-              Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+              (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
                 interfaceUpdate4;
-              Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+              (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
                 interfaceUpdate4;
             }
           }
@@ -620,7 +629,7 @@ private:
     for (int i = 0; i < dim - 1; i++)
     {
       Coordinate tau(0.0);
-      for (int j = 0; j < dim - 1; j++)
+      for (int j = 0; j < dim; j++)
       {
         if (j == i)
         {
diff --git a/src/dg.cc b/src/dg.cc
index 744b48cede991d649b5568e08eff1cceb0ff58ab..d2bba7c8ae1d818bd11e0908969cbeb332c1b8cf 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -39,7 +39,9 @@ int main(int argc, char** argv)
     //use YaspGrid if dim = 1 and MMesh otherwise
     using Grid = std::conditional<dim == 1, Dune::YaspGrid<dim>,
       Dune::MovingMesh<dim>>::type;
-    using GridFactory = Dune::DGFGridFactory<Grid>;
+    //using GridFactory = Dune::DGFGridFactory<Grid>;
+    using GridFactory = Dune::GmshGridFactory<Grid, /*implicit=*/false>;
+
     using GridView = typename Grid::LeafGridView;
     using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
     using Coordinate = Dune::FieldVector<double, dim>;
@@ -54,7 +56,9 @@ int main(int argc, char** argv)
     const double mu = std::stod(pt["mu"]);
 
     //create a grid from .dgf file
-    GridFactory gridFactory( "grids/cube" + std::to_string(dim) + "d.dgf" );
+    //GridFactory gridFactory( "grids/cube" + std::to_string(dim) + "d.dgf" );
+    GridFactory gridFactory( "grids/mmdg2.msh" );
+
     const Grid& grid = *gridFactory.grid();
     const GridView& gridView = grid.leafGridView();