diff --git a/test/averagedistanceassemblertest.cc b/test/averagedistanceassemblertest.cc
index 983b9707daf25fcac5cef1a211b83decdb4ea79c..03f97eac08b6ed3f0d99d94295ff3feb636123ea 100644
--- a/test/averagedistanceassemblertest.cc
+++ b/test/averagedistanceassemblertest.cc
@@ -9,6 +9,7 @@
#include <dune/gfe/rotation.hh>
#include <dune/gfe/realtuple.hh>
#include <dune/gfe/unitvector.hh>
+#include <dune/gfe/productmanifold.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
@@ -229,9 +230,29 @@ void testRotations()
}
+void testProductManifold()
+{
+ typedef Dune::GFE::ProductManifold<RealTuple<double,5>,UnitVector<double,3>, Rotation<double,3>> TargetSpace;
+
+ std::vector<TargetSpace> corners(dim+1);
+
+ std::generate(corners.begin(), corners.end(), []() {
+ return Dune::GFE::randomFieldVector<typename TargetSpace::field_type,TargetSpace::CoordinateType::dimension>(0.9,1.1);
+ });
+
+ TargetSpace argument = corners[0];
+ testWeightSet(corners, argument);
+ argument = corners[1];
+ testWeightSet(corners, argument);
+ argument = corners[2];
+ testWeightSet(corners, argument);
+}
+
+
int main()
{
testRealTuples();
testUnitVectors();
testRotations();
+ testProductManifold();
}
diff --git a/test/localgeodesicfefunctiontest.cc b/test/localgeodesicfefunctiontest.cc
index a2fa30ba11d9fd40d18180f39a36141ee3fef364..613fbda902bda20a398d7e7b313e40990108f89a 100644
--- a/test/localgeodesicfefunctiontest.cc
+++ b/test/localgeodesicfefunctiontest.cc
@@ -14,6 +14,7 @@
#include <dune/gfe/rotation.hh>
#include <dune/gfe/realtuple.hh>
#include <dune/gfe/unitvector.hh>
+#include <dune/gfe/productmanifold.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
#include "multiindex.hh"
@@ -303,6 +304,8 @@ int main()
test<UnitVector<double,3>,1>(GeometryTypes::simplex(1));
test<Rotation<double,3>,1>(GeometryTypes::simplex(1));
test<RigidBodyMotion<double,3>,1>(GeometryTypes::simplex(1));
+ typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2>> CrazyManifold;
+ test<CrazyManifold,1>(GeometryTypes::simplex(1));
////////////////////////////////////////////////////////////////
// Test functions on 2d simplex elements
@@ -313,6 +316,8 @@ int main()
test<UnitVector<double,3>,2>(GeometryTypes::simplex(2));
test<Rotation<double,3>,2>(GeometryTypes::simplex(2));
test<RigidBodyMotion<double,3>,2>(GeometryTypes::simplex(2));
+ typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2>> CrazyManifold;
+ test<CrazyManifold,2>(GeometryTypes::simplex(2));
////////////////////////////////////////////////////////////////
// Test functions on 2d quadrilateral elements
@@ -323,5 +328,7 @@ int main()
test<UnitVector<double,3>,2>(GeometryTypes::cube(2));
test<Rotation<double,3>,2>(GeometryTypes::cube(2));
test<RigidBodyMotion<double,3>,2>(GeometryTypes::cube(2));
+ typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2>> CrazyManifold;
+ test<CrazyManifold,2>(GeometryTypes::cube(2));
}
diff --git a/test/localgfetestfunctiontest.cc b/test/localgfetestfunctiontest.cc
index 4b9769c95f7348e66363ed16f9eaa88d8ec5ae38..2890f50e383039173b7d4581a44368bd93525866 100644
--- a/test/localgfetestfunctiontest.cc
+++ b/test/localgfetestfunctiontest.cc
@@ -10,6 +10,7 @@
#include <dune/gfe/rotation.hh>
#include <dune/gfe/realtuple.hh>
#include <dune/gfe/unitvector.hh>
+#include <dune/gfe/productmanifold.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/localgfetestfunctionbasis.hh>
@@ -79,6 +80,8 @@ int main() try
test<Rotation<double,3>, 1>();
test<Rotation<double,3>, 2>();
+ typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2>> CrazyManifold;
+ test<CrazyManifold, 2>();
} catch (Exception& e) {
std::cout << e << std::endl;
diff --git a/test/localprojectedfefunctiontest.cc b/test/localprojectedfefunctiontest.cc
index d2f6829398d13b3c60ba4052988e38ad8f701a7d..ec6164f4f3315f536ba722654b37de9c6338e1dd 100644
--- a/test/localprojectedfefunctiontest.cc
+++ b/test/localprojectedfefunctiontest.cc
@@ -15,6 +15,7 @@
#include <dune/gfe/rotation.hh>
#include <dune/gfe/realtuple.hh>
#include <dune/gfe/unitvector.hh>
+#include <dune/gfe/productmanifold.hh>
#include <dune/gfe/localprojectedfefunction.hh>
#include "multiindex.hh"
@@ -105,6 +106,20 @@ void testDerivativeTangentiality(const RigidBodyMotion<double,3>& x,
{
}
+// the columns of the derivative must be tangential to the manifold
+template <int domainDim, int vectorDim,typename ...TargetSpaces>
+void testDerivativeTangentiality(const Dune::GFE::ProductManifold<TargetSpaces...>& x,
+ const FieldMatrix<double,vectorDim,domainDim>& derivative)
+{
+ size_t posHelper=0;
+ using namespace Dune::Hybrid;
+ forEach(integralRange(Dune::Hybrid::size(x)), [&](auto&& i) {
+ using Manifold = std::remove_reference_t<decltype(x[i])>;
+ testDerivativeTangentiality(x[i],Dune::GFE::blockAt<Manifold::embeddedDim,domainDim>( derivative,posHelper,0));
+ posHelper +=Manifold::embeddedDim;
+ });
+}
+
/** \brief Test whether interpolation is invariant under permutation of the simplex vertices
* \todo Implement this for all dimensions
*/
@@ -255,6 +270,8 @@ int main()
test<UnitVector<double,3>,1>(GeometryTypes::line);
test<Rotation<double,3>,1>(GeometryTypes::line);
test<RigidBodyMotion<double,3>,1>(GeometryTypes::line);
+ typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2>> CrazyManifold;
+ test<CrazyManifold, 1>(GeometryTypes::line);
////////////////////////////////////////////////////////////////
// Test functions on 2d simplex elements
@@ -265,6 +282,8 @@ int main()
test<UnitVector<double,3>,2>(GeometryTypes::triangle);
test<Rotation<double,3>,2>(GeometryTypes::triangle);
test<RigidBodyMotion<double,3>,2>(GeometryTypes::triangle);
+ typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2>> CrazyManifold;
+ test<CrazyManifold, 2>(GeometryTypes::triangle);
////////////////////////////////////////////////////////////////
// Test functions on 2d quadrilateral elements
@@ -275,5 +294,7 @@ int main()
test<UnitVector<double,3>,2>(GeometryTypes::quadrilateral);
test<Rotation<double,3>,2>(GeometryTypes::quadrilateral);
test<RigidBodyMotion<double,3>,2>(GeometryTypes::quadrilateral);
+ typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2>> CrazyManifold;
+ test<CrazyManifold, 2>(GeometryTypes::quadrilateral);
}
diff --git a/test/targetspacetest.cc b/test/targetspacetest.cc
index cadfe4818e3f3c0c8bfce2ff69c7b0e22acaa9cb..b5e1c650bf1e91361fb1bf8fd7111ac40414bff8 100644
--- a/test/targetspacetest.cc
+++ b/test/targetspacetest.cc
@@ -3,6 +3,7 @@
#include <dune/gfe/unitvector.hh>
#include <dune/gfe/realtuple.hh>
#include <dune/gfe/rotation.hh>
+#include <dune/gfe/productmanifold.hh>
#include <dune/gfe/hyperbolichalfspacepoint.hh>
#include "valuefactory.hh"
@@ -138,8 +139,8 @@ void testDerivativeOfSquaredDistance(const TargetSpace& a, const TargetSpace& b)
// transform into embedded coordinates
typename TargetSpace::EmbeddedTangentVector d2_fd_embedded;
B.mtv(d2_fd,d2_fd_embedded);
-
- if ( (d2 - d2_fd_embedded).infinity_norm() > 100*eps ) {
+
+ if ( (d2 - d2_fd_embedded).infinity_norm() > 200*eps ) {
std::cout << className(a) << ": Analytical gradient does not match fd approximation." << std::endl;
std::cout << "d2 Analytical: " << d2 << std::endl;
std::cout << "d2 FD : " << d2_fd << std::endl;
@@ -156,14 +157,14 @@ void testHessianOfSquaredDistance(const TargetSpace& a, const TargetSpace& b)
///////////////////////////////////////////////////////////////////
// Test second derivative with respect to second argument
///////////////////////////////////////////////////////////////////
- FieldMatrix<double,embeddedDim,embeddedDim> d2d2 = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(a, b);
-
+ FieldMatrix<double,embeddedDim,embeddedDim> d2d2 = (TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(a, b)).matrix();
+
// 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() > 100*eps) {
+ if ( (d2d2_diff).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;
@@ -207,7 +208,8 @@ void testMixedDerivativesOfSquaredDistance(const TargetSpace& a, const TargetSpa
FieldMatrix<double,embeddedDim,embeddedDim> d1d2_diff = d1d2;
d1d2_diff -= d1d2_fd;
- if ( (d1d2_diff).infinity_norm() > 100*eps ) {
+
+ if ( (d1d2_diff).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;
@@ -245,8 +247,8 @@ void testDerivativeOfHessianOfSquaredDistance(const TargetSpace& a, const Target
d2d2d2_fd[i] /= 2*eps;
}
-
- if ( (d2d2d2 - d2d2d2_fd).infinity_norm() > 100*eps) {
+
+ if ( (d2d2d2 - d2d2d2_fd).infinity_norm() > 200*eps) {
std::cout << className(a) << ": Analytical third derivative does not match fd approximation." << std::endl;
std::cout << "d2d2d2 Analytical:" << std::endl << d2d2d2 << std::endl;
std::cout << "d2d2d2 FD :" << std::endl << d2d2d2_fd << std::endl;
@@ -284,8 +286,8 @@ void testMixedDerivativeOfHessianOfSquaredDistance(const TargetSpace& a, const T
d1d2d2_fd[i] /= 2*eps;
}
-
- if ( (d1d2d2 - d1d2d2_fd).infinity_norm() > 100*eps ) {
+
+ if ( (d1d2d2 - d1d2d2_fd).infinity_norm() > 200*eps ) {
std::cout << className(a) << ": Analytical mixed third derivative does not match fd approximation." << std::endl;
std::cout << "d1d2d2 Analytical:" << std::endl << d1d2d2 << std::endl;
std::cout << "d1d2d2 FD :" << std::endl << d1d2d2_fd << std::endl;
@@ -360,7 +362,7 @@ void test()
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]);
+ testDerivativesOfSquaredDistance<TargetSpace>(testPoints[i], testPoints[j]);
}
@@ -378,6 +380,9 @@ int main() try
test<UnitVector<double,3> >();
test<UnitVector<double,4> >();
+ 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> >();
diff --git a/test/valuefactory.hh b/test/valuefactory.hh
index 5a1a795696a85a38d6edaf2d3df1968416dcebb1..9d23579177b486e56e17bb168efda579c454c7bf 100644
--- a/test/valuefactory.hh
+++ b/test/valuefactory.hh
@@ -6,6 +6,7 @@
#include <dune/gfe/unitvector.hh>
#include <dune/gfe/rotation.hh>
#include <dune/gfe/rigidbodymotion.hh>
+#include <dune/gfe/productmanifold.hh>
#include <dune/gfe/orthogonalmatrix.hh>
#include <dune/gfe/hyperbolichalfspacepoint.hh>
@@ -328,4 +329,30 @@ public:
+
+
+/** \brief A class that creates sets of values of various types, to be used in unit tests
+*
+* This is the specialization for ProducManifold<...>
+*/
+template <typename ...TargetSpaces>
+class ValueFactory<Dune::GFE::ProductManifold<TargetSpaces...> >
+{
+using TargetSpace = Dune::GFE::ProductManifold<TargetSpaces...>;
+
+public:
+ static void get(std::vector<TargetSpace >& values) {
+
+ std::vector<typename TargetSpace::CoordinateType > testPoints(10);
+
+ std::generate(testPoints.begin(), testPoints.end(), [](){
+ return Dune::GFE::randomFieldVector<typename TargetSpace::field_type,TargetSpace::CoordinateType::dimension>(0.9,1.1) ; });
+
+ values.resize(testPoints.size());
+
+ std::transform(testPoints.cbegin(),testPoints.cend(),values.begin(),[](const auto& testPoint){return TargetSpace(testPoint);});
+ }
+};
+
+
#endif
\ No newline at end of file