add problem mmdg6, allow arbitrary aperture for problem mmdg2, add default values for xi

dune/mmdg/problems/coupleddgproblem.hh

@@ -32,7 +32,7 @@ public:
     return Scalar(1.0);
   }
 
-  //tangential gradient of the aperture d of the fracture at position pos 
+  //tangential gradient of the aperture d of the fracture at position pos
   virtual Coordinate gradAperture (const Coordinate& pos) const
   {
     return Coordinate(0.0);
@@ -73,7 +73,7 @@ public:
   //returns the recommended quadrature order to compute an integral
   //over x * interfaceBoundary(x)
   virtual int quadratureOrderInterfaceBoundary () const
-  { //NOTE: only relevant for 3d, not implemented!
+  {
     return 1;
   }
 

dune/mmdg/problems/mmdgproblem2.hh

@@ -14,6 +14,9 @@ public:
   using Base = CoupledDGProblem<Coordinate, Scalar>;
   using Matrix = typename Base::Matrix;
 
+  //constructor
+  MMDGProblem2 (const Scalar d) : d_(d) {}
+
   //the exact bulk solution at position pos
   Scalar exactSolution (const Coordinate& pos) const
   {
@@ -47,15 +50,15 @@ public:
   //interface source term at position pos
   Scalar qInterface (const Coordinate& pos) const
   {
-    constexpr Scalar sourcefactor =
-      (std::cos(2.0) + std::sin(2.0)) * (4.0 + 3.0 / 64.0 * M_PI * M_PI);
-    return sourcefactor * std::cos(M_PI * pos[1]);
+    constexpr Scalar sourcefactor = std::cos(2.0) + std::sin(2.0);
+    return sourcefactor * (4.0 + 0.75 * d_ * d_ * M_PI * M_PI)
+      * std::cos(M_PI * pos[1]);
   }
 
   //aperture d of the fracture at position pos
   Scalar aperture (const Coordinate& pos) const
   {
-    return Scalar(0.25);
+    return d_;
   }
 
   //permeability per aperture  of the fracture per aperture in normal direction
@@ -82,7 +85,7 @@ public:
   //returns the recommended quadrature order to compute an integral
   //over x * interfaceBoundary(x)
   int quadratureOrderInterfaceBoundary () const
-  { //NOTE: only relevant for 3d, not implemented!
+  {
     return 10;
   }
 
@@ -92,6 +95,9 @@ public:
   {
     return 10;
   }
+
+private:
+  const Scalar d_; //aperture
 };
 
 #endif

dune/mmdg/problems/mmdgproblem6.hh

+#ifndef __DUNE_MMDG_MMDGPROBLEM6_HH__
+#define __DUNE_MMDG_MMDGPROBLEM6_HH__
+
+#include <cmath>
+#include <dune/mmdg/problems/coupleddgproblem.hh>
+
+//a coupled dg problem
+template<class Coordinate, class Scalar = double>
+class MMDGProblem6 : public CoupledDGProblem<Coordinate, Scalar>
+{
+public:
+  static constexpr int dim = Coordinate::dimension;
+
+  //constructor
+  MMDGProblem6 (const Scalar d0) : d0_(d0) {}
+
+  //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)
+      : sqrt(2.0) * (2.0 * pos[0] + 1.0);
+
+    return std::cos(4.0 * pos[1]) * std::cosh(xArg);
+  }
+
+  //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]);
+  }
+
+  //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
+  {
+    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]);
+
+    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]);
+  }
+
+  //aperture d of the fracture at position pos
+  Scalar aperture (const Coordinate& pos) const
+  {
+    return d0_ * std::exp(-4.0 * 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 * d0_ * std::exp(-4.0 * pos[1]);
+    return gradD;
+  }
+
+  //permeability per aperture of the fracture in normal direction
+  //at position pos
+  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
+  //over x * boundary(x)
+  int quadratureOrderBoundary () const
+  {
+    return 10;
+  }
+
+  //returns the recommended quadrature order to compute an integral
+  //over x * interfaceBoundary(x)
+  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;
+  }
+
+private:
+  Scalar d0_; //prefactor of the aperture
+
+};
+
+#endif

src/mmdg.cc

@@ -25,6 +25,7 @@
 #include <dune/mmdg/problems/mmdgproblem3.hh>
 #include <dune/mmdg/problems/mmdgproblem4.hh>
 #include <dune/mmdg/problems/mmdgproblem5.hh>
+#include <dune/mmdg/problems/mmdgproblem6.hh>
 
 int main(int argc, char** argv)
 {
@@ -55,8 +56,8 @@ int main(int argc, char** argv)
 
     //assume constant penalty parameter //TODO
     const double mu0 = std::stod(pt["mu0"]);
-    const double xi = std::stod(pt["xi"]); //coupling parameter
     const double aperture = std::stod(pt["d"]);
+    double xi; //coupling parameter
 
     Problem* problem; //problem to be solved
 
@@ -67,32 +68,43 @@ int main(int argc, char** argv)
     {
       problem = new MMDGProblem1<Coordinate>(aperture);
       gridType = "mmdg2_5e-2";
+      xi = 0.75;
     }
     else if (pt["problem"] == "mmdg2")
     {
-      problem = new MMDGProblem2<Coordinate>();
+      problem = new MMDGProblem2<Coordinate>(aperture);
       gridType = "mmdg2_5e-2";
+      xi = 0.75;
     }
     else if (pt["problem"] == "mmdg3")
     {
       problem = new MMDGProblem3<Coordinate>(aperture);
       gridType = "mmdg3_5e-2";
+      xi = 1.0;
     }
     else if (pt["problem"] == "mmdg4")
     {
       problem = new MMDGProblem4<Coordinate>(aperture);
       gridType = "mmdg2_5e-2";
+      xi = 0.75;
     }
     else if (pt["problem"] == "mmdg5")
     {
       problem = new MMDGProblem5<Coordinate>(aperture);
       gridType = "mmdg5_5e-2";
+      xi = 0.75;
+    }
+    else if (pt["problem"] == "mmdg6")
+    {
+      problem = new MMDGProblem6<Coordinate>(aperture);
+      gridType = "mmdg2_5e-2";
+      xi = 1.0;
     }
     else
     {
       DUNE_THROW(Dune::Exception,
         "Invalid problem type or parameter 'problem' not specified in file "
-        "parameterDG.ini (use 'mmdg2' or TODO)");
+        "parameterDG.ini (use 'mmdg1' to 'mmdg6')");
     }
 
     if (pt.hasKey("gridfile")) //non-default grid type
@@ -100,6 +112,11 @@ int main(int argc, char** argv)
       gridType = pt["gridfile"];
     }
 
+    if (pt.hasKey("xi")) //non-default coupling parameter
+    {
+      xi = std::stod(pt["xi"]);
+    }
+
     //create a grid from .dgf file
     GridFactory gridFactory( "grids/" + gridType + "_" + std::to_string(dim)
       + "d.msh" );