Tensor data structure based on mdarray

dune/common/CMakeLists.txt

@@ -100,6 +100,9 @@ install(FILES
         stdthread.hh
         streamoperators.hh
         stringutility.hh
+        tensor.hh
+        tensorinterface.hh
+        tensorspan.hh
         timer.hh
         transpose.hh
         tupleutility.hh

dune/common/initializerlist.hh

+// -*- 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_INITIALIZER_LIST_HH
+#define DUNE_COMMON_INITIALIZER_LIST_HH
+
+#include <cassert>
+#include <initializer_list>
+
+namespace Dune {
+namespace Impl {
+
+template <class Value, int rank>
+struct NestedInitializerList
+{
+  using type = std::initializer_list<typename NestedInitializerList<Value,rank-1>::type>;
+};
+
+template <class Value>
+struct NestedInitializerList<Value,1>
+{
+  using type = std::initializer_list<Value>;
+};
+
+template <class Value>
+struct NestedInitializerList<Value,0>
+{
+  using type = Value;
+};
+
+} // end namespace Impl
+
+/// \brief A nested `std::initializer_list<std::initializer_list<...>>` up to depth `rank`
+template <class Value, int rank>
+using NestedInitializerList_t = typename Impl::NestedInitializerList<Value,rank>::type;
+
+
+/// \brief A utility to recursively unpack nested initializer lists
+template <class Value, class Extents, int I = Extents::rank()>
+class InitializerList
+{
+public:
+  using value_type = Value;
+  using extents_type = Extents;
+
+  template <class F>
+  static void apply (NestedInitializerList_t<Value,I> values, const extents_type& extents, F&& set_value)
+  {
+    assert(values.size() == std::size_t(extents.extent(extents_type::rank()-I)));
+
+    // process all the sub lists
+    for (auto&& sub : values)
+      InitializerList<value_type,extents_type,I-1>::apply(sub, extents, set_value);
+  }
+};
+
+#ifndef DOXYGEN
+template <class Value, class Extents>
+class InitializerList<Value,Extents,1>
+{
+public:
+  using value_type = Value;
+  using extents_type = Extents;
+
+  template <class F>
+  static void apply (std::initializer_list<Value> values, const extents_type& extents, F&& set_value)
+  {
+    assert(values.size() == std::size_t(extents.extent(extents_type::rank()-1)));
+
+    for (auto&& value : values)
+      set_value(value);
+  }
+};
+
+template <class Value, class Extents>
+class InitializerList<Value,Extents,0>
+{
+public:
+  using value_type = Value;
+  using extents_type = Extents;
+
+  template <class F>
+  static void apply (const Value& value, const extents_type& extents, F&& set_value)
+  {
+    set_value(value);
+  }
+};
+#endif // DOXYGEN
+
+} // end namespace Dune
+
+#endif // DUNE_COMMON_INITIALIZER_LIST_HH

dune/common/tensor.hh

+// -*- 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_TENSOR_HH
+#define DUNE_COMMON_TENSOR_HH
+
+#include <array>
+#include <type_traits>
+#include <vector>
+
+#include <dune/common/initializerlist.hh>
+#include <dune/common/tensormixin.hh>
+#include <dune/common/tensorspan.hh>
+#include <dune/common/std/extents.hh>
+#include <dune/common/std/layout_right.hh>
+#include <dune/common/std/mdarray.hh>
+#include <dune/common/std/span.hh>
+
+namespace Dune {
+
+//! A special value representing dynamic extents in a tensor
+static constexpr inline std::size_t dynamic = Std::dynamic_extent;
+
+namespace Impl {
+
+template <class Element, std::size_t... extents>
+struct TensorTraits
+{
+  using container_type = std::conditional_t<((extents == Dune::dynamic) || ...),
+    std::vector<Element>, std::array<Element, (extents * ... * std::size_t(1))>>;
+  using extents_type = Std::extents<typename container_type::size_type,extents...>;
+  using storage_type = Std::mdarray<Element, extents_type, Std::layout_right, container_type>;
+};
+
+} // end namespace Impl
+
+
+/**
+ * \brief A tensor of arbitrary rank and individual dimension extents, static
+ * or dynamic, storing field values in a container of corresponding static or
+ * dynamic size.
+ *
+ * \ingroup Tensors
+ * \nosubgrouping
+ *
+ * \tparam Element  The element type stored in the tensor
+ * \tparam extents  Individual static extents or Dune::dynamic
+ **/
+template <class Element, std::size_t... extents>
+class Tensor
+    : public TensorMixin<Tensor<Element,extents...>,
+        typename Impl::TensorTraits<Element,extents...>::storage_type>
+{
+  using self_type = Tensor;
+  using storage_type = typename Impl::TensorTraits<Element,extents...>::storage_type;
+  using base_type = TensorMixin<self_type, storage_type>;
+
+public:
+  using element_type = typename base_type::element_type;
+  using extents_type = typename base_type::extents_type;
+  using value_type = typename base_type::value_type;
+  using index_type = typename base_type::index_type;
+  using layout_type = typename base_type::layout_type;
+  using mapping_type = typename base_type::mapping_type;
+
+  using tensorspan_type = TensorSpan<element_type,extents_type,layout_type>;
+  using const_tensorspan_type = TensorSpan<const element_type,extents_type,layout_type>;
+
+public:
+  /// \brief Inherit constructor from TensorMixin
+  using base_type::base_type;
+
+  /// \name Additional Tensor constructors
+  /// @{
+
+  explicit constexpr Tensor (const value_type& value)
+    : base_type{extents_type{}, value}
+  {}
+
+  /**
+   * \brief Constructor from a brace-init list of values
+   *
+   * \b Example:
+   * \code{.cpp}
+   Tensor<double,2,2> matrix{
+     {1.0,2.0},
+     {2.0, 3.0}
+   };
+   * \endcode
+   **/
+  template <class E = extents_type,
+    std::enable_if_t<(E::rank_dynamic() == 0), int> = 0,
+    std::enable_if_t<std::is_default_constructible_v<E>, int> = 0>
+  constexpr Tensor (NestedInitializerList_t<value_type,extents_type::rank()> init)
+    : Tensor{extents_type{}, init}
+  {}
+
+  /**
+   * \brief Constructor from an extents and brace-init list of values
+   *
+   * \b Example:
+   * \code{.cpp}
+   using Extents = typename Tensor<double,dynamic,dynamic>::extents_type;
+   Tensor<double,dynamic,dynamic> matrix(Extents{2,2},
+   {
+     {1.0,2.0},
+     {2.0, 3.0}
+   });
+   * \endcode
+   **/
+  template <class M = mapping_type,
+    std::enable_if_t<std::is_constructible_v<M,extents_type>, int> = 0>
+  constexpr Tensor (const extents_type& e, NestedInitializerList_t<value_type,extents_type::rank()> init)
+    : Tensor{mapping_type{e}, init}
+  {}
+
+  /**
+   * \brief Constructor from a span of extents and brace-init list of values
+   *
+   * \b Example:
+   * \code{.cpp}
+   Tensor<double,dynamic,dynamic> matrix(std::array{2,2},
+   {
+     {1.0,2.0},
+     {2.0, 3.0}
+   });
+   * \endcode
+   **/
+  template <class Otherindex_type, std::size_t N,
+    std::enable_if_t<std::is_convertible_v<const Otherindex_type&, index_type>, int> = 0,
+    std::enable_if_t<std::is_nothrow_constructible_v<index_type,const Otherindex_type&>, int> = 0,
+    std::enable_if_t<(N == extents_type::rank_dynamic() || N == extents_type::rank()), int> = 0>
+  constexpr Tensor (Std::span<Otherindex_type,N> e,
+                    NestedInitializerList_t<value_type,extents_type::rank()> init)
+    : Tensor{mapping_type{extents_type{e}}, init}
+  {}
+
+  /**
+   * \brief Constructor from a layout-mapping and brace-init list of values
+   *
+   * \b Example:
+   * \code{.cpp}
+   using Extents = typename Tensor<double,dynamic,dynamic>::extents_type;
+   using Mapping = typename Tensor<double,dynamic,dynamic>::mapping_type;
+   Tensor<double,dynamic,dynamic> matrix(Mapping{Extents{2,2}},
+   {
+     {1.0,2.0},
+     {2.0, 3.0}
+   });
+   * \endcode
+   **/
+  constexpr Tensor (const mapping_type& m, NestedInitializerList_t<value_type,extents_type::rank()> init)
+    : base_type{m}
+  {
+    auto it = this->container_data();
+    InitializerList<value_type,extents_type>::apply(init,this->extents(),
+      [&it](value_type value) { *it++ = value; });
+  }
+
+  /// \brief Converting constructor from another TensorMixin
+  template <class D, class B,
+    std::enable_if_t<std::is_constructible_v<base_type, B>, int> = 0>
+  constexpr Tensor (const TensorMixin<D,B>& other)
+    : base_type{Std::mdspan(other)}
+  {}
+
+  /// @}
+
+
+  /// \name Multi index access
+  /// @{
+
+  /**
+   * \brief Subscript operator to access the tensor components using an array of indices.
+   * \b Examples:
+   * \code{c++}
+     Tensor<double,3,3> matrix;
+     matrix[std::array{0,1}] = 7.0;
+     \endcode
+   **/
+  using base_type::operator[];
+
+  /**
+   * \brief Access an element of the tensor using a variadic list of indices.
+   * \b Examples:
+   * \code{c++}
+     Tensor<double,3,3,3> tensor;
+     tensor(0,1,2) = 42.0;
+     \endcode
+   **/
+  using base_type::operator();
+
+  /// @}
+
+
+  /// \name Modifiers
+  /// @{
+
+  /// \brief Change the extents of the tensor and resize the underlying container with given default value
+  template <class V,
+    std::enable_if_t<std::is_convertible_v<V,value_type>, int> = 0>
+  void resize (const extents_type& e, const V& value)
+  {
+    auto container = std::move(*this).extract_container();
+    auto m = mapping_type{e};
+    container.resize(m.required_span_size(), value);
+    static_cast<base_type&>(*this) = base_type{m, std::move(container)};
+  }
+
+  /// \brief Change the extents of the tensor and resize the underlying container
+  void resize (const extents_type& e)
+  {
+    resize(e, value_type(0));
+  }
+
+  /// \brief Change the extents of the tensor by the given individual extents
+  template <class... index_types, class V,
+    std::enable_if_t<(sizeof...(index_types) == extents_type::rank() ||
+                      sizeof...(index_types) == extents_type::rank_dynamic()), int> = 0,
+    std::enable_if_t<(... && std::is_convertible_v<index_types,index_type>), int> = 0,
+    std::enable_if_t<std::is_constructible_v<extents_type,index_types...>, int> = 0,
+    std::enable_if_t<std::is_convertible_v<V,value_type>, int> = 0>
+  void resize (index_types... exts, const V& v)
+  {
+    resize(extents_type{exts...}, v);
+  }
+
+  /// \brief Change the extents of the tensor by the given individual extents
+  template <class... index_types,
+    std::enable_if_t<(sizeof...(index_types) == extents_type::rank() ||
+                      sizeof...(index_types) == extents_type::rank_dynamic()), int> = 0,
+    std::enable_if_t<(... && std::is_convertible_v<index_types,index_type>), int> = 0,
+    std::enable_if_t<std::is_constructible_v<extents_type,index_types...>, int> = 0>
+  void resize (index_types... exts)
+  {
+    resize(extents_type{exts...}, value_type(0));
+  }
+
+  /// @}
+
+
+  /// \name Conversion into mdspan
+  /// @{
+
+  /// \brief Conversion operator to TensorSpan
+  template <class V, class E, class L, class A,
+    std::enable_if_t<std::is_assignable_v<TensorSpan<V,E,L,A>, tensorspan_type>, int> = 0>
+  constexpr operator TensorSpan<V,E,L,A> ()
+  {
+    return tensorspan_type(this->container_data(), this->mapping());
+  }
+
+  /// \brief Conversion operator to TensorSpan
+  template <class V, class E, class L, class A,
+    std::enable_if_t<std::is_assignable_v<TensorSpan<V,E,L,A>, const_tensorspan_type>, int> = 0>
+  constexpr operator TensorSpan<V,E,L,A> () const
+  {
+    return const_tensorspan_type(this->container_data(), this->mapping());
+  }
+
+  /// \brief Conversion function to TensorSpan
+  template <class A = Std::default_accessor<element_type>,
+    std::enable_if_t<
+      std::is_assignable_v<tensorspan_type, TensorSpan<element_type,extents_type,layout_type,A>>, int> = 0>
+  constexpr TensorSpan<element_type,extents_type,layout_type,A>
+  toTensorSpan (const A& accessor = A{})
+  {
+    return TensorSpan<element_type,extents_type,layout_type,A>(this->container_data(), this->mapping(), accessor);
+  }
+
+  /// \brief Conversion function to TensorSpan
+  template <class A = Std::default_accessor<const element_type>,
+    std::enable_if_t<
+      std::is_assignable_v<const_tensorspan_type, TensorSpan<const element_type,extents_type,layout_type,A>>, int> = 0>
+  constexpr TensorSpan<const element_type,extents_type,layout_type,A>
+  toTensorSpan (const A& accessor = A{}) const
+  {
+    return TensorSpan<const element_type,extents_type,layout_type,A>(this->container_data(), this->mapping(), accessor);
+  }
+
+  /// @}
+};
+
+namespace Impl {
+
+template <class M>
+using IsLayoutMapping = std::is_same<M,
+    typename M::layout_type::template mapping<typename M::extents_type>>;
+
+template <class V, class E>
+struct TensorFromExtents;
+
+template <class V, class I, std::size_t... exts>
+struct TensorFromExtents<V, Dune::Std::extents<I, exts...>>
+{
+  using type = Dune::Tensor<V, exts...>;
+};
+
+} // end namespace Impl
+
+
+/// \name Deduction guides
+/// \relates Tensor
+/// @{
+
+template <class index_type, std::size_t... exts, class value_type>
+Tensor (Std::extents<index_type, exts...>, value_type)
+  -> Tensor<value_type, exts...>;
+
+template <class Mapping, class value_type,
+  std::enable_if_t<Impl::IsLayoutMapping<Mapping>::value, int> = 0,
+  std::enable_if_t<(Mapping::extents_type::rank() == 0), int> = 0>
+Tensor (Mapping, value_type)
+  -> Tensor<value_type>;
+
+template <class Mapping, class value_type,
+  std::enable_if_t<Impl::IsLayoutMapping<Mapping>::value, int> = 0,
+  std::enable_if_t<(Mapping::extents_type::rank() == 1), int> = 0>
+Tensor (Mapping, value_type)
+  -> Tensor<value_type, Mapping::extents_type::static_extent(0)>;
+
+template <class Mapping, class value_type,
+  std::enable_if_t<Impl::IsLayoutMapping<Mapping>::value, int> = 0,
+  std::enable_if_t<(Mapping::extents_type::rank() == 2), int> = 0>
+Tensor (Mapping, value_type)
+  -> Tensor<value_type,
+      Mapping::extents_type::static_extent(0),
+      Mapping::extents_type::static_extent(1)>;
+
+template <class Mapping, class value_type,
+  std::enable_if_t<Impl::IsLayoutMapping<Mapping>::value, int> = 0,
+  std::enable_if_t<(Mapping::extents_type::rank() == 3), int> = 0>
+Tensor (Mapping, value_type)
+  -> Tensor<value_type,
+      Mapping::extents_type::static_extent(0),
+      Mapping::extents_type::static_extent(1),
+      Mapping::extents_type::static_extent(2)>;
+
+template <class index_type, std::size_t... exts, class value_type, class Alloc>
+Tensor (Std::extents<index_type, exts...>, value_type, const Alloc&)
+  -> Tensor<value_type, exts...>;
+
+template <class Mapping, class value_type, class Alloc,
+  std::enable_if_t<Impl::IsLayoutMapping<Mapping>::value, int> = 0,
+  std::enable_if_t<(Mapping::extents_type::rank() == 0), int> = 0>
+Tensor (Mapping, value_type, const Alloc&)
+  -> Tensor<value_type>;
+
+template <class Mapping, class value_type, class Alloc,
+  std::enable_if_t<Impl::IsLayoutMapping<Mapping>::value, int> = 0,
+  std::enable_if_t<(Mapping::extents_type::rank() == 1), int> = 0>
+Tensor (Mapping, value_type, const Alloc&)
+  -> Tensor<value_type, Mapping::extents_type::static_extent(0)>;
+
+template <class Mapping, class value_type, class Alloc,
+  std::enable_if_t<Impl::IsLayoutMapping<Mapping>::value, int> = 0,
+  std::enable_if_t<(Mapping::extents_type::rank() == 2), int> = 0>
+Tensor (Mapping, value_type, const Alloc&)
+  -> Tensor<value_type,
+      Mapping::extents_type::static_extent(0),
+      Mapping::extents_type::static_extent(1)>;
+
+template <class Mapping, class value_type, class Alloc,
+  std::enable_if_t<Impl::IsLayoutMapping<Mapping>::value, int> = 0,
+  std::enable_if_t<(Mapping::extents_type::rank() == 3), int> = 0>
+Tensor (Mapping, value_type, const Alloc&)
+  -> Tensor<value_type,
+      Mapping::extents_type::static_extent(0),
+      Mapping::extents_type::static_extent(1),
+      Mapping::extents_type::static_extent(2)>;
+
+template <class V, class I, std::size_t... exts, class L, class C>
+Tensor (Std::mdarray<V,Dune::Std::extents<I,exts...>,L,C>)
+  -> Tensor<V, exts...>;
+
+template <class V, class I, std::size_t... exts, class L, class C>
+Tensor (Std::mdspan<V,Dune::Std::extents<I,exts...>,L,C>)
+  -> Tensor<std::remove_cv_t<V>, exts...>;
+
+/// @}
+
+
+template <class V, std::size_t... exts>
+struct FieldTraits< Tensor<V,exts...> >
+{
+  using field_type = typename FieldTraits<V>::field_type;
+  using real_type = typename FieldTraits<V>::real_type;
+};
+
+} // end namespace Dune
+
+#endif // DUNE_COMMON_TENSOR_HH

dune/common/tensormixin.hh

+// -*- 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_TENSORMIXIN_HH
+#define DUNE_COMMON_TENSORMIXIN_HH
+
+#include <cassert>
+#include <stdexcept>
+#include <type_traits>
+#include <utility>
+
+#include <dune/common/ftraits.hh>
+#include <dune/common/indices.hh>
+#include <dune/common/std/type_traits.hh>
+
+namespace Dune {
+
+/**
+ * \brief A tensor interface-class providing common functionality.
+ * \nosubgrouping
+ *
+ * The CRTP interface mixin-class for tensors extends the mdarray and mdspan implementations
+ * by common functionality for element access and size information and mathematical
+ * operations.
+ *
+ * \tparam Derived  The tensor type derived from this class.
+ * \tparam Base     Either mdarray or mdspan base type for storage.
+ **/
+template <class Derived, class Base>
+class TensorMixin
+    : public Base
+{
+  using self_type = TensorMixin;
+  using derived_type = Derived;
+  using base_type = Base;
+
+  template <class B>
+  using const_reference_t = typename B::const_reference;
+
+public:
+  using element_type = typename base_type::element_type;
+  using value_type = std::remove_const_t<element_type>;
+  using extents_type = typename base_type::extents_type;
+  using index_type = typename base_type::index_type;
+  using layout_type = typename base_type::layout_type;
+  using mapping_type = typename base_type::mapping_type;
+  using reference = typename base_type::reference;
+  using const_reference = Std::detected_or_t<reference,const_reference_t,base_type>;
+
+protected:
+  /// Derive constructors from base class
+  using base_type::base_type;
+
+
+  /// \brief base copy constructor
+  constexpr TensorMixin (const base_type& tensor)
+    : base_type{tensor}
+  {}
+
+  /// \brief base move constructor
+  constexpr TensorMixin (base_type&& tensor)
+    : base_type{std::move(tensor)}
+  {}
+
+public:
+  /// \name Multi index access
+  /// @{
+
+  /// \brief Access specified element at position (i0,i1,...) with mutable access
+  template <class... Indices,
+    std::enable_if_t<(sizeof...(Indices) == extents_type::rank()), int> = 0,
+    std::enable_if_t<(... && std::is_convertible_v<Indices, index_type>), int> = 0>
+  constexpr reference operator() (Indices... indices)
+  {
+    return base_type::operator[](std::array<index_type,extents_type::rank()>{index_type(indices)...});
+  }
+
+  /// \brief Access specified element at position (i0,i1,...) with const access
+  template <class... Indices,
+    std::enable_if_t<(sizeof...(Indices) == extents_type::rank()), int> = 0,
+    std::enable_if_t<(... && std::is_convertible_v<Indices, index_type>), int> = 0>
+  constexpr const_reference operator() (Indices... indices) const
+  {
+    return base_type::operator[](std::array<index_type,extents_type::rank()>{index_type(indices)...});
+  }
+
+  /**
+   * \brief Access specified element at position (i0,i1,...) with mutable access
+   * \throws std::out_of_range if the indices are out of the index space `[0,extent_0)x...x[0,extent_{r-1})`.
+   */
+  template <class... Indices,
+    std::enable_if_t<(sizeof...(Indices) == extents_type::rank()), int> = 0,
+    std::enable_if_t<(... && std::is_convertible_v<Indices, index_type>), int> = 0>
+  constexpr reference at (Indices... indices)
+  {
+    if (not indexInIndexSpace(indices...))
+      throw std::out_of_range("Indices out of index space");
+    return base_type::operator[](std::array<index_type,extents_type::rank()>{index_type(indices)...});
+  }
+
+  /**
+   * \brief Access specified element at position (i0,i1,...) with const access
+   * \throws std::out_of_range if the indices are out of the index space `[0,extent_0)x...x[0,extent_{r-1})`.
+   */
+  template <class... Indices,
+    std::enable_if_t<(sizeof...(Indices) == extents_type::rank()), int> = 0,
+    std::enable_if_t<(... && std::is_convertible_v<Indices, index_type>), int> = 0>
+  constexpr const_reference at (Indices... indices) const
+  {
+    if (not indexInIndexSpace(indices...))
+      throw std::out_of_range("Indices out of index space");
+    return base_type::operator[](std::array<index_type,extents_type::rank()>{index_type(indices)...});
+  }
+
+  /// \brief Access element at position [{i0,i1,...}]
+  template <class Index,
+    std::enable_if_t<std::is_convertible_v<const Index&, index_type>, int> = 0>
+  constexpr reference operator[] (const std::array<Index,extents_type::rank()>& indices)
+  {
+    return base_type::operator[](indices);
+  }
+
+  /// \brief Access element at position [{i0,i1,...}]
+  template <class Index,
+    std::enable_if_t<std::is_convertible_v<const Index&, index_type>, int> = 0>
+  constexpr const_reference operator[] (const std::array<Index,extents_type::rank()>& indices) const
+  {
+    return base_type::operator[](indices);
+  }
+
+  /// \brief Access vector-element at position [i0] with mutable access.
+  template <class E = extents_type,
+    std::enable_if_t<(E::rank() == 1), int> = 0>
+  constexpr reference operator[] (index_type index) noexcept
+  {
+    return base_type::operator[](std::array{index});
+  }
+
+  /// \brief Access vector-element at position [i0] with const access.
+  template <class E = extents_type,
+    std::enable_if_t<(E::rank() == 1), int> = 0>
+  constexpr const_reference operator[] (index_type index) const noexcept
+  {
+    return base_type::operator[](std::array{index});
+  }
+
+  /**
+   * \brief Return true when (i0,i1,...) is in pattern.
+   * This is always true for dense tensors.
+   **/
+  template <class... Indices,
+    std::enable_if_t<(sizeof...(Indices) == extents_type::rank()), int> = 0,
+    std::enable_if_t<(... && std::is_convertible_v<Indices, index_type>), int> = 0>
+  constexpr bool exists (Indices... indices) const noexcept
+  {
+    return indexInIndexSpace(indices...);
+  }
+
+  /// @}
+
+
+  /// \name Size information
+  /// @{
+
+  /// \brief Number of rows of a rank-2 tensor
+  template <class E = extents_type,
+    std::enable_if_t<(E::rank() == 2), int> = 0>
+  constexpr index_type rows () const noexcept
+  {
+    return asBase().extent(0);
+  }
+
+  /// \brief Number of columns of a rank-2 tensor
+  template <class E = extents_type,
+    std::enable_if_t<(E::rank() == 2), int> = 0>
+  constexpr index_type cols () const noexcept
+  {
+    return asBase().extent(1);
+  }
+
+  /// @}
+
+
+  /// \name Conversion to the underlying value if rank is zero
+  // @{
+
+  template <class ScalarType, class E = extents_type,
+    std::enable_if_t<std::is_convertible_v<value_type,ScalarType>, int> = 0,
+    std::enable_if_t<(E::rank() == 0), int> = 0>
+  constexpr operator ScalarType () const noexcept
+  {
+    return ScalarType(base_type::operator[](std::array<index_type,0>{}));
+  }
+
+  /// @}
+
+
+  /// \brief Comparison of two TensorMixins for equality
+  friend constexpr bool operator== (const TensorMixin& lhs, const TensorMixin& rhs) noexcept
+  {
+    return static_cast<const base_type&>(lhs) == static_cast<const base_type&>(rhs);
+  }
+
+private:
+  // Check whether a tuple of indices is in the index space `[0,extent_0)x...x[0,extent_{r-1})`.
+  template <class... Indices>
+  [[nodiscard]] constexpr bool indexInIndexSpace (Indices... indices) const noexcept
+  {
+    return unpackIntegerSequence([&](auto... i) {
+      return ( (0 <= indices && index_type(indices) < asBase().extent(i)) && ... );
+    }, std::make_index_sequence<sizeof...(Indices)>{});
+  }
+
+private:
+
+  derived_type const& asDerived () const
+  {
+    return static_cast<derived_type const&>(*this);
+  }
+
+  derived_type& asDerived ()
+  {
+    return static_cast<derived_type&>(*this);
+  }
+
+  base_type const& asBase () const
+  {
+    return static_cast<base_type const&>(*this);
+  }
+
+  base_type& asBase ()
+  {
+    return static_cast<base_type&>(*this);
+  }
+};
+
+// specialization for rank-0 tensor and comparison with scalar
+template <class D, class B, class V, class E = typename B::extents_type,
+  std::enable_if_t<(E::rank() == 0), int> = 0,
+  std::enable_if_t<std::is_convertible_v<V, typename B::value_type>, int> = 0,
+  std::enable_if_t<std::is_constructible_v<typename B::value_type, V>, int> = 0>
+constexpr bool operator== (const TensorMixin<D,B>& lhs, const V& rhs) noexcept
+{
+  return lhs() == rhs;
+}
+
+// specialization for rank-0 tensor and comparison with scalar
+template <class V, class D, class B, class E = typename B::extents_type,
+  std::enable_if_t<(E::rank() == 0), int> = 0,
+  std::enable_if_t<std::is_convertible_v<V, typename B::value_type>, int> = 0,
+  std::enable_if_t<std::is_constructible_v<typename B::value_type, V>, int> = 0>
+constexpr bool operator== (const V& lhs, const TensorMixin<D,B>& rhs) noexcept
+{
+  return lhs == rhs();
+}
+
+template <class D, class B>
+struct FieldTraits< TensorMixin<D,B> >
+{
+  using field_type = typename FieldTraits<D>::field_type;
+  using real_type = typename FieldTraits<D>::real_type;
+};
+
+} // end namespace Dune
+
+#endif // DUNE_COMMON_TENSORMIXIN_HH

dune/common/tensorspan.hh

+// -*- 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_TENSORSPAN_HH
+#define DUNE_COMMON_TENSORSPAN_HH
+
+#include <array>
+#include <type_traits>
+
+#include <dune/common/ftraits.hh>
+#include <dune/common/tensormixin.hh>
+#include <dune/common/std/default_accessor.hh>
+#include <dune/common/std/extents.hh>
+#include <dune/common/std/layout_right.hh>
+#include <dune/common/std/mdspan.hh>
+
+namespace Dune {
+
+/**
+ * \brief A multi-dimensional span with tensor interface.
+ * \nosubgrouping
+ *
+ * \tparam Element  The element type stored in the tensor.
+ * \tparam Extents  The type representing the tensor extents. Must be compatible
+ *                  to `Dune::Std::extents`.
+ * \tparam Layout   A class providing the index mapping from the multi-dimensional
+ *                  index space to a single index usable in the `Container`. See
+ *                  the `Dune::Concept::Layout` concept.
+ * \tparam Accessor A class that defines types and operations to create a reference
+ *                  to a single element stored by a data pointer and index.
+ **/
+template <class Element, class Extents,
+          class Layout = Std::layout_right,
+          class Accessor= Std::default_accessor<Element>>
+class TensorSpan
+    : public TensorMixin<TensorSpan<Element,Extents,Layout,Accessor>,
+        Std::mdspan<Element,Extents,Layout,Accessor>>
+{
+  using self_type = TensorSpan;
+  using storage_type = Std::mdspan<Element,Extents,Layout,Accessor>;
+  using base_type = TensorMixin<self_type,storage_type>;
+
+public:
+  using element_type =	Element;
+  using extents_type = Extents;
+  using layout_type = Layout;
+  using accessor_type = Accessor;
+  using mapping_type = typename base_type::mapping_type;
+
+public:
+  /// \name TensorSpan constructors
+  /// @{
+
+  using base_type::base_type;
+
+  /// \brief Converting constructor
+  template <class V, class E, class L, class A,
+    class M = typename L::template mapping<E>,
+    std::enable_if_t<std::is_constructible_v<mapping_type, const M&>, int> = 0,
+    std::enable_if_t<std::is_constructible_v<accessor_type, const A&>, int> = 0>
+  constexpr TensorSpan (const TensorSpan<V,E,L,A>& other)
+    : base_type{other}
+  {}
+
+  /// \brief Converting move constructor
+  template <class V, class E, class L, class A,
+    class M = typename L::template mapping<E>,
+    std::enable_if_t<std::is_constructible_v<mapping_type, const M&>, int> = 0,
+    std::enable_if_t<std::is_constructible_v<accessor_type, const A&>, int> = 0>
+  constexpr TensorSpan (TensorSpan<V,E,L,A>&& tensor)
+    : base_type{std::move(tensor)}
+  {}
+
+  /// \brief base copy constructor
+  constexpr TensorSpan (const storage_type& tensor)
+    : base_type{tensor}
+  {}
+
+  /// \brief base move constructor
+  constexpr TensorSpan (storage_type&& tensor)
+    : base_type{std::move(tensor)}
+  {}
+
+  /// @}
+};
+
+// deduction guides
+// @{
+
+template <class CArray,
+  std::enable_if_t<std::is_array_v<CArray>, int> = 0,
+  std::enable_if_t<(std::rank_v<CArray> == 1), int> = 0>
+TensorSpan (CArray&)
+  -> TensorSpan<std::remove_all_extents_t<CArray>, Std::extents<std::size_t, std::extent_v<CArray,0>>>;
+
+template <class Pointer,
+  std::enable_if_t<std::is_pointer_v<std::remove_reference_t<Pointer>>, int> = 0>
+TensorSpan (Pointer&&)
+  -> TensorSpan<std::remove_pointer_t<std::remove_reference_t<Pointer>>, Std::extents<std::size_t>>;
+
+template <class element_type, class... II,
+  std::enable_if_t<(... && std::is_convertible_v<II,std::size_t>), int> = 0,
+  std::enable_if_t<(sizeof...(II) > 0), int> = 0>
+TensorSpan (element_type*, II...)
+  -> TensorSpan<element_type, Std::dextents<std::size_t, sizeof...(II)>>;
+
+template <class element_type, class SizeType, std::size_t N>
+TensorSpan (element_type*, Std::span<SizeType,N>&)
+  -> TensorSpan<element_type, Std::dextents<std::size_t, N>>;
+
+template <class element_type, class SizeType, std::size_t N>
+TensorSpan (element_type*, const std::array<SizeType,N>&)
+  -> TensorSpan<element_type, Std::dextents<std::size_t, N>>;
+
+template <class element_type, class IndexType, std::size_t... exts>
+TensorSpan (element_type*, const Std::extents<IndexType,exts...>&)
+  -> TensorSpan<element_type, Std::extents<IndexType,exts...>>;
+
+template <class element_type, class Mapping,
+  class Extents = typename Mapping::extents_type,
+  class Layout = typename Mapping::layout_type>
+TensorSpan (element_type*, const Mapping&)
+  -> TensorSpan<element_type, Extents, Layout>;
+
+template <class Mapping, class Accessor,
+  class DataHandle = typename Accessor::data_handle_type,
+  class Element = typename Accessor::element_type,
+  class Extents = typename Mapping::extents_type,
+  class Layout = typename Mapping::layout_type>
+TensorSpan (const DataHandle&, const Mapping&, const Accessor&)
+  -> TensorSpan<Element, Extents, Layout, Accessor>;
+
+template <class V, class E, class L, class A>
+TensorSpan (Std::mdspan<V,E,L,A>)
+  -> TensorSpan<V,E,L,A>;
+
+/// @}
+
+
+template <class Element, class Extents, class Layout, class Accessor>
+struct FieldTraits< TensorSpan<Element,Extents,Layout,Accessor> >
+{
+  using field_type = typename FieldTraits<Element>::field_type;
+  using real_type = typename FieldTraits<Element>::real_type;
+};
+
+} // end namespace Dune
+
+#endif // DUNE_COMMON_TENSORSPAN_HH

dune/common/test/CMakeLists.txt

@@ -127,6 +127,10 @@ dune_add_test(SOURCES dunethrowtest.cc
 dune_add_test(SOURCES dynmatrixtest.cc
               LABELS quick)
 
+dune_add_test(SOURCES dyntensortest.cc
+              LINK_LIBRARIES dunecommon
+              LABELS quick)
+
 dune_add_test(SOURCES dynvectortest.cc
               LABELS quick)
 
@@ -146,6 +150,9 @@ add_dune_vc_flags(fmatrixtest)
 dune_add_test(SOURCES forceinlinetest.cc
               LABELS quick)
 
+dune_add_test(SOURCES ftensortest.cc
+              LABELS quick)
+
 dune_add_test(SOURCES fvectortest.cc
               LABELS quick)
 

dune/common/test/dyntensortest.cc

+// -*- 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
+
+#include <array>
+
+#include <dune/common/filledarray.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/tensor.hh>
+#include <dune/common/tensorspan.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/test/testsuite.hh>
+#include <dune/common/test/foreachindex.hh>
+#include <type_traits>
+
+using namespace Dune;
+
+template <class Tensor1, class Tensor2>
+void checkEqual(Dune::TestSuite& testSuite, Tensor1 const& a, Tensor2 const& b)
+{
+  Dune::TestSuite subTestSuite("checkEqual");
+  subTestSuite.require(a.extents() == b.extents());
+  Dune::forEachIndex(a.extents(), [&](auto index) {
+    subTestSuite.check(a[index] == b[index]);
+  });
+  testSuite.subTest(subTestSuite);
+}
+
+template <class Tensor>
+void checkEqualValue(Dune::TestSuite& testSuite, Tensor const& a, typename Tensor::value_type const& value)
+{
+  Dune::TestSuite subTestSuite("checkEqualValue");
+  Dune::forEachIndex(a.extents(), [&](auto index) {
+    subTestSuite.check(a[index] == value);
+  });
+  testSuite.subTest(subTestSuite);
+}
+
+
+template <class Tensor>
+void checkConstructors(Dune::TestSuite& testSuite)
+{
+  Dune::TestSuite subTestSuite("checkConstructors");
+
+  using T = typename Tensor::element_type;
+  using extents_type = typename Tensor::extents_type;
+  extents_type ext{Dune::filledArray<Tensor::rank()>(3)};
+
+  // default constructor
+  Tensor tensor0;
+  subTestSuite.check(tensor0.size() == 0 || tensor0.rank() == 0);
+
+  Tensor tensor(ext);
+  checkEqualValue(subTestSuite, tensor, T(0));
+
+  if constexpr(Tensor::rank() > 0) {
+    tensor0.resize(ext, T(0));
+    subTestSuite.check(tensor0.size() == tensor.size());
+  }
+
+#if 0 // not yet implemented
+  tensor = 1.0;
+  checkEqualValue(subTestSuite, tensor, 1.0);
+#endif
+
+  // constructor with a default value
+  Tensor tensor1(ext, T(1));
+  checkEqualValue(subTestSuite, tensor1, T(1));
+
+  // copy/move constructors
+  Tensor tensor2{tensor1};
+  checkEqual(subTestSuite, tensor1,tensor2);
+
+  Tensor tensor3 = tensor1;
+  checkEqual(subTestSuite, tensor1,tensor3);
+
+  Tensor tensor4{std::move(tensor2)};
+  checkEqual(subTestSuite, tensor1,tensor4);
+
+  Tensor tensor5 = std::move(tensor3);
+  checkEqual(subTestSuite, tensor1,tensor5);
+
+  // copy/move assignment operators
+  tensor4 = tensor1;
+  checkEqual(subTestSuite, tensor1,tensor4);
+
+  tensor5 = std::move(tensor4);
+  checkEqual(subTestSuite, tensor1,tensor5);
+
+  extents_type ext2{Dune::filledArray<Tensor::rank()>(2)};
+  if constexpr(Tensor::rank() == 1) {
+    if constexpr (std::is_floating_point_v<T>) {
+      Tensor tensor6(ext2, {6,7});
+      subTestSuite.check(tensor6(0) == T(6.0));
+      subTestSuite.check(tensor6(1) == T(7.0));
+    }
+
+    Tensor tensor7{ext2, {T(6.0),T(7.0)}};
+    subTestSuite.check(tensor7(0) == T(6.0));
+    subTestSuite.check(tensor7(1) == T(7.0));
+  }
+  else
+  if constexpr(Tensor::rank() == 2) {
+    if constexpr (std::is_floating_point_v<T>) {
+      Tensor tensor6(ext2, {{6,7},{8,9}});
+      subTestSuite.check(tensor6(0,0) == T(6.0));
+      subTestSuite.check(tensor6(0,1) == T(7.0));
+      subTestSuite.check(tensor6(1,0) == T(8.0));
+      subTestSuite.check(tensor6(1,1) == T(9.0));
+    }
+
+    Tensor tensor7(ext2, {{T(6.0),T(7.0)},{T(8.0),T(9.0)}});
+    Tensor tensor8a(std::array<int,2>{2,2}, {{T(6.0),T(7.0)},{T(8.0),T(9.0)}});
+    Tensor tensor8b(std::array<unsigned int,2>{2u,2u}, {{T(6.0),T(7.0)},{T(8.0),T(9.0)}});
+
+    subTestSuite.check(tensor7(0,0) == T(6.0));
+    subTestSuite.check(tensor7(0,1) == T(7.0));
+    subTestSuite.check(tensor7(1,0) == T(8.0));
+    subTestSuite.check(tensor7(1,1) == T(9.0));
+  }
+
+
+  Tensor tensor9(tensor1.toTensorSpan());
+
+  { // check deduction guides
+
+    Dune::Tensor t1{tensor1.extents(), T{}};
+    Dune::Tensor t2{tensor1.mapping(), T{}};
+    Dune::Tensor t3{tensor1.toTensorSpan()};
+  }
+
+  testSuite.subTest(subTestSuite);
+}
+
+
+template <class Tensor>
+void checkAccess(Dune::TestSuite& testSuite)
+{
+  Dune::TestSuite subTestSuite("checkAccess");
+
+  using T = typename Tensor::element_type;
+  using extents_type = typename Tensor::extents_type;
+  using index_type = typename Tensor::index_type;
+  extents_type ext{Dune::filledArray<Tensor::rank()>(3)};
+
+  Tensor tensor(ext, T(42.0));
+  checkEqualValue(subTestSuite, tensor, T(42.0));
+
+  if constexpr(Tensor::rank() == 0) {
+    subTestSuite.check(tensor[std::array<std::size_t,0>{}] == T(42.0));
+    subTestSuite.check(tensor() == T(42.0));
+    subTestSuite.check(tensor.at() == T(42.0));
+    subTestSuite.check(tensor == T(42.0));
+    double value = tensor;
+    subTestSuite.check(value == T(42.0));
+  }
+  else if constexpr(Tensor::rank() == 1) {
+    for (index_type i = 0; i < tensor.extent(0); ++i) {
+      subTestSuite.check(tensor[std::array{i}] == T(42.0));
+      subTestSuite.check(tensor[i] == T(42.0));
+      subTestSuite.check(tensor(i) == T(42.0));
+      subTestSuite.check(tensor.at(i) == T(42.0));
+    }
+  }
+  else if constexpr(Tensor::rank() == 2) {
+    for (index_type i = 0; i < tensor.extent(0); ++i) {
+      for (index_type j = 0; j < tensor.extent(1); ++j) {
+        subTestSuite.check(tensor[std::array{i,j}] == T(42.0));
+        subTestSuite.check(tensor(i,j) == T(42.0));
+        subTestSuite.check(tensor.at(i,j) == T(42.0));
+      }
+    }
+  }
+  else if constexpr(Tensor::rank() == 3) {
+    for (index_type i = 0; i < tensor.extent(0); ++i) {
+      for (index_type j = 0; j < tensor.extent(1); ++j) {
+        for (index_type k = 0; k < tensor.extent(2); ++k) {
+          subTestSuite.check(tensor[std::array{i,j,k}] == T(42.0));
+          subTestSuite.check(tensor(i,j,k) == T(42.0));
+          subTestSuite.check(tensor.at(i,j,k) == T(42.0));
+        }
+      }
+    }
+  }
+
+  testSuite.subTest(subTestSuite);
+}
+
+
+template <class Tensor>
+void checkArithmetic(Dune::TestSuite& testSuite)
+{
+  Dune::TestSuite subTestSuite("checkArithmetic");
+
+  using T = typename Tensor::element_type;
+  using extents_type = typename Tensor::extents_type;
+  extents_type ext{Dune::filledArray<Tensor::rank()>(3)};
+
+  Tensor tensor(ext, T(1));
+  Tensor tensor2(ext, T(2));
+  checkEqualValue(subTestSuite, tensor, T(1));
+  checkEqualValue(subTestSuite, tensor2, T(2));
+
+#if 0 // not yet implemented
+  if constexpr(std::is_floating_point_v<T>) {
+    tensor *= 2.0;
+    checkEqualValue(subTestSuite, tensor, T(2));
+
+    tensor += tensor2;
+    checkEqualValue(subTestSuite, tensor, T(4));
+
+    tensor.axpy(4.0, tensor2); // 12
+    checkEqualValue(subTestSuite, tensor, T(12));
+
+    tensor.aypx(4.0, tensor2); // 50
+    checkEqualValue(subTestSuite, tensor, T(50.0));
+
+    tensor -= tensor2;
+    checkEqualValue(subTestSuite, tensor, T(48.0));
+  }
+#endif
+
+  testSuite.subTest(subTestSuite);
+}
+
+
+
+int main(int argc, char** argv)
+{
+  MPIHelper::instance(argc, argv);
+
+  TestSuite testSuite;
+
+  using Tensor0 = Dune::Tensor<double>;
+  using Tensor1 = Dune::Tensor<double,Dune::dynamic>;
+  using Tensor2 = Dune::Tensor<double,Dune::dynamic,Dune::dynamic>;
+  using Tensor3 = Dune::Tensor<double,Dune::dynamic,Dune::dynamic,Dune::dynamic>;
+
+  // using Tensor4 = Dune::Tensor<bool,Dune::dynamic>;
+  // using Tensor5 = Dune::Tensor<bool,Dune::dynamic,Dune::dynamic>;
+
+  checkConstructors<Tensor0>(testSuite);
+  checkConstructors<Tensor1>(testSuite);
+  checkConstructors<Tensor2>(testSuite);
+  checkConstructors<Tensor3>(testSuite);
+  // checkConstructors<Tensor4>(testSuite);
+  // checkConstructors<Tensor5>(testSuite);
+
+  checkAccess<Tensor0>(testSuite);
+  checkAccess<Tensor1>(testSuite);
+  checkAccess<Tensor2>(testSuite);
+  checkAccess<Tensor3>(testSuite);
+  // checkAccess<Tensor4>(testSuite);
+  // checkAccess<Tensor5>(testSuite);
+
+  checkArithmetic<Tensor0>(testSuite);
+  checkArithmetic<Tensor1>(testSuite);
+  checkArithmetic<Tensor2>(testSuite);
+  checkArithmetic<Tensor3>(testSuite);
+  // checkArithmetic<Tensor4>(testSuite);
+  // checkArithmetic<Tensor5>(testSuite);
+
+  return testSuite.exit();
+}

dune/common/test/foreachindex.hh

+// -*- 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_TEST_FOREACHINDEX_HH
+#define DUNE_COMMON_TEST_FOREACHINDEX_HH
+
+#include <array>
+#include <functional>
+#include <type_traits>
+#include <utility>
+
+#include <dune/common/std/extents.hh>
+#include <dune/common/std/span.hh>
+#include <dune/common/std/functional.hh>
+
+namespace Dune {
+namespace Impl {
+
+template <class Extents, class Fun, class... Indices>
+inline __attribute__((always_inline))
+void forEachIndexImpl (const Extents& extents, Fun f, Indices... ii)
+{
+  constexpr typename Extents::rank_type pos = sizeof...(Indices);
+  if constexpr(pos < Extents::rank()) {
+    if constexpr(Extents::static_extent(pos) == Std::dynamic_extent)
+      for (typename Extents::index_type i = 0; i < extents.extent(pos); ++i)
+        forEachIndexImpl(extents, std::move(f), ii...,i);
+    else
+      for (std::size_t i = 0; i < Extents::static_extent(pos); ++i)
+        forEachIndexImpl(extents, std::move(f), ii...,typename Extents::index_type(i));
+  }
+  else {
+    using I = std::array<typename Extents::index_type,Extents::rank()>;
+    std::invoke(f, I{ii...});
+  }
+}
+
+} // end namespace Impl
+
+
+/// \brief Invoke the function `f` on all index-tuples in the multi dimensional index-space given by `extents`.
+template <class Extents, class Fun,
+  class I = std::array<typename Extents::index_type,Extents::rank()>,
+  std::enable_if_t<std::is_invocable_v<Fun,I>, int> = 0>
+inline __attribute__((always_inline))
+void forEachIndex (const Extents& extents, Fun f)
+{
+  Impl::forEachIndexImpl(extents, std::move(f));
+}
+
+} // end namespace Dune
+
+#endif // DUNE_COMMON_TEST_FOREACHINDEX_HH

dune/common/test/ftensortest.cc

+// -*- 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
+
+#include <iostream>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/tensor.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/test/foreachindex.hh>
+#include <dune/common/test/testsuite.hh>
+#include <stdexcept>
+
+using namespace Dune;
+
+
+template <class Tensor1, class Tensor2>
+void checkEqual(Dune::TestSuite& testSuite, Tensor1 const& a, Tensor2 const& b)
+{
+  Dune::TestSuite subTestSuite("checkEqual");
+  subTestSuite.require(a.extents() == b.extents());
+  Dune::forEachIndex(a.extents(), [&](auto index) {
+    subTestSuite.check(a[index] == b[index]);
+  });
+  testSuite.subTest(subTestSuite);
+}
+
+template <class Tensor>
+void checkEqualValue(Dune::TestSuite& testSuite, Tensor const& a, typename Tensor::value_type const& value)
+{
+  Dune::TestSuite subTestSuite("checkEqualValue");
+  Dune::forEachIndex(a.extents(), [&](auto index) {
+    subTestSuite.check(a[index] == value);
+  });
+  testSuite.subTest(subTestSuite);
+}
+
+template <class D, class B>
+void call (const TensorMixin<D,B>& tensor)
+{
+  // it is possible to call a function that expects the TensorInterface base class
+}
+
+template <class Tensor>
+void checkConstructors(Dune::TestSuite& testSuite)
+{
+  Dune::TestSuite subTestSuite("checkConstructors");
+
+  // default constructor
+  std::cout << "default constructor..." << std::endl;
+  Tensor tensor;
+  checkEqualValue(subTestSuite, tensor, 0.0);
+
+#if 0 // not yet implemented
+  std::cout << "assignment of value..." << std::endl;
+  tensor = 1.0;
+  checkEqualValue(subTestSuite, tensor, 1.0);
+#endif
+
+  // constructor with a default value
+  std::cout << "value constructor..." << std::endl;
+  Tensor tensor1(1.0);
+  checkEqualValue(subTestSuite, tensor1, 1.0);
+
+  // copy/move constructors
+  std::cout << "copy constructor..." << std::endl;
+  Tensor tensor2{tensor1};
+  checkEqual(subTestSuite, tensor1,tensor2);
+
+  Tensor tensor3 = tensor1;
+  checkEqual(subTestSuite, tensor1,tensor3);
+
+  std::cout << "move constructor..." << std::endl;
+  Tensor tensor4{std::move(tensor2)};
+  checkEqual(subTestSuite, tensor1,tensor4);
+
+  Tensor tensor5 = std::move(tensor3);
+  checkEqual(subTestSuite, tensor1,tensor5);
+
+  // copy/move assignment operators
+  std::cout << "copy assignment..." << std::endl;
+  tensor4 = tensor1;
+  checkEqual(subTestSuite, tensor1,tensor4);
+
+  std::cout << "move assignment..." << std::endl;
+  tensor5 = std::move(tensor4);
+  checkEqual(subTestSuite, tensor1,tensor5);
+
+  std::cout << "initializerlist constructor..." << std::endl;
+  if constexpr(Tensor::rank() == 1) {
+    Tensor tensor6{6,7};
+    Tensor tensor7{6.0,7.0};
+    Tensor tensor8({6.0,7.0});
+    Tensor tensor9{{6.0,7.0}};
+    Tensor tensor10 = {6.0,7.0};
+
+    subTestSuite.check(tensor6(0) == 6.0);
+    subTestSuite.check(tensor6(1) == 7.0);
+    subTestSuite.check(tensor7(0) == 6.0);
+    subTestSuite.check(tensor7(1) == 7.0);
+    subTestSuite.check(tensor8(0) == 6.0);
+    subTestSuite.check(tensor8(1) == 7.0);
+    subTestSuite.check(tensor9(0) == 6.0);
+    subTestSuite.check(tensor9(1) == 7.0);
+    subTestSuite.check(tensor10(0) == 6.0);
+    subTestSuite.check(tensor10(1) == 7.0);
+  }
+  else
+  if constexpr(Tensor::rank() == 2) {
+    Tensor tensor6{{6,7},{8,9}};
+    Tensor tensor7{{6.0,7.0},{8.0,9.0}};
+    Tensor tensor8({{6.0,7.0},{8.0,9.0}});
+    Tensor tensor9{{{6.0,7.0},{8.0,9.0}}};
+    Tensor tensor10 = {{6.0,7.0},{8.0,9.0}};
+
+    subTestSuite.check(tensor6(0,0) == 6.0);
+    subTestSuite.check(tensor6(0,1) == 7.0);
+    subTestSuite.check(tensor6(1,0) == 8.0);
+    subTestSuite.check(tensor6(1,1) == 9.0);
+    subTestSuite.check(tensor7(0,0) == 6.0);
+    subTestSuite.check(tensor7(0,1) == 7.0);
+    subTestSuite.check(tensor7(1,0) == 8.0);
+    subTestSuite.check(tensor7(1,1) == 9.0);
+    subTestSuite.check(tensor8(0,0) == 6.0);
+    subTestSuite.check(tensor8(0,1) == 7.0);
+    subTestSuite.check(tensor8(1,0) == 8.0);
+    subTestSuite.check(tensor8(1,1) == 9.0);
+    subTestSuite.check(tensor9(0,0) == 6.0);
+    subTestSuite.check(tensor9(0,1) == 7.0);
+    subTestSuite.check(tensor9(1,0) == 8.0);
+    subTestSuite.check(tensor9(1,1) == 9.0);
+    subTestSuite.check(tensor10(0,0) == 6.0);
+    subTestSuite.check(tensor10(0,1) == 7.0);
+    subTestSuite.check(tensor10(1,0) == 8.0);
+    subTestSuite.check(tensor10(1,1) == 9.0);
+  }
+
+  call(tensor);
+  testSuite.subTest(subTestSuite);
+}
+
+
+template <class Tensor>
+void checkAccess(Dune::TestSuite& testSuite)
+{
+  Dune::TestSuite subTestSuite("checkAccess");
+
+  Tensor tensor(42.0);
+  checkEqualValue(subTestSuite, tensor, 42.0);
+
+  if constexpr(Tensor::rank() == 0) {
+    subTestSuite.check(tensor[std::array<int,0>{}] == 42.0);
+    subTestSuite.check(tensor() == 42.0);
+    subTestSuite.check(tensor == 42.0);
+    double value = tensor;
+    subTestSuite.check(value == 42.0);
+  }
+  else if constexpr(Tensor::rank() == 1) {
+    for (std::size_t i = 0; i < Tensor::static_extent(0); ++i) {
+      subTestSuite.check(tensor[std::array{i}] == 42.0);
+      subTestSuite.check(tensor(i) == 42.0);
+      subTestSuite.check(tensor[i] == 42.0);
+      subTestSuite.check(tensor.at(i) == 42.0);
+    }
+    subTestSuite.checkThrow<std::out_of_range>([&]{tensor.at(-1);});
+    subTestSuite.checkThrow<std::out_of_range>([&]{tensor.at(tensor.extent(0));});
+  }
+  else if constexpr(Tensor::rank() == 2) {
+    for (std::size_t i = 0; i < Tensor::static_extent(0); ++i) {
+      for (std::size_t j = 0; j < Tensor::static_extent(1); ++j) {
+        subTestSuite.check(tensor[std::array{i,j}] == 42.0);
+        subTestSuite.check(tensor(i,j) == 42.0);
+        subTestSuite.check(tensor.at(i,j) == 42.0);
+      }
+      subTestSuite.checkThrow<std::out_of_range>([&]{tensor.at(-1,-2);});
+      subTestSuite.checkThrow<std::out_of_range>([&]{tensor.at(tensor.extent(0),tensor.extent(1));});
+    }
+  }
+  else if constexpr(Tensor::rank() == 3) {
+    for (std::size_t i = 0; i < Tensor::static_extent(0); ++i) {
+      for (std::size_t j = 0; j < Tensor::static_extent(1); ++j) {
+        for (std::size_t k = 0; k < Tensor::static_extent(2); ++k) {
+          subTestSuite.check(tensor[std::array{i,j,k}] == 42.0);
+          subTestSuite.check(tensor(i,j,k) == 42.0);
+          subTestSuite.check(tensor.at(i,j,k) == 42.0);
+        }
+        subTestSuite.checkThrow<std::out_of_range>([&]{tensor.at(-1,-2,-3);});
+        subTestSuite.checkThrow<std::out_of_range>([&]{tensor.at(tensor.extent(0),tensor.extent(1),tensor.extent(2));});
+      }
+    }
+  }
+
+  testSuite.subTest(subTestSuite);
+}
+
+
+template <class Tensor>
+void checkArithmetic(Dune::TestSuite& testSuite)
+{
+  Dune::TestSuite subTestSuite("checkArithmetic");
+
+  Tensor tensor(1.0);
+  Tensor tensor2(2.0);
+  checkEqualValue(subTestSuite, tensor, 1.0);
+  checkEqualValue(subTestSuite, tensor2, 2.0);
+
+#if 0 // not yet implemented
+  tensor *= 2.0;
+  checkEqualValue(subTestSuite, tensor, 2.0);
+
+  tensor += tensor2;
+  checkEqualValue(subTestSuite, tensor, 4.0);
+
+  tensor.axpy(4.0, tensor2); // 12
+  checkEqualValue(subTestSuite, tensor, 12.0);
+
+  tensor.aypx(4.0, tensor2); // 50
+  checkEqualValue(subTestSuite, tensor, 50.0);
+
+  tensor -= tensor2;
+  checkEqualValue(subTestSuite, tensor, 48.0);
+#endif
+
+  testSuite.subTest(subTestSuite);
+}
+
+
+
+int main(int argc, char** argv)
+{
+  MPIHelper::instance(argc, argv);
+
+  TestSuite testSuite;
+
+  using Tensor0 = Dune::Tensor<double>;
+  using Tensor1 = Dune::Tensor<double,2>;
+  using Tensor2 = Dune::Tensor<double,2,2>;
+  using Tensor3 = Dune::Tensor<double,2,2,2>;
+
+  checkConstructors<Tensor0>(testSuite);
+  checkConstructors<Tensor1>(testSuite);
+  checkConstructors<Tensor2>(testSuite);
+  checkConstructors<Tensor3>(testSuite);
+
+  checkAccess<Tensor0>(testSuite);
+  checkAccess<Tensor1>(testSuite);
+  checkAccess<Tensor2>(testSuite);
+  checkAccess<Tensor3>(testSuite);
+
+  checkArithmetic<Tensor0>(testSuite);
+  checkArithmetic<Tensor1>(testSuite);
+  checkArithmetic<Tensor2>(testSuite);
+  checkArithmetic<Tensor3>(testSuite);
+
+  return testSuite.exit();
+}