diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
index 0000000000000000000000000000000000000000..b7a71f58d5c52cc4112fddb1bc74a5f80b2965f9
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,18 @@
+# Master
+
+- Fix bug in the `RealTuple::log` method: Calling `log(a,b)`returned `a-b`
+ instead of `b-a`.
+
+- Fix the return value of `ProductManifold::log`: It was `TangentVector`,
+ but now it is `EmbeddedTangentVector`.
+
+- The `RigidBodyMotion` class has a `log` method now.
+
+- The method `Rotation<3>::log` now returns an `EmbeddedTangentVector`
+ instead of a `SkewMatrix`. This is consistent with the other manifold
+ implementations.
+
+- Deprecate the method `RigidBodyMotion::difference`; the method
+ `RigidBodyMotion::log`. Watch out: The `difference` method was buggy!
+ See https://gitlab.mn.tu-dresden.de/osander/dune-gfe/-/merge_requests/2
+ for details.
diff --git a/dune/gfe/polardecomposition.hh b/dune/gfe/polardecomposition.hh
index 6f3ecb862113d5e241bc267d0bdd701ec3e43406..7d48b9bee1698cecd6d67a064ce71d87ad8112f0 100644
--- a/dune/gfe/polardecomposition.hh
+++ b/dune/gfe/polardecomposition.hh
@@ -25,7 +25,7 @@ namespace Dune::GFE {
template <class field_type>
FieldMatrix<field_type,3,3> operator() (const FieldMatrix<field_type,3,3>& matrix, double tol = 0.001) const
{
- int maxIterations = 100;
+ size_t maxIterations = 100;
// Use Higham's method
auto polar = matrix;
for (size_t i=0; i<maxIterations; i++)
@@ -244,13 +244,13 @@ namespace Dune::GFE {
mini = i;
Pleft[3][mini] = 1;
- Pright[mini][3] = 1; // Smalest element to the last position
+ Pright[mini][3] = 1; // Smallest element to the last position
for (int i = 0; i < 4; ++i) {
- if ( i != maxi & i != mini ) {
+ if ( i != maxi && i != mini ) {
for (int j = 0; j < 4; ++j) {
- if ( j != maxi & j != mini & j != i ) {
+ if ( j != maxi && j != mini && j != i ) {
if ( Bdiag[i] < Bdiag[j] ) {
Pleft[1][j] = 1; // Second largest element at the second position
Pright[j][1] = 1;
@@ -323,4 +323,4 @@ namespace Dune::GFE {
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/dune/gfe/productmanifold.hh b/dune/gfe/productmanifold.hh
index 48e86df8a9d44ce355f9ab6715a91b9c8b5f8d3f..ebecd8d414115f266ea05b71792e999a6293f00c 100644
--- a/dune/gfe/productmanifold.hh
+++ b/dune/gfe/productmanifold.hh
@@ -161,9 +161,9 @@ namespace Dune::GFE
}
/** \brief Compute difference vector from a to b on the tangent space of a */
- static TangentVector log(const ProductManifold<TS,TargetSpaces...>& a, const ProductManifold<TS,TargetSpaces...>& b)
+ static EmbeddedTangentVector log(const ProductManifold<TS,TargetSpaces...>& a, const ProductManifold<TS,TargetSpaces...>& b)
{
- TangentVector diff;
+ EmbeddedTangentVector diff;
auto logFunctor =[] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
{
auto& res = std::get<0>(argsTuple);
@@ -172,7 +172,7 @@ namespace Dune::GFE
using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(a)> >;
const auto diffLoc = Manifold::log(a,b);
std::copy(diffLoc.begin(),diffLoc.end(),res.begin()+posHelper[0]);
- posHelper[0] += Manifold::dim;
+ posHelper[0] += Manifold::embeddedDim;
};
foreachManifold(logFunctor,diff,a, b);
return diff;
diff --git a/dune/gfe/realtuple.hh b/dune/gfe/realtuple.hh
index e4bc4df8561ab9c5fff83bd529b29c7ffe97f61d..7f28cc851254a8ec6981ffda2d3ebbcf067ebb50 100644
--- a/dune/gfe/realtuple.hh
+++ b/dune/gfe/realtuple.hh
@@ -110,10 +110,11 @@ public:
}
/** \brief The logarithmic map
- * Simply the difference vector for RealTuple
+ *
+ * \result A vector in the tangent space of a, viz: b-a
* */
static auto log(const RealTuple& a, const RealTuple& b) {
- return a.data_ - b.data_;
+ return b.data_ - a.data_;
}
#if ADOLC_ADOUBLE_H
diff --git a/dune/gfe/rigidbodymotion.hh b/dune/gfe/rigidbodymotion.hh
index b4541e3e722dac0d7ec67098771043d1e2ee863b..4a4704ff7bd9072981eeaf6535dda7eb477422ab 100644
--- a/dune/gfe/rigidbodymotion.hh
+++ b/dune/gfe/rigidbodymotion.hh
@@ -99,6 +99,26 @@ public:
return result;
}
+ /** \brief Compute difference vector from a to b on the tangent space of a */
+ static EmbeddedTangentVector log(const RigidBodyMotion<ctype,N>& a,
+ const RigidBodyMotion<ctype,N>& b)
+ {
+ EmbeddedTangentVector result;
+
+ // Usual linear difference
+ for (int i=0; i<N; i++)
+ result[i] = b.r[i] - a.r[i];
+
+ // Subtract orientations on the tangent space of 'a'
+ typename Rotation<ctype,N>::EmbeddedTangentVector v = Rotation<ctype,N>::log(a.q, b.q);
+
+ // Compute difference on T_a SO(3)
+ for (int i=0; i<Rotation<ctype,N>::EmbeddedTangentVector::dimension; i++)
+ result[i+N] = v[i];
+
+ return result;
+ }
+
/** \brief Compute geodesic distance from a to b */
static T distance(const RigidBodyMotion<ctype,N>& a, const RigidBodyMotion<ctype,N>& b) {
@@ -109,7 +129,12 @@ public:
return std::sqrt(euclideanDistanceSquared + rotationDistance*rotationDistance);
}
- /** \brief Compute difference vector from a to b on the tangent space of a */
+ /** \brief Compute difference vector from a to b on the tangent space of a
+
+ * \warning The method is buggy! See https://gitlab.mn.tu-dresden.de/osander/dune-gfe/-/merge_requests/2
+ * \deprecated Use the log method instead!
+ */
+ [[deprecated("Use RigidBodyMotion::log instead of RigidBodyMotion::difference!")]]
static TangentVector difference(const RigidBodyMotion<ctype,N>& a,
const RigidBodyMotion<ctype,N>& b) {
@@ -145,22 +170,26 @@ public:
}
/** \brief Compute the Hessian of the squared distance function keeping the first argument fixed */
- static Dune::FieldMatrix<T,embeddedDim,embeddedDim> secondDerivativeOfDistanceSquaredWRTSecondArgument(const RigidBodyMotion<ctype,N> & p, const RigidBodyMotion<ctype,N> & q)
+ static Dune::SymmetricMatrix<T,embeddedDim> secondDerivativeOfDistanceSquaredWRTSecondArgument(const RigidBodyMotion<ctype,N> & p, const RigidBodyMotion<ctype,N> & q)
{
- Dune::FieldMatrix<T,embeddedDim,embeddedDim> result(0);
+ Dune::SymmetricMatrix<T,embeddedDim> result;
// The linear part
- Dune::FieldMatrix<T,N,N> linearPart = RealTuple<T,N>::secondDerivativeOfDistanceSquaredWRTSecondArgument(p.r,q.r);
+ Dune::SymmetricMatrix<T,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];
+ for (int j=0; j<=i; j++)
+ result(i,j) = linearPart(i,j);
// The rotation part
- Dune::FieldMatrix<T,Rotation<T,N>::embeddedDim,Rotation<T,N>::embeddedDim> rotationPart
+ Dune::SymmetricMatrix<T,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];
+ {
+ for (int j=0; j<N; j++)
+ result(N+i,j) = 0;
+ for (int j=0; j<=i; j++)
+ result(N+i,N+j) = rotationPart(i,j);
+ }
return result;
}
diff --git a/dune/gfe/rotation.hh b/dune/gfe/rotation.hh
index 6a1cc01a69ab3be0d888c0981fbdf0c28ac2d0c3..139a30344bb755601c479688ccdce6b43233434d 100644
--- a/dune/gfe/rotation.hh
+++ b/dune/gfe/rotation.hh
@@ -305,7 +305,10 @@ public:
}
- /** \brief The exponential map from a given point $p \in SO(3)$. */
+ /** \brief The exponential map from a given point $p \in SO(3)$.
+
+ * \param v A tangent vector *at the identity*!
+ */
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);
@@ -403,6 +406,12 @@ public:
return skew;
}
+ /** \brief Derivative of the exponential map at the identity
+ *
+ * The exponential map at the identity is a map from a neighborhood of zero to the neighborhood of the identity rotation.
+ *
+ * \param v Where to evaluate the derivative of the (exponential map at the identity)
+ */
static Dune::FieldMatrix<T,4,3> Dexp(const SkewMatrix<T,3>& v) {
using std::cos;
@@ -613,7 +622,10 @@ 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> log(const Rotation<T,3>& a, const Rotation<T,3>& b) {
+ static EmbeddedTangentVector log(const Rotation<T,3>& a, const Rotation<T,3>& b) {
+
+ // embedded tangent vector at identity
+ Quaternion<T> v;
Quaternion<T> diff = a;
diff.invert();
@@ -622,7 +634,6 @@ public:
// Compute the geodesical distance between a and b on SO(3)
// Due to numerical dirt, diff[3] may be larger than 1.
// In that case, use 1 instead of diff[3].
- Dune::FieldVector<T,3> v;
if (diff[3] > 1.0) {
v = 0;
@@ -645,13 +656,14 @@ public:
T invSinc = 1/sincHalf(dist);
// Compute log on T_a SO(3)
- v[0] = diff[0] * invSinc;
- v[1] = diff[1] * invSinc;
- v[2] = diff[2] * invSinc;
-
+ v[0] = 0.5 * diff[0] * invSinc;
+ v[1] = 0.5 * diff[1] * invSinc;
+ v[2] = 0.5 * diff[2] * invSinc;
+ v[3] = 0;
}
- return SkewMatrix<T,3>(v);
+ // multiply with base point to get real embedded tangent vector
+ return ((Quaternion<T>) a).mult(v);
}
/** \brief Compute the derivatives of the director vectors with respect to the quaternion coordinates
@@ -917,11 +929,8 @@ public:
static Rotation<T,3> interpolate(const Rotation<T,3>& a, const Rotation<T,3>& b, T omega) {
// Compute log on T_a SO(3)
- SkewMatrix<T,3> v = log(a,b);
-
- v *= omega;
-
- return a.mult(exp(v));
+ EmbeddedTangentVector v = log(a,b);
+ return exp(a, omega*v);
}
/** \brief Interpolate between two rotations
diff --git a/dune/gfe/unitvector.hh b/dune/gfe/unitvector.hh
index 7122d1cda48f90638b8072f0211a6f1482a40baf..7a41e8435a4711db222ed55705e2cd4a57ff3717 100644
--- a/dune/gfe/unitvector.hh
+++ b/dune/gfe/unitvector.hh
@@ -219,6 +219,10 @@ public:
return result;
}
+ /** \brief The inverse of the exponential map
+ *
+ * \results A vector in the tangent space of p
+ */
static EmbeddedTangentVector log(const UnitVector& p, const UnitVector& q)
{
EmbeddedTangentVector result = p.projectOntoTangentSpace(q.data_-p.data_);
diff --git a/test/cosseratenergytest.cc b/test/cosseratenergytest.cc
index cdfc0f38c9d95afdf49e1aadd9cb380d69296214..9f4b5bc40e667f1355bf4531f6c25b5d889ecd44 100644
--- a/test/cosseratenergytest.cc
+++ b/test/cosseratenergytest.cc
@@ -7,8 +7,6 @@
#include <dune/functions/functionspacebases/lagrangebasis.hh>
-#include <dune/fufem/functions/constantfunction.hh>
-
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/cosseratenergystiffness.hh>
@@ -44,7 +42,7 @@ std::unique_ptr<GridType> makeSingleSimplexGrid()
v[i] = i;
factory.insertElement(GeometryTypes::simplex(domainDim), v);
- return std::unique_ptr<GridType>(factory.createGrid());
+ return factory.createGrid();
}
diff --git a/test/cosseratrodenergytest.cc b/test/cosseratrodenergytest.cc
index ff65271b6c6ed816cded185956c0158e4b1a1fd9..6a485ac5704538ca516071b8111c21702b149d0d 100644
--- a/test/cosseratrodenergytest.cc
+++ b/test/cosseratrodenergytest.cc
@@ -80,7 +80,7 @@ int main (int argc, char *argv[]) try
std::vector<RigidBodyMotion<double,3> > referenceConfiguration(gridView.size(1));
- for (const auto vertex : vertices(gridView))
+ for (const auto& vertex : vertices(gridView))
{
auto idx = gridView.indexSet().index(vertex);
diff --git a/test/polardecompositiontest.cc b/test/polardecompositiontest.cc
index cc10676d24424d59fc4316dace27efa2a9fb9720..779d929c07a21faa610704e10dab8f52b73abc02 100644
--- a/test/polardecompositiontest.cc
+++ b/test/polardecompositiontest.cc
@@ -168,7 +168,6 @@ static bool test4dDeterminant() {
static bool testdominantEVNewton() {
field_type detB = 0.992928;
- field_type detM = 0.000620104;
FieldVector<field_type,3> minpol = {0, -1.99999, -0.00496083};
double correctDominantEV = 1.05397;
diff --git a/test/surfacecosseratstressassemblertest.cc b/test/surfacecosseratstressassemblertest.cc
index 66cbddd561d31042533a57896b35b9730d5481e3..9e2ea9609c73b85646a9be8032e827ab63ab3754 100644
--- a/test/surfacecosseratstressassemblertest.cc
+++ b/test/surfacecosseratstressassemblertest.cc
@@ -92,7 +92,7 @@ static bool symmetryTest(std::unordered_map<std::string, FieldVector<double,d>>
int main (int argc, char *argv[])
{
- MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
+ MPIHelper::instance(argc, argv);
/////////////////////////////////////////////////////////////
// CREATE THE GRID
@@ -186,7 +186,7 @@ int main (int argc, char *argv[])
Functions::interpolate(basisOrderD, x, [](FieldVector<double,dim> x){ return x; });
Functions::interpolate(basisOrderD, xInitial, [](FieldVector<double,dim> x){ return x; });
- for (int i = 0; i < basisOrderD.size(); i++) {
+ for (std::size_t i = 0; i < basisOrderD.size(); i++) {
std::stringstream stream;
stream << x[i];
//Look up the displacement for this vertex in the deformationMap
@@ -201,7 +201,7 @@ int main (int argc, char *argv[])
DisplacementVector xOrderR;
xOrderR.resize(basisOrderR.size());
Functions::interpolate(basisOrderR, xOrderR, [](FieldVector<double,dim> x){ return x; });
- for (int i = 0; i < basisOrderR.size(); i++) {
+ for (std::size_t i = 0; i < basisOrderR.size(); i++) {
std::stringstream stream;
stream << xOrderR[i];
Rotation<double,dim> rotation(rotationMap.at(stream.str()));
diff --git a/test/targetspacetest.cc b/test/targetspacetest.cc
index b5e1c650bf1e91361fb1bf8fd7111ac40414bff8..cea7382c03b2361660487a408401ba17e9c755b5 100644
--- a/test/targetspacetest.cc
+++ b/test/targetspacetest.cc
@@ -29,23 +29,33 @@ double diameter(const std::vector<TargetSpace>& v)
const double eps = 1e-4;
template <class TargetSpace>
-double energy(const TargetSpace& a, const TargetSpace& b)
+void testExpLog(const TargetSpace& a, const TargetSpace& b)
{
- return TargetSpace::distance(a,b) * TargetSpace::distance(a,b);
+ // Check whether exp and log are mutually inverse
+ typename TargetSpace::EmbeddedTangentVector logarithm = TargetSpace::log(a,b);
+ TargetSpace exponential = TargetSpace::exp(a, logarithm);
+
+ if (TargetSpace::distance(b, exponential) > eps)
+ {
+ std::cout << className(a) << ": Exp and log are not mutually inverse." << std::endl;
+ std::cout << "exp(a,log(a,b)): " << exponential << std::endl;
+ std::cout << "b : " << b << std::endl;
+ assert(false);
+ }
}
-template <class TargetSpace, int dim>
-double energy(const TargetSpace& a, const FieldVector<double,dim>& b)
+template <class TargetSpace>
+double distanceSquared(const TargetSpace& a, const TargetSpace& b)
{
-#warning Cast where there should not be one
- return TargetSpace::distance(a,TargetSpace(b)) * TargetSpace::distance(a,TargetSpace(b));
+ return Dune::power(TargetSpace::distance(a,b), 2);
}
+// Squared distance between two points slightly off the manifold.
+// This is required for finite difference approximations.
template <class TargetSpace, int dim>
-double energy(const FieldVector<double,dim>& a, const FieldVector<double,dim>& b)
+double distanceSquared(const FieldVector<double,dim>& a, const FieldVector<double,dim>& b)
{
-#warning Cast where there should not be one
- return TargetSpace::distance(TargetSpace(a),TargetSpace(b)) * TargetSpace::distance(TargetSpace(a),TargetSpace(b));
+ return Dune::power(TargetSpace::distance(TargetSpace(a),TargetSpace(b)), 2);
}
/** \brief Compute the Riemannian Hessian of the squared distance function in global coordinates
@@ -67,15 +77,15 @@ FieldMatrix<double,worldDim,worldDim> getSecondDerivativeOfSecondArgumentFD(cons
for (size_t i=0; i<spaceDim; i++) {
for (size_t j=0; j<spaceDim; j++) {
- FieldVector<double,worldDim> epsXi = B[i]; epsXi *= eps;
- FieldVector<double,worldDim> epsEta = B[j]; epsEta *= eps;
+ FieldVector<double,worldDim> epsXi = eps * B[i];
+ FieldVector<double,worldDim> epsEta = eps * B[j];
- FieldVector<double,worldDim> minusEpsXi = epsXi; minusEpsXi *= -1;
- FieldVector<double,worldDim> minusEpsEta = epsEta; minusEpsEta *= -1;
+ FieldVector<double,worldDim> minusEpsXi = -1 * epsXi;
+ FieldVector<double,worldDim> minusEpsEta = -1 * epsEta;
- double forwardValue = energy(a,TargetSpace::exp(b,epsXi+epsEta)) - energy(a, TargetSpace::exp(b,epsXi)) - energy(a,TargetSpace::exp(b,epsEta));
- double centerValue = energy(a,b) - energy(a,b) - energy(a,b);
- double backwardValue = energy(a,TargetSpace::exp(b,minusEpsXi + minusEpsEta)) - energy(a, TargetSpace::exp(b,minusEpsXi)) - energy(a,TargetSpace::exp(b,minusEpsEta));
+ double forwardValue = distanceSquared(a,TargetSpace::exp(b,epsXi+epsEta)) - distanceSquared(a, TargetSpace::exp(b,epsXi)) - distanceSquared(a,TargetSpace::exp(b,epsEta));
+ double centerValue = distanceSquared(a,b) - distanceSquared(a,b) - distanceSquared(a,b);
+ double backwardValue = distanceSquared(a,TargetSpace::exp(b,minusEpsXi + minusEpsEta)) - distanceSquared(a, TargetSpace::exp(b,minusEpsXi)) - distanceSquared(a,TargetSpace::exp(b,minusEpsEta));
d2d2_fd[i][j] = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps);
@@ -111,7 +121,7 @@ void testOrthonormalFrame(const TargetSpace& a)
}
template <class TargetSpace>
-void testDerivativeOfSquaredDistance(const TargetSpace& a, const TargetSpace& b)
+void testDerivativeOfDistanceSquared(const TargetSpace& a, const TargetSpace& b)
{
static const size_t embeddedDim = TargetSpace::embeddedDim;
@@ -126,14 +136,12 @@ void testDerivativeOfSquaredDistance(const TargetSpace& a, const TargetSpace& b)
typename TargetSpace::TangentVector d2_fd;
for (size_t i=0; i<TargetSpace::TangentVector::dimension; i++) {
- typename TargetSpace::EmbeddedTangentVector fwVariation = B[i];
- typename TargetSpace::EmbeddedTangentVector bwVariation = B[i];
- fwVariation *= eps;
- bwVariation *= -eps;
+ typename TargetSpace::EmbeddedTangentVector fwVariation = eps * B[i];
+ typename TargetSpace::EmbeddedTangentVector bwVariation = -eps * B[i];
TargetSpace bPlus = TargetSpace::exp(b,fwVariation);
TargetSpace bMinus = TargetSpace::exp(b,bwVariation);
- d2_fd[i] = (energy(a,bPlus) - energy(a,bMinus)) / (2*eps);
+ d2_fd[i] = (distanceSquared(a,bPlus) - distanceSquared(a,bMinus)) / (2*eps);
}
// transform into embedded coordinates
@@ -150,7 +158,7 @@ void testDerivativeOfSquaredDistance(const TargetSpace& a, const TargetSpace& b)
}
template <class TargetSpace>
-void testHessianOfSquaredDistance(const TargetSpace& a, const TargetSpace& b)
+void testHessianOfDistanceSquared(const TargetSpace& a, const TargetSpace& b)
{
static const int embeddedDim = TargetSpace::embeddedDim;
@@ -162,9 +170,7 @@ void testHessianOfSquaredDistance(const TargetSpace& a, const TargetSpace& b)
// finite-difference approximation
FieldMatrix<double,embeddedDim,embeddedDim> d2d2_fd = getSecondDerivativeOfSecondArgumentFD<TargetSpace,embeddedDim>(a,b);
- FieldMatrix<double,embeddedDim,embeddedDim> d2d2_diff = d2d2;
- d2d2_diff -= d2d2_fd;
- if ( (d2d2_diff).infinity_norm() > 200*eps) {
+ if ( (d2d2 - d2d2_fd).infinity_norm() > 200*eps) {
std::cout << className(a) << ": Analytical second derivative does not match fd approximation." << std::endl;
std::cout << "d2d2 Analytical:" << std::endl << d2d2 << std::endl;
std::cout << "d2d2 FD :" << std::endl << d2d2_fd << std::endl;
@@ -174,7 +180,7 @@ void testHessianOfSquaredDistance(const TargetSpace& a, const TargetSpace& b)
}
template <class TargetSpace>
-void testMixedDerivativesOfSquaredDistance(const TargetSpace& a, const TargetSpace& b)
+void testMixedDerivativesOfDistanceSquared(const TargetSpace& a, const TargetSpace& b)
{
static const size_t embeddedDim = TargetSpace::embeddedDim;
@@ -200,16 +206,13 @@ void testMixedDerivativesOfSquaredDistance(const TargetSpace& a, const TargetSpa
bPlus[j] += eps;
bMinus[j] -= eps;
- d1d2_fd[i][j] = (energy<TargetSpace>(aPlus,bPlus) + energy<TargetSpace>(aMinus,bMinus)
- - energy<TargetSpace>(aPlus,bMinus) - energy<TargetSpace>(aMinus,bPlus)) / (4*eps*eps);
+ d1d2_fd[i][j] = (distanceSquared<TargetSpace>(aPlus,bPlus) + distanceSquared<TargetSpace>(aMinus,bMinus)
+ - distanceSquared<TargetSpace>(aPlus,bMinus) - distanceSquared<TargetSpace>(aMinus,bPlus)) / (4*eps*eps);
}
}
- FieldMatrix<double,embeddedDim,embeddedDim> d1d2_diff = d1d2;
- d1d2_diff -= d1d2_fd;
-
- if ( (d1d2_diff).infinity_norm() > 200*eps ) {
+ if ( (d1d2 - d1d2_fd).infinity_norm() > 200*eps ) {
std::cout << className(a) << ": Analytical mixed second derivative does not match fd approximation." << std::endl;
std::cout << "d1d2 Analytical:" << std::endl << d1d2 << std::endl;
std::cout << "d1d2 FD :" << std::endl << d1d2_fd << std::endl;
@@ -220,7 +223,7 @@ void testMixedDerivativesOfSquaredDistance(const TargetSpace& a, const TargetSpa
template <class TargetSpace>
-void testDerivativeOfHessianOfSquaredDistance(const TargetSpace& a, const TargetSpace& b)
+void testDerivativeOfHessianOfDistanceSquared(const TargetSpace& a, const TargetSpace& b)
{
static const size_t embeddedDim = TargetSpace::embeddedDim;
@@ -259,7 +262,7 @@ void testDerivativeOfHessianOfSquaredDistance(const TargetSpace& a, const Target
template <class TargetSpace>
-void testMixedDerivativeOfHessianOfSquaredDistance(const TargetSpace& a, const TargetSpace& b)
+void testMixedDerivativeOfHessianOfDistanceSquared(const TargetSpace& a, const TargetSpace& b)
{
static const size_t embeddedDim = TargetSpace::embeddedDim;
@@ -298,38 +301,38 @@ void testMixedDerivativeOfHessianOfSquaredDistance(const TargetSpace& a, const T
template <class TargetSpace>
-void testDerivativesOfSquaredDistance(const TargetSpace& a, const TargetSpace& b)
+void testDerivativesOfDistanceSquared(const TargetSpace& a, const TargetSpace& b)
{
///////////////////////////////////////////////////////////////////
// Test derivative with respect to second argument
///////////////////////////////////////////////////////////////////
- testDerivativeOfSquaredDistance<TargetSpace>(a,b);
+ testDerivativeOfDistanceSquared<TargetSpace>(a,b);
///////////////////////////////////////////////////////////////////
// Test second derivative with respect to second argument
///////////////////////////////////////////////////////////////////
- testHessianOfSquaredDistance<TargetSpace>(a,b);
+ testHessianOfDistanceSquared<TargetSpace>(a,b);
//////////////////////////////////////////////////////////////////////////////
// Test mixed second derivative with respect to first and second argument
//////////////////////////////////////////////////////////////////////////////
- testMixedDerivativesOfSquaredDistance<TargetSpace>(a,b);
+ testMixedDerivativesOfDistanceSquared<TargetSpace>(a,b);
/////////////////////////////////////////////////////////////////////////////////////////////
// Test third derivative with respect to second argument
/////////////////////////////////////////////////////////////////////////////////////////////
- testDerivativeOfHessianOfSquaredDistance<TargetSpace>(a,b);
+ testDerivativeOfHessianOfDistanceSquared<TargetSpace>(a,b);
/////////////////////////////////////////////////////////////////////////////////////////////
// Test mixed third derivative with respect to first (once) and second (twice) argument
/////////////////////////////////////////////////////////////////////////////////////////////
- testMixedDerivativeOfHessianOfSquaredDistance<TargetSpace>(a,b);
+ testMixedDerivativeOfHessianOfDistanceSquared<TargetSpace>(a,b);
}
@@ -357,12 +360,11 @@ void test()
if (diameter(testPointPair) > TargetSpace::convexityRadius)
continue;
- TargetSpace p = testPointPair[0];
- TargetSpace q = testPointPair[1];
- std::cout << "p: " << testPointPair[0] << ", q: " << testPointPair[1] << std::endl;
- std::cout << TargetSpace::exp(p, TargetSpace::log(p,q)) << std::endl;
-
- testDerivativesOfSquaredDistance<TargetSpace>(testPoints[i], testPoints[j]);
+ // Test the exponential map and the logarithm
+ testExpLog(testPoints[i], testPoints[j]);
+
+ // Test the various derivatives of the squared distance
+ testDerivativesOfDistanceSquared<TargetSpace>(testPoints[i], testPoints[j]);
}
@@ -373,25 +375,30 @@ void test()
int main() try
{
-// test<RealTuple<double,1> >();
-// test<RealTuple<double,3> >();
+ // Test the RealTuple class
+ test<RealTuple<double,1> >();
+ test<RealTuple<double,3> >();
+ // Test the UnitVector class
test<UnitVector<double,2> >();
test<UnitVector<double,3> >();
test<UnitVector<double,4> >();
+ // Test the rotation class
+ test<Rotation<double,3> >();
+
+ // Test the ProductManifold class
test<Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2>>>();
test<Dune::GFE::ProductManifold<Rotation<double,3>,UnitVector<double,5>>>();
-// test<Rotation<double,3> >();
-//
-// test<RigidBodyMotion<double,3> >();
+ // Test the RigidBodyMotion class
+ test<RigidBodyMotion<double,3> >();
//
// test<HyperbolicHalfspacePoint<double,2> >();
} catch (Exception& e) {
- std::cout << e << std::endl;
+ std::cout << e.what() << std::endl;
}