diff --git a/dune/gfe/linearalgebra.hh b/dune/gfe/linearalgebra.hh
index 1a08a863cbfb905f8dd8f2edd568bc7ea5485edb..57a1c66ae8016329459649e28c2781f39711615b 100644
--- a/dune/gfe/linearalgebra.hh
+++ b/dune/gfe/linearalgebra.hh
@@ -1,6 +1,8 @@
#ifndef DUNE_GFE_LINEARALGEBRA_HH
#define DUNE_GFE_LINEARALGEBRA_HH
+#include <random>
+
#include <dune/common/fmatrix.hh>
#include <dune/common/version.hh>
#include <dune/istl/scaledidmatrix.hh>
@@ -134,6 +136,66 @@ namespace Dune {
}
#endif
+ /** \brief Return a segment of a FieldVector from lower up to lower+size-1 */
+ template< int lower, int size,typename field_type,int n>
+ static FieldVector<field_type,size> segment(const FieldVector<field_type,n>& v)
+ {
+ FieldVector<field_type,size> res;
+ std::copy(v.begin()+lower,v.begin()+lower+size,res.begin());
+ return res;
+ }
+
+ /** \brief Return a segment of a FieldVector from lower up to lower+size-1
+ * lower is unkown at compile time*/
+ template< int size,typename field_type,int n>
+ static FieldVector<field_type,size> segmentAt(const FieldVector<field_type,n>& v,const size_t lower)
+ {
+ FieldVector<field_type,size> res;
+ std::copy(v.begin()+lower,v.begin()+lower+size,res.begin());
+ return res;
+ }
+
+ /** \brief Return a block of a FieldMatrix (lower1...lower1+size1-1,lower2...lower2+size2-1 */
+ template< int lower1, int size1, int lower2,int size2,typename field_type,int n,int m>
+ static auto block(const FieldMatrix<field_type,n,m>& v)
+ {
+ static_assert(lower1+size1<=n && lower2+size2<=m, "Size mismatch for Block!");
+ FieldMatrix<field_type,size1,size2> res;
+
+ for(int i=lower1; i<lower1+size1; ++i)
+ for(int j=lower2; j<lower2+size2; ++j)
+ res[i-lower1][j-lower2] = v[i][j];
+ return res;
+ }
+
+ /** \brief Return a block of a FieldMatrix (lower1...lower1+size1-1,lower2...lower2+size2-1
+ * * lower1 and lower2 is unkown at compile time*/
+ template< int size1,int size2,typename field_type,int n,int m>
+ static auto blockAt(const FieldMatrix<field_type,n,m>& v, const size_t& lower1, const size_t& lower2)
+ {
+ assert(lower1+size1<=n && lower2+size2<=m);
+ FieldMatrix<field_type,size1,size2> res;
+
+ for(size_t i=lower1; i<lower1+size1; ++i)
+ for(size_t j=lower2; j<lower2+size2; ++j)
+ res[i-lower1][j-lower2] = v[i][j];
+ return res;
+ }
+
+ /** \brief Generates FieldVector with random entries in the range -1..1 */
+ template<typename field_type,int n>
+ auto randomFieldVector(field_type lower=-1, field_type upper=1)
+ {
+ std::random_device rd;
+ std::mt19937 mt(rd());
+ std::uniform_real_distribution<field_type> dist(lower, upper);
+ auto rand = [&dist,&mt](){
+ return dist(mt);
+ };
+ FieldVector<field_type,n> vec;
+ std::generate(vec.begin(), vec.end(), rand);
+ return vec;
+ }
}
}
diff --git a/dune/gfe/productmanifold.hh b/dune/gfe/productmanifold.hh
new file mode 100644
index 0000000000000000000000000000000000000000..48e86df8a9d44ce355f9ab6715a91b9c8b5f8d3f
--- /dev/null
+++ b/dune/gfe/productmanifold.hh
@@ -0,0 +1,503 @@
+#ifndef DUNE_GFE_PRODUCTMANIFOLD_HH
+#define DUNE_GFE_PRODUCTMANIFOLD_HH
+
+#include <iostream>
+#include <tuple>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/math.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/tuplevector.hh>
+#include <dune/common/power.hh>
+
+#include <dune/gfe/linearalgebra.hh>
+
+namespace Dune::GFE
+{
+ namespace Impl
+ {
+ template<typename T, typename ... Ts>
+ constexpr auto sumDim()
+ {
+ return (T::dim + ... + sumDim<Ts>());
+ }
+
+ template<typename T, typename ... Ts>
+ constexpr auto sumEmbeddedDim()
+ {
+ return (T::embeddedDim + ... + sumEmbeddedDim<Ts>());
+ }
+
+ template<typename T, typename ... Ts>
+ constexpr auto variadicConvexityRadiusMin()
+ {
+ if constexpr (sizeof...(Ts)!=0)
+ return std::min(T::convexityRadius , variadicConvexityRadiusMin<Ts ...>());
+ else
+ return T::convexityRadius;
+ }
+
+ template <typename TS, typename ... TargetSpaces>
+ class ProductManifold;
+
+ template<class U,typename Tfirst,typename ... TargetSpaces2>
+ struct rebindHelper
+ {
+ typedef ProductManifold<typename Tfirst::template rebind<U>::other ,typename rebindHelper<U,TargetSpaces2...>::other> other;
+ };
+
+ template<class U,typename Tlast>
+ struct rebindHelper<U,Tlast>
+ {
+ typedef typename Tlast::template rebind<U>::other other;
+ };
+ }
+
+ /** \brief A Product manifold */
+ template <typename TS, typename ... TargetSpaces>
+ class ProductManifold
+ {
+ public:
+ template<std::size_t I>
+ using IC = Dune::index_constant<I>;
+ /** \brief The type used for coordinates. */
+ using ctype = typename std::common_type<typename TS::ctype, typename TargetSpaces::ctype...>::type ;
+ using field_type = typename std::common_type<typename TS::ctype, typename TargetSpaces::ctype...>::type ;
+
+ /** \brief Number of factors */
+ static constexpr int numTS = 1 + sizeof...(TargetSpaces);
+
+ /** \brief Dimension of manifold */
+ static constexpr int dim = TS::dim + Impl::sumDim<TargetSpaces ...>();
+
+ /** \brief Dimension of the embedding space */
+ static constexpr int embeddedDim = TS::embeddedDim + Impl::sumEmbeddedDim<TargetSpaces ...>();
+
+ /** \brief Type of a tangent of the ProductManifold with inner dimensions*/
+ typedef Dune::FieldVector<field_type, dim> TangentVector;
+
+ /** \brief Type of a tangent of the ProductManifold represented in the embedding space*/
+ typedef Dune::FieldVector<field_type, embeddedDim> EmbeddedTangentVector;
+
+ /** \brief The global convexity radius of the Product space */
+ static constexpr double convexityRadius = Impl::variadicConvexityRadiusMin<TS,TargetSpaces ...>();
+
+ /** \brief The type used for global coordinates */
+ typedef Dune::FieldVector<field_type ,embeddedDim> CoordinateType;
+
+ /** \brief Default constructor */
+ ProductManifold() = default;
+
+ /** \brief Constructor from another ProductManifold */
+ ProductManifold(const ProductManifold<TS,TargetSpaces ...>& productManifold)
+ : data_(productManifold.data_)
+ {}
+
+ /** \brief Constructor from a coordinates vector of the embedding space */
+ explicit ProductManifold(const CoordinateType& globalCoordinates)
+ {
+ DUNE_ASSERT_BOUNDS(globalCoordinates.size()== sumEmbeddedDim)
+ auto constructorFunctor =[] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& res = std::get<0>(argsTuple)[manifoldInt];
+ const auto& globalCoords =std::get<1>(argsTuple);
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(res)> >;
+ res = Manifold(Dune::GFE::segmentAt<Manifold::embeddedDim>(globalCoords,posHelper[0]));
+ posHelper[0] += Manifold::embeddedDim;
+ };
+ foreachManifold(constructorFunctor,*this,globalCoordinates);
+ }
+
+ /** \brief Assignment from a coordinates vector of the embedding space */
+ ProductManifold<TS,TargetSpaces ...>& operator=(const CoordinateType& globalCoordinates)
+ {
+ ProductManifold<TS,TargetSpaces ...> res(globalCoordinates);
+ data_ = res.data_;
+ return *this;
+ }
+
+ /** \brief Assignment from ProductManifold with different type -- used for automatic differentiation with ADOL-C */
+ template <typename TS2, typename ... TargetSpaces2>
+ ProductManifold<TS,TargetSpaces ...>& operator <<=(const ProductManifold<TS2,TargetSpaces2 ...>& other)
+ {
+ forEach(integralRange(IC<numTS>()), [&](auto&& i) {
+ (*this)[i] <<= other[i];
+ });
+ return *this;
+ }
+
+ /** \brief Rebind the ProductManifold to another coordinate type using an embedded tangent vector*/
+ template<class U,typename ... TargetSpaces2>
+ struct rebind
+ {
+ typedef typename Impl::rebindHelper<U,TS,TargetSpaces...>::other other;
+ };
+
+ /** \brief The exponential map from a given point. */
+ static ProductManifold<TS,TargetSpaces ...> exp(const ProductManifold<TS,TargetSpaces ...>& p, const EmbeddedTangentVector& v)
+ {
+ ProductManifold<TS,TargetSpaces ...> res;
+ auto expFunctor =[] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& res = std::get<0>(argsTuple)[manifoldInt];
+ const auto& p = std::get<1>(argsTuple)[manifoldInt];
+ const auto& tang = std::get<2>(argsTuple);
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(res)> >;
+ auto currentEmbeddedTangentVector = Dune::GFE::segmentAt<Manifold::embeddedDim>(tang,posHelper[0]);
+ res = Manifold::exp(p,currentEmbeddedTangentVector);
+ posHelper[0] += Manifold::embeddedDim;
+ };
+ foreachManifold(expFunctor,res,p,v);
+ return res;
+ }
+
+ /** \brief The exponential map from a given point using an intrinsic tangent vector*/
+ static ProductManifold<TS,TargetSpaces ...> exp(const ProductManifold<TS,TargetSpaces ...>& p, const TangentVector& v)
+ {
+ auto basis = p.orthonormalFrame();
+ EmbeddedTangentVector embeddedTangent;
+ basis.mtv(v, embeddedTangent);
+ return exp(p,embeddedTangent);
+ }
+
+ /** \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)
+ {
+ TangentVector diff;
+ auto logFunctor =[] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& res = std::get<0>(argsTuple);
+ const auto& a = std::get<1>(argsTuple)[manifoldInt];
+ const auto& b = std::get<2>(argsTuple)[manifoldInt];
+ 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;
+ };
+ foreachManifold(logFunctor,diff,a, b);
+ return diff;
+ }
+
+ /** \brief Compute geodesic distance from a to b */
+ static field_type distance(const ProductManifold<TS,TargetSpaces...>& a, const ProductManifold<TS,TargetSpaces...>& b)
+ {
+ field_type dist=0.0;
+ auto distanceFunctor =[] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& res = std::get<0>(argsTuple);
+ const auto& a = std::get<1>(argsTuple)[manifoldInt];
+ const auto& b = std::get<2>(argsTuple)[manifoldInt];
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(a)> >;
+ res += Dune::Power<2>().eval(Manifold::distance(a,b));
+ };
+ foreachManifold(distanceFunctor,dist,a, b);
+ return sqrt(dist);
+ }
+
+ /** \brief Compute the gradient of the squared distance function keeping the first argument fixed */
+ static EmbeddedTangentVector derivativeOfDistanceSquaredWRTSecondArgument(const ProductManifold<TS,TargetSpaces...>& a,
+ const ProductManifold<TS,TargetSpaces...>& b)
+ {
+ EmbeddedTangentVector derivative;
+ derivative= 0.0;
+
+ auto derivOfDistSqdWRTSecArgFunctor = [] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& derivative = std::get<0>(argsTuple);
+ const auto& a = std::get<1>(argsTuple)[manifoldInt];
+ const auto& b = std::get<2>(argsTuple)[manifoldInt];
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(a)> >;
+ const auto diffLoc = Manifold::derivativeOfDistanceSquaredWRTSecondArgument(a,b);
+ std::copy(diffLoc.begin(),diffLoc.end(),derivative.begin()+posHelper[0]);
+ posHelper[0] += Manifold::embeddedDim;
+ };
+ foreachManifold(derivOfDistSqdWRTSecArgFunctor, derivative,a, b);
+ return derivative;
+ }
+
+ /** \brief Compute the Hessian of the squared distance function keeping the first argument fixed */
+ static auto secondDerivativeOfDistanceSquaredWRTSecondArgument(const ProductManifold<TS,TargetSpaces...>& a,
+ const ProductManifold<TS, TargetSpaces...>& b)
+ {
+ Dune::SymmetricMatrix<field_type,embeddedDim> result;
+ auto secDerivOfDistSqWRTSecArgFunctor = [] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& deriv = std::get<0>(argsTuple);
+ const auto& a = std::get<1>(argsTuple)[manifoldInt];
+ const auto& b = std::get<2>(argsTuple)[manifoldInt];
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(a)> >;
+ const auto diffLoc = Manifold::secondDerivativeOfDistanceSquaredWRTSecondArgument(a,b);
+ for(size_t i=posHelper[0]; i<posHelper[0]+Manifold::embeddedDim; ++i )
+ for(size_t j=posHelper[0]; j<=i; ++j )
+ deriv(i,j) = diffLoc(i - posHelper[0], j - posHelper[0]);
+ posHelper[0] += Manifold::embeddedDim;
+ };
+ foreachManifold(secDerivOfDistSqWRTSecArgFunctor,result,a, b);
+ return result;
+ }
+
+ /** \brief Compute the mixed second derivate \partial d^2 / \partial da db */
+ static Dune::FieldMatrix<field_type ,embeddedDim,embeddedDim>
+ secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(const ProductManifold<TS,TargetSpaces ...>& a,
+ const ProductManifold<TS, TargetSpaces ...>& b)
+ {
+ Dune::FieldMatrix<field_type,embeddedDim,embeddedDim> result(0);
+ auto secDerivOfDistSqWRTFirstAndSecArgFunctor = [] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& deriv = std::get<0>(argsTuple);
+ const auto& a = std::get<1>(argsTuple)[manifoldInt];
+ const auto& b = std::get<2>(argsTuple)[manifoldInt];
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(a)> >;
+ const auto diffLoc = Manifold::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(a,b);
+ for(size_t i=posHelper[0]; i<posHelper[0]+Manifold::embeddedDim; ++i )
+ for(size_t j=posHelper[0]; j<posHelper[0]+Manifold::embeddedDim; ++j )
+ deriv[i][j] = diffLoc[i - posHelper[0]][ j - posHelper[0]];
+ posHelper[0] += Manifold::embeddedDim;
+ };
+ foreachManifold(secDerivOfDistSqWRTFirstAndSecArgFunctor,result,a, b);
+ return result;
+ }
+
+ /** \brief Compute the third derivative \partial d^3 / \partial dq^3 */
+ static Tensor3<field_type ,embeddedDim,embeddedDim,embeddedDim>
+ thirdDerivativeOfDistanceSquaredWRTSecondArgument(const ProductManifold<TS,TargetSpaces ...>& a,
+ const ProductManifold<TS, TargetSpaces ...>& b)
+ {
+ Tensor3<field_type,embeddedDim,embeddedDim,embeddedDim> result(0);
+ auto thirdDerivOfDistSqWRTSecArgFunctor =[](auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& deriv = std::get<0>(argsTuple);
+ const auto& a = std::get<1>(argsTuple)[manifoldInt];
+ const auto& b = std::get<2>(argsTuple)[manifoldInt];
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(a)>>;
+ const auto diffLoc = Manifold::thirdDerivativeOfDistanceSquaredWRTSecondArgument(a,b);
+ for(size_t i=posHelper[0]; i<posHelper[0]+Manifold::embeddedDim; ++i )
+ for(size_t j=posHelper[0]; j<posHelper[0]+Manifold::embeddedDim; ++j )
+ for(size_t k=posHelper[0]; k<posHelper[0]+Manifold::embeddedDim; ++k )
+ deriv[i][j][k] = diffLoc[i - posHelper[0]][ j - posHelper[0]][k - posHelper[0]];
+ posHelper[0] += Manifold::embeddedDim;
+ };
+ foreachManifold(thirdDerivOfDistSqWRTSecArgFunctor,result,a, b);
+ return result;
+ }
+
+ /** \brief Compute the mixed third derivative \partial d^3 / \partial da db^2 */
+ static Tensor3<field_type ,embeddedDim,embeddedDim,embeddedDim>
+ thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(const ProductManifold<TS,TargetSpaces ...>& a,
+ const ProductManifold<TS, TargetSpaces ...>& b)
+ {
+ Tensor3<field_type,embeddedDim,embeddedDim,embeddedDim> result(0);
+ auto thirdDerivOfDistSqWRT1And2ArgFunctor =[] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& deriv = std::get<0>(argsTuple);
+ const auto& a = std::get<1>(argsTuple)[manifoldInt];
+ const auto& b = std::get<2>(argsTuple)[manifoldInt];
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(a)>>;
+ const auto diffLoc = Manifold::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(a,b);
+ for(size_t i=posHelper[0]; i<posHelper[0]+Manifold::embeddedDim; ++i )
+ for(size_t j=posHelper[0]; j<posHelper[0]+Manifold::embeddedDim; ++j )
+ for(size_t k=posHelper[0]; k<posHelper[0]+Manifold::embeddedDim; ++k )
+ deriv[i][j][k] = diffLoc[i - posHelper[0]][ j - posHelper[0]][k - posHelper[0]];
+ posHelper[0] += Manifold::embeddedDim;
+ };
+ foreachManifold(thirdDerivOfDistSqWRT1And2ArgFunctor,result,a, b);
+ return result;
+ }
+
+ /** \brief Project tangent vector of R^n onto the tangent space */
+ EmbeddedTangentVector projectOntoTangentSpace(const EmbeddedTangentVector& v) const
+ {
+ EmbeddedTangentVector result {v};
+ auto projectOntoTangentSpaceFunctor =[] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& v = std::get<0>(argsTuple);
+ const auto& a = std::get<1>(argsTuple)[manifoldInt];
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(a)> >;
+ const auto vLoc = a.projectOntoTangentSpace(Dune::GFE::segmentAt<Manifold::embeddedDim>(v,posHelper[0]));
+ std::copy(vLoc.begin(),vLoc.end(),v.begin()+posHelper[0]);
+ posHelper[0] += Manifold::embeddedDim;
+ };
+ foreachManifold(projectOntoTangentSpaceFunctor,result,*this);
+ return result;
+ }
+
+ /** \brief Project tangent vector of R^n onto the normal space space */
+ EmbeddedTangentVector projectOntoNormalSpace(const EmbeddedTangentVector& v) const
+ {
+ EmbeddedTangentVector result {v};
+ auto projectOntoNormalSpaceFunctor =[] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& v = std::get<0>(argsTuple);
+ const auto& a = std::get<1>(argsTuple)[manifoldInt];
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(a)> >;
+ const auto vLoc = a.projectOntoNormalSpace(Dune::GFE::segmentAt<Manifold::embeddedDim>(v,posHelper[0]));
+ std::copy(vLoc.begin(),vLoc.end(),v.begin()+posHelper[0]);
+ posHelper[0] += Manifold::embeddedDim;
+ };
+ foreachManifold(projectOntoNormalSpaceFunctor,result,*this);
+ return result;
+ }
+
+ /** \brief Project tangent vector of R^n onto the normal space space */
+ static ProductManifold<TS,TargetSpaces ...> projectOnto(const CoordinateType& v)
+ {
+ ProductManifold<TS,TargetSpaces ...> result {v};
+ auto projectOntoFunctor =[] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& res = std::get<0>(argsTuple)[manifoldInt];
+ auto& v = std::get<1>(argsTuple);
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(res)> >;
+ res = Manifold::projectOnto(Dune::GFE::segmentAt<Manifold::embeddedDim>(v,posHelper[0]));
+ posHelper[0] += Manifold::embeddedDim;
+ };
+ foreachManifold(projectOntoFunctor,result, v);
+ return result;
+ }
+
+ /** \brief Project tangent vector of R^n onto the normal space space */
+ static auto derivativeOfProjection(const CoordinateType& v)
+ {
+ Dune::FieldMatrix<typename CoordinateType::value_type, embeddedDim, embeddedDim> result;
+ auto derivativeOfProjectionFunctor =[] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& res = std::get<0>(argsTuple);
+ const auto& v = std::get<1>(argsTuple);
+ const Dune::TupleVector<TS,TargetSpaces...> ManifoldTuple;
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(ManifoldTuple[manifoldInt])> >;
+ const auto vLoc = Manifold::derivativeOfProjection(Dune::GFE::segmentAt<Manifold::embeddedDim>(v,posHelper[0]));
+ for(size_t k=posHelper[0]; k<posHelper[0]+Manifold::embeddedDim; ++k )
+ for(size_t j=posHelper[0]; j<posHelper[0]+Manifold::embeddedDim; ++j )
+ res[k][j] = vLoc[k - posHelper[0]][j-posHelper[0]];
+ posHelper[0] += Manifold::embeddedDim;
+ };
+ foreachManifold(derivativeOfProjectionFunctor,result, v);
+ return result;
+ }
+
+ /** \brief The Weingarten map */
+ EmbeddedTangentVector weingarten(const EmbeddedTangentVector& z, const EmbeddedTangentVector& v) const
+ {
+ EmbeddedTangentVector result;
+ auto weingartenFunctor =[] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& res = std::get<0>(argsTuple);
+ const auto& a = std::get<1>(argsTuple)[manifoldInt];
+ const auto& z = std::get<2>(argsTuple);
+ const auto& v = std::get<3>(argsTuple);
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(a)> >;
+ const auto resLoc = a.weingarten(Dune::GFE::segmentAt<Manifold::embeddedDim>(z,posHelper[0]),
+ Dune::GFE::segmentAt<Manifold::embeddedDim>(v,posHelper[0]));
+ std::copy(resLoc.begin(),resLoc.end(),res.begin()+posHelper[0]);
+ posHelper[0] +=Manifold::embeddedDim;
+ };
+ foreachManifold(weingartenFunctor,result, *this,z, v);
+ return result;
+ }
+
+ /** \brief Compute an orthonormal basis of the tangent space of the ProductManifold
+ This basis may not be globally continuous. */
+ Dune::FieldMatrix<field_type ,dim,embeddedDim> orthonormalFrame() const
+ {
+ Dune::FieldMatrix<field_type,dim,embeddedDim> result(0);
+ auto orthonormalFrameFunctor =[] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& manifoldInt)
+ {
+ auto& res = std::get<0>(argsTuple);
+ const auto& a = std::get<1>(argsTuple)[manifoldInt];
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(a)> >;
+ const auto resLoc = a.orthonormalFrame();
+ for(size_t i=posHelper[0]; i<posHelper[0]+Manifold::dim; ++i )
+ for(size_t j=posHelper[1]; j<posHelper[1]+Manifold::embeddedDim; ++j )
+ res[i][j] = resLoc[i - posHelper[0]][j-posHelper[1]];
+ posHelper[0] += Manifold::dim;
+ posHelper[1] += Manifold::embeddedDim;
+ };
+ foreachManifold(orthonormalFrameFunctor,result,*this);
+ return result;
+ }
+
+ /** \brief The global coordinates */
+ CoordinateType globalCoordinates() const
+ {
+ CoordinateType returnValue;
+ auto globalCoordinatesFunctor =[] (auto& argsTuple, std::array<std::size_t,2>& posHelper, const auto& i)
+ {
+ auto& res = std::get<0>(argsTuple);
+ const auto& p = std::get<1>(argsTuple)[i];
+ const auto vLoc = p.globalCoordinates();
+ using Manifold = std::remove_const_t<typename std::remove_reference_t<decltype(p)> >;
+ std::copy(vLoc.begin(),vLoc.end(),res.begin()+posHelper[0]);
+ posHelper[0] += Manifold::embeddedDim;
+ };
+ foreachManifold(globalCoordinatesFunctor,returnValue,*this);
+ return returnValue;
+ }
+
+ /** \brief Const access to the tuple elements */
+ template<std::size_t i>
+ constexpr decltype(auto) operator[] (const Dune::index_constant<i>&) const
+ {
+ return std::get<i>(data_);
+ }
+
+ /** \brief Non-const access to the tuple elements */
+ template<std::size_t i>
+ decltype(auto) operator[] (const Dune::index_constant<i>&)
+ {
+ return std::get<i>(data_);
+ }
+
+ /** \brief Number of elements of the tuple */
+ static constexpr std::size_t size()
+ {
+ return numTS;
+ }
+
+ template<class TS2, class... TargetSpaces2>
+ friend std::ostream& operator<<(std::ostream& s, const ProductManifold<TS2, TargetSpaces2 ...>& c);
+
+ private:
+ /**
+ * \brief Range based for loop over all Manifolds in ProductManifold
+ *
+ * \tparam FunctorType the functor which should be applied for all manifolds
+ * \tparam Args... The arguments which are passed to the functor
+ *
+ * The functor has to extract the current manifold from its arguments.
+ * Therefore, the current index i is passed to the functor.
+ *
+ * Furthermore, the functor gets two integers to deal with FieldVectors
+ * or FieldMatrices which span the whole (embedding) coordinate range of the ProductManifold.
+ *
+ * Example:
+ * If a argument of the functor is
+ * ProductManifold<Realtuple<double,3>,UnitVector<double,3>>::EmbeddedTangentVector,
+ * the embedded coordinate vector has 6 entries. The functor needs to know where
+ * to insert the next subvector of the submanifold. Therefore, posHelper[0] should be 0 at the first
+ * iteration and 3 at the second one. Then first indices [0..2] stores the result of
+ * Realtuple<double,3> and [3..5] stores
+ * the result of UnitVector<double,3>.
+ * See e.g. the globalCoordinates() function
+ */
+ template<typename FunctorType, typename ... Args>
+ static void foreachManifold( FunctorType&& functor,Args&&... args)
+ {
+ auto argsTuple = std::forward_as_tuple(args ...);
+ std::array<std::size_t,2> posHelper({0,0});
+ Dune::Hybrid::forEach(Dune::Hybrid::integralRange(IC<numTS>()),[&](auto&& i) {
+ functor(argsTuple,posHelper,i);
+ });
+ }
+
+ std::tuple<TS,TargetSpaces ...> data_;
+};
+
+ template< typename TS,typename ... TargetSpaces>
+ std::ostream& operator<<(std::ostream& s, const ProductManifold<TS, TargetSpaces ...>& c)
+ {
+ Dune::Hybrid::forEach(Dune::Hybrid::integralRange(Dune::index_constant<ProductManifold<TS, TargetSpaces ...>::numTS>()), [&](auto&& i) {
+ s<<Dune::className<decltype(c[i])>()<<" "<< c[i]<<"\n";
+ });
+ return s;
+ }
+}
+#endif
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..3ad4023be07c66fc0941471f44128c46d7a54799 100644
--- a/dune/gfe/rotation.hh
+++ b/dune/gfe/rotation.hh
@@ -146,7 +146,7 @@ class Rotation<T,3> : public Quaternion<T>
/** \brief Computes sin(x/2) / x without getting unstable for small x */
static T sincHalf(const T& x) {
using std::sin;
- return (x < 1e-4) ? 0.5 - (x*x/48) + (x*x*x*x)/3840 : sin(x/2)/x;
+ return (x < 1e-1) ? 0.5 - (x*x/48) + Dune::power(x,4)/3840 - Dune::power(x,6)/645120 : sin(x/2)/x;
}
/** \brief Computes sin(sqrt(x)/2) / sqrt(x) without getting unstable for small x
@@ -168,7 +168,15 @@ class Rotation<T,3> : public Quaternion<T>
static T sincOfSquare(const T& x) {
using std::sin;
using std::sqrt;
- return (x < 1e-15) ? 1 - (x/6) + (x*x)/120 : sin(sqrt(x))/sqrt(x);
+ // we need here lots of terms to be sure that the numerical derivatives are also within maschine precission
+ return (x < 1e-2) ?
+ 1-x/6
+ +x*x/120
+ -Dune::power(x,3)/5040
+ +Dune::power(x,4)/362880
+ -Dune::power(x,5)/39916800
+ +Dune::power(x,6)/6227020800
+ -Dune::power(x,7)/1307674368000: sin(sqrt(x))/sqrt(x);
}
public:
@@ -289,7 +297,7 @@ public:
// hence the series of cos(x/2) is
// 1 - x*x/8 + x*x*x*x/384 - ...
q[3] = (normV2 < 1e-4)
- ? 1 - normV2/8 + normV2*normV2 / 384
+ ? 1 - normV2/8 + normV2*normV2 / 384-Dune::power(normV2,3)/46080 + Dune::power(normV2,4)/10321920
: cos(sqrt(normV2)/2);
return q;
@@ -321,8 +329,8 @@ public:
// The series expansion of cos(x) at x=0 is
// 1 - x*x/2 + x*x*x*x/24 - ...
- T cosValue = (norm2 < 1e-4)
- ? 1 - norm2/2 + norm2*norm2 / 24
+ T cosValue = (norm2 < 1e-5)
+ ? 1 - norm2/2 + norm2*norm2 / 24 - Dune::power(norm2,3)/720+Dune::power(norm2,4)/40320
: cos(sqrt(norm2));
result *= cosValue;
@@ -604,7 +612,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 +630,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 +643,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 +915,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 +930,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
{
diff --git a/dune/gfe/unitvector.hh b/dune/gfe/unitvector.hh
index 90e498086693e6df381a8e0b9ccc6c9b4ba6f295..7122d1cda48f90638b8072f0211a6f1482a40baf 100644
--- a/dune/gfe/unitvector.hh
+++ b/dune/gfe/unitvector.hh
@@ -26,15 +26,24 @@ class UnitVector
/** \brief Computes sin(x) / x without getting unstable for small x */
static T sinc(const T& x) {
using std::sin;
- return (x < 1e-4) ? 1 - (x*x/6) : sin(x)/x;
+ return (x < 1e-2) ? 1.0-x*x/6.0+ Dune::power(x,4)/120.0-Dune::power(x,6)/5040.0+Dune::power(x,8)/362880.0 : sin(x)/x;
}
/** \brief Compute arccos^2 without using the (at 1) nondifferentiable function acos x close to 1 */
static T arcCosSquared(const T& x) {
using std::acos;
- const T eps = 1e-4;
- if (x > 1-eps) { // acos is not differentiable, use the series expansion instead
- return -2 * (x-1) + 1.0/3 * (x-1)*(x-1) - 4.0/45 * (x-1)*(x-1)*(x-1);
+ const T eps = 1e-2;
+ if (x > 1-eps) { // acos is not differentiable, use the series expansion instead,
+ // we need here lots of terms to be sure that the numerical derivatives are also within maschine precission
+ //return -2 * (x-1) + 1.0/3 * (x-1)*(x-1) - 4.0/45 * (x-1)*(x-1)*(x-1);
+ return 11665028.0/4729725.0
+ -141088.0/45045.0*x
+ + 413.0/429.0*x*x
+ - 5344.0/12285.0*Dune::power(x,3)
+ + 245.0/1287.0*Dune::power(x,4)
+ - 1632.0/25025.0*Dune::power(x,5)
+ + 56.0/3861.0*Dune::power(x,6)
+ - 32.0/21021.0*Dune::power(x,7);
} else {
return Dune::power(acos(x),2);
}
@@ -44,9 +53,19 @@ class UnitVector
static T derivativeOfArcCosSquared(const T& x) {
using std::acos;
using std::sqrt;
- const T eps = 1e-4;
+ const T eps = 1e-2;
if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
- return -2 + 2*(x-1)/3 - 4/15*(x-1)*(x-1);
+ // we need here lots of terms to be sure that the numerical derivatives are also within maschine precission
+ //return -2 + 2*(x-1)/3 - 4/15*(x-1)*(x-1);
+ return -47104.0/15015.0
+ +12614.0/6435.0*x
+ -63488.0/45045.0*x*x
+ + 1204.0/1287.0*Dune::power(x,3)
+ - 2048.0/4095.0*Dune::power(x,4)
+ + 112.0/585.0*Dune::power(x,5)
+ -2048.0/45045.0*Dune::power(x,6)
+ + 32.0/6435.0*Dune::power(x,7);
+
} else if (x < -1+eps) { // The function is not differentiable
DUNE_THROW(Dune::Exception, "arccos^2 is not differentiable at x==-1!");
} else
@@ -57,9 +76,20 @@ class UnitVector
static T secondDerivativeOfArcCosSquared(const T& x) {
using std::acos;
using std::pow;
- const T eps = 1e-4;
+ const T eps = 1e-2;
if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
- return 2.0/3 - 8*(x-1)/15;
+ // we need here lots of terms to be sure that the numerical derivatives are also within maschine precission
+ //return 2.0/3 - 8*(x-1)/15;
+ return 1350030.0/676039.0+5632.0/2028117.0*Dune::power(x,10)
+ -1039056896.0/334639305.0*x
+ +150876.0/39767.0*x*x
+ -445186048.0/111546435.0*Dune::power(x,3)
+ + 343728.0/96577.0*Dune::power(x,4)
+ - 57769984.0/22309287.0*Dune::power(x,5)
+ + 710688.0/482885.0*Dune::power(x,6)
+ - 41615360.0/66927861.0*Dune::power(x,7)
+ + 616704.0/3380195.0*Dune::power(x,8)
+ - 245760.0/7436429.0*Dune::power(x,9);
} else if (x < -1+eps) { // The function is not differentiable
DUNE_THROW(Dune::Exception, "arccos^2 is not differentiable at x==-1!");
} else
@@ -70,10 +100,22 @@ class UnitVector
static T thirdDerivativeOfArcCosSquared(const T& x) {
using std::acos;
using std::sqrt;
- const T eps = 1e-4;
+ const T eps = 1e-2;
if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
- return -8.0/15 + 24*(x-1)/35;
+ // we need here lots of terms to be sure that the numerical derivatives are also within maschine precission
+ //return -8.0/15 + 24*(x-1)/35;
+ return -1039056896.0/334639305.0
+ +301752.0/39767.0*x
+ -445186048.0/37182145.0*x*x
+ +1374912.0/96577.0*Dune::power(x,3)
+ -288849920.0/22309287.0*Dune::power(x,4)
+ +4264128.0/482885.0*Dune::power(x,5)
+ -41615360.0/9561123.0*Dune::power(x,6)
+ +4933632.0/3380195.0*Dune::power(x,7)
+ -2211840.0/7436429.0*Dune::power(x,8)
+ +56320.0/2028117.0*Dune::power(x,9);
} else if (x < -1+eps) { // The function is not differentiable
+
DUNE_THROW(Dune::Exception, "arccos^2 is not differentiable at x==-1!");
} else {
T d = 1-x*x;
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 5d1fa21fe65589709880d0aab13e31543cc5304f..e9e7a2e52a5fe2991ba47d352ed23973c8de096f 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
*/
@@ -256,6 +271,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
@@ -267,6 +284,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
@@ -277,5 +296,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