add problem dg4

dune/mmdg/problems/dgproblem4.hh

+#ifndef __DUNE_MMDG_DGPROBLEM4_HH__
+#define __DUNE_MMDG_DGPROBLEM4_HH__
+
+#include <dune/mmdg/problems/dgproblem.hh>
+#include <cmath>
+
+template<class Coordinate, class Scalar = double>
+class DGProblem4 : public DGProblem<Coordinate, Scalar>
+{
+  public:
+    static constexpr int dim = Coordinate::dimension;
+    using Base = DGProblem<Coordinate, Scalar>;
+    using Matrix = typename Base::Matrix;
+
+    // static constexpr Scalar dmin_ = 0.05; //minimum aperture of the fracture
+    // static constexpr Scalar dmax_ = 0.2; //maximum aperture of the fracture
+
+    //require dim == 2
+    DGProblem4 (const Scalar dmin, const Scalar dmax)
+      : dmin_(dmin), dmax_(dmax) {}
+    //{
+      // if (dim != 2)
+      // {
+      //   DUNE_THROW(Dune::Exception,
+      //     "Problem dg4 is not implemented for dimension = " +
+      //     std::to_string(dim) + "!");
+      // }
+    //}
+
+    //the exact solution at position pos
+    Scalar exactSolution (const Coordinate& pos) const
+    {
+      const Scalar d = aperture(pos);
+
+      if ( pos[0] < 0.5 * ( 1.0 - d ) )
+      { //pos in Omega_1 (left of the fracture)
+        return 2.0 * ( pos[0] + d ) * exp(-d);
+      }
+
+      if ( pos[0] > 0.5 * ( 1.0 + d ) )
+      { //pos in Omega_2 (right of the fracture)
+        return 2.0 * pos[0] * exp(d);
+      }
+
+      //pos in Omega_R (inside the fracture)
+      return ( 1.0 + d ) * exp( 2.0 * pos[0] - 1.0 );
+    }
+
+    //exact solution is implemented
+    bool hasExactSolution() const
+    {
+      return true;
+    };
+
+    //source term at position pos
+    Scalar q (const Coordinate& pos) const
+    {
+      const Scalar d = aperture(pos);
+      Scalar dPrime2 = aperturePrime(pos);
+      dPrime2 *= dPrime2;
+      const Scalar dPrimePrime = aperturePrimePrime(pos);
+
+
+      if ( pos[0] < 0.5 * ( 1.0 - d ) )
+      { //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 );
+      }
+
+      if ( pos[0] > 0.5 * ( 1.0 + d ) )
+      { //pos in Omega_2 (right of the fracture)
+        return -2.0 * pos[0] * exp(d) / ( 1.0 + d )
+          * ( dPrimePrime + dPrime2 * d / ( 1.0 + d ) );
+      }
+
+      //pos in Omega_R (inside the fracture)
+      return
+        -exp(2.0 * pos[0] - 1) * ( 4.0 + dPrimePrime );
+    }
+
+    //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 * ( 1.0 - d ) )
+      { //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 if ( pos[0] > 0.5 * ( 1.0 + d ) )
+      { //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 );
+        }
+      }
+      else
+      { //pos in Omega_R (inside the fracture)
+        permeability[0][0] = 1.0 / ( 1.0 + d );
+
+        for (int i = 1; i < dim; i++)
+        {
+          permeability[i][i] = 1.0;
+        }
+      }
+
+      return permeability;
+    }
+
+    //returns the recommended quadrature order to compute an integral
+    //over x * q(x)
+    virtual int quadratureOrder_q () const
+    {
+      return 10;
+    }
+
+    //returns the recommended quadrature order to compute an integral
+    //over x * boundary(x) and K(x) * boundary(x)
+    virtual int quadratureOrderBoundary () const
+    {
+      return 10;
+    }
+
+    //returns the recommended quadrature order to compute an integral
+    //over an entry K[i][j] of the permeability tensor
+    virtual int quadratureOrder_K () const
+    {
+      return 10;
+    }
+
+  private:
+    //aperture of the fracture at position pos
+    const Scalar aperture (const Coordinate& pos) const
+    {
+      return dmin_ + 0.5 * ( dmax_ - dmin_ ) *
+        ( 1.0 + std::cos(8.0 * M_PI * pos[1]) );
+    }
+
+    //derivative of aperture at position pos
+    const Scalar aperturePrime (const Coordinate& pos) const
+    {
+      return -4.0 * M_PI * ( dmax_ - dmin_ ) * std::sin(8.0 * M_PI * pos[1]);
+    }
+
+    //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]);
+    }
+
+    const Scalar dmin_; //minimum aperture of the fracture
+    const Scalar dmax_; //maximum aperture of the fracture
+
+};
+
+#endif

src/dg.cc

@@ -24,6 +24,7 @@
 #include <dune/mmdg/problems/dgproblem1.hh>
 #include <dune/mmdg/problems/dgproblem2.hh>
 #include <dune/mmdg/problems/dgproblem3.hh>
+#include <dune/mmdg/problems/dgproblem4.hh>
 
 int main(int argc, char** argv)
 {
@@ -85,11 +86,26 @@ int main(int argc, char** argv)
           "an infinite series; increase the cut-off if necessary!\n\n";
       }
     }
+    else if (pt["problem"] == "dg4")
+    {
+      if ( pt.hasKey("dmin") && pt.hasKey("dmax") )
+      {
+        const double dmin = std::stod(pt["dmin"]);
+        const double dmax = std::stod(pt["dmax"]);
+        problem = new DGProblem4<Coordinate>(dmin, dmax);
+      }
+      else
+      {
+        DUNE_THROW(Dune::Exception,
+          "Problem dg4 requires the parameters 'dmin' and 'dmax'. Check the "
+          "file parameterDG.ini and make sure that they are specified.");
+      }
+    }
     else
     {
       DUNE_THROW(Dune::Exception,
         "Invalid problem type or parameter 'problem' not specified in file "
-        "parameterDG.ini (use 'dg1' to 'dg3')");
+        "parameterDG.ini (use 'dg1' to 'dg4')");
     }
 
     DG dg(gridView, mapper, *problem);

src/parameterDG.ini

 mu0 = 1000
-problem = dg2
+problem = dg4
 fastEOC = true
+dmin = 1e-1
+dmax = 3e-1