diff --git a/dune/common/CMakeLists.txt b/dune/common/CMakeLists.txt
index afecda540c139cdc4b666c2fc84e65ef492581dd..55ded18a7404f9c46b0eabfcc51bd901fc445d3f 100644
--- a/dune/common/CMakeLists.txt
+++ b/dune/common/CMakeLists.txt
@@ -104,6 +104,7 @@ install(FILES
tensordot.hh
tensormixin.hh
tensorspan.hh
+ tensortraits.hh
timer.hh
transpose.hh
tupleutility.hh
diff --git a/dune/common/concepts/archetypes/tensor.hh b/dune/common/concepts/archetypes/tensor.hh
index fe9e799c7a870b90df1c451876c87ef0050301c9..9934e5f803b67bcba246cc66bf2c194a4fec2827 100644
--- a/dune/common/concepts/archetypes/tensor.hh
+++ b/dune/common/concepts/archetypes/tensor.hh
@@ -12,6 +12,7 @@ namespace Dune::Concept::Archetypes {
template <std::size_t r>
struct Extents
{
+ using rank_type = std::size_t;
using index_type = std::size_t;
static constexpr std::size_t rank () { return r; }
diff --git a/dune/common/concepts/tensor.hh b/dune/common/concepts/tensor.hh
index 4b69c1a6873ca4ae91f22a7d7b228f33f37a25ea..15d991f0a9ec9346c9e52064c0c3c8038b685e0f 100644
--- a/dune/common/concepts/tensor.hh
+++ b/dune/common/concepts/tensor.hh
@@ -8,6 +8,7 @@
#include <array>
#include <concepts>
+#include <dune/common/tensortraits.hh>
#include <dune/common/concepts/archetypes/tensor.hh>
namespace Dune::Concept {
@@ -42,15 +43,15 @@ concept Extents = requires(E extents, std::size_t i)
* - `Dune::TensorSpan`
*/
template <class T>
-concept Tensor = Extents<T> && requires(T tensor)
+concept Tensor = requires(T tensor)
{
- requires Extents<typename T::extents_type>;
- { tensor.extents() } -> std::convertible_to<typename T::extents_type>;
+ requires Extents<typename TensorTraits<T>::extents_type>;
+ { TensorTraits<T>::extents(tensor) } -> std::convertible_to<typename TensorTraits<T>::extents_type>;
};
//! A `TensorWithRank` is a `Tensor` with given tensor-rank `rank`.
template <class T, std::size_t rank>
-concept TensorWithRank = Tensor<T> && T::rank() == rank;
+concept TensorWithRank = Tensor<T> && TensorTraits<T>::rank() == rank;
//! A `Vector` is a `Tensor` of rank 1.
template <class T>
@@ -75,14 +76,14 @@ static_assert(Concept::Matrix<Archetypes::Tensor<double,2>>);
*/
template <class T>
concept RandomAccessTensor = Tensor<T> &&
-requires(T tensor, std::array<typename T::index_type, T::rank()> indices)
+requires(T tensor, std::array<typename TensorTraits<T>::index_type, TensorTraits<T>::rank()> indices)
{
tensor[indices];
};
//! A `RandomAccessTensorWithRank` is a `RandomAccessTensor` with given tensor-rank `rank`.
template <class T, std::size_t rank>
-concept RandomAccessTensorWithRank = RandomAccessTensor<T> && T::rank() == rank;
+concept RandomAccessTensorWithRank = RandomAccessTensor<T> && TensorTraits<T>::rank() == rank;
//! A `RandomAccessVector` is a `RandomAccessTensor` of rank 1.
template <class T>
diff --git a/dune/common/densematrix.hh b/dune/common/densematrix.hh
index c696c85d3f94eace7c72e1c7099e01e5941d2261..7747368b001af88c20668249c2d3bcd6254d1719 100644
--- a/dune/common/densematrix.hh
+++ b/dune/common/densematrix.hh
@@ -201,6 +201,17 @@ namespace Dune
return asImp().mat_access(i);
}
+ //! random access with array of indices
+ value_type & operator[] (std::array<size_type,2> i)
+ {
+ return asImp().mat_access(i[0])[i[1]];
+ }
+
+ const value_type & operator[] (std::array<size_type,2> i) const
+ {
+ return asImp().mat_access(i[0])[i[1]];
+ }
+
//! size method (number of rows)
constexpr size_type size() const
{
diff --git a/dune/common/densevector.hh b/dune/common/densevector.hh
index cea5d13eeb8c7c99546c00f9c3fdc73bd2f26d27..26f538c4ddc2cb4145211686c450ba665ccab4e6 100644
--- a/dune/common/densevector.hh
+++ b/dune/common/densevector.hh
@@ -303,6 +303,21 @@ namespace Dune {
return asImp()[i];
}
+ //! random access with array of indices
+ template <class SizeType>
+ requires std::convertible_to<SizeType, size_type>
+ value_type & operator[] (std::array<SizeType,1> i)
+ {
+ return asImp()[i[0]];
+ }
+
+ template <class SizeType>
+ requires std::convertible_to<SizeType, size_type>
+ const value_type & operator[] (std::array<SizeType,1> i) const
+ {
+ return asImp()[i[0]];
+ }
+
//! return reference to first element
constexpr value_type& front()
{
diff --git a/dune/common/dynmatrix.hh b/dune/common/dynmatrix.hh
index ad2bd1809938861de9b2a2765d9d23990db2a76e..95815278f52c40b1b7fcc674b6feb695dafcd798 100644
--- a/dune/common/dynmatrix.hh
+++ b/dune/common/dynmatrix.hh
@@ -52,6 +52,30 @@ namespace Dune
typedef typename FieldTraits<K>::real_type real_type;
};
+ template< class K >
+ struct TensorTraits< DynamicMatrix<K> >
+ {
+ using index_type = typename DenseMatVecTraits<DynamicMatrix<K>>::size_type;
+ using extents_type = Std::extents<index_type, Std::dynamic_extent, Std::dynamic_extent>;
+ using rank_type = typename extents_type::rank_type;
+
+
+ /// \brief Number of elements in all dimensions of the array, \related extents
+ static constexpr extents_type extents (const DynamicMatrix<K>& tensor) noexcept { return extents_type{tensor.N(), tensor.M()}; }
+
+ /// \brief Number of dimensions of the array
+ static constexpr rank_type rank () noexcept { return 2; }
+
+ /// \brief Number of dimension with dynamic size
+ static constexpr rank_type rank_dynamic () noexcept { return 2; }
+
+ /// \brief Number of elements in the r'th dimension of the tensor
+ static constexpr std::size_t static_extent (rank_type /*r*/) noexcept { return Std::dynamic_extent; }
+
+ /// \brief Number of elements in the r'th dimension of the tensor
+ static constexpr index_type extent (const DynamicMatrix<K>& tensor, rank_type r) noexcept { return r == 0 ? tensor.N() : tensor.M(); }
+ };
+
/** \brief Construct a matrix with a dynamic size.
*
* \tparam K is the field type (use float, double, complex, etc)
diff --git a/dune/common/dynvector.hh b/dune/common/dynvector.hh
index ad54a52e9a44f567172c66c4c8b3c92f8aa32c35..773010b22e201bf45734c5885e201e5942a783c4 100644
--- a/dune/common/dynvector.hh
+++ b/dune/common/dynvector.hh
@@ -13,13 +13,16 @@
#include <initializer_list>
#include <limits>
#include <utility>
+#include <vector>
-#include "boundschecking.hh"
-#include "exceptions.hh"
-#include "genericiterator.hh"
+#include <dune/common/boundschecking.hh>
+#include <dune/common/densevector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/genericiterator.hh>
+#include <dune/common/tensortraits.hh>
+#include <dune/common/std/extents.hh>
+#include <dune/common/std/span.hh>
-#include <vector>
-#include "densevector.hh"
namespace Dune {
@@ -48,6 +51,31 @@ namespace Dune {
typedef typename FieldTraits< K >::real_type real_type;
};
+ template< class K, class Allocator >
+ struct TensorTraits< DynamicVector< K, Allocator > >
+ {
+ using index_type = typename DenseMatVecTraits<DynamicVector<K,Allocator>>::size_type;
+ using extents_type = Std::extents<index_type, Std::dynamic_extent>;
+ using rank_type = typename extents_type::rank_type;
+
+
+ /// \brief Number of elements in all dimensions of the array, \related extents
+ static constexpr extents_type extents (const DynamicVector<K,Allocator>& tensor) noexcept { return extents_type{tensor.size()}; }
+
+ /// \brief Number of dimensions of the array
+ static constexpr rank_type rank () noexcept { return 1; }
+
+ /// \brief Number of dimension with dynamic size
+ static constexpr rank_type rank_dynamic () noexcept { return 1; }
+
+ /// \brief Number of elements in the r'th dimension of the tensor
+ static constexpr std::size_t static_extent (rank_type /*r*/) noexcept { return Std::dynamic_extent; }
+
+ /// \brief Number of elements in the r'th dimension of the tensor
+ static constexpr index_type extent (const DynamicVector<K,Allocator>& tensor, rank_type /*r*/) noexcept { return tensor.size(); }
+ };
+
+
/** \brief Construct a vector with a dynamic size.
*
* \tparam K is the field type (use float, double, complex, etc)
@@ -158,6 +186,21 @@ namespace Dune {
return _data[i];
}
+ //! random access with array of indices
+ template <class SizeType>
+ requires std::convertible_to<SizeType, size_type>
+ K & operator[] (std::array<SizeType,1> i)
+ {
+ return _data[i[0]];
+ }
+
+ template <class SizeType>
+ requires std::convertible_to<SizeType, size_type>
+ const K & operator[] (std::array<SizeType,1> i) const
+ {
+ return _data[i[0]];
+ }
+
//! return pointer to underlying array
K* data() noexcept
{
diff --git a/dune/common/fmatrix.hh b/dune/common/fmatrix.hh
index 0ac53c7156c860300b42fb3436e02abeb843e890..d0c44eddce8972cddc1352226c1186f3b73e73cf 100644
--- a/dune/common/fmatrix.hh
+++ b/dune/common/fmatrix.hh
@@ -19,6 +19,8 @@
#include <dune/common/promotiontraits.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/matrixconcepts.hh>
+#include <dune/common/tensortraits.hh>
+#include <dune/common/std/extents.hh>
namespace Dune
{
@@ -104,6 +106,30 @@ namespace Dune
typedef typename FieldTraits<K>::real_type real_type;
};
+ template< class K, int ROWS, int COLS >
+ struct TensorTraits< FieldMatrix<K,ROWS,COLS> >
+ {
+ using index_type = typename DenseMatVecTraits<FieldMatrix<K,ROWS,COLS>>::size_type;
+ using extents_type = Std::extents<std::size_t, std::size_t(ROWS),std::size_t(COLS)>;
+ using rank_type = typename extents_type::rank_type;
+
+
+ /// \brief Number of elements in all dimensions of the array, \related extents
+ static constexpr extents_type extents (const FieldMatrix<K,ROWS,COLS>& /*tensor*/) noexcept { return extents_type{}; }
+
+ /// \brief Number of dimensions of the array
+ static constexpr rank_type rank () noexcept { return 2; }
+
+ /// \brief Number of dimension with dynamic size
+ static constexpr rank_type rank_dynamic () noexcept { return 0; }
+
+ /// \brief Number of elements in the r'th dimension of the tensor
+ static constexpr std::size_t static_extent (rank_type r) noexcept { return r == 0 ? ROWS : COLS; }
+
+ /// \brief Number of elements in the r'th dimension of the tensor
+ static constexpr index_type extent (const FieldMatrix<K,ROWS,COLS>& /*tensor*/, rank_type r) noexcept { return r == 0 ? ROWS : COLS; }
+ };
+
/**
@brief A dense n x m matrix.
diff --git a/dune/common/fvector.hh b/dune/common/fvector.hh
index 605988abe311df2d0b8afa46a5c34e8e3bc5e010..fa41dae5db323ac697c73b348e8fdd20750bd754 100644
--- a/dune/common/fvector.hh
+++ b/dune/common/fvector.hh
@@ -19,8 +19,10 @@
#include <dune/common/ftraits.hh>
#include <dune/common/math.hh>
#include <dune/common/promotiontraits.hh>
+#include <dune/common/tensortraits.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/typeutilities.hh>
+#include <dune/common/std/extents.hh>
namespace Dune {
@@ -50,6 +52,30 @@ namespace Dune {
typedef typename FieldTraits<K>::real_type real_type;
};
+ template< class K, int SIZE >
+ struct TensorTraits< FieldVector<K,SIZE> >
+ {
+ using index_type = typename DenseMatVecTraits<FieldVector<K,SIZE>>::size_type;
+ using extents_type = Std::extents<std::size_t, std::size_t(SIZE)>;
+ using rank_type = typename extents_type::rank_type;
+
+
+ /// \brief Number of elements in all dimensions of the array, \related extents
+ static constexpr extents_type extents (const FieldVector<K,SIZE>& /*tensor*/) noexcept { return extents_type{}; }
+
+ /// \brief Number of dimensions of the array
+ static constexpr rank_type rank () noexcept { return 1; }
+
+ /// \brief Number of dimension with dynamic size
+ static constexpr rank_type rank_dynamic () noexcept { return 0; }
+
+ /// \brief Number of elements in the r'th dimension of the tensor
+ static constexpr std::size_t static_extent (rank_type /*r*/) noexcept { return SIZE; }
+
+ /// \brief Number of elements in the r'th dimension of the tensor
+ static constexpr index_type extent (const FieldVector<K,SIZE>& /*tensor*/, rank_type /*r*/) noexcept { return SIZE; }
+ };
+
/**
* @brief TMP to check the size of a DenseVectors statically, if possible.
*
@@ -209,6 +235,21 @@ namespace Dune {
return _data[i];
}
+ //! random access with array of indices
+ template <class SizeType>
+ requires std::convertible_to<SizeType, size_type>
+ reference operator[] (std::array<SizeType,1> i)
+ {
+ return _data[i[0]];
+ }
+
+ template <class SizeType>
+ requires std::convertible_to<SizeType, size_type>
+ const_reference operator[] (std::array<SizeType,1> i) const
+ {
+ return _data[i[0]];
+ }
+
//! Return pointer to underlying array
constexpr K* data () noexcept
{
@@ -402,6 +443,21 @@ namespace Dune {
return _data;
}
+ //! random access with array of indices
+ template <class SizeType>
+ requires std::convertible_to<SizeType, size_type>
+ reference operator[] (std::array<SizeType,1> i)
+ {
+ return _data;
+ }
+
+ template <class SizeType>
+ requires std::convertible_to<SizeType, size_type>
+ const_reference operator[] (std::array<SizeType,1> i) const
+ {
+ return _data;
+ }
+
//! return pointer to underlying array
constexpr K* data () noexcept
{
diff --git a/dune/common/tensordot.hh b/dune/common/tensordot.hh
index ef38b6dcb72ea7031d5ff9df2eaf8c4b5aca94b7..850300bba1e8f8616b84f217468d635ce1718619 100644
--- a/dune/common/tensordot.hh
+++ b/dune/common/tensordot.hh
@@ -14,6 +14,7 @@
#include <dune/common/integersequence.hh>
#include <dune/common/std/extents.hh>
#include <dune/common/rangeutilities.hh>
+#include <dune/common/tensortraits.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/concepts/tensor.hh>
#include <dune/common/std/extents.hh>
@@ -103,20 +104,23 @@ for (int l = 0; l < a.extent(aSeq[0]); ++l)
template <std::size_t K = 0,
class A, class ASeq, class ASeqInv,
class B, class BSeq, class BSeqInv,
- class C, class BinaryOp1, class BinaryOp2>
+ class C, class BinaryOp1, class BinaryOp2,
+ class ATraits = TensorTraits<A>,
+ class BTraits = TensorTraits<B>,
+ class CTraits = TensorTraits<C>>
constexpr DUNE_FORCE_INLINE
void tensorDotImpl (const A& a, ASeq aSeq, ASeqInv aSeqInv,
const B& b, BSeq bSeq, BSeqInv bSeqInv,
C& c, BinaryOp1 op1, BinaryOp2 op2,
- std::array<typename A::index_type,A::rank()> aIndices = {},
- std::array<typename B::index_type,B::rank()> bIndices = {},
- std::array<typename C::index_type,C::rank()> cIndices = {})
+ std::array<typename ATraits::index_type,ATraits::rank()> aIndices = {},
+ std::array<typename BTraits::index_type,BTraits::rank()> bIndices = {},
+ std::array<typename CTraits::index_type,CTraits::rank()> cIndices = {})
{
if constexpr(aSeq.size() > 0 && bSeq.size() > 0) {
// first, loop over the contraction indices of A and B
constexpr std::size_t I = head(aSeq);
constexpr std::size_t J = head(bSeq);
- for (typename A::index_type k = 0; k < a.extent(I); ++k) {
+ for (typename ATraits::index_type k = 0; k < ATraits::extent(a,I); ++k) {
aIndices[I] = k;
bIndices[J] = k;
tensorDotImpl<K>(a,tail(aSeq),aSeqInv,b,tail(bSeq),bSeqInv,c,op1,op2,aIndices,bIndices,cIndices);
@@ -125,7 +129,7 @@ void tensorDotImpl (const A& a, ASeq aSeq, ASeqInv aSeqInv,
else if constexpr(aSeqInv.size() > 0) {
// second, loop over the remaining indices of tensor A
constexpr std::size_t I = head(aSeqInv);
- for (typename A::index_type i = 0; i < a.extent(I); ++i) {
+ for (typename ATraits::index_type i = 0; i < ATraits::extent(a,I); ++i) {
aIndices[I] = i;
cIndices[K] = i;
tensorDotImpl<K+1>(a,aSeq,tail(aSeqInv),b,bSeq,bSeqInv,c,op1,op2,aIndices,bIndices,cIndices);
@@ -134,7 +138,7 @@ void tensorDotImpl (const A& a, ASeq aSeq, ASeqInv aSeqInv,
else if constexpr(bSeqInv.size() > 0) {
// third, loop over the remaining indices of tensor B
constexpr std::size_t J = head(bSeqInv);
- for (typename B::index_type j = 0; j < b.extent(J); ++j) {
+ for (typename BTraits::index_type j = 0; j < BTraits::extent(b,J); ++j) {
bIndices[J] = j;
cIndices[K] = j;
tensorDotImpl<K+1>(a,aSeq,aSeqInv,b,bSeq,tail(bSeqInv),c,op1,op2,aIndices,bIndices,cIndices);
@@ -170,16 +174,20 @@ constexpr auto tensordotOut (const A& a, std::index_sequence<II...> aSeq,
{
static_assert(aSeq.size() == bSeq.size());
+ using ATraits = TensorTraits<A>;
+ using BTraits = TensorTraits<B>;
+ using CTraits = TensorTraits<C>;
+
// create integer sequences that do not include the contraction indices
- const auto aSeqInv = difference<A::rank()>(aSeq); // {0,1,...A::rank()-1} \ {II...}
- const auto bSeqInv = difference<B::rank()>(bSeq); // {0,1,...B::rank()-1} \ {JJ...}
- static_assert(aSeqInv.size() + bSeqInv.size() == C::rank());
+ const auto aSeqInv = difference<ATraits::rank()>(aSeq); // {0,1,...A::rank()-1} \ {II...}
+ const auto bSeqInv = difference<BTraits::rank()>(bSeq); // {0,1,...B::rank()-1} \ {JJ...}
+ static_assert(aSeqInv.size() + bSeqInv.size() == CTraits::rank());
// the extents of a and the extents of b must be compatible to c
- using EA = typename A::extents_type;
- using EB = typename B::extents_type;
+ using EA = typename ATraits::extents_type;
+ using EB = typename BTraits::extents_type;
static_assert(Impl::checkStaticExtents<EA,EB>(aSeq, bSeq));
- assert((Impl::checkExtents(a.extents(), aSeq, b.extents(), bSeq)));
+ assert((Impl::checkExtents(ATraits::extents(a), aSeq, BTraits::extents(b), bSeq)));
// Objects of `A` and `B` must be different from object `C`
assert((void*)(&a) != (void*)(&c) && (void*)(&b) != (void*)(&c));
@@ -208,17 +216,21 @@ constexpr void tensordotOut (const A& a, const B& b, C& c,
std::integral_constant<std::size_t,N> axes = {},
BinaryOp1 op1 = {}, BinaryOp2 op2 = {})
{
- using SeqI = typename StaticIntegralRange<std::size_t,A::rank(),A::rank()-N>::integer_sequence;
- using InvSeqI = std::make_index_sequence<A::rank()-N>;
+ using ATraits = TensorTraits<A>;
+ using BTraits = TensorTraits<B>;
+ using CTraits = TensorTraits<C>;
+
+ using SeqI = typename StaticIntegralRange<std::size_t,ATraits::rank(),ATraits::rank()-N>::integer_sequence;
+ using InvSeqI = std::make_index_sequence<ATraits::rank()-N>;
using SeqJ = std::make_index_sequence<N>;
- using InvSeqJ = typename StaticIntegralRange<std::size_t,B::rank(),N>::integer_sequence;
- static_assert(InvSeqI::size() + InvSeqJ::size() == C::rank());
+ using InvSeqJ = typename StaticIntegralRange<std::size_t,BTraits::rank(),N>::integer_sequence;
+ static_assert(InvSeqI::size() + InvSeqJ::size() == CTraits::rank());
// the extents of a and the extents of b must be compatible to c
- using EA = typename A::extents_type;
- using EB = typename B::extents_type;
+ using EA = typename ATraits::extents_type;
+ using EB = typename BTraits::extents_type;
static_assert(Impl::checkStaticExtents<EA,EB>(SeqI{}, SeqJ{}));
- assert((Impl::checkExtents(a.extents(), SeqI{}, b.extents(), SeqJ{})));
+ assert((Impl::checkExtents(ATraits::extents(a), SeqI{}, BTraits::extents(b), SeqJ{})));
// Objects of `A` and `B` must be different from object `C`
assert((void*)(&a) != (void*)(&c) && (void*)(&b) != (void*)(&c));
@@ -252,23 +264,26 @@ constexpr auto tensordot (const A& a, std::index_sequence<II...> aSeq,
const B& b, std::index_sequence<JJ...> bSeq,
BinaryOp1 op1 = {}, BinaryOp2 op2 = {})
{
+ using ATraits = TensorTraits<A>;
+ using BTraits = TensorTraits<B>;
+
// the extents(II) of a and the extents(JJ) of b must match
- using EA = typename A::extents_type;
- using EB = typename B::extents_type;
+ using EA = typename ATraits::extents_type;
+ using EB = typename BTraits::extents_type;
static_assert(Impl::checkStaticExtents<EA,EB>(aSeq, bSeq));
- assert((Impl::checkExtents(a.extents(), aSeq, b.extents(), bSeq)));
+ assert((Impl::checkExtents(ATraits::extents(a), aSeq, BTraits::extents(b), bSeq)));
// create integer sequences that do not include the contraction indices
- const auto aSeqInv = difference<A::rank()>(aSeq); // {0,1,...A::rank()-1} \ {II...}
- const auto bSeqInv = difference<B::rank()>(bSeq); // {0,1,...B::rank()-1} \ {JJ...}
+ const auto aSeqInv = difference<ATraits::rank()>(aSeq); // {0,1,...A::rank()-1} \ {II...}
+ const auto bSeqInv = difference<BTraits::rank()>(bSeq); // {0,1,...B::rank()-1} \ {JJ...}
// create result extents by collecting the extents of a and b that are not contracted
auto cExtents = Impl::concatExtents(
- Impl::sliceExtents(a.extents(), aSeqInv),
- Impl::sliceExtents(b.extents(), bSeqInv));
+ Impl::sliceExtents(ATraits::extents(a), aSeqInv),
+ Impl::sliceExtents(BTraits::extents(b), bSeqInv));
- using VA = decltype(std::declval<A>()[std::array<typename A::index_type,A::rank()>{}]);
- using VB = decltype(std::declval<B>()[std::array<typename B::index_type,B::rank()>{}]);
+ using VA = decltype(std::declval<A>()[std::array<typename ATraits::index_type,ATraits::rank()>{}]);
+ using VB = decltype(std::declval<B>()[std::array<typename BTraits::index_type,BTraits::rank()>{}]);
using V = std::invoke_result_t<BinaryOp2,VA,VB>;
auto c = Tensor{cExtents, V(0)};
Impl::tensorDotImpl(a,aSeq,aSeqInv,b,bSeq,bSeqInv,c,std::ref(op1),std::ref(op2));
@@ -296,7 +311,8 @@ constexpr auto tensordot (const A& a, const B& b,
std::integral_constant<std::size_t,N> axes = {},
BinaryOp1 op1 = {}, BinaryOp2 op2 = {})
{
- using SeqI = typename StaticIntegralRange<std::size_t,A::rank(),A::rank()-N>::integer_sequence;
+ using ATraits = TensorTraits<A>;
+ using SeqI = typename StaticIntegralRange<std::size_t,ATraits::rank(),ATraits::rank()-N>::integer_sequence;
using SeqJ = std::make_index_sequence<N>;
return tensordot(a,SeqI{},b,SeqJ{},op1,op2);
}
diff --git a/dune/common/tensormixin.hh b/dune/common/tensormixin.hh
index 4a05ea27248cd3482ed550ca0849aae6600caf0b..93734663a6ce929b85182a09b4ce40597de7bb74 100644
--- a/dune/common/tensormixin.hh
+++ b/dune/common/tensormixin.hh
@@ -651,6 +651,49 @@ constexpr bool operator== (const S& number, const TensorMixin<D,B>& rhs) noexcep
return number == rhs();
}
+
+/** \brief Output stream overload for tensor types */
+template <class D, class B>
+std::ostream& operator<< (std::ostream& out, const Dune::TensorMixin<D,B>& tensor)
+{
+ using extents_type = typename Dune::TensorMixin<D,B>::extents_type;
+ using index_type = typename Dune::TensorMixin<D,B>::index_type;
+ if constexpr(extents_type::rank() == 0) {
+ out << tensor();
+ } else if constexpr(extents_type::rank() == 1) {
+ out << "[";
+ for (index_type i = 0; i < tensor.extent(0); ++i)
+ out << tensor(i) << (i < tensor.extent(0)-1 ? ", " : "]");
+ } else if constexpr(extents_type::rank() == 2) {
+ out << "[\n";
+ for (index_type i = 0; i < tensor.extent(0); ++i) {
+ out << " [";
+ for (index_type j = 0; j < tensor.extent(1); ++j)
+ out << tensor(i,j) << (j < tensor.extent(1)-1 ? ", " : "]");
+ out << (i < tensor.extent(0)-1 ? ",\n" : "\n");
+ }
+ out << ']';
+ } else if constexpr(extents_type::rank() == 3) {
+ out << "[\n";
+ for (index_type i = 0; i < tensor.extent(0); ++i) {
+ out << " [\n";
+ for (index_type j = 0; j < tensor.extent(1); ++j) {
+ out << " [";
+ for (index_type k = 0; k < tensor.extent(2); ++k)
+ out << tensor(i,j,k) << (k < tensor.extent(2)-1 ? ", " : "]");
+ out << (j < tensor.extent(1)-1 ? ",\n" : "\n");
+ }
+ out << " ]";
+ out << (i < tensor.extent(0)-1 ? ",\n" : "\n");
+ }
+ out << ']';
+ } else {
+ out << "Tensor<" << extents_type::rank() << ">";
+ }
+ return out;
+}
+
+
template <class D, class B>
struct FieldTraits< TensorMixin<D,B> >
{
diff --git a/dune/common/tensortraits.hh b/dune/common/tensortraits.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0bd0d5f13310fbaad402a9fc950eb7ccb5aea27d
--- /dev/null
+++ b/dune/common/tensortraits.hh
@@ -0,0 +1,38 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
+// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
+#ifndef DUNE_COMMON_TENSORTRAITS_HH
+#define DUNE_COMMON_TENSORTRAITS_HH
+
+#include <cstddef>
+
+namespace Dune {
+
+template <class T>
+struct TensorTraits
+{
+ using extents_type = typename T::extents_type;
+ using index_type = typename extents_type::index_type;
+ using rank_type = typename extents_type::rank_type;
+
+
+ /// \brief Number of elements in all dimensions of the array, \related extents
+ static constexpr const extents_type& extents (const T& tensor) noexcept { return tensor.extents(); }
+
+ /// \brief Number of dimensions of the array
+ static constexpr rank_type rank () noexcept { return extents_type::rank(); }
+
+ /// \brief Number of dimension with dynamic size
+ static constexpr rank_type rank_dynamic () noexcept { return extents_type::rank_dynamic(); }
+
+ /// \brief Number of elements in the r'th dimension of the tensor
+ static constexpr std::size_t static_extent (rank_type r) noexcept { return extents_type::static_extent(r); }
+
+ /// \brief Number of elements in the r'th dimension of the tensor
+ static constexpr index_type extent (const T& tensor, rank_type r) noexcept { return extents(tensor).extent(r); }
+};
+
+} // end namespace Dune
+
+#endif // DUNE_COMMON_TENSORTRAITS_HH
diff --git a/dune/common/test/tensordottest.cc b/dune/common/test/tensordottest.cc
index 3a00235a47a1c5cb27ba0502cc31097f726baf92..89bb43df942e5e53b9ad3f37200fa70782626e35 100644
--- a/dune/common/test/tensordottest.cc
+++ b/dune/common/test/tensordottest.cc
@@ -3,6 +3,10 @@
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/dynvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
#include <dune/common/tensor.hh>
#include <dune/common/tensordot.hh>
#include <dune/common/tensorspan.hh>
@@ -23,6 +27,11 @@ int main(int argc, char** argv)
auto dTensor32 = Dune::Tensor<double,Dune::dynamic,Dune::dynamic>{3,2};
auto dTensor234 = Dune::Tensor<double,Dune::dynamic,Dune::dynamic,Dune::dynamic>{2,3,4};
+ auto dVector2 = Dune::DynamicVector<double>{2};
+ auto dVector3 = Dune::DynamicVector<double>{3};
+ auto dMatrix23 = Dune::DynamicMatrix<double>{2,3};
+ auto dMatrix32 = Dune::DynamicMatrix<double>{3,2};
+
auto fTensor = Dune::Tensor<double>{};
auto fTensor2 = Dune::Tensor<double,2>{};
auto fTensor3 = Dune::Tensor<double,3>{};
@@ -30,6 +39,11 @@ int main(int argc, char** argv)
auto fTensor32 = Dune::Tensor<double,3,2>{};
auto fTensor234 = Dune::Tensor<double,2,3,4>{};
+ auto fVector2 = Dune::FieldVector<double,2>{};
+ auto fVector3 = Dune::FieldVector<double,3>{};
+ auto fMatrix23 = Dune::FieldMatrix<double,2,3>{};
+ auto fMatrix32 = Dune::FieldMatrix<double,3,2>{};
+
// test dynamic tensors
{
auto d = tensordot<0>(dTensor,dTensor);
@@ -54,6 +68,11 @@ int main(int argc, char** argv)
testSuite.check(d23.extent(1) == 3);
}
+ {
+ auto d = tensordot<1>(dVector2,dVector2);
+ testSuite.check(d.rank() == 0);
+ }
+
{
auto d = tensordot<2>(dTensor23,dTensor23);
testSuite.check(d.rank() == 0);
@@ -101,6 +120,43 @@ int main(int argc, char** argv)
testSuite.check(d2332.extent(3) == 2);
}
+ {
+ auto d = tensordot<2>(dMatrix23,dMatrix23);
+ testSuite.check(d.rank() == 0);
+
+ auto d2 = tensordot<1>(dMatrix23,dVector3);
+ testSuite.check(d2.rank() == 1);
+ testSuite.check(d2.extent(0) == 2);
+
+ auto d3 = tensordot<1>(dVector2,dMatrix23);
+ testSuite.check(d3.rank() == 1);
+ testSuite.check(d3.extent(0) == 3);
+
+ auto d22 = tensordot<1>(dMatrix23,dMatrix32);
+ testSuite.check(d22.rank() == 2);
+ testSuite.check(d22.extent(0) == 2);
+ testSuite.check(d22.extent(1) == 2);
+
+ auto d223 = tensordot<0>(dVector2,dMatrix23);
+ testSuite.check(d223.rank() == 3);
+ testSuite.check(d223.extent(0) == 2);
+ testSuite.check(d223.extent(1) == 2);
+ testSuite.check(d223.extent(2) == 3);
+
+ auto d233 = tensordot<0>(dMatrix23,dVector3);
+ testSuite.check(d233.rank() == 3);
+ testSuite.check(d233.extent(0) == 2);
+ testSuite.check(d233.extent(1) == 3);
+ testSuite.check(d233.extent(2) == 3);
+
+ auto d2332 = tensordot<0>(dMatrix23,dMatrix32);
+ testSuite.check(d2332.rank() == 4);
+ testSuite.check(d2332.extent(0) == 2);
+ testSuite.check(d2332.extent(1) == 3);
+ testSuite.check(d2332.extent(2) == 3);
+ testSuite.check(d2332.extent(3) == 2);
+ }
+
// test mixed static/dynamic tensors
{
auto d = tensordot<0>(fTensor,dTensor);
@@ -125,6 +181,16 @@ int main(int argc, char** argv)
testSuite.check(d23.extent(1) == 3);
}
+ {
+ auto d = tensordot<1>(fVector2,dVector2);
+ testSuite.check(d.rank() == 0);
+
+ auto d23 = tensordot<0>(fVector2,dVector3);
+ testSuite.check(d23.rank() == 2);
+ testSuite.check(d23.extent(0) == 2);
+ testSuite.check(d23.extent(1) == 3);
+ }
+
{
auto d = tensordot<2>(fTensor23,dTensor23);
testSuite.check(d.rank() == 0);
@@ -172,6 +238,24 @@ int main(int argc, char** argv)
testSuite.check(d2332.extent(3) == 2);
}
+ {
+ auto d = tensordot<2>(fMatrix23,dTensor23);
+ testSuite.check(d.rank() == 0);
+
+ auto d2 = tensordot<1>(fMatrix23,dTensor3);
+ testSuite.check(d2.rank() == 1);
+ testSuite.check(d2.extent(0) == 2);
+
+ auto d3 = tensordot<1>(fVector2,dTensor23);
+ testSuite.check(d3.rank() == 1);
+ testSuite.check(d3.extent(0) == 3);
+
+ auto d22 = tensordot<1>(fMatrix23,dTensor32);
+ testSuite.check(d22.rank() == 2);
+ testSuite.check(d22.extent(0) == 2);
+ testSuite.check(d22.extent(1) == 2);
+ }
+
// test interaction with TensorSpan
{
auto dMat0 = tensordot<1>(fTensor23,dTensor32.toTensorSpan());