diff --git a/dune/mmdg/problems/mmdgproblem7.hh b/dune/mmdg/problems/mmdgproblem7.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8477a75b04bc3b8e11073db2fb87ec4c0285238c
--- /dev/null
+++ b/dune/mmdg/problems/mmdgproblem7.hh
@@ -0,0 +1,198 @@
+#ifndef __DUNE_MMDG_MMDGPROBLEM7_HH__
+#define __DUNE_MMDG_MMDGPROBLEM7_HH__
+
+#include <cmath>
+#include <dune/mmdg/problems/coupleddgproblem.hh>
+
+//a coupled dg problem
+template<class Coordinate, class Scalar = double>
+class MMDGProblem7 : public CoupledDGProblem<Coordinate, Scalar>
+{
+public:
+ static constexpr int dim = Coordinate::dimension;
+ using Base = CoupledDGProblem<Coordinate, Scalar>;
+ using Matrix = typename Base::Matrix;
+
+ static constexpr Scalar dmin = 0.001; //minimum aperture of the fracture
+ static constexpr Scalar dmax = 0.002; //maximum aperture of the fracture
+
+ //constructor
+ MMDGProblem7 (const Scalar d0) : d0_(d0), xi_(1.0) {}
+
+ //constructor
+ MMDGProblem7 (const Scalar d0, const Scalar xi) : d0_(d0), xi_(xi) {}
+
+ //the exact bulk solution at position pos
+ Scalar exactSolution (const Coordinate& pos) const
+ {
+ const Scalar d = aperture(pos);
+
+ if ( pos[0] < 0.5 )
+ { //pos in Omega_1 (left of the fracture)
+ return 2.0 * ( pos[0] + d ) * exp(-d);
+ }
+
+ //pos in Omega_2 (right of the fracture)
+ return 2.0 * pos[0] * exp(d);
+ }
+
+ //the exact solution on the interface at position pos
+ Scalar exactInterfaceSolution (const Coordinate& pos) const
+ {
+
+ const Scalar d = aperture(pos);
+ return ( 1.0 + d ) / d * std::sinh(d);
+ }
+
+ //indicates whether an exact solution is implemented for the problem
+ bool hasExactSolution () const
+ {
+ return true;
+ };
+
+ //source term at position pos
+ Scalar q (const Coordinate& pos) const
+ {
+ const Scalar d = aperture(pos);
+ const Scalar dPrime2 = gradAperture(pos) * gradAperture(pos);
+ const Scalar dPrimePrime = aperturePrimePrime(pos);
+
+ if ( pos[0] < 0.5 )
+ { //pos in Omega_1 (left of the fracture)
+ return -2.0 * exp(-d) / ( 1.0 - d ) * ( ( 1.0 - pos[0] - d )
+ * ( dPrimePrime + dPrime2 * d / ( 1.0 - d ) ) - dPrime2 );
+ }
+
+ //pos in Omega_2 (right of the fracture)
+ return -2.0 * pos[0] * exp(d) / ( 1.0 + d )
+ * ( dPrimePrime + dPrime2 * d / ( 1.0 + d ) );
+
+ }
+
+ //interface source term at position pos
+ Scalar qInterface (const Coordinate& pos) const
+ {
+ return -( 4.0 + aperturePrimePrime(pos) ) * std::sinh(aperture(pos)) * aperture(pos);
+ }
+
+ //aperture d of the fracture at position pos
+ Scalar aperture (const Coordinate& pos) const
+ {
+ return dmin + 0.5 * ( dmax - dmin ) *
+ ( 1.0 + std::cos(8.0 * M_PI * pos[1]) );
+ }
+
+ //tangential gradient of the aperture d of the fracture at position pos
+ Coordinate gradAperture (const Coordinate& pos) const
+ {
+ Coordinate gradD(0.0);
+ gradD[1] = -4.0 * M_PI * ( dmax - dmin ) * std::sin(8.0 * M_PI * pos[1]);
+ return gradD;
+ }
+
+ //bulk permeability tensor at position pos
+ Matrix K (const Coordinate& pos) const
+ {
+ const Scalar d = aperture(pos);
+
+ Matrix permeability(0.0);
+
+ if ( pos[0] < 0.5 )
+ { //pos in Omega_1 (left of the fracture)
+ permeability[0][0] = 1.0;
+
+ for (int i = 1; i < dim; i++)
+ {
+ permeability[i][i] = 1.0 / ( 1.0 - d );
+ }
+ }
+ else
+ { //pos in Omega_2 (right of the fracture)
+ permeability[0][0] = 1.0;
+
+ for (int i = 1; i < dim; i++)
+ {
+ permeability[i][i] = 1.0 / ( 1.0 + d );
+ }
+ }
+
+ return permeability;
+ }
+
+ //permeability per aperture of the fracture in normal direction
+ //at position pos
+ Scalar Kperp (const Coordinate& pos) const
+ {
+ const Scalar d = aperture(pos);
+
+ return 1.0 / ( d * ( 1.0 + d ) );
+ }
+
+ //tangential permeability tensor of the interface at position pos
+ Matrix Kparallel (const Coordinate& pos) const
+ {
+ Matrix permeability(0.0);
+
+ for (int i = 0; i < dim; i++)
+ {
+ permeability[i][i] = 1.0;// / aperture(pos);
+ }
+
+ return permeability;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * boundary(x)
+ int quadratureOrderBoundary () const
+ {
+ return 10;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * d^2 * interfaceBoundary(x)
+ //and x * d * interfaceBoundary(x) * (Kpar * grad(d))
+ //and d^2 * interfaceBoundary(x) * Kpar
+ int quadratureOrderInterfaceBoundary () const
+ {
+ return 10;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * qInterface(x)
+ int quadratureOrder_qInterface () const
+ {
+ return 10;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * q(x)
+ int quadratureOrder_q () const
+ {
+ return 10;
+ }
+
+ //returns the recommended quadrature order to compute integrals
+ //over x^2 * (Kpar * grad(d) * grad(d))
+ //and x^2 * d^2
+ //and x^2 * d * (Kpar * grad(d))
+ //and x * d^2 * Kpar
+ int quadratureOrder_Kparallel () const
+ {
+ return 10;
+ }
+
+private:
+ Scalar d0_; //prefactor of the aperture
+ Scalar xi_; //coupling constant
+
+ //second derivative of aperture at position pos
+ const Scalar aperturePrimePrime (const Coordinate& pos) const
+ {
+ return -32.0 * M_PI * M_PI * ( dmax - dmin )
+ * std::cos(8.0 * M_PI * pos[1]);
+ }
+
+
+};
+
+#endif
diff --git a/src/mmdg.cc b/src/mmdg.cc
index 7a2bca12179fef4e816757b7f015b7776fc5d5c0..244d92f275eae709d490cd1088b5ee527d318bf4 100644
--- a/src/mmdg.cc
+++ b/src/mmdg.cc
@@ -26,6 +26,7 @@
#include <dune/mmdg/problems/mmdgproblem4.hh>
#include <dune/mmdg/problems/mmdgproblem5.hh>
#include <dune/mmdg/problems/mmdgproblem6.hh>
+#include <dune/mmdg/problems/mmdgproblem7.hh>
int main(int argc, char** argv)
{
@@ -101,6 +102,12 @@ int main(int argc, char** argv)
problem = new MMDGProblem6<Coordinate>(aperture, xi);
gridType = "mmdg2_C";
}
+ else if (pt["problem"] == "mmdg7")
+ {
+ xi = (pt.hasKey("xi")) ? std::stod(pt["xi"]) : 1.0;
+ problem = new MMDGProblem7<Coordinate>(aperture, xi);
+ gridType = "mmdg2_C";
+ }
else
{
DUNE_THROW(Dune::Exception,