From fbef589c49f2d7ca53beebbc891cfc4cf2eb9bcc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
 <maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Thu, 19 Mar 2020 12:25:00 +0100
Subject: [PATCH] add method to mmdg to compute L2 error

---
 dune/mmdg/dg.hh   |  1 +
 dune/mmdg/mmdg.hh | 45 +++++++++++++++++++++++++++++++++++++++++++++
 src/mmdg.cc       |  2 +-
 3 files changed, 47 insertions(+), 1 deletion(-)

diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 4390d21..f352e23 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -63,6 +63,7 @@ public:
     writeVTKoutput();
   }
 
+  //compute L2 error using a quadrature rule
   Scalar computeL2error(const int quadratureOrder = 5) const
   {
     if (!problem_.hasExactSolution())
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 69a48fe..62d2d0b 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -44,6 +44,51 @@ public:
     writeVTKoutput();
   }
 
+  //compute L2 error using quadrature rules
+  // norm = sqrt( ||.||^2_L2(bulk) + ||.||^2_L2(interface) )
+  Scalar computeL2error(const int quadratureOrder = 5) const
+  {
+    const Scalar bulkError = Base::computeL2error(quadratureOrder);
+    Scalar interfaceError(0.0);
+
+    for (const auto& iElem : elements(iGridView_))
+    {
+      const auto& iGeo = iElem.geometry();
+      const int iElemIdxSLE =
+        (dim + 1) * Base::gridView_.size(0) + dim * iMapper_.index(iElem);
+      const QRuleInterface& qrule =
+        QRulesInterface::rule(iSimplex, quadratureOrder);
+      const InterfaceBasis iFrame = interfaceFrame(
+        Base::grid_.asIntersection(iElem).centerUnitOuterNormal());
+
+      //determine L2-error on the interface element iElem using a
+      //quadrature rule
+      for (const auto& ip : qrule)
+      {
+        const auto& qpLocal = ip.position();
+        const auto& qpGlobal = iGeo.global(qpLocal);
+
+        const Scalar quadratureFactor =
+          ip.weight() * iGeo.integrationElement(qpLocal);
+
+        Scalar numericalSolutionAtQP = Base::d[iElemIdxSLE];
+
+        for (int i = 0; i < dim - 1; i++)
+        {
+          numericalSolutionAtQP +=
+            Base::d[iElemIdxSLE + i + 1] * (qpGlobal * iFrame[i]);
+        }
+
+        const Scalar errorEvaluation = numericalSolutionAtQP
+          - Base::problem_.exactInterfaceSolution(qpGlobal);
+
+        interfaceError += quadratureFactor * errorEvaluation * errorEvaluation;
+      }
+    }
+
+    return sqrt(bulkError * bulkError + interfaceError);
+  }
+
   const InterfaceGridView& iGridView_;
   const InterfaceMapper& iMapper_;
 
diff --git a/src/mmdg.cc b/src/mmdg.cc
index 4d567e2..e2d4348 100644
--- a/src/mmdg.cc
+++ b/src/mmdg.cc
@@ -83,7 +83,7 @@ int main(int argc, char** argv)
     //print total run time
     const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);
     std::cout << "Finished with a total run time of " << timeTotal <<
-      " Seconds.\n";
+      " Seconds.\nL2-error: " << mmdg.computeL2error() << ".\n";
 
     return 0;
   }
-- 
GitLab