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