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 2a6f0f9b21f1a2075cb9fd25b8538102deee5639..e4f75601a9ebbc537006ad852a022b3a8c3d35e1 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) { @@ -299,8 +295,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_);