add problem poisson, move exact solution to problem, add parameter problem in ini file

dune/mmdg/dg.hh

@@ -40,40 +40,11 @@ public:
     //NOTE: what would be an appropiate solver here?
     A.solve(d, b);
 
-    for (int i = 0; i < dof; i++)
-      for (int j = 0; j < i; j++)
-        /*assert*/if(std::abs(A[i][j] - A[j][i]) >
-          std::numeric_limits<Scalar>::epsilon())
-            std::cout << i << ", " << j << std::endl;
-
-    //NOTE: dim = 1
-//    auto exactSolution = [](const auto& pos) {return 0.5*pos[0]*(pos[0] - 1);};
-
-    //NOTE: dim = 2
-    auto exactSolution = [](const auto& pos)
-      {
-        double solution = 0.0;
-        for (int m = 1; m <= 21; m = m + 2)
-        {
-          for (int n = 1; n <= 21; 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;
-      };
-
-
-
-
     //storage for pressure data, for each element we store the
     //pressure at the corners of the element, the VTKFunction will
     //create the output using a linear interpolation for each element
-    Dune::DynamicVector<Scalar> pressure((dim + 1)*gridView_.size(0), 0.0);
-    Dune::DynamicVector<double>
-      exactPressure((dim + 1)*gridView_.size(0), 0.0);
+    Dune::DynamicVector<Scalar> pressure(dof, 0.0);
+    Dune::DynamicVector<double> exactPressure(dof, 0.0);
 
     for (const auto& elem : elements(gridView_))
     {
@@ -82,7 +53,7 @@ public:
 
       for (int k = 0; k < geo.corners(); k++)
       {
-        exactPressure[elemIdxSLE + k] = exactSolution(geo.corner(k));
+        exactPressure[elemIdxSLE + k] = problem_.exactSolution(geo.corner(k));
 
         //contribution of the basis function
         // phi_elem,0 (x) = indicator(elem);
@@ -401,6 +372,13 @@ private:
         }
       }
     }
+
+    //NOTE: check if A is symmetric
+    for (int i = 0; i < dof; i++)
+      for (int j = 0; j < i; j++)
+        /*assert*/if(std::abs(A[i][j] - A[j][i]) >
+          std::numeric_limits<Scalar>::epsilon())
+            std::cout << i << ", " << j << std::endl;
   }
 
 

dune/mmdg/problems/dgproblem.hh

@@ -3,7 +3,7 @@
 
 //abstract class providing an interface for a problem that is to be solved with
 //a discontinuous Galerkin scheme
-template<class Coordinate, class ctype>
+template<class Coordinate, class ctype = double>
 class DGProblem
 {
   public:

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 ctype = double>
+class Poisson : public DGProblem<Dune::FieldVector<ctype, dim>, ctype>
+{
+  public:
+    using Scalar = ctype;
+    using Vector = Dune::FieldVector<Scalar, dim>;
+
+    //the exact solution at position pos
+    virtual Scalar exactSolution (const Vector& pos) const
+    {
+      if (dim == 1)
+      {
+        return 0.5*pos[0]*(pos[0] - 1);
+      }
+
+      if (dim == 2)
+      {
+        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;
+      }
+
+      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 Vector& pos) const
+    {
+      return Scalar(-1.0);
+    }
+
+    //NOTE: necessary or already default?
+    //homogeneous Dirichlet boundary conditions
+    virtual Scalar boundary (const Vector& pos) const
+    {
+      return Scalar(0.0);
+    }
+};
+
+#endif

src/dg.cc

@@ -21,6 +21,7 @@
 #include <dune/mmdg/dg.hh>
 #include <dune/mmdg/clockhelper.hh>
 #include <dune/mmdg/problems/dgproblem.hh>
+#include <dune/mmdg/problems/poisson.hh>
 
 int main(int argc, char** argv)
 {
@@ -39,6 +40,7 @@ int main(int argc, char** argv)
     using GridFactory = Dune::DGFGridFactory<Grid>;
     using GridView = typename Grid::LeafGridView;
     using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
+    using Coordinate = Dune::FieldVector<double, dim>;
 
     //Message Passing Interface (MPI) for parallel computation
     Dune::MPIHelper::instance(argc, argv);
@@ -62,9 +64,23 @@ int main(int argc, char** argv)
     //element mapper for the leaf grid view
     const Mapper mapper(gridView, Dune::mcmgElementLayout());
 
-    const DGProblem<Dune::FieldVector<double, dim>, double> dummyProblem;
+//    const DGProblem<Dune::FieldVector<double, dim>, double> dummyProblem;
 
-    const DG dg(gridView, mapper, dummyProblem);
+    DGProblem<Coordinate>* problem; //problem to be solved
+
+    //determine problem type
+    if (pt["problem"] == "poisson")
+    {
+      problem = new Poisson<dim>();
+    }
+    else
+    {
+      DUNE_THROW(Dune::Exception,
+        "Invalid problem type or parameter 'problem' not specified in file "
+        "parameterDG.ini (use 'poisson' or TODO)");
+    }
+
+    const DG dg(gridView, mapper, *problem);
 
     dg(K, mu);
 

src/parameterDG.ini

 K = 1
 mu = 1000
+problem = poisson