add new problem mmdg3

dune/mmdg/problems/mmdgproblem3.hh

+#ifndef __DUNE_MMDG_MMDGPROBLEM3_HH__
+#define __DUNE_MMDG_MMDGPROBLEM3_HH__
+
+#include <cmath>
+#include <dune/mmdg/problems/coupleddgproblem.hh>
+
+//a coupled dg problem for dim = 2,
+//taken from [Antonietti et al. (2019): Example 6.3]
+template<class Coordinate, class Scalar = double>
+class MMDGProblem3 : public CoupledDGProblem<Coordinate, Scalar>
+{
+  public:
+    static constexpr int dim = Coordinate::dimension;
+
+    //the exact bulk solution at position pos
+    Scalar exactSolution (const Coordinate& pos) const
+    {
+      Scalar xPlusy = pos * Coordinate(1.0);
+      Scalar solution = std::exp(xPlusy);
+      return (xPlusy < 1.0) ? solution :
+        0.5 * solution + std::exp(1.0) * (0.5 + 1.5 * sqrt(2) * this->aperture(pos));
+    }
+
+    //the exact solution on the interface at position pos
+    Scalar exactInterfaceSolution (const Coordinate& pos) const
+    {
+      return std::exp(1.0) * (1.0 + sqrt(2) * this->aperture(pos));
+    }
+
+    //indicates whether an exact solution is implemented for the problem
+    bool hasExactSolution () const
+    {
+      return true;
+    };
+
+    //bulk source term at position pos
+    Scalar q (const Coordinate& pos) const
+    {
+      Scalar xPlusy = pos * Coordinate(1.0);
+      Scalar source = -std::exp(xPlusy);
+      return (xPlusy < 1.0) ? -2.0 * source : source;
+    }
+
+    //interface source term at position pos
+    Scalar qInterface (const Coordinate& pos) const
+    {
+      return std::exp(1.0) / sqrt(2);
+    }
+
+    // //aperture d of the fracture at position pos
+    // Scalar aperture (const Coordinate& pos) const
+    // {
+    //   return Scalar(1.0);
+    // }
+
+    //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 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
+    { //NOTE: only relevant for 3d, not implemented!
+      return 10;
+    }
+
+    //returns the recommended quadrature order to compute an integral
+    //over x * qInterface(x)
+    int quadratureOrder_qInterface () const
+    {
+      return 10;
+    }
+};
+
+#endif

src/grids/mmdg3.geo

+//characteristic length
+lc = 5e-2;
+
+//domain corners, domain = [0,1] x [0,1]
+Point(1) = {0, 0, 0, lc};
+Point(2) = {1, 0, 0, lc};
+Point(3) = {1, 1, 0, lc};
+Point(4) = {0, 1, 0, lc};
+
+//domain outline
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+
+//interface outline
+Line(5) = {2,4};
+
+//curve loops
+Curve Loop(1) = {1:4}; //domain boundary
+
+//surfaces
+Plane Surface(1) = {1};
+Physical Surface(1) = {1};
+
+Curve{5} In Surface{1};
+Physical Curve(2) = {5};

src/mmdg.cc

@@ -22,6 +22,7 @@
 #include <dune/mmdg/problems/coupleddgproblem.hh>
 #include <dune/mmdg/problems/mmdgproblem1.hh>
 #include <dune/mmdg/problems/mmdgproblem2.hh>
+#include <dune/mmdg/problems/mmdgproblem3.hh>
 
 int main(int argc, char** argv)
 {
@@ -69,6 +70,11 @@ int main(int argc, char** argv)
       problem = new MMDGProblem2<Coordinate>();
       gridType = "mmdg2";
     }
+    else if (pt["problem"] == "mmdg3")
+    {
+      problem = new MMDGProblem3<Coordinate>();
+      gridType = "mmdg3";
+    }
     else
     {
       DUNE_THROW(Dune::Exception,