diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index e76112b2b7a3b7357c1b87c643bc72a1a532fb57..f208821357a1d5c628ab74e1f79542738070aa9e 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -40,8 +40,6 @@ public:
     //NOTE: what would be an appropiate solver here?
     A.solve(d, b);
 
-    std::cout << b << std::endl;
-
     for (int i = 0; i < dof; i++)
       for (int j = 0; j < i; j++)
         /*assert*/if(std::abs(A[i][j] - A[j][i]) >
@@ -127,377 +125,288 @@ private:
   //assemble stiffness matrix A and load vector b
   void assembleSLE (const Scalar K, const Scalar mu)
   {
-  //NOTE:
-//const int order = 3; //order of the quadrature rule
+    //NOTE:
+    //const int order = 3; //order of the quadrature rule
 
-  for (const auto& elem : elements(gridView_))
-  {
-    const auto& geo = elem.geometry();
-    std::cout << mapper_.index(elem) << "\t";
-    for (int k = 0; k < geo.corners(); k++)
-    {
-      std::cout << geo.corner(k) << "\t";
-    }
-    std::cout << "\n\n";
-  }
-
-  //we use the basis
-  // phi_elem,0 (x) = indicator(elem);
-  // phi_elem,i (x) = x[i]*indicator(elem);
-  //for elem in elements(gridView) and i = 1,...,dim
-  for (const auto& elem : elements(gridView_))
-  {
-    const int elemIdx = mapper_.index(elem);
-    const auto& geo = elem.geometry();
-    const double elemVol = geo.volume();
-    const auto& center = geo.center();
-
-    std::cout << "===============================\nElement " << elemIdx;
-    for (int k = 0; k < geo.corners(); k++)
+    //we use the basis
+    // phi_elem,0 (x) = indicator(elem);
+    // phi_elem,i (x) = x[i]*indicator(elem);
+    //for elem in elements(gridView) and i = 1,...,dim
+    for (const auto& elem : elements(gridView_))
     {
-      std::cout << "\t" << geo.corner(k);
-    }
-    std::cout << "\n\n";
-
-    //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;
+      const int elemIdx = mapper_.index(elem);
+      const auto& geo = elem.geometry();
+      const double elemVol = geo.volume();
+      const auto& center = geo.center();
 
-/*    //TODO: can be done outside of the loop?
-    const Dune::QuadratureRule<double,dim>& rule =
-      Dune::QuadratureRules<double,dim>::rule(geo.type(),order);
+      //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;
 
-    //NOTE: how are quadrature rules in Dune applied correctly?
-    for (const auto& ip : rule)
-    {
-      const auto qp = ip.position();
-      //NOTE: is the volume of the element taken into account
-      //automatically?
-      const auto weight = ip.weight();
+  /*    //TODO: can be done outside of the loop?
+      const Dune::QuadratureRule<double,dim>& rule =
+        Dune::QuadratureRules<double,dim>::rule(geo.type(),order);
 
-      //quadrature for int_elem q*phi_elem,0 dV
-      b[elemIdxSLE] += weight * q(qp);
-
-      //quadrature for int_elem q*phi_elem,i dV
-      for (int i = 0; i < dim; i++)
+      //NOTE: how are quadrature rules in Dune applied correctly?
+      for (const auto& ip : rule)
       {
-        b[elemIdxSLE + i + 1] += weight * qp[i] * q(qp);
-      }
-    }
-*/
-
-    //NOTE: makeshift solution for source term q = -1
-    b[elemIdxSLE] += -elemVol;
-    for (int i = 0; i < dim; i++)
-    {
-      b[elemIdxSLE + i + 1] += -elemVol * center[i];
-    }
+        const auto qp = ip.position();
+        //NOTE: is the volume of the element taken into account
+        //automatically?
+        const auto weight = ip.weight();
 
-    for (int i = 0; i < dim; i++)
-    {
-      //exact evaluation of
-      // int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV
-      A[elemIdxSLE + i + 1][elemIdxSLE + i + 1] += K * elemVol;
-    }
+        //quadrature for int_elem q*phi_elem,0 dV
+        b[elemIdxSLE] += weight * q(qp);
 
-    //iterate over all intersection with the boundary of elem
-    for (const auto& intersct : intersections(gridView_, elem))
-    {
-      const auto& normal = intersct.centerUnitOuterNormal();
-      const auto& intersctGeo = intersct.geometry();
-      const double intersctVol = intersctGeo.volume();
-      const auto& intersctCenter = intersctGeo.center();
+        //quadrature for int_elem q*phi_elem,i dV
+        for (int i = 0; i < dim; i++)
+        {
+          b[elemIdxSLE + i + 1] += weight * qp[i] * q(qp);
+        }
+      }
+  */
 
-      std::cout << "+++++++++++++++++++++++++++\nIntersection ";
-      for (int k = 0; k < intersctGeo.corners(); k++)
+      //NOTE: makeshift solution for source term q = -1
+      b[elemIdxSLE] += -elemVol;
+      for (int i = 0; i < dim; i++)
       {
-        std::cout << "\t" << intersctGeo.corner(k);
+        b[elemIdxSLE + i + 1] += -elemVol * center[i];
       }
-      std::cout << "\tnormal: " << normal << "\n\n";
-
-      //TODO: quadrature rule cannot be used for dim = 1!
-//        const Dune::QuadratureRule<double,dim-1>& secondOrderRule =
-//          Dune::QuadratureRules<double,dim-1>::rule(
-//          intersct.type(), 2);
-
-      //storage for multiply used integral values,
-      //note that only the lower diagonal and diagonal entries of
-      //quadraticIntregrals are used
-      Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);
-      Dune::FieldMatrix<Scalar, dim, dim> quadraticIntregrals(0.0);
-      //NOTE: is there a better type for symmetric matrix?
 
       for (int i = 0; i < dim; i++)
       {
-        //use midpoint rule for exact evaluation of
-        // int_intersct x_i ds
-        linearIntegrals[i] = intersctVol * intersctCenter[i];
-
-        for (int j = 0; j <= i; j++)
-        {
-          const auto& leftCorner = intersctGeo.corner(0);
-          const auto& rightCorner = intersctGeo.corner(1);
-          std::cout << "left corner: " << leftCorner << "\tright corner: " << rightCorner << "\n";
-
-          quadraticIntregrals[i][j] = intersctVol / 3 *
-            ( leftCorner[i] * leftCorner[j] +
-              rightCorner[i] * rightCorner[j]
-            + 0.5 * (leftCorner[i] * rightCorner[j] +
-              leftCorner[j] * rightCorner[i]) );
-
-          std::cout << i << ", " << j << "\t" << leftCorner[i] * leftCorner[j] << "\t" <<
-            rightCorner[i] * rightCorner[j] << "\t" <<
-            0.5 * (leftCorner[i] * rightCorner[j] +
-            leftCorner[j] * rightCorner[i]) << "\t" << intersctVol / 3 *
-              ( leftCorner[i] * leftCorner[j] +
-                rightCorner[i] * rightCorner[j]
-              + 0.5 * (leftCorner[i] * rightCorner[j] +
-                leftCorner[j] * rightCorner[i]) ) << "\t" << quadraticIntregrals[i][j] << "\n\n";
-/*
-          //use second order quadrature rule for exact evaluation of
-          // int_intersct x_i*x_j ds
-          //NOTE: use Simpson's rule instead manually?
-          for (const auto& ip : secondOrderRule)
-          {
-            const auto& qp = ip.position();
-            quadraticIntregrals[i][j] += //NOTE: volume necessary?
-              ip.weight() * qp[i] * qp[j] * intersctVol;
-          }
-*/
-          //NOTE: probably unnecessary
-          quadraticIntregrals[j][i] = quadraticIntregrals[i][j];
-        }
+        //exact evaluation of
+        // int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV
+        A[elemIdxSLE + i + 1][elemIdxSLE + i + 1] += K * elemVol;
       }
-      std::cout << "linearIntegrals:\n" << linearIntegrals << "\n\n";
-      std::cout << "quadraticIntregrals:\n" << quadraticIntregrals << "\n\n";
-
-      //exact evaluation of
-      // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
-      A[elemIdxSLE][elemIdxSLE] += mu * intersctVol;
-      std::cout << elemIdxSLE << ", " << elemIdxSLE <<"\t"  << intersctCenter<< ":\n"
-        << mu * intersctVol << "\t" << A[elemIdx][elemIdx] << "\n\n";
 
-      if (intersct.neighbor()) //intersct has neighboring element
+      //iterate over all intersection with the boundary of elem
+      for (const auto& intersct : intersections(gridView_, elem))
       {
-        std::cout << "------------------------------\nNeighbor\t";
-        for (int k = 0; k < intersct.outside().geometry().corners(); k++)
-        {
-          std::cout << "\t" << intersct.outside().geometry().corner(k);
-        }
-        std::cout << "\n\n";
-
-
-        //index of the neighboring element
-        const int neighborIdx = mapper_.index(intersct.outside());
-        const int neighborIdxSLE = (dim + 1)*neighborIdx;
+        const auto& normal = intersct.centerUnitOuterNormal();
+        const auto& intersctGeo = intersct.geometry();
+        const double intersctVol = intersctGeo.volume();
+        const auto& intersctCenter = intersctGeo.center();
+
+        //TODO: quadrature rule cannot be used for dim = 1!
+  //        const Dune::QuadratureRule<double,dim-1>& secondOrderRule =
+  //          Dune::QuadratureRules<double,dim-1>::rule(
+  //          intersct.type(), 2);
+
+        //storage for multiply used integral values,
+        //note that only the lower diagonal and diagonal entries of
+        //quadraticIntregrals are used
+        Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);
+        Dune::FieldMatrix<Scalar, dim, dim> quadraticIntregrals(0.0);
+        //NOTE: is there a better type for a symmetric matrix?
 
         for (int i = 0; i < dim; i++)
-        { //we use the relations
-          // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
-          // = mu * int_intersct x_i ds
-          //and
-          // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds
-          // = 0.5 * K * normal[i] * vol(intersct)
-          A[elemIdxSLE + i + 1][elemIdxSLE] +=
-            mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
-          std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE << ":\n"
-            << mu * linearIntegrals[i] << "\t" <<
-            - 0.5 * K * normal[i] * intersctVol << "\t"
-            << A[elemIdxSLE + i + 1][elemIdxSLE] << "\n\n";
+        {
+          //use midpoint rule for exact evaluation of
+          // int_intersct x_i ds
+          linearIntegrals[i] = intersctVol * intersctCenter[i];
 
           for (int j = 0; j <= i; j++)
           {
-            //we use the relations
-            // int_intersct mu*jump(phi_elem,i)*jump(phi_elem,j) ds
-            // = mu * int_intersct x_i*x_j ds
-            //and
-            // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,j) ds
-            // = 0.5 * K * normal[i] * int_intersct x_j ds
-            A[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
-              += mu * quadraticIntregrals[i][j]
-                - 0.5 * K * (normal[i] * linearIntegrals[j]
-                  + normal[j] * linearIntegrals[i]);
-            std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE + j + 1 << ":\n"
-              << mu * quadraticIntregrals[i][j] << "\t" <<
-              - 0.5 * K * (normal[i] * linearIntegrals[j]
-                + normal[j] * linearIntegrals[i]) << "\t" <<
-              A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] << "\n\n";
-          }
-        }
+            const auto& leftCorner = intersctGeo.corner(0);
+            const auto& rightCorner = intersctGeo.corner(1);
 
-        if (neighborIdx > elemIdx)
-        { //make sure that each facet is considered only once
-          continue;
+            quadraticIntregrals[i][j] = intersctVol / 3 *
+              ( leftCorner[i] * leftCorner[j] + rightCorner[i] * rightCorner[j]
+              + 0.5 * (leftCorner[i] * rightCorner[j] +
+                leftCorner[j] * rightCorner[i]) );
+  /*
+            //use second order quadrature rule for exact evaluation of
+            // int_intersct x_i*x_j ds
+            //NOTE: use Simpson's rule instead manually?
+            for (const auto& ip : secondOrderRule)
+            {
+              const auto& qp = ip.position();
+              quadraticIntregrals[i][j] += //NOTE: volume necessary?
+                ip.weight() * qp[i] * qp[j] * intersctVol;
+            }
+  */
+            //NOTE: probably unnecessary
+            quadraticIntregrals[j][i] = quadraticIntregrals[i][j];
+          }
         }
 
         //exact evaluation of
-        // int_intersct mu*jump(phi_elem,0)*jump(phi_neighbor,0) ds
-        A[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;
-
-        //stiffness matrix A is symmetric
-        A[neighborIdxSLE][elemIdxSLE] += A[elemIdxSLE][neighborIdxSLE];
+        // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
+        A[elemIdxSLE][elemIdxSLE] += mu * intersctVol;
 
-        for (int i = 0; i < dim; i++)
+        if (intersct.neighbor()) //intersct has neighboring element
         {
-          //we use the relations
-          // int_intersct mu*jump(phi_elem,i)
-          //  *jump(phi_neighbor,0) ds
-          // = -mu*int_intersct x_i ds
-          //and
-          // int_intersct avg(K*grad(phi_neighbor,i))
-          //  *jump(phi_elem,0) ds
-          // = 0.5 * K * normal[i] * vol(intersct)
-          A[elemIdxSLE + i + 1][neighborIdxSLE] +=
-            -mu * linearIntegrals[i] + 0.5 * K * normal[i] * intersctVol;
-          std::cout << elemIdxSLE + i + 1 << ", " << neighborIdxSLE << ":\n"
-            <<   -mu * linearIntegrals[i] << "\t" <<
-            0.5 * K * normal[i] * intersctVol << "\t" <<
-            A[elemIdxSLE + i + 1][neighborIdxSLE] << "\n\n";
-
-          //we use the relations
-          // int_intersct mu*jump(phi_elem,0)
-          //  *jump(phi_neighbor,i) ds
-          // = -mu*int_intersct x_i ds
-          //and
-          // int_intersct avg(K*grad(phi_neighbor,i))
-          //  *jump(phi_elem,0) ds
-          // = 0.5 * K * normal[i] * vol(intersct)
-          A[elemIdxSLE][neighborIdxSLE + i + 1] +=
-            -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
-          std::cout << elemIdxSLE << ", " << neighborIdxSLE + i + 1<< ":\n"
-            << -mu * linearIntegrals[i] << "\t" <<
-            -0.5 * K * normal[i] * intersctVol << "\t" <<
-            A[elemIdxSLE][neighborIdxSLE + i + 1] << "\n\n";
+          //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
+            // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
+            // = mu * int_intersct x_i ds
+            //and
+            // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds
+            // = 0.5 * K * normal[i] * vol(intersct)
+            A[elemIdxSLE + i + 1][elemIdxSLE] +=
+              mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
+
+            for (int j = 0; j <= i; j++)
+            {
+              //we use the relations
+              // int_intersct mu*jump(phi_elem,i)*jump(phi_elem,j) ds
+              // = mu * int_intersct x_i*x_j ds
+              //and
+              // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,j) ds
+              // = 0.5 * K * normal[i] * int_intersct x_j ds
+              A[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
+                += mu * quadraticIntregrals[i][j]
+                  - 0.5 * K * (normal[i] * linearIntegrals[j]
+                    + normal[j] * linearIntegrals[i]);
+            }
+          }
+
+          if (neighborIdx > elemIdx)
+          { //make sure that each facet is considered only once
+            continue;
+          }
+
+          //exact evaluation of
+          // int_intersct mu*jump(phi_elem,0)*jump(phi_neighbor,0) ds
+          A[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;
 
           //stiffness matrix A is symmetric
-          A[neighborIdxSLE][elemIdxSLE + i + 1] +=
-            A[elemIdxSLE + i + 1][neighborIdxSLE];
-          A[neighborIdxSLE + i + 1][elemIdxSLE] +=
-            A[elemIdxSLE][neighborIdxSLE + i + 1];
+          A[neighborIdxSLE][elemIdxSLE] += A[elemIdxSLE][neighborIdxSLE];
 
-          for (int j = 0; j <= i; j++)
+          for (int i = 0; i < dim; i++)
           {
             //we use the relations
             // int_intersct mu*jump(phi_elem,i)
-            //  *jump(phi_neighbor,j) ds
-            // = -mu*int_intersct x_i*x_j ds
+            //  *jump(phi_neighbor,0) ds
+            // = -mu*int_intersct x_i ds
             //and
-            // int_intersct avg(K*grad(phi_neighbor,j))
-            //  *jump(phi_elem,i) ds
-            // = 0.5 * K * normal[j] * int_intersct x_i ds
-            //as well as
-            // int_intersct avg(K*grad(phi_elem,i))
-            //  *jump(phi_neighbor,j) ds
-            // = -0.5 * K * normal[i] * int_intersct x_j ds
-            A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
-              -mu * quadraticIntregrals[i][j]
-              - 0.5 * K * (normal[j] * linearIntegrals[i]
-                - normal[i] * linearIntegrals[j]);
-            std::cout << elemIdxSLE + i + 1 << ", " << neighborIdxSLE + j + 1<< ":\n"
-              << -mu * quadraticIntregrals[i][j] << "\t" <<
-              - 0.5 * K * (normal[j] * linearIntegrals[i]
-                - normal[i] * linearIntegrals[j]) << "\t" <<
-              A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] << "\n\n";
+            // int_intersct avg(K*grad(phi_neighbor,i))
+            //  *jump(phi_elem,0) ds
+            // = 0.5 * K * normal[i] * vol(intersct)
+            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)
+            //  *jump(phi_neighbor,i) ds
+            // = -mu*int_intersct x_i ds
+            //and
+            // int_intersct avg(K*grad(phi_neighbor,i))
+            //  *jump(phi_elem,0) ds
+            // = 0.5 * K * normal[i] * vol(intersct)
+            A[elemIdxSLE][neighborIdxSLE + i + 1] +=
+              -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
 
             //stiffness matrix A is symmetric
-            A[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
-              A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1];
+            A[neighborIdxSLE][elemIdxSLE + i + 1] +=
+              A[elemIdxSLE + i + 1][neighborIdxSLE];
+            A[neighborIdxSLE + i + 1][elemIdxSLE] +=
+              A[elemIdxSLE][neighborIdxSLE + i + 1];
 
-            if (i != j)
+            for (int j = 0; j <= i; j++)
             {
-              // int_intersct mu*jump(phi_elem,j)
-              //  *jump(phi_neighbor,i) ds
+              //we use the relations
+              // int_intersct mu*jump(phi_elem,i)
+              //  *jump(phi_neighbor,j) 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
+              // int_intersct avg(K*grad(phi_neighbor,j))
+              //  *jump(phi_elem,i) ds
+              // = 0.5 * K * normal[j] * int_intersct x_i 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] +=
+              // int_intersct avg(K*grad(phi_elem,i))
+              //  *jump(phi_neighbor,j) ds
+              // = -0.5 * K * normal[i] * int_intersct x_j ds
+              A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
                 -mu * quadraticIntregrals[i][j]
-                - 0.5 * K * (normal[i] * linearIntegrals[j]
-                  - normal[j] * linearIntegrals[i]);
-              std::cout << elemIdxSLE + j + 1 << ", " << neighborIdxSLE + i + 1<< ":\n"
-                << -mu * quadraticIntregrals[i][j] << "\t" <<
-                - 0.5 * K * (normal[i] * linearIntegrals[j]
-                  - normal[j] * linearIntegrals[i]) << "\t" <<
-                A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] << "\n\n";
+                - 0.5 * K * (normal[j] * linearIntegrals[i]
+                  - normal[i] * linearIntegrals[j]);
 
               //stiffness matrix A is symmetric
-              A[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] +=
-                A[elemIdxSLE + j + 1][neighborIdxSLE + 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];
+              }
             }
           }
         }
-      }
-      else //boundary facet
-      {
-        std::cout << "------------------------------\nBoundary\n\n";
-        for (int i = 0; i < dim; i++)
-        { //we use the relations
-          // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
-          // = mu * int_intersct x_i ds
-          //and for boundary facets
-          // int_intersct avg(K*grad(phi_elem,i))
-          //  *jump(phi_elem,0) ds
-          // = K * normal[i] * vol(intersct)
-          A[elemIdxSLE + i + 1][elemIdxSLE] +=
-            mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
-          std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE << "\t"
-          << intersctCenter << ":\n" << mu * linearIntegrals[i] << "\t" <<
-            - 0.5 * K * normal[i] * intersctVol << "\t"
-          << A[elemIdxSLE + i + 1][elemIdxSLE] << "\n\n";
-
-
-          for (int j = 0; j <= i; j++)
-          {
-            //we use the relations
-            // int_intersct mu*jump(phi_elem,i)*jump(phi_elem,j) ds
-            // = mu * int_intersct x_i*x_j ds
+        else //boundary facet
+        {
+          for (int i = 0; i < dim; i++)
+          { //we use the relations
+            // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
+            // = mu * int_intersct x_i ds
             //and for boundary facets
             // int_intersct avg(K*grad(phi_elem,i))
-            //  *jump(phi_elem,j) ds
-            // = 0.5 * K * normal[i] * int_intersct x_j ds
-            A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
-              mu * quadraticIntregrals[i][j]
-              - 0.5 * K * (normal[i] * linearIntegrals[j]
-                + normal[j] * linearIntegrals[i]);
-            std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE + j+ 1 << "\t"
-              << intersctCenter << ":\n" << mu * quadraticIntregrals[i][j] << "\t" <<
-              - 0.5 * K * (normal[i] * linearIntegrals[j]
-                + normal[j] * linearIntegrals[i]) << "\t"
-              << A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] << "\n\n";
+            //  *jump(phi_elem,0) ds
+            // = K * normal[i] * vol(intersct)
+            A[elemIdxSLE + i + 1][elemIdxSLE] +=
+              mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
+
+            for (int j = 0; j <= i; j++)
+            {
+              //we use the relations
+              // int_intersct mu*jump(phi_elem,i)*jump(phi_elem,j) ds
+              // = mu * int_intersct x_i*x_j ds
+              //and for boundary facets
+              // int_intersct avg(K*grad(phi_elem,i))
+              //  *jump(phi_elem,j) ds
+              // = 0.5 * K * normal[i] * int_intersct x_j ds
+              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
-    for (int i = 0; i < dim; i++)
-    {
-      A[elemIdxSLE][elemIdxSLE + i + 1] = A[elemIdxSLE + i + 1][elemIdxSLE];
-
-      for(int j = 0; j < i; j++)
+      //stiffness matrix A is symmetric
+      for (int i = 0; i < dim; i++)
       {
-        A[elemIdxSLE + j + 1][elemIdxSLE + i + 1] =
-          A[elemIdxSLE + i + 1][elemIdxSLE + j + 1];
+        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;
-  }
-
 
   Matrix A; //stiffness matrix
   Vector b; //load vector
   Vector d; //solution vector
-
 };
 
 #endif