Reference documentation for deal.II version GIT 7b2de2f2f9 2023-09-24 11:00:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tensor_accessors.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_tensor_accessors_h
17 #define dealii_tensor_accessors_h
18 
19 #include <deal.II/base/config.h>
20 
23 
24 
26 
69 namespace TensorAccessors
70 {
71  // forward declarations
72  namespace internal
73  {
74  template <int index, int rank, typename T>
75  class ReorderedIndexView;
76  template <int position, int rank>
77  struct ExtractHelper;
78  template <int no_contr, int rank_1, int rank_2, int dim>
79  class Contract;
80  template <int rank_1, int rank_2, int dim>
81  class Contract3;
82  } // namespace internal
83 
84 
101  template <typename T>
102  struct ValueType
103  {
104  using value_type = typename T::value_type;
105  };
106 
107  template <typename T>
108  struct ValueType<const T>
109  {
110  using value_type = const typename T::value_type;
111  };
112 
113  template <typename T, std::size_t N>
114  struct ValueType<T[N]>
115  {
116  using value_type = T;
117  };
118 
119  template <typename T, std::size_t N>
120  struct ValueType<const T[N]>
121  {
122  using value_type = const T;
123  };
124 
125 
133  template <int deref_steps, typename T>
134  struct ReturnType
135  {
136  using value_type =
137  typename ReturnType<deref_steps - 1,
139  };
140 
141  template <typename T>
142  struct ReturnType<0, T>
143  {
144  using value_type = T;
145  };
146 
147 
185  template <int index, int rank, typename T>
188  {
189  static_assert(0 <= index && index < rank,
190  "The specified index must lie within the range [0,rank)");
191 
193  }
194 
195 
217  template <int rank, typename T, typename ArrayType>
219  extract(T &t, const ArrayType &indices)
220  {
221  return internal::ExtractHelper<0, rank>::template extract<T, ArrayType>(
222  t, indices);
223  }
224 
225 
264  template <int no_contr,
265  int rank_1,
266  int rank_2,
267  int dim,
268  typename T1,
269  typename T2,
270  typename T3>
271  constexpr inline DEAL_II_ALWAYS_INLINE void
272  contract(T1 &result, const T2 &left, const T3 &right)
273  {
274  static_assert(rank_1 >= no_contr,
275  "The rank of the left tensor must be "
276  "equal or greater than the number of "
277  "contractions");
278  static_assert(rank_2 >= no_contr,
279  "The rank of the right tensor must be "
280  "equal or greater than the number of "
281  "contractions");
282 
284  template contract<T1, T2, T3>(result, left, right);
285  }
286 
287 
316  template <int rank_1,
317  int rank_2,
318  int dim,
319  typename T1,
320  typename T2,
321  typename T3,
322  typename T4>
323  constexpr T1
324  contract3(const T2 &left, const T3 &middle, const T4 &right)
325  {
327  template contract3<T1, T2, T3, T4>(left, middle, right);
328  }
329 
330 
331  namespace internal
332  {
333  // -------------------------------------------------------------------------
334  // Forward declarations and type traits
335  // -------------------------------------------------------------------------
336 
337  template <int rank, typename S>
338  class StoreIndex;
339  template <typename T>
340  class Identity;
341  template <int no_contr, int dim>
342  class Contract2;
343 
353  template <typename T>
355  {
356  using type = T &;
357  };
358 
359  template <int rank, typename S>
360  struct ReferenceType<StoreIndex<rank, S>>
361  {
363  };
364 
365  template <int index, int rank, typename T>
366  struct ReferenceType<ReorderedIndexView<index, rank, T>>
367  {
369  };
370 
371 
372  // TODO: Is there a possibility to just have the following block of
373  // explanation on an internal page in doxygen? If, yes. Doxygen
374  // wizards, your call!
375 
376  // -------------------------------------------------------------------------
377  // Implementation of helper classes for reordered_index_view
378  // -------------------------------------------------------------------------
379 
380  // OK. This is utterly brutal template magic. Therefore, we will not
381  // comment on the individual internal helper classes, because this is
382  // of not much value, but explain the general recursion procedure.
383  //
384  // (In order of appearance)
385  //
386  // Our task is to reorder access to a tensor object where a specified
387  // index is moved to the end. Thus we want to construct an object
388  // <code>reordered</code> out of a <code>tensor</code> where the
389  // following access patterns are equivalent:
390  // @code
391  // tensor [i_0]...[i_index-1][i_index][i_index+1]...[i_n]
392  // reordered [i_0]...[i_index_1][i_index+1]...[i_n][i_index]
393  // @endcode
394  //
395  // The first task is to get rid of the application of
396  // [i_0]...[i_index-1]. This is a classical recursion pattern - relay
397  // the task from <index, rank> to <index-1, rank-1> by accessing the
398  // subtensor object:
399 
400  template <int index, int rank, typename T>
402  {
403  public:
405  : t_(t)
406  {}
407 
408  using value_type = ReorderedIndexView<index - 1,
409  rank - 1,
410  typename ValueType<T>::value_type>;
411 
412  // Recurse by applying index j directly:
414  operator[](unsigned int j) const
415  {
416  return value_type(t_[j]);
417  }
418 
419  private:
421  };
422 
423  // At some point we hit the condition index == 0 and rank > 1, i.e.,
424  // the first index should be reordered to the end.
425  //
426  // At this point we cannot be lazy any more and have to start storing
427  // indices because we get them in the wrong order. The user supplies
428  // [i_0][i_1]...[i_{rank - 1}]
429  // but we have to call the subtensor object with
430  // [i_{rank - 1}[i_0][i_1]...[i_{rank-2}]
431  //
432  // So give up and relay the task to the StoreIndex class:
433 
434  template <int rank, typename T>
435  class ReorderedIndexView<0, rank, T>
436  {
437  public:
439  : t_(t)
440  {}
441 
443 
445  operator[](unsigned int j) const
446  {
447  return value_type(Identity<T>(t_), j);
448  }
449 
450  private:
452  };
453 
454  // Sometimes, we're lucky and don't have to do anything. In this case
455  // just return the original tensor.
456 
457  template <typename T>
458  class ReorderedIndexView<0, 1, T>
459  {
460  public:
462  : t_(t)
463  {}
464 
465  using value_type =
467 
469  operator[](unsigned int j) const
470  {
471  return t_[j];
472  }
473 
474  private:
476  };
477 
478  // Here, Identity is a helper class to ground the recursion in
479  // StoreIndex. Its implementation is easy - we haven't stored any
480  // indices yet. So, we just provide a function apply that returns the
481  // application of an index j to the stored tensor t_:
482 
483  template <typename T>
484  class Identity
485  {
486  public:
487  constexpr Identity(typename ReferenceType<T>::type t)
488  : t_(t)
489  {}
490 
492 
494  apply(unsigned int j) const
495  {
496  return t_[j];
497  }
498 
499  private:
501  };
502 
503  // StoreIndex is a class that stores an index recursively with every
504  // invocation of operator[](unsigned int j): We do this by recursively
505  // creating a new StoreIndex class of lower rank that stores the
506  // supplied index j and holds a copy of the current class (with all
507  // other stored indices). Again, we provide an apply member function
508  // that knows how to apply an index on the highest rank and all
509  // subsequently stored indices:
510 
511  template <int rank, typename S>
513  {
514  public:
515  constexpr StoreIndex(S s, int i)
516  : s_(s)
517  , i_(i)
518  {}
519 
521 
523  operator[](unsigned int j) const
524  {
525  return value_type(*this, j);
526  }
527 
528  using return_type =
530 
531  constexpr typename ReferenceType<return_type>::type
532  apply(unsigned int j) const
533  {
534  return s_.apply(j)[i_];
535  }
536 
537  private:
538  const S s_;
539  const int i_;
540  };
541 
542  // We have to store indices until we hit rank == 1. Then, upon the next
543  // invocation of operator[](unsigned int j) we have all necessary
544  // information available to return the actual object.
545 
546  template <typename S>
547  class StoreIndex<1, S>
548  {
549  public:
550  constexpr StoreIndex(S s, int i)
551  : s_(s)
552  , i_(i)
553  {}
554 
555  using return_type =
558 
560  operator[](unsigned int j) const
561  {
562  return s_.apply(j)[i_];
563  }
564 
565  private:
566  const S s_;
567  const int i_;
568  };
569 
570 
571  // -------------------------------------------------------------------------
572  // Implementation of helper classes for extract
573  // -------------------------------------------------------------------------
574 
575  // Straightforward recursion implemented by specializing ExtractHelper
576  // for position == rank. We use the type trait ReturnType<rank, T> to
577  // have an idea what the final type will be.
578  template <int position, int rank>
580  {
581  template <typename T, typename ArrayType>
582  constexpr static typename ReturnType<rank - position, T>::value_type &
583  extract(T &t, const ArrayType &indices)
584  {
586  typename ValueType<T>::value_type,
587  ArrayType>(t[indices[position]], indices);
588  }
589  };
590 
591  // For position == rank there is nothing to extract, just return the
592  // object.
593  template <int rank>
594  struct ExtractHelper<rank, rank>
595  {
596  template <typename T, typename ArrayType>
597  constexpr static T &
598  extract(T &t, const ArrayType &)
599  {
600  return t;
601  }
602  };
603 
604 
605  // -------------------------------------------------------------------------
606  // Implementation of helper classes for contract
607  // -------------------------------------------------------------------------
608 
609  // Straightforward recursive pattern:
610  //
611  // As long as rank_1 > no_contr, assign indices from the left tensor to
612  // result. This builds up the first part of the nested outer loops:
613  //
614  // for(unsigned int i_0; i_0 < dim; ++i_0)
615  // ...
616  // for(i_; i_ < dim; ++i_)
617  // [...]
618  // result[i_0]..[i_] ... left[i_0]..[i_] ...
619 
620  template <int no_contr, int rank_1, int rank_2, int dim>
621  class Contract
622  {
623  public:
624  template <typename T1, typename T2, typename T3>
625  constexpr inline DEAL_II_ALWAYS_INLINE static void
626  contract(T1 &result, const T2 &left, const T3 &right)
627  {
628  for (unsigned int i = 0; i < dim; ++i)
630  left[i],
631  right);
632  }
633  };
634 
635  // If rank_1 == no_contr leave out the remaining no_contr indices for
636  // the contraction and assign indices from the right tensor to the
637  // result. This builds up the second part of the nested loops:
638  //
639  // for(unsigned int i_0 = 0; i_0 < dim; ++i_0)
640  // ...
641  // for(unsigned int i_ = 0; i_ < dim; ++i_)
642  // for(unsigned int j_0 = 0; j_0 < dim; ++j_0)
643  // ...
644  // for(unsigned int j_ = 0; j_ < dim; ++j_)
645  // [...]
646  // result[i_0]..[i_][j_0]..[j_] ... left[i_0]..[i_] ...
647  // right[j_0]..[j_]
648  //
649 
650  template <int no_contr, int rank_2, int dim>
651  class Contract<no_contr, no_contr, rank_2, dim>
652  {
653  public:
654  template <typename T1, typename T2, typename T3>
655  constexpr inline DEAL_II_ALWAYS_INLINE static void
656  contract(T1 &result, const T2 &left, const T3 &right)
657  {
658  for (unsigned int i = 0; i < dim; ++i)
660  left,
661  right[i]);
662  }
663  };
664 
665  // If rank_1 == rank_2 == no_contr we have built up all of the outer
666  // loop. Now, it is time to do the actual contraction:
667  //
668  // [...]
669  // {
670  // result[i_0]..[i_][j_0]..[j_] = 0.;
671  // for(unsigned int k_0 = 0; k_0 < dim; ++k_0)
672  // ...
673  // for(unsigned int k_ = 0; k_ < dim; ++k_)
674  // result[i_0]..[i_][j_0]..[j_] += left[i_0]..[i_][k_0]..[k_] *
675  // right[j_0]..[j_][k_0]..[k_];
676  // }
677  //
678  // Relay this summation to another helper class.
679 
680  template <int no_contr, int dim>
681  class Contract<no_contr, no_contr, no_contr, dim>
682  {
683  public:
684  template <typename T1, typename T2, typename T3>
685  constexpr inline DEAL_II_ALWAYS_INLINE static void
686  contract(T1 &result, const T2 &left, const T3 &right)
687  {
688  result = Contract2<no_contr, dim>::template contract2<T1>(left, right);
689  }
690  };
691 
692  // Straightforward recursion:
693  //
694  // Contract leftmost index and recurse one down.
695 
696  template <int no_contr, int dim>
697  class Contract2
698  {
699  public:
700  template <typename T1, typename T2, typename T3>
701  constexpr inline DEAL_II_ALWAYS_INLINE static T1
702  contract2(const T2 &left, const T3 &right)
703  {
704  // Some auto-differentiable numbers need explicit
705  // zero initialization.
706  if (dim == 0)
707  {
708  T1 result = ::internal::NumberType<T1>::value(0.0);
709  return result;
710  }
711  else
712  {
713  T1 result =
714  Contract2<no_contr - 1, dim>::template contract2<T1>(left[0],
715  right[0]);
716  for (unsigned int i = 1; i < dim; ++i)
717  result +=
718  Contract2<no_contr - 1, dim>::template contract2<T1>(left[i],
719  right[i]);
720  return result;
721  }
722  }
723  };
724 
725  // A contraction of two objects of order 0 is just a scalar
726  // multiplication:
727 
728  template <int dim>
729  class Contract2<0, dim>
730  {
731  public:
732  template <typename T1, typename T2, typename T3>
733  constexpr DEAL_II_ALWAYS_INLINE static T1
734  contract2(const T2 &left, const T3 &right)
735  {
736  return left * right;
737  }
738  };
739 
740 
741  // -------------------------------------------------------------------------
742  // Implementation of helper classes for contract3
743  // -------------------------------------------------------------------------
744 
745  // Fully contract three tensorial objects
746  //
747  // As long as rank_1 > 0, recurse over left and middle:
748  //
749  // for(unsigned int i_0; i_0 < dim; ++i_0)
750  // ...
751  // for(i_; i_ < dim; ++i_)
752  // [...]
753  // left[i_0]..[i_] ... middle[i_0]..[i_] ... right
754 
755  template <int rank_1, int rank_2, int dim>
756  class Contract3
757  {
758  public:
759  template <typename T1, typename T2, typename T3, typename T4>
760  constexpr static inline T1
761  contract3(const T2 &left, const T3 &middle, const T4 &right)
762  {
763  // Some auto-differentiable numbers need explicit
764  // zero initialization.
765  T1 result = ::internal::NumberType<T1>::value(0.0);
766  for (unsigned int i = 0; i < dim; ++i)
767  result += Contract3<rank_1 - 1, rank_2, dim>::template contract3<T1>(
768  left[i], middle[i], right);
769  return result;
770  }
771  };
772 
773  // If rank_1 ==0, continue to recurse over middle and right:
774  //
775  // for(unsigned int i_0; i_0 < dim; ++i_0)
776  // ...
777  // for(i_; i_ < dim; ++i_)
778  // for(unsigned int j_0; j_0 < dim; ++j_0)
779  // ...
780  // for(j_; j_ < dim; ++j_)
781  // [...]
782  // left[i_0]..[i_] ... middle[i_0]..[i_][j_0]..[j_] ...
783  // right[j_0]..[j_]
784 
785  template <int rank_2, int dim>
786  class Contract3<0, rank_2, dim>
787  {
788  public:
789  template <typename T1, typename T2, typename T3, typename T4>
790  constexpr static inline T1
791  contract3(const T2 &left, const T3 &middle, const T4 &right)
792  {
793  // Some auto-differentiable numbers need explicit
794  // zero initialization.
795  T1 result = ::internal::NumberType<T1>::value(0.0);
796  for (unsigned int i = 0; i < dim; ++i)
797  result +=
798  Contract3<0, rank_2 - 1, dim>::template contract3<T1>(left,
799  middle[i],
800  right[i]);
801  return result;
802  }
803  };
804 
805  // Contraction of three tensorial objects of rank 0 is just a scalar
806  // multiplication.
807 
808  template <int dim>
809  class Contract3<0, 0, dim>
810  {
811  public:
812  template <typename T1, typename T2, typename T3, typename T4>
813  constexpr static T1
814  contract3(const T2 &left, const T3 &middle, const T4 &right)
815  {
816  return left * middle * right;
817  }
818  };
819 
820  // -------------------------------------------------------------------------
821 
822  } /* namespace internal */
823 } /* namespace TensorAccessors */
824 
826 
827 #endif /* dealii_tensor_accessors_h */
constexpr static T1 contract2(const T2 &left, const T3 &right)
constexpr static T1 contract2(const T2 &left, const T3 &right)
constexpr static T1 contract3(const T2 &left, const T3 &middle, const T4 &right)
constexpr static T1 contract3(const T2 &left, const T3 &middle, const T4 &right)
constexpr static T1 contract3(const T2 &left, const T3 &middle, const T4 &right)
constexpr static void contract(T1 &result, const T2 &left, const T3 &right)
constexpr static void contract(T1 &result, const T2 &left, const T3 &right)
constexpr static void contract(T1 &result, const T2 &left, const T3 &right)
constexpr Identity(typename ReferenceType< T >::type t)
constexpr ReferenceType< return_type >::type apply(unsigned int j) const
typename ValueType< T >::value_type return_type
typename ReferenceType< typename ValueType< T >::value_type >::type value_type
constexpr value_type operator[](unsigned int j) const
constexpr ReorderedIndexView(typename ReferenceType< T >::type t)
constexpr value_type operator[](unsigned int j) const
constexpr ReorderedIndexView(typename ReferenceType< T >::type t)
ReorderedIndexView< index - 1, rank - 1, typename ValueType< T >::value_type > value_type
constexpr value_type operator[](unsigned int j) const
constexpr ReorderedIndexView(typename ReferenceType< T >::type t)
constexpr return_type & operator[](unsigned int j) const
typename ValueType< typename S::return_type >::value_type return_type
constexpr ReferenceType< return_type >::type apply(unsigned int j) const
typename ValueType< typename S::return_type >::value_type return_type
constexpr value_type operator[](unsigned int j) const
StoreIndex< rank - 1, StoreIndex< rank, S > > value_type
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static const char T
static const char N
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
constexpr void contract(T1 &result, const T2 &left, const T3 &right)
constexpr T1 contract3(const T2 &left, const T3 &middle, const T4 &right)
constexpr internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
typename ReturnType< deref_steps - 1, typename ValueType< T >::value_type >::value_type value_type
const typename T::value_type value_type
typename T::value_type value_type
constexpr static T & extract(T &t, const ArrayType &)
constexpr static ReturnType< rank - position, T >::value_type & extract(T &t, const ArrayType &indices)
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition: numbers.h:703