New TargetSpace implementation for points in hyperbolic half-spaces of any dimension

This is just an initial commit.  It is a copy of the corresponding code for
unit vectors, with just a few names exchanged.  The actual functionality
will come in subsequent commits.

[[Imported from SVN: r9076]]

dune/gfe/Makefile.am

@@ -15,6 +15,7 @@ srcinclude_HEADERS = averagedistanceassembler.hh \
                      globalgfetestfunctionbasis.hh \
                      globalgfetestfunction.hh \
                      harmonicenergystiffness.hh \
+                     hyperbolichalfspacepoint.hh \
                      linearalgebra.hh \
                      localgeodesicfefunction.hh \
                      localgeodesicfestiffness.hh \

dune/gfe/hyperbolichalfspacepoint.hh

+#ifndef HYPERBOLIC_HALF_SPACE_POINT_HH
+#define HYPERBOLIC_HALF_SPACE_POINT_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/gfe/tensor3.hh>
+
+/** \brief A point in the hyperbolic half-space H^N
+
+    \tparam N Dimension of the hyperbolic space space
+    \tparam T The type used for individual coordinates
+*/
+template <class T, int N>
+class HyperbolicHalfspacePoint
+{
+    dune_static_assert(N>=2, "A hyperbolic half-space needs to be at least two-dimensional!");
+    
+    /** \brief Computes sin(x) / x without getting unstable for small x */
+    static T sinc(const T& x) {
+        return (x < 1e-4) ? 1 - (x*x/6) : std::sin(x)/x;
+    }
+
+    /** \brief Compute the derivative of arccos^2 without getting unstable for x close to 1 */
+    static T derivativeOfArcCosSquared(const T& x) {
+        const T eps = 1e-4;
+        if (x > 1-eps) {  // regular expression is unstable, use the series expansion instead
+            return -2 + 2*(x-1)/3 - 4/15*(x-1)*(x-1);
+        } else if (x < -1+eps) {  // The function is not differentiable
+            DUNE_THROW(Dune::Exception, "arccos^2 is not differentiable at x==-1!");
+        } else
+            return -2*std::acos(x) / std::sqrt(1-x*x);
+    }
+
+    /** \brief Compute the second derivative of arccos^2 without getting unstable for x close to 1 */
+    static T secondDerivativeOfArcCosSquared(const T& x) {
+        const T eps = 1e-4;
+        if (x > 1-eps) {  // regular expression is unstable, use the series expansion instead
+            return 2.0/3 - 8*(x-1)/15;
+        } else if (x < -1+eps) {  // The function is not differentiable
+            DUNE_THROW(Dune::Exception, "arccos^2 is not differentiable at x==-1!");
+        } else
+            return 2/(1-x*x) - 2*x*std::acos(x) / std::pow(1-x*x,1.5);
+    }
+
+    /** \brief Compute the third derivative of arccos^2 without getting unstable for x close to 1 */
+    static T thirdDerivativeOfArcCosSquared(const T& x) {
+        const T eps = 1e-4;
+        if (x > 1-eps) {  // regular expression is unstable, use the series expansion instead
+            return -8.0/15 + 24*(x-1)/35;
+        } else if (x < -1+eps) {  // The function is not differentiable
+            DUNE_THROW(Dune::Exception, "arccos^2 is not differentiable at x==-1!");
+        } else {
+            T d = 1-x*x;
+            return 6*x/(d*d) - 6*x*x*std::acos(x)/(d*d*std::sqrt(d)) - 2*std::acos(x)/(d*std::sqrt(d));
+        }
+    }
+
+public:
+
+    /** \brief The type used for coordinates */
+    typedef T ctype;
+
+    /** \brief The type used for global coordinates */
+    typedef Dune::FieldVector<T,N> CoordinateType;
+    
+    /** \brief Dimension of the manifold formed by unit vectors */
+    static const int dim = N-1;
+
+    /** \brief Dimension of the Euclidean space the manifold is embedded in */
+    static const int embeddedDim = N;
+
+    typedef Dune::FieldVector<T,N-1> TangentVector;
+
+    typedef Dune::FieldVector<T,N> EmbeddedTangentVector;
+    
+    /** \brief Default constructor */
+    HyperbolicHalfspacePoint()
+    {}
+    
+    /** \brief Constructor from a vector.  The vector gets normalized */
+    HyperbolicHalfspacePoint(const Dune::FieldVector<T,N>& vector)
+        : data_(vector)
+    {
+        data_ /= data_.two_norm();
+    }
+    
+    /** \brief Constructor from an array.  The array gets normalized */
+    HyperbolicHalfspacePoint(const Dune::array<T,N>& vector)
+    {
+        for (int i=0; i<N; i++)
+            data_[i] = vector[i];
+        data_ /= data_.two_norm();
+    }
+
+    UnitVector<T,N>& operator=(const Dune::FieldVector<T,N>& vector)
+    {
+        data_ = vector;
+        data_ /= data_.two_norm();
+        return *this;
+    }
+
+     /** \brief The exponential map */
+    static HyperbolicHalfspacePoint exp(const HyperbolicHalfspacePoint& p, const TangentVector& v) {
+
+        Dune::FieldMatrix<T,N-1,N> frame = p.orthonormalFrame();
+
+        EmbeddedTangentVector ev;
+        frame.mtv(v,ev);
+            
+        return exp(p,ev);
+    }
+
+     /** \brief The exponential map */
+    static HyperbolicHalfspacePoint exp(const HyperbolicHalfspacePoint& p, const EmbeddedTangentVector& v) {
+
+        assert( std::abs(p.data_*v) < 1e-5 );
+
+        const T norm = v.two_norm();
+        HyperbolicHalfspacePoint result = p;
+        result.data_ *= std::cos(norm);
+        result.data_.axpy(sinc(norm), v);
+        return result;
+    }
+
+    /** \brief Length of the great arc connecting the two points */
+     static T distance(const HyperbolicHalfspacePoint& a, const HyperbolicHalfspacePoint& b) {
+
+         // Not nice: we are in a class for unit vectors, but the class is actually
+         // supposed to handle perturbations of unit vectors as well.  Therefore
+         // we normalize here.
+         T x = a.data_ * b.data_/a.data_.two_norm()/b.data_.two_norm();
+         
+         // paranoia:  if the argument is just eps larger than 1 acos returns NaN
+         x = std::min(x,1.0);
+         
+         return std::acos(x);
+    }
+
+    /** \brief Compute the gradient of the squared distance function keeping the first argument fixed
+
+    Unlike the distance itself the squared distance is differentiable at zero
+     */
+    static EmbeddedTangentVector derivativeOfDistanceSquaredWRTSecondArgument(const HyperbolicHalfspacePoint& a, const HyperbolicHalfspacePoint& b) {
+        T x = a.data_ * b.data_;
+
+        EmbeddedTangentVector result = a.data_;
+
+        result *= derivativeOfArcCosSquared(x);
+
+        // Project gradient onto the tangent plane at b in order to obtain the surface gradient
+        result = b.projectOntoTangentSpace(result);
+
+        // Gradient must be a tangent vector at b, in other words, orthogonal to it
+        assert( std::abs(b.data_ * result) < 1e-5);
+
+        return result;
+    }
+
+    /** \brief Compute the Hessian of the squared distance function keeping the first argument fixed
+
+    Unlike the distance itself the squared distance is differentiable at zero
+     */
+    static Dune::FieldMatrix<T,N,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const HyperbolicHalfspacePoint& p, const HyperbolicHalfspacePoint& q) {
+
+        T sp = p.data_ * q.data_;
+        
+        Dune::FieldVector<T,N> pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
+
+        Dune::FieldMatrix<T,N,N> A;
+        for (int i=0; i<N; i++)
+            for (int j=0; j<N; j++)
+                A[i][j] = pProjected[i]*pProjected[j];
+        
+        A *= secondDerivativeOfArcCosSquared(sp);
+
+        // Compute matrix B (see notes)
+        Dune::FieldMatrix<T,N,N> Pq;
+        for (int i=0; i<N; i++)
+            for (int j=0; j<N; j++)
+                Pq[i][j] = (i==j) - q.data_[i]*q.data_[j];
+
+        // Bring it all together
+        A.axpy(-1*derivativeOfArcCosSquared(sp)*sp, Pq);
+
+        return A;
+    }
+
+    /** \brief Compute the mixed second derivate \partial d^2 / \partial da db
+
+    Unlike the distance itself the squared distance is differentiable at zero
+     */
+    static Dune::FieldMatrix<T,N,N> secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(const HyperbolicHalfspacePoint& a, const HyperbolicHalfspacePoint& b) {
+
+        T sp = a.data_ * b.data_;
+
+        // Compute vector A (see notes)
+        Dune::FieldMatrix<T,1,N> row;
+        row[0] = b.projectOntoTangentSpace(a.globalCoordinates());
+
+        Dune::FieldVector<T,N> tmp = a.projectOntoTangentSpace(b.globalCoordinates());
+        Dune::FieldMatrix<T,N,1> column;
+        for (int i=0; i<N; i++)  // turn row vector into column vector
+            column[i] = tmp[i];
+
+        Dune::FieldMatrix<T,N,N> A;
+        // A = row * column
+        Dune::FMatrixHelp::multMatrix(column,row,A);
+        A *= secondDerivativeOfArcCosSquared(sp);
+
+        // Compute matrix B (see notes)
+        Dune::FieldMatrix<T,N,N> Pp, Pq;
+        for (int i=0; i<N; i++)
+            for (int j=0; j<N; j++) {
+                Pp[i][j] = (i==j) - a.data_[i]*a.data_[j];
+                Pq[i][j] = (i==j) - b.data_[i]*b.data_[j];
+            }
+            
+        Dune::FieldMatrix<T,N,N> B;
+        Dune::FMatrixHelp::multMatrix(Pp,Pq,B);
+
+        // Bring it all together
+        A.axpy(derivativeOfArcCosSquared(sp), B);
+
+        return A;
+    }
+    
+    
+    /** \brief Compute the third derivative \partial d^3 / \partial dq^3
+
+    Unlike the distance itself the squared distance is differentiable at zero
+     */
+    static Tensor3<T,N,N,N> thirdDerivativeOfDistanceSquaredWRTSecondArgument(const HyperbolicHalfspacePoint& p, const HyperbolicHalfspacePoint& q) {
+
+        Tensor3<T,N,N,N> result;
+
+        T sp = p.data_ * q.data_;
+        
+        // The projection matrix onto the tangent space at p and q
+        Dune::FieldMatrix<T,N,N> Pq;
+        for (int i=0; i<N; i++)
+            for (int j=0; j<N; j++)
+                Pq[i][j] = (i==j) - q.globalCoordinates()[i]*q.globalCoordinates()[j];
+            
+        Dune::FieldVector<T,N> pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
+
+        for (int i=0; i<N; i++)
+            for (int j=0; j<N; j++)
+                for (int k=0; k<N; k++) {
+
+                    result[i][j][k] = thirdDerivativeOfArcCosSquared(sp) * pProjected[i] * pProjected[j] * pProjected[k]
+                                    - secondDerivativeOfArcCosSquared(sp) * ((i==j)*sp + p.globalCoordinates()[i]*q.globalCoordinates()[j])*pProjected[k]
+                                    - secondDerivativeOfArcCosSquared(sp) * ((i==k)*sp + p.globalCoordinates()[i]*q.globalCoordinates()[k])*pProjected[j]
+                                    - secondDerivativeOfArcCosSquared(sp) * pProjected[i] * Pq[j][k] * sp
+                                    + derivativeOfArcCosSquared(sp) * ((i==j)*q.globalCoordinates()[k] + (i==k)*q.globalCoordinates()[j]) * sp
+                                    - derivativeOfArcCosSquared(sp) * p.globalCoordinates()[i] * Pq[j][k];
+                }
+                
+        result = Pq * result;
+                
+        return result;
+    }    
+        
+    /** \brief Compute the mixed third derivative \partial d^3 / \partial da db^2
+
+    Unlike the distance itself the squared distance is differentiable at zero
+     */
+    static Tensor3<T,N,N,N> thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(const HyperbolicHalfspacePoint& p, const HyperbolicHalfspacePoint& q) {
+
+        Tensor3<T,N,N,N> result;
+
+        T sp = p.data_ * q.data_;
+        
+        // The projection matrix onto the tangent space at p and q
+        Dune::FieldMatrix<T,N,N> Pp, Pq;
+        for (int i=0; i<N; i++)
+            for (int j=0; j<N; j++) {
+                Pp[i][j] = (i==j) - p.globalCoordinates()[i]*p.globalCoordinates()[j];
+                Pq[i][j] = (i==j) - q.globalCoordinates()[i]*q.globalCoordinates()[j];
+            }
+            
+        Dune::FieldVector<T,N> pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
+        Dune::FieldVector<T,N> qProjected = p.projectOntoTangentSpace(q.globalCoordinates());
+        
+        Tensor3<T,N,N,N> derivativeOfPqOTimesPq;
+        for (int i=0; i<N; i++)
+            for (int j=0; j<N; j++)
+                for (int k=0; k<N; k++) {
+                    derivativeOfPqOTimesPq[i][j][k] = 0;
+                    for (int l=0; l<N; l++)
+                        derivativeOfPqOTimesPq[i][j][k] += Pp[i][l] * (Pq[j][l]*pProjected[k] + pProjected[j]*Pq[k][l]);
+                }
+                
+        result = thirdDerivativeOfArcCosSquared(sp)         * Tensor3<T,N,N,N>::product(qProjected,pProjected,pProjected)
+                 + secondDerivativeOfArcCosSquared(sp)      * derivativeOfPqOTimesPq
+                 - secondDerivativeOfArcCosSquared(sp) * sp * Tensor3<T,N,N,N>::product(qProjected,Pq)
+                 - derivativeOfArcCosSquared(sp)            * Tensor3<T,N,N,N>::product(qProjected,Pq);
+               
+        return result;
+    }
+    
+    
+    /** \brief Project tangent vector of R^n onto the tangent space */
+    EmbeddedTangentVector projectOntoTangentSpace(const EmbeddedTangentVector& v) const {
+        EmbeddedTangentVector result = v;
+        result.axpy(-1*(data_*result), data_);
+        return result;
+    }
+
+    /** \brief The global coordinates, if you really want them */
+    const CoordinateType& globalCoordinates() const {
+        return data_;
+    }
+
+    /** \brief Compute an orthonormal basis of the tangent space of S^n.
+
+    This basis is of course not globally continuous.
+    */
+    Dune::FieldMatrix<T,N-1,N> orthonormalFrame() const {
+
+        Dune::FieldMatrix<T,N-1,N> result;
+
+        // Coordinates of the stereographic projection
+        Dune::FieldVector<T,N-1> X;
+        
+        if (data_[N-1] <= 0) {
+            
+            // Stereographic projection from the north pole onto R^{N-1}
+            for (size_t i=0; i<N-1; i++)
+                X[i] = data_[i] / (1-data_[N-1]);
+            
+        } else {
+            
+            // Stereographic projection from the south pole onto R^{N-1}
+            for (size_t i=0; i<N-1; i++)
+                X[i] = data_[i] / (1+data_[N-1]);
+            
+        }
+            
+        T RSquared = X.two_norm2();
+            
+        for (size_t i=0; i<N-1; i++)
+            for (size_t j=0; j<N-1; j++)
+                // Note: the matrix is the transpose of the one in the paper
+                result[j][i] = 2*(i==j)*(1+RSquared) - 4*X[i]*X[j];
+                
+        for (size_t j=0; j<N-1; j++)
+            result[j][N-1] = 4*X[j];
+            
+        // Upper hemisphere: adapt formulas so it is the stereographic projection from the south pole
+        if (data_[N-1] > 0) 
+            for (size_t j=0; j<N-1; j++)
+                result[j][N-1] *= -1;
+            
+        // normalize the rows to make the orthogonal basis orthonormal
+        for (size_t i=0; i<N-1; i++)
+            result[i] /= result[i].two_norm();
+        
+        return result;
+    }
+
+    /** \brief Write unit vector object to output stream */
+    friend std::ostream& operator<< (std::ostream& s, const HyperbolicHalfspacePoint& unitVector)
+    {
+        return s << unitVector.data_;
+    }
+
+
+private:
+
+    Dune::FieldVector<T,N> data_;
+};
+
+#endif