diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 3b567de885524ab9cfb893ababe42e12390d2fbf..c2040e3d811ba8a29363b9fcb3628431e6aa6754 100644
--- a/dune/mmdg/mmdg.hh
+++ b/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))
diff --git a/dune/mmdg/problems/mmdgproblem6.hh b/dune/mmdg/problems/mmdgproblem6.hh
index ea581e748a53399d1dbdcfdee6dd6d6cbdd57919..f7b4ce1d86cb3c62fc280f4ef8b32a681e386065 100644
--- a/dune/mmdg/problems/mmdgproblem6.hh
+++ b/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