From 567d2adda542b87c371a1614856f4d2347ad44c9 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Fri, 16 May 2014 09:23:13 +0000
Subject: [PATCH] Do not modify the stiffness matrix to account for Dirichlet
 nodes

The Dirichlet nodes are handled directly via the ignoreNodes_ field.
I don't know how the matrix modification code ever got in here.

Removing the code doesn't seem to make any difference.  I find that
a bit surprising: I would expect multigrid and IPOpt to work better
if the matrix is truly symmetric.

[[Imported from SVN: r9742]]
---
 dune/gfe/riemanniantrsolver.cc | 32 --------------------------------
 1 file changed, 32 deletions(-)

diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index 78fa3d23..d3c02028 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -303,38 +303,6 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
 
             rhs *= -1;        // The right hand side is the _negative_ gradient
 
-            ////////////////////////////////////////////////////////////////////////
-            //   Modify matrix and right-hand side to account for Dirichlet values
-            ////////////////////////////////////////////////////////////////////////
-
-            typedef typename MatrixType::row_type::Iterator ColumnIterator;
-
-            for (size_t j=0; j<ignoreNodes_->size(); j++) {
-
-                if (ignoreNodes_->operator[](j).count() > 0) {
-
-                    // make matrix row an identity row
-                    ColumnIterator cIt    = (*hessianMatrix_)[j].begin();
-                    ColumnIterator cEndIt = (*hessianMatrix_)[j].end();
-
-                    for (; cIt!=cEndIt; ++cIt) {
-                        for (int k=0; k<blocksize; k++) {
-                            if (ignoreNodes_->operator[](j)[k]) {
-                                (*cIt)[k] = 0;
-                                if (j==cIt.index())
-                                    (*cIt)[k][k] = 1;
-                            }
-                        }
-                    }
-
-                    // Dirichlet value.  Zero, because we are solving defect problems
-                    for (int k=0; k<blocksize; k++)
-                        if (ignoreNodes_->operator[](j)[k])
-                             rhs[j][k] = 0;
-                }
-
-            }
-
             if (this->verbosity_ == Solver::FULL)
               std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
 
-- 
GitLab