From 1fc5de53e2a7f234e7b9bb1e1aece7ed62550280 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 30 Apr 2009 15:50:08 +0000
Subject: [PATCH] measure L2 and H1 errors as well
[[Imported from SVN: r4135]]
---
rod-eoc.cc | 32 +++++++++++++++++++++++++++++---
1 file changed, 29 insertions(+), 3 deletions(-)
diff --git a/rod-eoc.cc b/rod-eoc.cc
index 4c5e1f36..67b5f2ee 100644
--- a/rod-eoc.cc
+++ b/rod-eoc.cc
@@ -7,6 +7,8 @@
#include <dune/istl/io.hh>
+#include <dune/disc/miscoperators/massmatrix.hh>
+#include <dune/disc/miscoperators/laplace.hh>
#include <dune-solvers/solvers/iterativesolver.hh>
#include <dune-solvers/norms/energynorm.hh>
@@ -173,12 +175,23 @@ int main (int argc, char *argv[]) try
solve(referenceGrid, referenceSolution, numLevels, dirichletValue, parameterSet);
+ // //////////////////////////////////////////////////////////////////////
+ // Compute mass matrix and laplace matrix to emulate L2 and H1 norms
+ // //////////////////////////////////////////////////////////////////////
+
+ Dune::LeafP1Function<GridType,double> u(referenceGrid),f(referenceGrid);
+
+ Dune::MassMatrixLocalStiffness<GridType::LeafGridView,double,1> massMatrixStiffness;
+ Dune::LeafP1OperatorAssembler<GridType,double,1> massMatrix(referenceGrid);
+ massMatrix.assemble(massMatrixStiffness,u,f);
+
+ Dune::LaplaceLocalStiffness<GridType::LeafGridView,double> laplaceStiffness;
+ Dune::LeafP1OperatorAssembler<GridType,double,1> laplace(referenceGrid);
+ laplace.assemble(laplaceStiffness,u,f);
// ///////////////////////////////////////////////////////////
// Compute on all coarser levels, and compare
// ///////////////////////////////////////////////////////////
-
-
for (int i=1; i<=numLevels; i++) {
@@ -199,9 +212,22 @@ int main (int argc, char *argv[]) try
assert(referenceSolution.size() == solution.size());
+ BlockVector<TargetSpace::TangentVector> difference = computeGeodesicDifference(solution,referenceSolution);
+
+ H1SemiNorm< BlockVector<TargetSpace::TangentVector> > h1Norm(*laplace);
+ H1SemiNorm< BlockVector<TargetSpace::TangentVector> > l2Norm(*massMatrix);
+
// Compute max-norm difference
std::cout << "Level: " << i-1
- << ", max-norm error: " << computeGeodesicDifference(solution,referenceSolution).infinity_norm()
+ << ", max-norm error: " << difference.infinity_norm()
+ << std::endl;
+
+ std::cout << "Level: " << i-1
+ << ", L2 error: " << l2Norm(difference)
+ << std::endl;
+
+ std::cout << "Level: " << i-1
+ << ", H1 error: " << h1Norm(difference)
<< std::endl;
}
--
GitLab