Reference documentation for deal.II version Git a7d922cd5a 2021-01-15 21:19:32 -0700
\(\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 
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_ >= 0,
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 
556  operator[](const unsigned int i);
557 
563  constexpr DEAL_II_CUDA_HOST_DEV const value_type &
564  operator[](const unsigned int i) const;
565 
569  DEAL_II_CONSTEXPR const Number &
570  operator[](const TableIndices<rank_> &indices) const;
571 
575  DEAL_II_CONSTEXPR Number &operator[](const TableIndices<rank_> &indices);
576 
580  Number *
581  begin_raw();
582 
586  const Number *
587  begin_raw() const;
588 
592  Number *
593  end_raw();
594 
598  const Number *
599  end_raw() const;
600 
608  template <typename OtherNumber>
610  operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
611 
618  DEAL_II_CONSTEXPR Tensor &
619  operator=(const Number &d);
620 
624  template <typename OtherNumber>
625  DEAL_II_CONSTEXPR bool
627 
631  template <typename OtherNumber>
632  constexpr bool
634 
640  template <typename OtherNumber>
642  operator+=(const Tensor<rank_, dim, OtherNumber> &);
643 
649  template <typename OtherNumber>
651  operator-=(const Tensor<rank_, dim, OtherNumber> &);
652 
659  template <typename OtherNumber>
661  operator*=(const OtherNumber &factor);
662 
668  template <typename OtherNumber>
670  operator/=(const OtherNumber &factor);
671 
678  operator-() const;
679 
692  DEAL_II_CONSTEXPR void
693  clear();
694 
704  norm() const;
705 
714  norm_square() const;
715 
723  template <typename OtherNumber>
724  void
725  unroll(Vector<OtherNumber> &result) const;
726 
731  static DEAL_II_CONSTEXPR unsigned int
732  component_to_unrolled_index(const TableIndices<rank_> &indices);
733 
740  unrolled_to_component_indices(const unsigned int i);
741 
746  static constexpr std::size_t
748 
754  template <class Archive>
755  void
756  serialize(Archive &ar, const unsigned int version);
757 
763 
764 private:
768  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
769  // ... avoid a compiler warning in case of dim == 0 and ensure that the
770  // array always has positive size.
771 
775  template <typename OtherNumber>
776  void
777  unroll_recursion(Vector<OtherNumber> &result,
778  unsigned int & start_index) const;
779 
786  template <typename ArrayLike, std::size_t... Indices>
787  constexpr DEAL_II_CUDA_HOST_DEV
788  Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
789 
790  // Allow an arbitrary Tensor to access the underlying values.
791  template <int, int, typename>
792  friend class Tensor;
793 
794  // Point is allowed access to the coordinates. This is supposed to improve
795  // speed.
796  friend class Point<dim, Number>;
797 };
798 
799 
800 #ifndef DOXYGEN
801 namespace internal
802 {
803  // Workaround: The following 4 overloads are necessary to be able to
804  // compile the library with Apple Clang 8 and older. We should remove
805  // these overloads again when we bump the minimal required version to
806  // something later than clang-3.6 / Apple Clang 6.3.
807  template <int rank, int dim, typename T, typename U>
808  struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
809  {
810  using type =
812  };
813 
814  template <int rank, int dim, typename T, typename U>
815  struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
816  {
817  using type =
819  };
820 
821  template <typename T, int rank, int dim, typename U>
822  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
823  {
824  using type =
826  };
827 
828  template <int rank, int dim, typename T, typename U>
829  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
830  {
831  using type =
833  };
834  // end workaround
835 
840  template <int rank, int dim, typename T>
841  struct NumberType<Tensor<rank, dim, T>>
842  {
843  static constexpr DEAL_II_ALWAYS_INLINE const Tensor<rank, dim, T> &
844  value(const Tensor<rank, dim, T> &t)
845  {
846  return t;
847  }
848 
850  value(const T &t)
851  {
853  tmp = t;
854  return tmp;
855  }
856  };
857 } // namespace internal
858 
859 
860 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
861 
862 
863 template <int dim, typename Number>
866  // Some auto-differentiable numbers need explicit
867  // zero initialization such as adtl::adouble.
868  : Tensor{0.0}
869 {}
870 
871 
872 
873 template <int dim, typename Number>
874 template <typename OtherNumber>
876  Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
877  : value(internal::NumberType<Number>::value(initializer))
878 {}
879 
880 
881 
882 template <int dim, typename Number>
883 template <typename OtherNumber>
886  : Tensor{p.value}
887 {}
888 
889 
890 
891 template <int dim, typename Number>
892 inline Number *
894 {
895  return std::addressof(value);
896 }
897 
898 
899 
900 template <int dim, typename Number>
901 inline const Number *
903 {
904  return std::addressof(value);
905 }
906 
907 
908 
909 template <int dim, typename Number>
910 inline Number *
912 {
914 }
915 
916 
917 
918 template <int dim, typename Number>
919 const Number *
921 {
923 }
924 
925 
926 
927 template <int dim, typename Number>
930 {
931  // We cannot use Assert inside a CUDA kernel
932 # ifndef __CUDA_ARCH__
933  Assert(dim != 0,
934  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
935 # endif
936  return value;
937 }
938 
939 
940 template <int dim, typename Number>
943 {
944  // We cannot use Assert inside a CUDA kernel
945 # ifndef __CUDA_ARCH__
946  Assert(dim != 0,
947  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
948 # endif
949  return value;
950 }
951 
952 
953 template <int dim, typename Number>
954 template <typename OtherNumber>
958 {
960  return *this;
961 }
962 
963 
964 # ifdef __INTEL_COMPILER
965 template <int dim, typename Number>
969 {
970  value = p.value;
971  return *this;
972 }
973 # endif
974 
975 
976 template <int dim, typename Number>
977 template <typename OtherNumber>
980  Tensor<0, dim, Number>::operator=(const OtherNumber &d)
981 {
983  return *this;
984 }
985 
986 
987 template <int dim, typename Number>
988 template <typename OtherNumber>
989 DEAL_II_CONSTEXPR inline bool
991 {
992 # if defined(DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING)
993  Assert(!(std::is_same<Number, adouble>::value ||
994  std::is_same<OtherNumber, adouble>::value),
995  ExcMessage(
996  "The Tensor equality operator for ADOL-C taped numbers has not yet "
997  "been extended to support advanced branching."));
998 # endif
999 
1000  return numbers::values_are_equal(value, p.value);
1001 }
1002 
1003 
1004 template <int dim, typename Number>
1005 template <typename OtherNumber>
1006 constexpr bool
1008 {
1009  return !((*this) == p);
1010 }
1011 
1012 
1013 template <int dim, typename Number>
1014 template <typename OtherNumber>
1018 {
1019  value += p.value;
1020  return *this;
1021 }
1022 
1023 
1024 template <int dim, typename Number>
1025 template <typename OtherNumber>
1029 {
1030  value -= p.value;
1031  return *this;
1032 }
1033 
1034 
1035 
1036 namespace internal
1037 {
1038  namespace ComplexWorkaround
1039  {
1040  template <typename Number, typename OtherNumber>
1042  multiply_assign_scalar(Number &val, const OtherNumber &s)
1043  {
1044  val *= s;
1045  }
1046 
1047 # ifdef __CUDA_ARCH__
1048  template <typename Number, typename OtherNumber>
1050  multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
1051  {
1052  printf("This function is not implemented for std::complex<Number>!\n");
1053  assert(false);
1054  }
1055 # endif
1056  } // namespace ComplexWorkaround
1057 } // namespace internal
1058 
1059 
1060 template <int dim, typename Number>
1061 template <typename OtherNumber>
1064  Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1065 {
1066  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1067  return *this;
1068 }
1069 
1070 
1071 
1072 template <int dim, typename Number>
1073 template <typename OtherNumber>
1075 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1076 {
1077  value /= s;
1078  return *this;
1079 }
1080 
1081 
1082 template <int dim, typename Number>
1085 {
1086  return -value;
1087 }
1088 
1089 
1090 template <int dim, typename Number>
1093 {
1094  Assert(dim != 0,
1095  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1096  return numbers::NumberTraits<Number>::abs(value);
1097 }
1098 
1099 
1100 template <int dim, typename Number>
1104 {
1105  // We cannot use Assert inside a CUDA kernel
1106 # ifndef __CUDA_ARCH__
1107  Assert(dim != 0,
1108  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1109 # endif
1111 }
1112 
1113 
1114 template <int dim, typename Number>
1115 template <typename OtherNumber>
1116 inline void
1117 Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1118  unsigned int & index) const
1119 {
1120  Assert(dim != 0,
1121  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1122  result[index] = value;
1123  ++index;
1124 }
1125 
1126 
1127 template <int dim, typename Number>
1128 DEAL_II_CONSTEXPR inline void
1130 {
1131  // Some auto-differentiable numbers need explicit
1132  // zero initialization.
1134 }
1135 
1136 
1137 template <int dim, typename Number>
1138 template <class Archive>
1139 inline void
1140 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1141 {
1142  ar &value;
1143 }
1144 
1145 
1146 template <int dim, typename Number>
1148 
1149 
1150 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1151 
1152 template <int rank_, int dim, typename Number>
1153 template <typename ArrayLike, std::size_t... indices>
1155  Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1156  std::index_sequence<indices...>)
1157  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
1158 {
1159  static_assert(sizeof...(indices) == dim,
1160  "dim should match the number of indices");
1161 }
1162 
1163 
1164 
1165 template <int rank_, int dim, typename Number>
1168  // We would like to use =default, but this causes compile errors with some
1169  // MSVC versions and internal compiler errors with -O1 in gcc 5.4.
1170  : values{}
1171 {}
1172 
1173 
1174 
1175 template <int rank_, int dim, typename Number>
1178  : Tensor(initializer, std::make_index_sequence<dim>{})
1179 {}
1180 
1181 
1182 
1183 template <int rank_, int dim, typename Number>
1184 template <typename ElementType, typename MemorySpace>
1187  const ArrayView<ElementType, MemorySpace> &initializer)
1188 {
1190 
1191  for (unsigned int i = 0; i < n_independent_components; ++i)
1192  (*this)[unrolled_to_component_indices(i)] = initializer[i];
1193 }
1194 
1195 
1196 
1197 template <int rank_, int dim, typename Number>
1198 template <typename OtherNumber>
1201  const Tensor<rank_, dim, OtherNumber> &initializer)
1202  : Tensor(initializer, std::make_index_sequence<dim>{})
1203 {}
1204 
1205 
1206 template <int rank_, int dim, typename Number>
1207 template <typename OtherNumber>
1208 constexpr DEAL_II_ALWAYS_INLINE
1210  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1211  : Tensor(initializer, std::make_index_sequence<dim>{})
1212 {}
1213 
1214 
1215 template <int rank_, int dim, typename Number>
1216 template <typename OtherNumber>
1218  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1219 {
1220  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
1221 }
1222 
1223 
1224 
1225 namespace internal
1226 {
1227  namespace TensorSubscriptor
1228  {
1229  template <typename ArrayElementType, int dim>
1231  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1232  subscript(ArrayElementType * values,
1233  const unsigned int i,
1234  std::integral_constant<int, dim>)
1235  {
1236  // We cannot use Assert in a CUDA kernel
1237 # ifndef __CUDA_ARCH__
1238  AssertIndexRange(i, dim);
1239 # endif
1240  return values[i];
1241  }
1242 
1243  // The variables within this struct will be referenced in the next function.
1244  // It is a workaround that allows returning a reference to a static variable
1245  // while allowing constexpr evaluation of the function.
1246  // It has to be defined outside the function because constexpr functions
1247  // cannot define static variables
1248  template <typename ArrayElementType>
1249  struct Uninitialized
1250  {
1251  static ArrayElementType value;
1252  };
1253 
1254  template <typename Type>
1255  Type Uninitialized<Type>::value;
1256 
1257  template <typename ArrayElementType>
1259  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1260  subscript(ArrayElementType *,
1261  const unsigned int,
1262  std::integral_constant<int, 0>)
1263  {
1264  // We cannot use Assert in a CUDA kernel
1265 # ifndef __CUDA_ARCH__
1266  Assert(
1267  false,
1268  ExcMessage(
1269  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1270 # endif
1271  return Uninitialized<ArrayElementType>::value;
1272  }
1273  } // namespace TensorSubscriptor
1274 } // namespace internal
1275 
1276 
1277 template <int rank_, int dim, typename Number>
1280  operator[](const unsigned int i)
1281 {
1282  return ::internal::TensorSubscriptor::subscript(
1283  values, i, std::integral_constant<int, dim>());
1284 }
1285 
1286 
1287 template <int rank_, int dim, typename Number>
1288 constexpr DEAL_II_ALWAYS_INLINE
1290  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1291 {
1292 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1293  AssertIndexRange(i, dim);
1294 # endif
1295 
1296  return values[i];
1297 }
1298 
1299 
1300 template <int rank_, int dim, typename Number>
1301 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1303  operator[](const TableIndices<rank_> &indices) const
1304 {
1305 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1306  Assert(dim != 0,
1307  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1308 # endif
1309 
1310  return TensorAccessors::extract<rank_>(*this, indices);
1311 }
1312 
1313 
1314 
1315 template <int rank_, int dim, typename Number>
1318 {
1319 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1320  Assert(dim != 0,
1321  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1322 # endif
1323 
1324  return TensorAccessors::extract<rank_>(*this, indices);
1325 }
1326 
1327 
1328 
1329 template <int rank_, int dim, typename Number>
1330 inline Number *
1332 {
1333  return std::addressof(
1334  this->operator[](this->unrolled_to_component_indices(0)));
1335 }
1336 
1337 
1338 
1339 template <int rank_, int dim, typename Number>
1340 inline const Number *
1342 {
1343  return std::addressof(
1344  this->operator[](this->unrolled_to_component_indices(0)));
1345 }
1346 
1347 
1348 
1349 template <int rank_, int dim, typename Number>
1350 inline Number *
1352 {
1354 }
1355 
1356 
1357 
1358 template <int rank_, int dim, typename Number>
1359 inline const Number *
1361 {
1363 }
1364 
1365 
1366 
1367 template <int rank_, int dim, typename Number>
1368 template <typename OtherNumber>
1371 {
1372  // The following loop could be written more concisely using std::copy, but
1373  // that function is only constexpr from C++20 on.
1374  for (unsigned int i = 0; i < dim; ++i)
1375  values[i] = t.values[i];
1376  return *this;
1377 }
1378 
1379 
1380 template <int rank_, int dim, typename Number>
1383 {
1385  ExcMessage("Only assignment with zero is allowed"));
1386  (void)d;
1387 
1388  for (unsigned int i = 0; i < dim; ++i)
1390  return *this;
1391 }
1392 
1393 
1394 template <int rank_, int dim, typename Number>
1395 template <typename OtherNumber>
1396 DEAL_II_CONSTEXPR inline bool
1399 {
1400  for (unsigned int i = 0; i < dim; ++i)
1401  if (values[i] != p.values[i])
1402  return false;
1403  return true;
1404 }
1405 
1406 
1407 // At some places in the library, we have Point<0> for formal reasons
1408 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1409 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1410 // in the above function that the loop end check always fails, we
1411 // implement this function here
1412 template <>
1413 template <>
1414 DEAL_II_CONSTEXPR inline bool
1416 {
1417  return true;
1418 }
1419 
1420 
1421 template <int rank_, int dim, typename Number>
1422 template <typename OtherNumber>
1423 constexpr bool
1426 {
1427  return !((*this) == p);
1428 }
1429 
1430 
1431 template <int rank_, int dim, typename Number>
1432 template <typename OtherNumber>
1437 {
1438  for (unsigned int i = 0; i < dim; ++i)
1439  values[i] += p.values[i];
1440  return *this;
1441 }
1442 
1443 
1444 template <int rank_, int dim, typename Number>
1445 template <typename OtherNumber>
1450 {
1451  for (unsigned int i = 0; i < dim; ++i)
1452  values[i] -= p.values[i];
1453  return *this;
1454 }
1455 
1456 
1457 template <int rank_, int dim, typename Number>
1458 template <typename OtherNumber>
1461  Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1462 {
1463  for (unsigned int i = 0; i < dim; ++i)
1464  values[i] *= s;
1465  return *this;
1466 }
1467 
1468 
1469 namespace internal
1470 {
1471  namespace TensorImplementation
1472  {
1473  template <int rank,
1474  int dim,
1475  typename Number,
1476  typename OtherNumber,
1477  typename std::enable_if<
1478  !std::is_integral<
1479  typename ProductType<Number, OtherNumber>::type>::value &&
1480  !std::is_same<Number, Differentiation::SD::Expression>::value,
1481  int>::type = 0>
1484  const OtherNumber &factor)
1485  {
1486  const Number inverse_factor = Number(1.) / factor;
1487  // recurse over the base objects
1488  for (unsigned int d = 0; d < dim; ++d)
1489  t[d] *= inverse_factor;
1490  }
1491 
1492 
1493  template <int rank,
1494  int dim,
1495  typename Number,
1496  typename OtherNumber,
1497  typename std::enable_if<
1498  std::is_integral<
1499  typename ProductType<Number, OtherNumber>::type>::value ||
1500  std::is_same<Number, Differentiation::SD::Expression>::value,
1501  int>::type = 0>
1504  const OtherNumber &factor)
1505  {
1506  // recurse over the base objects
1507  for (unsigned int d = 0; d < dim; ++d)
1508  t[d] /= factor;
1509  }
1510  } // namespace TensorImplementation
1511 } // namespace internal
1512 
1513 
1514 template <int rank_, int dim, typename Number>
1515 template <typename OtherNumber>
1518  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1519 {
1521  return *this;
1522 }
1523 
1524 
1525 template <int rank_, int dim, typename Number>
1529 {
1531 
1532  for (unsigned int i = 0; i < dim; ++i)
1533  tmp.values[i] = -values[i];
1534 
1535  return tmp;
1536 }
1537 
1538 
1539 template <int rank_, int dim, typename Number>
1542 {
1543  return std::sqrt(norm_square());
1544 }
1545 
1546 
1547 template <int rank_, int dim, typename Number>
1551 {
1553  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1554  for (unsigned int i = 0; i < dim; ++i)
1555  s += values[i].norm_square();
1556 
1557  return s;
1558 }
1559 
1560 
1561 template <int rank_, int dim, typename Number>
1562 template <typename OtherNumber>
1563 inline void
1564 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
1565 {
1566  AssertDimension(result.size(),
1567  (Utilities::fixed_power<rank_, unsigned int>(dim)));
1568 
1569  unsigned int index = 0;
1570  unroll_recursion(result, index);
1571 }
1572 
1573 
1574 template <int rank_, int dim, typename Number>
1575 template <typename OtherNumber>
1576 inline void
1577 Tensor<rank_, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1578  unsigned int & index) const
1579 {
1580  for (unsigned int i = 0; i < dim; ++i)
1581  values[i].unroll_recursion(result, index);
1582 }
1583 
1584 
1585 template <int rank_, int dim, typename Number>
1586 DEAL_II_CONSTEXPR inline unsigned int
1588  const TableIndices<rank_> &indices)
1589 {
1590  unsigned int index = 0;
1591  for (int r = 0; r < rank_; ++r)
1592  index = index * dim + indices[r];
1593 
1594  return index;
1595 }
1596 
1597 
1598 
1599 namespace internal
1600 {
1601  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1602  // and rank=2. Make sure we don't have compiler warnings.
1603 
1604  template <int dim>
1605  inline DEAL_II_CONSTEXPR unsigned int
1606  mod(const unsigned int x)
1607  {
1608  return x % dim;
1609  }
1610 
1611  template <>
1612  inline unsigned int
1613  mod<0>(const unsigned int x)
1614  {
1615  Assert(false, ExcInternalError());
1616  return x;
1617  }
1618 
1619  template <int dim>
1620  inline DEAL_II_CONSTEXPR unsigned int
1621  div(const unsigned int x)
1622  {
1623  return x / dim;
1624  }
1625 
1626  template <>
1627  inline unsigned int
1628  div<0>(const unsigned int x)
1629  {
1630  Assert(false, ExcInternalError());
1631  return x;
1632  }
1633 
1634 } // namespace internal
1635 
1636 
1637 
1638 template <int rank_, int dim, typename Number>
1641 {
1642  AssertIndexRange(i, n_independent_components);
1643 
1644  TableIndices<rank_> indices;
1645 
1646  unsigned int remainder = i;
1647  for (int r = rank_ - 1; r >= 0; --r)
1648  {
1649  indices[r] = internal::mod<dim>(remainder);
1650  remainder = internal::div<dim>(remainder);
1651  }
1652  Assert(remainder == 0, ExcInternalError());
1653 
1654  return indices;
1655 }
1656 
1657 
1658 template <int rank_, int dim, typename Number>
1659 DEAL_II_CONSTEXPR inline void
1661 {
1662  for (unsigned int i = 0; i < dim; ++i)
1664 }
1665 
1666 
1667 template <int rank_, int dim, typename Number>
1668 constexpr std::size_t
1670 {
1671  return sizeof(Tensor<rank_, dim, Number>);
1672 }
1673 
1674 
1675 template <int rank_, int dim, typename Number>
1676 template <class Archive>
1677 inline void
1678 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1679 {
1680  ar &values;
1681 }
1682 
1683 
1684 template <int rank_, int dim, typename Number>
1686 
1687 #endif // DOXYGEN
1688 
1689 /* ----------------- Non-member functions operating on tensors. ------------ */
1690 
1695 
1703 template <int rank_, int dim, typename Number>
1704 inline std::ostream &
1705 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1706 {
1707  for (unsigned int i = 0; i < dim; ++i)
1708  {
1709  out << p[i];
1710  if (i != dim - 1)
1711  out << ' ';
1712  }
1713 
1714  return out;
1715 }
1716 
1717 
1724 template <int dim, typename Number>
1725 inline std::ostream &
1726 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1727 {
1728  out << static_cast<const Number &>(p);
1729  return out;
1730 }
1731 
1732 
1734 
1738 
1739 
1750 template <int dim, typename Number, typename Other>
1753  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1754 {
1755  return object * static_cast<const Number &>(t);
1756 }
1757 
1758 
1759 
1770 template <int dim, typename Number, typename Other>
1773  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1774 {
1775  return static_cast<const Number &>(t) * object;
1776 }
1777 
1778 
1790 template <int dim, typename Number, typename OtherNumber>
1794  const Tensor<0, dim, OtherNumber> &src2)
1795 {
1796  return static_cast<const Number &>(src1) *
1797  static_cast<const OtherNumber &>(src2);
1798 }
1799 
1800 
1808 template <int dim, typename Number, typename OtherNumber>
1810  Tensor<0,
1811  dim,
1812  typename ProductType<Number,
1813  typename EnableIfScalar<OtherNumber>::type>::type>
1814  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
1815 {
1816  return static_cast<const Number &>(t) / factor;
1817 }
1818 
1819 
1827 template <int dim, typename Number, typename OtherNumber>
1831  const Tensor<0, dim, OtherNumber> &q)
1832 {
1833  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1834 }
1835 
1836 
1844 template <int dim, typename Number, typename OtherNumber>
1848  const Tensor<0, dim, OtherNumber> &q)
1849 {
1850  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1851 }
1852 
1853 
1866 template <int rank, int dim, typename Number, typename OtherNumber>
1868  Tensor<rank,
1869  dim,
1870  typename ProductType<Number,
1871  typename EnableIfScalar<OtherNumber>::type>::type>
1872  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1873 {
1874  // recurse over the base objects
1876  for (unsigned int d = 0; d < dim; ++d)
1877  tt[d] = t[d] * factor;
1878  return tt;
1879 }
1880 
1881 
1894 template <int rank, int dim, typename Number, typename OtherNumber>
1896  Tensor<rank,
1897  dim,
1899  OtherNumber>::type>
1900  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
1901 {
1902  // simply forward to the operator above
1903  return t * factor;
1904 }
1905 
1906 
1907 namespace internal
1908 {
1909  namespace TensorImplementation
1910  {
1911  template <int rank,
1912  int dim,
1913  typename Number,
1914  typename OtherNumber,
1915  typename std::enable_if<
1916  !std::is_integral<
1917  typename ProductType<Number, OtherNumber>::type>::value,
1918  int>::type = 0>
1922  const OtherNumber & factor)
1923  {
1925  const Number inverse_factor = Number(1.) / factor;
1926  // recurse over the base objects
1927  for (unsigned int d = 0; d < dim; ++d)
1928  tt[d] = t[d] * inverse_factor;
1929  return tt;
1930  }
1931 
1932 
1933  template <int rank,
1934  int dim,
1935  typename Number,
1936  typename OtherNumber,
1937  typename std::enable_if<
1938  std::is_integral<
1939  typename ProductType<Number, OtherNumber>::type>::value,
1940  int>::type = 0>
1944  const OtherNumber & factor)
1945  {
1947  // recurse over the base objects
1948  for (unsigned int d = 0; d < dim; ++d)
1949  tt[d] = t[d] / factor;
1950  return tt;
1951  }
1952  } // namespace TensorImplementation
1953 } // namespace internal
1954 
1955 
1965 template <int rank, int dim, typename Number, typename OtherNumber>
1967  Tensor<rank,
1968  dim,
1969  typename ProductType<Number,
1970  typename EnableIfScalar<OtherNumber>::type>::type>
1971  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1972 {
1974 }
1975 
1976 
1986 template <int rank, int dim, typename Number, typename OtherNumber>
1991 {
1993 
1994  for (unsigned int i = 0; i < dim; ++i)
1995  tmp[i] += q[i];
1996 
1997  return tmp;
1998 }
1999 
2000 
2010 template <int rank, int dim, typename Number, typename OtherNumber>
2015 {
2017 
2018  for (unsigned int i = 0; i < dim; ++i)
2019  tmp[i] -= q[i];
2020 
2021  return tmp;
2022 }
2023 
2030 template <int dim, typename Number, typename OtherNumber>
2034  const Tensor<0, dim, OtherNumber> &src2)
2035 {
2037 
2038  tmp *= src2;
2039 
2040  return tmp;
2041 }
2042 
2059 template <int rank, int dim, typename Number, typename OtherNumber>
2063  const Tensor<rank, dim, OtherNumber> &src2)
2064 {
2066 
2067  for (unsigned int i = 0; i < dim; ++i)
2068  tmp[i] = schur_product(Tensor<rank - 1, dim, Number>(src1[i]),
2070 
2071  return tmp;
2072 }
2073 
2075 
2079 
2080 
2103 template <int rank_1,
2104  int rank_2,
2105  int dim,
2106  typename Number,
2107  typename OtherNumber,
2108  typename = typename std::enable_if<rank_1 >= 1 && rank_2 >= 1>::type>
2110  typename Tensor<rank_1 + rank_2 - 2,
2111  dim,
2115 {
2116  typename Tensor<rank_1 + rank_2 - 2,
2117  dim,
2119  result{};
2120 
2121  TensorAccessors::internal::
2122  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2123  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2124  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2125 
2126  return result;
2127 }
2128 
2129 
2158 template <int index_1,
2159  int index_2,
2160  int rank_1,
2161  int rank_2,
2162  int dim,
2163  typename Number,
2164  typename OtherNumber>
2166  typename Tensor<rank_1 + rank_2 - 2,
2167  dim,
2171 {
2172  Assert(0 <= index_1 && index_1 < rank_1,
2173  ExcMessage(
2174  "The specified index_1 must lie within the range [0,rank_1)"));
2175  Assert(0 <= index_2 && index_2 < rank_2,
2176  ExcMessage(
2177  "The specified index_2 must lie within the range [0,rank_2)"));
2178 
2179  using namespace TensorAccessors;
2180  using namespace TensorAccessors::internal;
2181 
2182  // Reorder index_1 to the end of src1:
2184  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2185 
2186  // Reorder index_2 to the end of src2:
2188  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2189 
2190  typename Tensor<rank_1 + rank_2 - 2,
2191  dim,
2193  result{};
2194  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2195  return result;
2196 }
2197 
2198 
2229 template <int index_1,
2230  int index_2,
2231  int index_3,
2232  int index_4,
2233  int rank_1,
2234  int rank_2,
2235  int dim,
2236  typename Number,
2237  typename OtherNumber>
2238 DEAL_II_CONSTEXPR inline
2239  typename Tensor<rank_1 + rank_2 - 4,
2240  dim,
2244 {
2245  Assert(0 <= index_1 && index_1 < rank_1,
2246  ExcMessage(
2247  "The specified index_1 must lie within the range [0,rank_1)"));
2248  Assert(0 <= index_3 && index_3 < rank_1,
2249  ExcMessage(
2250  "The specified index_3 must lie within the range [0,rank_1)"));
2251  Assert(index_1 != index_3,
2252  ExcMessage("index_1 and index_3 must not be the same"));
2253  Assert(0 <= index_2 && index_2 < rank_2,
2254  ExcMessage(
2255  "The specified index_2 must lie within the range [0,rank_2)"));
2256  Assert(0 <= index_4 && index_4 < rank_2,
2257  ExcMessage(
2258  "The specified index_4 must lie within the range [0,rank_2)"));
2259  Assert(index_2 != index_4,
2260  ExcMessage("index_2 and index_4 must not be the same"));
2261 
2262  using namespace TensorAccessors;
2263  using namespace TensorAccessors::internal;
2264 
2265  // Reorder index_1 to the end of src1:
2267  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2268 
2269  // Reorder index_2 to the end of src2:
2271  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2272 
2273  // Now, reorder index_3 to the end of src1. We have to make sure to
2274  // preserve the original ordering: index_1 has been removed. If
2275  // index_3 > index_1, we have to use (index_3 - 1) instead:
2277  (index_3 < index_1 ? index_3 : index_3 - 1),
2278  rank_1,
2279  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2280  reord_3 =
2281  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2282  index_3 - 1,
2283  rank_1 > (reord_1);
2284 
2285  // Now, reorder index_4 to the end of src2. We have to make sure to
2286  // preserve the original ordering: index_2 has been removed. If
2287  // index_4 > index_2, we have to use (index_4 - 1) instead:
2288  ReorderedIndexView<
2289  (index_4 < index_2 ? index_4 : index_4 - 1),
2290  rank_2,
2291  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2292  reord_4 =
2293  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2294  index_4 - 1,
2295  rank_2 > (reord_2);
2296 
2297  typename Tensor<rank_1 + rank_2 - 4,
2298  dim,
2300  result{};
2301  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2302  return result;
2303 }
2304 
2305 
2318 template <int rank, int dim, typename Number, typename OtherNumber>
2322  const Tensor<rank, dim, OtherNumber> &right)
2323 {
2324  typename ProductType<Number, OtherNumber>::type result{};
2325  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2326  return result;
2327 }
2328 
2329 
2347 template <template <int, int, typename> class TensorT1,
2348  template <int, int, typename> class TensorT2,
2349  template <int, int, typename> class TensorT3,
2350  int rank_1,
2351  int rank_2,
2352  int dim,
2353  typename T1,
2354  typename T2,
2355  typename T3>
2358  contract3(const TensorT1<rank_1, dim, T1> & left,
2359  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2360  const TensorT3<rank_2, dim, T3> & right)
2361 {
2362  using return_type =
2364  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2365  middle,
2366  right);
2367 }
2368 
2369 
2380 template <int rank_1,
2381  int rank_2,
2382  int dim,
2383  typename Number,
2384  typename OtherNumber>
2389 {
2390  typename Tensor<rank_1 + rank_2,
2391  dim,
2393  result{};
2394  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2395  return result;
2396 }
2397 
2398 
2400 
2404 
2405 
2416 template <int dim, typename Number>
2419 {
2420  Assert(dim == 2, ExcInternalError());
2421 
2422  Tensor<1, dim, Number> result;
2423 
2424  result[0] = src[1];
2425  result[1] = -src[0];
2426 
2427  return result;
2428 }
2429 
2430 
2440 template <int dim, typename Number1, typename Number2>
2444  const Tensor<1, dim, Number2> &src2)
2445 {
2446  Assert(dim == 3, ExcInternalError());
2447 
2449 
2450  // avoid compiler warnings
2451  constexpr int s0 = 0 % dim;
2452  constexpr int s1 = 1 % dim;
2453  constexpr int s2 = 2 % dim;
2454 
2455  result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2456  result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2457  result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2458 
2459  return result;
2460 }
2461 
2462 
2464 
2468 
2469 
2475 template <int dim, typename Number>
2478 {
2479  // Compute the determinant using the Laplace expansion of the
2480  // determinant. We expand along the last row.
2481  Number det = internal::NumberType<Number>::value(0.0);
2482 
2483  for (unsigned int k = 0; k < dim; ++k)
2484  {
2485  Tensor<2, dim - 1, Number> minor;
2486  for (unsigned int i = 0; i < dim - 1; ++i)
2487  for (unsigned int j = 0; j < dim - 1; ++j)
2488  minor[i][j] = t[i][j < k ? j : j + 1];
2489 
2490  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2491 
2492  det += t[dim - 1][k] * cofactor;
2493  }
2494 
2495  return ((dim % 2 == 0) ? 1. : -1.) * det;
2496 }
2497 
2503 template <typename Number>
2504 constexpr DEAL_II_ALWAYS_INLINE Number
2506 {
2507  return t[0][0];
2508 }
2509 
2515 template <typename Number>
2516 constexpr DEAL_II_ALWAYS_INLINE Number
2518 {
2519  // hard-coded for efficiency reasons
2520  return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2521 }
2522 
2528 template <typename Number>
2529 constexpr DEAL_II_ALWAYS_INLINE Number
2531 {
2532  // hard-coded for efficiency reasons
2533  const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2534  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2535  const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2536  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2537  const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2538  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2539  return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2540 }
2541 
2542 
2549 template <int dim, typename Number>
2552 {
2553  Number t = d[0][0];
2554  for (unsigned int i = 1; i < dim; ++i)
2555  t += d[i][i];
2556  return t;
2557 }
2558 
2559 
2568 template <int dim, typename Number>
2571 {
2572  Number return_tensor[dim][dim];
2573 
2574  // if desired, take over the
2575  // inversion of a 4x4 tensor
2576  // from the FullMatrix
2577  AssertThrow(false, ExcNotImplemented());
2578 
2579  return Tensor<2, dim, Number>(return_tensor);
2580 }
2581 
2582 
2583 #ifndef DOXYGEN
2584 
2585 template <typename Number>
2587  invert(const Tensor<2, 1, Number> &t)
2588 {
2589  Tensor<2, 1, Number> return_tensor;
2590 
2591  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2592 
2593  return return_tensor;
2594 }
2595 
2596 
2597 template <typename Number>
2599  invert(const Tensor<2, 2, Number> &t)
2600 {
2601  Tensor<2, 2, Number> return_tensor;
2602 
2603  const Number inv_det_t = internal::NumberType<Number>::value(
2604  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2605  return_tensor[0][0] = t[1][1];
2606  return_tensor[0][1] = -t[0][1];
2607  return_tensor[1][0] = -t[1][0];
2608  return_tensor[1][1] = t[0][0];
2609  return_tensor *= inv_det_t;
2610 
2611  return return_tensor;
2612 }
2613 
2614 
2615 template <typename Number>
2617  invert(const Tensor<2, 3, Number> &t)
2618 {
2619  Tensor<2, 3, Number> return_tensor;
2620 
2621  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2622  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2623  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2624  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2625  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2626  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2627  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2628  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2629  return_tensor[1][1] = internal::NumberType<Number>::value(t[0][0] * t[2][2]) -
2630  internal::NumberType<Number>::value(t[0][2] * t[2][0]);
2631  return_tensor[1][2] = internal::NumberType<Number>::value(t[0][2] * t[1][0]) -
2632  internal::NumberType<Number>::value(t[0][0] * t[1][2]);
2633  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2634  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2635  return_tensor[2][1] = internal::NumberType<Number>::value(t[0][1] * t[2][0]) -
2636  internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2637  return_tensor[2][2] = internal::NumberType<Number>::value(t[0][0] * t[1][1]) -
2638  internal::NumberType<Number>::value(t[0][1] * t[1][0]);
2639  const Number inv_det_t = internal::NumberType<Number>::value(
2640  1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2641  t[0][2] * return_tensor[2][0]));
2642  return_tensor *= inv_det_t;
2643 
2644  return return_tensor;
2645 }
2646 
2647 #endif /* DOXYGEN */
2648 
2649 
2655 template <int dim, typename Number>
2658 {
2660  for (unsigned int i = 0; i < dim; ++i)
2661  {
2662  tt[i][i] = t[i][i];
2663  for (unsigned int j = i + 1; j < dim; ++j)
2664  {
2665  tt[i][j] = t[j][i];
2666  tt[j][i] = t[i][j];
2667  };
2668  }
2669  return tt;
2670 }
2671 
2672 
2686 template <int dim, typename Number>
2687 constexpr Tensor<2, dim, Number>
2689 {
2690  return determinant(t) * invert(t);
2691 }
2692 
2693 
2707 template <int dim, typename Number>
2708 constexpr Tensor<2, dim, Number>
2710 {
2711  return transpose(adjugate(t));
2712 }
2713 
2714 
2778 template <int dim, typename Number>
2781 
2782 
2790 template <int dim, typename Number>
2791 inline Number
2793 {
2795  for (unsigned int j = 0; j < dim; ++j)
2796  {
2798  for (unsigned int i = 0; i < dim; ++i)
2799  sum += std::fabs(t[i][j]);
2800 
2801  if (sum > max)
2802  max = sum;
2803  }
2804 
2805  return max;
2806 }
2807 
2808 
2816 template <int dim, typename Number>
2817 inline Number
2819 {
2821  for (unsigned int i = 0; i < dim; ++i)
2822  {
2824  for (unsigned int j = 0; j < dim; ++j)
2825  sum += std::fabs(t[i][j]);
2826 
2827  if (sum > max)
2828  max = sum;
2829  }
2830 
2831  return max;
2832 }
2833 
2835 
2836 
2837 #ifndef DOXYGEN
2838 
2839 
2840 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2841 
2842 // Specialization of functions for ADOL-C number types when
2843 // the advanced branching feature is used
2844 template <int dim>
2845 inline adouble
2847 {
2849  for (unsigned int j = 0; j < dim; ++j)
2850  {
2852  for (unsigned int i = 0; i < dim; ++i)
2853  sum += std::fabs(t[i][j]);
2854 
2855  condassign(max, (sum > max), sum, max);
2856  }
2857 
2858  return max;
2859 }
2860 
2861 
2862 template <int dim>
2863 inline adouble
2865 {
2867  for (unsigned int i = 0; i < dim; ++i)
2868  {
2870  for (unsigned int j = 0; j < dim; ++j)
2871  sum += std::fabs(t[i][j]);
2872 
2873  condassign(max, (sum > max), sum, max);
2874  }
2875 
2876  return max;
2877 }
2878 
2879 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2880 
2881 
2882 #endif // DOXYGEN
2883 
2885 
2886 #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:2570
Number * begin_raw()
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
constexpr ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:2321
constexpr Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2709
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:2657
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:1753
constexpr Tensor operator-() const
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2418
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2792
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
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:2443
constexpr void clear()
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2818
constexpr Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:768
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:1576
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:762
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:1466
constexpr Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2551
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:380
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:94
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:2242
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:2688
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:2358
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:1814
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > division_operator(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1921
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:1830
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:379
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:2033
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:2387
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:2169
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:155
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:2477