From 6318caa8ff7265657ebbf8e88e8141cccef40a10 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Tue, 27 Jan 2015 20:33:36 +0000
Subject: [PATCH] Use new differentiable Python functions to also compute the
H^1 error
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
Kudos to Carsten Gräser for those new Python functions.
[[Imported from SVN: r10028]]
---
...monicmaps-stereographic-measurement.parset | 5 ++---
src/harmonicmaps.cc | 20 ++++++++++++++-----
2 files changed, 17 insertions(+), 8 deletions(-)
diff --git a/src/harmonicmaps-stereographic-measurement.parset b/src/harmonicmaps-stereographic-measurement.parset
index 272ff5cd..294942fa 100644
--- a/src/harmonicmaps-stereographic-measurement.parset
+++ b/src/harmonicmaps-stereographic-measurement.parset
@@ -61,6 +61,5 @@ initialIterate = "[2*x[0] / (x[0]*x[0]+x[1]*x[1]+1), 2*x[1] / (x[0]*x[0]+x[1]*x[
# none / analytical / gridfunction
discretizationErrorMode = analytical
-# Again, the inverse stereographic projection
-referenceSolution = "[2*x[0] / (x[0]*x[0]+x[1]*x[1]+1), 2*x[1] / (x[0]*x[0]+x[1]*x[1]+1), (x[0]*x[0]+x[1]*x[1]-1)/ (x[0]*x[0]+x[1]*x[1]+1)]"
-
+# The python file implementing the reference solution and its derivative
+referenceSolution = inverse-stereographic-projection
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 3187cf4a..52bd3b8e 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -238,10 +238,13 @@ int main (int argc, char *argv[]) try
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));
+ // Read reference solution and its derivative into a PythonFunction
+ typedef VirtualDifferentiableFunction<FieldVector<double, dim>, TargetSpace::CoordinateType> FBase;
+ Python::Module module = Python::import(parameterSet.get<std::string>("referenceSolution"));
+ auto referenceSolution = module.get("fdf").toC<std::shared_ptr<FBase>>();
+
+ // The numerical solution, as a grid function
GFE::EmbeddedGlobalGFEFunction<FEBasis, TargetSpace> numericalSolution(feBasis, x);
// QuadratureRule for the integral of the L^2 error
@@ -249,11 +252,18 @@ int main (int argc, char *argv[]) try
// Compute the embedded L^2 error
double l2Error = DiscretizationError<GridType::LeafGridView>::computeL2Error(&numericalSolution,
- &pythonReferenceSolution,
+ referenceSolution.get(),
quadKey);
- std::cout << "L^2 error: " << l2Error << std::endl;
+ // Compute the embedded H^1 error
+ double h1Error = DiscretizationError<GridType::LeafGridView>::computeH1HalfNormDifferenceSquared(grid->leafGridView(),
+ &numericalSolution,
+ referenceSolution.get(),
+ quadKey);
+ std::cout << "L^2 error: " << l2Error
+ << " ";
+ std::cout << "H^1 error: " << std::sqrt(l2Error*l2Error + h1Error) << std::endl;
}
--
GitLab