[bugfix] fix index error in an interface entry of the load vector

dune/mmdg/mmdg.hh

@@ -12,7 +12,7 @@ public:
   static constexpr auto iSimplex = Dune::GeometryTypes::simplex(dim-1);
 
   using Base = DG<GridView, Mapper, Problem>;
-  using Scalar = typename Base::Scalar; //NOTE using Base::Scalar
+  using Scalar = typename Base::Scalar;
   using Matrix = typename Base::Matrix;
   using Vector = typename Base::Vector;
   using Coordinate = Dune::FieldVector<Scalar, dim>;
@@ -203,22 +203,25 @@ private:
 
         const Scalar quadratureFactor =
           ip.weight() * iGeo.integrationElement(qpLocal);
-        const Scalar betaEvaluation = beta(qpGlobal);
+        const Scalar interfaceUpdate1 = quadratureFactor * beta(qpGlobal);
 
         for (int i = 0; i < dim - 1; i++)
         {
-          for (int j = 0; j <= i; j++)
-          {
-            (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
-              quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal)
-              * (iFrame[j] * qpGlobal);
+          const Scalar qpI = iFrame[i] * qpGlobal;
+          const Scalar interfaceUpdate2 = interfaceUpdate1 * qpI;
 
-            if (i != j)
+          (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
+            interfaceUpdate2 * qpI;
+
+          for (int j = 0; j < i; j++)
           {
+            const Scalar interfaceUpdate3 =
+              interfaceUpdate2 * (iFrame[j] * qpGlobal);
+
+            (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
+              interfaceUpdate3;
             (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
-                quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal)
-                * (iFrame[j] * qpGlobal);
-            }
+              interfaceUpdate3;
           }
         }
       }
@@ -571,9 +574,10 @@ private:
           // int_intersct d * g * (mu * phi_iElem,i +
           //      (Kparallel * grad(phi_i,iElem)) * normal) dr
           //for i = 0,...,dim-1
-          const Scalar dgGamma = Base::problem_.aperture(center)
-            * Base::problem_.boundary(center);
-          Base::b[iElemIdxSLE] += mu * dgGamma;
+          const Scalar dEvaluation = Base::problem_.aperture(center);
+          const Scalar dgGamma = dEvaluation * Base::problem_.boundary(center);
+
+          Base::b[iElemIdxSLE] += mu * dgGamma * dEvaluation;
 
           for (int i = 0; i < dim - 1; i++)
           {
@@ -584,7 +588,8 @@ private:
             const Scalar interfaceUpdate3 =
               interfaceUpdate1 - interfaceUpdate2;
 
-            Base::b[iElemIdxSLE] += dgGamma * (mu * centerI - interfaceUpdate2);
+            Base::b[iElemIdxSLE + i + 1] +=
+              dgGamma * dEvaluation * (mu * centerI - interfaceUpdate2);
 
             (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
               interfaceUpdate3;