add exact 3d solution to problem dg3

README.md

@@ -23,7 +23,6 @@ Dependencies
 ------------
 
 `dune-mmdg` requires the DUNE core modules, version 2.7 or later, and a current version of the `dune-mmesh` module.
-
 Additionally, for the `dune-mmesh` module, you will need a current version of `CGAL` installed on your system.
 
 

dune/mmdg/problems/dgproblem3.hh

@@ -13,7 +13,9 @@ class DGProblem3 : public DGProblem<Coordinate, Scalar>
   public:
     static constexpr int dim = Coordinate::dimension;
 
-    //the exact solution at position pos
+    //the exact solution at position pos,
+    //in 2D and 3D the exact solution is obtained by a Fourier series expansion
+    //of the rectangular function
     Scalar exactSolution (const Coordinate& pos) const
     {
       if constexpr (dim == 1)
@@ -21,19 +23,42 @@ class DGProblem3 : public DGProblem<Coordinate, Scalar>
         return 0.5 * pos[0] * (pos[0] - 1);
       }
 
-      if constexpr (dim == 2) //approximate series
+      if constexpr (dim == 2) //Fourier series with cut-off
       {
-        const int cutoff = 21;
+        constexpr int cutoff = 31;
         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));
+              / (m * n * (m * m + n * n));
           }
         }
-        solution *= -16 / pow(M_PI, 4);
+
+        solution *= -16.0 / pow(M_PI, 4);
+        return solution;
+      }
+
+      if constexpr (dim == 3) //Fourier series with cut-off
+      {
+        constexpr int cutoff = 21;
+        double solution = 0.0;
+
+        for (int k = 1; k <= cutoff; k = k + 2)
+        {
+          for (int l = 1; l <= cutoff; l = l + 2)
+          {
+            for (int m = 1; m <= cutoff; m = m + 2)
+            {
+              solution += sin(M_PI * k * pos[0]) * sin(M_PI * l * pos[1])
+                * sin(M_PI * m* pos[2]) / (k * l * m * (k * k + l * l + m * m));
+            }
+          }
+        }
+
+        solution *= -64.0 / pow(M_PI, 5);
         return solution;
       }
 
@@ -42,10 +67,10 @@ class DGProblem3 : public DGProblem<Coordinate, Scalar>
         std::to_string(dim) + "!");
     }
 
-    //exact solution is implemented for dim = 1, 2
+    //exact solution is implemented
     bool hasExactSolution() const
     {
-      return (dim == 1 || dim == 2);
+      return true;
     };
 
     //source term

src/dg.cc

@@ -77,7 +77,7 @@ int main(int argc, char** argv)
     {
       problem = new DGProblem3<Coordinate>();
 
-      if constexpr (dim == 2)
+      if constexpr (dim >= 2)
       {
         std::cout << "NOTE: The exact solution of your problem is given by "
           "an infinite series; increase the cut-off if necessary!\n\n";