Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • lisa_julia.nebel_at_tu-dresden.de/dune-solvers
  • patrick.jaap_at_tu-dresden.de/dune-solvers
  • burchardt_at_igpm.rwth-aachen.de/dune-solvers
  • agnumpde/dune-solvers
4 results
Select Git revision
Show changes
Showing
with 1370 additions and 757 deletions
SUBDIRS =
iterationstepsdir = $(includedir)/dune/solvers/iterationsteps
iterationsteps_HEADERS = amgstep.hh \
blockgsstep.hh \
blockgsstep.cc \
cgstep.hh \
cgstep.cc \
istlseqilu0step.hh \
iterationstep.hh \
lineariterationstep.hh \
linegsstep.hh \
linegsstep.cc \
minimalpolynomialextrapolationstep.hh \
mmgstep.hh mmgstep.cc \
multigridstep.hh \
multigridstep.cc \
obstacletnnmgstep.hh \
projectedblockgsstep.hh \
projectedblockgsstep.cc \
projectedlinegsstep.cc projectedlinegsstep.hh \
richardsonstep.hh \
truncatedblockgsstep.hh \
truncatedsaddlepointgsstep.hh \
trustregiongsstep.cc \
trustregiongsstep.hh
EXTRA_DIST = CMakeLists.txt
include $(top_srcdir)/am/global-rules
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef ISTL_AMG_STEP_HH #ifndef ISTL_AMG_STEP_HH
#define ISTL_AMG_STEP_HH #define ISTL_AMG_STEP_HH
/** \file /** \file
\brief A wrapper class for the ISTL AMG implementation \brief A wrapper class for the ISTL AMG implementation
*/ */
#include <memory>
#include <dune/solvers/iterationsteps/lineariterationstep.hh> #include <dune/solvers/iterationsteps/lineariterationstep.hh>
#include <dune/istl/paamg/amg.hh> #include <dune/istl/paamg/amg.hh>
namespace Dune {
namespace Solvers {
/** \brief A wrapper class for the ISTL AMG implementation /** \brief A wrapper class for the ISTL AMG implementation
*/ */
template <class MatrixType, class VectorType> template <class MatrixType, class VectorType>
...@@ -21,58 +28,53 @@ class AMGStep ...@@ -21,58 +28,53 @@ class AMGStep
/** \brief Use a sequential SSOR for smoothing */ /** \brief Use a sequential SSOR for smoothing */
typedef Dune::SeqSSOR<MatrixType,VectorType,VectorType> Smoother; typedef Dune::SeqSSOR<MatrixType,VectorType,VectorType> Smoother;
typedef typename Dune::Amg::template SmootherTraits<Smoother>::Arguments SmootherArgs; typedef typename Dune::Amg::template SmootherTraits<Smoother>::Arguments SmootherArgs;
typedef Dune::Amg::AMG<Operator,VectorType,Smoother> AMG; typedef Dune::Amg::AMG<Operator,VectorType,Smoother> AMG;
public: public:
/** \brief Default constructor */ /** \brief Default constructor */
AMGStep () {} AMGStep () = default;
/** \brief Constructor which initializes and sets up an algebraic hierarchy /** \brief Constructor which initializes and sets up an algebraic hierarchy
\param smootherArgs Arguments for the smoother. See the dune-istl documentation for details \param smootherArgs Arguments for the smoother. See the dune-istl documentation for details
\param coarseningCriterion Arguments for the coarsening. See the dune-istl documentation for details \param coarseningCriterion Arguments for the coarsening. See the dune-istl documentation for details
*/ */
AMGStep (const MatrixType* stiffnessMatrix, AMGStep (const MatrixType* stiffnessMatrix,
VectorType& x, VectorType& x,
VectorType& rhs, VectorType& rhs,
const SmootherArgs& smootherArgs, const SmootherArgs& smootherArgs,
const Criterion& coarseningCriterion) const Criterion& coarseningCriterion)
{ {
setProblem(*stiffnessMatrix, x, rhs); setProblem(*stiffnessMatrix, x, rhs);
fop_ = std::auto_ptr<Operator>(new Operator(*stiffnessMatrix));
amg_ = std::auto_ptr<AMG>(new AMG(*fop_, coarseningCriterion, smootherArgs, 1, 1, 1, false));
amg_->pre(*this->x_,*this->rhs_);
} }
/** \brief Initialize the iteration step but don't actually build the matrix hierarchy yet */ /** \brief Initialize the iteration step but don't actually build the matrix hierarchy yet */
AMGStep (const MatrixType* stiffnessMatrix, AMGStep (const MatrixType* stiffnessMatrix,
VectorType& x, VectorType& x,
VectorType& rhs) VectorType& rhs)
{ {
setProblem(*stiffnessMatrix, x, rhs); setProblem(*stiffnessMatrix, x, rhs);
} }
/** \brief Initialize the iteration step but don't actually build the matrix hierarchy yet */ /** \brief Initialize the iteration step but don't actually build the matrix hierarchy yet */
virtual void setProblem(const MatrixType& stiffnessMatrix, void setProblem(const MatrixType& stiffnessMatrix,
VectorType& x, VectorType& x,
VectorType& rhs) const VectorType& rhs) override
{ {
LinearIterationStep<MatrixType, VectorType>::setProblem(stiffnessMatrix, x, rhs); LinearIterationStep<MatrixType, VectorType>::setProblem(stiffnessMatrix, x, rhs);
fop_ = std::auto_ptr<Operator>(new Operator(stiffnessMatrix)); residual_ = rhs;
} }
/** \brief Sets up an algebraic hierarchy /** \brief Sets up an algebraic hierarchy
*/ */
virtual void preprocess(); void preprocess() override;
/** \brief Perform one iteration */ /** \brief Perform one iteration */
virtual void iterate(); void iterate() override;
/** \brief Get the solution */ void apply(VectorType& x, const VectorType& rhs) override;
virtual VectorType getSol();
/** \brief Arguments for the smoother. See the dune-istl documentation for details */ /** \brief Arguments for the smoother. See the dune-istl documentation for details */
SmootherArgs smootherArgs_; SmootherArgs smootherArgs_;
...@@ -83,29 +85,51 @@ public: ...@@ -83,29 +85,51 @@ public:
private: private:
/** \brief A matrix wrapper demanded by istl */ /** \brief A matrix wrapper demanded by istl */
std::auto_ptr<Operator> fop_; std::unique_ptr<Operator> fop_;
/** \brief The dune-istl AMG step */ /** \brief The dune-istl AMG step */
std::auto_ptr<AMG> amg_; std::unique_ptr<AMG> amg_;
VectorType residual_;
void setupCoarseningCriterion()
{
coarseningCriterion_.setNoPreSmoothSteps(1);
coarseningCriterion_.setNoPostSmoothSteps(1);
coarseningCriterion_.setGamma(1);
coarseningCriterion_.setAdditive(false);
}
}; };
template <class MatrixType, class VectorType> template <class MatrixType, class VectorType>
void AMGStep<MatrixType,VectorType>::preprocess() void AMGStep<MatrixType,VectorType>::preprocess()
{ {
amg_ = std::auto_ptr<AMG>(new AMG(*fop_, coarseningCriterion_, smootherArgs_, 1, 1, 1, false)); setupCoarseningCriterion();
amg_->pre(*this->x_,*this->rhs_); fop_ = std::make_unique<Operator>(*this->getMatrix());
amg_ = std::make_unique<AMG>(*fop_, coarseningCriterion_, smootherArgs_);
amg_->pre(*this->x_, residual_);
} }
template <class MatrixType, class VectorType> template <class MatrixType, class VectorType>
void AMGStep<MatrixType,VectorType>::iterate() void AMGStep<MatrixType,VectorType>::iterate()
{ {
amg_->apply(*this->x_,*this->rhs_); amg_->apply(*this->x_, residual_);
} }
template <class MatrixType, class VectorType> template <class MatrixType, class VectorType>
VectorType AMGStep<MatrixType,VectorType>::getSol() void AMGStep<MatrixType,VectorType>::apply(VectorType& x, const VectorType& rhs)
{ {
return *this->x_; x = 0;
this->x_ = &x;
residual_ = rhs;
preprocess();
iterate();
} }
} /* namespace Solvers */
} /* namespace Dune */
// For backward compatibility: will be removed eventually
using Dune::Solvers::AMGStep;
#endif #endif
#include <cassert>
template<class OperatorType, class DiscFuncType, class BitVectorType>
inline
DiscFuncType BlockGSStep<OperatorType, DiscFuncType, BitVectorType>::getSol()
{
return *(this->x_);
}
template<class OperatorType, class DiscFuncType, class BitVectorType>
inline
void BlockGSStep<OperatorType, DiscFuncType, BitVectorType>::
residual(int index, VectorBlock& r) const
{
const OperatorType& mat = *this->mat_;
typedef typename OperatorType::row_type RowType;
const RowType& row = mat[index];
typedef typename RowType::ConstIterator ColumnIterator;
r = (*this->rhs_)[index];
/* The following loop subtracts
* \f[ sum_i = \sum_j A_{ij}w_j \f]
*/
ColumnIterator cIt = row.begin();
ColumnIterator cEndIt = row.end();
for (; cIt!=cEndIt; ++cIt) {
// r_i -= A_ij x_j
cIt->mmv((*this->x_)[cIt.index()], r);
}
}
template<class OperatorType, class DiscFuncType, class BitVectorType>
inline
void BlockGSStep<OperatorType, DiscFuncType, BitVectorType>::iterate()
{
const OperatorType& mat = *this->mat_;
assert(this->ignoreNodes_ != NULL);
for (size_t i=0; i<this->x_->size(); i++) {
size_t const count = (*this->ignoreNodes_)[i].count();
if (count == BlockSize)
continue;
VectorBlock r;
residual(i, r);
// Compute x_i += A_{i,i}^{-1} r[i]
VectorBlock v;
VectorBlock& x = (*this->x_)[i];
if (count == 0) {
// No degree of freedom shall be ignored --> solve linear problem
mat[i][i].solve(v, r);
} else {
// Copy the matrix and adjust rhs and matrix so that dofs given in ignoreNodes[i]
// are not touched
typename OperatorType::block_type matRes;
for (int j = 0; j < BlockSize; ++j) {
if ((*this->ignoreNodes_)[i][j]) {
r[j] = 0;
for (int k = 0; k < BlockSize; ++k)
matRes[j][k] = (k == j);
} else
matRes[j] = mat[i][i][j];
}
matRes.solve(v, r);
}
// Add correction;
x += v;
}
}
#ifndef DUNE_BLOCK_GAUSS_SEIDEL_STEP_HH
#define DUNE_BLOCK_GAUSS_SEIDEL_STEP_HH
#include <dune/common/bitsetvector.hh>
#include "lineariterationstep.hh"
/** \brief A linear Gauß-Seidel iteration step for convex problems which have a block-structure.
*
* \tparam OperatorType The linear operator type
* \tparam DiscFuncType The block vector type of the right hand side and the iterates
* \tparam BitVectorType The type of the bit-vector specifying degrees of freedom that shall be ignored.
*
*/
template<class OperatorType,
class DiscFuncType,
class BitVectorType = Dune::BitSetVector<DiscFuncType::block_type::dimension> >
class BlockGSStep : public LinearIterationStep<OperatorType, DiscFuncType, BitVectorType>
{
typedef typename DiscFuncType::block_type VectorBlock;
enum {BlockSize = VectorBlock::dimension};
public:
//! Default constructor. Doesn't init anything
BlockGSStep() {}
//! Constructor with a linear problem
BlockGSStep(const OperatorType& mat, DiscFuncType& x, const DiscFuncType& rhs)
: LinearIterationStep<OperatorType,DiscFuncType>(mat, x, rhs)
{}
//! Perform one iteration
virtual void iterate();
//! Return the solution
virtual DiscFuncType getSol();
/** \brief Compute the residual of the current iterate of a (block) degree of freedom
*
* \param index Global index of the (block) degree of freedom
* \param r Write residual in this vector
*/
virtual void residual(int index, VectorBlock& r) const;
};
#include "blockgsstep.cc"
#endif
#include <algorithm>
#include "blockgssteps.hh"
namespace Dune {
namespace Solvers {
namespace BlockGS {
std::istream& operator>>(std::istream& lhs, Direction& t) {
std::string s;
lhs >> s;
std::transform(s.begin(), s.end(), s.begin(), ::tolower);
if (s == "forward")
t = Direction::FORWARD;
else if (s == "backward")
t = Direction::BACKWARD;
else if (s == "symmetric")
t = Direction::SYMMETRIC;
else
lhs.setstate(std::ios_base::failbit);
return lhs;
}
}
std::istream& operator>>(std::istream& lhs, BlockGSType& t) {
std::string s;
lhs >> s;
std::transform(s.begin(), s.end(), s.begin(), ::tolower);
if (s == "direct")
t = BlockGSType::Direct;
else if (s == "ldlt")
t = BlockGSType::LDLt;
else if (s == "cg")
t = BlockGSType::CG;
else if (s == "gs")
t = BlockGSType::GS;
else
lhs.setstate(std::ios_base::failbit);
return lhs;
}
}
}
#ifndef DUNE_SOLVERS_ITERATIONSTEPS_BLOCKGSSTEPS_HH
#define DUNE_SOLVERS_ITERATIONSTEPS_BLOCKGSSTEPS_HH
#include <cassert>
#include <functional>
#include <type_traits>
#include <dune/common/parametertree.hh>
#include <dune/matrix-vector/ldlt.hh>
#include <dune/matrix-vector/algorithm.hh>
#include <dune/solvers/common/defaultbitvector.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
namespace Dune {
namespace Solvers {
namespace BlockGS {
enum class Direction { FORWARD, BACKWARD, SYMMETRIC };
/**
* @brief Provides support to parse the @Direction enum.
*/
std::istream& operator>>(std::istream& lhs, Direction& t);
/**
* \brief Iterates over all rows i and updates the i-th component by solving
* the i-th correction problem by means of localSolver.
* \param m Matrix
* \param x Solution vector
* \param b Inhomogeneity
* \param ignore Ignored DOFs
* \param localSolver
* \param direction Iteration direction (e.g. forward, backward, ...)
*/
template <class M, class V, class BitVector, class LocalSolver>
void linearStep(const M& m, V& x, const V& b, const BitVector* ignore,
LocalSolver&& localSolver,
Direction direction = Direction::FORWARD) {
// Note: move-capture requires C++14
auto blockStep = [&, localSolver = std::move(localSolver) ](size_t i) {
const auto& row_i = m[i];
// Compute residual
auto ri = b[i];
using Block = typename M::block_type;
const Block* diag = nullptr;
MatrixVector::sparseRangeFor(row_i, [&](const auto& mij, auto j) {
mij.mmv(x[j], ri);
if (j == i)
diag = &mij;
});
using Ignore = std::bitset<V::block_type::dimension>;
// Update iterate with correction
x[i] += localSolver(diag != nullptr ? *diag : Block(0.0), std::move(ri),
ignore != nullptr ? (*ignore)[i] : Ignore(0));
};
if (direction != Direction::BACKWARD)
for (size_t i = 0; i < x.size(); ++i)
blockStep(i);
if (direction != Direction::FORWARD)
for (size_t i = x.size(); i > 0; --i)
blockStep(i - 1);
}
/// \brief Convenience method for LinearIterationSteps.
template <class M, class V, class BitVector, class LocalSolver>
void linearStep(LinearIterationStep<M, V, BitVector>& lis,
LocalSolver&& localSolver,
Direction direction = Direction::FORWARD) {
BitVector const* ignore = lis.hasIgnore() ? &lis.ignore() : nullptr;
linearStep(*lis.mat_, *lis.x_, *lis.rhs_, ignore,
std::forward<LocalSolver>(localSolver), direction);
}
/**
* \brief All functions in this namespace are assumed to compute the solution
* to a correction problem. Hence they do not require an initial guess
* (and assume it to be zero) and return the (approximate) solution.
*/
namespace LocalSolvers {
namespace LinearSolvers {
/** \brief Solves a block problem using the MBlock-inherent solve method. See
* e.g. FieldMatrix::solve.
*/
template <class MBlock, class VBlock>
VBlock direct(const MBlock& m, const VBlock& b) {
VBlock x;
m.solve(x, b);
return x;
}
//! \brief Solves a block problem using LDL^t decomposition.
template <class MBlock, class VBlock>
VBlock ldlt(MBlock m, const VBlock& b) {
VBlock x;
Dune::MatrixVector::ldlt(m, m, m);
Dune::MatrixVector::ldltSolve(m, m, b, x);
return x;
}
static constexpr size_t defaultCgMaxIter = 100;
static constexpr double defaultCgTol = 1e-15;
//! \brief Solves a block problem using the conjugate gradient method.
template <class MBlock, class VBlock>
VBlock cg(const MBlock& m, VBlock r, size_t maxIter = defaultCgMaxIter,
double tol = defaultCgTol) {
VBlock x(0);
VBlock p = r;
auto q = r;
auto rOldSquared = r * r;
auto tolSquared = tol * tol;
for (size_t j = 0; j < maxIter; ++j) {
m.mv(p, q);
if (rOldSquared < tolSquared)
break;
auto alpha = rOldSquared / (q * p);
x.axpy(alpha, p);
r.axpy(-alpha, q);
auto rSquared = r * r;
auto beta = rSquared / rOldSquared;
p *= beta;
p += r;
rOldSquared = rSquared;
}
return x;
}
static constexpr double defaultGsTol = 0.0;
/**
* \brief Solves a block problem using a single Gauss--Seidel iteration.
* \param tol Lines with diagonal entries below this value are skipped and
* their respective solutions are set to 0.
*/
template <class MBlock, class VBlock>
VBlock gs(const MBlock& m, const VBlock& b, double tol = defaultGsTol) {
VBlock x(0);
for (size_t i = 0; i < m.N(); ++i) {
const auto& mi = m[i];
const auto& mii = mi[i];
if (std::abs(mii) <= tol)
continue;
x[i] = b[i];
MatrixVector::sparseRangeFor(mi, [&](const auto& mij, auto j) {
if (j != i)
x[i] -= mij * x[j];
});
x[i] /= mii;
}
return x;
}
} // end namespace LinearSolvers
namespace LocalSolverFromLinearSolver {
// Note: move-capture, auto-return and auto-lambda-arguments require C++14
/**
* \brief Applies truncation to the block problem according to the local ignore
* nodes. The system is modified in such a way that for nodes that are to be
* ignored, the corresponding rhs is set to 0 while the corresponding matrix row
* is set to be the appropriate euclidean basis vector, essentially enforcing an
* exact solution for this node to be zero. This leads to zero correction
* contributions for the ignored DOFs.
* \param ignore The global ignore nodes where the local ignore nodes can be
* obtained from.
* \param localSolver The block solver used to solve the (modified) local
* problem.
*/
template <class LinearSolver>
auto truncateSymmetrically(LinearSolver&& linearSolver) {
return [linearSolver = std::move(linearSolver) ](
const auto& m, const auto& b, const auto& ignore) {
using Return = std::invoke_result_t<LinearSolver, decltype(m), decltype(b)>;
if (ignore.all())
return Return(0);
if (ignore.none())
return linearSolver(m, b);
auto mTruncated = m;
auto bTruncated = b;
assert(b.size() == m.N());
assert(m.N() == m.M());
size_t blockSize = b.size();
for (size_t j = 0; j < blockSize; ++j) {
if (not ignore[j])
continue;
bTruncated[j] = 0;
for (size_t k = 0; k < blockSize; ++k)
mTruncated[j][k] = mTruncated[k][j] = (k == j);
}
return linearSolver(std::move(mTruncated), std::move(bTruncated));
};
}
} // end namespace LocalSolverFromLinearSolver
namespace LocalSolverRegularizer {
// a suitable value for defaultDiagRegularization could be 1e-10
// nevertheless, the default is 0 to avoid unintentional regularization
static constexpr double defaultDiagRegularizeParameter = 0;
// Note: move-capture, auto-return and auto-lambda-arguments require C++14
/**
* \brief Applies regularization to the diagonal of the block problem, to
* render possibly indefinite problems definite (and hence enable the
* application of more block problem solve techniques).
* \param p Diagonal entries that are non-zero are scaled by (1+p) while
* zero diagonal entries are set to be p.
*/
template <class LocalSolver>
auto diagRegularize(double p, LocalSolver&& localSolver) {
return [ p, localSolver = std::move(localSolver) ](
const auto& m, const auto& b, const auto& ignore) {
if (p == 0)
return localSolver(m, b, ignore);
auto m_regularized = m;
for (size_t j = 0; j < m_regularized.N(); ++j)
if (!ignore[j]) {
if (m_regularized[j][j] == 0)
m_regularized[j][j] = p;
else
m_regularized[j][j] *= 1.0 + p;
}
return localSolver(std::move(m_regularized), b, ignore);
};
}
} // end namespace LocalSolverRegularizer
//! \brief is @LinearSolvers::direct with ignore nodes and regularization.
//! \param r Regularization parameter. Set to 0.0 (default value)
//! to switch off regularization.
struct Direct {
Direct(double r = LocalSolverRegularizer::defaultDiagRegularizeParameter)
: r_(r) {}
template <class M, class V, class I>
V operator()(const M& m, const V& b, const I& ignore) const {
auto directSolver = LinearSolvers::direct<M, V>;
auto directTruncatedSolver =
LocalSolverFromLinearSolver::truncateSymmetrically(directSolver);
auto directTruncatedRegularizedSolver =
LocalSolverRegularizer::diagRegularize(r_, directTruncatedSolver);
return directTruncatedRegularizedSolver(m, b, ignore);
}
private:
double r_;
};
template <class... Args>
auto direct(Args&&... args) {
return Direct(std::forward<Args>(args)...);
}
//! \brief is @LinearSolvers::ldlt with ignore nodes and regularization.
//! \param r Regularization parameter. Set to 0.0 (default value)
//! to switch off regularization.
struct LDLt {
LDLt(double r = LocalSolverRegularizer::defaultDiagRegularizeParameter)
: r_(r) {}
template <class M, class V, class I>
V operator()(const M& m, const V& b, const I& ignore) const {
auto ldltSolver = LinearSolvers::ldlt<M, V>;
auto ldltTruncatedSolver =
LocalSolverFromLinearSolver::truncateSymmetrically(ldltSolver);
auto ldltTruncatedRegularizedSolver =
LocalSolverRegularizer::diagRegularize(r_, ldltTruncatedSolver);
return ldltTruncatedRegularizedSolver(m, b, ignore);
}
private:
double r_;
};
template <class... Args>
auto ldlt(Args&&... args) {
return LDLt(std::forward<Args>(args)...);
}
//! \brief is @LinearSolvers::cg with ignore nodes and regularization.
//! \param r Regularization parameter. Set to 0.0 (default value)
//! to switch off regularization.
struct CG {
CG(size_t maxIter = LinearSolvers::defaultCgMaxIter,
double tol = LinearSolvers::defaultCgTol,
double r = LocalSolverRegularizer::defaultDiagRegularizeParameter)
: maxIter_(maxIter)
, tol_(tol)
, r_(r) {}
template <class M, class V, class I>
V operator()(const M& m, const V& b, const I& ignore) const {
using namespace std::placeholders;
auto cgSolver = std::bind(LinearSolvers::cg<M, V>, _1, _2, maxIter_, tol_);
auto cgTruncatedSolver =
LocalSolverFromLinearSolver::truncateSymmetrically(cgSolver);
auto cgTruncatedRegularizedSolver =
LocalSolverRegularizer::diagRegularize(r_, cgTruncatedSolver);
return cgTruncatedRegularizedSolver(m, b, ignore);
}
private:
size_t maxIter_;
double tol_;
double r_;
};
template <class... Args>
auto cg(Args&&... args) {
return CG(std::forward<Args>(args)...);
}
//! \brief is @LinearSolvers::gs with ignore nodes and regularization.
//! \param r Regularization parameter. Set to 0.0 (default value)
//! to switch off regularization.
struct GS {
GS(double tol = LinearSolvers::defaultGsTol,
double r = LocalSolverRegularizer::defaultDiagRegularizeParameter)
: tol_(tol)
, r_(r) {}
template <class M, class V, class I>
V operator()(const M& m, const V& b, const I& ignore) const {
using namespace std::placeholders;
auto gsSolver = std::bind(LinearSolvers::gs<M, V>, _1, _2, tol_);
auto gsTruncatedSolver =
LocalSolverFromLinearSolver::truncateSymmetrically(gsSolver);
auto gsTruncatedRegularizedSolver =
LocalSolverRegularizer::diagRegularize(r_, gsTruncatedSolver);
return gsTruncatedRegularizedSolver(m, b, ignore);
}
private:
double tol_;
double r_;
};
template <class... Args>
auto gs(Args&&... args) {
return GS(std::forward<Args>(args)...);
}
} // end namespace LocalSolvers
} // end namespace BlockGS
/**
* \brief A Gauss--Seidel-type linear iteration step.
* \param localSolver The solver for the linear block correction problems.
*/
template <class LocalSolver, class Matrix, class Vector,
class BitVector = Dune::Solvers::DefaultBitVector_t<Vector>>
struct BlockGSStep : public LinearIterationStep<Matrix, Vector, BitVector> {
//! \brief Implicitly default-construct local solver
BlockGSStep(BlockGS::Direction direction = BlockGS::Direction::FORWARD)
: localSolver_()
, direction_(direction) {}
//! \brief Use given local solver instance
template <class LS>
BlockGSStep(LS&& localSolver,
BlockGS::Direction direction = BlockGS::Direction::FORWARD)
: localSolver_(std::forward<LS>(localSolver))
, direction_(direction) {}
//! Implicitly construct local solver with forwarded arguments
//! \note Gauss--Seidel direction is hardwired to 'forward' here!
template <class... Args>
BlockGSStep(Args&&... args)
: localSolver_(std::forward<Args>(args)...)
, direction_(BlockGS::Direction::FORWARD) {}
void iterate() {
assert(this->mat_ != nullptr);
assert(this->x_ != nullptr);
assert(this->rhs_ != nullptr);
BlockGS::linearStep(*this, localSolver_, direction_);
}
private:
LocalSolver localSolver_;
BlockGS::Direction direction_;
};
/**
* @brief The Type enum provides representatives for (new) block Gauss-Seidel
* steps that can be constructed via the @BlockGSStepFactory.
*/
enum class BlockGSType { Direct, LDLt, CG, GS };
/**
* @brief Provides support to parse the @BlockGSType enum.
*/
std::istream& operator>>(std::istream& lhs, BlockGSType& t);
/**
* @brief The BlockGSStepFactory struct allows to conveniently construct various
* (new) block Gauss-Seidel steps through a ParameterTree config.
*/
template <class Matrix, class Vector,
class BitVector = DefaultBitVector_t<Vector>>
struct BlockGSStepFactory {
using BlockGSStepPtr =
std::shared_ptr<LinearIterationStep<Matrix, Vector, BitVector>>;
// TODO add ProjectedBlockGSStep (not linear!)
// TODO add CachedLDLtBlockGSStep
//! \brief Create a BlockGSStep with custom local solver.
template <class LocalSolver>
static auto create(
LocalSolver&& localSolver,
BlockGS::Direction direction = BlockGS::Direction::FORWARD) {
return BlockGSStep<LocalSolver, Matrix, Vector, BitVector>(
std::forward<LocalSolver>(localSolver), direction);
}
template <class LocalSolver>
static auto createPtr(
LocalSolver&& localSolver,
BlockGS::Direction direction = BlockGS::Direction::FORWARD) {
return std::make_shared<
BlockGSStep<LocalSolver, Matrix, Vector, BitVector>>(
std::forward<LocalSolver>(localSolver), direction);
}
static BlockGSStepPtr createPtrFromConfig(const Dune::ParameterTree& config) {
auto type = config.get<BlockGSType>("type");
auto dir = config.get<BlockGS::Direction>("direction",
BlockGS::Direction::FORWARD);
auto reg = config.get("regularize_diag",
BlockGS::LocalSolvers::LocalSolverRegularizer::
defaultDiagRegularizeParameter);
switch (type) {
case BlockGSType::Direct:
return createPtr(BlockGS::LocalSolvers::direct(reg), dir);
case BlockGSType::LDLt:
return createPtr(BlockGS::LocalSolvers::ldlt(reg), dir);
case BlockGSType::CG: {
auto maxIter = config.get(
"maxIter", BlockGS::LocalSolvers::LinearSolvers::defaultCgMaxIter);
auto tol =
config.get("tol", BlockGS::LocalSolvers::LinearSolvers::defaultCgTol);
return createPtr(BlockGS::LocalSolvers::cg(maxIter, tol, reg), dir);
}
case BlockGSType::GS: {
auto tol =
config.get("tol", BlockGS::LocalSolvers::LinearSolvers::defaultGsTol);
return createPtr(BlockGS::LocalSolvers::gs(tol, reg), dir);
}
default:
DUNE_THROW(Dune::NotImplemented,
"Factory cannot create this BlockGSType.");
}
}
};
} // end namespace Solvers
} // end namespace Dune
#endif // DUNE_SOLVERS_ITERATIONSTEPS_BLOCKGSSTEPS_HH
template <class MatrixType, class VectorType> // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
void CGStep<MatrixType, VectorType>::check() const // vi: set et ts=8 sw=4 sts=4:
#include <dune/solvers/common/canignore.hh>
#include <dune/matrix-vector/genericvectortools.hh>
namespace Dune {
namespace Solvers {
template <class MatrixType, class VectorType, class Ignore>
void CGStep<MatrixType, VectorType, Ignore>::check() const
{ {
if (preconditioner_)
preconditioner_->check();
} }
template <class MatrixType, class VectorType> template <class MatrixType, class VectorType, class Ignore>
void CGStep<MatrixType, VectorType>::preprocess() void CGStep<MatrixType, VectorType, Ignore>::preprocess()
{ {
// Compute the residual (r starts out as the rhs) // Compute the residual (r starts out as the rhs)
matrix_.mmv(x_,r_); this->mat_->mmv(*x_,r_);
if (this->hasIgnore())
Dune::MatrixVector::Generic::truncate(r_, this->ignore());
if (preconditioner_) { if (preconditioner_) {
preconditioner_->setMatrix(matrix_); using CanIgnore_t = CanIgnore<Ignore>;
CanIgnore_t *asCanIgnore = dynamic_cast<CanIgnore_t*>(preconditioner_.get());
if (asCanIgnore != nullptr and this->hasIgnore())
asCanIgnore->setIgnore(this->ignore());
preconditioner_->setMatrix(*this->mat_);
preconditioner_->apply(p_, r_); preconditioner_->apply(p_, r_);
} else } else
p_ = r_; p_ = r_;
...@@ -20,16 +34,23 @@ void CGStep<MatrixType, VectorType>::preprocess() ...@@ -20,16 +34,23 @@ void CGStep<MatrixType, VectorType>::preprocess()
r_squared_old_ = p_*r_; r_squared_old_ = p_*r_;
} }
template <class MatrixType, class VectorType> template <class MatrixType, class VectorType, class Ignore>
void CGStep<MatrixType, VectorType>::iterate() void CGStep<MatrixType, VectorType, Ignore>::iterate()
{ {
VectorType q(x_); // Avoid divide-by-zero. If r_squared was zero, we're done anyway.
if (r_squared_old_ <= 0)
return;
VectorType q(*x_);
matrix_.mv(p_, q); // q_0 = Ap_0 this->mat_->mv(p_, q); // q_0 = Ap_0
const double alpha = r_squared_old_ / (p_ * q); // alpha_0 = r_0*r_0/p_0*Ap_0 const double alpha = r_squared_old_ / (p_ * q); // alpha_0 = r_0*r_0/p_0*Ap_0
x_.axpy(alpha, p_); // x_1 = x_0 + alpha_0 p_0 x_->axpy(alpha, p_); // x_1 = x_0 + alpha_0 p_0
r_.axpy(-alpha, q); // r_1 = r_0 - alpha_0 Ap_0 r_.axpy(-alpha, q); // r_1 = r_0 - alpha_0 Ap_0
if (this->hasIgnore())
Dune::MatrixVector::Generic::truncate(r_, this->ignore());
if (preconditioner_) if (preconditioner_)
preconditioner_->apply(q, r_); preconditioner_->apply(q, r_);
else else
...@@ -41,3 +62,6 @@ void CGStep<MatrixType, VectorType>::iterate() ...@@ -41,3 +62,6 @@ void CGStep<MatrixType, VectorType>::iterate()
p_ += q; p_ += q;
r_squared_old_ = r_squared; r_squared_old_ = r_squared;
} }
} /* namespace Solvers */
} /* namespace Dune */
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_SOLVERS_ITERATIONSTEPS_CGSTEP_HH #ifndef DUNE_SOLVERS_ITERATIONSTEPS_CGSTEP_HH
#define DUNE_SOLVERS_ITERATIONSTEPS_CGSTEP_HH #define DUNE_SOLVERS_ITERATIONSTEPS_CGSTEP_HH
#include <memory>
#include <dune/solvers/common/preconditioner.hh> #include <dune/solvers/common/preconditioner.hh>
#include <dune/solvers/iterationsteps/iterationstep.hh> #include <dune/solvers/common/defaultbitvector.hh>
#include <dune/solvers/common/wrapownshare.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
namespace Dune { namespace Dune {
namespace Solvers { namespace Solvers {
//! A conjugate gradient solver //! A conjugate gradient solver
template <class MatrixType, class VectorType> template <class MatrixType, class VectorType, class Ignore = DefaultBitVector_t<VectorType>>
class CGStep : public IterationStep<VectorType> class CGStep : public LinearIterationStep<MatrixType,VectorType,Ignore>
{ {
using Base = LinearIterationStep<MatrixType,VectorType,Ignore>;
using Preconditioner = Dune::Solvers::Preconditioner<MatrixType, VectorType, Ignore>;
public: public:
CGStep() CGStep() = default;
{}
CGStep(const MatrixType& matrix, CGStep(const MatrixType& matrix,
VectorType& x, VectorType& x,
const VectorType& rhs) const VectorType& rhs)
: p_(rhs.size()), r_(rhs), x_(x), matrix_(matrix), : Base(matrix,x), p_(rhs.size()), r_(rhs)
preconditioner_(nullptr)
{} {}
template<typename P>
CGStep(const MatrixType& matrix, CGStep(const MatrixType& matrix,
VectorType& x, VectorType& x,
const VectorType& rhs, const VectorType& rhs,
Preconditioner<MatrixType, VectorType>& preconditioner) P&& preconditioner)
: p_(rhs.size()), r_(rhs), x_(x), matrix_(matrix), : Base(matrix,x), p_(rhs.size()), r_(rhs),
preconditioner_(&preconditioner) preconditioner_(wrap_own_share<Preconditioner>(std::forward<P>(preconditioner)))
{} {}
//! Set linear operator, solution and right hand side //! Set linear operator, solution and right hand side
virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) override
{ {
matrix_ = mat; this->setMatrix(mat);
x_ = x; setProblem(x);
r_ = rhs; r_ = rhs;
p_.resize(r_.size());
}
template<typename P>
void setPreconditioner(P&& preconditioner)
{
preconditioner_ = wrap_own_share<Preconditioner>(std::forward<P>(preconditioner));
} }
void check() const; void check() const override;
void preprocess(); void preprocess() override;
void iterate(); void iterate() override;
// Q: do we really want this interface?
VectorType getSol() { return x_; } using Base::setProblem;
private: private:
VectorType p_; // search direction VectorType p_; // search direction
VectorType r_; // residual VectorType r_; // residual
VectorType& x_; using Base::x_;
const MatrixType& matrix_;
double r_squared_old_; double r_squared_old_;
Preconditioner<MatrixType, VectorType>* preconditioner_; std::shared_ptr<Preconditioner> preconditioner_;
}; };
#include "cgstep.cc"
} }
} }
#include "cgstep.cc"
#endif #endif
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef ISTL_SEQILU0_STEP_HH #ifndef ISTL_SEQILU0_STEP_HH
#define ISTL_SEQILU0_STEP_HH #define ISTL_SEQILU0_STEP_HH
/** \file /** \file
\brief A wrapper class for the ISTL SeqILU0 implementation \brief A wrapper class for the ISTL SeqILU(0) implementation
*/ */
#include <memory> #include <memory>
...@@ -10,80 +12,56 @@ ...@@ -10,80 +12,56 @@
#include <dune/solvers/iterationsteps/lineariterationstep.hh> #include <dune/solvers/iterationsteps/lineariterationstep.hh>
#include <dune/istl/preconditioners.hh> #include <dune/istl/preconditioners.hh>
#include <dune/istl/operators.hh>
/** \brief A wrapper class for the ISTL SeqILU0 implementation /** \brief A wrapper class for the ISTL SeqILU(0) implementation
*/ */
template <class MatrixType, class VectorType> template <class MatrixType, class VectorType>
class ISTLSeqILU0Step class ISTLSeqILU0Step
// FIXME: ignoreNodes are not handled
: public LinearIterationStep<MatrixType, VectorType> : public LinearIterationStep<MatrixType, VectorType>
{ {
typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator; typedef Dune::SeqILU<MatrixType,VectorType,VectorType> SeqILU0;
typedef Dune::SeqILU0<MatrixType,VectorType,VectorType> SeqILU0;
public: public:
ISTLSeqILU0Step (double relaxationFactor)
/** \brief Constructor which initializes and sets up an algebraic hierarchy
\param smootherArgs Arguments for the smoother. See the dune-istl documentation for details
\param coarseningCriterion Arguments for the coarsening. See the dune-istl documentation for details
*/
ISTLSeqILU0Step (const MatrixType* stiffnessMatrix,
VectorType& x,
VectorType& rhs,
double relaxationFactor)
: relaxationFactor_(relaxationFactor) : relaxationFactor_(relaxationFactor)
{ {}
setProblem(*stiffnessMatrix, x, rhs);
seqILU0_ = std::auto_ptr<SeqILU0>(new SeqILU0(*stiffnessMatrix, relaxationFactor));
seqILU0_->pre(*this->x_,*this->rhs_);
}
/** \brief Initialize the iteration step but don't actually build the matrix hierarchy yet */
ISTLSeqILU0Step (const MatrixType* stiffnessMatrix,
VectorType& x,
VectorType& rhs)
{
setProblem(*stiffnessMatrix, x, rhs);
}
/** \brief Sets up an algebraic hierarchy /** \brief Sets up an algebraic hierarchy
*/ */
virtual void preprocess(); virtual void preprocess() override {
if (needReconstruction_) {
seqILU0_ = std::make_unique<SeqILU0>(*this->mat_, 0, relaxationFactor_);
needReconstruction_ = false;
}
// Note: as of now, pre() is a dummy
mutable_rhs = *this->rhs_;
seqILU0_->pre(*this->x_,mutable_rhs);
}
/** \brief Perform one iteration */ /** \brief Perform one iteration */
virtual void iterate(); virtual void iterate() override {
mutable_rhs = *this->rhs_;
seqILU0_->apply(*this->x_, mutable_rhs);
}
/** \brief Get the solution */ virtual void setMatrix(const MatrixType& mat) override {
virtual VectorType getSol(); this->mat_ = Dune::stackobject_to_shared_ptr(mat);
needReconstruction_ = true;
}
private: private:
/** \brief The dune-istl sequential ILU0 preconditioner */ /** \brief The dune-istl sequential ILU0 preconditioner */
std::auto_ptr<SeqILU0> seqILU0_; std::unique_ptr<SeqILU0> seqILU0_;
/** \brief Some magic relaxation factor */ VectorType mutable_rhs;
double relaxationFactor_;
};
template <class MatrixType, class VectorType>
void ISTLSeqILU0Step<MatrixType,VectorType>::preprocess()
{
seqILU0_ = std::auto_ptr<SeqILU0>(new SeqILU0(*this->mat_, relaxationFactor_));
seqILU0_->pre(*this->x_,*this->rhs_);
}
template <class MatrixType, class VectorType>
void ISTLSeqILU0Step<MatrixType,VectorType>::iterate()
{
seqILU0_->apply(*this->x_,*this->rhs_);
}
template <class MatrixType, class VectorType> double const relaxationFactor_;
VectorType ISTLSeqILU0Step<MatrixType,VectorType>::getSol()
{
return *this->x_;
}
/** \brief Calling the Dune::SeqILU0 constructor is expensive; we
thus remember if setMatrix has been called. And only
call the constructor if necessary */
bool needReconstruction_;
};
#endif #endif
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef ISTL_SEQSSOR_STEP_HH
#define ISTL_SEQSSOR_STEP_HH
/** \file
\brief A wrapper class for the ISTL SeqSSOR implementation
*/
#include <memory>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
#include <dune/istl/preconditioners.hh>
/** \brief A wrapper class for the ISTL SeqSSOR implementation
*/
template <class MatrixType, class VectorType>
class ISTLSeqSSORStep
// FIXME: ignoreNodes are not handled
: public LinearIterationStep<MatrixType, VectorType>
{
typedef Dune::SeqSSOR<MatrixType,VectorType,VectorType> SeqSSOR;
public:
ISTLSeqSSORStep (int iterations, double relaxationFactor)
: iterations_(iterations), relaxationFactor_(relaxationFactor)
{}
/** \brief Sets up an algebraic hierarchy
*/
virtual void preprocess() override {
seqSSOR_ = std::make_unique<SeqSSOR>(*this->mat_, iterations_,
relaxationFactor_);
// Note: as of now, pre() is a dummy
mutable_rhs = *this->rhs_;
seqSSOR_->pre(*this->x_, mutable_rhs);
}
/** \brief Perform one iteration */
virtual void iterate() override {
mutable_rhs = *this->rhs_;
seqSSOR_->apply(*this->x_, mutable_rhs);
}
private:
/** \brief The dune-istl sequential SSOR preconditioner */
std::unique_ptr<SeqSSOR> seqSSOR_;
VectorType mutable_rhs;
int const iterations_;
double const relaxationFactor_;
};
#endif
#ifndef ITERATIONSTEP_HH // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
#define ITERATIONSTEP_HH // vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_SOLVERS_ITERATIONSTEPS_ITERATIONSTEP_HH
#define DUNE_SOLVERS_ITERATIONSTEPS_ITERATIONSTEP_HH
#include <vector> #include <vector>
#include <string> #include <string>
#include <dune/common/bitsetvector.hh>
#include <dune/solvers/common/canignore.hh> #include <dune/solvers/common/canignore.hh>
#include <dune/solvers/common/defaultbitvector.hh>
#include <dune/solvers/common/numproc.hh> #include <dune/solvers/common/numproc.hh>
namespace Dune {
namespace Solvers {
//! Base class for iteration steps being called by an iterative solver //! Base class for iteration steps being called by an iterative solver
template<class VectorType, class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> > template<class VectorType, class BitVectorType = Solvers::DefaultBitVector_t<VectorType> >
class IterationStep : virtual public NumProc, public CanIgnore<BitVectorType> class IterationStep : virtual public NumProc, public CanIgnore<BitVectorType>
{ {
public: public:
typedef VectorType Vector;
typedef BitVectorType BitVector;
//! Default constructor //! Default constructor
IterationStep() IterationStep() {
{} x_ = nullptr;
}
/** \brief Destructor */ /** \brief Destructor */
virtual ~IterationStep() {} virtual ~IterationStep() {}
...@@ -37,8 +46,23 @@ class IterationStep : virtual public NumProc, public CanIgnore<BitVectorType> ...@@ -37,8 +46,23 @@ class IterationStep : virtual public NumProc, public CanIgnore<BitVectorType>
//! Do the actual iteration //! Do the actual iteration
virtual void iterate() = 0; virtual void iterate() = 0;
//! Access the stored pointer to iterate
virtual const Vector* getIterate() const
{
return x_;
}
//! Access the stored pointer to iterate
virtual Vector* getIterate()
{
return x_;
}
//! Return solution object //! Return solution object
virtual VectorType getSol() = 0; virtual VectorType getSol()
{
return *getIterate();
}
/** \brief Checks whether all relevant member variables are set /** \brief Checks whether all relevant member variables are set
* \exception SolverError if the iteration step is not set up properly * \exception SolverError if the iteration step is not set up properly
...@@ -58,10 +82,14 @@ class IterationStep : virtual public NumProc, public CanIgnore<BitVectorType> ...@@ -58,10 +82,14 @@ class IterationStep : virtual public NumProc, public CanIgnore<BitVectorType>
//! The solution container //! The solution container
VectorType* x_; VectorType* x_;
};
} // namespace Solvers
using CanIgnore<BitVectorType>::ignoreNodes_; } // namespace Dune
}; // For backward compatibility: will be removed eventually
using Dune::Solvers::IterationStep;
#endif #endif
#ifndef LINEAR_ITERATIONSTEP_HH // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
#define LINEAR_ITERATIONSTEP_HH // vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_SOLVERS_ITERATIONSTEPS_LINEARITERATIONSTEP_HH
#define DUNE_SOLVERS_ITERATIONSTEPS_LINEARITERATIONSTEP_HH
#include <memory>
#include <vector> #include <vector>
#include <dune/common/bitsetvector.hh> #include <dune/common/bitsetvector.hh>
#include <dune/common/shared_ptr.hh> #include <dune/common/shared_ptr.hh>
#include <dune/solvers/common/preconditioner.hh> #include <dune/solvers/common/preconditioner.hh>
#include <dune/solvers/common/defaultbitvector.hh>
#include "iterationstep.hh" #include "iterationstep.hh"
//! Base class for iteration steps being called by an iterative linear solver namespace Dune {
template<class MatrixType, class VectorType, class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> >
namespace Solvers {
/** \brief Base class for iteration steps being called by an iterative solver for linear problems
*
* This class is called Linear, because all iteration step implementations that derive from it
* are supposed to solve linear problems. In particular, this means that the LinearIterationStep
* class knows about matrix data types, and it has the setProblem and setMatrix methods that
* allow to specify the matrix of the linear system to solve.
*
* The 'Linear' in the name is not about the type of algorithm. For example, one step of a
* CG solver would qualify as a LinearIterationStep: CG is not a linear algorithm, but it is
* a solver for solving linear problems.
*/
template<class MatrixType, class VectorType, class BitVectorType = Dune::Solvers::DefaultBitVector_t<VectorType> >
class LinearIterationStep : public IterationStep<VectorType, BitVectorType>, class LinearIterationStep : public IterationStep<VectorType, BitVectorType>,
public Dune::Solvers::Preconditioner<MatrixType, public Dune::Solvers::Preconditioner<MatrixType,
VectorType, VectorType,
BitVectorType> BitVectorType>
{ {
using Base = IterationStep<VectorType, BitVectorType>;
public: public:
//! Default constructor //! Default constructor
LinearIterationStep() {} LinearIterationStep() {
rhs_ = nullptr;
}
/** \brief Destructor */ /** \brief Destructor */
virtual ~LinearIterationStep() {} virtual ~LinearIterationStep() {}
//! Constructor being given linear operator, solution and right hand side //! Constructor with given matrix, solution and right hand side
LinearIterationStep(const MatrixType& mat, VectorType& x, const VectorType& rhs) { LinearIterationStep(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
mat_ = Dune::stackobject_to_shared_ptr(mat); setProblem(mat, x, rhs);
this->x_ = &x;
rhs_ = &rhs;
} }
//! Constructor with given matrix and solution vector
LinearIterationStep(const MatrixType& matrix, VectorType& x)
: IterationStep<VectorType,BitVectorType>(x),
mat_(stackobject_to_shared_ptr(matrix))
{}
//! Set linear operator, solution and right hand side //! Set linear operator, solution and right hand side
virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) { virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
this->x_ = &x; setProblem(x);
rhs_ = &rhs; rhs_ = &rhs;
mat_ = Dune::stackobject_to_shared_ptr(mat); setMatrix(mat);
} }
//! Set linear operator //! Set linear operator
virtual void setMatrix(const MatrixType& mat) { virtual void setMatrix(const MatrixType& mat) {
mat_ = Dune::stackobject_to_shared_ptr(mat); mat_ = Dune::stackobject_to_shared_ptr(mat);
} }
//! Do the actual iteration //! Do the actual iteration
virtual void iterate() = 0; virtual void iterate() = 0;
//! Return solution object
virtual VectorType getSol()
{
return *(this->x_);
}
//! Return linear operator //! Return linear operator
virtual const MatrixType* getMatrix() {return mat_.get();} virtual const MatrixType* getMatrix() {return mat_.get();}
/** \brief Checks whether all relevant member variables are set /** \brief Checks whether all relevant member variables are set
* \exception SolverError if the iteration step is not set up properly * \exception SolverError if the iteration step is not set up properly
*/ */
...@@ -75,16 +94,24 @@ public: ...@@ -75,16 +94,24 @@ public:
x = 0; x = 0;
this->x_ = &x; this->x_ = &x;
rhs_ = &r; rhs_ = &r;
this->preprocess();
iterate(); iterate();
x = getSol(); x = this->getSol();
} }
using Base::setProblem;
//! The container for the right hand side //! The container for the right hand side
const VectorType* rhs_; const VectorType* rhs_;
//! The linear operator //! The linear operator
Dune::shared_ptr<const MatrixType> mat_; std::shared_ptr<const MatrixType> mat_;
}; };
} // namespace Solvers
} // namespace Dune
// For backward compatibility: will be removed eventually
using Dune::Solvers::LinearIterationStep;
#endif #endif
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#include <dune/istl/btdmatrix.hh> #include <dune/istl/btdmatrix.hh>
#include <dune/istl/scaledidmatrix.hh> #include <dune/istl/scaledidmatrix.hh>
template<class OperatorType, class DiscFuncType, class BitVectorType>
inline
DiscFuncType LineGSStep<OperatorType, DiscFuncType, BitVectorType>::getSol()
{
return *(this->x_);
}
template<class OperatorType, class DiscFuncType, class BitVectorType> template<class MatrixType, class DiscFuncType, class BitVectorType>
void LineGSStep<OperatorType, DiscFuncType, BitVectorType >::iterate() void LineGSStep<MatrixType, DiscFuncType, BitVectorType >::iterate()
{ {
// input of this method: x^(k) (not the permuted version of x^(k)!) // input of this method: x^(k) (not the permuted version of x^(k)!)
const OperatorType& mat = *this->mat_; const MatrixType& mat = *this->mat_;
int number_of_blocks = blockStructure_.size(); int number_of_blocks = blockStructure_.size();
...@@ -29,13 +23,13 @@ void LineGSStep<OperatorType, DiscFuncType, BitVectorType >::iterate() ...@@ -29,13 +23,13 @@ void LineGSStep<OperatorType, DiscFuncType, BitVectorType >::iterate()
const int current_block_size = blockStructure_[b_num].size(); const int current_block_size = blockStructure_[b_num].size();
//! compute and save the residuals for the curent block: //! compute and save the residuals for the current block:
// determine the (permuted) residuals r[p(i)],..., r[p(i+current_block_size-1)] // determine the (permuted) residuals r[p(i)],..., r[p(i+current_block_size-1)]
// this means that we determine the residuals for the current block // this means that we determine the residuals for the current block
DiscFuncType permuted_r_i(current_block_size); DiscFuncType permuted_r_i(current_block_size);
for (int k=0; k<current_block_size; k++) for (int k=0; k<current_block_size; k++)
{ {
if ( (*this->ignoreNodes_)[blockStructure_[b_num][k]][0]) if ( this->ignore()[blockStructure_[b_num][k]][0])
permuted_r_i[k] = 0.0; permuted_r_i[k] = 0.0;
else else
this->residual( blockStructure_[b_num][k], permuted_r_i[k]); // get r[p(i+k)] this->residual( blockStructure_[b_num][k], permuted_r_i[k]); // get r[p(i+k)]
...@@ -48,22 +42,22 @@ void LineGSStep<OperatorType, DiscFuncType, BitVectorType >::iterate() ...@@ -48,22 +42,22 @@ void LineGSStep<OperatorType, DiscFuncType, BitVectorType >::iterate()
// (Later: x_now[p(i+k)] = x_now[p(i+k)] + v[p(i+k)] ) // (Later: x_now[p(i+k)] = x_now[p(i+k)] + v[p(i+k)] )
// ////////////////////////////////////////////////////////////////////////////////////////////////// // //////////////////////////////////////////////////////////////////////////////////////////////////
// Copy the linear system for the current line/block into a tridiagonal matrix // Copy the linear system for the current line/block into a tridiagonal matrix
// ////////////////////////////////////////////////////////////////////////////////////////////////// // //////////////////////////////////////////////////////////////////////////////////////////////////
Dune::BTDMatrix<typename OperatorType::block_type> tridiagonalMatrix(current_block_size); Dune::BTDMatrix<typename MatrixType::block_type> tridiagonalMatrix(current_block_size);
for (int j=0; j<current_block_size; j++) { for (int j=0; j<current_block_size; j++) {
if ( (*this->ignoreNodes_)[blockStructure_[b_num][j]][0] ) { if ( this->ignore()[blockStructure_[b_num][j]][0] ) {
// left off-diagonal: // left off-diagonal:
if (j>0) if (j>0)
tridiagonalMatrix[j][j-1] = 0; tridiagonalMatrix[j][j-1] = 0;
// diagonal // diagonal
tridiagonalMatrix[j][j] = Dune::ScaledIdentityMatrix<double,BlockSize>(1); tridiagonalMatrix[j][j] = Dune::ScaledIdentityMatrix<double,BlockSize>(1);
// right off-diagonal: // right off-diagonal:
if (j<current_block_size-1) if (j<current_block_size-1)
tridiagonalMatrix[j][j+1] = 0; tridiagonalMatrix[j][j+1] = 0;
...@@ -73,10 +67,10 @@ void LineGSStep<OperatorType, DiscFuncType, BitVectorType >::iterate() ...@@ -73,10 +67,10 @@ void LineGSStep<OperatorType, DiscFuncType, BitVectorType >::iterate()
// left off-diagonal: // left off-diagonal:
if (j>0) if (j>0)
tridiagonalMatrix[j][j-1] = mat[blockStructure_[b_num][j]][blockStructure_[b_num][j-1]]; tridiagonalMatrix[j][j-1] = mat[blockStructure_[b_num][j]][blockStructure_[b_num][j-1]];
// diagonal // diagonal
tridiagonalMatrix[j][j] = mat[blockStructure_[b_num][j]][blockStructure_[b_num][j]]; tridiagonalMatrix[j][j] = mat[blockStructure_[b_num][j]][blockStructure_[b_num][j]];
// right off-diagonal: // right off-diagonal:
if (j<current_block_size-1) if (j<current_block_size-1)
tridiagonalMatrix[j][j+1] = mat[blockStructure_[b_num][j]][blockStructure_[b_num][j+1]]; tridiagonalMatrix[j][j+1] = mat[blockStructure_[b_num][j]][blockStructure_[b_num][j+1]];
......
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_LINE_GAUSS_SEIDEL_STEP_HH #ifndef DUNE_LINE_GAUSS_SEIDEL_STEP_HH
#define DUNE_LINE_GAUSS_SEIDEL_STEP_HH #define DUNE_LINE_GAUSS_SEIDEL_STEP_HH
#include <dune/common/bitsetvector.hh> #include <dune/matrix-vector/axpy.hh>
#include <dune/solvers/iterationsteps/blockgsstep.hh> #include <dune/common/bitsetvector.hh>
template<class OperatorType, template<class MatrixType,
class DiscFuncType, class DiscFuncType,
class BitVectorType = Dune::BitSetVector<DiscFuncType::block_type::dimension> > class BitVectorType = Dune::BitSetVector<DiscFuncType::block_type::dimension> >
class LineGSStep : public BlockGSStep<OperatorType, DiscFuncType, BitVectorType> class LineGSStep : public LinearIterationStep<MatrixType, DiscFuncType, BitVectorType>
{ {
using Base = public LinearIterationStep<MatrixType, DiscFuncType, BitVectorType>;
typedef typename DiscFuncType::block_type VectorBlock; typedef typename DiscFuncType::block_type VectorBlock;
enum {BlockSize = VectorBlock::dimension}; enum {BlockSize = VectorBlock::dimension};
...@@ -21,7 +23,7 @@ template<class OperatorType, ...@@ -21,7 +23,7 @@ template<class OperatorType,
std::vector<std::vector<unsigned int> > blockStructure_; std::vector<std::vector<unsigned int> > blockStructure_;
public: public:
//! Default constructor. Doesn't init anything //! Default constructor. Doesn't init anything
LineGSStep() {} LineGSStep() {}
...@@ -31,15 +33,27 @@ template<class OperatorType, ...@@ -31,15 +33,27 @@ template<class OperatorType,
{} {}
//! Constructor with a linear problem //! Constructor with a linear problem
LineGSStep(const OperatorType& mat, DiscFuncType& x, DiscFuncType& rhs) LineGSStep(const MatrixType& mat, DiscFuncType& x, DiscFuncType& rhs)
: BlockGSStep<OperatorType,DiscFuncType>(mat, x, rhs) : Base(mat, x, rhs)
{} {}
//! Perform one iteration //! Perform one iteration
virtual void iterate(); virtual void iterate();
virtual DiscFuncType getSol(); /** \brief Compute the residual of the current iterate of a (block) degree of freedom
*
* \param index Global index of the (block) degree of freedom
* \param r Write residual in this vector
*/
virtual void residual(int index, VectorBlock& r) const {
r = (*this->rhs_)[index];
const auto& row = (*this->mat_)[index];
for (auto cIt = row.begin(); cIt!=row.end(); ++cIt) {
// r_i -= A_ij x_j
Dune::MatrixVector::subtractProduct(r, *cIt, (*this->x_)[cIt.index()]);
}
}
}; };
......
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_MINIMAL_POLYNOMIAL_EXTRAPOLATION_STEP_HH #ifndef DUNE_MINIMAL_POLYNOMIAL_EXTRAPOLATION_STEP_HH
#define DUNE_MINIMAL_POLYNOMIAL_EXTRAPOLATION_STEP_HH #define DUNE_MINIMAL_POLYNOMIAL_EXTRAPOLATION_STEP_HH
...@@ -12,23 +14,23 @@ ...@@ -12,23 +14,23 @@
#include <dune/solvers/iterationsteps/iterationstep.hh> #include <dune/solvers/iterationsteps/iterationstep.hh>
/** \brief A single step of the Minimal Polynomial Extrapolation method /** \brief A single step of the Minimal Polynomial Extrapolation method
* *
* Minimal Polynomial Extrapolation (MPE) is a general method for sequence acceleration. * Minimal Polynomial Extrapolation (MPE) is a general method for sequence acceleration.
* It takes a converging sequence (actually, it even works for some nonconverging ones) * It takes a converging sequence (actually, it even works for some nonconverging ones)
* and interpolates the iterates such as to obtain a second sequence which converges * and interpolates the iterates such as to obtain a second sequence which converges
* faster. * faster.
* *
* In this implementation, the original sequence is produced by an IterationStep object. * In this implementation, the original sequence is produced by an IterationStep object.
*/ */
template<class VectorType, template<class VectorType,
class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> > class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> >
class MinimalPolynomialExtrapolationStep class MinimalPolynomialExtrapolationStep
: public IterationStep<VectorType, BitVectorType> : public IterationStep<VectorType, BitVectorType>
{ {
public: public:
/** \brief Constructor /** \brief Constructor
* \param iterationStep The iteration step object that creates the sequence being accelerated * \param iterationStep The iteration step object that creates the sequence being accelerated
* \param restart Restart the method after this number of iterations * \param restart Restart the method after this number of iterations
*/ */
...@@ -41,43 +43,33 @@ public: ...@@ -41,43 +43,33 @@ public:
//! Perform one iteration //! Perform one iteration
virtual void iterate(); virtual void iterate();
/** \brief Retrieve the solution */ /** \brief To be called before starting to iterate
virtual VectorType getSol(); \note This calls the preprocess method for the dependent iteration step class, too!
/** \brief To be called before starting to iterate
\note This calls the preprocess method for the dependant iteration step class, too!
*/ */
virtual void preprocess(); virtual void preprocess();
IterationStep<VectorType,BitVectorType>* iterationStep_; IterationStep<VectorType,BitVectorType>* iterationStep_;
/** \brief Restart the method after this number of iterations */ /** \brief Restart the method after this number of iterations */
int restart_; int restart_;
/** \brief The history of iterates of the original sequence */ /** \brief The history of iterates of the original sequence */
std::vector<VectorType> xHistory_; std::vector<VectorType> xHistory_;
/** \brief The history of differences between the iterates of the original sequence */ /** \brief The history of differences between the iterates of the original sequence */
std::vector<VectorType> U_; std::vector<VectorType> U_;
}; };
template<class VectorType, class BitVectorType>
inline
VectorType MinimalPolynomialExtrapolationStep<VectorType, BitVectorType>::getSol()
{
return *(this->x_);
}
template<class VectorType, class BitVectorType> template<class VectorType, class BitVectorType>
inline inline
void MinimalPolynomialExtrapolationStep<VectorType, BitVectorType>::preprocess() void MinimalPolynomialExtrapolationStep<VectorType, BitVectorType>::preprocess()
{ {
iterationStep_->preprocess(); iterationStep_->preprocess();
xHistory_.resize(1); xHistory_.resize(1);
xHistory_[0] = *this->x_; xHistory_[0] = *this->x_;
U_.resize(0); U_.resize(0);
} }
...@@ -87,13 +79,13 @@ void MinimalPolynomialExtrapolationStep<VectorType, BitVectorType>::iterate() ...@@ -87,13 +79,13 @@ void MinimalPolynomialExtrapolationStep<VectorType, BitVectorType>::iterate()
{ {
// append a copy of the last element to the history, for the starting value // append a copy of the last element to the history, for the starting value
xHistory_.push_back(xHistory_.back()); xHistory_.push_back(xHistory_.back());
// iterate once // iterate once
iterationStep_->x_ = &xHistory_.back(); iterationStep_->x_ = &xHistory_.back();
iterationStep_->iterate(); iterationStep_->iterate();
//std::cout << "x[" << xHistory_.size()-1 << "]:\n" << xHistory_.back() << std::endl; //std::cout << "x[" << xHistory_.size()-1 << "]:\n" << xHistory_.back() << std::endl;
// update the history of sequence differences // update the history of sequence differences
VectorType diff = xHistory_.back(); VectorType diff = xHistory_.back();
diff -= xHistory_[xHistory_.size()-2]; diff -= xHistory_[xHistory_.size()-2];
...@@ -101,35 +93,35 @@ void MinimalPolynomialExtrapolationStep<VectorType, BitVectorType>::iterate() ...@@ -101,35 +93,35 @@ void MinimalPolynomialExtrapolationStep<VectorType, BitVectorType>::iterate()
// current number of iterates used for interpolation // current number of iterates used for interpolation
const size_t k= U_.size()-1; const size_t k= U_.size()-1;
/** Compute weights c by solving the system /** Compute weights c by solving the system
* \f[ U^T_{k-1} U_{k-1} c = -U_{k-1} u_k \f] * \f[ U^T_{k-1} U_{k-1} c = -U_{k-1} u_k \f]
*/ */
typedef Dune::Matrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType; typedef Dune::Matrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
typedef Dune::BlockVector<Dune::FieldVector<double,1> > ScalarVectorType; typedef Dune::BlockVector<Dune::FieldVector<double,1> > ScalarVectorType;
// set up the matrix // set up the matrix
ScalarMatrixType UTU(k,k); ScalarMatrixType UTU(k,k);
for (size_t i=0; i<k; i++) for (size_t i=0; i<k; i++)
for (size_t j=0; j<k; j++) for (size_t j=0; j<k; j++)
UTU[i][j] = U_[i] * U_[j]; UTU[i][j] = U_[i] * U_[j];
// set up the right hand side // set up the right hand side
ScalarVectorType rhs(k); ScalarVectorType rhs(k);
for (size_t i=0; i<k; i++) for (size_t i=0; i<k; i++)
rhs[i] = U_[i] * U_[k]; rhs[i] = U_[i] * U_[k];
rhs *= -1; rhs *= -1;
// solve the system // solve the system
ScalarVectorType c(k); ScalarVectorType c(k);
c = 0; c = 0;
// Technicality: turn the matrix into a linear operator // Technicality: turn the matrix into a linear operator
Dune::MatrixAdapter<ScalarMatrixType,ScalarVectorType,ScalarVectorType> op(UTU); Dune::MatrixAdapter<ScalarMatrixType,ScalarVectorType,ScalarVectorType> op(UTU);
// A preconditioner // A preconditioner
Dune::SeqILU0<ScalarMatrixType,ScalarVectorType,ScalarVectorType> ilu0(UTU,1.0); Dune::SeqILU<ScalarMatrixType,ScalarVectorType,ScalarVectorType> ilu0(UTU, 0, 1.0);
// A preconditioned conjugate-gradient solver // A preconditioned conjugate-gradient solver
Dune::CGSolver<ScalarVectorType> cg(op,ilu0,1e-6, // tolerance Dune::CGSolver<ScalarVectorType> cg(op,ilu0,1e-6, // tolerance
...@@ -145,32 +137,32 @@ void MinimalPolynomialExtrapolationStep<VectorType, BitVectorType>::iterate() ...@@ -145,32 +137,32 @@ void MinimalPolynomialExtrapolationStep<VectorType, BitVectorType>::iterate()
// The last coefficient of the c array is '1' // The last coefficient of the c array is '1'
c.resize(c.size()+1, true); c.resize(c.size()+1, true);
c[c.size()-1] = 1; c[c.size()-1] = 1;
//std::cout << "c:\n" << c << std::endl; //std::cout << "c:\n" << c << std::endl;
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////
// Compute the new accelerated iterate // Compute the new accelerated iterate
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////
VectorType& newIterate = *this->x_; VectorType& newIterate = *this->x_;
newIterate = 0; newIterate = 0;
double cSum = std::accumulate(c.begin(), c.end(), double(0)); double cSum = std::accumulate(c.begin(), c.end(), double(0));
for (size_t i=1; i<=k+1; i++) for (size_t i=1; i<=k+1; i++)
newIterate.axpy(c[i-1]/cSum, xHistory_[i]); newIterate.axpy(c[i-1]/cSum, xHistory_[i]);
//std::cout << "y:\n" << newIterate << std::endl; //std::cout << "y:\n" << newIterate << std::endl;
// Acceleration methods for nonlinear problems should be restarted // Acceleration methods for nonlinear problems should be restarted
// every now and then // every now and then
if (k==restart_) { if (k==restart_) {
std::cout << "Restarting MPE..." << std::endl; std::cout << "Restarting MPE..." << std::endl;
U_.resize(0); U_.resize(0);
xHistory_.resize(1); xHistory_.resize(1);
xHistory_[0] = *this->x_; xHistory_[0] = *this->x_;
} }
} }
#endif #endif
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#include <dune/solvers/transferoperators/truncatedmgtransfer.hh> #include <dune/solvers/transferoperators/truncatedmgtransfer.hh>
#include <dune/solvers/iterationsteps/projectedblockgsstep.hh> #include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
...@@ -12,41 +14,71 @@ template <class MatrixType, class VectorType> ...@@ -12,41 +14,71 @@ template <class MatrixType, class VectorType>
void MonotoneMGStep<MatrixType, VectorType>:: void MonotoneMGStep<MatrixType, VectorType>::
preprocess() preprocess()
{ {
// Unset the recompute bitfields, so we compute the full stiffness matrix hierarchy at the beginning
for (size_t i=0; i<this->mgTransfer_.size(); i++) {
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[i])->setRecomputeBitField(nullptr);
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[i])->setCriticalBitField(nullptr);
}
// Call preprocess of the base class // Call preprocess of the base class
MultigridStep<MatrixType,VectorType>::preprocess(); MultigridStep<MatrixType,VectorType>::preprocess();
// Then specify the subset of entries to be reassembled at each iteration // Then setup the obstacle flags which specify the subset of entries to be reassembled at each iteration
recompute_.resize(this->numLevels()); hasObstacleHierarchy_.resize(this->numLevels());
recompute_[recompute_.size()-1] = (*hasObstacle_)[recompute_.size()-1]; hasObstacleHierarchy_.back() = hasObstacle_;
for (int i=this->mgTransfer_.size()-1; i>=0; i--)
this->mgTransfer_[i]->restrictScalarBitField(recompute_[i+1], recompute_[i]); for (int i=this->mgTransfer_.size()-1; i>=0; i--) {
hasObstacleHierarchy_[i] = std::make_shared<Dune::BitSetVector<dim> >();
this->mgTransfer_[i]->restrict(*hasObstacleHierarchy_[i+1], *hasObstacleHierarchy_[i]);
for (size_t j=0; j<hasObstacleHierarchy_[i]->size(); j++)
if ((*this->ignoreNodesHierarchy_[i])[j].any())
(*hasObstacleHierarchy_[i])[j][0]=false;
}
recompute_.resize(this->mgTransfer_.size());
for (size_t i=0; i<this->mgTransfer_.size(); i++) for (size_t i=0; i<this->mgTransfer_.size(); i++)
dynamic_cast<TruncatedMGTransfer<VectorType>*>(this->mgTransfer_[i])->recompute_ = &recompute_[i]; {
recompute_[i].resize(hasObstacleHierarchy_[i]->size());
recompute_[i].unsetAll();
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[i])->setRecomputeBitField(&recompute_[i]);
}
oldCritical_.resize(this->numLevels());
for (size_t i=0; i<this->numLevels(); i++)
{
oldCritical_[i].resize(hasObstacleHierarchy_[i]->size());
oldCritical_[i].unsetAll();
}
// ///////////////////////////////////////////// // /////////////////////////////////////////////
// Set up base solver // Set up base solver
// ///////////////////////////////////////////// // /////////////////////////////////////////////
if (typeid(*this->basesolver_) == typeid(LoopSolver<VectorType>)) {
LoopSolver<VectorType>* loopBaseSolver = dynamic_cast<LoopSolver<VectorType>* > (this->basesolver_); obstacleHierarchy_.resize(this->numLevels());
obstacleHierarchy_.back() = obstacles_;
for (size_t i=0; i<obstacleHierarchy_.size()-1; i++)
obstacleHierarchy_[i] = std::make_shared<ObstacleVectorType>();
if (typeid(*this->basesolver_) == typeid(::LoopSolver<VectorType>)) {
::LoopSolver<VectorType>* loopBaseSolver = dynamic_cast< ::LoopSolver<VectorType>* > (this->basesolver_.get());
typedef ProjectedBlockGSStep<MatrixType, VectorType> SmootherType; typedef ProjectedBlockGSStep<MatrixType, VectorType> SmootherType;
dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->hasObstacle_ = &(*hasObstacle_)[0]; dynamic_cast<SmootherType*>(&loopBaseSolver->getIterationStep())->hasObstacle_ = hasObstacleHierarchy_[0].get();
dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->obstacles_ = &(*obstacles_)[0]; dynamic_cast<SmootherType*>(&loopBaseSolver->getIterationStep())->obstacles_ = obstacleHierarchy_[0].get();
#if HAVE_IPOPT #if HAVE_IPOPT
} else if (typeid(*this->basesolver_) == typeid(QuadraticIPOptSolver<MatrixType,VectorType>)) { } else if (typeid(*this->basesolver_) == typeid(QuadraticIPOptSolver<MatrixType,VectorType>)) {
QuadraticIPOptSolver<MatrixType,VectorType>* ipoptBaseSolver = dynamic_cast<QuadraticIPOptSolver<MatrixType,VectorType>*> (this->basesolver_); QuadraticIPOptSolver<MatrixType,VectorType>* ipoptBaseSolver = dynamic_cast<QuadraticIPOptSolver<MatrixType,VectorType>*> (this->basesolver_.get());
ipoptBaseSolver->obstacles_ = &(*obstacles_)[0]; ipoptBaseSolver->setObstacles(*obstacleHierarchy_[0]);
#endif #endif
} else { } else {
DUNE_THROW(SolverError, "You can't use " << typeid(*this->basesolver_).name() DUNE_THROW(SolverError, "You can't use " << typeid(*this->basesolver_).name()
<< " as a base solver for contact problems!"); << " as a base solver for obstacle problems!");
} }
} }
...@@ -60,12 +92,9 @@ void MonotoneMGStep<MatrixType, VectorType>::nestedIteration() ...@@ -60,12 +92,9 @@ void MonotoneMGStep<MatrixType, VectorType>::nestedIteration()
// ///////////////////////////// // /////////////////////////////
// Restrict obstacles // Restrict obstacles
// ////////////////////////////// // //////////////////////////////
obstacleRestrictor_->restrict((*obstacles_)[i], obstacleRestrictor_->restrict(*obstacleHierarchy_[i], *obstacleHierarchy_[i-1],
(*obstacles_)[i-1], *hasObstacleHierarchy_[i], *hasObstacleHierarchy_[i-1],
(*hasObstacle_)[i], *this->mgTransfer_[i-1], dummy);
(*hasObstacle_)[i-1],
*this->mgTransfer_[i-1],
dummy);
// Restrict right hand side // Restrict right hand side
this->mgTransfer_[i-1]->restrict(this->rhsHierarchy_[i], this->rhsHierarchy_[i-1]); this->mgTransfer_[i-1]->restrict(this->rhsHierarchy_[i], this->rhsHierarchy_[i-1]);
...@@ -76,21 +105,21 @@ void MonotoneMGStep<MatrixType, VectorType>::nestedIteration() ...@@ -76,21 +105,21 @@ void MonotoneMGStep<MatrixType, VectorType>::nestedIteration()
// If we start from an infeasible configuration, the restricted // If we start from an infeasible configuration, the restricted
// obstacles may be inconsistent. We do an ad hoc correction here. // obstacles may be inconsistent. We do an ad hoc correction here.
// The true obstacles on the maxlevel are not touched. // The true obstacles on the maxlevel are not touched.
for (size_t i=0; i<(*obstacles_)[this->level_].size(); i++) { for (size_t i=0; i<obstacleHierarchy_[this->level_]->size(); i++) {
BoxConstraint<field_type,dim>& cO = (*obstacles_)[this->level_][i]; BoxConstraint<field_type,dim>& cO = (*obstacleHierarchy_[this->level_])[i];
for (int j=0; j<dim; j++) for (int j=0; j<dim; j++)
if (cO.lower(j) > cO.upper(j)) if (cO.lower(j) > cO.upper(j))
cO.lower(j) = cO.upper(j) = (cO.lower(j) + cO.upper(j)) / 2; cO.lower(j) = cO.upper(j) = (cO.lower(j) + cO.upper(j)) / 2;
} }
std::cout << "Nested iteration on level " << this->level_ << std::endl; std::cout << "Nested iteration on level " << this->level_ << std::endl;
iterate(); iterate();
iterate(); iterate();
//std::cout << this->xHierarchy_[this->level_] << std::endl; //std::cout << this->xHierarchy_[this->level_] << std::endl;
this->mgTransfer_[this->level_]->prolong(*this->xHierarchy_[this->level_], *this->xHierarchy_[this->level_+1]); this->mgTransfer_[this->level_]->prolong(*this->xHierarchy_[this->level_], *this->xHierarchy_[this->level_+1]);
} }
} }
...@@ -102,10 +131,10 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate() ...@@ -102,10 +131,10 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
int& level = this->level_; int& level = this->level_;
// Define references just for ease of notation // Define references just for ease of notation
std::vector<Dune::shared_ptr<const MatrixType> >& mat = this->matrixHierarchy_; std::vector<std::shared_ptr<const MatrixType> >& mat = this->matrixHierarchy_;
std::vector<Dune::shared_ptr<VectorType> >& x = this->xHierarchy_; std::vector<std::shared_ptr<VectorType> >& x = this->xHierarchy_;
std::vector<VectorType>& rhs = this->rhsHierarchy_; std::vector<VectorType>& rhs = this->rhsHierarchy_;
std::vector<std::vector<BoxConstraint<field_type,dim> > >& obstacles = *this->obstacles_; auto& obstacles = obstacleHierarchy_;
Dune::BitSetVector<dim> critical(x[level]->size(), false); Dune::BitSetVector<dim> critical(x[level]->size(), false);
...@@ -116,67 +145,89 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate() ...@@ -116,67 +145,89 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
// Determine critical nodes (only for debugging) // Determine critical nodes (only for debugging)
const double eps = 1e-12; const double eps = 1e-12;
for (size_t i=0; i<obstacles[level].size(); i++) { for (size_t i=0; i<obstacles[level]->size(); i++) {
for (int j=0; j<dim; j++) { for (int j=0; j<dim; j++) {
if (obstacles[level][i].lower(j) >= (*x[level])[i][j] - eps if ((*obstacles[level])[i].lower(j) >= (*x[level])[i][j] - eps
|| obstacles[level][i].upper(j) <= (*x[level])[i][j] + eps) { || (*obstacles[level])[i].upper(j) <= (*x[level])[i][j] + eps) {
critical[i][j] = true; critical[i][j] = true;
} }
} }
} }
} else { } else {
// Presmoothing // Presmoothing
ProjectedBlockGSStep<MatrixType,VectorType>* presmoother = dynamic_cast<ProjectedBlockGSStep<MatrixType, VectorType>*>(this->presmoother_[level].get()); ProjectedBlockGSStep<MatrixType,VectorType>* presmoother = dynamic_cast<ProjectedBlockGSStep<MatrixType, VectorType>*>(this->presmoother_[level].get());
assert(presmoother); assert(presmoother);
presmoother->setProblem(*(mat[level]), *x[level], rhs[level]); presmoother->setProblem(*(mat[level]), *x[level], rhs[level]);
presmoother->obstacles_ = &(obstacles[level]); presmoother->obstacles_ = obstacles[level].get();
presmoother->hasObstacle_ = &((*hasObstacle_)[level]); presmoother->hasObstacle_ = hasObstacleHierarchy_[level].get();
presmoother->ignoreNodes_ = this->ignoreNodesHierarchy_[level]; presmoother->ignoreNodes_ = this->ignoreNodesHierarchy_[level];
for (int i=0; i<this->nu1_; i++) for (int i=0; i<this->nu1_; i++)
presmoother->iterate(); presmoother->iterate();
// First, a backup of the obstacles for postsmoothing // First, a backup of the obstacles for postsmoothing
std::vector<BoxConstraint<field_type,dim> > obstacleBackup = obstacles[level]; std::vector<BoxConstraint<field_type,dim> > obstacleBackup = *obstacles[level];
// Compute defect obstacles // Compute defect obstacles
for (size_t i=0; i<(*obstacles_)[level].size(); i++) for (size_t i=0; i<obstacles[level]->size(); i++)
obstacles[level][i] -= (*x[level])[i]; (*obstacles[level])[i] -= (*x[level])[i];
// /////////////////////// // ///////////////////////
// Truncation // Truncation
// /////////////////////// // ///////////////////////
// Determine critical nodes // Determine critical nodes
const double eps = 1e-12; const double eps = 1e-12;
for (size_t i=0; i<obstacles[level].size(); i++) { for (size_t i=0; i<obstacles[level]->size(); i++) {
for (int j=0; j<dim; j++) { for (int j=0; j<dim; j++) {
if (obstacles[level][i].lower(j) >= -eps if ((*obstacles[level])[i].lower(j) >= -eps
|| obstacles[level][i].upper(j) <= eps) { || (*obstacles[level])[i].upper(j) <= eps) {
critical[i][j] = true; critical[i][j] = true;
} }
} }
} }
///////////////////////////////////////////////////////////////////////////////
// Compute the part of the coarse grid matrix that needs to be recomputed.
// There are two reasons why a matrix entry may need to be recomputed:
// 1) A corresponding fine grid vertex switched from critical to non-critical or vice versa
// 2) A corresponding fine grid matrix entry got recomputed
///////////////////////////////////////////////////////////////////////////////
Dune::BitSetVector<dim> changed(critical.size());
for (size_t i=0; i<changed.size(); i++)
for (int j=0; j<dim; j++)
changed[i][j] = (critical[i][j] != oldCritical_[level][i][j]);
if (level < (int) this->numLevels()-1 )
for (size_t i=0; i<changed.size(); i++)
for (int j=0; j<dim; j++)
changed[i][j] = (changed[i][j] || recompute_[level][i][j]);
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[level-1])->restrict(changed, recompute_[level-1]);
oldCritical_[level] = critical;
// Set bitfield of nodes that will be truncated
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[level-1])->setCriticalBitField(&critical);
// Restrict stiffness matrix // Restrict stiffness matrix
dynamic_cast<TruncatedMGTransfer<VectorType>*>(this->mgTransfer_[level-1]) this->mgTransfer_[level-1]->galerkinRestrict(*(mat[level]), *(const_cast<MatrixType*>(mat[level-1].get())));
->galerkinRestrict(*(mat[level]), *(const_cast<MatrixType*>(mat[level-1].get())), critical);
// Restrict obstacles // Restrict obstacles
obstacleRestrictor_->restrict(obstacles[level], obstacles[level-1], obstacleRestrictor_->restrict(*obstacles[level], *obstacles[level-1],
(*hasObstacle_)[level], (*hasObstacle_)[level-1], *hasObstacleHierarchy_[level], *hasObstacleHierarchy_[level-1],
*this->mgTransfer_[level-1], *this->mgTransfer_[level-1],
critical); critical);
// ////////////////////////////////////////////////////////////////////// // //////////////////////////////////////////////////////////////////////
// Restriction: fineResiduum = rhs[level] - mat[level] * x[level]; // Restriction: fineResiduum = rhs[level] - mat[level] * x[level];
// ////////////////////////////////////////////////////////////////////// // //////////////////////////////////////////////////////////////////////
...@@ -184,11 +235,11 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate() ...@@ -184,11 +235,11 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
mat[level]->mmv(*x[level], fineResidual); mat[level]->mmv(*x[level], fineResidual);
// restrict residual // restrict residual
dynamic_cast<TruncatedMGTransfer<VectorType>*>(this->mgTransfer_[level-1])->restrict(fineResidual, rhs[level-1], critical); this->mgTransfer_[level-1]->restrict(fineResidual, rhs[level-1]);
// Choose all zeros as the initial correction // Choose all zeros as the initial correction
*x[level-1] = 0; *x[level-1] = 0;
// /////////////////////////////////////// // ///////////////////////////////////////
// Recursively solve the coarser system // Recursively solve the coarser system
// /////////////////////////////////////// // ///////////////////////////////////////
...@@ -196,66 +247,75 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate() ...@@ -196,66 +247,75 @@ void MonotoneMGStep<MatrixType, VectorType>::iterate()
for (int i=0; i<this->mu_; i++) for (int i=0; i<this->mu_; i++)
iterate(); iterate();
level++; level++;
// //////////////////////////////////////// // ////////////////////////////////////////
// Prolong // Prolong
// //////////////////////////////////////// // ////////////////////////////////////////
// add correction to the presmoothed solution // add correction to the presmoothed solution
VectorType tmp; VectorType tmp;
dynamic_cast<TruncatedMGTransfer<VectorType>*>(this->mgTransfer_[level-1])->prolong(*x[level-1], tmp, critical); this->mgTransfer_[level-1]->prolong(*x[level-1], tmp);
*x[level] += tmp; *x[level] += tmp;
// Remove pointer to the temporary critical bitfield
// this avoids a memory problem when the same mmg step is reused
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[level-1])->setCriticalBitField(nullptr);
// restore the true (non-defect) obstacles // restore the true (non-defect) obstacles
obstacles[level] = obstacleBackup; *obstacles[level] = obstacleBackup;
#ifndef NDEBUG #ifndef NDEBUG
// Debug: is the current iterate really admissible? // Debug: is the current iterate really admissible?
for (size_t i=0; i<obstacles[level].size(); i++) for (size_t i=0; i<obstacles[level]->size(); i++)
for (int j=0; j<VectorType::block_type::dimension; j++) for (int j=0; j<VectorType::block_type::dimension; j++)
if ((*x[level])[i][j] < obstacles[level][i].lower(j) if (((*x[level])[i][j] - (*obstacles[level])[i].lower(j)<-1e-14
|| (*x[level])[i][j] > obstacles[level][i].upper(j)) { || (*x[level])[i][j] - (*obstacles[level])[i].upper(j) >1e-14 )
&& (!(*this->ignoreNodesHierarchy_[level])[i][j])) {
std::cout << "Obstacle disregarded!\n"; std::cout << "Obstacle disregarded!\n";
std::cout << (*x[level])[i] << std::endl << obstacles[level][i] << std::endl; std::cout << (*x[level])[i] << std::endl << (*obstacles[level])[i] << std::endl;
std::cout << "level: " << level << " index: " << i << " komponent: " << j << std::endl; std::cout << "level: " << level << " index: " << i << " komponent: " << j << std::endl;
std::cout << "is " << ((critical[i][j]) ? "" : "not") << "critical" << std::endl; std::cout << "is " << ((critical[i][j]) ? "" : "not") << "critical" << std::endl;
} }
#endif #endif
// Postsmoothing // Postsmoothing
ProjectedBlockGSStep<MatrixType,VectorType>* postsmoother = dynamic_cast<ProjectedBlockGSStep<MatrixType, VectorType>*>(this->postsmoother_[level].get()); ProjectedBlockGSStep<MatrixType,VectorType>* postsmoother = dynamic_cast<ProjectedBlockGSStep<MatrixType, VectorType>*>(this->postsmoother_[level].get());
assert(postsmoother); assert(postsmoother);
postsmoother->setProblem(*(mat[level]), *x[level], rhs[level]); postsmoother->setProblem(*(mat[level]), *x[level], rhs[level]);
postsmoother->obstacles_ = &obstacles[level]; postsmoother->obstacles_ = obstacles[level].get();
postsmoother->hasObstacle_ = &((*hasObstacle_)[level]); postsmoother->hasObstacle_ = hasObstacleHierarchy_[level].get();
postsmoother->ignoreNodes_ = this->ignoreNodesHierarchy_[level]; postsmoother->ignoreNodes_ = this->ignoreNodesHierarchy_[level];
for (int i=0; i<this->nu2_; i++) for (int i=0; i<this->nu2_; i++)
postsmoother->iterate(); postsmoother->iterate();
} }
// //////////////////////////////////////////////////////////////////// // ////////////////////////////////////////////////////////////////////
// Track the number of critical nodes found during this iteration // Track the number of critical nodes found during this iteration
// //////////////////////////////////////////////////////////////////// // ////////////////////////////////////////////////////////////////////
if (level==(int) this->numLevels()-1 && this->verbosity_==NumProc::FULL) { if (level==(int) this->numLevels()-1 && this->verbosity_==NumProc::FULL) {
std::cout << critical.count() << " critical nodes found on level " << level; std::cout << critical.count() << " critical nodes found on level " << level;
int changes = 0; int changes = 0;
for (unsigned int i=0; i<oldCritical.size(); i++) for (unsigned int i=0; i<oldCritical.size(); i++)
for (int j=0; j<dim; j++) for (int j=0; j<dim; j++)
if (oldCritical[i][j] !=critical[i][j]) if (oldCritical[i][j] !=critical[i][j])
changes++; changes++;
std::cout << ", and " << changes << " changes." << std::endl; std::cout << ", and " << changes << " changes." << std::endl;
oldCritical = critical; oldCritical = critical;
} }
// Debug: output energy // Debug: output energy
if (level==(int) this->numLevels()-1 && this->verbosity_==NumProc::FULL) if (level==(int) this->numLevels()-1 && this->verbosity_==NumProc::FULL) {
std::cout << "Total energy: " std::streamsize const oldPrecision = std::cout.precision();
<< std::setprecision(10) << computeEnergy(*mat[level], *x[level], rhs[level]) << std::endl; std::cout << "Total energy: "
<< std::setprecision(10)
<< computeEnergy(*mat[level], *x[level], rhs[level]) << std::endl
<< std::setprecision(oldPrecision);
}
} }
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_MONOTONE_MULTIGRID_STEP_HH #ifndef DUNE_MONOTONE_MULTIGRID_STEP_HH
#define DUNE_MONOTONE_MULTIGRID_STEP_HH #define DUNE_MONOTONE_MULTIGRID_STEP_HH
...@@ -6,6 +8,7 @@ ...@@ -6,6 +8,7 @@
#include <dune/solvers/transferoperators/multigridtransfer.hh> #include <dune/solvers/transferoperators/multigridtransfer.hh>
#include <dune/solvers/common/boxconstraint.hh> #include <dune/solvers/common/boxconstraint.hh>
#include <dune/solvers/common/wrapownshare.hh>
#include <dune/solvers/transferoperators/obstaclerestrictor.hh> #include <dune/solvers/transferoperators/obstaclerestrictor.hh>
/** \brief The general monotone multigrid solver /** \brief The general monotone multigrid solver
...@@ -13,73 +16,104 @@ ...@@ -13,73 +16,104 @@
template<class MatrixType, class VectorType> template<class MatrixType, class VectorType>
class MonotoneMGStep : public MultigridStep<MatrixType, VectorType> class MonotoneMGStep : public MultigridStep<MatrixType, VectorType>
{ {
using Base = MultigridStep<MatrixType,VectorType>;
static const int dim = VectorType::block_type::dimension; static const int dim = VectorType::block_type::dimension;
typedef typename VectorType::field_type field_type; typedef typename VectorType::field_type field_type;
typedef std::vector<BoxConstraint<field_type,dim> > ObstacleVectorType;
public: public:
MonotoneMGStep() {} MonotoneMGStep() : hasObstacleHierarchy_(0), obstacleHierarchy_(0)
MonotoneMGStep(int numLevels)
DUNE_DEPRECATED_MSG("The number of levels is no longer set explicitely, but instead inferred from the number of transfer operators. Just erase the number of levels from the argument list in your constructor call.")
{} {}
MonotoneMGStep(const MatrixType& mat,
VectorType& x,
VectorType& rhs, int numLevels)
DUNE_DEPRECATED_MSG("The number of levels is no longer set explicitely, but instead inferred from the number of transfer operators. Just erase the number of levels from the argument list in your constructor call.")
: MultigridStep<MatrixType, VectorType>(mat, x, rhs, numLevels)
{
oldCritical.resize(x.size(), false);
}
MonotoneMGStep(const MatrixType& mat, MonotoneMGStep(const MatrixType& mat,
VectorType& x, VectorType& x,
VectorType& rhs) VectorType& rhs)
: MultigridStep<MatrixType, VectorType>(mat, x, rhs) : Base(mat, x, rhs),
hasObstacleHierarchy_(0), obstacleHierarchy_(0)
{ {
oldCritical.resize(x.size(), false); oldCritical.resize(x.size(), false);
} }
virtual ~MonotoneMGStep() {} virtual ~MonotoneMGStep() {}
virtual void setProblem(const MatrixType& mat,
VectorType& x,
VectorType& rhs,
int numLevels)
DUNE_DEPRECATED_MSG("Use setProblem(const MatrixType& mat,VectorType& x,const VectorType& rhs). The number of levels is no longer set explicitely, but instead inferred from the number of transfer operators.")
{
MultigridStep<MatrixType, VectorType>::setProblem(mat,x,rhs);
oldCritical.resize(x.size(), false);
}
virtual void setProblem(const MatrixType& mat, virtual void setProblem(const MatrixType& mat,
VectorType& x, VectorType& x,
VectorType& rhs) const VectorType& rhs)
{ {
MultigridStep<MatrixType, VectorType>::setProblem(mat,x,rhs); Base::setProblem(mat,x,rhs);
oldCritical.resize(x.size(), false); oldCritical.resize(x.size(), false);
oldCritical.unsetAll();
} }
virtual void iterate(); virtual void iterate();
virtual void preprocess(); virtual void preprocess();
virtual void nestedIteration(); virtual void nestedIteration();
ObstacleRestrictor<VectorType>* obstacleRestrictor_; //! Set the hasObstacle bitfield
[[deprecated("Setting by raw pointer is deprecated. Use l-value or r-value or a shared_ptr")]]
std::vector<Dune::BitSetVector<1> >* hasObstacle_; void setHasObstacles(Dune::BitSetVector<dim>* hasObstacle) {
hasObstacle_ = Dune::stackobject_to_shared_ptr(*hasObstacle);
std::vector<Dune::BitSetVector<1> > recompute_; }
std::vector<std::vector<BoxConstraint<field_type,dim> > >* obstacles_; //! Set the hasObstacle bitfield
template <class BitVector, class Enable = std::enable_if_t<not std::is_pointer<BitVector>::value> >
void setHasObstacles(BitVector&& hasObstacles) {
hasObstacle_ = Dune::Solvers::wrap_own_share<Dune::BitSetVector<dim> >(std::forward<BitVector>(hasObstacles));
}
//! Set the obstacle field
[[deprecated("Setting by raw pointer is deprecated. Use l-value or r-value or a shared_ptr")]]
void setObstacles(ObstacleVectorType* obstacles) {
obstacles_ = Dune::stackobject_to_shared_ptr(*obstacles);
}
//! Set the obstacle field
template <class ObstacleVector, class Enable = std::enable_if_t<not std::is_pointer<ObstacleVector>::value> >
void setObstacles(ObstacleVector&& obstacles) {
obstacles_ = Dune::Solvers::wrap_own_share<ObstacleVectorType>(std::forward<ObstacleVector>(obstacles));
}
//! Get the obstacle field
ObstacleVectorType& getObstacles() { return *obstacles_; }
//! Set the obstacle restrictor
template <class Restrictor>
void setObstacleRestrictor(Restrictor&& restrictor) {
obstacleRestrictor_ = Dune::Solvers::wrap_own_share<ObstacleRestrictor<VectorType> >(std::forward<Restrictor>(restrictor));
}
protected:
//! Bitfield determining which fine grid nodes have an obstacle
std::shared_ptr<Dune::BitSetVector<dim> > hasObstacle_;
//! Vector containing the obstacle values of the fine grid nodes
std::shared_ptr<ObstacleVectorType> obstacles_;
// Needed to track changes in the set of critical bits, and allows // Needed to track changes in the set of critical bits, and allows
// to check which dofs where critical after the last iteration // to check which dofs where critical after the last iteration
Dune::BitSetVector<dim> oldCritical; Dune::BitSetVector<dim> oldCritical;
//! The restrictor used to construct the coarse obstacles
std::shared_ptr<ObstacleRestrictor<VectorType> > obstacleRestrictor_;
//! Bitfield hierarchy containing the coarse obstacle nodes
std::vector<std::shared_ptr<Dune::BitSetVector<dim>> > hasObstacleHierarchy_;
//! For each level specify the matrix lines that need reassembly
std::vector<Dune::BitSetVector<dim> > recompute_;
//! For each level specify dofs that changed status
std::vector<Dune::BitSetVector<dim> > oldCritical_;
//! Hierarchy containing the obstacle values of the coarse obstacles
std::vector<std::shared_ptr<ObstacleVectorType> > obstacleHierarchy_;
private:
using Base::setProblem;
}; };
#include "mmgstep.cc" #include "mmgstep.cc"
......