diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index e5b4adbb0a8c0524605cb99e9e83ff3f4f3ea16c..a2959980416278b316410c6373886b2eb9ee1684 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -1,6 +1,13 @@
 #ifndef DUNE_MMDG_DG_HH
 #define DUNE_MMDG_DG_HH
 
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/dynvector.hh>
+
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
 #include <dune/mmdg/nonconformingp1vtkfunction.hh>
 
 template<class GridView, class Mapper, class Problem>
@@ -38,6 +45,12 @@ public:
       const auto& geo = elem.geometry();
       const double elemVol = geo.volume();
 
+      //in the system of linear equations (SLE) Ad = b,
+      //the index elemIdxSLE refers to the basis function phi_elem,0
+      //and the indices elemIdxSLE + i + 1, i = 1,...,dim, refer to the
+      //basis function phi_elem,i
+      const int elemIdxSLE = (dim + 1)*elemIdx;
+
   /*    //TODO: can be done outside of the loop?
       const Dune::QuadratureRule<double,dim>& rule =
         Dune::QuadratureRules<double,dim>::rule(geo.type(),order);
@@ -51,29 +64,28 @@ public:
         const auto weight = ip.weight();
 
         //quadrature for int_elem q*phi_elem,0 dV
-        b[(1 + dim)*elemIdx] += weight * q(qp);
+        b[elemIdxSLE] += weight * q(qp);
 
         //quadrature for int_elem q*phi_elem,i dV
         for (int i = 0; i < dim; i++)
         {
-          b[(1 + dim)*elemIdx + i + 1] += weight * qp[i] * q(qp);
+          b[elemIdxSLE + i + 1] += weight * qp[i] * q(qp);
         }
       }
   */
 
       //NOTE: makeshift solution for source term q = -1
-      b[(1 + dim)*elemIdx] = -elemVol;
+      b[elemIdxSLE] = -elemVol;
       for (int i = 0; i < dim; i++)
       {
-        b[(1 + dim)*elemIdx + i + 1] = -elemVol * geo.center()[i];
+        b[elemIdxSLE + i + 1] = -elemVol * geo.center()[i];
       }
 
       for (int i = 0; i < dim; i++)
       {
         //exact evaluation of
         // int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV
-        A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + i + 1] =
-          K * elemVol;
+        A[elemIdxSLE + i + 1][elemIdxSLE + i + 1] = K * elemVol;
       }
 
       //iterate over all intersection with the boundary of elem
@@ -122,18 +134,20 @@ public:
                 ip.weight() * qp[i] * qp[j] * intersctVol;
             }
 */
+            //NOTE: probably unnecessary
+            quadraticIntregrals[i][j] = quadraticIntregrals[j][i];
           }
         }
 
         //exact evaluation of
         // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
-        A[(1 + dim)*elemIdx][(1 + dim)*elemIdx] +=
-          mu * intersctGeo.volume();
+        A[elemIdxSLE][elemIdxSLE] += mu * intersctGeo.volume();
 
         if (intersct.neighbor()) //intersct has neighboring element
         {
           //index of the neighboring element
           const int neighborIdx = mapper_.index(intersct.outside());
+          const int neighborIdxSLE = (dim + 1)*neighborIdx;
 
           for (int i = 0; i < dim; i++)
           { //we use the relations
@@ -143,13 +157,8 @@ public:
             // int_intersct avg(K*grad(phi_elem,i))
             //   *jump(phi_elem,0) ds
             // = 0.5 * K * normal[i] * vol(intersct)
-            A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] +=
-              mu * linearIntegrals[i]
-              - 0.5 * K * normal[i] * intersctVol;
-
-            //stiffness matrix A is symmetric
-            A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] +=
-              A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx];
+            A[elemIdxSLE + i + 1][elemIdxSLE] +=
+              mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
 
             for (int j = 0; j <= i; j++)
             {
@@ -160,15 +169,10 @@ public:
               // int_intersct avg(K*grad(phi_elem,i))
               //   *jump(phi_elem,j) ds
               // = 0.5 * K * normal[i] * int_intersct x_j ds
-              A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1]
+              A[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
                 += mu * quadraticIntregrals[i][j]
                   - 0.5 * K * (normal[i] * linearIntegrals[j]
                     + normal[j] * linearIntegrals[i]);
-
-              //stiffness matrix A is symmetric
-              A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1]
-                += A[(1 + dim)*elemIdx + i + 1]
-                  [(1 + dim)*elemIdx + j + 1];
             }
           }
 
@@ -179,12 +183,10 @@ public:
 
           //exact evaluation of
           // int_intersct mu*jump(phi_elem,0)*jump(phi_neighbor,0) ds
-          A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx] +=
-            -mu * intersctVol;
+          A[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;
 
           //stiffness matrix A is symmetric
-          A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx] +=
-            A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx];
+          A[neighborIdxSLE][elemIdxSLE] += A[elemIdxSLE][neighborIdxSLE];
 
           for (int i = 0; i < dim; i++)
           {
@@ -196,9 +198,8 @@ public:
             // int_intersct avg(K*grad(phi_neighbor,i))
             //  *jump(phi_elem,0) ds
             // = 0.5 * K * normal[i] * vol(intersct)
-            A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx] +=
-              -mu * linearIntegrals[i]
-              - 0.5 * K * normal[i] * intersctVol;
+            A[elemIdxSLE + i + 1][neighborIdxSLE] +=
+              -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
 
             //we use the relations
             // int_intersct mu*jump(phi_elem,0)
@@ -208,15 +209,14 @@ public:
             // int_intersct avg(K*grad(phi_neighbor,i))
             //  *jump(phi_elem,0) ds
             // = 0.5 * K * normal[i] * vol(intersct)
-            A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + i + 1] +=
-              -mu * linearIntegrals[i]
-              - 0.5 * K * normal[i] * intersctVol;
+            A[elemIdxSLE][neighborIdxSLE + i + 1] +=
+              -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
 
             //stiffness matrix A is symmetric
-            A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx + i + 1] +=
-              A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx];
-            A[(1 + dim)*neighborIdx + i + 1][(1 + dim)*elemIdx] +=
-              A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + 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++)
             {
@@ -232,38 +232,37 @@ public:
               // int_intersct avg(K*grad(phi_elem,i))
               //  *jump(phi_neighbor,j) ds
               // = -0.5 * K * normal[i] * int_intersct x_j ds
-              A[(1 + dim)*elemIdx + i + 1]
-                [(1 + dim)*neighborIdx + j + 1] +=
+              A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
                 -mu * quadraticIntregrals[i][j]
                 - 0.5 * K * (normal[j] * linearIntegrals[i]
                   - normal[i] * linearIntegrals[j]);
 
-              // int_intersct mu*jump(phi_elem,j)
-              //  *jump(phi_neighbor,i) ds
-              // = -mu*int_intersct x_i*x_j ds
-              //and
-              // int_intersct avg(K*grad(phi_neighbor,i))
-              //  *jump(phi_elem,j) ds
-              // = 0.5 * K * normal[i] * int_intersct x_j ds
-              //as well as
-              // int_intersct avg(K*grad(phi_elem,j))
-              //  *jump(phi_neighbor,i) ds
-              // = -0.5 * K * normal[j] * int_intersct x_i ds
-              A[(1 + dim)*elemIdx + j + 1]
-                [(1 + dim)*neighborIdx + i + 1] +=
-                -mu * quadraticIntregrals[i][j]
-                - 0.5 * K * (normal[i] * linearIntegrals[j]
-                  - normal[j] * linearIntegrals[i]);
-
               //stiffness matrix A is symmetric
-              A[(1 + dim)*neighborIdx + j + 1]
-                [(1 + dim)*elemIdx + i + 1] +=
-                A[(1 + dim)*elemIdx + i + 1]
-                  [(1 + dim)*neighborIdx + j + 1];
-              A[(1 + dim)*neighborIdx + i + 1]
-                [(1 + dim)*elemIdx + j + 1] +=
-                A[(1 + dim)*elemIdx + j + 1]
-                  [(1 + dim)*neighborIdx + i + 1];
+              A[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
+                A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1];
+
+              if (i != j)
+              {
+                // int_intersct mu*jump(phi_elem,j)
+                //  *jump(phi_neighbor,i) ds
+                // = -mu*int_intersct x_i*x_j ds
+                //and
+                // int_intersct avg(K*grad(phi_neighbor,i))
+                //  *jump(phi_elem,j) ds
+                // = 0.5 * K * normal[i] * int_intersct x_j ds
+                //as well as
+                // int_intersct avg(K*grad(phi_elem,j))
+                //  *jump(phi_neighbor,i) ds
+                // = -0.5 * K * normal[j] * int_intersct x_i ds
+                A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] +=
+                  -mu * quadraticIntregrals[i][j]
+                  - 0.5 * K * (normal[i] * linearIntegrals[j]
+                    - normal[j] * linearIntegrals[i]);
+
+                //stiffness matrix A is symmetric
+                A[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] +=
+                  A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1];
+              }
             }
           }
         }
@@ -277,13 +276,8 @@ public:
             // int_intersct avg(K*grad(phi_elem,i))
             //  *jump(phi_elem,0) ds
             // = K * normal[i] * vol(intersct)
-            A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] +=
-              mu * linearIntegrals[i]
-              - 0.5 * K * normal[i] * intersctVol;
-
-            //stiffness matrix A is symmetric
-            A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] +=
-              A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx];
+            A[elemIdxSLE + i + 1][elemIdxSLE] +=
+              mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
 
             for (int j = 0; j <= i; j++)
             {
@@ -294,20 +288,26 @@ public:
               // int_intersct avg(K*grad(phi_elem,i))
               //  *jump(phi_elem,j) ds
               // = 0.5 * K * normal[i] * int_intersct x_j ds
-              A[(1 + dim)*elemIdx + i + 1]
-                [(1 + dim)*elemIdx + j + 1] +=
+              A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
                 mu * quadraticIntregrals[i][j]
                 - 0.5 * K * (normal[i] * linearIntegrals[j]
                   + normal[j] * linearIntegrals[i]);
-
-              //stiffness matrix A is symmetric
-              A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1]
-                += A[(1 + dim)*elemIdx + i + 1]
-                  [(1 + dim)*elemIdx + j + 1];
             }
           }
         }
       }
+
+      //stiffness matrix A is symmetric
+      for (int i = 0; i < dim; i++)
+      {
+        A[elemIdxSLE][elemIdxSLE + i + 1] = A[elemIdxSLE + i + 1][elemIdxSLE];
+
+        for(int j = 0; j < i; j++)
+        {
+          A[elemIdxSLE + j + 1][elemIdxSLE + i + 1] =
+            A[elemIdxSLE + i + 1][elemIdxSLE + j + 1];
+        }
+      }
     }
 
     std::cout << A << std::endl;
@@ -315,6 +315,12 @@ public:
     //NOTE: what would be an appropiate solver here?
     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;
+
     //storage for pressure data, for each element we store the
     //pressure at the corners of the element, the VTKFunction will
     //create the output using a linear interpolation for each element
@@ -323,7 +329,7 @@ public:
 
     for (const auto& elem : elements(gridView_))
     {
-      const int elemIdx = mapper_.index(elem);
+      const int elemIdxSLE = (dim + 1)*mapper_.index(elem);
       const auto& geo = elem.geometry();
 
       for (int k = 0; k < geo.corners(); k++)
@@ -331,15 +337,14 @@ public:
         //contribution of the basis function
         // phi_elem,0 (x) = indicator(elem);
         //at the kth corner of elem
-        pressure[(dim + 1)*elemIdx + k] = d[(1 + dim)*elemIdx];
+        pressure[elemIdxSLE + k] = d[elemIdxSLE];
 
         for (int i = 0; i < dim; i++)
         {
           //contribution of the basis function
           // phi_elem,i (x) = x[i]*indicator(elem);
           //at the kth corner of elem
-          pressure[(dim + 1)*elemIdx + k] +=
-            d[(1 + dim)*elemIdx + i + 1] * geo.corner(k)[i];
+          pressure[elemIdxSLE + k] += d[elemIdxSLE + i + 1] * geo.corner(k)[i];
         }
       }
     }
diff --git a/src/dg.cc b/src/dg.cc
index 58bfe30efc8fd91c216d9bd53118880f795912a9..945503e3a4571f9366b2c7cf15a35ee0126cb576 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -10,17 +10,12 @@
 #include <dune/common/exceptions.hh>
 #include <dune/common/parametertree.hh>
 #include <dune/common/parametertreeparser.hh>
-#include <dune/common/dynmatrix.hh>
-#include <dune/common/dynvector.hh>
 
 #include <dune/grid/common/mcmgmapper.hh>
 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
 #include <dune/grid/io/file/dgfparser/dgfug.hh>
-#include <dune/grid/io/file/vtk.hh>
 #include <dune/grid/yaspgrid.hh>
 
-#include <dune/geometry/quadraturerules.hh>
-
 #include <dune/mmesh/mmesh.hh>
 
 #include <dune/mmdg/dg.hh>
@@ -38,12 +33,12 @@ int main(int argc, char** argv)
 
     static constexpr int dim = GRIDDIM;
 
+    //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 GridView = typename Grid::LeafGridView;
     using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
-    using Coordinate = Dune::FieldVector<double, dim>;
 
     //Message Passing Interface (MPI) for parallel computation
     Dune::MPIHelper::instance(argc, argv);