start implementation interpolation with barycentric coordinates

[[Imported from SVN: r4138]]

src/averagedistanceassembler.hh

+#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

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>

src/targetspacertrsolver.cc

+
+// 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;
+
+    }
+
+}

src/targetspacertrsolver.hh

+#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