From 8b80c38706b48fbe70a2a8a40df0bfd2e07e1398 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Tue, 12 May 2009 12:39:20 +0000
Subject: [PATCH] use mass matrix norm to define the trust region
[[Imported from SVN: r4176]]
---
src/riemanniantrsolver.cc | 17 ++++++++++++++---
1 file changed, 14 insertions(+), 3 deletions(-)
diff --git a/src/riemanniantrsolver.cc b/src/riemanniantrsolver.cc
index 32389f07..683e8dd7 100644
--- a/src/riemanniantrsolver.cc
+++ b/src/riemanniantrsolver.cc
@@ -3,8 +3,9 @@
#include <dune/istl/io.hh>
#include <dune/disc/functions/p1function.hh>
-#include <dune/disc/miscoperators/laplace.hh>
#include <dune/disc/operators/p1operator.hh>
+#include <dune/disc/miscoperators/laplace.hh>
+#include <dune/disc/miscoperators/massmatrix.hh>
// For using a monotone multigrid as the inner solver
#include <dune-solvers/iterationsteps/trustregiongsstep.hh>
@@ -180,19 +181,28 @@ setupTCG(const GridType& grid,
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
+ // This is used to measure convergence of the inner solver
// //////////////////////////////////////////////////////////////////////////////////////
Dune::LeafP1Function<GridType,double> u(grid),f(grid);
Dune::LaplaceLocalStiffness<typename GridType::LeafGridView,double> laplaceStiffness;
Dune::LeafP1OperatorAssembler<GridType,double,1>* A = new Dune::LeafP1OperatorAssembler<GridType,double,1>(grid);
A->assemble(laplaceStiffness,u,f);
- typedef typename Dune::LeafP1OperatorAssembler<GridType,double,1>::RepresentationType LaplaceMatrixType;
-
if (h1SemiNorm_)
delete h1SemiNorm_;
h1SemiNorm_ = new H1SemiNorm<CorrectionType>(**A);
+ // //////////////////////////////////////////////////////////////////////////////////////
+ // Assemble a mass matrix to create a norm that's equivalent to the L2-norm
+ // This is used to to define the trust region
+ // //////////////////////////////////////////////////////////////////////////////////////
+
+ Dune::LeafP1Function<GridType,double,blocksize> uvv(grid),fvv(grid);
+ Dune::MassMatrixLocalStiffness<typename GridType::LeafGridView,double,blocksize> massMatrixStiffness;
+ Dune::LeafP1OperatorAssembler<GridType,double,blocksize>* B = new Dune::LeafP1OperatorAssembler<GridType,double,blocksize>(grid);
+ B->assemble(massMatrixStiffness,uvv,fvv);
+
// ////////////////////////////////////////////////////
// Create a truncated conjugate gradient solver
// ////////////////////////////////////////////////////
@@ -200,6 +210,7 @@ setupTCG(const GridType& grid,
innerSolver_ = new TruncatedCGSolver<MatrixType,CorrectionType>(innerIterations_,
innerTolerance_,
h1SemiNorm_,
+ &**B, // the norm of the trust region
Solver::FULL);
// Write all intermediate solutions, if requested
--
GitLab