add mmdg problems

dune/mmdg/problems/coupleddgproblem.hh

@@ -6,34 +6,39 @@
 //abstract class providing an interface for a coupled flow problem in a
 //fractured porous medium that is to be solved with a coupled discontinuous
 //Galerkin scheme
-template<class Vector, class Scalar = double>
-class CoupledDGProblem : public DGProblem<Vector, Scalar>
+template<class Coordinate, class Scalar = double>
+class CoupledDGProblem : public DGProblem<Coordinate, Scalar>
 {
   public:
-    static constexpr int dim = Vector::dimension;
-    using Base = DGProblem<Vector, Scalar>;
+    static constexpr int dim = Coordinate::dimension;
+    using Base = DGProblem<Coordinate, Scalar>;
     using Matrix = typename Base::Matrix;
 
     //the exact solution on the interface at position pos
-    virtual Scalar exactInterfaceSolution (const Vector& pos) const
+    virtual Scalar exactInterfaceSolution (const Coordinate& pos) const
     {
       return Scalar(0.0);
     }
 
     //interface source term at position pos
-    virtual Scalar qInterface (const Vector& pos) const
+    virtual Scalar qInterface (const Coordinate& pos) const
     {
       return Scalar(0.0);
     }
 
     //aperture d of the fracture at position pos
-    virtual Scalar aperture (const Vector& pos) const
+    virtual Scalar aperture (const Coordinate& pos) const
     {
       return Scalar(1.0);
     }
 
+    virtual Scalar interfaceBoundary(const Coordinate& pos) const
+    {
+      return exactInterfaceSolution(pos);
+    }
+
     //tangential permeability tensor of the interface at position pos
-    virtual Matrix Kparallel (const Vector& pos) const
+    virtual Matrix Kparallel (const Coordinate& pos) const
     {
       Matrix permeability(0.0);
 
@@ -47,12 +52,11 @@ class CoupledDGProblem : public DGProblem<Vector, Scalar>
     }
 
     //permeability of the fracture in normal direction at position pos
-    virtual Scalar Kperp (const Vector& pos) const
+    virtual Scalar Kperp (const Coordinate& pos) const
     {
       return Scalar(1.0);
     }
 
-
     //returns the recommended quadrature order to compute an integral
     //over x * qInterface(x)
     virtual int quadratureOrder_qInterface () const
@@ -60,6 +64,13 @@ class CoupledDGProblem : public DGProblem<Vector, Scalar>
       return 1;
     }
 
+    //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;
+    }
+
     //returns the recommended quadrature order to compute an integral
     //over Kparallel(x)
     virtual int quadratureOrder_Kparallel () const

dune/mmdg/problems/dgproblem.hh

@@ -31,7 +31,7 @@ class DGProblem
     //boundary condition at position pos
     virtual Scalar boundary (const Coordinate& pos) const
     {
-      return Scalar(0.0);
+      return exactSolution(pos);
     }
 
     //permeability tensor at position pos

dune/mmdg/problems/inhomogeneousbulkproblem.hh

@@ -32,12 +32,6 @@ class InhomogeneousBulkProblem : public DGProblem<Coordinate, Scalar>
       return pos * Coordinate(2.0);
     }
 
-    //boundary condition
-    virtual Scalar boundary (const Coordinate& pos) const
-    {
-      return exactSolution(pos);
-    }
-
     //permeability tensor at position pos
     virtual Matrix K (const Coordinate& pos) const
     {

dune/mmdg/problems/mmdgproblem1.hh

+#ifndef __DUNE_MMDG_MMDGPROBLEM1_HH__
+#define __DUNE_MMDG_MMDGPROBLEM1_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 MMDGProblem1 : 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 solution = pos * Coordinate(1.0);
+
+      if (pos[0] > 0.5)
+      {
+        solution += 1.0;
+      }
+
+      return solution;
+    }
+
+    //the exact solution on the interface at position pos
+    Scalar exactInterfaceSolution (const Coordinate& pos) const
+    {
+      return pos * Coordinate(1.0) + 0.5;
+    }
+
+    //indicates whether an exact solution is implemented for the problem
+    bool hasExactSolution () const
+    {
+      return true;
+    };
+
+    //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 * boundary(x)
+    int quadratureOrderBoundary () const
+    {
+      return 2;
+    }
+
+    //returns the recommended quadrature order to compute an integral
+    //over x * interfaceBoundary(x)
+    int quadratureOrderInterfaceBoundary () const
+    { //NOTE: only relevant for 3d, not implemented!
+      return 2;
+    }
+};
+
+#endif

dune/mmdg/problems/mmdgproblem2.hh

+#ifndef __DUNE_MMDG_MMDGPROBLEM2_HH__
+#define __DUNE_MMDG_MMDGPROBLEM2_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 MMDGProblem2 : 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 solution =
+        (pos[0] < 0.5) ? std::sin(4 * pos[0]) : std::cos(4 * pos[0]);
+      return solution * std::cos(M_PI * pos[1]);
+    }
+
+    //the exact solution on the interface at position pos
+    Scalar exactInterfaceSolution (const Coordinate& pos) const
+    {
+      constexpr Scalar cos2 = std::cos(2.0);
+      constexpr Scalar sin2 = std::sin(2.0);
+      return 0.75 * (cos2 - sin2) * std::cos(M_PI * pos[1]);
+    }
+
+    //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
+    {
+      constexpr Scalar sourcefactor = 16.0 + M_PI * M_PI;
+      Scalar source =
+        (pos[0] < 0.5) ? std::sin(4 * pos[0]) : std::cos(4 * pos[0]);
+      return sourcefactor * source * std::cos(M_PI * pos[1]);
+    }
+
+    //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 + 0.1875 * M_PI * M_PI);
+      return sourcefactor * aperture(pos) * std::cos(M_PI * pos[1]);
+    }
+
+    //aperture d of the fracture at position pos
+    Scalar aperture (const Coordinate& pos) const
+    {
+      return Scalar(0.25);
+    }
+
+    //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

dune/mmdg/problems/poisson.hh

@@ -12,7 +12,7 @@ class Poisson : public DGProblem<Coordinate, Scalar>
     //1d: f(x)     = x,
     //2d: f(x,y)   = x*y,
     //3d: f(x,y,z) = x*y*z
-    virtual Scalar exactSolution (const Coordinate& pos) const
+    Scalar exactSolution (const Coordinate& pos) const
     {
       Scalar solution(1.0);
 
@@ -29,12 +29,6 @@ class Poisson : public DGProblem<Coordinate, Scalar>
     {
       return true;
     };
-
-    //boundary condition
-    virtual Scalar boundary (const Coordinate& pos) const
-    {
-      return exactSolution(pos);
-    }
 };
 
 #endif

dune/mmdg/problems/zeropoisson.hh

@@ -15,7 +15,7 @@ class ZeroPoisson : public DGProblem<Coordinate, Scalar>
     static constexpr int dim = Coordinate::dimension;
 
     //the exact solution at position pos
-    virtual Scalar exactSolution (const Coordinate& pos) const
+    Scalar exactSolution (const Coordinate& pos) const
     {
       if constexpr (dim == 1)
       {
@@ -54,6 +54,12 @@ class ZeroPoisson : public DGProblem<Coordinate, Scalar>
     {
       return Scalar(-1.0);
     }
+
+    //boundary condition at position pos
+    virtual Scalar boundary (const Coordinate& pos) const
+    {
+      return Scalar(0.0);
+    }
 };
 
 #endif

src/grids/mmdg2.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};
+
+//points for interface x = 0.5
+Point(5) = {0.5, 0, 0, lc};
+Point(6) = {0.5, 1, 0, lc};
+Point(7) = {0.5, 0.5, 0, lc};
+
+//domain outline
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+
+//interface outline
+Line(6) = {5, 7};
+Line(7) = {7, 6};
+
+//curve loops
+Curve Loop(1) = {1:4}; //domain boundary
+
+//surfaces
+Plane Surface(1) = {1};
+Physical Surface(1) = {1};
+
+Curve{6:7} In Surface{1};
+Physical Curve(2) = {6:7};

src/mmdg.cc

@@ -5,6 +5,7 @@
 #include <sys/stat.h>
 #include <sys/types.h>
 #include <type_traits>
+#include <string>
 
 #include <dune/common/exceptions.hh>
 #include <dune/common/parametertree.hh>
@@ -19,7 +20,8 @@
 #include <dune/mmdg/mmdg.hh>
 #include <dune/mmdg/clockhelper.hh>
 #include <dune/mmdg/problems/coupleddgproblem.hh>
-
+#include <dune/mmdg/problems/mmdgproblem1.hh>
+#include <dune/mmdg/problems/mmdgproblem2.hh>
 
 int main(int argc, char** argv)
 {
@@ -52,31 +54,39 @@ int main(int argc, char** argv)
     const double mu = std::stod(pt["mu"]);
     const double xi = std::stod(pt["xi"]); //coupling parameter
 
-    //create a grid from .dgf file
-    GridFactory gridFactory( "grids/sphere" + std::to_string(dim) + "d.msh" );
-    const Grid& grid = *gridFactory.grid();
-    const GridView& gridView = grid.leafGridView();
-    const IGrid& iGrid = grid.interfaceGrid();
-    const IGridView& iGridView = iGrid.leafGridView();
-
-    //element mapper for the leaf grid views
-    const Mapper mapper(gridView, Dune::mcmgElementLayout());
-    const IMapper iMapper(iGridView, Dune::mcmgElementLayout());
-
     Problem* problem; //problem to be solved
 
+    std::string gridType;
+
     //determine problem type
-    if (pt["problem"] == "poisson")
-    { //TODO
-      problem = new CoupledDGProblem<Coordinate>();
+    if (pt["problem"] == "mmdg1")
+    {
+      problem = new MMDGProblem1<Coordinate>();
+      gridType = "mmdg2";
+    }
+    else if (pt["problem"] == "mmdg2")
+    {
+      problem = new MMDGProblem2<Coordinate>();
+      gridType = "mmdg2";
     }
     else
     {
       DUNE_THROW(Dune::Exception,
         "Invalid problem type or parameter 'problem' not specified in file "
-        "parameterDG.ini (use 'poisson' or TODO)");
+        "parameterDG.ini (use 'mmdg2' or TODO)");
     }
 
+    //create a grid from .dgf file
+    GridFactory gridFactory( "grids/" + gridType + ".msh" );
+    const Grid& grid = *gridFactory.grid();
+    const GridView& gridView = grid.leafGridView();
+    const IGrid& iGrid = grid.interfaceGrid();
+    const IGridView& iGridView = iGrid.leafGridView();
+
+    //element mapper for the leaf grid views
+    const Mapper mapper(gridView, Dune::mcmgElementLayout());
+    const IMapper iMapper(iGridView, Dune::mcmgElementLayout());
+
     MMDG mmdg(grid, gridView, mapper, iGridView, iMapper, *problem);
     mmdg(mu, xi);
 

src/parameterMMDG.ini

 mu = 1000
 xi = 0.75
-problem = poisson
+problem = mmdg1