diff --git a/dune/gfe/realtuple.hh b/dune/gfe/realtuple.hh index c1f52db2e40acabdc6345ae94e88ed3beddf0616..e4bc4df8561ab9c5fff83bd529b29c7ffe97f61d 100644 --- a/dune/gfe/realtuple.hh +++ b/dune/gfe/realtuple.hh @@ -106,7 +106,14 @@ public: Simply the Euclidean distance */ static T distance(const RealTuple& a, const RealTuple& b) { - return (a.data_ - b.data_).two_norm(); + return log(a.data_,b.data_).two_norm(); + } + + /** \brief The logarithmic map + * Simply the difference vector for RealTuple + * */ + static auto log(const RealTuple& a, const RealTuple& b) { + return a.data_ - b.data_; } #if ADOLC_ADOUBLE_H @@ -215,6 +222,11 @@ public: return data_; } + Dune::FieldVector<T,N>& globalCoordinates() { + return data_; + } + + /** \brief Compute an orthonormal basis of the tangent space of R^n. In general this frame field, may of course not be continuous, but for RealTuples it is. diff --git a/dune/gfe/rotation.hh b/dune/gfe/rotation.hh index 96ea13db09cfbd34112410fdb530091c8ede82f0..809138eeab775d3daadebd64c07c4f6826248065 100644 --- a/dune/gfe/rotation.hh +++ b/dune/gfe/rotation.hh @@ -604,7 +604,7 @@ public: /** \brief Compute the vector in T_aSO(3) that is mapped by the exponential map to the geodesic from a to b */ - static SkewMatrix<T,3> difference(const Rotation<T,3>& a, const Rotation<T,3>& b) { + static SkewMatrix<T,3> log(const Rotation<T,3>& a, const Rotation<T,3>& b) { Quaternion<T> diff = a; diff.invert(); @@ -622,7 +622,7 @@ public: // TODO: ADOL-C does not like this part of the code, // because arccos is not differentiable at -1 and 1. - // (Even though the overall 'difference' function is differentiable.) + // (Even though the overall 'log' function is differentiable.) using std::acos; T dist = 2*acos( diff[3] ); @@ -635,7 +635,7 @@ public: T invSinc = 1/sincHalf(dist); - // Compute difference on T_a SO(3) + // Compute log on T_a SO(3) v[0] = diff[0] * invSinc; v[1] = diff[1] * invSinc; v[2] = diff[2] * invSinc; @@ -907,8 +907,8 @@ public: /** \brief Interpolate between two rotations */ static Rotation<T,3> interpolate(const Rotation<T,3>& a, const Rotation<T,3>& b, T omega) { - // Compute difference on T_a SO(3) - SkewMatrix<T,3> v = difference(a,b); + // Compute log on T_a SO(3) + SkewMatrix<T,3> v = log(a,b); v *= omega; @@ -922,8 +922,8 @@ public: T omega) { Quaternion<T> result(0); - // Compute difference on T_a SO(3) - SkewMatrix<T,3> xi = difference(a,b); + // Compute log on T_a SO(3) + SkewMatrix<T,3> xi = log(a,b); SkewMatrix<T,3> v = xi; v *= omega; diff --git a/dune/gfe/skewmatrix.hh b/dune/gfe/skewmatrix.hh index 083cd19152d2ca588af070715cf2bec0042a3ef4..76614d7cfefdcb15bac52e3de919378f50f0f457 100644 --- a/dune/gfe/skewmatrix.hh +++ b/dune/gfe/skewmatrix.hh @@ -48,7 +48,34 @@ public: { return data_; } - + + typedef typename Dune::FieldVector<T,3>::Iterator Iterator; + + //! begin iterator + Iterator begin () + { + return Iterator(data_,0); + } + + //! end iterator + Iterator end () + { + return Iterator(data_,3); + } + + typedef typename Dune::FieldVector<T,3>::ConstIterator ConstIterator; + //! begin iterator + ConstIterator begin () const + { + return ConstIterator(data_,0); + } + + //! end iterator + ConstIterator end () const + { + return ConstIterator(data_,3); + } + /** \brief Embedd the skey-symmetric matrix in R^3x3 */ Dune::FieldMatrix<T,3,3> toMatrix() const {