Reference documentation for deal.II version Git 1ecc23629d 2021-05-18 09:57:04 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tensor.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2021 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 
22 #include <deal.II/base/numbers.h>
26 #include <deal.II/base/utilities.h>
27 
28 #ifdef DEAL_II_WITH_ADOLC
29 # include <adolc/adouble.h> // Taped double
30 #endif
31 
32 #include <cmath>
33 #include <ostream>
34 #include <utility>
35 #include <vector>
36 
37 
39 
40 // Forward declarations:
41 #ifndef DOXYGEN
42 template <typename ElementType, typename MemorySpace>
43 class ArrayView;
44 template <int dim, typename Number>
45 class Point;
46 template <int rank_, int dim, typename Number = double>
47 class Tensor;
48 template <typename Number>
49 class Vector;
50 template <typename number>
51 class FullMatrix;
52 namespace Differentiation
53 {
54  namespace SD
55  {
56  class Expression;
57  }
58 } // namespace Differentiation
59 #endif
60 
61 
91 template <int dim, typename Number>
92 class Tensor<0, dim, Number>
93 {
94 public:
95  static_assert(dim >= 0,
96  "Tensors must have a dimension greater than or equal to one.");
97 
106  static constexpr unsigned int dimension = dim;
107 
111  static constexpr unsigned int rank = 0;
112 
116  static constexpr unsigned int n_independent_components = 1;
117 
127 
132  using value_type = Number;
133 
139  using array_type = Number;
140 
146  constexpr DEAL_II_CUDA_HOST_DEV
147  Tensor();
148 
156  template <typename OtherNumber>
157  constexpr DEAL_II_CUDA_HOST_DEV
158  Tensor(const Tensor<0, dim, OtherNumber> &initializer);
159 
165  template <typename OtherNumber>
166  constexpr DEAL_II_CUDA_HOST_DEV
167  Tensor(const OtherNumber &initializer);
168 
172  Number *
173  begin_raw();
174 
178  const Number *
179  begin_raw() const;
180 
184  Number *
185  end_raw();
186 
191  const Number *
192  end_raw() const;
193 
203  constexpr DEAL_II_CUDA_HOST_DEV operator Number &();
204 
213  constexpr DEAL_II_CUDA_HOST_DEV operator const Number &() const;
214 
222  template <typename OtherNumber>
223  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
224  operator=(const Tensor<0, dim, OtherNumber> &rhs);
225 
226 #ifdef __INTEL_COMPILER
227 
235  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
236  operator=(const Tensor<0, dim, Number> &rhs);
237 #endif
238 
245  template <typename OtherNumber>
246  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
247  operator=(const OtherNumber &d);
248 
252  template <typename OtherNumber>
253  constexpr bool
254  operator==(const Tensor<0, dim, OtherNumber> &rhs) const;
255 
259  template <typename OtherNumber>
260  constexpr bool
261  operator!=(const Tensor<0, dim, OtherNumber> &rhs) const;
262 
268  template <typename OtherNumber>
269  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
270  operator+=(const Tensor<0, dim, OtherNumber> &rhs);
271 
277  template <typename OtherNumber>
278  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
279  operator-=(const Tensor<0, dim, OtherNumber> &rhs);
280 
286  template <typename OtherNumber>
287  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
288  operator*=(const OtherNumber &factor);
289 
295  template <typename OtherNumber>
296  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
297  operator/=(const OtherNumber &factor);
298 
304  constexpr DEAL_II_CUDA_HOST_DEV Tensor
305  operator-() const;
306 
319  constexpr void
320  clear();
321 
327  real_type
328  norm() const;
329 
337  norm_square() const;
338 
344  template <class Archive>
345  void
346  serialize(Archive &ar, const unsigned int version);
347 
352  using tensor_type = Number;
353 
354 private:
358  Number value;
359 
363  template <typename OtherNumber>
364  void
365  unroll_recursion(Vector<OtherNumber> &result,
366  unsigned int & start_index) const;
367 
368  // Allow an arbitrary Tensor to access the underlying values.
369  template <int, int, typename>
370  friend class Tensor;
371 };
372 
373 
374 
448 template <int rank_, int dim, typename Number>
449 class Tensor
450 {
451 public:
452  static_assert(rank_ >= 1,
453  "Tensors must have a rank greater than or equal to one.");
454  static_assert(dim >= 0,
455  "Tensors must have a dimension greater than or equal to one.");
464  static constexpr unsigned int dimension = dim;
465 
469  static constexpr unsigned int rank = rank_;
470 
475  static constexpr unsigned int n_independent_components =
476  Tensor<rank_ - 1, dim>::n_independent_components * dim;
477 
483  using value_type = typename Tensor<rank_ - 1, dim, Number>::tensor_type;
484 
489  using array_type =
490  typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1];
491 
498  Tensor();
499 
505  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
506  const array_type &initializer);
507 
521  template <typename ElementType, typename MemorySpace>
522  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
523  const ArrayView<ElementType, MemorySpace> &initializer);
524 
532  template <typename OtherNumber>
533  constexpr DEAL_II_CUDA_HOST_DEV
534  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
535 
539  template <typename OtherNumber>
540  constexpr Tensor(
541  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
542 
546  template <typename OtherNumber>
547  constexpr
548  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
549 
555  constexpr DEAL_II_CUDA_HOST_DEV value_type &operator[](const unsigned int i);
556 
562  constexpr DEAL_II_CUDA_HOST_DEV const value_type &
563  operator[](const unsigned int i) const;
564 
568  constexpr const Number &operator[](const TableIndices<rank_> &indices) const;
569 
573  constexpr Number &operator[](const TableIndices<rank_> &indices);
574 
578  Number *
579  begin_raw();
580 
584  const Number *
585  begin_raw() const;
586 
590  Number *
591  end_raw();
592 
596  const Number *
597  end_raw() const;
598 
606  template <typename OtherNumber>
607  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
608  operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
609 
616  constexpr Tensor &
617  operator=(const Number &d);
618 
622  template <typename OtherNumber>
623  constexpr bool
625 
629  template <typename OtherNumber>
630  constexpr bool
632 
638  template <typename OtherNumber>
639  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
640  operator+=(const Tensor<rank_, dim, OtherNumber> &);
641 
647  template <typename OtherNumber>
648  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
649  operator-=(const Tensor<rank_, dim, OtherNumber> &);
650 
657  template <typename OtherNumber>
658  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
659  operator*=(const OtherNumber &factor);
660 
666  template <typename OtherNumber>
667  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
668  operator/=(const OtherNumber &factor);
669 
675  constexpr DEAL_II_CUDA_HOST_DEV Tensor
676  operator-() const;
677 
690  constexpr void
691  clear();
692 
702  norm() const;
703 
710  constexpr DEAL_II_CUDA_HOST_DEV
712  norm_square() const;
713 
721  template <typename OtherNumber>
722  void
723  unroll(Vector<OtherNumber> &result) const;
724 
729  static constexpr unsigned int
730  component_to_unrolled_index(const TableIndices<rank_> &indices);
731 
737  static constexpr TableIndices<rank_>
738  unrolled_to_component_indices(const unsigned int i);
739 
744  static constexpr std::size_t
746 
752  template <class Archive>
753  void
754  serialize(Archive &ar, const unsigned int version);
755 
761 
762 private:
766  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
767  // ... avoid a compiler warning in case of dim == 0 and ensure that the
768  // array always has positive size.
769 
773  template <typename OtherNumber>
774  void
775  unroll_recursion(Vector<OtherNumber> &result,
776  unsigned int & start_index) const;
777 
784  template <typename ArrayLike, std::size_t... Indices>
785  constexpr DEAL_II_CUDA_HOST_DEV
786  Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
787 
788  // Allow an arbitrary Tensor to access the underlying values.
789  template <int, int, typename>
790  friend class Tensor;
791 
792  // Point is allowed access to the coordinates. This is supposed to improve
793  // speed.
794  friend class Point<dim, Number>;
795 };
796 
797 
798 #ifndef DOXYGEN
799 namespace internal
800 {
801  // Workaround: The following 4 overloads are necessary to be able to
802  // compile the library with Apple Clang 8 and older. We should remove
803  // these overloads again when we bump the minimal required version to
804  // something later than clang-3.6 / Apple Clang 6.3.
805  template <int rank, int dim, typename T, typename U>
806  struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
807  {
808  using type =
810  };
811 
812  template <int rank, int dim, typename T, typename U>
813  struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
814  {
815  using type =
817  };
818 
819  template <typename T, int rank, int dim, typename U>
820  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
821  {
822  using type =
824  };
825 
826  template <int rank, int dim, typename T, typename U>
827  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
828  {
829  using type =
831  };
832  // end workaround
833 
838  template <int rank, int dim, typename T>
839  struct NumberType<Tensor<rank, dim, T>>
840  {
841  static constexpr DEAL_II_ALWAYS_INLINE const Tensor<rank, dim, T> &
842  value(const Tensor<rank, dim, T> &t)
843  {
844  return t;
845  }
846 
848  value(const T &t)
849  {
851  tmp = t;
852  return tmp;
853  }
854  };
855 } // namespace internal
856 
857 
858 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
859 
860 
861 template <int dim, typename Number>
864  // Some auto-differentiable numbers need explicit
865  // zero initialization such as adtl::adouble.
866  : Tensor{0.0}
867 {}
868 
869 
870 
871 template <int dim, typename Number>
872 template <typename OtherNumber>
874  Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
875  : value(internal::NumberType<Number>::value(initializer))
876 {}
877 
878 
879 
880 template <int dim, typename Number>
881 template <typename OtherNumber>
884  : Tensor{p.value}
885 {}
886 
887 
888 
889 template <int dim, typename Number>
890 inline Number *
892 {
893  return std::addressof(value);
894 }
895 
896 
897 
898 template <int dim, typename Number>
899 inline const Number *
901 {
902  return std::addressof(value);
903 }
904 
905 
906 
907 template <int dim, typename Number>
908 inline Number *
910 {
912 }
913 
914 
915 
916 template <int dim, typename Number>
917 const Number *
919 {
921 }
922 
923 
924 
925 template <int dim, typename Number>
926 constexpr inline DEAL_II_ALWAYS_INLINE
928 {
929  // We cannot use Assert inside a CUDA kernel
930 # ifndef __CUDA_ARCH__
931  Assert(dim != 0,
932  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
933 # endif
934  return value;
935 }
936 
937 
938 template <int dim, typename Number>
939 constexpr inline DEAL_II_ALWAYS_INLINE
941 {
942  // We cannot use Assert inside a CUDA kernel
943 # ifndef __CUDA_ARCH__
944  Assert(dim != 0,
945  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
946 # endif
947  return value;
948 }
949 
950 
951 template <int dim, typename Number>
952 template <typename OtherNumber>
953 constexpr inline DEAL_II_ALWAYS_INLINE
956 {
958  return *this;
959 }
960 
961 
962 # ifdef __INTEL_COMPILER
963 template <int dim, typename Number>
964 constexpr inline DEAL_II_ALWAYS_INLINE
967 {
968  value = p.value;
969  return *this;
970 }
971 # endif
972 
973 
974 template <int dim, typename Number>
975 template <typename OtherNumber>
976 constexpr inline DEAL_II_ALWAYS_INLINE
978  Tensor<0, dim, Number>::operator=(const OtherNumber &d)
979 {
981  return *this;
982 }
983 
984 
985 template <int dim, typename Number>
986 template <typename OtherNumber>
987 constexpr inline bool
989 {
990 # if defined(DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING)
991  Assert(!(std::is_same<Number, adouble>::value ||
992  std::is_same<OtherNumber, adouble>::value),
993  ExcMessage(
994  "The Tensor equality operator for ADOL-C taped numbers has not yet "
995  "been extended to support advanced branching."));
996 # endif
997 
998  return numbers::values_are_equal(value, p.value);
999 }
1000 
1001 
1002 template <int dim, typename Number>
1003 template <typename OtherNumber>
1004 constexpr bool
1006 {
1007  return !((*this) == p);
1008 }
1009 
1010 
1011 template <int dim, typename Number>
1012 template <typename OtherNumber>
1013 constexpr inline DEAL_II_ALWAYS_INLINE
1016 {
1017  value += p.value;
1018  return *this;
1019 }
1020 
1021 
1022 template <int dim, typename Number>
1023 template <typename OtherNumber>
1024 constexpr inline DEAL_II_ALWAYS_INLINE
1027 {
1028  value -= p.value;
1029  return *this;
1030 }
1031 
1032 
1033 
1034 namespace internal
1035 {
1036  namespace ComplexWorkaround
1037  {
1038  template <typename Number, typename OtherNumber>
1039  constexpr inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV void
1040  multiply_assign_scalar(Number &val, const OtherNumber &s)
1041  {
1042  val *= s;
1043  }
1044 
1045 # ifdef __CUDA_ARCH__
1046  template <typename Number, typename OtherNumber>
1047  constexpr inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV void
1048  multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
1049  {
1050  printf("This function is not implemented for std::complex<Number>!\n");
1051  assert(false);
1052  }
1053 # endif
1054  } // namespace ComplexWorkaround
1055 } // namespace internal
1056 
1057 
1058 template <int dim, typename Number>
1059 template <typename OtherNumber>
1060 constexpr inline DEAL_II_ALWAYS_INLINE
1062  Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1063 {
1064  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1065  return *this;
1066 }
1067 
1068 
1069 
1070 template <int dim, typename Number>
1071 template <typename OtherNumber>
1073 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1074 {
1075  value /= s;
1076  return *this;
1077 }
1078 
1079 
1080 template <int dim, typename Number>
1083 {
1084  return -value;
1085 }
1086 
1087 
1088 template <int dim, typename Number>
1091 {
1092  Assert(dim != 0,
1093  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1094  return numbers::NumberTraits<Number>::abs(value);
1095 }
1096 
1097 
1098 template <int dim, typename Number>
1102 {
1103  // We cannot use Assert inside a CUDA kernel
1104 # ifndef __CUDA_ARCH__
1105  Assert(dim != 0,
1106  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1107 # endif
1109 }
1110 
1111 
1112 template <int dim, typename Number>
1113 template <typename OtherNumber>
1114 inline void
1115 Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1116  unsigned int & index) const
1117 {
1118  Assert(dim != 0,
1119  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1120  result[index] = value;
1121  ++index;
1122 }
1123 
1124 
1125 template <int dim, typename Number>
1126 constexpr inline void
1128 {
1129  // Some auto-differentiable numbers need explicit
1130  // zero initialization.
1132 }
1133 
1134 
1135 template <int dim, typename Number>
1136 template <class Archive>
1137 inline void
1138 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1139 {
1140  ar &value;
1141 }
1142 
1143 
1144 template <int dim, typename Number>
1146 
1147 
1148 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1149 
1150 template <int rank_, int dim, typename Number>
1151 template <typename ArrayLike, std::size_t... indices>
1153  Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1154  std::index_sequence<indices...>)
1155  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
1156 {
1157  static_assert(sizeof...(indices) == dim,
1158  "dim should match the number of indices");
1159 }
1160 
1161 
1162 
1163 template <int rank_, int dim, typename Number>
1166  // We would like to use =default, but this causes compile errors with some
1167  // MSVC versions and internal compiler errors with -O1 in gcc 5.4.
1168  : values{}
1169 {}
1170 
1171 
1172 
1173 template <int rank_, int dim, typename Number>
1176  : Tensor(initializer, std::make_index_sequence<dim>{})
1177 {}
1178 
1179 
1180 
1181 template <int rank_, int dim, typename Number>
1182 template <typename ElementType, typename MemorySpace>
1185  const ArrayView<ElementType, MemorySpace> &initializer)
1186 {
1188 
1189  for (unsigned int i = 0; i < n_independent_components; ++i)
1190  (*this)[unrolled_to_component_indices(i)] = initializer[i];
1191 }
1192 
1193 
1194 
1195 template <int rank_, int dim, typename Number>
1196 template <typename OtherNumber>
1199  const Tensor<rank_, dim, OtherNumber> &initializer)
1200  : Tensor(initializer, std::make_index_sequence<dim>{})
1201 {}
1202 
1203 
1204 template <int rank_, int dim, typename Number>
1205 template <typename OtherNumber>
1206 constexpr DEAL_II_ALWAYS_INLINE
1208  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1209  : Tensor(initializer, std::make_index_sequence<dim>{})
1210 {}
1211 
1212 
1213 template <int rank_, int dim, typename Number>
1214 template <typename OtherNumber>
1216  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1217 {
1218  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
1219 }
1220 
1221 
1222 
1223 namespace internal
1224 {
1225  namespace TensorSubscriptor
1226  {
1227  template <typename ArrayElementType, int dim>
1228  constexpr inline DEAL_II_ALWAYS_INLINE
1229  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1230  subscript(ArrayElementType * values,
1231  const unsigned int i,
1232  std::integral_constant<int, dim>)
1233  {
1234  // We cannot use Assert in a CUDA kernel
1235 # ifndef __CUDA_ARCH__
1236  AssertIndexRange(i, dim);
1237 # endif
1238  return values[i];
1239  }
1240 
1241  // The variables within this struct will be referenced in the next function.
1242  // It is a workaround that allows returning a reference to a static variable
1243  // while allowing constexpr evaluation of the function.
1244  // It has to be defined outside the function because constexpr functions
1245  // cannot define static variables
1246  template <typename ArrayElementType>
1247  struct Uninitialized
1248  {
1249  static ArrayElementType value;
1250  };
1251 
1252  template <typename Type>
1253  Type Uninitialized<Type>::value;
1254 
1255  template <typename ArrayElementType>
1256  constexpr inline DEAL_II_ALWAYS_INLINE
1257  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1258  subscript(ArrayElementType *,
1259  const unsigned int,
1260  std::integral_constant<int, 0>)
1261  {
1262  // We cannot use Assert in a CUDA kernel
1263 # ifndef __CUDA_ARCH__
1264  Assert(
1265  false,
1266  ExcMessage(
1267  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1268 # endif
1269  return Uninitialized<ArrayElementType>::value;
1270  }
1271  } // namespace TensorSubscriptor
1272 } // namespace internal
1273 
1274 
1275 template <int rank_, int dim, typename Number>
1278  operator[](const unsigned int i)
1279 {
1280  return ::internal::TensorSubscriptor::subscript(
1281  values, i, std::integral_constant<int, dim>());
1282 }
1283 
1284 
1285 template <int rank_, int dim, typename Number>
1286 constexpr DEAL_II_ALWAYS_INLINE
1288  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1289 {
1290 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1291  AssertIndexRange(i, dim);
1292 # endif
1293 
1294  return values[i];
1295 }
1296 
1297 
1298 template <int rank_, int dim, typename Number>
1299 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
1301  operator[](const TableIndices<rank_> &indices) const
1302 {
1303 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1304  Assert(dim != 0,
1305  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1306 # endif
1307 
1308  return TensorAccessors::extract<rank_>(*this, indices);
1309 }
1310 
1311 
1312 
1313 template <int rank_, int dim, typename Number>
1314 constexpr inline DEAL_II_ALWAYS_INLINE Number &Tensor<rank_, dim, Number>::
1315  operator[](const TableIndices<rank_> &indices)
1316 {
1317 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1318  Assert(dim != 0,
1319  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1320 # endif
1321 
1322  return TensorAccessors::extract<rank_>(*this, indices);
1323 }
1324 
1325 
1326 
1327 template <int rank_, int dim, typename Number>
1328 inline Number *
1330 {
1331  return std::addressof(
1332  this->operator[](this->unrolled_to_component_indices(0)));
1333 }
1334 
1335 
1336 
1337 template <int rank_, int dim, typename Number>
1338 inline const Number *
1340 {
1341  return std::addressof(
1342  this->operator[](this->unrolled_to_component_indices(0)));
1343 }
1344 
1345 
1346 
1347 template <int rank_, int dim, typename Number>
1348 inline Number *
1350 {
1352 }
1353 
1354 
1355 
1356 template <int rank_, int dim, typename Number>
1357 inline const Number *
1359 {
1361 }
1362 
1363 
1364 
1365 template <int rank_, int dim, typename Number>
1366 template <typename OtherNumber>
1369 {
1370  // The following loop could be written more concisely using std::copy, but
1371  // that function is only constexpr from C++20 on.
1372  for (unsigned int i = 0; i < dim; ++i)
1373  values[i] = t.values[i];
1374  return *this;
1375 }
1376 
1377 
1378 template <int rank_, int dim, typename Number>
1381 {
1383  (void)d;
1384 
1385  for (unsigned int i = 0; i < dim; ++i)
1387  return *this;
1388 }
1389 
1390 
1391 template <int rank_, int dim, typename Number>
1392 template <typename OtherNumber>
1393 constexpr inline bool
1396 {
1397  for (unsigned int i = 0; i < dim; ++i)
1398  if (values[i] != p.values[i])
1399  return false;
1400  return true;
1401 }
1402 
1403 
1404 // At some places in the library, we have Point<0> for formal reasons
1405 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1406 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1407 // in the above function that the loop end check always fails, we
1408 // implement this function here
1409 template <>
1410 template <>
1411 constexpr inline bool
1413 {
1414  return true;
1415 }
1416 
1417 
1418 template <int rank_, int dim, typename Number>
1419 template <typename OtherNumber>
1420 constexpr bool
1423 {
1424  return !((*this) == p);
1425 }
1426 
1427 
1428 template <int rank_, int dim, typename Number>
1429 template <typename OtherNumber>
1430 constexpr inline DEAL_II_ALWAYS_INLINE
1434 {
1435  for (unsigned int i = 0; i < dim; ++i)
1436  values[i] += p.values[i];
1437  return *this;
1438 }
1439 
1440 
1441 template <int rank_, int dim, typename Number>
1442 template <typename OtherNumber>
1443 constexpr inline DEAL_II_ALWAYS_INLINE
1447 {
1448  for (unsigned int i = 0; i < dim; ++i)
1449  values[i] -= p.values[i];
1450  return *this;
1451 }
1452 
1453 
1454 template <int rank_, int dim, typename Number>
1455 template <typename OtherNumber>
1456 constexpr inline DEAL_II_ALWAYS_INLINE
1458  Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1459 {
1460  for (unsigned int i = 0; i < dim; ++i)
1461  values[i] *= s;
1462  return *this;
1463 }
1464 
1465 
1466 namespace internal
1467 {
1468  namespace TensorImplementation
1469  {
1470  template <int rank,
1471  int dim,
1472  typename Number,
1473  typename OtherNumber,
1474  typename std::enable_if<
1475  !std::is_integral<
1476  typename ProductType<Number, OtherNumber>::type>::value &&
1477  !std::is_same<Number, Differentiation::SD::Expression>::value,
1478  int>::type = 0>
1479  constexpr DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
1481  const OtherNumber &factor)
1482  {
1483  const Number inverse_factor = Number(1.) / factor;
1484  // recurse over the base objects
1485  for (unsigned int d = 0; d < dim; ++d)
1486  t[d] *= inverse_factor;
1487  }
1488 
1489 
1490  template <int rank,
1491  int dim,
1492  typename Number,
1493  typename OtherNumber,
1494  typename std::enable_if<
1495  std::is_integral<
1496  typename ProductType<Number, OtherNumber>::type>::value ||
1497  std::is_same<Number, Differentiation::SD::Expression>::value,
1498  int>::type = 0>
1499  constexpr DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
1501  const OtherNumber &factor)
1502  {
1503  // recurse over the base objects
1504  for (unsigned int d = 0; d < dim; ++d)
1505  t[d] /= factor;
1506  }
1507  } // namespace TensorImplementation
1508 } // namespace internal
1509 
1510 
1511 template <int rank_, int dim, typename Number>
1512 template <typename OtherNumber>
1513 constexpr inline DEAL_II_ALWAYS_INLINE
1515  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1516 {
1518  return *this;
1519 }
1520 
1521 
1522 template <int rank_, int dim, typename Number>
1523 constexpr inline DEAL_II_ALWAYS_INLINE
1526 {
1528 
1529  for (unsigned int i = 0; i < dim; ++i)
1530  tmp.values[i] = -values[i];
1531 
1532  return tmp;
1533 }
1534 
1535 
1536 template <int rank_, int dim, typename Number>
1539 {
1540  return std::sqrt(norm_square());
1541 }
1542 
1543 
1544 template <int rank_, int dim, typename Number>
1548 {
1550  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1551  for (unsigned int i = 0; i < dim; ++i)
1552  s += values[i].norm_square();
1553 
1554  return s;
1555 }
1556 
1557 
1558 template <int rank_, int dim, typename Number>
1559 template <typename OtherNumber>
1560 inline void
1561 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
1562 {
1563  AssertDimension(result.size(),
1564  (Utilities::fixed_power<rank_, unsigned int>(dim)));
1565 
1566  unsigned int index = 0;
1567  unroll_recursion(result, index);
1568 }
1569 
1570 
1571 template <int rank_, int dim, typename Number>
1572 template <typename OtherNumber>
1573 inline void
1574 Tensor<rank_, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1575  unsigned int & index) const
1576 {
1577  for (unsigned int i = 0; i < dim; ++i)
1578  values[i].unroll_recursion(result, index);
1579 }
1580 
1581 
1582 template <int rank_, int dim, typename Number>
1583 constexpr inline unsigned int
1585  const TableIndices<rank_> &indices)
1586 {
1587  unsigned int index = 0;
1588  for (int r = 0; r < rank_; ++r)
1589  index = index * dim + indices[r];
1590 
1591  return index;
1592 }
1593 
1594 
1595 
1596 namespace internal
1597 {
1598  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1599  // and rank=2. Make sure we don't have compiler warnings.
1600 
1601  template <int dim>
1602  inline constexpr unsigned int
1603  mod(const unsigned int x)
1604  {
1605  return x % dim;
1606  }
1607 
1608  template <>
1609  inline unsigned int
1610  mod<0>(const unsigned int x)
1611  {
1612  Assert(false, ExcInternalError());
1613  return x;
1614  }
1615 
1616  template <int dim>
1617  inline constexpr unsigned int
1618  div(const unsigned int x)
1619  {
1620  return x / dim;
1621  }
1622 
1623  template <>
1624  inline unsigned int
1625  div<0>(const unsigned int x)
1626  {
1627  Assert(false, ExcInternalError());
1628  return x;
1629  }
1630 
1631 } // namespace internal
1632 
1633 
1634 
1635 template <int rank_, int dim, typename Number>
1636 constexpr inline TableIndices<rank_>
1638 {
1639  AssertIndexRange(i, n_independent_components);
1640 
1641  TableIndices<rank_> indices;
1642 
1643  unsigned int remainder = i;
1644  for (int r = rank_ - 1; r >= 0; --r)
1645  {
1646  indices[r] = internal::mod<dim>(remainder);
1647  remainder = internal::div<dim>(remainder);
1648  }
1649  Assert(remainder == 0, ExcInternalError());
1650 
1651  return indices;
1652 }
1653 
1654 
1655 template <int rank_, int dim, typename Number>
1656 constexpr inline void
1658 {
1659  for (unsigned int i = 0; i < dim; ++i)
1661 }
1662 
1663 
1664 template <int rank_, int dim, typename Number>
1665 constexpr std::size_t
1667 {
1668  return sizeof(Tensor<rank_, dim, Number>);
1669 }
1670 
1671 
1672 template <int rank_, int dim, typename Number>
1673 template <class Archive>
1674 inline void
1675 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1676 {
1677  ar &values;
1678 }
1679 
1680 
1681 template <int rank_, int dim, typename Number>
1683 
1684 #endif // DOXYGEN
1685 
1686 /* ----------------- Non-member functions operating on tensors. ------------ */
1687 
1692 
1700 template <int rank_, int dim, typename Number>
1701 inline std::ostream &
1702 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1703 {
1704  for (unsigned int i = 0; i < dim; ++i)
1705  {
1706  out << p[i];
1707  if (i != dim - 1)
1708  out << ' ';
1709  }
1710 
1711  return out;
1712 }
1713 
1714 
1721 template <int dim, typename Number>
1722 inline std::ostream &
1723 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1724 {
1725  out << static_cast<const Number &>(p);
1726  return out;
1727 }
1728 
1729 
1731 
1735 
1736 
1747 template <int dim, typename Number, typename Other>
1750  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1751 {
1752  return object * static_cast<const Number &>(t);
1753 }
1754 
1755 
1756 
1767 template <int dim, typename Number, typename Other>
1770  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1771 {
1772  return static_cast<const Number &>(t) * object;
1773 }
1774 
1775 
1787 template <int dim, typename Number, typename OtherNumber>
1791  const Tensor<0, dim, OtherNumber> &src2)
1792 {
1793  return static_cast<const Number &>(src1) *
1794  static_cast<const OtherNumber &>(src2);
1795 }
1796 
1797 
1805 template <int dim, typename Number, typename OtherNumber>
1807  Tensor<0,
1808  dim,
1809  typename ProductType<Number,
1810  typename EnableIfScalar<OtherNumber>::type>::type>
1811  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
1812 {
1813  return static_cast<const Number &>(t) / factor;
1814 }
1815 
1816 
1824 template <int dim, typename Number, typename OtherNumber>
1828  const Tensor<0, dim, OtherNumber> &q)
1829 {
1830  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1831 }
1832 
1833 
1841 template <int dim, typename Number, typename OtherNumber>
1845  const Tensor<0, dim, OtherNumber> &q)
1846 {
1847  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1848 }
1849 
1850 
1863 template <int rank, int dim, typename Number, typename OtherNumber>
1865  Tensor<rank,
1866  dim,
1867  typename ProductType<Number,
1868  typename EnableIfScalar<OtherNumber>::type>::type>
1869  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1870 {
1871  // recurse over the base objects
1873  for (unsigned int d = 0; d < dim; ++d)
1874  tt[d] = t[d] * factor;
1875  return tt;
1876 }
1877 
1878 
1891 template <int rank, int dim, typename Number, typename OtherNumber>
1893  Tensor<rank,
1894  dim,
1896  OtherNumber>::type>
1897  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
1898 {
1899  // simply forward to the operator above
1900  return t * factor;
1901 }
1902 
1903 
1904 namespace internal
1905 {
1906  namespace TensorImplementation
1907  {
1908  template <int rank,
1909  int dim,
1910  typename Number,
1911  typename OtherNumber,
1912  typename std::enable_if<
1913  !std::is_integral<
1914  typename ProductType<Number, OtherNumber>::type>::value,
1915  int>::type = 0>
1919  const OtherNumber & factor)
1920  {
1922  const Number inverse_factor = Number(1.) / factor;
1923  // recurse over the base objects
1924  for (unsigned int d = 0; d < dim; ++d)
1925  tt[d] = t[d] * inverse_factor;
1926  return tt;
1927  }
1928 
1929 
1930  template <int rank,
1931  int dim,
1932  typename Number,
1933  typename OtherNumber,
1934  typename std::enable_if<
1935  std::is_integral<
1936  typename ProductType<Number, OtherNumber>::type>::value,
1937  int>::type = 0>
1941  const OtherNumber & factor)
1942  {
1944  // recurse over the base objects
1945  for (unsigned int d = 0; d < dim; ++d)
1946  tt[d] = t[d] / factor;
1947  return tt;
1948  }
1949  } // namespace TensorImplementation
1950 } // namespace internal
1951 
1952 
1962 template <int rank, int dim, typename Number, typename OtherNumber>
1964  Tensor<rank,
1965  dim,
1966  typename ProductType<Number,
1967  typename EnableIfScalar<OtherNumber>::type>::type>
1968  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1969 {
1971 }
1972 
1973 
1983 template <int rank, int dim, typename Number, typename OtherNumber>
1988 {
1990 
1991  for (unsigned int i = 0; i < dim; ++i)
1992  tmp[i] += q[i];
1993 
1994  return tmp;
1995 }
1996 
1997 
2007 template <int rank, int dim, typename Number, typename OtherNumber>
2012 {
2014 
2015  for (unsigned int i = 0; i < dim; ++i)
2016  tmp[i] -= q[i];
2017 
2018  return tmp;
2019 }
2020 
2027 template <int dim, typename Number, typename OtherNumber>
2028 inline constexpr DEAL_II_ALWAYS_INLINE
2031  const Tensor<0, dim, OtherNumber> &src2)
2032 {
2034 
2035  tmp *= src2;
2036 
2037  return tmp;
2038 }
2039 
2056 template <int rank, int dim, typename Number, typename OtherNumber>
2057 inline constexpr DEAL_II_ALWAYS_INLINE
2060  const Tensor<rank, dim, OtherNumber> &src2)
2061 {
2063 
2064  for (unsigned int i = 0; i < dim; ++i)
2065  tmp[i] = schur_product(Tensor<rank - 1, dim, Number>(src1[i]),
2067 
2068  return tmp;
2069 }
2070 
2072 
2076 
2077 
2100 template <int rank_1,
2101  int rank_2,
2102  int dim,
2103  typename Number,
2104  typename OtherNumber,
2105  typename = typename std::enable_if<rank_1 >= 1 && rank_2 >= 1>::type>
2106 constexpr inline DEAL_II_ALWAYS_INLINE
2107  typename Tensor<rank_1 + rank_2 - 2,
2108  dim,
2112 {
2113  typename Tensor<rank_1 + rank_2 - 2,
2114  dim,
2116  result{};
2117 
2118  TensorAccessors::internal::
2119  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2120  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2121  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2122 
2123  return result;
2124 }
2125 
2126 
2155 template <int index_1,
2156  int index_2,
2157  int rank_1,
2158  int rank_2,
2159  int dim,
2160  typename Number,
2161  typename OtherNumber>
2162 constexpr inline DEAL_II_ALWAYS_INLINE
2163  typename Tensor<rank_1 + rank_2 - 2,
2164  dim,
2168 {
2169  Assert(0 <= index_1 && index_1 < rank_1,
2170  ExcMessage(
2171  "The specified index_1 must lie within the range [0,rank_1)"));
2172  Assert(0 <= index_2 && index_2 < rank_2,
2173  ExcMessage(
2174  "The specified index_2 must lie within the range [0,rank_2)"));
2175 
2176  using namespace TensorAccessors;
2177  using namespace TensorAccessors::internal;
2178 
2179  // Reorder index_1 to the end of src1:
2181  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2182 
2183  // Reorder index_2 to the end of src2:
2184  const ReorderedIndexView<index_2,
2185  rank_2,
2187  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2188 
2189  typename Tensor<rank_1 + rank_2 - 2,
2190  dim,
2192  result{};
2193  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2194  return result;
2195 }
2196 
2197 
2228 template <int index_1,
2229  int index_2,
2230  int index_3,
2231  int index_4,
2232  int rank_1,
2233  int rank_2,
2234  int dim,
2235  typename Number,
2236  typename OtherNumber>
2237 constexpr inline
2238  typename Tensor<rank_1 + rank_2 - 4,
2239  dim,
2243 {
2244  Assert(0 <= index_1 && index_1 < rank_1,
2245  ExcMessage(
2246  "The specified index_1 must lie within the range [0,rank_1)"));
2247  Assert(0 <= index_3 && index_3 < rank_1,
2248  ExcMessage(
2249  "The specified index_3 must lie within the range [0,rank_1)"));
2250  Assert(index_1 != index_3,
2251  ExcMessage("index_1 and index_3 must not be the same"));
2252  Assert(0 <= index_2 && index_2 < rank_2,
2253  ExcMessage(
2254  "The specified index_2 must lie within the range [0,rank_2)"));
2255  Assert(0 <= index_4 && index_4 < rank_2,
2256  ExcMessage(
2257  "The specified index_4 must lie within the range [0,rank_2)"));
2258  Assert(index_2 != index_4,
2259  ExcMessage("index_2 and index_4 must not be the same"));
2260 
2261  using namespace TensorAccessors;
2262  using namespace TensorAccessors::internal;
2263 
2264  // Reorder index_1 to the end of src1:
2266  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2267 
2268  // Reorder index_2 to the end of src2:
2270  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2271 
2272  // Now, reorder index_3 to the end of src1. We have to make sure to
2273  // preserve the original ordering: index_1 has been removed. If
2274  // index_3 > index_1, we have to use (index_3 - 1) instead:
2276  (index_3 < index_1 ? index_3 : index_3 - 1),
2277  rank_1,
2278  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2279  reord_3 =
2280  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2281  index_3 - 1,
2282  rank_1 > (reord_1);
2283 
2284  // Now, reorder index_4 to the end of src2. We have to make sure to
2285  // preserve the original ordering: index_2 has been removed. If
2286  // index_4 > index_2, we have to use (index_4 - 1) instead:
2287  ReorderedIndexView<
2288  (index_4 < index_2 ? index_4 : index_4 - 1),
2289  rank_2,
2290  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2291  reord_4 =
2292  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2293  index_4 - 1,
2294  rank_2 > (reord_2);
2295 
2296  typename Tensor<rank_1 + rank_2 - 4,
2297  dim,
2299  result{};
2300  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2301  return result;
2302 }
2303 
2304 
2317 template <int rank, int dim, typename Number, typename OtherNumber>
2318 constexpr inline DEAL_II_ALWAYS_INLINE
2321  const Tensor<rank, dim, OtherNumber> &right)
2322 {
2323  typename ProductType<Number, OtherNumber>::type result{};
2324  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2325  return result;
2326 }
2327 
2328 
2346 template <template <int, int, typename> class TensorT1,
2347  template <int, int, typename> class TensorT2,
2348  template <int, int, typename> class TensorT3,
2349  int rank_1,
2350  int rank_2,
2351  int dim,
2352  typename T1,
2353  typename T2,
2354  typename T3>
2355 constexpr inline DEAL_II_ALWAYS_INLINE
2357  contract3(const TensorT1<rank_1, dim, T1> & left,
2358  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2359  const TensorT3<rank_2, dim, T3> & right)
2360 {
2361  using return_type =
2363  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2364  middle,
2365  right);
2366 }
2367 
2368 
2379 template <int rank_1,
2380  int rank_2,
2381  int dim,
2382  typename Number,
2383  typename OtherNumber>
2384 constexpr inline DEAL_II_ALWAYS_INLINE
2388 {
2389  typename Tensor<rank_1 + rank_2,
2390  dim,
2392  result{};
2393  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2394  return result;
2395 }
2396 
2397 
2399 
2403 
2404 
2415 template <int dim, typename Number>
2418 {
2419  Assert(dim == 2, ExcInternalError());
2420 
2421  Tensor<1, dim, Number> result;
2422 
2423  result[0] = src[1];
2424  result[1] = -src[0];
2425 
2426  return result;
2427 }
2428 
2429 
2439 template <int dim, typename Number1, typename Number2>
2440 constexpr inline DEAL_II_ALWAYS_INLINE
2443  const Tensor<1, dim, Number2> &src2)
2444 {
2445  Assert(dim == 3, ExcInternalError());
2446 
2448 
2449  // avoid compiler warnings
2450  constexpr int s0 = 0 % dim;
2451  constexpr int s1 = 1 % dim;
2452  constexpr int s2 = 2 % dim;
2453 
2454  result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2455  result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2456  result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2457 
2458  return result;
2459 }
2460 
2461 
2463 
2467 
2468 
2474 template <int dim, typename Number>
2475 constexpr inline DEAL_II_ALWAYS_INLINE Number
2477 {
2478  // Compute the determinant using the Laplace expansion of the
2479  // determinant. We expand along the last row.
2480  Number det = internal::NumberType<Number>::value(0.0);
2481 
2482  for (unsigned int k = 0; k < dim; ++k)
2483  {
2484  Tensor<2, dim - 1, Number> minor;
2485  for (unsigned int i = 0; i < dim - 1; ++i)
2486  for (unsigned int j = 0; j < dim - 1; ++j)
2487  minor[i][j] = t[i][j < k ? j : j + 1];
2488 
2489  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2490 
2491  det += t[dim - 1][k] * cofactor;
2492  }
2493 
2494  return ((dim % 2 == 0) ? 1. : -1.) * det;
2495 }
2496 
2502 template <typename Number>
2503 constexpr DEAL_II_ALWAYS_INLINE Number
2505 {
2506  return t[0][0];
2507 }
2508 
2514 template <typename Number>
2515 constexpr DEAL_II_ALWAYS_INLINE Number
2517 {
2518  // hard-coded for efficiency reasons
2519  return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2520 }
2521 
2527 template <typename Number>
2528 constexpr DEAL_II_ALWAYS_INLINE Number
2530 {
2531  // hard-coded for efficiency reasons
2532  const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2533  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2534  const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2535  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2536  const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2537  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2538  return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2539 }
2540 
2541 
2548 template <int dim, typename Number>
2549 constexpr inline DEAL_II_ALWAYS_INLINE Number
2551 {
2552  Number t = d[0][0];
2553  for (unsigned int i = 1; i < dim; ++i)
2554  t += d[i][i];
2555  return t;
2556 }
2557 
2558 
2567 template <int dim, typename Number>
2568 constexpr inline Tensor<2, dim, Number>
2570 {
2571  Number return_tensor[dim][dim];
2572 
2573  // if desired, take over the
2574  // inversion of a 4x4 tensor
2575  // from the FullMatrix
2576  AssertThrow(false, ExcNotImplemented());
2577 
2578  return Tensor<2, dim, Number>(return_tensor);
2579 }
2580 
2581 
2582 #ifndef DOXYGEN
2583 
2584 template <typename Number>
2586  invert(const Tensor<2, 1, Number> &t)
2587 {
2588  Tensor<2, 1, Number> return_tensor;
2589 
2590  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2591 
2592  return return_tensor;
2593 }
2594 
2595 
2596 template <typename Number>
2598  invert(const Tensor<2, 2, Number> &t)
2599 {
2600  Tensor<2, 2, Number> return_tensor;
2601 
2602  const Number inv_det_t = internal::NumberType<Number>::value(
2603  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2604  return_tensor[0][0] = t[1][1];
2605  return_tensor[0][1] = -t[0][1];
2606  return_tensor[1][0] = -t[1][0];
2607  return_tensor[1][1] = t[0][0];
2608  return_tensor *= inv_det_t;
2609 
2610  return return_tensor;
2611 }
2612 
2613 
2614 template <typename Number>
2616  invert(const Tensor<2, 3, Number> &t)
2617 {
2618  Tensor<2, 3, Number> return_tensor;
2619 
2620  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2621  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2622  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2623  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2624  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2625  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2626  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2627  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2628  return_tensor[1][1] = internal::NumberType<Number>::value(t[0][0] * t[2][2]) -
2629  internal::NumberType<Number>::value(t[0][2] * t[2][0]);
2630  return_tensor[1][2] = internal::NumberType<Number>::value(t[0][2] * t[1][0]) -
2631  internal::NumberType<Number>::value(t[0][0] * t[1][2]);
2632  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2633  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2634  return_tensor[2][1] = internal::NumberType<Number>::value(t[0][1] * t[2][0]) -
2635  internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2636  return_tensor[2][2] = internal::NumberType<Number>::value(t[0][0] * t[1][1]) -
2637  internal::NumberType<Number>::value(t[0][1] * t[1][0]);
2638  const Number inv_det_t = internal::NumberType<Number>::value(
2639  1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2640  t[0][2] * return_tensor[2][0]));
2641  return_tensor *= inv_det_t;
2642 
2643  return return_tensor;
2644 }
2645 
2646 #endif /* DOXYGEN */
2647 
2648 
2654 template <int dim, typename Number>
2657 {
2659  for (unsigned int i = 0; i < dim; ++i)
2660  {
2661  tt[i][i] = t[i][i];
2662  for (unsigned int j = i + 1; j < dim; ++j)
2663  {
2664  tt[i][j] = t[j][i];
2665  tt[j][i] = t[i][j];
2666  };
2667  }
2668  return tt;
2669 }
2670 
2671 
2685 template <int dim, typename Number>
2686 constexpr Tensor<2, dim, Number>
2688 {
2689  return determinant(t) * invert(t);
2690 }
2691 
2692 
2706 template <int dim, typename Number>
2707 constexpr Tensor<2, dim, Number>
2709 {
2710  return transpose(adjugate(t));
2711 }
2712 
2713 
2777 template <int dim, typename Number>
2780 
2781 
2789 template <int dim, typename Number>
2790 inline Number
2792 {
2794  for (unsigned int j = 0; j < dim; ++j)
2795  {
2797  for (unsigned int i = 0; i < dim; ++i)
2798  sum += std::fabs(t[i][j]);
2799 
2800  if (sum > max)
2801  max = sum;
2802  }
2803 
2804  return max;
2805 }
2806 
2807 
2815 template <int dim, typename Number>
2816 inline Number
2818 {
2820  for (unsigned int i = 0; i < dim; ++i)
2821  {
2823  for (unsigned int j = 0; j < dim; ++j)
2824  sum += std::fabs(t[i][j]);
2825 
2826  if (sum > max)
2827  max = sum;
2828  }
2829 
2830  return max;
2831 }
2832 
2834 
2835 
2836 #ifndef DOXYGEN
2837 
2838 
2839 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2840 
2841 // Specialization of functions for ADOL-C number types when
2842 // the advanced branching feature is used
2843 template <int dim>
2844 inline adouble
2846 {
2848  for (unsigned int j = 0; j < dim; ++j)
2849  {
2851  for (unsigned int i = 0; i < dim; ++i)
2852  sum += std::fabs(t[i][j]);
2853 
2854  condassign(max, (sum > max), sum, max);
2855  }
2856 
2857  return max;
2858 }
2859 
2860 
2861 template <int dim>
2862 inline adouble
2864 {
2866  for (unsigned int i = 0; i < dim; ++i)
2867  {
2869  for (unsigned int j = 0; j < dim; ++j)
2870  sum += std::fabs(t[i][j]);
2871 
2872  condassign(max, (sum > max), sum, max);
2873  }
2874 
2875  return max;
2876 }
2877 
2878 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2879 
2880 
2881 #endif // DOXYGEN
2882 
2884 
2885 #endif
constexpr Tensor & operator+=(const Tensor< rank_, dim, OtherNumber > &)
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
constexpr Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:2569
Number * begin_raw()
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
constexpr ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:2320
constexpr Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2708
static constexpr const T & value(const T &t)
Definition: numbers.h:693
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:915
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
constexpr Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2656
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
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:577
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1750
constexpr Tensor operator-() const
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2417
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2791
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
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:2442
constexpr void clear()
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2817
constexpr Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:766
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
STL namespace.
typename Tensor< rank_ - 1, dim, Number >::array_type[(dim !=0) ? dim :1] array_type
Definition: tensor.h:490
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
static real_type abs(const number &x)
Definition: numbers.h:599
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: point.h:110
constexpr Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
std::size_t size() const
Definition: array_view.h:574
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:760
Number * end_raw()
static ::ExceptionBase & ExcMessage(std::string arg1)
constexpr value_type & operator[](const unsigned int i)
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
static const char T
#define Assert(cond, exc)
Definition: exceptions.h:1465
constexpr Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2550
constexpr internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
void serialize(Archive &ar, const unsigned int version)
constexpr bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:96
typename Tensor< rank_ - 1, dim, Number >::tensor_type value_type
Definition: tensor.h:483
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:2241
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
Expression fabs(const Expression &x)
static constexpr unsigned int rank
Definition: tensor.h:469
constexpr Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2687
constexpr Tensor & operator*=(const OtherNumber &factor)
constexpr Tensor & operator/=(const OtherNumber &factor)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static constexpr std::size_t memory_consumption()
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:2357
static const char A
constexpr bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
typename numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:126
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:1811
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > division_operator(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1918
Definition: tensor.h:449
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:1827
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:931
#define DEAL_II_CUDA_HOST_DEV
Definition: numbers.h:34
Expression operator-(Expression lhs, const Expression &rhs)
static ::ExceptionBase & ExcNotImplemented()
void unroll(Vector< OtherNumber > &result) const
Tensor< 2, dim, Number > project_onto_orthogonal_tensors(const Tensor< 2, dim, Number > &A)
numbers::NumberTraits< Number >::real_type norm() const
static constexpr unsigned int n_independent_components
Definition: tensor.h:475
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:2030
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:2386
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:2166
Expression operator!=(const Expression &lhs, const Expression &rhs)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
inline ::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
static ::ExceptionBase & ExcInternalError()
constexpr Tensor()
constexpr Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2476