Reference documentation for deal.II version Git d2d482dcec 2020-10-26 10:29:03 -0400
\(\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 - 2020 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  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV operator Number &();
204 
213  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV operator const Number &() const;
214 
222  template <typename OtherNumber>
224  operator=(const Tensor<0, dim, OtherNumber> &rhs);
225 
226 #ifdef __INTEL_COMPILER
227 
236  operator=(const Tensor<0, dim, Number> &rhs);
237 #endif
238 
245  template <typename OtherNumber>
247  operator=(const OtherNumber &d);
248 
252  template <typename OtherNumber>
253  DEAL_II_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>
270  operator+=(const Tensor<0, dim, OtherNumber> &rhs);
271 
277  template <typename OtherNumber>
279  operator-=(const Tensor<0, dim, OtherNumber> &rhs);
280 
286  template <typename OtherNumber>
288  operator*=(const OtherNumber &factor);
289 
295  template <typename OtherNumber>
297  operator/=(const OtherNumber &factor);
298 
304  constexpr DEAL_II_CUDA_HOST_DEV Tensor
305  operator-() const;
306 
319  DEAL_II_CONSTEXPR void
320  clear();
321 
327  real_type
328  norm() const;
329 
337  norm_square() const;
338 
343  template <class Archive>
344  void
345  serialize(Archive &ar, const unsigned int version);
346 
351  using tensor_type = Number;
352 
353 private:
357  Number value;
358 
362  template <typename OtherNumber>
363  void
364  unroll_recursion(Vector<OtherNumber> &result,
365  unsigned int & start_index) const;
366 
367  // Allow an arbitrary Tensor to access the underlying values.
368  template <int, int, typename>
369  friend class Tensor;
370 };
371 
372 
373 
447 template <int rank_, int dim, typename Number>
448 class Tensor
449 {
450 public:
451  static_assert(rank_ >= 0,
452  "Tensors must have a rank greater than or equal to one.");
453  static_assert(dim >= 0,
454  "Tensors must have a dimension greater than or equal to one.");
463  static constexpr unsigned int dimension = dim;
464 
468  static constexpr unsigned int rank = rank_;
469 
474  static constexpr unsigned int n_independent_components =
475  Tensor<rank_ - 1, dim>::n_independent_components * dim;
476 
482  using value_type = typename Tensor<rank_ - 1, dim, Number>::tensor_type;
483 
488  using array_type =
489  typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1];
490 
497  Tensor();
498 
504  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
505  const array_type &initializer);
506 
520  template <typename ElementType, typename MemorySpace>
521  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
522  const ArrayView<ElementType, MemorySpace> &initializer);
523 
531  template <typename OtherNumber>
532  constexpr DEAL_II_CUDA_HOST_DEV
533  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
534 
538  template <typename OtherNumber>
539  constexpr Tensor(
540  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
541 
545  template <typename OtherNumber>
546  constexpr
547  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
548 
555  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  DEAL_II_CONSTEXPR const Number &
569  operator[](const TableIndices<rank_> &indices) const;
570 
574  DEAL_II_CONSTEXPR Number &operator[](const TableIndices<rank_> &indices);
575 
579  Number *
580  begin_raw();
581 
585  const Number *
586  begin_raw() const;
587 
591  Number *
592  end_raw();
593 
597  const Number *
598  end_raw() const;
599 
607  template <typename OtherNumber>
609  operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
610 
617  DEAL_II_CONSTEXPR Tensor &
618  operator=(const Number &d);
619 
623  template <typename OtherNumber>
624  DEAL_II_CONSTEXPR bool
626 
630  template <typename OtherNumber>
631  constexpr bool
633 
639  template <typename OtherNumber>
641  operator+=(const Tensor<rank_, dim, OtherNumber> &);
642 
648  template <typename OtherNumber>
650  operator-=(const Tensor<rank_, dim, OtherNumber> &);
651 
658  template <typename OtherNumber>
660  operator*=(const OtherNumber &factor);
661 
667  template <typename OtherNumber>
669  operator/=(const OtherNumber &factor);
670 
677  operator-() const;
678 
691  DEAL_II_CONSTEXPR void
692  clear();
693 
703  norm() const;
704 
713  norm_square() const;
714 
722  template <typename OtherNumber>
723  void
724  unroll(Vector<OtherNumber> &result) const;
725 
730  static DEAL_II_CONSTEXPR unsigned int
731  component_to_unrolled_index(const TableIndices<rank_> &indices);
732 
739  unrolled_to_component_indices(const unsigned int i);
740 
745  static constexpr std::size_t
747 
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>
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>
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>
956 {
958  return *this;
959 }
960 
961 
962 # ifdef __INTEL_COMPILER
963 template <int dim, typename Number>
967 {
968  value = p.value;
969  return *this;
970 }
971 # endif
972 
973 
974 template <int dim, typename Number>
975 template <typename OtherNumber>
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 DEAL_II_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>
1016 {
1017  value += p.value;
1018  return *this;
1019 }
1020 
1021 
1022 template <int dim, typename Number>
1023 template <typename OtherNumber>
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>
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>
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>
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 DEAL_II_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>
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>
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 DEAL_II_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>
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  ExcMessage("Only assignment with zero is allowed"));
1384  (void)d;
1385 
1386  for (unsigned int i = 0; i < dim; ++i)
1388  return *this;
1389 }
1390 
1391 
1392 template <int rank_, int dim, typename Number>
1393 template <typename OtherNumber>
1394 DEAL_II_CONSTEXPR inline bool
1397 {
1398  for (unsigned int i = 0; i < dim; ++i)
1399  if (values[i] != p.values[i])
1400  return false;
1401  return true;
1402 }
1403 
1404 
1405 // At some places in the library, we have Point<0> for formal reasons
1406 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1407 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1408 // in the above function that the loop end check always fails, we
1409 // implement this function here
1410 template <>
1411 template <>
1412 DEAL_II_CONSTEXPR inline bool
1414 {
1415  return true;
1416 }
1417 
1418 
1419 template <int rank_, int dim, typename Number>
1420 template <typename OtherNumber>
1421 constexpr bool
1424 {
1425  return !((*this) == p);
1426 }
1427 
1428 
1429 template <int rank_, int dim, typename Number>
1430 template <typename OtherNumber>
1435 {
1436  for (unsigned int i = 0; i < dim; ++i)
1437  values[i] += p.values[i];
1438  return *this;
1439 }
1440 
1441 
1442 template <int rank_, int dim, typename Number>
1443 template <typename OtherNumber>
1448 {
1449  for (unsigned int i = 0; i < dim; ++i)
1450  values[i] -= p.values[i];
1451  return *this;
1452 }
1453 
1454 
1455 template <int rank_, int dim, typename Number>
1456 template <typename OtherNumber>
1459  Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1460 {
1461  for (unsigned int i = 0; i < dim; ++i)
1462  values[i] *= s;
1463  return *this;
1464 }
1465 
1466 
1467 namespace internal
1468 {
1469  namespace TensorImplementation
1470  {
1471  template <int rank,
1472  int dim,
1473  typename Number,
1474  typename OtherNumber,
1475  typename std::enable_if<
1476  !std::is_integral<
1477  typename ProductType<Number, OtherNumber>::type>::value &&
1478  !std::is_same<Number, Differentiation::SD::Expression>::value,
1479  int>::type = 0>
1482  const OtherNumber &factor)
1483  {
1484  const Number inverse_factor = Number(1.) / factor;
1485  // recurse over the base objects
1486  for (unsigned int d = 0; d < dim; ++d)
1487  t[d] *= inverse_factor;
1488  }
1489 
1490 
1491  template <int rank,
1492  int dim,
1493  typename Number,
1494  typename OtherNumber,
1495  typename std::enable_if<
1496  std::is_integral<
1497  typename ProductType<Number, OtherNumber>::type>::value ||
1498  std::is_same<Number, Differentiation::SD::Expression>::value,
1499  int>::type = 0>
1502  const OtherNumber &factor)
1503  {
1504  // recurse over the base objects
1505  for (unsigned int d = 0; d < dim; ++d)
1506  t[d] /= factor;
1507  }
1508  } // namespace TensorImplementation
1509 } // namespace internal
1510 
1511 
1512 template <int rank_, int dim, typename Number>
1513 template <typename OtherNumber>
1516  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1517 {
1519  return *this;
1520 }
1521 
1522 
1523 template <int rank_, int dim, typename Number>
1527 {
1529 
1530  for (unsigned int i = 0; i < dim; ++i)
1531  tmp.values[i] = -values[i];
1532 
1533  return tmp;
1534 }
1535 
1536 
1537 template <int rank_, int dim, typename Number>
1540 {
1541  return std::sqrt(norm_square());
1542 }
1543 
1544 
1545 template <int rank_, int dim, typename Number>
1549 {
1551  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1552  for (unsigned int i = 0; i < dim; ++i)
1553  s += values[i].norm_square();
1554 
1555  return s;
1556 }
1557 
1558 
1559 template <int rank_, int dim, typename Number>
1560 template <typename OtherNumber>
1561 inline void
1562 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
1563 {
1564  AssertDimension(result.size(),
1565  (Utilities::fixed_power<rank_, unsigned int>(dim)));
1566 
1567  unsigned int index = 0;
1568  unroll_recursion(result, index);
1569 }
1570 
1571 
1572 template <int rank_, int dim, typename Number>
1573 template <typename OtherNumber>
1574 inline void
1575 Tensor<rank_, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1576  unsigned int & index) const
1577 {
1578  for (unsigned int i = 0; i < dim; ++i)
1579  values[i].unroll_recursion(result, index);
1580 }
1581 
1582 
1583 template <int rank_, int dim, typename Number>
1584 DEAL_II_CONSTEXPR inline unsigned int
1586  const TableIndices<rank_> &indices)
1587 {
1588  unsigned int index = 0;
1589  for (int r = 0; r < rank_; ++r)
1590  index = index * dim + indices[r];
1591 
1592  return index;
1593 }
1594 
1595 
1596 
1597 namespace internal
1598 {
1599  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1600  // and rank=2. Make sure we don't have compiler warnings.
1601 
1602  template <int dim>
1603  inline DEAL_II_CONSTEXPR unsigned int
1604  mod(const unsigned int x)
1605  {
1606  return x % dim;
1607  }
1608 
1609  template <>
1610  inline unsigned int
1611  mod<0>(const unsigned int x)
1612  {
1613  Assert(false, ExcInternalError());
1614  return x;
1615  }
1616 
1617  template <int dim>
1618  inline DEAL_II_CONSTEXPR unsigned int
1619  div(const unsigned int x)
1620  {
1621  return x / dim;
1622  }
1623 
1624  template <>
1625  inline unsigned int
1626  div<0>(const unsigned int x)
1627  {
1628  Assert(false, ExcInternalError());
1629  return x;
1630  }
1631 
1632 } // namespace internal
1633 
1634 
1635 
1636 template <int rank_, int dim, typename Number>
1639 {
1640  AssertIndexRange(i, n_independent_components);
1641 
1642  TableIndices<rank_> indices;
1643 
1644  unsigned int remainder = i;
1645  for (int r = rank_ - 1; r >= 0; --r)
1646  {
1647  indices[r] = internal::mod<dim>(remainder);
1648  remainder = internal::div<dim>(remainder);
1649  }
1650  Assert(remainder == 0, ExcInternalError());
1651 
1652  return indices;
1653 }
1654 
1655 
1656 template <int rank_, int dim, typename Number>
1657 DEAL_II_CONSTEXPR inline void
1659 {
1660  for (unsigned int i = 0; i < dim; ++i)
1662 }
1663 
1664 
1665 template <int rank_, int dim, typename Number>
1666 constexpr std::size_t
1668 {
1669  return sizeof(Tensor<rank_, dim, Number>);
1670 }
1671 
1672 
1673 template <int rank_, int dim, typename Number>
1674 template <class Archive>
1675 inline void
1676 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1677 {
1678  ar &values;
1679 }
1680 
1681 
1682 template <int rank_, int dim, typename Number>
1684 
1685 #endif // DOXYGEN
1686 
1687 /* ----------------- Non-member functions operating on tensors. ------------ */
1688 
1693 
1701 template <int rank_, int dim, typename Number>
1702 inline std::ostream &
1703 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1704 {
1705  for (unsigned int i = 0; i < dim; ++i)
1706  {
1707  out << p[i];
1708  if (i != dim - 1)
1709  out << ' ';
1710  }
1711 
1712  return out;
1713 }
1714 
1715 
1722 template <int dim, typename Number>
1723 inline std::ostream &
1724 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1725 {
1726  out << static_cast<const Number &>(p);
1727  return out;
1728 }
1729 
1730 
1732 
1736 
1737 
1748 template <int dim, typename Number, typename Other>
1751  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1752 {
1753  return object * static_cast<const Number &>(t);
1754 }
1755 
1756 
1757 
1768 template <int dim, typename Number, typename Other>
1771  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1772 {
1773  return static_cast<const Number &>(t) * object;
1774 }
1775 
1776 
1788 template <int dim, typename Number, typename OtherNumber>
1792  const Tensor<0, dim, OtherNumber> &src2)
1793 {
1794  return static_cast<const Number &>(src1) *
1795  static_cast<const OtherNumber &>(src2);
1796 }
1797 
1798 
1806 template <int dim, typename Number, typename OtherNumber>
1808  Tensor<0,
1809  dim,
1810  typename ProductType<Number,
1811  typename EnableIfScalar<OtherNumber>::type>::type>
1812  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
1813 {
1814  return static_cast<const Number &>(t) / factor;
1815 }
1816 
1817 
1825 template <int dim, typename Number, typename OtherNumber>
1829  const Tensor<0, dim, OtherNumber> &q)
1830 {
1831  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1832 }
1833 
1834 
1842 template <int dim, typename Number, typename OtherNumber>
1846  const Tensor<0, dim, OtherNumber> &q)
1847 {
1848  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1849 }
1850 
1851 
1864 template <int rank, int dim, typename Number, typename OtherNumber>
1866  Tensor<rank,
1867  dim,
1868  typename ProductType<Number,
1869  typename EnableIfScalar<OtherNumber>::type>::type>
1870  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1871 {
1872  // recurse over the base objects
1874  for (unsigned int d = 0; d < dim; ++d)
1875  tt[d] = t[d] * factor;
1876  return tt;
1877 }
1878 
1879 
1892 template <int rank, int dim, typename Number, typename OtherNumber>
1894  Tensor<rank,
1895  dim,
1897  OtherNumber>::type>
1898  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
1899 {
1900  // simply forward to the operator above
1901  return t * factor;
1902 }
1903 
1904 
1905 namespace internal
1906 {
1907  namespace TensorImplementation
1908  {
1909  template <int rank,
1910  int dim,
1911  typename Number,
1912  typename OtherNumber,
1913  typename std::enable_if<
1914  !std::is_integral<
1915  typename ProductType<Number, OtherNumber>::type>::value,
1916  int>::type = 0>
1920  const OtherNumber & factor)
1921  {
1923  const Number inverse_factor = Number(1.) / factor;
1924  // recurse over the base objects
1925  for (unsigned int d = 0; d < dim; ++d)
1926  tt[d] = t[d] * inverse_factor;
1927  return tt;
1928  }
1929 
1930 
1931  template <int rank,
1932  int dim,
1933  typename Number,
1934  typename OtherNumber,
1935  typename std::enable_if<
1936  std::is_integral<
1937  typename ProductType<Number, OtherNumber>::type>::value,
1938  int>::type = 0>
1942  const OtherNumber & factor)
1943  {
1945  // recurse over the base objects
1946  for (unsigned int d = 0; d < dim; ++d)
1947  tt[d] = t[d] / factor;
1948  return tt;
1949  }
1950  } // namespace TensorImplementation
1951 } // namespace internal
1952 
1953 
1963 template <int rank, int dim, typename Number, typename OtherNumber>
1965  Tensor<rank,
1966  dim,
1967  typename ProductType<Number,
1968  typename EnableIfScalar<OtherNumber>::type>::type>
1969  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1970 {
1972 }
1973 
1974 
1984 template <int rank, int dim, typename Number, typename OtherNumber>
1989 {
1991 
1992  for (unsigned int i = 0; i < dim; ++i)
1993  tmp[i] += q[i];
1994 
1995  return tmp;
1996 }
1997 
1998 
2008 template <int rank, int dim, typename Number, typename OtherNumber>
2013 {
2015 
2016  for (unsigned int i = 0; i < dim; ++i)
2017  tmp[i] -= q[i];
2018 
2019  return tmp;
2020 }
2021 
2028 template <int dim, typename Number, typename OtherNumber>
2032  const Tensor<0, dim, OtherNumber> &src2)
2033 {
2035 
2036  tmp *= src2;
2037 
2038  return tmp;
2039 }
2040 
2057 template <int rank, int dim, typename Number, typename OtherNumber>
2061  const Tensor<rank, dim, OtherNumber> &src2)
2062 {
2064 
2065  for (unsigned int i = 0; i < dim; ++i)
2066  tmp[i] = schur_product(Tensor<rank - 1, dim, Number>(src1[i]),
2068 
2069  return tmp;
2070 }
2071 
2073 
2077 
2078 
2101 template <int rank_1,
2102  int rank_2,
2103  int dim,
2104  typename Number,
2105  typename OtherNumber,
2106  typename = typename std::enable_if<rank_1 >= 1 && rank_2 >= 1>::type>
2108  typename Tensor<rank_1 + rank_2 - 2,
2109  dim,
2113 {
2114  typename Tensor<rank_1 + rank_2 - 2,
2115  dim,
2117  result{};
2118 
2119  TensorAccessors::internal::
2120  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2121  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2122  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2123 
2124  return result;
2125 }
2126 
2127 
2156 template <int index_1,
2157  int index_2,
2158  int rank_1,
2159  int rank_2,
2160  int dim,
2161  typename Number,
2162  typename OtherNumber>
2164  typename Tensor<rank_1 + rank_2 - 2,
2165  dim,
2169 {
2170  Assert(0 <= index_1 && index_1 < rank_1,
2171  ExcMessage(
2172  "The specified index_1 must lie within the range [0,rank_1)"));
2173  Assert(0 <= index_2 && index_2 < rank_2,
2174  ExcMessage(
2175  "The specified index_2 must lie within the range [0,rank_2)"));
2176 
2177  using namespace TensorAccessors;
2178  using namespace TensorAccessors::internal;
2179 
2180  // Reorder index_1 to the end of src1:
2182  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2183 
2184  // Reorder index_2 to the end of src2:
2186  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2187 
2188  typename Tensor<rank_1 + rank_2 - 2,
2189  dim,
2191  result{};
2192  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2193  return result;
2194 }
2195 
2196 
2227 template <int index_1,
2228  int index_2,
2229  int index_3,
2230  int index_4,
2231  int rank_1,
2232  int rank_2,
2233  int dim,
2234  typename Number,
2235  typename OtherNumber>
2236 DEAL_II_CONSTEXPR inline
2237  typename Tensor<rank_1 + rank_2 - 4,
2238  dim,
2242 {
2243  Assert(0 <= index_1 && index_1 < rank_1,
2244  ExcMessage(
2245  "The specified index_1 must lie within the range [0,rank_1)"));
2246  Assert(0 <= index_3 && index_3 < rank_1,
2247  ExcMessage(
2248  "The specified index_3 must lie within the range [0,rank_1)"));
2249  Assert(index_1 != index_3,
2250  ExcMessage("index_1 and index_3 must not be the same"));
2251  Assert(0 <= index_2 && index_2 < rank_2,
2252  ExcMessage(
2253  "The specified index_2 must lie within the range [0,rank_2)"));
2254  Assert(0 <= index_4 && index_4 < rank_2,
2255  ExcMessage(
2256  "The specified index_4 must lie within the range [0,rank_2)"));
2257  Assert(index_2 != index_4,
2258  ExcMessage("index_2 and index_4 must not be the same"));
2259 
2260  using namespace TensorAccessors;
2261  using namespace TensorAccessors::internal;
2262 
2263  // Reorder index_1 to the end of src1:
2265  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2266 
2267  // Reorder index_2 to the end of src2:
2269  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2270 
2271  // Now, reorder index_3 to the end of src1. We have to make sure to
2272  // preserve the original ordering: index_1 has been removed. If
2273  // index_3 > index_1, we have to use (index_3 - 1) instead:
2275  (index_3 < index_1 ? index_3 : index_3 - 1),
2276  rank_1,
2277  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2278  reord_3 =
2279  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2280  index_3 - 1,
2281  rank_1 > (reord_1);
2282 
2283  // Now, reorder index_4 to the end of src2. We have to make sure to
2284  // preserve the original ordering: index_2 has been removed. If
2285  // index_4 > index_2, we have to use (index_4 - 1) instead:
2286  ReorderedIndexView<
2287  (index_4 < index_2 ? index_4 : index_4 - 1),
2288  rank_2,
2289  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2290  reord_4 =
2291  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2292  index_4 - 1,
2293  rank_2 > (reord_2);
2294 
2295  typename Tensor<rank_1 + rank_2 - 4,
2296  dim,
2298  result{};
2299  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2300  return result;
2301 }
2302 
2303 
2316 template <int rank, int dim, typename Number, typename OtherNumber>
2320  const Tensor<rank, dim, OtherNumber> &right)
2321 {
2322  typename ProductType<Number, OtherNumber>::type result{};
2323  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2324  return result;
2325 }
2326 
2327 
2345 template <template <int, int, typename> class TensorT1,
2346  template <int, int, typename> class TensorT2,
2347  template <int, int, typename> class TensorT3,
2348  int rank_1,
2349  int rank_2,
2350  int dim,
2351  typename T1,
2352  typename T2,
2353  typename T3>
2356  contract3(const TensorT1<rank_1, dim, T1> & left,
2357  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2358  const TensorT3<rank_2, dim, T3> & right)
2359 {
2360  using return_type =
2362  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2363  middle,
2364  right);
2365 }
2366 
2367 
2378 template <int rank_1,
2379  int rank_2,
2380  int dim,
2381  typename Number,
2382  typename OtherNumber>
2387 {
2388  typename Tensor<rank_1 + rank_2,
2389  dim,
2391  result{};
2392  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2393  return result;
2394 }
2395 
2396 
2398 
2402 
2403 
2414 template <int dim, typename Number>
2417 {
2418  Assert(dim == 2, ExcInternalError());
2419 
2420  Tensor<1, dim, Number> result;
2421 
2422  result[0] = src[1];
2423  result[1] = -src[0];
2424 
2425  return result;
2426 }
2427 
2428 
2438 template <int dim, typename Number1, typename Number2>
2442  const Tensor<1, dim, Number2> &src2)
2443 {
2444  Assert(dim == 3, ExcInternalError());
2445 
2447 
2448  // avoid compiler warnings
2449  constexpr int s0 = 0 % dim;
2450  constexpr int s1 = 1 % dim;
2451  constexpr int s2 = 2 % dim;
2452 
2453  result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2454  result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2455  result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2456 
2457  return result;
2458 }
2459 
2460 
2462 
2466 
2467 
2473 template <int dim, typename Number>
2476 {
2477  // Compute the determinant using the Laplace expansion of the
2478  // determinant. We expand along the last row.
2479  Number det = internal::NumberType<Number>::value(0.0);
2480 
2481  for (unsigned int k = 0; k < dim; ++k)
2482  {
2483  Tensor<2, dim - 1, Number> minor;
2484  for (unsigned int i = 0; i < dim - 1; ++i)
2485  for (unsigned int j = 0; j < dim - 1; ++j)
2486  minor[i][j] = t[i][j < k ? j : j + 1];
2487 
2488  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2489 
2490  det += t[dim - 1][k] * cofactor;
2491  }
2492 
2493  return ((dim % 2 == 0) ? 1. : -1.) * det;
2494 }
2495 
2501 template <typename Number>
2502 constexpr DEAL_II_ALWAYS_INLINE Number
2504 {
2505  return t[0][0];
2506 }
2507 
2513 template <typename Number>
2514 constexpr DEAL_II_ALWAYS_INLINE Number
2516 {
2517  // hard-coded for efficiency reasons
2518  return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2519 }
2520 
2526 template <typename Number>
2527 constexpr DEAL_II_ALWAYS_INLINE Number
2529 {
2530  // hard-coded for efficiency reasons
2531  const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2532  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2533  const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2534  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2535  const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2536  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2537  return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2538 }
2539 
2540 
2547 template <int dim, typename Number>
2550 {
2551  Number t = d[0][0];
2552  for (unsigned int i = 1; i < dim; ++i)
2553  t += d[i][i];
2554  return t;
2555 }
2556 
2557 
2566 template <int dim, typename Number>
2569 {
2570  Number return_tensor[dim][dim];
2571 
2572  // if desired, take over the
2573  // inversion of a 4x4 tensor
2574  // from the FullMatrix
2575  AssertThrow(false, ExcNotImplemented());
2576 
2577  return Tensor<2, dim, Number>(return_tensor);
2578 }
2579 
2580 
2581 #ifndef DOXYGEN
2582 
2583 template <typename Number>
2585  invert(const Tensor<2, 1, Number> &t)
2586 {
2587  Tensor<2, 1, Number> return_tensor;
2588 
2589  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2590 
2591  return return_tensor;
2592 }
2593 
2594 
2595 template <typename Number>
2597  invert(const Tensor<2, 2, Number> &t)
2598 {
2599  Tensor<2, 2, Number> return_tensor;
2600 
2601  const Number inv_det_t = internal::NumberType<Number>::value(
2602  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2603  return_tensor[0][0] = t[1][1];
2604  return_tensor[0][1] = -t[0][1];
2605  return_tensor[1][0] = -t[1][0];
2606  return_tensor[1][1] = t[0][0];
2607  return_tensor *= inv_det_t;
2608 
2609  return return_tensor;
2610 }
2611 
2612 
2613 template <typename Number>
2615  invert(const Tensor<2, 3, Number> &t)
2616 {
2617  Tensor<2, 3, Number> return_tensor;
2618 
2619  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2620  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2621  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2622  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2623  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2624  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2625  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2626  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2627  return_tensor[1][1] = internal::NumberType<Number>::value(t[0][0] * t[2][2]) -
2628  internal::NumberType<Number>::value(t[0][2] * t[2][0]);
2629  return_tensor[1][2] = internal::NumberType<Number>::value(t[0][2] * t[1][0]) -
2630  internal::NumberType<Number>::value(t[0][0] * t[1][2]);
2631  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2632  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2633  return_tensor[2][1] = internal::NumberType<Number>::value(t[0][1] * t[2][0]) -
2634  internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2635  return_tensor[2][2] = internal::NumberType<Number>::value(t[0][0] * t[1][1]) -
2636  internal::NumberType<Number>::value(t[0][1] * t[1][0]);
2637  const Number inv_det_t = internal::NumberType<Number>::value(
2638  1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2639  t[0][2] * return_tensor[2][0]));
2640  return_tensor *= inv_det_t;
2641 
2642  return return_tensor;
2643 }
2644 
2645 #endif /* DOXYGEN */
2646 
2647 
2653 template <int dim, typename Number>
2656 {
2658  for (unsigned int i = 0; i < dim; ++i)
2659  {
2660  tt[i][i] = t[i][i];
2661  for (unsigned int j = i + 1; j < dim; ++j)
2662  {
2663  tt[i][j] = t[j][i];
2664  tt[j][i] = t[i][j];
2665  };
2666  }
2667  return tt;
2668 }
2669 
2670 
2684 template <int dim, typename Number>
2685 constexpr Tensor<2, dim, Number>
2687 {
2688  return determinant(t) * invert(t);
2689 }
2690 
2691 
2705 template <int dim, typename Number>
2706 constexpr Tensor<2, dim, Number>
2708 {
2709  return transpose(adjugate(t));
2710 }
2711 
2712 
2776 template <int dim, typename Number>
2779 
2780 
2788 template <int dim, typename Number>
2789 inline Number
2791 {
2793  for (unsigned int j = 0; j < dim; ++j)
2794  {
2796  for (unsigned int i = 0; i < dim; ++i)
2797  sum += std::fabs(t[i][j]);
2798 
2799  if (sum > max)
2800  max = sum;
2801  }
2802 
2803  return max;
2804 }
2805 
2806 
2814 template <int dim, typename Number>
2815 inline Number
2817 {
2819  for (unsigned int i = 0; i < dim; ++i)
2820  {
2822  for (unsigned int j = 0; j < dim; ++j)
2823  sum += std::fabs(t[i][j]);
2824 
2825  if (sum > max)
2826  max = sum;
2827  }
2828 
2829  return max;
2830 }
2831 
2833 
2834 
2835 #ifndef DOXYGEN
2836 
2837 
2838 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2839 
2840 // Specialization of functions for ADOL-C number types when
2841 // the advanced branching feature is used
2842 template <int dim>
2843 inline adouble
2845 {
2847  for (unsigned int j = 0; j < dim; ++j)
2848  {
2850  for (unsigned int i = 0; i < dim; ++i)
2851  sum += std::fabs(t[i][j]);
2852 
2853  condassign(max, (sum > max), sum, max);
2854  }
2855 
2856  return max;
2857 }
2858 
2859 
2860 template <int dim>
2861 inline adouble
2863 {
2865  for (unsigned int i = 0; i < dim; ++i)
2866  {
2868  for (unsigned int j = 0; j < dim; ++j)
2869  sum += std::fabs(t[i][j]);
2870 
2871  condassign(max, (sum > max), sum, max);
2872  }
2873 
2874  return max;
2875 }
2876 
2877 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2878 
2879 
2880 #endif // DOXYGEN
2881 
2883 
2884 #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:2568
Number * begin_raw()
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1580
constexpr ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:2319
constexpr Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2707
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:2655
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:1751
constexpr Tensor operator-() const
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2416
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2790
#define AssertIndexRange(index, range)
Definition: exceptions.h:1648
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:2441
constexpr void clear()
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2816
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:489
#define AssertThrow(cond, exc)
Definition: exceptions.h:1533
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:542
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:1423
constexpr Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2549
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:369
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:94
typename Tensor< rank_ - 1, dim, Number >::tensor_type value_type
Definition: tensor.h:482
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:2240
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
Expression fabs(const Expression &x)
static constexpr unsigned int rank
Definition: tensor.h:468
constexpr Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2686
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:2356
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:1812
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > division_operator(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1919
Definition: tensor.h:448
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:1828
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
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:474
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:2031
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:2385
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:2167
Expression operator!=(const Expression &lhs, const Expression &rhs)
T max(const T &t, const MPI_Comm &mpi_communicator)
#define DEAL_II_CONSTEXPR
Definition: config.h:154
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:2475