use different problem for poisson

dune/mmdg/problems/poisson.hh

@@ -9,18 +9,18 @@ class Poisson : public DGProblem<Coordinate, Scalar>
 {
   public:
     static constexpr int dim = Coordinate::dimension;
-    
+
     //the exact solution at position pos
-    //1d: f(x)     = x,
-    //2d: f(x,y)   = x*y,
-    //3d: f(x,y,z) = x*y*z
+    //1d: f(x)     = 0.5*x^2,
+    //2d: f(x,y)   = 0.5*x^2*y^2,
+    //3d: f(x,y,z) = 0.5*x^2*y^2*z^2
     Scalar exactSolution (const Coordinate& pos) const
     {
-      Scalar solution(1.0);
+      Scalar solution(0.5);
 
-      for (int i = 0; i < pos.dim(); i++)
+      for (int i = 0; i < dim; i++)
       {
-        solution *= pos[i];
+        solution *= pos[i] * pos[i];
       }
 
       return solution;
@@ -32,11 +32,43 @@ class Poisson : public DGProblem<Coordinate, Scalar>
       return true;
     };
 
+    //source term at position pos
+    Scalar q (const Coordinate& pos) const
+    {
+      if constexpr (dim == 1)
+      { //q(x) = -1
+        return Scalar(-1.0);
+      }
+
+      if constexpr (dim == 2)
+      { //q(x,y) = -x^2 - y^2
+        return -(pos[0] * pos[0] + pos[1] * pos[1]);
+      }
+
+      if constexpr (dim == 3)
+      { //q(x,y,z) = -x^2*y^2 - x^2*z^2 - y^2*z^2
+        Scalar pos1_2 = pos[1] * pos[1];
+        Scalar pos2_2 = pos[2] * pos[2];
+        return -pos[0] * pos[0] * (pos1_2 + pos2_2) - pos1_2 * pos2_2;
+      }
+
+      DUNE_THROW(Dune::Exception,
+        "Exact Solution not implemented for dimension = " +
+        std::to_string(dim) + "!");
+    }
+
+    //returns the recommended quadrature order to compute an integral
+    //over x * q(x)
+    virtual int quadratureOrder_q () const
+    {
+      return 2 * dim - 1;
+    }
+
     //returns the recommended quadrature order to compute an integral
     //over x * boundary(x) and K(x) * boundary(x)
     int quadratureOrderBoundary () const
     {
-      return dim;
+      return 2 * dim - 1;
     }
 };