diff --git a/dune/gfe/realtuple.hh b/dune/gfe/realtuple.hh
index f572f96e145d9a48c7fcd24ebaffa4455806c197..99fd8dbbf0737ab2990f1d91b7ba955165cebf51 100644
--- a/dune/gfe/realtuple.hh
+++ b/dune/gfe/realtuple.hh
@@ -11,7 +11,7 @@
Currently this class only exists for testing purposes.
*/
-template <int N, class T=double>
+template <class T, int N>
class RealTuple
{
public:
@@ -39,7 +39,7 @@ public:
}
/** \brief Copy constructor */
- RealTuple(const RealTuple<N>& other)
+ RealTuple(const RealTuple<T,N>& other)
: data_(other.data_)
{}
diff --git a/dune/gfe/rigidbodymotion.hh b/dune/gfe/rigidbodymotion.hh
index 2d1194f44b9f711b52a03891e1ea3430bf232f47..8945d17861146c33499da0d67520a64ae8b4fb55 100644
--- a/dune/gfe/rigidbodymotion.hh
+++ b/dune/gfe/rigidbodymotion.hh
@@ -7,16 +7,16 @@
#include "rotation.hh"
/** \brief A rigid-body motion in R^N, i.e., a member of SE(N) */
-template <int N, class T=double>
+template <class T, int N>
struct RigidBodyMotion
{
public:
/** \brief Dimension of manifold */
- static const int dim = N + Rotation<N,T>::dim;
+ static const int dim = N + Rotation<T,N>::dim;
/** \brief Dimension of the embedding space */
- static const int embeddedDim = N + Rotation<N,T>::embeddedDim;
+ static const int embeddedDim = N + Rotation<T,N>::embeddedDim;
/** \brief Type of an infinitesimal rigid body motion */
typedef Dune::FieldVector<T, dim> TangentVector;
@@ -36,7 +36,7 @@ public:
/** \brief Constructor from a translation and a rotation */
RigidBodyMotion(const Dune::FieldVector<ctype, N>& translation,
- const Rotation<N,ctype>& rotation)
+ const Rotation<ctype,N>& rotation)
: r(translation), q(rotation)
{}
@@ -59,9 +59,9 @@ public:
and then the compiler complains about having two methods with the same signature.
*/
template <class TVector>
- static RigidBodyMotion<N,ctype> exp(const RigidBodyMotion<N,ctype>& p, const TVector& v) {
+ static RigidBodyMotion<ctype,N> exp(const RigidBodyMotion<ctype,N>& p, const TVector& v) {
- RigidBodyMotion<N,ctype> result;
+ RigidBodyMotion<ctype,N> result;
// Add translational correction
for (int i=0; i<N; i++)
@@ -69,29 +69,29 @@ public:
// Add rotational correction
typedef typename Dune::SelectType<Dune::is_same<TVector,TangentVector>::value,
- typename Rotation<N,ctype>::TangentVector,
- typename Rotation<N,ctype>::EmbeddedTangentVector>::Type RotationTangentVector;
+ typename Rotation<ctype,N>::TangentVector,
+ typename Rotation<ctype,N>::EmbeddedTangentVector>::Type RotationTangentVector;
RotationTangentVector qCorr;
for (int i=0; i<RotationTangentVector::dimension; i++)
qCorr[i] = v[N+i];
- result.q = Rotation<N,ctype>::exp(p.q, qCorr);
+ result.q = Rotation<ctype,N>::exp(p.q, qCorr);
return result;
}
/** \brief Compute geodesic distance from a to b */
- static T distance(const RigidBodyMotion<N,ctype>& a, const RigidBodyMotion<N,ctype>& b) {
+ static T distance(const RigidBodyMotion<ctype,N>& a, const RigidBodyMotion<ctype,N>& b) {
T euclideanDistanceSquared = (a.r - b.r).two_norm2();
- T rotationDistance = Rotation<N,ctype>::distance(a.q, b.q);
+ T rotationDistance = Rotation<ctype,N>::distance(a.q, b.q);
return std::sqrt(euclideanDistanceSquared + rotationDistance*rotationDistance);
}
/** \brief Compute difference vector from a to b on the tangent space of a */
- static TangentVector difference(const RigidBodyMotion<N,ctype>& a,
- const RigidBodyMotion<N,ctype>& b) {
+ static TangentVector difference(const RigidBodyMotion<ctype,N>& a,
+ const RigidBodyMotion<ctype,N>& b) {
TangentVector result;
@@ -100,17 +100,17 @@ public:
result[i] = a.r[i] - b.r[i];
// Subtract orientations on the tangent space of 'a'
- typename Rotation<N,ctype>::TangentVector v = Rotation<N,ctype>::difference(a.q, b.q).axial();
+ typename Rotation<ctype,N>::TangentVector v = Rotation<ctype,N>::difference(a.q, b.q).axial();
// Compute difference on T_a SO(3)
- for (int i=0; i<Rotation<N,ctype>::TangentVector::dimension; i++)
+ for (int i=0; i<Rotation<ctype,N>::TangentVector::dimension; i++)
result[i+N] = v[i];
return result;
}
- static EmbeddedTangentVector derivativeOfDistanceSquaredWRTSecondArgument(const RigidBodyMotion<N,ctype>& a,
- const RigidBodyMotion<N,ctype>& b) {
+ static EmbeddedTangentVector derivativeOfDistanceSquaredWRTSecondArgument(const RigidBodyMotion<ctype,N>& a,
+ const RigidBodyMotion<ctype,N>& b) {
// linear part
Dune::FieldVector<ctype,N> linearDerivative = a.r;
@@ -118,28 +118,28 @@ public:
linearDerivative *= -2;
// rotation part
- typename Rotation<N,ctype>::EmbeddedTangentVector rotationDerivative
- = Rotation<N,ctype>::derivativeOfDistanceSquaredWRTSecondArgument(a.q, b.q);
+ typename Rotation<ctype,N>::EmbeddedTangentVector rotationDerivative
+ = Rotation<ctype,N>::derivativeOfDistanceSquaredWRTSecondArgument(a.q, b.q);
return concat(linearDerivative, rotationDerivative);
}
/** \brief Compute the Hessian of the squared distance function keeping the first argument fixed */
- static Dune::FieldMatrix<T,embeddedDim,embeddedDim> secondDerivativeOfDistanceSquaredWRTSecondArgument(const RigidBodyMotion<N,ctype> & p, const RigidBodyMotion<N,ctype> & q)
+ static Dune::FieldMatrix<T,embeddedDim,embeddedDim> secondDerivativeOfDistanceSquaredWRTSecondArgument(const RigidBodyMotion<ctype,N> & p, const RigidBodyMotion<ctype,N> & q)
{
Dune::FieldMatrix<T,embeddedDim,embeddedDim> result(0);
// The linear part
- Dune::FieldMatrix<T,N,N> linearPart = RealTuple<N>::secondDerivativeOfDistanceSquaredWRTSecondArgument(p.r,q.r);
+ Dune::FieldMatrix<T,N,N> linearPart = RealTuple<T,N>::secondDerivativeOfDistanceSquaredWRTSecondArgument(p.r,q.r);
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
result[i][j] = linearPart[i][j];
// The rotation part
- Dune::FieldMatrix<T,Rotation<N,T>::embeddedDim,Rotation<N,T>::embeddedDim> rotationPart
- = Rotation<N,ctype>::secondDerivativeOfDistanceSquaredWRTSecondArgument(p.q,q.q);
- for (int i=0; i<Rotation<N,T>::embeddedDim; i++)
- for (int j=0; j<Rotation<N,T>::embeddedDim; j++)
+ Dune::FieldMatrix<T,Rotation<T,N>::embeddedDim,Rotation<T,N>::embeddedDim> rotationPart
+ = Rotation<ctype,N>::secondDerivativeOfDistanceSquaredWRTSecondArgument(p.q,q.q);
+ for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
+ for (int j=0; j<Rotation<T,N>::embeddedDim; j++)
result[N+i][N+j] = rotationPart[i][j];
return result;
@@ -149,21 +149,21 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
- static Dune::FieldMatrix<T,embeddedDim,embeddedDim> secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(const RigidBodyMotion<N,ctype> & p, const RigidBodyMotion<N,ctype> & q)
+ static Dune::FieldMatrix<T,embeddedDim,embeddedDim> secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(const RigidBodyMotion<ctype,N> & p, const RigidBodyMotion<ctype,N> & q)
{
Dune::FieldMatrix<T,embeddedDim,embeddedDim> result(0);
// The linear part
- Dune::FieldMatrix<T,N,N> linearPart = RealTuple<N>::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(p.r,q.r);
+ Dune::FieldMatrix<T,N,N> linearPart = RealTuple<T,N>::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(p.r,q.r);
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
result[i][j] = linearPart[i][j];
// The rotation part
- Dune::FieldMatrix<T,Rotation<N,T>::embeddedDim,Rotation<N,T>::embeddedDim> rotationPart
- = Rotation<N,ctype>::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(p.q,q.q);
- for (int i=0; i<Rotation<N,T>::embeddedDim; i++)
- for (int j=0; j<Rotation<N,T>::embeddedDim; j++)
+ Dune::FieldMatrix<T,Rotation<T,N>::embeddedDim,Rotation<T,N>::embeddedDim> rotationPart
+ = Rotation<ctype,N>::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(p.q,q.q);
+ for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
+ for (int j=0; j<Rotation<T,N>::embeddedDim; j++)
result[N+i][N+j] = rotationPart[i][j];
return result;
@@ -173,24 +173,24 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
- static Tensor3<T,embeddedDim,embeddedDim,embeddedDim> thirdDerivativeOfDistanceSquaredWRTSecondArgument(const RigidBodyMotion<N,ctype> & p, const RigidBodyMotion<N,ctype> & q)
+ static Tensor3<T,embeddedDim,embeddedDim,embeddedDim> thirdDerivativeOfDistanceSquaredWRTSecondArgument(const RigidBodyMotion<ctype,N> & p, const RigidBodyMotion<ctype,N> & q)
{
Tensor3<T,embeddedDim,embeddedDim,embeddedDim> result(0);
// The linear part
- Tensor3<T,N,N,N> linearPart = RealTuple<N>::thirdDerivativeOfDistanceSquaredWRTSecondArgument(p.r,q.r);
+ Tensor3<T,N,N,N> linearPart = RealTuple<T,N>::thirdDerivativeOfDistanceSquaredWRTSecondArgument(p.r,q.r);
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] = linearPart[i][j][k];
// The rotation part
- Tensor3<T,Rotation<N,T>::embeddedDim,Rotation<N,T>::embeddedDim,Rotation<N,T>::embeddedDim> rotationPart
- = Rotation<N,ctype>::thirdDerivativeOfDistanceSquaredWRTSecondArgument(p.q,q.q);
+ Tensor3<T,Rotation<T,N>::embeddedDim,Rotation<T,N>::embeddedDim,Rotation<T,N>::embeddedDim> rotationPart
+ = Rotation<ctype,N>::thirdDerivativeOfDistanceSquaredWRTSecondArgument(p.q,q.q);
- for (int i=0; i<Rotation<N,T>::embeddedDim; i++)
- for (int j=0; j<Rotation<N,T>::embeddedDim; j++)
- for (int k=0; k<Rotation<N,T>::embeddedDim; k++)
+ for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
+ for (int j=0; j<Rotation<T,N>::embeddedDim; j++)
+ for (int k=0; k<Rotation<T,N>::embeddedDim; k++)
result[N+i][N+j][N+k] = rotationPart[i][j][k];
return result;
@@ -200,22 +200,22 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
- static Tensor3<T,embeddedDim,embeddedDim,embeddedDim> thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(const RigidBodyMotion<N,ctype> & p, const RigidBodyMotion<N,ctype> & q)
+ static Tensor3<T,embeddedDim,embeddedDim,embeddedDim> thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(const RigidBodyMotion<ctype,N> & p, const RigidBodyMotion<ctype,N> & q)
{
Tensor3<T,embeddedDim,embeddedDim,embeddedDim> result(0);
// The linear part
- Tensor3<T,N,N,N> linearPart = RealTuple<N>::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(p.r,q.r);
+ Tensor3<T,N,N,N> linearPart = RealTuple<T,N>::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(p.r,q.r);
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] = linearPart[i][j][k];
// The rotation part
- Tensor3<T,Rotation<N,T>::embeddedDim,Rotation<N,T>::embeddedDim,Rotation<N,T>::embeddedDim> rotationPart = Rotation<N,ctype>::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(p.q,q.q);
- for (int i=0; i<Rotation<N,T>::embeddedDim; i++)
- for (int j=0; j<Rotation<N,T>::embeddedDim; j++)
- for (int k=0; k<Rotation<N,T>::embeddedDim; k++)
+ Tensor3<T,Rotation<T,N>::embeddedDim,Rotation<T,N>::embeddedDim,Rotation<T,N>::embeddedDim> rotationPart = Rotation<ctype,N>::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(p.q,q.q);
+ for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
+ for (int j=0; j<Rotation<T,N>::embeddedDim; j++)
+ for (int k=0; k<Rotation<T,N>::embeddedDim; k++)
result[N+i][N+j][N+k] = rotationPart[i][j][k];
return result;
@@ -240,10 +240,10 @@ public:
for (int i=0; i<N; i++)
result[i][i] = 1;
- Dune::FieldMatrix<T,Rotation<N>::dim,Rotation<N>::embeddedDim> SO3Part = q.orthonormalFrame();
+ Dune::FieldMatrix<T,Rotation<T,N>::dim,Rotation<T,N>::embeddedDim> SO3Part = q.orthonormalFrame();
- for (int i=0; i<Rotation<N>::dim; i++)
- for (int j=0; j<Rotation<N>::embeddedDim; j++)
+ for (int i=0; i<Rotation<T,N>::dim; i++)
+ for (int j=0; j<Rotation<T,N>::embeddedDim; j++)
result[N+i][N+j] = SO3Part[i][j];
return result;
@@ -260,7 +260,7 @@ public:
Dune::FieldVector<ctype, N> r;
// Rotational part
- Rotation<N,ctype> q;
+ Rotation<ctype,N> q;
private:
@@ -281,7 +281,7 @@ private:
//! Send configuration to output stream
template <int N, class ctype>
-std::ostream& operator<< (std::ostream& s, const RigidBodyMotion<N,ctype>& c)
+std::ostream& operator<< (std::ostream& s, const RigidBodyMotion<ctype,N>& c)
{
s << "(" << c.r << ") (" << c.q << ")";
return s;
diff --git a/dune/gfe/rotation.hh b/dune/gfe/rotation.hh
index 4b72a1d44b4f08a32fd8379300995d0932102541..398f9d78370c97328ff0383f3571995e90b1089c 100644
--- a/dune/gfe/rotation.hh
+++ b/dune/gfe/rotation.hh
@@ -15,7 +15,7 @@
#include <dune/gfe/unitvector.hh>
#include <dune/gfe/skewmatrix.hh>
-template <int dim, class T=double>
+template <class T, int dim>
class Rotation
{
@@ -25,7 +25,7 @@ class Rotation
\tparam T The type used for coordinates
*/
template <class T>
-class Rotation<2,T>
+class Rotation<T,2>
{
public:
/** \brief The type used for coordinates */
@@ -55,13 +55,13 @@ public:
{}
/** \brief Return the identity element */
- static Rotation<2,T> identity() {
+ static Rotation<T,2> identity() {
// Default constructor creates an identity
- Rotation<2,T> id;
+ Rotation<T,2> id;
return id;
}
- static T distance(const Rotation<2,T>& a, const Rotation<2,T>& b) {
+ static T distance(const Rotation<T,2>& a, const Rotation<T,2>& b) {
T dist = a.angle_ - b.angle_;
while (dist < 0)
dist += 2*M_PI;
@@ -72,36 +72,36 @@ public:
}
/** \brief The exponential map from a given point $p \in SO(3)$. */
- static Rotation<2,T> exp(const Rotation<2,T>& p, const TangentVector& v) {
- Rotation<2,T> result = p;
+ static Rotation<T,2> exp(const Rotation<T,2>& p, const TangentVector& v) {
+ Rotation<T,2> result = p;
result.angle_ += v;
return result;
}
/** \brief The exponential map from \f$ \mathfrak{so}(2) \f$ to \f$ SO(2) \f$
*/
- static Rotation<2,T> exp(const Dune::FieldVector<T,1>& v) {
- Rotation<2,T> result;
+ static Rotation<T,2> exp(const Dune::FieldVector<T,1>& v) {
+ Rotation<T,2> result;
result.angle_ = v[0];
return result;
}
- static TangentVector derivativeOfDistanceSquaredWRTSecondArgument(const Rotation<2,T>& a,
- const Rotation<2,T>& b) {
+ static TangentVector derivativeOfDistanceSquaredWRTSecondArgument(const Rotation<T,2>& a,
+ const Rotation<T,2>& b) {
// This assertion is here to remind me of the following laziness:
// The difference has to be computed modulo 2\pi
assert( std::fabs(a.angle_ - b.angle_) <= M_PI );
return -2 * (a.angle_ - b.angle_);
}
- static Dune::FieldMatrix<T,1,1> secondDerivativeOfDistanceSquaredWRTSecondArgument(const Rotation<2,T>& a,
- const Rotation<2,T>& b) {
+ static Dune::FieldMatrix<T,1,1> secondDerivativeOfDistanceSquaredWRTSecondArgument(const Rotation<T,2>& a,
+ const Rotation<T,2>& b) {
return 2;
}
/** \brief Right multiplication */
- Rotation<2,T> mult(const Rotation<2,T>& other) const {
- Rotation<2,T> q = *this;
+ Rotation<T,2> mult(const Rotation<T,2>& other) const {
+ Rotation<T,2> q = *this;
q.angle_ += other.angle_;
return q;
}
@@ -122,7 +122,7 @@ public:
//! Send configuration to output stream
template <class T>
-std::ostream& operator<< (std::ostream& s, const Rotation<2,T>& c)
+std::ostream& operator<< (std::ostream& s, const Rotation<T,2>& c)
{
return s << "[" << c.angle_ << " (" << std::sin(c.angle_) << " " << std::cos(c.angle_) << ") ]";
}
@@ -133,7 +133,7 @@ std::ostream& operator<< (std::ostream& s, const Rotation<2,T>& c)
Uses unit quaternion coordinates.
*/
template <class T>
-class Rotation<3,T> : public Quaternion<T>
+class Rotation<T,3> : public Quaternion<T>
{
/** \brief Computes sin(x/2) / x without getting unstable for small x */
@@ -177,7 +177,7 @@ public:
: Quaternion<T>(0,0,0,1)
{}
- Rotation<3,T>(const Dune::array<T,4>& c)
+ Rotation<T,3>(const Dune::array<T,4>& c)
{
for (int i=0; i<4; i++)
(*this)[i] = c[i];
@@ -185,26 +185,26 @@ public:
*this /= this->two_norm();
}
- Rotation<3,T>(const Dune::FieldVector<T,4>& c)
+ Rotation<T,3>(const Dune::FieldVector<T,4>& c)
: Quaternion<T>(c)
{
*this /= this->two_norm();
}
- Rotation<3,T>(Dune::FieldVector<T,3> axis, T angle)
+ Rotation<T,3>(Dune::FieldVector<T,3> axis, T angle)
: Quaternion<T>(axis, angle)
{}
/** \brief Return the identity element */
- static Rotation<3,T> identity() {
+ static Rotation<T,3> identity() {
// Default constructor creates an identity
- Rotation<3,T> id;
+ Rotation<T,3> id;
return id;
}
/** \brief Right multiplication */
- Rotation<3,T> mult(const Rotation<3,T>& other) const {
- Rotation<3,T> q;
+ Rotation<T,3> mult(const Rotation<T,3>& other) const {
+ Rotation<T,3> q;
q[0] = (*this)[3]*other[0] - (*this)[2]*other[1] + (*this)[1]*other[2] + (*this)[0]*other[3];
q[1] = (*this)[2]*other[0] + (*this)[3]*other[1] - (*this)[0]*other[2] + (*this)[1]*other[3];
q[2] = - (*this)[1]*other[0] + (*this)[0]*other[1] + (*this)[3]*other[2] + (*this)[2]*other[3];
@@ -215,8 +215,8 @@ public:
/** \brief The exponential map from \f$ \mathfrak{so}(3) \f$ to \f$ SO(3) \f$
*/
- static Rotation<3,T> exp(const SkewMatrix<T,3>& v) {
- Rotation<3,T> q;
+ static Rotation<T,3> exp(const SkewMatrix<T,3>& v) {
+ Rotation<T,3> q;
Dune::FieldVector<T,3> vAxial = v.axial();
T normV = vAxial.two_norm();
@@ -237,8 +237,8 @@ public:
/** \brief The exponential map from a given point $p \in SO(3)$. */
- static Rotation<3,T> exp(const Rotation<3,T>& p, const SkewMatrix<T,3>& v) {
- Rotation<3,T> corr = exp(v);
+ static Rotation<T,3> exp(const Rotation<T,3>& p, const SkewMatrix<T,3>& v) {
+ Rotation<T,3> corr = exp(v);
return p.mult(corr);
}
@@ -248,7 +248,7 @@ public:
\param v A tangent vector in quaternion coordinates
*/
- static Rotation<3,T> exp(const Rotation<3,T>& p, const EmbeddedTangentVector& v) {
+ static Rotation<T,3> exp(const Rotation<T,3>& p, const EmbeddedTangentVector& v) {
assert( std::fabs(p*v) < 1e-8 );
@@ -272,7 +272,7 @@ public:
\param v A tangent vector.
*/
- static Rotation<3,T> exp(const Rotation<3,T>& p, const TangentVector& v) {
+ static Rotation<T,3> exp(const Rotation<T,3>& p, const TangentVector& v) {
// embedded tangent vector
Dune::FieldMatrix<T,3,4> basis = p.orthonormalFrame();
@@ -284,7 +284,7 @@ public:
}
/** \brief Compute tangent vector from given basepoint and skew symmetric matrix. */
- static TangentVector skewToTangentVector(const Rotation<3,T>& p, const SkewMatrix<T,3>& v ) {
+ static TangentVector skewToTangentVector(const Rotation<T,3>& p, const SkewMatrix<T,3>& v ) {
// embedded tangent vector at identity
Quaternion<T> vAtIdentity(0);
@@ -306,7 +306,7 @@ public:
}
/** \brief Compute skew matrix from given basepoint and tangent vector. */
- static SkewMatrix<T,3> tangentToSkew(const Rotation<3,T>& p, const TangentVector& tangent) {
+ static SkewMatrix<T,3> tangentToSkew(const Rotation<T,3>& p, const TangentVector& tangent) {
// embedded tangent vector
Dune::FieldMatrix<T,3,4> basis = p.orthonormalFrame();
@@ -317,7 +317,7 @@ public:
}
/** \brief Compute skew matrix from given basepoint and an embedded tangent vector. */
- static SkewMatrix<T,3> tangentToSkew(const Rotation<3,T>& p, const EmbeddedTangentVector& q) {
+ static SkewMatrix<T,3> tangentToSkew(const Rotation<T,3>& p, const EmbeddedTangentVector& q) {
// left multiplication by the inverse base point yields a tangent vector at the identity
Quaternion<T> vAtIdentity = p.inverse().mult(q);
@@ -331,7 +331,7 @@ public:
return skew;
}
- static Rotation<3,T> exp(const Rotation<3,T>& p, const Dune::FieldVector<T,4>& v) {
+ static Rotation<T,3> exp(const Rotation<T,3>& p, const Dune::FieldVector<T,4>& v) {
assert( std::fabs(p*v) < 1e-8 );
@@ -415,7 +415,7 @@ public:
}
/** \brief The inverse of the exponential map */
- static Dune::FieldVector<T,3> expInv(const Rotation<3,T>& q) {
+ static Dune::FieldVector<T,3> expInv(const Rotation<T,3>& q) {
// Compute v = exp^{-1} q
// Due to numerical dirt, q[3] may be larger than 1.
// In that case, use 1 instead of q[3].
@@ -437,7 +437,7 @@ public:
}
/** \brief The derivative of the inverse of the exponential map, evaluated at q */
- static Dune::FieldMatrix<T,3,4> DexpInv(const Rotation<3,T>& q) {
+ static Dune::FieldMatrix<T,3,4> DexpInv(const Rotation<T,3>& q) {
// Compute v = exp^{-1} q
Dune::FieldVector<T,3> v = expInv(q);
@@ -474,8 +474,8 @@ public:
* three-dimensional rods:Exact energy and momentum conserving algorithms'
* (but the corrected version with 0.25 instead of 0.5 in the denominator)
*/
- static Rotation<3,T> cayley(const SkewMatrix<T,3>& s) {
- Rotation<3,T> q;
+ static Rotation<T,3> cayley(const SkewMatrix<T,3>& s) {
+ Rotation<T,3> q;
Dune::FieldVector<T,3> vAxial = s.axial();
T norm = 0.25*vAxial.two_norm2() + 1;
@@ -498,7 +498,7 @@ public:
*
* The formula is taken from J.M.Selig - Cayley Maps for SE(3).
*/
- static SkewMatrix<T,3> cayleyInv(const Rotation<3,T> q) {
+ static SkewMatrix<T,3> cayleyInv(const Rotation<T,3> q) {
Dune::FieldMatrix<T,3,3> mat;
@@ -509,7 +509,7 @@ public:
q.matrix(mat);
Dune::FieldMatrix<T,3,3> matT;
- Rotation<3,T>(q.inverse()).matrix(matT);
+ Rotation<T,3>(q.inverse()).matrix(matT);
mat -= matT;
mat *= 2/(1+trace);
}
@@ -537,7 +537,7 @@ public:
}
- static T distance(const Rotation<3,T>& a, const Rotation<3,T>& b) {
+ static T distance(const Rotation<T,3>& a, const Rotation<T,3>& b) {
Quaternion<T> diff = a;
diff.invert();
@@ -559,7 +559,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<3,T>& a, const Rotation<3,T>& b) {
+ static SkewMatrix<T,3> difference(const Rotation<T,3>& a, const Rotation<T,3>& b) {
Quaternion<T> diff = a;
diff.invert();
@@ -620,10 +620,10 @@ public:
}
- static EmbeddedTangentVector derivativeOfDistanceSquaredWRTSecondArgument(const Rotation<3,T>& p,
- const Rotation<3,T>& q) {
+ static EmbeddedTangentVector derivativeOfDistanceSquaredWRTSecondArgument(const Rotation<T,3>& p,
+ const Rotation<T,3>& q) {
- Rotation<3,T> pInv = p;
+ Rotation<T,3> pInv = p;
pInv.invert();
// the forth component of pInv times q
@@ -650,9 +650,9 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
- static Dune::FieldMatrix<T,4,4> secondDerivativeOfDistanceSquaredWRTSecondArgument(const Rotation<3,T>& p, const Rotation<3,T>& q) {
+ static Dune::FieldMatrix<T,4,4> secondDerivativeOfDistanceSquaredWRTSecondArgument(const Rotation<T,3>& p, const Rotation<T,3>& q) {
// use the functionality from the unitvector class
- Dune::FieldMatrix<T,4,4> result = UnitVector<4>::secondDerivativeOfDistanceSquaredWRTSecondArgument(p.globalCoordinates(),
+ Dune::FieldMatrix<T,4,4> result = UnitVector<T,4>::secondDerivativeOfDistanceSquaredWRTSecondArgument(p.globalCoordinates(),
q.globalCoordinates());
// for some reason that I don't really understand, the distance we have defined for the rotations (== Unit quaternions)
// is twice the corresponding distance on the unit quaternions seen as a sphere. Hence the derivative of the
@@ -665,9 +665,9 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
- static Dune::FieldMatrix<T,4,4> secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(const Rotation<3,T>& p, const Rotation<3,T>& q) {
+ static Dune::FieldMatrix<T,4,4> secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(const Rotation<T,3>& p, const Rotation<T,3>& q) {
// use the functionality from the unitvector class
- Dune::FieldMatrix<T,4,4> result = UnitVector<4>::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(p.globalCoordinates(),
+ Dune::FieldMatrix<T,4,4> result = UnitVector<T,4>::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(p.globalCoordinates(),
q.globalCoordinates());
// for some reason that I don't really understand, the distance we have defined for the rotations (== Unit quaternions)
// is twice the corresponding distance on the unit quaternions seen as a sphere. Hence the derivative of the
@@ -680,9 +680,9 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
- static Tensor3<T,4,4,4> thirdDerivativeOfDistanceSquaredWRTSecondArgument(const Rotation<3,T>& p, const Rotation<3,T>& q) {
+ static Tensor3<T,4,4,4> thirdDerivativeOfDistanceSquaredWRTSecondArgument(const Rotation<T,3>& p, const Rotation<T,3>& q) {
// use the functionality from the unitvector class
- Tensor3<T,4,4,4> result = UnitVector<4>::thirdDerivativeOfDistanceSquaredWRTSecondArgument(p.globalCoordinates(),
+ Tensor3<T,4,4,4> result = UnitVector<T,4>::thirdDerivativeOfDistanceSquaredWRTSecondArgument(p.globalCoordinates(),
q.globalCoordinates());
// for some reason that I don't really understand, the distance we have defined for the rotations (== Unit quaternions)
// is twice the corresponding distance on the unit quaternions seen as a sphere. Hence the derivative of the
@@ -695,9 +695,9 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
- static Tensor3<T,4,4,4> thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(const Rotation<3,T>& p, const Rotation<3,T>& q) {
+ static Tensor3<T,4,4,4> thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(const Rotation<T,3>& p, const Rotation<T,3>& q) {
// use the functionality from the unitvector class
- Tensor3<T,4,4,4> result = UnitVector<4>::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(p.globalCoordinates(),
+ Tensor3<T,4,4,4> result = UnitVector<T,4>::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(p.globalCoordinates(),
q.globalCoordinates());
// for some reason that I don't really understand, the distance we have defined for the rotations (== Unit quaternions)
// is twice the corresponding distance on the unit quaternions seen as a sphere. Hence the derivative of the
@@ -710,7 +710,7 @@ public:
/** \brief Interpolate between two rotations */
- static Rotation<3,T> interpolate(const Rotation<3,T>& a, const Rotation<3,T>& b, T omega) {
+ 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);
@@ -723,7 +723,7 @@ public:
/** \brief Interpolate between two rotations
\param omega must be between 0 and 1
*/
- static Quaternion<T> interpolateDerivative(const Rotation<3,T>& a, const Rotation<3,T>& b,
+ static Quaternion<T> interpolateDerivative(const Rotation<T,3>& a, const Rotation<T,3>& b,
T omega) {
Quaternion<T> result(0);
diff --git a/dune/gfe/unitvector.hh b/dune/gfe/unitvector.hh
index edce122d33876dd594652cdcd613b75d65e888bd..7be1a3daee603ea7c369d19f12348fc467fff659 100644
--- a/dune/gfe/unitvector.hh
+++ b/dune/gfe/unitvector.hh
@@ -11,7 +11,7 @@
\tparam N Dimension of the embedding space
\tparam T The type used for individual coordinates
*/
-template <int N, class T=double>
+template <class T, int N>
class UnitVector
{
/** \brief Computes sin(x) / x without getting unstable for small x */
@@ -91,7 +91,7 @@ public:
data_ /= data_.two_norm();
}
- UnitVector<N>& operator=(const Dune::FieldVector<T,N>& vector)
+ UnitVector<T,N>& operator=(const Dune::FieldVector<T,N>& vector)
{
data_ = vector;
data_ /= data_.two_norm();