diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 102d0e254d6286ee7eaa6ed4ee3328a3a9d38823..0983c462f89e97d8d17228071061f29e10611bd6 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -40,7 +40,7 @@ public:
       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 << "NON-symmetric!: " << i << ", " << j << std::endl;
 
     std::cout << *Base::A << std::endl;
 
@@ -232,6 +232,10 @@ private:
         }
       }
 
+      std::cout << iFrame[0] << "\t" << (*Base::A)[iElemIdxSLE][iElemIdxSLE] << "\t" <<
+        (*Base::A)[iElemIdxSLE+1][iElemIdxSLE] << "\t" << (*Base::A)[iElemIdxSLE][iElemIdxSLE+1] << "\t" <<
+        (*Base::A)[iElemIdxSLE+1][iElemIdxSLE+1] << "\t" << std::endl;
+
       //get neighboring bulk grid elements of iELem
       const auto elemIn = intersct.inside();
       const auto elemOut = intersct.outside();
@@ -550,6 +554,16 @@ private:
           const int neighborIdxSLE =
             (dim + 1) * Base::gridView_.size(0) + dim * neighborIdx;
 
+          //get basis of the interface coordinate system of the neighboring
+          //interface element, for uncurved interfaces the sign of normals
+          //and hence the frames of iElem and its neighbor can differ
+          const InterfaceBasis neighborFrame =
+            interfaceFrame(Base::gridView_.grid().asIntersection(
+              intersct.outside()).centerUnitOuterNormal());
+
+          //storage for K_par * neighborFrame[i]
+          Coordinate KparDotNeighborTau_i;
+
           //evaluate the interface terms //TODO: terms with grad(d)
           // int_intersct mu^Gamma * jump(phi_iElem,i) * jump(phi_neighbor,j) dr
           //and
@@ -565,49 +579,50 @@ private:
           for (int i = 0; i < dim - 1; i++)
           {
             Base::problem_.Kparallel(center).mv(iFrame[i], KparDotTau_i);
+            Base::problem_.Kparallel(center).mv(
+              neighborFrame[i], KparDotNeighborTau_i);
             const Scalar centerI = center * iFrame[i];
-            const Scalar interfaceUpdate1 = mu * centerI;
-            const Scalar interfaceUpdate2 = KparDotTau_i * normal;
-            const Scalar interfaceUpdate3 =
-              dEvaluation2 * (-interfaceUpdate1 + 0.5 * interfaceUpdate2);
-            const Scalar interfaceUpdate4 =
-              dEvaluation2 * (-interfaceUpdate1 - 0.5 * interfaceUpdate2);
+            const Scalar neighborCenterI = center * neighborFrame[i];
+            const Scalar interfaceUpdate1a = mu * centerI;
+            const Scalar interfaceUpdate1b = mu * neighborCenterI;
+            const Scalar interfaceUpdate2a = KparDotTau_i * normal;
+            const Scalar interfaceUpdate2b = KparDotNeighborTau_i * normal;
+            const Scalar interfaceUpdate3a =
+              dEvaluation2 * (-interfaceUpdate1a + 0.5 * interfaceUpdate2a);
+            const Scalar interfaceUpdate3b =
+              dEvaluation2 * (-interfaceUpdate1b - 0.5 * interfaceUpdate2b);
 
             (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE] +=
-              interfaceUpdate3;
+              interfaceUpdate3a;
             (*Base::A)[neighborIdxSLE][iElemIdxSLE + i + 1] +=
-              interfaceUpdate3;
+              interfaceUpdate3a;
             (*Base::A)[neighborIdxSLE + i + 1][iElemIdxSLE] +=
-              interfaceUpdate4;
+              interfaceUpdate3b;
             (*Base::A)[iElemIdxSLE][neighborIdxSLE + i + 1] +=
-              interfaceUpdate4;
+              interfaceUpdate3b;
 
-            (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE + i + 1] +=
-              -centerI * dEvaluation2 * interfaceUpdate1;
+            for (int j = 0; j < dim - 1; j++)
+            {
+              Coordinate KparDotNeighborTau_j;
+
+              Base::problem_.Kparallel(center).mv(
+                neighborFrame[j], KparDotNeighborTau_j);
 
-            for (int j = 0; j < i; j++)
-            { //NOTE: only relevant for n=3, requires quadrature rule
-              Coordinate KparDotTau_j;
-              Base::problem_.Kparallel(center).mv(iFrame[j], KparDotTau_j);
               const Scalar interfaceUpdate5 =
-                (center * iFrame[j]) * interfaceUpdate1;
-              const Scalar interfaceUpdate6 = interfaceUpdate2 *
-                (center * iFrame[j]) - (KparDotTau_j * normal) * centerI;
-              const Scalar interfaceUpdate7 =
-                dEvaluation2 * (-interfaceUpdate5 - 0.5 * interfaceUpdate6);
-              const Scalar interfaceUpdate8 =
-                dEvaluation2 * (-interfaceUpdate5 + 0.5 * interfaceUpdate6);
+                -dEvaluation2 * ((center * neighborFrame[j]) * interfaceUpdate1a
+                + 0.5 * (-interfaceUpdate2a * (center * neighborFrame[j])
+                + (KparDotNeighborTau_j * normal) * centerI));
 
               (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
-                interfaceUpdate7;
+                interfaceUpdate5;
               (*Base::A)[neighborIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
-                interfaceUpdate7;
-              (*Base::A)[iElemIdxSLE + j + 1][neighborIdxSLE + i + 1] +=
-                interfaceUpdate8;
-              (*Base::A)[neighborIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
-                interfaceUpdate8;
+                interfaceUpdate5;
             }
           }
+
+          std::cout << "neighborTerms: " << (*Base::A)[neighborIdxSLE][iElemIdxSLE] << "\t" <<
+          (*Base::A)[neighborIdxSLE][iElemIdxSLE+1] << "\t" << (*Base::A)[neighborIdxSLE+1][iElemIdxSLE] << "\t" <<
+          (*Base::A)[neighborIdxSLE+1][iElemIdxSLE+1] << "\t\n\n";
         }
         else //boundary interface edge
         {
@@ -663,6 +678,9 @@ private:
           }
         }
       }
+      std::cout << iFrame[0] << "\t" << (*Base::A)[iElemIdxSLE][iElemIdxSLE] << "\t" <<
+        (*Base::A)[iElemIdxSLE+1][iElemIdxSLE] << "\t" << (*Base::A)[iElemIdxSLE][iElemIdxSLE+1] << "\t" <<
+        (*Base::A)[iElemIdxSLE+1][iElemIdxSLE+1] << "\t" << std::endl;
     }
   }