tensor.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2019-2023 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
#ifndef TENSOR_HH
#define TENSOR_HH
 
 
#include "dune/common/fvector.hh"
#include "dune/common/fmatrix.hh"
 
#include "utilities/scalar.hh"
 
namespace Kaskade
{
 
  template <int first, int last>
  struct StaticIndexRange {};
 
  // forward declaration
  template <class Entry, int Size0, int ... Sizes>
  class Tensor;
 
  namespace TensorDetail
  {
    template <class Entry, int Size>
    struct Prepend { using type = Tensor<Entry,Size>; };
 
    template <class Entry, int SizePre, int Size0, int ... Sizes>
    struct Prepend<Tensor<Entry,Size0,Sizes...>,SizePre> { using type = Tensor<Entry,SizePre,Size0,Sizes...>; };
  }
 
  // Provide the rank of tensors with the same interface as for vectors and matrices.
  namespace ScalarDetail
  {
      template <class K, int Size0, int ... Sizes>
      struct Rank<Tensor<K,Size0,Sizes...>>
      {
        static constexpr int value = Tensor<K,Size0,Sizes...>::rank;
      };
  }
  namespace TensorAccess
  {
    constexpr StaticIndexRange<0,std::numeric_limits<int>::max()> _;
  }
 
  template <class Entry, int Size0, int ... Sizes>
  class Tensor: public Dune::FieldVector<Tensor<Entry,Sizes...>,Size0>
  {
    using Base = Dune::FieldVector<Tensor<Entry,Sizes...>,Size0>;
    static_assert(Size0>0);
    static constexpr int Size1 = Tensor<Entry,Sizes...>::dimension;
  public:
 
    static constexpr int rank = Tensor<Entry,Sizes...>::rank + 1;
 
    Tensor() = default;
 
    Tensor(Entry const& e)
    {
      (*this) = e;
    }
 
    Tensor& operator=(Entry const& e)
    {
      for (auto& t: *this)
        t = e;
      return *this;
    }
 
    Tensor& operator=(Dune::FieldMatrix<Entry,Size0,Size1> const& a)
    {
      static_assert(rank==2);
      if constexpr (rank==2)
        for (int i=0; i<Size0; ++i)
            (*this)[i] = a[i];
      return *this;
    }
    template <class ... Access>
    auto operator()(int i0, Access... is) const
    {
      return (*this)[i0](is...);
    }
 
    template <int first, int last, class ... Access>
    auto operator()(StaticIndexRange<first,last> , Access... is) const
    {
      constexpr int s = std::min(Size0,last-first);
      using SubT = decltype( (*this)[0](is...) );
      using R = typename TensorDetail::Prepend<SubT,s>::type;
      R r;
      for (int i=0; i<s; ++i)
        r[i] = (*this)[i+first](is...);
      return r;
    }
 
    operator std::conditional_t<rank==2,
                                Dune::FieldMatrix<Entry,Size0,Base::block_type::dimension>,
                                Base&> () const
    {
      if constexpr(rank==2)
      {
        Dune::FieldMatrix<Entry,Size0,Base::value_type::dimension> A;
        // We want that A[i][j] = (*this)[i][j]. A[i] are the rows of A (of type FieldVector),
        // and (*this)[i] are Tensors of rank 1, and thus derived from FieldVector. So we can
        // assign the matrix row by row.
        for (int i=0; i<Size0; ++i)
          A[i] = (*this)[i];
        return A;
      }
      else
        return *this;
    }
 
  };
 
  // ----------------------------------------------------------------------------------------------
 
  // Terminal of parameter pack recursion
  template <class Entry, int Size0>
  class Tensor<Entry,Size0>: public Dune::FieldVector<Entry,Size0>
  {
    using Base = Dune::FieldVector<Entry,Size0>;
  public:
 
    using field_type = typename EntryTraits<Entry>::field_type;
 
    static constexpr int rank = 1;
 
    Tensor() = default;
 
    Tensor(Entry const& e)
    {
      (*this) = e;
    }
 
    Tensor& operator=(Entry const& e)
    {
      for (auto& t: *this)
        t = e;
      return *this;
    }
 
    Tensor& operator=(Base const& v)
    {
      static_cast<Base&>(*this) = v;
      return *this;
    }
 
    Entry operator()(int i) const
    {
      return (*this)[i];
    }
 
    template <int first, int last>
    auto operator()(StaticIndexRange<first,last>) const
    {
      constexpr int s = std::min(Size0,last-first);
 
      Tensor<Entry,s> r;
      for (int i=0; i<s; ++i)
        r[i] = (*this)[i+first];
      return r;
    }
 
    Tensor& operator*=(field_type const& x)
    {
      static_cast<Base&>(*this) *= x;
      return *this;
    }
  };
 
 
  // Template specialization for accessing type information.
  template <class Entry, int ... Sizes>
  struct EntryTraits<Tensor<Entry,Sizes...>>
  {
    using field_type = typename EntryTraits<Entry>::field_type;
    using real_type = typename ScalarTraits<field_type>::Real;
    static int const lapackLayout = false;
  };
 
  // ----------------------------------------------------------------------------------------------
 
  template <class Entry, int Size0, int ...Sizes>
  Tensor<Entry,Size0,Sizes...> operator * (typename EntryTraits<Entry>::field_type a, Tensor<Entry,Size0,Sizes...> A)
  {
    A *= a;
    return A;
  }
 
  // ----------------------------------------------------------------------------------------------
 
  template <class Entry, int n>
  Entry trace(Tensor<Entry,n,n> const& A)
  {
    Entry t = 0;
    for (int i=0; i<n; ++i)
      t += A[i][i];
    return t;
  }
 
  // ----------------------------------------------------------------------------------------------
 
  // Tensor output
  template <class Entry, int ... Sizes>
  std::ostream& operator<<(std::ostream& out, Tensor<Entry,Sizes...> const& t)
  {
    out << "[" << t[0];
    for (int i=1; i<std::size(t); ++i)
      out << "," << t[i];
    out << "]";
    return out;
  }
 
}
 
// ------------------------------------------------------------------------------------------------
 
// ------------------------------------------------------------------------------------------------
 
namespace Dune
{
  // For interoperability with Dune (and as we put Tensors as entries into FieldVectors,
  // we need to specify the underlying field type.
  template <class K, int Size0, int ... Sizes>
  struct FieldTraits<Kaskade::Tensor<K,Size0,Sizes...>>
  {
    using field_type = typename FieldTraits<K>::field_type;
    using real_type = typename FieldTraits<K>::real_type;
  };
}
#endif