diff --git a/dune/solvers/common/boxconstraint.hh b/dune/solvers/common/boxconstraint.hh
index 0e812a5bd6418974e7869bf8ddb0ccbd4427ae12..ebd8ce427f8d3b14f78fb43976d911f212b73621 100644
--- a/dune/solvers/common/boxconstraint.hh
+++ b/dune/solvers/common/boxconstraint.hh
@@ -8,63 +8,67 @@
#include <dune/common/fvector.hh>
+#include "interval.hh"
+
template <class T, int dim>
class BoxConstraint {
public:
+ using Interval = typename Dune::Solvers::Interval<T>;
/** \brief Default constructor, makes the obstacle as large as possible */
BoxConstraint() {
- for (int i=0; i<dim; i++) {
- val[2*i] = -std::numeric_limits<T>::max();
- val[2*i+1] = std::numeric_limits<T>::max();
- }
+ clear();
}
/** \brief Set obstacle values to +/- std::numeric_limits<T>::max() */
void clear() {
- for (int i=0; i<dim; i++) {
- val[2*i] = -std::numeric_limits<T>::max();
- val[2*i+1] = std::numeric_limits<T>::max();
- }
+ for (int i=0; i<dim; i++)
+ val[i] = { -std::numeric_limits<T>::max(),
+ std::numeric_limits<T>::max() };
}
//! Subtract vector from box
BoxConstraint<T,dim>& operator-= (const Dune::FieldVector<T, dim>& v)
{
for (int i=0; i<dim; i++) {
- val[2*i] -= v[i];
- val[2*i+1] -= v[i];
+ val[i][0] -= v[i];
+ val[i][1] -= v[i];
}
return *this;
}
+ Interval const &operator[](size_t i) const {
+ return val[i];
+ }
+ Interval &operator[](size_t i) {
+ return val[i];
+ }
+
//! Access to the lower obstacles
- T& lower(size_t i) {return val[2*i];}
+ T& lower(size_t i) {return val[i][0];}
//! Const access to the lower obstacles
- const T& lower(size_t i) const {return val[2*i];}
+ const T& lower(size_t i) const {return val[i][0];}
//! Access to the upper obstacles
- T& upper(size_t i) {return val[2*i+1];}
+ T& upper(size_t i) {return val[i][1];}
//! Const access to the upper obstacles
- const T& upper(size_t i) const {return val[2*i+1];}
+ const T& upper(size_t i) const {return val[i][1];}
//! Send box constraint to an output stream
friend
std::ostream& operator<< (std::ostream& s, const BoxConstraint<T,dim>& v)
{
for (int i=0; i<dim; i++)
- s << "Direction: " << i <<", val[0] = " << v.val[2*i]
- << ", val[1] = " << v.val[2*i+1] << std::endl;
+ s << "Direction: " << i <<", val[0] = " << v[i][0]
+ << ", val[1] = " << v[i][1] << std::endl;
return s;
}
protected:
-
- std::array<T, 2*dim> val;
-
+ std::array<Interval, dim> val;
};
#endif
diff --git a/dune/solvers/common/interval.hh b/dune/solvers/common/interval.hh
index 6530dac10d92641cfaf66d4e655f447695736fd4..8c68f19e9a0186e4de66a71ee809768537417178 100644
--- a/dune/solvers/common/interval.hh
+++ b/dune/solvers/common/interval.hh
@@ -24,8 +24,9 @@ public:
/** \brief Construct from an initializer list */
Interval(std::initializer_list<field_type> const &input)
- : data_(input)
- {}
+ {
+ std::copy(input.begin(), input.end(), data_.begin());
+ }
/** \brief Array-like access
*/
diff --git a/dune/solvers/iterationsteps/obstacletnnmgstep.hh b/dune/solvers/iterationsteps/obstacletnnmgstep.hh
index 80824e16f69e77382b80c281cd8c0eb8f063444e..ebe8eac2c8cb061cb34ae10e2f61ed2ad52755c5 100644
--- a/dune/solvers/iterationsteps/obstacletnnmgstep.hh
+++ b/dune/solvers/iterationsteps/obstacletnnmgstep.hh
@@ -265,12 +265,8 @@ class ObstacleTNNMGStep
// project correction into local defect constraints
for (int j = 0; j < blockSize; ++j)
- {
- if (coarseCorrection_[i][j] < defectConstraint.lower(j))
- coarseCorrection_[i][j] = defectConstraint.lower(j);
- if (coarseCorrection_[i][j] > defectConstraint.upper(j))
- coarseCorrection_[i][j] = defectConstraint.upper(j);
- }
+ coarseCorrection_[i][j]
+ = defectConstraint[j].projectIn(coarseCorrection_[i][j]);
for (int j = 0; j < blockSize; ++j)
{
@@ -302,8 +298,7 @@ class ObstacleTNNMGStep
damping = (residual_*coarseCorrection_) / coarseCorrectionNorm2;
else
damping = 0;
- damping = std::max(damping, directionalConstraints.lower(0));
- damping = std::min(damping, directionalConstraints.upper(0));
+ damping = directionalConstraints[0].projectIn(damping);
// apply correction
x.axpy(damping, coarseCorrection_);
diff --git a/dune/solvers/iterationsteps/projectedlinegsstep.cc b/dune/solvers/iterationsteps/projectedlinegsstep.cc
index 542754ce6b963867eff0c12b45f9fa7b921da1c9..6b806b2ba728ba0400af74fe1086fa1e04b27dae 100755
--- a/dune/solvers/iterationsteps/projectedlinegsstep.cc
+++ b/dune/solvers/iterationsteps/projectedlinegsstep.cc
@@ -136,7 +136,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 914d80402bc3f573366de4495d958741057809a6..ef95d3cf0075a5e325f704e37143a2451f27194c 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 2feb30593812e6a80449af8c02e3d44f9819dae3..3f65815b8199387026fd6bf878e8f3e8f73123fb 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 9438b23ea2ad07aae79bfbdbbece0cd393817345..04aa5cb16cd0d5d3d243956cd97a829cbf0a5b15 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 5cb6525fb9e11b7e55a3f95ce11138e489e9e7fd..3b534924ff53ee8485871fa9358ffcad887e19a8 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 373c0da183af343c050e166814e52a7ff581800f..6321ca035355eec91a46d3e2710b2e057fc2ab37 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 dcc7cdea6f9cd15a67d90ab1d9ae12981ae1d46a..df1fe3468f5a9a27c84c06b95a3965650dd78303 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 e5f19898ed9fa6cf242c9336c6a8ca80db4bff7d..7024362aed532429621e4a9a3e3ec52a4a9b7394 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);