From cd1bed328a8f1f5e731f7e0c39df97a5da7b4d89 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Tue, 7 Jun 2016 13:30:46 +0200
Subject: [PATCH] Box constraints: Use new constructor, projectIn()

---
 dune/solvers/iterationsteps/projectedlinegsstep.cc |  2 +-
 dune/solvers/iterationsteps/trustregiongsstep.cc   |  7 +------
 dune/solvers/solvers/trustregionsolver.cc          |  6 ++----
 dune/solvers/solvers/trustregionsolver.hh          |  6 ++----
 dune/solvers/test/mmgtest.cc                       |  7 +------
 dune/solvers/test/obstacletnnmgtest.cc             |  7 +------
 dune/solvers/test/projectedgradienttest.cc         | 10 ++++------
 dune/solvers/test/quadraticipoptsolvertest.cc      |  7 +------
 8 files changed, 13 insertions(+), 39 deletions(-)

diff --git a/dune/solvers/iterationsteps/projectedlinegsstep.cc b/dune/solvers/iterationsteps/projectedlinegsstep.cc
index b852c1af..2714f34b 100755
--- a/dune/solvers/iterationsteps/projectedlinegsstep.cc
+++ b/dune/solvers/iterationsteps/projectedlinegsstep.cc
@@ -145,7 +145,7 @@ solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
 
         // Line search
         field_type alpha = (residual*linearCorrection) / correctionNormSquared;
-        alpha = std::max(std::min(alpha, lineSearchObs.upper(0)), lineSearchObs.lower(0));
+        alpha = lineSearchObs[0].projectIn(alpha);
 
         if (isnan(alpha))
             alpha = 1;
diff --git a/dune/solvers/iterationsteps/trustregiongsstep.cc b/dune/solvers/iterationsteps/trustregiongsstep.cc
index 914d8040..ef95d3cf 100644
--- a/dune/solvers/iterationsteps/trustregiongsstep.cc
+++ b/dune/solvers/iterationsteps/trustregiongsstep.cc
@@ -41,12 +41,7 @@ void TrustRegionGSStep<MatrixType, VectorType>::iterate()
                 // 1d problem is convex
                 x[i][j] = r / diag;
 
-                if ( x[i][j] < obstacles[i].lower(j)) {
-                    x[i][j] = obstacles[i].lower(j);
-                }
-                else if ( x[i][j] > obstacles[i].upper(j)) {
-                    x[i][j] = obstacles[i].upper(j);
-                }
+                x[i][j] = obstacles[i][j].projectIn(x[i][j]);
 
             } else  {
 
diff --git a/dune/solvers/solvers/trustregionsolver.cc b/dune/solvers/solvers/trustregionsolver.cc
index 2feb3059..3f65815b 100644
--- a/dune/solvers/solvers/trustregionsolver.cc
+++ b/dune/solvers/solvers/trustregionsolver.cc
@@ -23,10 +23,8 @@ void TrustRegionSolver<ProblemType,VectorType,MatrixType>::solve()
     // setup the trust-region in the maximums norm
     std::vector<BoxConstraint<field_type,blocksize> > trustRegionObstacles(x_.size());
     for (size_t i=0; i<trustRegionObstacles.size(); i++)
-        for (int j=0; j<blocksize; j++) {
-            trustRegionObstacles[i].lower(j) = -initialRadius_;
-            trustRegionObstacles[i].upper(j) = initialRadius_;
-        }
+        for (int j=0; j<blocksize; j++)
+            trustRegionObstacles[i][j] = { -initialRadius_, initialRadius_ };
 
     mgStep->obstacles_ = &trustRegionObstacles;
 
diff --git a/dune/solvers/solvers/trustregionsolver.hh b/dune/solvers/solvers/trustregionsolver.hh
index 9438b23e..04aa5cb1 100644
--- a/dune/solvers/solvers/trustregionsolver.hh
+++ b/dune/solvers/solvers/trustregionsolver.hh
@@ -108,10 +108,8 @@ protected:
             initialRadius_ *= scaling;
 
         for (size_t i=0; i<obs.size(); i++)
-            for (int j=0; j<blocksize; j++) {
-                obs[i].lower(j) = -initialRadius_;
-                obs[i].upper(j) = initialRadius_;
-        }
+            for (int j=0; j<blocksize; j++)
+                obs[i][j] = { -initialRadius_, initialRadius_ };
     }
 
     /** \brief The problem type. */
diff --git a/dune/solvers/test/mmgtest.cc b/dune/solvers/test/mmgtest.cc
index 5cb6525f..3b534924 100644
--- a/dune/solvers/test/mmgtest.cc
+++ b/dune/solvers/test/mmgtest.cc
@@ -66,13 +66,8 @@ void solveObstacleProblemByMMGSolver(const GridType& grid, const MatrixType& mat
 
     std::vector<BoxConstraint<typename VectorType::field_type,blockSize> > boxConstraints(rhs.size());
     for (size_t i = 0; i < boxConstraints.size(); ++i)
-    {
         for (int j = 0; j < blockSize; ++j)
-        {
-            boxConstraints[i].lower(j) = -1;
-            boxConstraints[i].upper(j) = +1;
-        }
-    }
+            boxConstraints[i][j] = { -1, +1 };
     mmgStep.obstacles_         = &boxConstraints;
 
     // we need a vector of pointers to the transfer operator base class
diff --git a/dune/solvers/test/obstacletnnmgtest.cc b/dune/solvers/test/obstacletnnmgtest.cc
index 373c0da1..6321ca03 100644
--- a/dune/solvers/test/obstacletnnmgtest.cc
+++ b/dune/solvers/test/obstacletnnmgtest.cc
@@ -66,13 +66,8 @@ void solveObstacleProblemByTNNMG(const GridType& grid, const MatrixType& mat, Ve
     // create double obstacle constraints
     BoxConstraintVector boxConstraints(rhs.size());
     for (size_t i = 0; i < boxConstraints.size(); ++i)
-    {
         for (int j = 0; j < blockSize; ++j)
-        {
-            boxConstraints[i].lower(j) = -1;
-            boxConstraints[i].upper(j) = +1;
-        }
-    }
+            boxConstraints[i][j] = { -1, +1 };
 
     // we want to control errors with respect to the energy norm induced by the matrix
     Norm norm(mat);
diff --git a/dune/solvers/test/projectedgradienttest.cc b/dune/solvers/test/projectedgradienttest.cc
index dcc7cdea..df1fe346 100644
--- a/dune/solvers/test/projectedgradienttest.cc
+++ b/dune/solvers/test/projectedgradienttest.cc
@@ -30,13 +30,11 @@ bool solveObstacleProblemByProjectedGradientSolver(const MatrixType& mat, Vector
     typedef std::vector<BoxConstraint<double, blockSize> > BoxConstraintVector;
     BoxConstraintVector boxConstraints(rhs.size());
     for (size_t i = 0; i < boxConstraints.size(); ++i)
-    {
         for (size_t j = 0; j < blockSize; ++j)
-        {
-            boxConstraints[i].lower(j) = (ignore[i][j]==true) ? 0 : -1;
-            boxConstraints[i].upper(j) = (ignore[i][j]==true) ? 0 : 1;
-        }
-    }
+          boxConstraints[i][j] = {
+            (ignore[i][j]==true) ? 0.0 : -1.0,
+            (ignore[i][j]==true) ? 0.0 : +1.0
+          };
 
     ProjectedGradientStep<MatrixType, VectorType> projGradStep(mat, x, rhs);
     projGradStep.obstacles_ = &boxConstraints;
diff --git a/dune/solvers/test/quadraticipoptsolvertest.cc b/dune/solvers/test/quadraticipoptsolvertest.cc
index e5f19898..7024362a 100644
--- a/dune/solvers/test/quadraticipoptsolvertest.cc
+++ b/dune/solvers/test/quadraticipoptsolvertest.cc
@@ -35,13 +35,8 @@ void solveObstacleProblemByQuadraticIPOptSolver(const GridType& grid, const Matr
     // create double obstacle constraints
     std::vector<BoxConstraint<typename VectorType::field_type,blockSize> > boxConstraints(rhs.size());
     for (size_t i = 0; i < boxConstraints.size(); ++i)
-    {
         for (int j = 0; j < blockSize; ++j)
-        {
-            boxConstraints[i].lower(j) = -1;
-            boxConstraints[i].upper(j) = +1;
-        }
-    }
+            boxConstraints[i][j] = { -1, +1 };
 
     // create solver
     Solver solver(mat,x,rhs, NumProc::REDUCED);
-- 
GitLab