[WIP] add problem dg5

dune/mmdg/problems/dgproblem5.hh

+#ifndef __DUNE_MMDG_DGPROBLEM5_HH__
+#define __DUNE_MMDG_DGPROBLEM5_HH__
+
+#include <dune/mmdg/problems/dgproblem.hh>
+#include <cmath>
+
+template<class Coordinate, class Scalar = double>
+class DGProblem5 : 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.1; //maximum aperture of the fracture
+    static constexpr Scalar Kbulk = 1.0; //bulk permeability
+    static constexpr Scalar K_par = 3.0; //tangential fracture permeability
+    static constexpr Scalar K_perp = 0.5; //normal fracture permeability
+
+    //require dim == 2
+    DGProblem5 ()
+    {
+      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
+    {
+      if ( pos[0] < 0.5 * (1.0 - aperture(pos)) )
+      { //pos in Omega_1 (left of the fracture)
+        return -0.5 * ( pos[0] - 0.5 ) * pos[1] * aperture(pos);
+      }
+
+      if ( pos[0] > 0.5 * (1.0 + aperture(pos)) )
+      { //pos in Omega_2 (right of the fracture)
+        return 0.5 * ( pos[0] - 0.5 ) * pos[1] * aperture(pos);
+      }
+
+      //pos in Omega_R (inside the fracture)
+      return ( pos[0] - 0.5 ) * ( pos[0] - 0.5 ) * pos[1];
+    }
+
+    //exact solution is implemented
+    bool hasExactSolution() const
+    {
+      return true;
+    };
+
+    //source term at position pos
+    Scalar q (const Coordinate& pos) const
+    {
+      if ( pos[0] < 0.5 * (1.0 - aperture(pos)) )
+      { //pos in Omega_1 (left of the fracture)
+        return -4.0 * M_PI * Kbulk * ( dmax - dmin ) *
+          ( pos[0] - 0.5 ) * ( std::sin(8.0 * M_PI * pos[1])
+            + 4.0 * M_PI * pos[1] * std::cos(8.0 * M_PI * pos[1]) );
+      }
+
+      if ( pos[0] > 0.5 * (1.0 + aperture(pos)) )
+      { //pos in Omega_2 (right of the fracture)
+        return 4.0 * M_PI * Kbulk * ( dmax - dmin ) *
+          ( pos[0] - 0.5 ) * ( std::sin(8.0 * M_PI * pos[1])
+            + 4.0 * M_PI * pos[1] * std::cos(8.0 * M_PI * pos[1]) );
+      }
+
+      //pos in Omega_R (inside the fracture)
+      return -2.0 * K_perp * pos[1];
+      // return -2.0 * K_perp * pos[1] - ( pos[0] - 0.5 ) * ( pos[0] - 0.5 ) *
+      //   ( -1.0 + 4.0 * M_PI * ( dmax - dmin ) * ( 8.0 * M_PI *
+      //     std::cos(8.0 * M_PI * pos[1]) - std::sin(8.0 * M_PI * pos[1]) )
+      //     / aperture(pos) );
+    }
+
+    //permeability tensor at position pos
+    Matrix K (const Coordinate& pos) const
+    {
+      Matrix permeability(0.0);
+
+      if ( pos[0] < 0.5 * (1.0 - aperture(pos)) ||
+        pos[0] > 0.5 * (1.0 + aperture(pos)))
+      { //pos in Omega_1 or Omega_2 (outside the fracture)
+        for (int i = 0; i < dim; i++)
+        {
+          permeability[i][i] = Kbulk;
+          //  aperture(pos) / ( aperture(pos) - 4.0 * pos[1] * M_PI * (dmax - dmin) * std::sin(8.0 * M_PI * pos[1]) );
+
+          if (permeability[i][i] <= 0.0)
+          {
+            std::cout << "WARNING!" << std::endl;
+          }
+        }
+      }
+      else
+      { //pos in Omega_R (inside the fracture)
+        permeability[0][0] = K_perp;
+
+        for (int i = 1; i < dim; i++)
+        {
+          permeability[i][i] = K_par;
+            // 1.0 - pos[1] * 4.0 * M_PI * ( dmax - dmin ) *
+            //   std::sin(8.0 * M_PI * pos[1]) / aperture(pos);
+
+          if (permeability[i][i] <= 0.0)
+          {
+            std::cout << "WARNING" << std::endl;
+          }
+        }
+      }
+
+      return permeability;
+    }
+
+  private:
+    //aperture of the fracture
+    const Scalar aperture (const Coordinate& pos) const
+    {
+      return dmin + 0.5 * ( dmax - dmin ) *
+        ( 1.0 + std::cos(8.0 * M_PI * pos[1]) );
+      //return dmin + ( dmax - dmin ) * std::cos(8.0 * M_PI * pos[1]);
+
+      // const Scalar y = pos[1] - 0.25 * floor(4.0 * pos[1]);
+      //
+      // if ( y <= 0.125 )
+      // { //0 <= y <= 0.125
+      //   return dmin + ( 1.0 - 8.0 * y ) * ( dmax - dmin );
+      // }
+      //
+      // //0.125 <= y <= 0.25
+      // return dmin + ( 8.0 * y - 1.0 ) * ( dmax - dmin );
+    }
+};
+
+#endif