diff --git a/src/averagedistanceassembler.hh b/src/averagedistanceassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..592df5b45f3081ac764d968953f1ee5ed4ceb40f
--- /dev/null
+++ b/src/averagedistanceassembler.hh
@@ -0,0 +1,55 @@
+#ifndef AVERAGE_DISTANCE_ASSEMBLER_HH
+#define AVERAGE_DISTANCE_ASSEMBLER_HH
+
+#include "rotation.hh"
+
+template <class TargetSpace>
+class AverageDistanceAssembler
+{};
+
+
+template <>
+class AverageDistanceAssembler<Rotation<3,double> >
+{
+ typedef Rotation<3,double> TargetSpace;
+
+ static const int size = TargetSpace::TangentVector::size;
+
+public:
+
+ AverageDistanceAssembler(const std::vector<TargetSpace> coefficients,
+ const std::vector<double> weights)
+ : coefficients_(coefficients),
+ weights_(weights)
+ {}
+
+ double value(const TargetSpace& x) {
+
+ double result = 0;
+ for (int i=0; i<coefficients_.size(); i++) {
+ double dist = TargetSpace::distance(coefficients_[i], x);
+ result += 0.5*weights_[i]*dist*dist;
+ }
+
+ return result;
+ }
+
+ void assembleGradient(const TargetSpace& x,
+ TargetSpace::TangentVector& gradient)
+ {
+ DUNE_THROW(Dune::NotImplemented, "assembleGradient");
+ }
+
+ void assembleMatrix(const TargetSpace& x,
+ Dune::FieldMatrix<double,size,size>& matrix)
+ {
+ DUNE_THROW(Dune::NotImplemented, "assembleMatrix");
+ }
+
+ const std::vector<TargetSpace> coefficients_;
+
+ const std::vector<double> weights_;
+
+};
+
+#endif
diff --git a/src/localgeodesicfefunction.hh b/src/localgeodesicfefunction.hh
index beb52d7e5a2203709acc7455ae281258359a2530..086af32838ec9ba02217af31665ff7991ce06679 100644
--- a/src/localgeodesicfefunction.hh
+++ b/src/localgeodesicfefunction.hh
@@ -6,6 +6,8 @@
#include <dune/common/fvector.hh>
#include <dune/common/geometrytype.hh>
+#include <dune/src/averagedistanceassembler.hh>
+#include <dune/src/targetspacertrsolver.hh>
/** \brief A geodesic function from the reference element to a manifold
@@ -66,6 +68,30 @@ evaluate(const Dune::FieldVector<ctype, dim>& local)
return result;
#endif
+
+ Dune::FieldVector<ctype, dim+1> barycentricCoordinates;
+
+ barycentricCoordinates[0] = 1;
+ for (int i=0; i<dim; i++) {
+ barycentricCoordinates[0] -= local[i];
+ barycentricCoordinates[i+1] = local[i];
+ }
+
+ AverageDistanceAssembler<TargetSpace> assembler(coefficients_, barycentricCoordinates);
+
+ TargetSpaceRiemannianTRSolver<TargetSpace> solver;
+
+ solver.setup(&assembler,
+ coefficients_[0], // initial iterate
+ 20, // maxTrustRegionSteps
+ 1, // initial trust region radius
+ 20, // inner iterations
+ 1e-8 // inner tolerance
+ );
+
+ solver.solve();
+
+ return solver.getSol();
}
template <int dim, class ctype, class TargetSpace>
diff --git a/src/targetspacertrsolver.cc b/src/targetspacertrsolver.cc
new file mode 100644
index 0000000000000000000000000000000000000000..fdf5b1ffc5e82b223acb942bbcf362a3370db74d
--- /dev/null
+++ b/src/targetspacertrsolver.cc
@@ -0,0 +1,186 @@
+
+// For using a monotone multigrid as the inner solver
+#include <dune-solvers/iterationsteps/trustregiongsstep.hh>
+#include <dune-solvers/solvers/iterativesolver.hh>
+#include "maxnormtrustregion.hh"
+
+template <class TargetSpace>
+void TargetSpaceRiemannianTRSolver<TargetSpace>::
+setup(const AverageDistanceAssembler<TargetSpace>* assembler,
+ const TargetSpace& x,
+ double tolerance,
+ int maxTrustRegionSteps,
+ double initialTrustRegionRadius,
+ int innerIterations,
+ double innerTolerance)
+{
+ x_ = x;
+ this->tolerance_ = tolerance;
+ maxTrustRegionSteps_ = maxTrustRegionSteps;
+ initialTrustRegionRadius_ = initialTrustRegionRadius;
+ innerIterations_ = innerIterations;
+ innerTolerance_ = innerTolerance;
+
+ // ////////////////////////////////
+ // Create a projected gauss-seidel solver
+ // ////////////////////////////////
+
+ // First create a Gauss-seidel base solver
+ TrustRegionGSStep<MatrixType, CorrectionType>* innerSolverStep = new TrustRegionGSStep<MatrixType, CorrectionType>;
+
+ EnergyNorm<MatrixType, CorrectionType>* energyNorm = new EnergyNorm<MatrixType, CorrectionType>(*baseSolverStep);
+
+ innerSolver_ = new ::LoopSolver<CorrectionType>(innerSolverStep,
+ innerIterations,
+ innerTolerance,
+ energyNorm,
+ Solver::QUIET);
+
+ // //////////////////////////////////////////////////////////
+ // Create obstacles
+ // //////////////////////////////////////////////////////////
+
+ innerSolverStep->obstacle_->resize(1);
+ innerSolverStep->obstacle_->setAll();
+
+}
+
+
+template <class TargetSpace>
+void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
+{
+ MaxNormTrustRegion<blocksize> trustRegion(x_.size(), initialTrustRegionRadius_);
+
+ std::vector<std::vector<BoxConstraint<field_type,blocksize> > > trustRegionObstacles((mgStep)
+ ? mgStep->numLevels_
+ : 0);
+
+ // /////////////////////////////////////////////////////
+ // Trust-Region Solver
+ // /////////////////////////////////////////////////////
+ for (int i=0; i<maxTrustRegionSteps_; i++) {
+
+ if (this->verbosity_ == Solver::FULL) {
+ std::cout << "----------------------------------------------------" << std::endl;
+ std::cout << " Trust-Region Step Number: " << i
+ << ", radius: " << trustRegion.radius()
+ << ", energy: " << assembler_->value(x_) << std::endl;
+ std::cout << "----------------------------------------------------" << std::endl;
+ }
+
+ CorrectionType rhs;
+ CorrectionType corr(x_.size());
+ corr = 0;
+
+ assembler_->assembleGradient(x_, rhs);
+ assembler_->assembleMatrix(x_, hesseMatrix);
+
+ //gradientFDCheck(x_, rhs, *rodAssembler_);
+ //hessianFDCheck(x_, *hessianMatrix_, *rodAssembler_);
+
+ // The right hand side is the _negative_ gradient
+ rhs *= -1;
+
+ mgStep->setProblem(*hessianMatrix_, corr, rhs, grid_->maxLevel()+1);
+
+ trustRegionObstacles.back() = trustRegion.obstacles();
+ mgStep->obstacles_ = &trustRegionObstacles;
+
+ innerSolver_->preprocess();
+
+ // /////////////////////////////
+ // Solve !
+ // /////////////////////////////
+
+ innerSolver_->solve();
+
+ corr = mgStep->getSol();
+
+ //std::cout << "Correction: " << std::endl << corr << std::endl;
+
+
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
+
+ if (corr.infinity_norm() < this->tolerance_) {
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
+
+ if (this->verbosity_ != NumProc::QUIET)
+ std::cout << i+1 << " trust-region steps were taken." << std::endl;
+ break;
+ }
+
+ // ////////////////////////////////////////////////////
+ // Check whether trust-region step can be accepted
+ // ////////////////////////////////////////////////////
+
+ TargetSpace newIterate = x_;
+ newIterate = TargetSpace::exp(newIterate, corr);
+
+ /** \todo Don't always recompute oldEnergy */
+ double oldEnergy = assembler_->value(x_);
+ double energy = assembler_->value(newIterate);
+
+ // compute the model decrease
+ // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
+ // Note that rhs = -g
+ CorrectionType tmp(corr.size());
+ tmp = 0;
+ hessianMatrix_->umv(corr, tmp);
+ double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
+
+ if (/* this->verbosity_ == NumProc::FULL */) {
+ std::cout << "Absolute model decrease: " << modelDecrease
+ << ", functional decrease: " << oldEnergy - energy << std::endl;
+ std::cout << "Relative model decrease: " << modelDecrease / energy
+ << ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
+ }
+
+ assert(modelDecrease >= 0);
+
+ if (energy >= oldEnergy) {
+ // if (this->verbosity_ == NumProc::FULL)
+ printf("Richtung ist keine Abstiegsrichtung!\n");
+ }
+
+ if (energy >= oldEnergy &&
+ (std::abs(oldEnergy-energy)/energy < 1e-9 || modelDecrease/energy < 1e-9)) {
+// if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Suspecting rounding problems" << std::endl;
+
+ // if (this->verbosity_ != NumProc::QUIET)
+ std::cout << i+1 << " trust-region steps were taken." << std::endl;
+
+ x_ = newIterate;
+ break;
+ }
+
+ // //////////////////////////////////////////////
+ // Check for acceptance of the step
+ // //////////////////////////////////////////////
+ if ( (oldEnergy-energy) / modelDecrease > 0.9) {
+ // very successful iteration
+
+ x_ = newIterate;
+ trustRegion.scale(2);
+
+ } else if ( (oldEnergy-energy) / modelDecrease > 0.01
+ || std::abs(oldEnergy-energy) < 1e-12) {
+ // successful iteration
+ x_ = newIterate;
+
+ } else {
+ // unsuccessful iteration
+ trustRegion.scale(0.5);
+ // if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Unsuccessful iteration!" << std::endl;
+ }
+
+ // Write current energy
+// if (this->verbosity_ == NumProc::FULL)
+ std::cout << "--- Current energy: " << energy << " ---" << std::endl;
+
+ }
+
+}
diff --git a/src/targetspacertrsolver.hh b/src/targetspacertrsolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7210319f22637914ff5eda2e9b31e755c1abaa70
--- /dev/null
+++ b/src/targetspacertrsolver.hh
@@ -0,0 +1,74 @@
+#ifndef TARGET_SPACE_RIEMANNIAN_TRUST_REGION_SOLVER_HH
+#define TARGET_SPACE_RIEMANNIAN_TRUST_REGION_SOLVER_HH
+
+#include <dune/istl/matrix.hh>
+
+#include <dune-solvers/boxconstraint.hh>
+#include <dune-solvers/solvers/loopsolver.hh>
+
+/** \brief Riemannian trust-region solver for geodesic finite-element problems */
+template <class TargetSpace>
+class TargetSpaceRiemannianTRSolver
+// : public IterativeSolver<std::vector<TargetSpace>,
+// Dune::BitSetVector<TargetSpace::TangentVector::size> >
+{
+ const static int blocksize = TargetSpace::TangentVector::size;
+
+ // Centralize the field type here
+ typedef double field_type;
+
+ // Some types that I need
+ typedef Dune::FieldMatrix<field_type, blocksize, blocksize> MatrixType;
+ typedef Dune::FieldVector<field_type, blocksize> CorrectionType;
+
+public:
+
+ TargetSpaceRiemannianTRSolver()
+ // : IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL)
+ {}
+
+ /** \brief Set up the solver using a monotone multigrid method as the inner solver */
+ void setup(const AverageDistanceAssembler<TargetSpace>* rodAssembler,
+ const TargetSpace& x,
+ double tolerance,
+ int maxTrustRegionSteps,
+ double initialTrustRegionRadius,
+ int innerIterations,
+ double innerTolerance);
+
+ void solve();
+
+ void setInitialSolution(const TargetSpace& x) {
+ x_ = x;
+ }
+
+ TargetSpace getSol() const {return x_;}
+
+protected:
+
+ /** \brief The solution vector */
+ TargetSpace x_;
+
+ /** \brief The initial trust-region radius in the maximum-norm */
+ double initialTrustRegionRadius_;
+
+ /** \brief Maximum number of trust-region steps */
+ int maxTrustRegionSteps_;
+
+ /** \brief Maximum number of multigrid iterations */
+ int innerIterations_;
+
+ /** \brief Error tolerance of the multigrid QP solver */
+ double innerTolerance_;
+
+ /** \brief The assembler for the material law */
+ const AverageDistanceAssembler<TargetSpace>* assembler_;
+
+ /** \brief The solver for the quadratic inner problems */
+ ::LoopSolver<Dune::array<CorrectionType,1> >* innerSolver_;
+
+};
+
+#include "targetspacertrsolver.cc"
+
+#endif