Reference documentation for deal.II version 9.1.1
\(\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  inline 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  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  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  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  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  inline DEAL_II_ALWAYS_INLINE value_type operator[](unsigned int j) const
425  {
426  return value_type(t_[j]);
427  }
428 
429  private:
430  typename ReferenceType<T>::type t_;
431  };
432 
433  // At some point we hit the condition index == 0 and rank > 1, i.e.,
434  // the first index should be reordered to the end.
435  //
436  // At this point we cannot be lazy any more and have to start storing
437  // indices because we get them in the wrong order. The user supplies
438  // [i_0][i_1]...[i_{rank - 1}]
439  // but we have to call the subtensor object with
440  // [i_{rank - 1}[i_0][i_1]...[i_{rank-2}]
441  //
442  // So give up and relay the task to the StoreIndex class:
443 
444  template <int rank, typename T>
445  class ReorderedIndexView<0, rank, T>
446  {
447  public:
448  ReorderedIndexView(typename ReferenceType<T>::type t)
449  : t_(t)
450  {}
451 
452  using value_type = StoreIndex<rank - 1, internal::Identity<T>>;
453 
454  inline DEAL_II_ALWAYS_INLINE value_type operator[](unsigned int j) const
455  {
456  return value_type(Identity<T>(t_), j);
457  }
458 
459  private:
460  typename ReferenceType<T>::type t_;
461  };
462 
463  // Sometimes, we're lucky and don't have to do anything. In this case
464  // just return the original tensor.
465 
466  template <typename T>
467  class ReorderedIndexView<0, 1, T>
468  {
469  public:
470  ReorderedIndexView(typename ReferenceType<T>::type t)
471  : t_(t)
472  {}
473 
474  using value_type =
476 
477  inline DEAL_II_ALWAYS_INLINE value_type operator[](unsigned int j) const
478  {
479  return t_[j];
480  }
481 
482  private:
483  typename ReferenceType<T>::type t_;
484  };
485 
486  // Here, Identity is a helper class to ground the recursion in
487  // StoreIndex. Its implementation is easy - we haven't stored any
488  // indices yet. So, we just provide a function apply that returns the
489  // application of an index j to the stored tensor t_:
490 
491  template <typename T>
492  class Identity
493  {
494  public:
495  Identity(typename ReferenceType<T>::type t)
496  : t_(t)
497  {}
498 
499  using return_type = typename ValueType<T>::value_type;
500 
501  inline DEAL_II_ALWAYS_INLINE typename ReferenceType<return_type>::type
502  apply(unsigned int j) const
503  {
504  return t_[j];
505  }
506 
507  private:
508  typename ReferenceType<T>::type t_;
509  };
510 
511  // StoreIndex is a class that stores an index recursively with every
512  // invocation of operator[](unsigned int j): We do this by recursively
513  // creating a new StoreIndex class of lower rank that stores the
514  // supplied index j and holds a copy of the current class (with all
515  // other stored indices). Again, we provide an apply member function
516  // that knows how to apply an index on the highest rank and all
517  // subsequently stored indices:
518 
519  template <int rank, typename S>
520  class StoreIndex
521  {
522  public:
523  StoreIndex(S s, int i)
524  : s_(s)
525  , i_(i)
526  {}
527 
528  using value_type = StoreIndex<rank - 1, StoreIndex<rank, S>>;
529 
530  inline DEAL_II_ALWAYS_INLINE value_type operator[](unsigned int j) const
531  {
532  return value_type(*this, j);
533  }
534 
535  using return_type =
536  typename ValueType<typename S::return_type>::value_type;
537 
538  inline typename ReferenceType<return_type>::type
539  apply(unsigned int j) const
540  {
541  return s_.apply(j)[i_];
542  }
543 
544  private:
545  const S s_;
546  const int i_;
547  };
548 
549  // We have to store indices until we hit rank == 1. Then, upon the next
550  // invocation of operator[](unsigned int j) we have all necessary
551  // information available to return the actual object.
552 
553  template <typename S>
554  class StoreIndex<1, S>
555  {
556  public:
557  StoreIndex(S s, int i)
558  : s_(s)
559  , i_(i)
560  {}
561 
562  using return_type =
563  typename ValueType<typename S::return_type>::value_type;
564  using value_type = return_type;
565 
566  inline DEAL_II_ALWAYS_INLINE return_type &operator[](unsigned int j) const
567  {
568  return s_.apply(j)[i_];
569  }
570 
571  private:
572  const S s_;
573  const int i_;
574  };
575 
576 
577  // -------------------------------------------------------------------------
578  // Implementation of helper classes for extract
579  // -------------------------------------------------------------------------
580 
581  // Straightforward recursion implemented by specializing ExtractHelper
582  // for position == rank. We use the type trait ReturnType<rank, T> to
583  // have an idea what the final type will be.
584  template <int position, int rank>
585  struct ExtractHelper
586  {
587  template <typename T, typename ArrayType>
588  inline static typename ReturnType<rank - position, T>::value_type &
589  extract(T &t, const ArrayType &indices)
590  {
591  return ExtractHelper<position + 1, rank>::template extract<
592  typename ValueType<T>::value_type,
593  ArrayType>(t[indices[position]], indices);
594  }
595  };
596 
597  // For position == rank there is nothing to extract, just return the
598  // object.
599  template <int rank>
600  struct ExtractHelper<rank, rank>
601  {
602  template <typename T, typename ArrayType>
603  inline static T &
604  extract(T &t, const ArrayType &)
605  {
606  return t;
607  }
608  };
609 
610 
611  // -------------------------------------------------------------------------
612  // Implementation of helper classes for contract
613  // -------------------------------------------------------------------------
614 
615  // Straightforward recursive pattern:
616  //
617  // As long as rank_1 > no_contr, assign indices from the left tensor to
618  // result. This builds up the first part of the nested outer loops:
619  //
620  // for(unsigned int i_0; i_0 < dim; ++i_0)
621  // ...
622  // for(i_; i_ < dim; ++i_)
623  // [...]
624  // result[i_0]..[i_] ... left[i_0]..[i_] ...
625 
626  template <int no_contr, int rank_1, int rank_2, int dim>
627  class Contract
628  {
629  public:
630  template <typename T1, typename T2, typename T3>
631  inline DEAL_II_ALWAYS_INLINE static void
632  contract(T1 &result, const T2 &left, const T3 &right)
633  {
634  for (unsigned int i = 0; i < dim; ++i)
636  left[i],
637  right);
638  }
639  };
640 
641  // If rank_1 == no_contr leave out the remaining no_contr indices for
642  // the contraction and assign indices from the right tensor to the
643  // result. This builds up the second part of the nested loops:
644  //
645  // for(unsigned int i_0 = 0; i_0 < dim; ++i_0)
646  // ...
647  // for(unsigned int i_ = 0; i_ < dim; ++i_)
648  // for(unsigned int j_0 = 0; j_0 < dim; ++j_0)
649  // ...
650  // for(unsigned int j_ = 0; j_ < dim; ++j_)
651  // [...]
652  // result[i_0]..[i_][j_0]..[j_] ... left[i_0]..[i_] ...
653  // right[j_0]..[j_]
654  //
655 
656  template <int no_contr, int rank_2, int dim>
657  class Contract<no_contr, no_contr, rank_2, dim>
658  {
659  public:
660  template <typename T1, typename T2, typename T3>
661  inline DEAL_II_ALWAYS_INLINE static void
662  contract(T1 &result, const T2 &left, const T3 &right)
663  {
664  for (unsigned int i = 0; i < dim; ++i)
666  left,
667  right[i]);
668  }
669  };
670 
671  // If rank_1 == rank_2 == no_contr we have built up all of the outer
672  // loop. Now, it is time to do the actual contraction:
673  //
674  // [...]
675  // {
676  // result[i_0]..[i_][j_0]..[j_] = 0.;
677  // for(unsigned int k_0 = 0; k_0 < dim; ++k_0)
678  // ...
679  // for(unsigned int k_ = 0; k_ < dim; ++k_)
680  // result[i_0]..[i_][j_0]..[j_] += left[i_0]..[i_][k_0]..[k_] *
681  // right[j_0]..[j_][k_0]..[k_];
682  // }
683  //
684  // Relay this summation to another helper class.
685 
686  template <int no_contr, int dim>
687  class Contract<no_contr, no_contr, no_contr, dim>
688  {
689  public:
690  template <typename T1, typename T2, typename T3>
691  inline DEAL_II_ALWAYS_INLINE static void
692  contract(T1 &result, const T2 &left, const T3 &right)
693  {
694  result = Contract2<no_contr, dim>::template contract2<T1>(left, right);
695  }
696  };
697 
698  // Straightforward recursion:
699  //
700  // Contract leftmost index and recurse one down.
701 
702  template <int no_contr, int dim>
703  class Contract2
704  {
705  public:
706  template <typename T1, typename T2, typename T3>
707  inline DEAL_II_ALWAYS_INLINE static T1
708  contract2(const T2 &left, const T3 &right)
709  {
710  // Some auto-differentiable numbers need explicit
711  // zero initialization.
712  T1 result = ::internal::NumberType<T1>::value(0.0);
713  for (unsigned int i = 0; i < dim; ++i)
714  result +=
715  Contract2<no_contr - 1, dim>::template contract2<T1>(left[i],
716  right[i]);
717  return result;
718  }
719  };
720 
721  // A contraction of two objects of order 0 is just a scalar
722  // multiplication:
723 
724  template <int dim>
725  class Contract2<0, dim>
726  {
727  public:
728  template <typename T1, typename T2, typename T3>
729  inline DEAL_II_ALWAYS_INLINE static T1
730  contract2(const T2 &left, const T3 &right)
731  {
732  return left * right;
733  }
734  };
735 
736 
737  // -------------------------------------------------------------------------
738  // Implementation of helper classes for contract3
739  // -------------------------------------------------------------------------
740 
741  // Fully contract three tensorial objects
742  //
743  // As long as rank_1 > 0, recurse over left and middle:
744  //
745  // for(unsigned int i_0; i_0 < dim; ++i_0)
746  // ...
747  // for(i_; i_ < dim; ++i_)
748  // [...]
749  // left[i_0]..[i_] ... middle[i_0]..[i_] ... right
750 
751  template <int rank_1, int rank_2, int dim>
752  class Contract3
753  {
754  public:
755  template <typename T1, typename T2, typename T3, typename T4>
756  static inline T1
757  contract3(const T2 &left, const T3 &middle, const T4 &right)
758  {
759  // Some auto-differentiable numbers need explicit
760  // zero initialization.
761  T1 result = ::internal::NumberType<T1>::value(0.0);
762  for (unsigned int i = 0; i < dim; ++i)
763  result += Contract3<rank_1 - 1, rank_2, dim>::template contract3<T1>(
764  left[i], middle[i], right);
765  return result;
766  }
767  };
768 
769  // If rank_1 ==0, continue to recurse over middle and right:
770  //
771  // for(unsigned int i_0; i_0 < dim; ++i_0)
772  // ...
773  // for(i_; i_ < dim; ++i_)
774  // for(unsigned int j_0; j_0 < dim; ++j_0)
775  // ...
776  // for(j_; j_ < dim; ++j_)
777  // [...]
778  // left[i_0]..[i_] ... middle[i_0]..[i_][j_0]..[j_] ...
779  // right[j_0]..[j_]
780 
781  template <int rank_2, int dim>
782  class Contract3<0, rank_2, dim>
783  {
784  public:
785  template <typename T1, typename T2, typename T3, typename T4>
786  static inline T1
787  contract3(const T2 &left, const T3 &middle, const T4 &right)
788  {
789  // Some auto-differentiable numbers need explicit
790  // zero initialization.
791  T1 result = ::internal::NumberType<T1>::value(0.0);
792  for (unsigned int i = 0; i < dim; ++i)
793  result +=
794  Contract3<0, rank_2 - 1, dim>::template contract3<T1>(left,
795  middle[i],
796  right[i]);
797  return result;
798  }
799  };
800 
801  // Contraction of three tensorial objects of rank 0 is just a scalar
802  // multiplication.
803 
804  template <int dim>
805  class Contract3<0, 0, dim>
806  {
807  public:
808  template <typename T1, typename T2, typename T3, typename T4>
809  static inline T1
810  contract3(const T2 &left, const T3 &middle, const T4 &right)
811  {
812  return left * middle * right;
813  }
814  };
815 
816  // -------------------------------------------------------------------------
817 
818  } /* namespace internal */
819 } /* namespace TensorAccessors */
820 
821 DEAL_II_NAMESPACE_CLOSE
822 
823 #endif /* dealii_tensor_accessors_h */
internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
T1 contract3(const T2 &left, const T3 &middle, const T4 &right)
void contract(T1 &result, const T2 &left, const T3 &right)
void contract(Tensor< 2, dim, Number > &dest, const Tensor< 2, dim, Number > &src1, const unsigned int index1, const Tensor< 2, dim, Number > &src2, const unsigned int index3)