diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index f208821357a1d5c628ab74e1f79542738070aa9e..67227159f95b204edecff0331c6507ee45f3715f 100644
--- a/dune/mmdg/dg.hh
+++ b/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;
}
diff --git a/dune/mmdg/problems/dgproblem.hh b/dune/mmdg/problems/dgproblem.hh
index 8f62e78740576005225a8ad1997e44b657afb999..f5374154a202163a44653160715ac1a10af845a3 100644
--- a/dune/mmdg/problems/dgproblem.hh
+++ b/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:
diff --git a/dune/mmdg/problems/poisson.hh b/dune/mmdg/problems/poisson.hh
new file mode 100644
index 0000000000000000000000000000000000000000..516a7c014a352181aa4ca96020b2eebc7449d3ec
--- /dev/null
+++ b/dune/mmdg/problems/poisson.hh
@@ -0,0 +1,67 @@
+#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
diff --git a/src/dg.cc b/src/dg.cc
index 945503e3a4571f9366b2c7cf15a35ee0126cb576..e087fd03ebec52cf0af6a3a6be49fa97cf5b8e92 100644
--- a/src/dg.cc
+++ b/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);
diff --git a/src/parameterDG.ini b/src/parameterDG.ini
index 860cc1eb24e9168e7f663903cda9173e31ae048a..808b309c788e7fb01abf6fe8127d74a91f0a991e 100644
--- a/src/parameterDG.ini
+++ b/src/parameterDG.ini
@@ -1,2 +1,3 @@
K = 1
mu = 1000
+problem = poisson