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