From 7f7e26b91927c84ccfbf55de3c7efe16f4714261 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sat, 12 Nov 2011 19:21:10 +0000
Subject: [PATCH] Remove the default parameter for the coordinate type. It
makes it difficult to be sure that an explicitly requested type is really
handed over to every place.
While we are at it: reverse the order of the template
parameters. Now first comes the type, then the dimension.
This is more like FieldVector and friends, and hence
should be less confusing.
[[Imported from SVN: r8148]]
---
dune/gfe/realtuple.hh | 4 +-
dune/gfe/rigidbodymotion.hh | 98 ++++++++++++++++----------------
dune/gfe/rotation.hh | 110 ++++++++++++++++++------------------
dune/gfe/unitvector.hh | 4 +-
4 files changed, 108 insertions(+), 108 deletions(-)
diff --git a/dune/gfe/realtuple.hh b/dune/gfe/realtuple.hh
index f572f96e..99fd8dbb 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 2d1194f4..8945d178 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 4b72a1d4..398f9d78 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 edce122d..7be1a3da 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();
--
GitLab