From 67d438e80227e22d6b90ce9180365724f1be043e Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Tue, 27 Jan 2015 13:53:01 +0000
Subject: [PATCH] Prototype: measure the L^2 error of the solution against a
reference function
This is only a prototype, because the GFE solution function is replaced by
on where the interpolation is done in the embedding space. That defeats the
purpose of measuring the GFE discretization error :-) but it is easier to
implement. The real thing will follow shortly.
[[Imported from SVN: r10023]]
---
src/harmonicmaps.cc | 30 ++++++++++++++++++++++++++++++
1 file changed, 30 insertions(+)
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 45decd8e..23e6d958 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -21,6 +21,7 @@
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functions/vtkbasisgridfunction.hh>
#include <dune/fufem/functiontools/basisinterpolator.hh>
+#include <dune/fufem/discretizationerror.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
@@ -230,6 +231,35 @@ int main (int argc, char *argv[]) try
vtkWriter.write(resultPath + "_" + energy + "_result");
+ /////////////////////////////////////////////////////////////////
+ // Measure the discretization error, if requested
+ /////////////////////////////////////////////////////////////////
+
+ if (parameterSet.get<std::string>("discretizationErrorMode")=="analytical")
+ {
+ // Read reference solution into a PythonFunction
+ std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("referenceSolution") + std::string(")");
+ PythonFunction<FieldVector<double,dim>, TargetSpace::CoordinateType > pythonReferenceSolution(Python::evaluate(lambda));
+
+ std::vector<TargetSpace::CoordinateType> xEmbedded(x.size());
+ for (size_t i=0; i<x.size(); i++)
+ xEmbedded[i] = x[i].globalCoordinates();
+
+ BasisGridFunction<FEBasis,std::vector<TargetSpace::CoordinateType> > numericalSolution(feBasis, xEmbedded);
+
+ // QuadratureRule for the integral of the L^2 error
+ QuadratureRuleKey quadKey(dim,3);
+
+ // Compute the embedded L^2 error
+ double l2Error = DiscretizationError<GridType::LeafGridView>::computeL2Error(&numericalSolution,
+ &pythonReferenceSolution,
+ quadKey);
+
+ std::cout << "L^2 error: " << l2Error << std::endl;
+
+
+ }
+
} catch (Exception e) {
std::cout << e << std::endl;
--
GitLab