[bugfix] fix problem mmdg6

dune/mmdg/mmdg.hh

@@ -236,16 +236,22 @@ private:
 
               //quadrature for
               // int_iElem beta * phi_iElem,0 * phi_iElem,i ds
-              Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) += interfaceUpdate;
-              Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) += interfaceUpdate;
+              Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
+                interfaceUpdate;
+              Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
+                interfaceUpdate;
 
               //quadrature for
               // -int_iElem beta * phi_iElem,i * avg(phi_elem,0) ds
               //for elem in {elemIn, elemOut} using a quadrature rule
-              Base::A->entry(elemInIdxSLE, iElemIdxSLE + i + 1) -= couplingUpdate3;
-              Base::A->entry(iElemIdxSLE + i + 1, elemInIdxSLE) -= couplingUpdate3;
-              Base::A->entry(elemOutIdxSLE, iElemIdxSLE + i + 1) -= couplingUpdate3;
-              Base::A->entry(iElemIdxSLE + i + 1, elemOutIdxSLE) -= couplingUpdate3;
+              Base::A->entry(elemInIdxSLE, iElemIdxSLE + i + 1) -=
+                couplingUpdate3;
+              Base::A->entry(iElemIdxSLE + i + 1, elemInIdxSLE) -=
+                couplingUpdate3;
+              Base::A->entry(elemOutIdxSLE, iElemIdxSLE + i + 1) -=
+                couplingUpdate3;
+              Base::A->entry(iElemIdxSLE + i + 1, elemOutIdxSLE) -=
+                couplingUpdate3;
             }
 
             //quadrature for
@@ -265,10 +271,14 @@ private:
             //quadrature for
             // -int_iElem beta * phi_iElem,0 * avg(phi_elem,i) ds
             //for elem in {elemIn, elemOut}
-            Base::A->entry(elemInIdxSLE + i + 1, iElemIdxSLE) -= couplingUpdate2;
-            Base::A->entry(iElemIdxSLE, elemInIdxSLE + i + 1) -= couplingUpdate2;
-            Base::A->entry(elemOutIdxSLE + i + 1, iElemIdxSLE) -= couplingUpdate2;
-            Base::A->entry(iElemIdxSLE, elemOutIdxSLE + i + 1) -= couplingUpdate2;
+            Base::A->entry(elemInIdxSLE + i + 1, iElemIdxSLE) -=
+              couplingUpdate2;
+            Base::A->entry(iElemIdxSLE, elemInIdxSLE + i + 1) -=
+              couplingUpdate2;
+            Base::A->entry(elemOutIdxSLE + i + 1, iElemIdxSLE) -=
+              couplingUpdate2;
+            Base::A->entry(iElemIdxSLE, elemOutIdxSLE + i + 1) -=
+              couplingUpdate2;
           }
         };
 
@@ -470,8 +480,10 @@ private:
             //quadrature for
             // int_iElem (Kpar * grad(d * phi_iElem,i))
             //    * grad(d * phi,iElem,0) ds
-            Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) += interfaceUpdate3;
-            Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) += interfaceUpdate3;
+            Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
+              interfaceUpdate3;
+            Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
+              interfaceUpdate3;
 
             //quadrature for
             // int_iElem (Kpar * grad(d * phi_iElem,i))

dune/mmdg/problems/mmdgproblem6.hh

@@ -17,7 +17,7 @@ public:
   //the exact bulk solution at position pos
   Scalar exactSolution (const Coordinate& pos) const
   {
-    Scalar xArg = (pos[0] > 0.5) ? sqrt(2.0) * (2.0 * pos[0] - 1.0)
+    Scalar xArg = (pos[0] < 0.5) ? sqrt(2.0) * (2.0 * pos[0] - 1.0)
       : sqrt(2.0) * (2.0 * pos[0] + 1.0);
 
     return std::cos(4.0 * pos[1]) * std::cosh(xArg);
@@ -26,14 +26,7 @@ public:
   //the exact solution on the interface at position pos
   Scalar exactInterfaceSolution (const Coordinate& pos) const
   {
-    constexpr Scalar coshsqrt2 = std::cosh(sqrt(2.0));
-
-    return ( coshsqrt2 * coshsqrt2 - 0.5 * std::sinh(2.0 * sqrt(2.0))
-      * std::tanh(2.0 * sqrt(2.0)) ) * std::cos(4.0 * pos[1]);
-
-    return 0.5 * coshsqrt2 * coshsqrt2 * (2.0 - sqrt(2.0)) * std::cos(4.0 * pos[1]);
-
-    return ( coshsqrt2 * coshsqrt2 - 1.0 ) * std::cos(4.0 * pos[1]);
+    return std::cos(4.0 * pos[1]);
   }
 
   //indicates whether an exact solution is implemented for the problem
@@ -45,33 +38,20 @@ public:
   //source term at position pos
   Scalar q (const Coordinate& pos) const
   {
-    Scalar xArg = (pos[0] > 0.5) ? sqrt(2.0) * (2.0 * pos[0] - 1.0)
+    Scalar xArg = (pos[0] < 0.5) ? sqrt(2.0) * (2.0 * pos[0] - 1.0)
       : sqrt(2.0) * (2.0 * pos[0] + 1.0);
 
     return 8.0 * std::cos(4.0 * pos[1]) * std::cosh(xArg);
-    return 8.0 * std::cos(4.0 * pos[1]) * (2.0 * std::cosh(xArg)
-      - std::sinh(xArg));
   }
 
   //interface source term at position pos
   Scalar qInterface (const Coordinate& pos) const
   {
-    constexpr Scalar coshsqrt2 = std::cosh(sqrt(2.0));
-
-    return -16.0 * d0_ * d0_ * ( coshsqrt2 * coshsqrt2 - 0.5 * std::sinh(2.0 * sqrt(2.0))
-      * std::tanh(2.0 * sqrt(2.0)) ) * std::exp(-8.0 * pos[1]) * ( std::cos(4.0 * pos[1])
-    + 3.0 * std::sin(4.0 * pos[1])) - 2.0 * sqrt(2.0) * std::sinh(2.0 * sqrt(2.0)) * std::cos(4.0 * pos[1]);
-
-
-    return 8.0 * d0_ * d0_ * coshsqrt2 * coshsqrt2 * (sqrt(2.0) - 2.0)
-      * std::exp(-8.0 * pos[1]) * ( std::cos(4.0 * pos[1])
-      + 3.0 * std::sin(4.0 * pos[1]) )
-      - 2.0 * sqrt(2.0) * std::sinh(2.0 * sqrt(2.0)) * std::cos(4.0 * pos[1]);
+    const Scalar cos4y = std::cos(4.0 * pos[1]);
+    constexpr Scalar sqrt2 = 2.0 * sqrt(2.0);
 
-    return 16.0 * d0_ * d0_ * ( 1.0 - coshsqrt2 * coshsqrt2 )
-      * std::exp(-8.0 * pos[1]) * ( std::cos(4.0 * pos[1])
-      + 3.0 * std::sin(4.0 * pos[1]) )
-      - 2.0 * sqrt(2.0) * std::sinh(2.0 * sqrt(2.0)) * std::cos(4.0 * pos[1]);
+    return -48.0 * d0_ * d0_ * std::exp(-8.0 * pos[1]) * ( cos4y
+      + 3.0 * std::sin(4.0 * pos[1]) ) - sqrt2 * std::sinh(sqrt2) * cos4y;
   }
 
   //aperture d of the fracture at position pos
@@ -93,7 +73,6 @@ public:
   Scalar Kperp (const Coordinate& pos) const
   {
     return sqrt(2.0) / std::tanh(sqrt(2.0));
-    return 2.0 * sqrt(2.0) * std::sinh(2.0 * sqrt(2.0));
   }
 
   //returns the recommended quadrature order to compute an integral