Newer
Older
#ifndef RIEMANNIAN_TRUST_REGION_SOLVER_HH
#define RIEMANNIAN_TRUST_REGION_SOLVER_HH
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/solvers/common/boxconstraint.hh>
#include <dune/solvers/norms/h1seminorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>

Oliver Sander
committed
#include "geodesicfeassembler.hh"
#include <dune/grid/utility/globalindexset.hh>
#include <dune/gfe/parallel/globalmapper.hh>
#include <dune/gfe/parallel/globalp1mapper.hh>

Oliver Sander
committed
#include <dune/gfe/parallel/globalp2mapper.hh>
#include <dune/gfe/parallel/p2mapper.hh>
/** \brief Assign GlobalMapper and LocalMapper types to a dune-fufem FunctionSpaceBasis */
template <typename Basis>
/** \brief Specialization for LagrangeBasis<1> */

Oliver Sander
committed
template <typename GridView>
struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,1> >
typedef Dune::GlobalP1Mapper<Dune::Functions::LagrangeBasis<GridView,1>> GlobalMapper;
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> LocalMapper;
static LocalMapper createLocalMapper(const GridView& gridView)
{
return LocalMapper(gridView, Dune::mcmgVertexLayout());
}

Oliver Sander
committed
template <typename GridView>
struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,2> >
typedef Dune::GlobalP2Mapper<Dune::Functions::LagrangeBasis<GridView,2>> GlobalMapper;
typedef P2BasisMapper<GridView> LocalMapper;
static LocalMapper createLocalMapper(const GridView& gridView)
{
return LocalMapper(gridView);
}
/** \brief Specialization for LagrangeBasis<3> */

Oliver Sander
committed
template <typename GridView>
struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,3> >
{
// Error: we don't currently have a global P3 mapper
};

Oliver Sander
committed
/** \brief Riemannian trust-region solver for geodesic finite-element problems */

Lisa Julia Nebel
committed
template <class Basis, class TargetSpace, class Assembler = GeodesicFEAssembler<Basis,TargetSpace>>
: public IterativeSolver<std::vector<TargetSpace>,
Dune::BitSetVector<TargetSpace::TangentVector::dimension> >
typedef typename Basis::GridView::Grid GridType;
const static int blocksize = TargetSpace::TangentVector::dimension;

Oliver Sander
committed
const static int gridDim = GridType::dimension;
// Centralize the field type here
typedef double field_type;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize, blocksize> > MatrixType;
typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > CorrectionType;

Oliver Sander
committed
typedef std::vector<TargetSpace> SolutionType;
#if HAVE_MPI
typedef typename MapperFactory<Basis>::GlobalMapper GlobalMapper;
typedef typename MapperFactory<Basis>::LocalMapper LocalMapper;
#endif
/** \brief Records information about the last run of the RiemannianTrustRegionSolver
*
* This is used primarily for unit testing.
*/
struct Statistics
{
std::size_t finalIteration;
field_type finalEnergy;
};
RiemannianTrustRegionSolver()

Oliver Sander
committed
: IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
hessianMatrix_(nullptr), h1SemiNorm_(NULL)
{
std::fill(scaling_.begin(), scaling_.end(), 1.0);
}

Oliver Sander
committed
/** \brief Set up the solver using a monotone multigrid method as the inner solver */

Lisa Julia Nebel
committed
const Assembler* assembler,
const Dune::BitSetVector<blocksize>& dirichletNodes,
int maxTrustRegionSteps,
double initialTrustRegionRadius,
int multigridIterations,
double mgTolerance,
int nu1,
int nu2,
int baseIterations,
double baseTolerance,
bool instrumented);
void setScaling(const Dune::FieldVector<double,blocksize>& scaling)
{
void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
{
ignoreNodes_ = &ignoreNodes;
std::shared_ptr<LoopSolver<CorrectionType> > loopSolver = std::dynamic_pointer_cast<LoopSolver<CorrectionType> >(innerSolver_);
assert(loopSolver);
loopSolver->iterationStep_->ignoreNodes_ = ignoreNodes_;
}
[[deprecated]]
void setInitialSolution(const SolutionType& x) {
x_ = x;
}
void setInitialIterate(const SolutionType& x) {
x_ = x;
}
SolutionType getSol() const {return x_;}
const Statistics& getStatistics() const {return statistics_;}
#if HAVE_MPI
std::unique_ptr<GlobalMapper> globalMapper_;
#endif

Oliver Sander
committed
/** \brief The grid */
const GridType* grid_;
/** \brief The solution vector */
SolutionType x_;
/** \brief The initial trust-region radius in the maximum-norm */
double initialTrustRegionRadius_;
/** \brief Trust-region norm scaling */
Dune::FieldVector<double,blocksize> scaling_;
/** \brief Maximum number of trust-region steps */
std::size_t maxTrustRegionSteps_;
/** \brief Maximum number of multigrid iterations */

Oliver Sander
committed
int innerIterations_;
/** \brief Error tolerance of the multigrid QP solver */

Oliver Sander
committed
double innerTolerance_;
std::unique_ptr<MatrixType> hessianMatrix_;
/** \brief The assembler for the material law */

Lisa Julia Nebel
committed
const Assembler* assembler_;

Oliver Sander
committed
/** \brief The solver for the quadratic inner problems */
std::shared_ptr<Solver> innerSolver_;
/** \brief Contains 'true' everywhere -- the trust-region is bounded */
Dune::BitSetVector<blocksize> hasObstacle_;
/** \brief The Dirichlet nodes */
const Dune::BitSetVector<blocksize>* ignoreNodes_;
/** \brief The norm used to measure multigrid convergence */
std::shared_ptr<H1SemiNorm<CorrectionType> > h1SemiNorm_;
/** \brief An L2-norm, really. The H1SemiNorm class is badly named */
std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_;
/** \brief If set to true we log convergence speed and other stuff */
bool instrumented_;
/** \brief Store information about solver runs for unit testing */
Statistics statistics_;

Oliver Sander
committed
#include "riemanniantrsolver.cc"