Reference documentation for deal.II version Git a0b41b6d0f 2020-02-26 20:08:13 -0600
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tensor.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 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_h
17 #define dealii_tensor_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/numbers.h>
23 #include <deal.II/base/std_cxx14/utility.h>
24 #include <deal.II/base/table_indices.h>
25 #include <deal.II/base/template_constraints.h>
26 #include <deal.II/base/tensor_accessors.h>
27 #include <deal.II/base/utilities.h>
28 
29 #ifdef DEAL_II_WITH_ADOLC
30 # include <adolc/adouble.h> // Taped double
31 #endif
32 
33 #include <cmath>
34 #include <ostream>
35 #include <vector>
36 
37 
38 DEAL_II_NAMESPACE_OPEN
39 
40 // Forward declarations:
41 #ifndef DOXYGEN
42 template <int dim, typename Number>
43 class Point;
44 template <int rank_, int dim, typename Number = double>
45 class Tensor;
46 template <typename Number>
47 class Vector;
48 namespace Differentiation
49 {
50  namespace SD
51  {
52  class Expression;
53  }
54 } // namespace Differentiation
55 #endif
56 
57 #ifndef DOXYGEN
58 // Overload invalid tensor types of negative rank that come up during
59 // overload resolution of operator* and related contraction variants.
60 template <int dim, typename Number>
61 class Tensor<-2, dim, Number>
62 {};
63 
64 template <int dim, typename Number>
65 class Tensor<-1, dim, Number>
66 {};
67 #endif /* DOXYGEN */
68 
69 
100 template <int dim, typename Number>
101 class Tensor<0, dim, Number>
102 {
103 public:
112  static constexpr unsigned int dimension = dim;
113 
117  static constexpr unsigned int rank = 0;
118 
122  static constexpr unsigned int n_independent_components = 1;
123 
133 
138  using value_type = Number;
139 
145  using array_type = Number;
146 
152  constexpr DEAL_II_CUDA_HOST_DEV
153  Tensor();
154 
162  template <typename OtherNumber>
163  constexpr DEAL_II_CUDA_HOST_DEV
164  Tensor(const Tensor<0, dim, OtherNumber> &initializer);
165 
171  template <typename OtherNumber>
172  constexpr DEAL_II_CUDA_HOST_DEV
173  Tensor(const OtherNumber &initializer);
174 
178  Number *
179  begin_raw();
180 
184  const Number *
185  begin_raw() const;
186 
190  Number *
191  end_raw();
192 
197  const Number *
198  end_raw() const;
199 
209  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV operator Number &();
210 
219  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV operator const Number &() const;
220 
228  template <typename OtherNumber>
229  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
230  operator=(const Tensor<0, dim, OtherNumber> &rhs);
231 
232 #ifdef __INTEL_COMPILER
233 
241  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
242  operator=(const Tensor<0, dim, Number> &rhs);
243 #endif
244 
251  template <typename OtherNumber>
252  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
253  operator=(const OtherNumber &d);
254 
258  template <typename OtherNumber>
259  DEAL_II_CONSTEXPR bool
260  operator==(const Tensor<0, dim, OtherNumber> &rhs) const;
261 
265  template <typename OtherNumber>
266  constexpr bool
267  operator!=(const Tensor<0, dim, OtherNumber> &rhs) const;
268 
274  template <typename OtherNumber>
275  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
276  operator+=(const Tensor<0, dim, OtherNumber> &rhs);
277 
283  template <typename OtherNumber>
284  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
285  operator-=(const Tensor<0, dim, OtherNumber> &rhs);
286 
292  template <typename OtherNumber>
293  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
294  operator*=(const OtherNumber &factor);
295 
301  template <typename OtherNumber>
302  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
303  operator/=(const OtherNumber &factor);
304 
310  constexpr DEAL_II_CUDA_HOST_DEV Tensor
311  operator-() const;
312 
325  DEAL_II_CONSTEXPR void
326  clear();
327 
333  real_type
334  norm() const;
335 
342  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV real_type
343  norm_square() const;
344 
349  template <class Archive>
350  void
351  serialize(Archive &ar, const unsigned int version);
352 
357  using tensor_type = Number;
358 
359 private:
363  Number value;
364 
368  template <typename OtherNumber>
369  void
370  unroll_recursion(Vector<OtherNumber> &result,
371  unsigned int & start_index) const;
372 
373  // Allow an arbitrary Tensor to access the underlying values.
374  template <int, int, typename>
375  friend class Tensor;
376 };
377 
378 
379 
421 template <int rank_, int dim, typename Number>
422 class Tensor
423 {
424 public:
433  static constexpr unsigned int dimension = dim;
434 
438  static constexpr unsigned int rank = rank_;
439 
444  static constexpr unsigned int n_independent_components =
445  Tensor<rank_ - 1, dim>::n_independent_components * dim;
446 
452  using value_type = typename Tensor<rank_ - 1, dim, Number>::tensor_type;
453 
458  using array_type =
459  typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1];
460 
466  constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
467  Tensor()
468 #ifdef DEAL_II_MSVC
469  : values{}
470  {}
471 #else
472  = default;
473 #endif
474 
480  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
481  const array_type &initializer);
482 
490  template <typename OtherNumber>
491  constexpr DEAL_II_CUDA_HOST_DEV
492  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
493 
497  template <typename OtherNumber>
498  constexpr Tensor(
499  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
500 
504  template <typename OtherNumber>
505  constexpr
506  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
507 
513  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV value_type &
514  operator[](const unsigned int i);
515 
521  constexpr DEAL_II_CUDA_HOST_DEV const value_type &
522  operator[](const unsigned int i) const;
523 
527  DEAL_II_CONSTEXPR const Number &
528  operator[](const TableIndices<rank_> &indices) const;
529 
533  DEAL_II_CONSTEXPR Number &operator[](const TableIndices<rank_> &indices);
534 
538  Number *
539  begin_raw();
540 
544  const Number *
545  begin_raw() const;
546 
550  Number *
551  end_raw();
552 
556  const Number *
557  end_raw() const;
558 
566  template <typename OtherNumber>
567  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
568  operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
569 
576  DEAL_II_CONSTEXPR Tensor &
577  operator=(const Number &d);
578 
582  template <typename OtherNumber>
583  DEAL_II_CONSTEXPR bool
585 
589  template <typename OtherNumber>
590  constexpr bool
592 
598  template <typename OtherNumber>
599  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
600  operator+=(const Tensor<rank_, dim, OtherNumber> &);
601 
607  template <typename OtherNumber>
608  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
609  operator-=(const Tensor<rank_, dim, OtherNumber> &);
610 
617  template <typename OtherNumber>
618  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
619  operator*=(const OtherNumber &factor);
620 
626  template <typename OtherNumber>
627  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
628  operator/=(const OtherNumber &factor);
629 
635  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor
636  operator-() const;
637 
650  DEAL_II_CONSTEXPR void
651  clear();
652 
660  DEAL_II_CUDA_HOST_DEV
662  norm() const;
663 
670  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV
672  norm_square() const;
673 
681  template <typename OtherNumber>
682  void
683  unroll(Vector<OtherNumber> &result) const;
684 
689  static DEAL_II_CONSTEXPR unsigned int
690  component_to_unrolled_index(const TableIndices<rank_> &indices);
691 
696  static DEAL_II_CONSTEXPR TableIndices<rank_>
697  unrolled_to_component_indices(const unsigned int i);
698 
703  static constexpr std::size_t
704  memory_consumption();
705 
710  template <class Archive>
711  void
712  serialize(Archive &ar, const unsigned int version);
713 
719 
720 private:
724  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
725  // ... avoid a compiler warning in case of dim == 0 and ensure that the
726  // array always has positive size.
727 
731  template <typename OtherNumber>
732  void
733  unroll_recursion(Vector<OtherNumber> &result,
734  unsigned int & start_index) const;
735 
742  template <typename ArrayLike, std::size_t... Indices>
743  constexpr DEAL_II_CUDA_HOST_DEV
744  Tensor(const ArrayLike &initializer, std_cxx14::index_sequence<Indices...>);
745 
746  // Allow an arbitrary Tensor to access the underlying values.
747  template <int, int, typename>
748  friend class Tensor;
749 
750  // Point is allowed access to the coordinates. This is supposed to improve
751  // speed.
752  friend class Point<dim, Number>;
753 };
754 
755 
756 namespace internal
757 {
758  template <int rank, int dim, typename T, typename U>
759  struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
760  {
761  using type =
763  };
764 
765  template <typename T, int rank, int dim, typename U>
766  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
767  {
768  using type =
770  };
771 
776  template <int rank, int dim, typename T>
777  struct NumberType<Tensor<rank, dim, T>>
778  {
779  static constexpr DEAL_II_ALWAYS_INLINE const Tensor<rank, dim, T> &
780  value(const Tensor<rank, dim, T> &t)
781  {
782  return t;
783  }
784 
785  static DEAL_II_CONSTEXPR DEAL_II_ALWAYS_INLINE Tensor<rank, dim, T>
786  value(const T &t)
787  {
789  tmp = t;
790  return tmp;
791  }
792  };
793 } // namespace internal
794 
795 
796 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
797 
798 
799 template <int dim, typename Number>
800 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
802  // Some auto-differentiable numbers need explicit
803  // zero initialization such as adtl::adouble.
804  : Tensor{0.0}
805 {}
806 
807 
808 
809 template <int dim, typename Number>
810 template <typename OtherNumber>
811 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
812  Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
813  : value(internal::NumberType<Number>::value(initializer))
814 {}
815 
816 
817 
818 template <int dim, typename Number>
819 template <typename OtherNumber>
820 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
822  : Tensor{p.value}
823 {}
824 
825 
826 
827 template <int dim, typename Number>
828 inline Number *
830 {
831  return std::addressof(value);
832 }
833 
834 
835 
836 template <int dim, typename Number>
837 inline const Number *
839 {
840  return std::addressof(value);
841 }
842 
843 
844 
845 template <int dim, typename Number>
846 inline Number *
848 {
850 }
851 
852 
853 
854 template <int dim, typename Number>
855 const Number *
857 {
859 }
860 
861 
862 
863 template <int dim, typename Number>
864 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
865  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::operator Number &()
866 {
867  // We cannot use Assert inside a CUDA kernel
868 #ifndef __CUDA_ARCH__
869  Assert(dim != 0,
870  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
871 #endif
872  return value;
873 }
874 
875 
876 template <int dim, typename Number>
877 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
878  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::operator const Number &() const
879 {
880  // We cannot use Assert inside a CUDA kernel
881 #ifndef __CUDA_ARCH__
882  Assert(dim != 0,
883  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
884 #endif
885  return value;
886 }
887 
888 
889 template <int dim, typename Number>
890 template <typename OtherNumber>
891 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
892  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
894 {
895  value = internal::NumberType<Number>::value(p);
896  return *this;
897 }
898 
899 
900 #ifdef __INTEL_COMPILER
901 template <int dim, typename Number>
902 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
903  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
905 {
906  value = p.value;
907  return *this;
908 }
909 #endif
910 
911 
912 template <int dim, typename Number>
913 template <typename OtherNumber>
914 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
915  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
916  Tensor<0, dim, Number>::operator=(const OtherNumber &d)
917 {
918  value = internal::NumberType<Number>::value(d);
919  return *this;
920 }
921 
922 
923 template <int dim, typename Number>
924 template <typename OtherNumber>
925 DEAL_II_CONSTEXPR inline bool
927 {
928 #if defined(DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING)
929  Assert(!(std::is_same<Number, adouble>::value ||
930  std::is_same<OtherNumber, adouble>::value),
931  ExcMessage(
932  "The Tensor equality operator for ADOL-C taped numbers has not yet "
933  "been extended to support advanced branching."));
934 #endif
935 
936  return numbers::values_are_equal(value, p.value);
937 }
938 
939 
940 template <int dim, typename Number>
941 template <typename OtherNumber>
942 constexpr bool
944 {
945  return !((*this) == p);
946 }
947 
948 
949 template <int dim, typename Number>
950 template <typename OtherNumber>
951 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
952  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
954 {
955  value += p.value;
956  return *this;
957 }
958 
959 
960 template <int dim, typename Number>
961 template <typename OtherNumber>
962 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
963  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
965 {
966  value -= p.value;
967  return *this;
968 }
969 
970 
971 
972 namespace internal
973 {
974  namespace ComplexWorkaround
975  {
976  template <typename Number, typename OtherNumber>
977  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV void
978  multiply_assign_scalar(Number &val, const OtherNumber &s)
979  {
980  val *= s;
981  }
982 
983 #ifdef __CUDA_ARCH__
984  template <typename Number, typename OtherNumber>
985  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV void
986  multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
987  {
988  printf("This function is not implemented for std::complex<Number>!\n");
989  assert(false);
990  }
991 #endif
992  } // namespace ComplexWorkaround
993 } // namespace internal
994 
995 
996 template <int dim, typename Number>
997 template <typename OtherNumber>
998 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
999  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
1000  Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1001 {
1002  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1003  return *this;
1004 }
1005 
1006 
1007 
1008 template <int dim, typename Number>
1009 template <typename OtherNumber>
1010 DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
1011 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1012 {
1013  value /= s;
1014  return *this;
1015 }
1016 
1017 
1018 template <int dim, typename Number>
1019 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>
1021 {
1022  return -value;
1023 }
1024 
1025 
1026 template <int dim, typename Number>
1027 inline DEAL_II_ALWAYS_INLINE typename Tensor<0, dim, Number>::real_type
1029 {
1030  Assert(dim != 0,
1031  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1032  return numbers::NumberTraits<Number>::abs(value);
1033 }
1034 
1035 
1036 template <int dim, typename Number>
1037 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1040 {
1041  // We cannot use Assert inside a CUDA kernel
1042 #ifndef __CUDA_ARCH__
1043  Assert(dim != 0,
1044  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1045 #endif
1047 }
1048 
1049 
1050 template <int dim, typename Number>
1051 template <typename OtherNumber>
1052 inline void
1053 Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1054  unsigned int & index) const
1055 {
1056  Assert(dim != 0,
1057  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1058  result[index] = value;
1059  ++index;
1060 }
1061 
1062 
1063 template <int dim, typename Number>
1064 DEAL_II_CONSTEXPR inline void
1066 {
1067  // Some auto-differentiable numbers need explicit
1068  // zero initialization.
1069  value = internal::NumberType<Number>::value(0.0);
1070 }
1071 
1072 
1073 template <int dim, typename Number>
1074 template <class Archive>
1075 inline void
1076 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1077 {
1078  ar &value;
1079 }
1080 
1081 
1082 template <int dim, typename Number>
1084 
1085 
1086 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1087 
1088 template <int rank_, int dim, typename Number>
1089 template <typename ArrayLike, std::size_t... indices>
1090 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1091  Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1092  std_cxx14::index_sequence<indices...>)
1093  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
1094 {
1095  static_assert(sizeof...(indices) == dim,
1096  "dim should match the number of indices");
1097 }
1098 
1099 
1100 template <int rank_, int dim, typename Number>
1101 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1103  : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
1104 {}
1105 
1106 
1107 template <int rank_, int dim, typename Number>
1108 template <typename OtherNumber>
1109 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1111  const Tensor<rank_, dim, OtherNumber> &initializer)
1112  : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
1113 {}
1114 
1115 
1116 template <int rank_, int dim, typename Number>
1117 template <typename OtherNumber>
1118 constexpr DEAL_II_ALWAYS_INLINE
1120  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1121  : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
1122 {}
1123 
1124 
1125 template <int rank_, int dim, typename Number>
1126 template <typename OtherNumber>
1127 constexpr DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number>::
1128  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1129 {
1130  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
1131 }
1132 
1133 
1134 
1135 namespace internal
1136 {
1137  namespace TensorSubscriptor
1138  {
1139  template <typename ArrayElementType, int dim>
1140  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1141  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1142  subscript(ArrayElementType * values,
1143  const unsigned int i,
1144  std::integral_constant<int, dim>)
1145  {
1146  // We cannot use Assert in a CUDA kernel
1147 #ifndef __CUDA_ARCH__
1148  AssertIndexRange(i, dim);
1149 #endif
1150  return values[i];
1151  }
1152 
1153  // The variables within this struct will be referenced in the next function.
1154  // It is a workaround that allows returning a reference to a static variable
1155  // while allowing constexpr evaluation of the function.
1156  // It has to be defined outside the function because constexpr functions
1157  // cannot define static variables
1158  template <typename ArrayElementType>
1159  struct Uninitialized
1160  {
1161  static ArrayElementType value;
1162  };
1163 
1164  template <typename Type>
1165  Type Uninitialized<Type>::value;
1166 
1167  template <typename ArrayElementType>
1168  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1169  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1170  subscript(ArrayElementType *,
1171  const unsigned int,
1172  std::integral_constant<int, 0>)
1173  {
1174  // We cannot use Assert in a CUDA kernel
1175 #ifndef __CUDA_ARCH__
1176  Assert(
1177  false,
1178  ExcMessage(
1179  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1180 #endif
1181  return Uninitialized<ArrayElementType>::value;
1182  }
1183  } // namespace TensorSubscriptor
1184 } // namespace internal
1185 
1186 
1187 template <int rank_, int dim, typename Number>
1188 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1190  operator[](const unsigned int i)
1191 {
1192  return ::internal::TensorSubscriptor::subscript(
1193  values, i, std::integral_constant<int, dim>());
1194 }
1195 
1196 
1197 template <int rank_, int dim, typename Number>
1198 constexpr DEAL_II_ALWAYS_INLINE
1199  DEAL_II_CUDA_HOST_DEV const typename Tensor<rank_, dim, Number>::value_type &
1200  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1201 {
1202  return values[i];
1203 }
1204 
1205 
1206 template <int rank_, int dim, typename Number>
1207 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1209  operator[](const TableIndices<rank_> &indices) const
1210 {
1211  Assert(dim != 0,
1212  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1213 
1214  return TensorAccessors::extract<rank_>(*this, indices);
1215 }
1216 
1217 
1218 
1219 template <int rank_, int dim, typename Number>
1220 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number &
1222 {
1223  Assert(dim != 0,
1224  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1225 
1226  return TensorAccessors::extract<rank_>(*this, indices);
1227 }
1228 
1229 
1230 
1231 template <int rank_, int dim, typename Number>
1232 inline Number *
1234 {
1235  return std::addressof(
1236  this->operator[](this->unrolled_to_component_indices(0)));
1237 }
1238 
1239 
1240 
1241 template <int rank_, int dim, typename Number>
1242 inline const Number *
1244 {
1245  return std::addressof(
1246  this->operator[](this->unrolled_to_component_indices(0)));
1247 }
1248 
1249 
1250 
1251 template <int rank_, int dim, typename Number>
1252 inline Number *
1254 {
1256 }
1257 
1258 
1259 
1260 template <int rank_, int dim, typename Number>
1261 inline const Number *
1263 {
1265 }
1266 
1267 
1268 
1269 template <int rank_, int dim, typename Number>
1270 template <typename OtherNumber>
1271 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
1273 {
1274  // The following loop could be written more concisely using std::copy, but
1275  // that function is only constexpr from C++20 on.
1276  for (unsigned int i = 0; i < dim; ++i)
1277  values[i] = t.values[i];
1278  return *this;
1279 }
1280 
1281 
1282 template <int rank_, int dim, typename Number>
1283 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
1285 {
1287  ExcMessage("Only assignment with zero is allowed"));
1288  (void)d;
1289 
1290  for (unsigned int i = 0; i < dim; ++i)
1291  values[i] = internal::NumberType<Number>::value(0.0);
1292  return *this;
1293 }
1294 
1295 
1296 template <int rank_, int dim, typename Number>
1297 template <typename OtherNumber>
1298 DEAL_II_CONSTEXPR inline bool
1301 {
1302  for (unsigned int i = 0; i < dim; ++i)
1303  if (values[i] != p.values[i])
1304  return false;
1305  return true;
1306 }
1307 
1308 
1309 // At some places in the library, we have Point<0> for formal reasons
1310 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1311 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1312 // in the above function that the loop end check always fails, we
1313 // implement this function here
1314 template <>
1315 template <>
1316 DEAL_II_CONSTEXPR inline bool
1318 {
1319  return true;
1320 }
1321 
1322 
1323 template <int rank_, int dim, typename Number>
1324 template <typename OtherNumber>
1325 constexpr bool
1328 {
1329  return !((*this) == p);
1330 }
1331 
1332 
1333 template <int rank_, int dim, typename Number>
1334 template <typename OtherNumber>
1335 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1336  DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
1339 {
1340  for (unsigned int i = 0; i < dim; ++i)
1341  values[i] += p.values[i];
1342  return *this;
1343 }
1344 
1345 
1346 template <int rank_, int dim, typename Number>
1347 template <typename OtherNumber>
1348 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1349  DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
1352 {
1353  for (unsigned int i = 0; i < dim; ++i)
1354  values[i] -= p.values[i];
1355  return *this;
1356 }
1357 
1358 
1359 template <int rank_, int dim, typename Number>
1360 template <typename OtherNumber>
1361 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1362  DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
1363  Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1364 {
1365  for (unsigned int i = 0; i < dim; ++i)
1366  values[i] *= s;
1367  return *this;
1368 }
1369 
1370 
1371 namespace internal
1372 {
1373  namespace TensorImplementation
1374  {
1375  template <int rank,
1376  int dim,
1377  typename Number,
1378  typename OtherNumber,
1379  typename std::enable_if<
1380  !std::is_integral<
1381  typename ProductType<Number, OtherNumber>::type>::value &&
1382  !std::is_same<Number, Differentiation::SD::Expression>::value,
1383  int>::type = 0>
1384  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
1385  division_operator(Tensor<rank, dim, Number> (&t)[dim],
1386  const OtherNumber &factor)
1387  {
1388  const Number inverse_factor = Number(1.) / factor;
1389  // recurse over the base objects
1390  for (unsigned int d = 0; d < dim; ++d)
1391  t[d] *= inverse_factor;
1392  }
1393 
1394 
1395  template <int rank,
1396  int dim,
1397  typename Number,
1398  typename OtherNumber,
1399  typename std::enable_if<
1400  std::is_integral<
1401  typename ProductType<Number, OtherNumber>::type>::value ||
1402  std::is_same<Number, Differentiation::SD::Expression>::value,
1403  int>::type = 0>
1404  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
1405  division_operator(Tensor<rank, dim, Number> (&t)[dim],
1406  const OtherNumber &factor)
1407  {
1408  // recurse over the base objects
1409  for (unsigned int d = 0; d < dim; ++d)
1410  t[d] /= factor;
1411  }
1412  } // namespace TensorImplementation
1413 } // namespace internal
1414 
1415 
1416 template <int rank_, int dim, typename Number>
1417 template <typename OtherNumber>
1418 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1419  DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
1420  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1421 {
1422  internal::TensorImplementation::division_operator(values, s);
1423  return *this;
1424 }
1425 
1426 
1427 template <int rank_, int dim, typename Number>
1428 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1429  DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number>
1431 {
1433 
1434  for (unsigned int i = 0; i < dim; ++i)
1435  tmp.values[i] = -values[i];
1436 
1437  return tmp;
1438 }
1439 
1440 
1441 template <int rank_, int dim, typename Number>
1444 {
1445  return std::sqrt(norm_square());
1446 }
1447 
1448 
1449 template <int rank_, int dim, typename Number>
1450 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1453 {
1454  typename numbers::NumberTraits<Number>::real_type s = internal::NumberType<
1455  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1456  for (unsigned int i = 0; i < dim; ++i)
1457  s += values[i].norm_square();
1458 
1459  return s;
1460 }
1461 
1462 
1463 template <int rank_, int dim, typename Number>
1464 template <typename OtherNumber>
1465 inline void
1466 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
1467 {
1468  AssertDimension(result.size(),
1469  (Utilities::fixed_power<rank_, unsigned int>(dim)));
1470 
1471  unsigned int index = 0;
1472  unroll_recursion(result, index);
1473 }
1474 
1475 
1476 template <int rank_, int dim, typename Number>
1477 template <typename OtherNumber>
1478 inline void
1480  unsigned int & index) const
1481 {
1482  for (unsigned int i = 0; i < dim; ++i)
1483  values[i].unroll_recursion(result, index);
1484 }
1485 
1486 
1487 template <int rank_, int dim, typename Number>
1488 DEAL_II_CONSTEXPR inline unsigned int
1490  const TableIndices<rank_> &indices)
1491 {
1492  unsigned int index = 0;
1493  for (int r = 0; r < rank_; ++r)
1494  index = index * dim + indices[r];
1495 
1496  return index;
1497 }
1498 
1499 
1500 
1501 namespace internal
1502 {
1503  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1504  // and rank=2. Make sure we don't have compiler warnings.
1505 
1506  template <int dim>
1507  inline DEAL_II_CONSTEXPR unsigned int
1508  mod(const unsigned int x)
1509  {
1510  return x % dim;
1511  }
1512 
1513  template <>
1514  inline unsigned int
1515  mod<0>(const unsigned int x)
1516  {
1517  Assert(false, ExcInternalError());
1518  return x;
1519  }
1520 
1521  template <int dim>
1522  inline DEAL_II_CONSTEXPR unsigned int
1523  div(const unsigned int x)
1524  {
1525  return x / dim;
1526  }
1527 
1528  template <>
1529  inline unsigned int
1530  div<0>(const unsigned int x)
1531  {
1532  Assert(false, ExcInternalError());
1533  return x;
1534  }
1535 
1536 } // namespace internal
1537 
1538 
1539 
1540 template <int rank_, int dim, typename Number>
1541 DEAL_II_CONSTEXPR inline TableIndices<rank_>
1543 {
1544  AssertIndexRange(i, n_independent_components);
1545 
1546  TableIndices<rank_> indices;
1547 
1548  unsigned int remainder = i;
1549  for (int r = rank_ - 1; r >= 0; --r)
1550  {
1551  indices[r] = internal::mod<dim>(remainder);
1552  remainder = internal::div<dim>(remainder);
1553  }
1554  Assert(remainder == 0, ExcInternalError());
1555 
1556  return indices;
1557 }
1558 
1559 
1560 template <int rank_, int dim, typename Number>
1561 DEAL_II_CONSTEXPR inline void
1563 {
1564  for (unsigned int i = 0; i < dim; ++i)
1565  values[i] = internal::NumberType<Number>::value(0.0);
1566 }
1567 
1568 
1569 template <int rank_, int dim, typename Number>
1570 constexpr std::size_t
1572 {
1573  return sizeof(Tensor<rank_, dim, Number>);
1574 }
1575 
1576 
1577 template <int rank_, int dim, typename Number>
1578 template <class Archive>
1579 inline void
1580 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1581 {
1582  ar &values;
1583 }
1584 
1585 
1586 template <int rank, int dim, typename Number>
1588 
1589 
1590 /* ----------------- Non-member functions operating on tensors. ------------ */
1591 
1596 
1604 template <int rank_, int dim, typename Number>
1605 inline std::ostream &
1606 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1607 {
1608  for (unsigned int i = 0; i < dim; ++i)
1609  {
1610  out << p[i];
1611  if (i != dim - 1)
1612  out << ' ';
1613  }
1614 
1615  return out;
1616 }
1617 
1618 
1625 template <int dim, typename Number>
1626 inline std::ostream &
1627 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1628 {
1629  out << static_cast<const Number &>(p);
1630  return out;
1631 }
1632 
1633 
1635 
1639 
1640 
1651 template <int dim, typename Number, typename Other>
1652 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1653  typename ProductType<Other, Number>::type
1654  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1655 {
1656  return object * static_cast<const Number &>(t);
1657 }
1658 
1659 
1660 
1671 template <int dim, typename Number, typename Other>
1672 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1673  typename ProductType<Number, Other>::type
1674  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1675 {
1676  return static_cast<const Number &>(t) * object;
1677 }
1678 
1679 
1691 template <int dim, typename Number, typename OtherNumber>
1692 DEAL_II_CUDA_HOST_DEV constexpr DEAL_II_ALWAYS_INLINE
1693  typename ProductType<Number, OtherNumber>::type
1695  const Tensor<0, dim, OtherNumber> &src2)
1696 {
1697  return static_cast<const Number &>(src1) *
1698  static_cast<const OtherNumber &>(src2);
1699 }
1700 
1701 
1709 template <int dim, typename Number, typename OtherNumber>
1710 DEAL_II_CUDA_HOST_DEV constexpr DEAL_II_ALWAYS_INLINE
1711  Tensor<0,
1712  dim,
1713  typename ProductType<Number,
1714  typename EnableIfScalar<OtherNumber>::type>::type>
1715  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
1716 {
1717  return static_cast<const Number &>(t) / factor;
1718 }
1719 
1720 
1728 template <int dim, typename Number, typename OtherNumber>
1729 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1732  const Tensor<0, dim, OtherNumber> &q)
1733 {
1734  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1735 }
1736 
1737 
1745 template <int dim, typename Number, typename OtherNumber>
1746 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1749  const Tensor<0, dim, OtherNumber> &q)
1750 {
1751  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1752 }
1753 
1754 
1767 template <int rank, int dim, typename Number, typename OtherNumber>
1768 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1769  Tensor<rank,
1770  dim,
1771  typename ProductType<Number,
1772  typename EnableIfScalar<OtherNumber>::type>::type>
1773  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1774 {
1775  // recurse over the base objects
1777  for (unsigned int d = 0; d < dim; ++d)
1778  tt[d] = t[d] * factor;
1779  return tt;
1780 }
1781 
1782 
1795 template <int rank, int dim, typename Number, typename OtherNumber>
1796 DEAL_II_CUDA_HOST_DEV DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1797  Tensor<rank,
1798  dim,
1800  OtherNumber>::type>
1801  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
1802 {
1803  // simply forward to the operator above
1804  return t * factor;
1805 }
1806 
1807 
1808 namespace internal
1809 {
1810  namespace TensorImplementation
1811  {
1812  template <int rank,
1813  int dim,
1814  typename Number,
1815  typename OtherNumber,
1816  typename std::enable_if<
1817  !std::is_integral<
1818  typename ProductType<Number, OtherNumber>::type>::value,
1819  int>::type = 0>
1820  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1822  division_operator(const Tensor<rank, dim, Number> &t,
1823  const OtherNumber & factor)
1824  {
1826  const Number inverse_factor = Number(1.) / factor;
1827  // recurse over the base objects
1828  for (unsigned int d = 0; d < dim; ++d)
1829  tt[d] = t[d] * inverse_factor;
1830  return tt;
1831  }
1832 
1833 
1834  template <int rank,
1835  int dim,
1836  typename Number,
1837  typename OtherNumber,
1838  typename std::enable_if<
1839  std::is_integral<
1840  typename ProductType<Number, OtherNumber>::type>::value,
1841  int>::type = 0>
1842  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1844  division_operator(const Tensor<rank, dim, Number> &t,
1845  const OtherNumber & factor)
1846  {
1848  // recurse over the base objects
1849  for (unsigned int d = 0; d < dim; ++d)
1850  tt[d] = t[d] / factor;
1851  return tt;
1852  }
1853  } // namespace TensorImplementation
1854 } // namespace internal
1855 
1856 
1866 template <int rank, int dim, typename Number, typename OtherNumber>
1867 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1868  Tensor<rank,
1869  dim,
1870  typename ProductType<Number,
1871  typename EnableIfScalar<OtherNumber>::type>::type>
1872  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1873 {
1874  return internal::TensorImplementation::division_operator(t, factor);
1875 }
1876 
1877 
1887 template <int rank, int dim, typename Number, typename OtherNumber>
1888 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1892 {
1894 
1895  for (unsigned int i = 0; i < dim; ++i)
1896  tmp[i] += q[i];
1897 
1898  return tmp;
1899 }
1900 
1901 
1911 template <int rank, int dim, typename Number, typename OtherNumber>
1912 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1916 {
1918 
1919  for (unsigned int i = 0; i < dim; ++i)
1920  tmp[i] -= q[i];
1921 
1922  return tmp;
1923 }
1924 
1931 template <int dim, typename Number, typename OtherNumber>
1932 inline DEAL_II_CONSTEXPR DEAL_II_ALWAYS_INLINE
1935  const Tensor<0, dim, OtherNumber> &src2)
1936 {
1938 
1939  tmp *= src2;
1940 
1941  return tmp;
1942 }
1943 
1949 template <int dim, typename Number, typename OtherNumber>
1950 inline DEAL_II_CONSTEXPR DEAL_II_ALWAYS_INLINE
1953  const Tensor<1, dim, OtherNumber> &src2)
1954 {
1956 
1957  for (unsigned int i = 0; i < dim; ++i)
1958  tmp[i] *= src2[i];
1959 
1960  return tmp;
1961 }
1962 
1979 template <int rank, int dim, typename Number, typename OtherNumber>
1980 inline DEAL_II_CONSTEXPR DEAL_II_ALWAYS_INLINE
1983  const Tensor<rank, dim, OtherNumber> &src2)
1984 {
1986 
1987  for (unsigned int i = 0; i < dim; ++i)
1988  tmp[i] = schur_product(src1[i], src2[i]);
1989 
1990  return tmp;
1991 }
1992 
1994 
1998 
1999 
2023 template <int rank_1,
2024  int rank_2,
2025  int dim,
2026  typename Number,
2027  typename OtherNumber>
2028 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2029  typename Tensor<rank_1 + rank_2 - 2,
2030  dim,
2031  typename ProductType<Number, OtherNumber>::type>::tensor_type
2034 {
2035  typename Tensor<rank_1 + rank_2 - 2,
2036  dim,
2037  typename ProductType<Number, OtherNumber>::type>::tensor_type
2038  result{};
2039 
2040  TensorAccessors::internal::
2041  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2042  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2043  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2044 
2045  return result;
2046 }
2047 
2048 
2078 template <int index_1,
2079  int index_2,
2080  int rank_1,
2081  int rank_2,
2082  int dim,
2083  typename Number,
2084  typename OtherNumber>
2085 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2086  typename Tensor<rank_1 + rank_2 - 2,
2087  dim,
2088  typename ProductType<Number, OtherNumber>::type>::tensor_type
2091 {
2092  Assert(0 <= index_1 && index_1 < rank_1,
2093  ExcMessage(
2094  "The specified index_1 must lie within the range [0,rank_1)"));
2095  Assert(0 <= index_2 && index_2 < rank_2,
2096  ExcMessage(
2097  "The specified index_2 must lie within the range [0,rank_2)"));
2098 
2099  using namespace TensorAccessors;
2100  using namespace TensorAccessors::internal;
2101 
2102  // Reorder index_1 to the end of src1:
2103  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>
2104  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2105 
2106  // Reorder index_2 to the end of src2:
2107  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2108  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2109 
2110  typename Tensor<rank_1 + rank_2 - 2,
2111  dim,
2112  typename ProductType<Number, OtherNumber>::type>::tensor_type
2113  result{};
2114  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2115  return result;
2116 }
2117 
2118 
2150 template <int index_1,
2151  int index_2,
2152  int index_3,
2153  int index_4,
2154  int rank_1,
2155  int rank_2,
2156  int dim,
2157  typename Number,
2158  typename OtherNumber>
2159 DEAL_II_CONSTEXPR inline
2160  typename Tensor<rank_1 + rank_2 - 4,
2161  dim,
2162  typename ProductType<Number, OtherNumber>::type>::tensor_type
2165 {
2166  Assert(0 <= index_1 && index_1 < rank_1,
2167  ExcMessage(
2168  "The specified index_1 must lie within the range [0,rank_1)"));
2169  Assert(0 <= index_3 && index_3 < rank_1,
2170  ExcMessage(
2171  "The specified index_3 must lie within the range [0,rank_1)"));
2172  Assert(index_1 != index_3,
2173  ExcMessage("index_1 and index_3 must not be the same"));
2174  Assert(0 <= index_2 && index_2 < rank_2,
2175  ExcMessage(
2176  "The specified index_2 must lie within the range [0,rank_2)"));
2177  Assert(0 <= index_4 && index_4 < rank_2,
2178  ExcMessage(
2179  "The specified index_4 must lie within the range [0,rank_2)"));
2180  Assert(index_2 != index_4,
2181  ExcMessage("index_2 and index_4 must not be the same"));
2182 
2183  using namespace TensorAccessors;
2184  using namespace TensorAccessors::internal;
2185 
2186  // Reorder index_1 to the end of src1:
2187  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>
2188  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2189 
2190  // Reorder index_2 to the end of src2:
2191  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2192  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2193 
2194  // Now, reorder index_3 to the end of src1. We have to make sure to
2195  // preserve the original ordering: index_1 has been removed. If
2196  // index_3 > index_1, we have to use (index_3 - 1) instead:
2197  ReorderedIndexView<
2198  (index_3 < index_1 ? index_3 : index_3 - 1),
2199  rank_1,
2200  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2201  reord_3 =
2202  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2203  index_3 - 1,
2204  rank_1 > (reord_1);
2205 
2206  // Now, reorder index_4 to the end of src2. We have to make sure to
2207  // preserve the original ordering: index_2 has been removed. If
2208  // index_4 > index_2, we have to use (index_4 - 1) instead:
2209  ReorderedIndexView<
2210  (index_4 < index_2 ? index_4 : index_4 - 1),
2211  rank_2,
2212  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2213  reord_4 =
2214  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2215  index_4 - 1,
2216  rank_2 > (reord_2);
2217 
2218  typename Tensor<rank_1 + rank_2 - 4,
2219  dim,
2220  typename ProductType<Number, OtherNumber>::type>::tensor_type
2221  result{};
2222  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2223  return result;
2224 }
2225 
2226 
2240 template <int rank, int dim, typename Number, typename OtherNumber>
2241 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2242  typename ProductType<Number, OtherNumber>::type
2244  const Tensor<rank, dim, OtherNumber> &right)
2245 {
2246  typename ProductType<Number, OtherNumber>::type result{};
2247  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2248  return result;
2249 }
2250 
2251 
2270 template <template <int, int, typename> class TensorT1,
2271  template <int, int, typename> class TensorT2,
2272  template <int, int, typename> class TensorT3,
2273  int rank_1,
2274  int rank_2,
2275  int dim,
2276  typename T1,
2277  typename T2,
2278  typename T3>
2279 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2281  contract3(const TensorT1<rank_1, dim, T1> & left,
2282  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2283  const TensorT3<rank_2, dim, T3> & right)
2284 {
2285  using return_type =
2287  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2288  middle,
2289  right);
2290 }
2291 
2292 
2304 template <int rank_1,
2305  int rank_2,
2306  int dim,
2307  typename Number,
2308  typename OtherNumber>
2309 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2313 {
2314  typename Tensor<rank_1 + rank_2,
2315  dim,
2316  typename ProductType<Number, OtherNumber>::type>::tensor_type
2317  result{};
2318  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2319  return result;
2320 }
2321 
2322 
2324 
2328 
2329 
2341 template <int dim, typename Number>
2342 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
2344 {
2345  Assert(dim == 2, ExcInternalError());
2346 
2347  Tensor<1, dim, Number> result;
2348 
2349  result[0] = src[1];
2350  result[1] = -src[0];
2351 
2352  return result;
2353 }
2354 
2355 
2366 template <int dim, typename Number1, typename Number2>
2367 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2370  const Tensor<1, dim, Number2> &src2)
2371 {
2372  Assert(dim == 3, ExcInternalError());
2373 
2375 
2376  result[0] = src1[1] * src2[2] - src1[2] * src2[1];
2377  result[1] = src1[2] * src2[0] - src1[0] * src2[2];
2378  result[2] = src1[0] * src2[1] - src1[1] * src2[0];
2379 
2380  return result;
2381 }
2382 
2383 
2385 
2389 
2390 
2397 template <int dim, typename Number>
2398 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
2400 {
2401  // Compute the determinant using the Laplace expansion of the
2402  // determinant. We expand along the last row.
2403  Number det = internal::NumberType<Number>::value(0.0);
2404 
2405  for (unsigned int k = 0; k < dim; ++k)
2406  {
2407  Tensor<2, dim - 1, Number> minor;
2408  for (unsigned int i = 0; i < dim - 1; ++i)
2409  for (unsigned int j = 0; j < dim - 1; ++j)
2410  minor[i][j] = t[i][j < k ? j : j + 1];
2411 
2412  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2413 
2414  det += t[dim - 1][k] * cofactor;
2415  }
2416 
2417  return ((dim % 2 == 0) ? 1. : -1.) * det;
2418 }
2419 
2425 template <typename Number>
2426 constexpr DEAL_II_ALWAYS_INLINE Number
2428 {
2429  return t[0][0];
2430 }
2431 
2432 
2440 template <int dim, typename Number>
2441 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
2443 {
2444  Number t = d[0][0];
2445  for (unsigned int i = 1; i < dim; ++i)
2446  t += d[i][i];
2447  return t;
2448 }
2449 
2450 
2460 template <int dim, typename Number>
2461 DEAL_II_CONSTEXPR inline Tensor<2, dim, Number>
2463 {
2464  Number return_tensor[dim][dim];
2465 
2466  // if desired, take over the
2467  // inversion of a 4x4 tensor
2468  // from the FullMatrix
2469  AssertThrow(false, ExcNotImplemented());
2470 
2471  return Tensor<2, dim, Number>(return_tensor);
2472 }
2473 
2474 
2475 #ifndef DOXYGEN
2476 
2477 template <typename Number>
2478 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<2, 1, Number>
2479  invert(const Tensor<2, 1, Number> &t)
2480 {
2481  Tensor<2, 1, Number> return_tensor;
2482 
2483  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2484 
2485  return return_tensor;
2486 }
2487 
2488 
2489 template <typename Number>
2490 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<2, 2, Number>
2491  invert(const Tensor<2, 2, Number> &t)
2492 {
2493  Tensor<2, 2, Number> return_tensor;
2494 
2495  // this is Maple output,
2496  // thus a bit unstructured
2497  const Number inv_det_t = internal::NumberType<Number>::value(
2498  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2499  return_tensor[0][0] = t[1][1];
2500  return_tensor[0][1] = -t[0][1];
2501  return_tensor[1][0] = -t[1][0];
2502  return_tensor[1][1] = t[0][0];
2503  return_tensor *= inv_det_t;
2504 
2505  return return_tensor;
2506 }
2507 
2508 
2509 template <typename Number>
2510 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<2, 3, Number>
2511  invert(const Tensor<2, 3, Number> &t)
2512 {
2513  Tensor<2, 3, Number> return_tensor;
2514 
2515  const Number t4 = internal::NumberType<Number>::value(t[0][0] * t[1][1]),
2516  t6 = internal::NumberType<Number>::value(t[0][0] * t[1][2]),
2517  t8 = internal::NumberType<Number>::value(t[0][1] * t[1][0]),
2518  t00 = internal::NumberType<Number>::value(t[0][2] * t[1][0]),
2519  t01 = internal::NumberType<Number>::value(t[0][1] * t[2][0]),
2520  t04 = internal::NumberType<Number>::value(t[0][2] * t[2][0]),
2521  inv_det_t = internal::NumberType<Number>::value(
2522  1.0 / (t4 * t[2][2] - t6 * t[2][1] - t8 * t[2][2] +
2523  t00 * t[2][1] + t01 * t[1][2] - t04 * t[1][1]));
2524  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2525  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2526  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2527  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2528  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2529  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2530  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2531  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2532  return_tensor[1][1] =
2533  internal::NumberType<Number>::value(t[0][0] * t[2][2]) - t04;
2534  return_tensor[1][2] = t00 - t6;
2535  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2536  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2537  return_tensor[2][1] =
2538  t01 - internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2539  return_tensor[2][2] = internal::NumberType<Number>::value(t4 - t8);
2540  return_tensor *= inv_det_t;
2541 
2542  return return_tensor;
2543 }
2544 
2545 #endif /* DOXYGEN */
2546 
2547 
2554 template <int dim, typename Number>
2555 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, Number>
2557 {
2559  for (unsigned int i = 0; i < dim; ++i)
2560  {
2561  tt[i][i] = t[i][i];
2562  for (unsigned int j = i + 1; j < dim; ++j)
2563  {
2564  tt[i][j] = t[j][i];
2565  tt[j][i] = t[i][j];
2566  };
2567  }
2568  return tt;
2569 }
2570 
2571 
2586 template <int dim, typename Number>
2587 constexpr Tensor<2, dim, Number>
2589 {
2590  return determinant(t) * invert(t);
2591 }
2592 
2593 
2608 template <int dim, typename Number>
2609 constexpr Tensor<2, dim, Number>
2611 {
2612  return transpose(adjugate(t));
2613 }
2614 
2615 
2623 template <int dim, typename Number>
2624 inline Number
2626 {
2627  Number max = internal::NumberType<Number>::value(0.0);
2628  for (unsigned int j = 0; j < dim; ++j)
2629  {
2630  Number sum = internal::NumberType<Number>::value(0.0);
2631  for (unsigned int i = 0; i < dim; ++i)
2632  sum += std::fabs(t[i][j]);
2633 
2634  if (sum > max)
2635  max = sum;
2636  }
2637 
2638  return max;
2639 }
2640 
2641 
2649 template <int dim, typename Number>
2650 inline Number
2652 {
2653  Number max = internal::NumberType<Number>::value(0.0);
2654  for (unsigned int i = 0; i < dim; ++i)
2655  {
2656  Number sum = internal::NumberType<Number>::value(0.0);
2657  for (unsigned int j = 0; j < dim; ++j)
2658  sum += std::fabs(t[i][j]);
2659 
2660  if (sum > max)
2661  max = sum;
2662  }
2663 
2664  return max;
2665 }
2666 
2668 
2669 
2670 #ifndef DOXYGEN
2671 
2672 
2673 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2674 
2675 // Specialization of functions for ADOL-C number types when
2676 // the advanced branching feature is used
2677 template <int dim>
2678 inline adouble
2680 {
2681  adouble max = internal::NumberType<adouble>::value(0.0);
2682  for (unsigned int j = 0; j < dim; ++j)
2683  {
2684  adouble sum = internal::NumberType<adouble>::value(0.0);
2685  for (unsigned int i = 0; i < dim; ++i)
2686  sum += std::fabs(t[i][j]);
2687 
2688  condassign(max, (sum > max), sum, max);
2689  }
2690 
2691  return max;
2692 }
2693 
2694 
2695 template <int dim>
2696 inline adouble
2698 {
2699  adouble max = internal::NumberType<adouble>::value(0.0);
2700  for (unsigned int i = 0; i < dim; ++i)
2701  {
2702  adouble sum = internal::NumberType<adouble>::value(0.0);
2703  for (unsigned int j = 0; j < dim; ++j)
2704  sum += std::fabs(t[i][j]);
2705 
2706  condassign(max, (sum > max), sum, max);
2707  }
2708 
2709  return max;
2710 }
2711 
2712 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2713 
2714 
2715 #endif // DOXYGEN
2716 
2717 DEAL_II_NAMESPACE_CLOSE
2718 
2719 #endif
constexpr Tensor()=default
DEAL_II_CONSTEXPR Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > schur_product(const Tensor< rank, dim, Number > &src1, const Tensor< rank, dim, OtherNumber > &src2)
Definition: tensor.h:1982
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
DEAL_II_CONSTEXPR ProductType< Number, Other >::type operator*(const Tensor< 0, dim, Number > &t, const Other &object)
Definition: tensor.h:1674
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1571
DEAL_II_CONSTEXPR ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:2243
DEAL_II_CONSTEXPR ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1654
DEAL_II_CONSTEXPR Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2556
constexpr Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2610
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:943
DEAL_II_CONSTEXPR Tensor< rank_1+rank_2 - 2, dim, typename ProductType< Number, OtherNumber >::type >::tensor_type contract(const Tensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
Definition: tensor.h:2089
DEAL_II_CONSTEXPR Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:2462
DEAL_II_CONSTEXPR Tensor & operator/=(const OtherNumber &factor)
static constexpr std::enable_if< std::is_same< Dummy, number >::value &&is_cuda_compatible< Dummy >::value, real_type >::type abs_square(const number &x)
Definition: numbers.h:605
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2625
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
DEAL_II_CONSTEXPR Tensor< rank, dim, typename ProductType< typename EnableIfScalar< Number >::type, OtherNumber >::type > operator*(const Number &factor, const Tensor< rank, dim, OtherNumber > &t)
Definition: tensor.h:1801
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2651
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:724
DEAL_II_CONSTEXPR Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
STL namespace.
typename Tensor< rank_ - 1, dim, Number >::array_type[(dim !=0) ? dim :1] array_type
Definition: tensor.h:459
DEAL_II_CONSTEXPR Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2369
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
static real_type abs(const number &x)
Definition: numbers.h:627
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1443
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: tensor.h:1479
Definition: point.h:110
LinearAlgebra::distributed::Vector< Number > Vector
DEAL_II_CONSTEXPR Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< rank, dim, Number > &p, const Tensor< rank, dim, OtherNumber > &q)
Definition: tensor.h:1890
DEAL_II_CONSTEXPR void clear()
Definition: tensor.h:1562
DEAL_II_CONSTEXPR Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< rank, dim, Number > &p, const Tensor< rank, dim, OtherNumber > &q)
Definition: tensor.h:1914
DEAL_II_CONSTEXPR value_type & operator[](const unsigned int i)
Definition: tensor.h:1190
DEAL_II_CONSTEXPR Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2399
static ::ExceptionBase & ExcMessage(std::string arg1)
DEAL_II_CONSTEXPR Tensor & operator*=(const OtherNumber &factor)
DEAL_II_CONSTEXPR Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2343
constexpr Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
Definition: tensor.h:1748
DEAL_II_CONSTEXPR Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1872
DEAL_II_CONSTEXPR Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
DEAL_II_CONSTEXPR ProductType< T1, typename ProductType< T2, T3 >::type >::type contract3(const TensorT1< rank_1, dim, T1 > &left, const TensorT2< rank_1+rank_2, dim, T2 > &middle, const TensorT3< rank_2, dim, T3 > &right)
Definition: tensor.h:2281
#define Assert(cond, exc)
Definition: exceptions.h:1411
constexpr internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
void serialize(Archive &ar, const unsigned int version)
Definition: tensor.h:1580
static DEAL_II_CONSTEXPR unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
Definition: tensor.h:1489
Number * end_raw()
Definition: tensor.h:1253
typename Tensor< rank_ - 1, dim, Number >::tensor_type value_type
Definition: tensor.h:452
DEAL_II_CONSTEXPR Tensor< rank_1+rank_2 - 2, dim, typename ProductType< Number, OtherNumber >::type >::tensor_type operator*(const Tensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
Definition: tensor.h:2032
static constexpr unsigned int rank
Definition: tensor.h:438
constexpr Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2588
constexpr ProductType< Number, OtherNumber >::type operator*(const Tensor< 0, dim, Number > &src1, const Tensor< 0, dim, OtherNumber > &src2)
Definition: tensor.h:1694
DEAL_II_CONSTEXPR Tensor< rank_1+rank_2, dim, typename ProductType< Number, OtherNumber >::type > outer_product(const Tensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
Definition: tensor.h:2311
DEAL_II_CONSTEXPR Tensor operator-() const
Definition: tensor.h:1430
constexpr bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1327
typename numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:132
DEAL_II_CONSTEXPR Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > schur_product(const Tensor< 0, dim, Number > &src1, const Tensor< 0, dim, OtherNumber > &src2)
Definition: tensor.h:1934
DEAL_II_CONSTEXPR numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1452
constexpr Tensor< 0, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const Tensor< 0, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1715
Definition: tensor.h:422
Expression operator==(const Expression &lhs, const Expression &rhs)
constexpr Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
Definition: tensor.h:1731
DEAL_II_CONSTEXPR Tensor & operator+=(const Tensor< rank_, dim, OtherNumber > &)
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:959
Number * begin_raw()
Definition: tensor.h:1233
DEAL_II_CONSTEXPR Tensor< rank_1+rank_2 - 4, dim, typename ProductType< Number, OtherNumber >::type >::tensor_type double_contract(const Tensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
Definition: tensor.h:2163
Expression operator-(Expression lhs, const Expression &rhs)
static ::ExceptionBase & ExcNotImplemented()
void unroll(Vector< OtherNumber > &result) const
Definition: tensor.h:1466
static DEAL_II_CONSTEXPR TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Definition: tensor.h:1542
static constexpr unsigned int n_independent_components
Definition: tensor.h:444
DEAL_II_CONSTEXPR bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1300
Expression operator!=(const Expression &lhs, const Expression &rhs)
DEAL_II_CONSTEXPR Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2442
static constexpr std::size_t memory_consumption()
Definition: tensor.h:1571
constexpr Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:2427
DEAL_II_CONSTEXPR Tensor< 1, dim, typename ProductType< Number, OtherNumber >::type > schur_product(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, OtherNumber > &src2)
Definition: tensor.h:1952
static ::ExceptionBase & ExcInternalError()
DEAL_II_CONSTEXPR Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1773