[bugfix] add missing coupling entries, [bugfix] fix index errors for coupling entries

dune/mmdg/mmdg.hh

@@ -384,11 +384,44 @@ private:
         }
 
         //i = dim - 1
-        const Scalar couplingUpdate4 = qpGlobal[dim] * couplingUpdate1;
-        (*Base::A)[elemInIdxSLE + dim + 1][iElemIdxSLE] -= couplingUpdate4;
-        (*Base::A)[iElemIdxSLE][elemInIdxSLE + dim + 1] -= couplingUpdate4;
-        (*Base::A)[elemOutIdxSLE + dim + 1][iElemIdxSLE] -= couplingUpdate4;
-        (*Base::A)[iElemIdxSLE][elemOutIdxSLE + dim + 1] -= couplingUpdate4;
+        const Scalar couplingUpdate4 = qpGlobal[dim - 1] * couplingUpdate1;
+
+        (*Base::A)[elemInIdxSLE + dim][iElemIdxSLE] -= couplingUpdate4;
+        (*Base::A)[iElemIdxSLE][elemInIdxSLE + dim] -= couplingUpdate4;
+        (*Base::A)[elemOutIdxSLE + dim][iElemIdxSLE] -= couplingUpdate4;
+        (*Base::A)[iElemIdxSLE][elemOutIdxSLE + dim] -= couplingUpdate4;
+      }
+
+      //evalute the coupling terms
+      // -int_iElem beta * phi_iElem,i * avg(phi_elem,j) ds
+      //for i = 1,...,dim-1 and j = 1,...,dim
+      //and for elem in {elemIn, elemOut} using a quadrature rule
+      for (const auto& ip : rule_KperpPlus1) //TODO: symmetrize more efficiently
+      //TODO: apply quadrature rule only once
+      { //TODO TODO TODO
+        const auto& qpLocal = ip.position();
+        const auto& qpGlobal = iGeo.global(qpLocal);
+
+        const Scalar quadratureFactor =
+          ip.weight() * iGeo.integrationElement(qpLocal);
+        const Scalar betaEvaluation2 = 0.5 * beta(qpGlobal);
+
+        for (int i = 0; i < dim; i++)
+        {
+          for (int j = 0; j < dim - 1; j++)
+          {
+            const Scalar couplingUpdate = quadratureFactor * betaEvaluation2
+              * (iFrame[j] * qpGlobal) * qpGlobal[i];
+            (*Base::A)[elemInIdxSLE + i + 1][iElemIdxSLE + j + 1] -=
+              couplingUpdate;
+            (*Base::A)[iElemIdxSLE + j + 1][elemInIdxSLE + i + 1] -=
+              couplingUpdate;
+            (*Base::A)[elemOutIdxSLE + i + 1][iElemIdxSLE + j + 1] -=
+              couplingUpdate;
+            (*Base::A)[iElemIdxSLE + j + 1][elemOutIdxSLE + i + 1] -=
+              couplingUpdate;
+          }
+        }
       }
 
       // === interface entries ===