add method to compute L2 error, bugfix for boundary terms

dune/mmdg/dg.hh

@@ -19,12 +19,13 @@
 
 template<class Grid, class GridView, class Mapper, class Problem>
 class DG
-{ //TODO: Ordnung (private / public)
+{
 public:
-  //using GridView = typename Grid::LeafGridView;
-  //using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
-  using Scalar = typename GridView::ctype;
   static constexpr int dim = GridView::dimension;
+  static constexpr auto simplex = Dune::GeometryTypes::simplex(dim);
+  static constexpr auto edge = Dune::GeometryTypes::simplex(dim-1);
+
+  using Scalar = typename GridView::ctype;
   using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
   using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
   using VTKFunction = Dune::VTKFunction<GridView>;
@@ -34,8 +35,6 @@ public:
   using QRulesElem = Dune::QuadratureRules<Scalar, dim>;
   using QRuleEdge = Dune::QuadratureRule<Scalar, dim-1>;
   using QRulesEdge = Dune::QuadratureRules<Scalar, dim-1>;
-  static constexpr auto simplex = Dune::GeometryTypes::simplex(dim);
-  static constexpr auto edge = Dune::GeometryTypes::simplex(dim-1);
 
   //constructor
   DG (const Grid& grid, const GridView& gridView, const Mapper& mapper, const Problem& problem) :
@@ -50,7 +49,7 @@ public:
       d = Vector(dof_);
     }
 
-  const void operator() (const Scalar mu)
+  void operator() (const Scalar mu)
   {
     //assemble stiffness matrix A and load vector b
     assembleSLE(mu);
@@ -64,6 +63,48 @@ public:
     writeVTKoutput();
   }
 
+  Scalar computeL2error(const int quadratureOrder = 5) const
+  {
+    if (!problem_.hasExactSolution())
+    {
+      DUNE_THROW(Dune::Exception, "ERROR: Problem has no exact solution or "
+        "exact solution is not implemented");
+    }
+
+    Scalar error(0.0);
+
+    for (const auto& elem : elements(gridView_))
+    {
+      const int elemIdxSLE = (dim + 1) * mapper_.index(elem);
+      const auto& geo = elem.geometry();
+      const QRuleElem& qrule = QRulesElem::rule(simplex, quadratureOrder);
+
+      //determine L2-error on element elem using a quadrature rule
+      for (const auto& ip : qrule)
+      {
+        const auto& qpLocal = ip.position();
+        const auto& qpGlobal = geo.global(qpLocal);
+
+        const Scalar quadratureFactor =
+          ip.weight() * geo.integrationElement(qpLocal);
+
+        Scalar numericalSolutionAtQP = d[elemIdxSLE];
+
+        for (int i = 0; i < dim; i++)
+        {
+          numericalSolutionAtQP += d[elemIdxSLE + i + 1] * qpGlobal[i];
+        }
+
+        const Scalar errorEvaluation =
+          numericalSolutionAtQP - problem_.exactSolution(qpGlobal);
+
+        error += quadratureFactor * errorEvaluation * errorEvaluation;
+      }
+    }
+
+    return std::sqrt(error);
+  }
+
   const Grid& grid_;
   const GridView& gridView_;
   const Mapper& mapper_; //element mapper for gridView_
@@ -388,15 +429,15 @@ protected:
             const auto& qpLocal = ip.position();
             const auto& qpGlobal = intersctGeo.global(qpLocal);
 
-            const Scalar weight(ip.weight());
-            const Scalar
-              integrationElem(intersctGeo.integrationElement(qpLocal));
+            const Scalar quadratureFactor =
+              ip.weight() * intersctGeo.integrationElement(qpLocal);
+
             const Scalar boundaryEvaluation(problem_.boundary(qpGlobal));
             const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
               problem_.K(qpGlobal);
 
             //quadrature for int_intersct mu * g * phi_elem,0 ds
-            b[elemIdxSLE] += weight * integrationElem * mu * boundaryEvaluation;
+            b[elemIdxSLE] += quadratureFactor * mu * boundaryEvaluation;
 
             //quadrature for
             // int_intersct (mu * phi_elem,i - K[:,i] * normal) * g ds
@@ -404,7 +445,7 @@ protected:
             for (int i = 0; i < dim; i++)
             {
               b[elemIdxSLE + i + 1] +=
-                weight * integrationElem * (mu * qpGlobal[i] -
+                quadratureFactor * (mu * qpGlobal[i] -
                   (KEvaluation[i] * normal)) * boundaryEvaluation;
             }
           }
@@ -417,7 +458,7 @@ protected:
             // int_intersct phi_elem,0 * K*grad(phi_elem,i) * normal ds
             //with K*grad(phi_elem,i) = K[:,i]
             A->entry(elemIdxSLE + i + 1, elemIdxSLE) +=
-              mu * linearIntegrals[i] - 0.5 * KIntegrals[i];
+              mu * linearIntegrals[i] - KIntegrals[i];
 
             for (int j = 0; j <= i; j++)
             {
@@ -429,7 +470,7 @@ protected:
                 //with K*grad(phi_elem,i) = K[:,i]
               A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
                 mu * quadraticIntregrals[i][j]
-                - 0.5 * ( KIntegralsX[i][j] + KIntegralsX[j][i] );
+                - ( KIntegralsX[i][j] + KIntegralsX[j][i] );
             }
           }
         }