add new bulk problems

dune/mmdg/problems/coupleddgproblem.hh

@@ -10,7 +10,7 @@ template<class Vector, class Scalar = double>
 class CoupledDGProblem : public DGProblem<Vector, Scalar>
 {
   public:
-    static constexpr int dimension = Vector::dimension;
+    static constexpr int dim = Vector::dimension;
     using Base = DGProblem<Vector, Scalar>;
     using Matrix = typename Base::Matrix;
 
@@ -37,7 +37,7 @@ class CoupledDGProblem : public DGProblem<Vector, Scalar>
     {
       Matrix permeability(0.0);
 
-      for (int i = 0; i < dimension; i++)
+      for (int i = 0; i < dim; i++)
       {
         permeability[i][i] = 1.0;
       }

dune/mmdg/problems/dgproblem.hh

@@ -3,15 +3,15 @@
 
 //abstract class providing an interface for a problem that is to be solved with
 //a discontinuous Galerkin scheme
-template<class Vector, class Scalar = double>
+template<class Coordinate, class Scalar = double>
 class DGProblem
 {
   public:
-    static constexpr int dimension = Vector::dimension;
-    using Matrix = Dune::FieldMatrix<Scalar, dimension, dimension>;
+    static constexpr int dim = Coordinate::dimension;
+    using Matrix = Dune::FieldMatrix<Scalar, dim, dim>;
 
     //the exact solution at position pos
-    virtual Scalar exactSolution (const Vector& pos) const
+    virtual Scalar exactSolution (const Coordinate& pos) const
     {
       return Scalar(0.0);
     }
@@ -23,23 +23,23 @@ class DGProblem
     };
 
     //source term at position pos
-    virtual Scalar q (const Vector& pos) const
+    virtual Scalar q (const Coordinate& pos) const
     {
       return Scalar(0.0);
     }
 
     //boundary condition at position pos
-    virtual Scalar boundary (const Vector& pos) const
+    virtual Scalar boundary (const Coordinate& pos) const
     {
       return Scalar(0.0);
     }
 
     //permeability tensor at position pos
-    virtual Matrix K (const Vector& pos) const
+    virtual Matrix K (const Coordinate& pos) const
     {
       Matrix permeability(0.0);
 
-      for (int i = 0; i < dimension; i++)
+      for (int i = 0; i < dim; i++)
       {
         permeability[i][i] = 1.0;
       }

dune/mmdg/problems/inhomogeneousbulkproblem.hh

+#ifndef __DUNE_MMDG_INHOMOGENEOUSBULKPROBLEM_HH__
+#define __DUNE_MMDG_INHOMOGENEOUSBULKPROBLEM_HH__
+
+#include <dune/mmdg/problems/dgproblem.hh>
+
+template<class Coordinate, class Scalar = double>
+class InhomogeneousBulkProblem : public DGProblem<Coordinate, Scalar>
+{
+  public:
+    static constexpr int dim = Coordinate::dimension;
+    using Base = DGProblem<Coordinate, Scalar>;
+    using Matrix = typename Base::Matrix;
+
+    //the exact solution at position pos
+    //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
+    {
+      return pos * Coordinate(1.0);
+    }
+
+    //exact solution is implemented
+    bool hasExactSolution() const
+    {
+      return true;
+    };
+
+    //source term at position pos
+    virtual Scalar q (const Coordinate& pos) const
+    {
+      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
+    {
+      Matrix permeability(0.0);
+
+      for (int i = 0; i < dim; i++)
+      {
+        permeability[i][i] = 1.0 + pos[i] * pos[i];
+      }
+
+      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)
+    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;
+    }
+};
+
+#endif

dune/mmdg/problems/poisson.hh

 #ifndef __DUNE_MMDG_POISSON_HH__
 #define __DUNE_MMDG_POISSON_HH__
 
-#include <cmath>
-#include <string>
-
 #include <dune/mmdg/problems/dgproblem.hh>
 
-
-//Poisson problem with right side q = -1
-template<int dim, class Scalar = double>
-class Poisson : public DGProblem<Dune::FieldVector<Scalar, dim>, Scalar>
+//a Poisson problem
+template<class Coordinate, class Scalar = double>
+class Poisson : public DGProblem<Coordinate, Scalar>
 {
   public:
-    using Vector = Dune::FieldVector<Scalar, dim>;
-
     //the exact solution at position pos
-    virtual Scalar exactSolution (const Vector& pos) const
+    //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
     {
-      if (dim == 1)
-      {
-        return 0.5*pos[0]*(pos[0] - 1);
-      }
+      Scalar solution(1.0);
 
-      if (dim == 2)
+      for (int i = 0; i < pos.dim(); i++)
       {
-        const int cutoff = 21; //NOTE: error?
-        double solution = 0.0;
-        for (int m = 1; m <= cutoff; m = m + 2)
-        {
-          for (int n = 1; n <= cutoff; n = n + 2)
-          {
-            solution += sin(M_PI * m * pos[0]) * sin(M_PI * n * pos[1])
-              / (m * n * (m*m + n*n));
-          }
-        }
-        solution *= -16 / pow(M_PI, 4);
-        return solution;
+        solution *= pos[i];
       }
 
-      DUNE_THROW(Dune::Exception,
-        "Exact Solution not implemented for dimension = " +
-        std::to_string(dim) + "!");
+      return solution;
     }
 
-    //exact solution is implemented for dim = 1, 2
+    //exact solution is implemented
     bool hasExactSolution() const
     {
-      return (dim == 1 || dim == 2);
+      return true;
     };
 
-    //source term
-    Scalar q (const Vector& pos) const
-    {
-      return Scalar(-1.0);
-    }
-
-    //NOTE: necessary or already default?
-    //homogeneous Dirichlet boundary conditions
-    virtual Scalar boundary (const Vector& pos) const
+    //boundary condition
+    virtual Scalar boundary (const Coordinate& pos) const
     {
-      return Scalar(0.0);
+      return exactSolution(pos);
     }
 };
 

dune/mmdg/problems/zeropoisson.hh

+#ifndef __DUNE_MMDG_ZEROPOISSON_HH__
+#define __DUNE_MMDG_ZEROPOISSON_HH__
+
+#include <cmath>
+#include <string>
+
+#include <dune/mmdg/problems/dgproblem.hh>
+
+
+//Poisson problem with right side q = -1
+template<class Coordinate, class Scalar = double>
+class ZeroPoisson : public DGProblem<Coordinate, Scalar>
+{
+  public:
+    static constexpr int dim = Coordinate::dimension;
+
+    //the exact solution at position pos
+    virtual Scalar exactSolution (const Coordinate& pos) const
+    {
+      if constexpr (dim == 1)
+      {
+        return 0.5*pos[0]*(pos[0] - 1);
+      }
+
+      if constexpr (dim == 2) //approximate series
+      {
+        const int cutoff = 21;
+        double solution = 0.0;
+        for (int m = 1; m <= cutoff; m = m + 2)
+        {
+          for (int n = 1; n <= cutoff; n = n + 2)
+          {
+            solution += sin(M_PI * m * pos[0]) * sin(M_PI * n * pos[1])
+              / (m * n * (m*m + n*n));
+          }
+        }
+        solution *= -16 / pow(M_PI, 4);
+        return solution;
+      }
+
+      DUNE_THROW(Dune::Exception,
+        "Exact Solution not implemented for dimension = " +
+        std::to_string(dim) + "!");
+    }
+
+    //exact solution is implemented for dim = 1, 2
+    bool hasExactSolution() const
+    {
+      return (dim == 1 || dim == 2);
+    };
+
+    //source term
+    Scalar q (const Coordinate& pos) const
+    {
+      return Scalar(-1.0);
+    }
+};
+
+#endif

src/CMakeLists.txt

@@ -6,6 +6,10 @@ add_executable("dg-2d" dg.cc)
 target_compile_definitions("dg-2d" PRIVATE GRIDDIM=2)
 target_link_dune_default_libraries("dg-2d")
 
+add_executable("dg-3d" dg.cc)
+target_compile_definitions("dg-3d" PRIVATE GRIDDIM=3)
+target_link_dune_default_libraries("dg-3d")
+
 add_executable("mmdg-2d" mmdg.cc)
 target_compile_definitions("mmdg-2d" PRIVATE GRIDDIM=2)
 target_link_dune_default_libraries("mmdg-2d")

src/dg.cc

@@ -20,8 +20,9 @@
 #include <dune/mmdg/dg.hh>
 #include <dune/mmdg/clockhelper.hh>
 #include <dune/mmdg/problems/dgproblem.hh>
+#include <dune/mmdg/problems/zeropoisson.hh>
 #include <dune/mmdg/problems/poisson.hh>
-
+#include <dune/mmdg/problems/inhomogeneousbulkproblem.hh>
 
 int main(int argc, char** argv)
 {
@@ -62,9 +63,17 @@ int main(int argc, char** argv)
     Problem* problem; //problem to be solved
 
     //determine problem type
-    if (pt["problem"] == "poisson")
+    if (pt["problem"] == "zeropoisson")
+    {
+      problem = new ZeroPoisson<Coordinate>();
+    }
+    else if (pt["problem"] == "poisson")
+    {
+      problem = new Poisson<Coordinate>();
+    }
+    else if (pt["problem"] == "inhomogeneousbulk")
     {
-      problem = new Poisson<dim>();
+      problem = new InhomogeneousBulkProblem<Coordinate>();
     }
     else
     {
@@ -81,7 +90,7 @@ int main(int argc, char** argv)
     //print total run time
     const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);
     std::cout << "Finished with a total run time of " << timeTotal <<
-      " Seconds.\n";
+      " Seconds.\nL2-error: " << dg.computeL2error() << ".\n";
 
     return 0;
   }

src/grids/cube1d.dgf

@@ -3,7 +3,7 @@ DGF
 Interval
 0	
 1	
-100	
+1000	
 #
 
 Simplex

src/grids/cube2d.dgf

@@ -3,7 +3,7 @@ DGF
 Interval
 0	0	
 1	1	
-20	20	
+100	100	
 #
 
 Simplex

src/grids/cube3d.dgf

@@ -3,7 +3,7 @@ DGF
 Interval
 0	0	0	
 1	1	1	
-20	20	20
+10	10	10
 #
 
 Simplex