From e1f94d765a0068880dbf63a17726ca42da54ccd2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
 <maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Sat, 18 Jan 2020 13:29:24 +0100
Subject: [PATCH] add quadrature order to problem, change gitignore

---
 .gitignore                      |  1 +
 dune/mmdg/dg.hh                 | 19 ++++++++-----------
 dune/mmdg/problems/dgproblem.hh |  9 ++++++++-
 src/dg.cc                       |  3 +++
 4 files changed, 20 insertions(+), 12 deletions(-)

diff --git a/.gitignore b/.gitignore
index d756843..6a99568 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,2 @@
 build-cmake
+.fuse_hidden*
diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 4ee978c..395ca14 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -57,9 +57,6 @@ private:
   //assemble stiffness matrix A and load vector b
   void assembleSLE (const Scalar K, const Scalar mu)
   {
-    //NOTE:
-    //const int order = 3; //order of the quadrature rule
-
     //we use the basis
     // phi_elem,0 (x) = indicator(elem);
     // phi_elem,i (x) = x[i]*indicator(elem);
@@ -77,29 +74,29 @@ private:
       //basis function phi_elem,i
       const int elemIdxSLE = (dim + 1)*elemIdx;
 
-  /*    //TODO: can be done outside of the loop?
+/*      //TODO: can be done outside of the loop?
       const Dune::QuadratureRule<double,dim>& rule =
-        Dune::QuadratureRules<double,dim>::rule(geo.type(),order);
+        Dune::QuadratureRules<double,dim>::rule(geo.type(), problem_.quadratureOrder());
+      Dune::FieldVector<double, dim+1> update(0.0);
 
       //NOTE: how are quadrature rules in Dune applied correctly?
       for (const auto& ip : rule)
       {
-        const auto qp = ip.position();
+        const auto& qp = ip.position();
         //NOTE: is the volume of the element taken into account
         //automatically?
-        const auto weight = ip.weight();
+        const double weight = ip.weight();
 
         //quadrature for int_elem q*phi_elem,0 dV
-        b[elemIdxSLE] += weight * q(qp);
+        b[elemIdxSLE] += weight * problem_.q(qp);
 
         //quadrature for int_elem q*phi_elem,i dV
         for (int i = 0; i < dim; i++)
         {
-          b[elemIdxSLE + i + 1] += weight * qp[i] * q(qp);
+          b[elemIdxSLE + i + 1] += weight * qp[i] * problem_.q(qp);
         }
       }
-  */
-
+*/
       //NOTE: makeshift solution for source term q = -1
       b[elemIdxSLE] += -elemVol;
       for (int i = 0; i < dim; i++)
diff --git a/dune/mmdg/problems/dgproblem.hh b/dune/mmdg/problems/dgproblem.hh
index f537415..c7ffb3d 100644
--- a/dune/mmdg/problems/dgproblem.hh
+++ b/dune/mmdg/problems/dgproblem.hh
@@ -17,7 +17,7 @@ class DGProblem
     }
 
     //indicates whether an exact solution is implemented for the problem
-    virtual bool hasExactSolution() const
+    virtual bool hasExactSolution () const
     {
       return false;
     };
@@ -33,6 +33,13 @@ class DGProblem
     {
       return Scalar(0.0);
     }
+
+    //returns the recommended quadrature order to compute an integral
+    //over x * q(x)
+    virtual int quadratureOrder () const
+    {
+      return 1;
+    }
 };
 
 #endif
diff --git a/src/dg.cc b/src/dg.cc
index e087fd0..8a50902 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -23,6 +23,9 @@
 #include <dune/mmdg/problems/dgproblem.hh>
 #include <dune/mmdg/problems/poisson.hh>
 
+
+//TODO: specify K in problem? exact solution depends on K!
+
 int main(int argc, char** argv)
 {
   try
-- 
GitLab