modify problem mmdg4

dune/mmdg/problems/mmdgproblem4.hh

@@ -10,6 +10,7 @@ class MMDGProblem4 : public CoupledDGProblem<Coordinate, Scalar>
 {
 public:
   static constexpr int dim = Coordinate::dimension;
+  static constexpr Scalar jump = 1.0;
   using Base = CoupledDGProblem<Coordinate, Scalar>;
   using Matrix = typename Base::Matrix;
 
@@ -23,7 +24,7 @@ public:
 
     if (pos[0] > 0.5)
     {
-      return solution += d_;
+      return solution += jump;
     }
 
     return solution;
@@ -32,7 +33,7 @@ public:
   //the exact solution on the interface at position pos
   Scalar exactInterfaceSolution (const Coordinate& pos) const
   {
-    return 0.25 - pos[1] * pos[1] + 0.5 * d_;
+    return 0.25 - pos[1] * pos[1] + 0.5 * jump;
   }
 
   //indicates whether an exact solution is implemented for the problem
@@ -44,19 +45,34 @@ public:
   //interface source term at position pos
   Scalar qInterface (const Coordinate& pos) const
   {
-    return 2.0 * d_;
+    // return 2.0 * d_ * d_;
+    const Scalar de = aperture(pos);
+    const Scalar dPrime = sqrt(gradAperture(pos) * gradAperture(pos));
+    const Scalar dPrimePrime = 0.;
+
+    return (d_ * d_ * ( 10. * pos[1] * 0.01 + 9. * pos[1] * pos[1] + 2. * 0.01 * 0.01 - 0.5 * (jump + 0.5))) * de;
   }
 
   //aperture d of the fracture at position pos
   Scalar aperture (const Coordinate& pos) const
   {
-    return d_;
+    return d_ * (pos[1] + 0.01);;
+  }
+
+  //tangential gradient of the aperture d of the fracture at position pos
+  Coordinate gradAperture (const Coordinate& pos) const
+  {
+    Coordinate gradD(0.0);
+    gradD[1] = d_;
+    return gradD;
+
+    return Coordinate(0.0);
   }
 
   //permeability of the fracture in normal direction at position pos
   Scalar Kperp (const Coordinate& pos) const
   {
-    return Scalar(1.0) / d_;
+    return Scalar(1.0) / jump;
   }
 
   //tangential permeability tensor of the interface at position pos
@@ -66,7 +82,7 @@ public:
 
     for (int i = 0; i < dim; i++)
     {
-      permeability[i][i] = 1.0 / d_;
+      permeability[i][i] = 1.0 * aperture(pos);//d_;
     }
 
     return permeability;
@@ -76,7 +92,7 @@ public:
   //over x * boundary(x)
   int quadratureOrderBoundary () const
   {
-    return 3;
+    return 10;//3;
   }
 
   //returns the recommended quadrature order to compute an integral
@@ -85,7 +101,7 @@ public:
   //and   d^2 * interfaceBoundary(x) * Kpar
   int quadratureOrderInterfaceBoundary () const
   {
-    return 3;
+    return 10;//3;
   }
 
 private: