Reference documentation for deal.II version Git 409ee4b167 2020-08-14 09:46:12 -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 #ifdef DEAL_II_MSVC
499  : values{}
500  {}
501 #else
502  = default;
503 #endif
504 
510  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
511  const array_type &initializer);
512 
526  template <typename ElementType, typename MemorySpace>
527  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
528  const ArrayView<ElementType, MemorySpace> &initializer);
529 
537  template <typename OtherNumber>
538  constexpr DEAL_II_CUDA_HOST_DEV
539  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
540 
544  template <typename OtherNumber>
545  constexpr Tensor(
546  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
547 
551  template <typename OtherNumber>
552  constexpr
553  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
554 
561  operator[](const unsigned int i);
562 
568  constexpr DEAL_II_CUDA_HOST_DEV const value_type &
569  operator[](const unsigned int i) const;
570 
574  DEAL_II_CONSTEXPR const Number &
575  operator[](const TableIndices<rank_> &indices) const;
576 
580  DEAL_II_CONSTEXPR Number &operator[](const TableIndices<rank_> &indices);
581 
585  Number *
586  begin_raw();
587 
591  const Number *
592  begin_raw() const;
593 
597  Number *
598  end_raw();
599 
603  const Number *
604  end_raw() const;
605 
613  template <typename OtherNumber>
615  operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
616 
623  DEAL_II_CONSTEXPR Tensor &
624  operator=(const Number &d);
625 
629  template <typename OtherNumber>
630  DEAL_II_CONSTEXPR bool
632 
636  template <typename OtherNumber>
637  constexpr bool
639 
645  template <typename OtherNumber>
647  operator+=(const Tensor<rank_, dim, OtherNumber> &);
648 
654  template <typename OtherNumber>
656  operator-=(const Tensor<rank_, dim, OtherNumber> &);
657 
664  template <typename OtherNumber>
666  operator*=(const OtherNumber &factor);
667 
673  template <typename OtherNumber>
675  operator/=(const OtherNumber &factor);
676 
683  operator-() const;
684 
697  DEAL_II_CONSTEXPR void
698  clear();
699 
709  norm() const;
710 
719  norm_square() const;
720 
728  template <typename OtherNumber>
729  void
730  unroll(Vector<OtherNumber> &result) const;
731 
736  static DEAL_II_CONSTEXPR unsigned int
737  component_to_unrolled_index(const TableIndices<rank_> &indices);
738 
745  unrolled_to_component_indices(const unsigned int i);
746 
751  static constexpr std::size_t
753 
758  template <class Archive>
759  void
760  serialize(Archive &ar, const unsigned int version);
761 
767 
768 private:
772  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
773  // ... avoid a compiler warning in case of dim == 0 and ensure that the
774  // array always has positive size.
775 
779  template <typename OtherNumber>
780  void
781  unroll_recursion(Vector<OtherNumber> &result,
782  unsigned int & start_index) const;
783 
790  template <typename ArrayLike, std::size_t... Indices>
791  constexpr DEAL_II_CUDA_HOST_DEV
792  Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
793 
794  // Allow an arbitrary Tensor to access the underlying values.
795  template <int, int, typename>
796  friend class Tensor;
797 
798  // Point is allowed access to the coordinates. This is supposed to improve
799  // speed.
800  friend class Point<dim, Number>;
801 };
802 
803 
804 #ifndef DOXYGEN
805 namespace internal
806 {
807  // Workaround: The following 4 overloads are necessary to be able to
808  // compile the library with Apple Clang 8 and older. We should remove
809  // these overloads again when we bump the minimal required version to
810  // something later than clang-3.6 / Apple Clang 6.3.
811  template <int rank, int dim, typename T, typename U>
812  struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
813  {
814  using type =
816  };
817 
818  template <int rank, int dim, typename T, typename U>
819  struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
820  {
821  using type =
823  };
824 
825  template <typename T, int rank, int dim, typename U>
826  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
827  {
828  using type =
830  };
831 
832  template <int rank, int dim, typename T, typename U>
833  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
834  {
835  using type =
837  };
838  // end workaround
839 
844  template <int rank, int dim, typename T>
845  struct NumberType<Tensor<rank, dim, T>>
846  {
847  static constexpr DEAL_II_ALWAYS_INLINE const Tensor<rank, dim, T> &
848  value(const Tensor<rank, dim, T> &t)
849  {
850  return t;
851  }
852 
854  value(const T &t)
855  {
857  tmp = t;
858  return tmp;
859  }
860  };
861 } // namespace internal
862 
863 
864 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
865 
866 
867 template <int dim, typename Number>
870  // Some auto-differentiable numbers need explicit
871  // zero initialization such as adtl::adouble.
872  : Tensor{0.0}
873 {}
874 
875 
876 
877 template <int dim, typename Number>
878 template <typename OtherNumber>
880  Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
881  : value(internal::NumberType<Number>::value(initializer))
882 {}
883 
884 
885 
886 template <int dim, typename Number>
887 template <typename OtherNumber>
890  : Tensor{p.value}
891 {}
892 
893 
894 
895 template <int dim, typename Number>
896 inline Number *
898 {
899  return std::addressof(value);
900 }
901 
902 
903 
904 template <int dim, typename Number>
905 inline const Number *
907 {
908  return std::addressof(value);
909 }
910 
911 
912 
913 template <int dim, typename Number>
914 inline Number *
916 {
918 }
919 
920 
921 
922 template <int dim, typename Number>
923 const Number *
925 {
927 }
928 
929 
930 
931 template <int dim, typename Number>
934 {
935  // We cannot use Assert inside a CUDA kernel
936 # ifndef __CUDA_ARCH__
937  Assert(dim != 0,
938  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
939 # endif
940  return value;
941 }
942 
943 
944 template <int dim, typename Number>
947 {
948  // We cannot use Assert inside a CUDA kernel
949 # ifndef __CUDA_ARCH__
950  Assert(dim != 0,
951  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
952 # endif
953  return value;
954 }
955 
956 
957 template <int dim, typename Number>
958 template <typename OtherNumber>
962 {
964  return *this;
965 }
966 
967 
968 # ifdef __INTEL_COMPILER
969 template <int dim, typename Number>
973 {
974  value = p.value;
975  return *this;
976 }
977 # endif
978 
979 
980 template <int dim, typename Number>
981 template <typename OtherNumber>
984  Tensor<0, dim, Number>::operator=(const OtherNumber &d)
985 {
987  return *this;
988 }
989 
990 
991 template <int dim, typename Number>
992 template <typename OtherNumber>
993 DEAL_II_CONSTEXPR inline bool
995 {
996 # if defined(DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING)
997  Assert(!(std::is_same<Number, adouble>::value ||
998  std::is_same<OtherNumber, adouble>::value),
999  ExcMessage(
1000  "The Tensor equality operator for ADOL-C taped numbers has not yet "
1001  "been extended to support advanced branching."));
1002 # endif
1003 
1004  return numbers::values_are_equal(value, p.value);
1005 }
1006 
1007 
1008 template <int dim, typename Number>
1009 template <typename OtherNumber>
1010 constexpr bool
1012 {
1013  return !((*this) == p);
1014 }
1015 
1016 
1017 template <int dim, typename Number>
1018 template <typename OtherNumber>
1022 {
1023  value += p.value;
1024  return *this;
1025 }
1026 
1027 
1028 template <int dim, typename Number>
1029 template <typename OtherNumber>
1033 {
1034  value -= p.value;
1035  return *this;
1036 }
1037 
1038 
1039 
1040 namespace internal
1041 {
1042  namespace ComplexWorkaround
1043  {
1044  template <typename Number, typename OtherNumber>
1046  multiply_assign_scalar(Number &val, const OtherNumber &s)
1047  {
1048  val *= s;
1049  }
1050 
1051 # ifdef __CUDA_ARCH__
1052  template <typename Number, typename OtherNumber>
1054  multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
1055  {
1056  printf("This function is not implemented for std::complex<Number>!\n");
1057  assert(false);
1058  }
1059 # endif
1060  } // namespace ComplexWorkaround
1061 } // namespace internal
1062 
1063 
1064 template <int dim, typename Number>
1065 template <typename OtherNumber>
1068  Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1069 {
1070  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1071  return *this;
1072 }
1073 
1074 
1075 
1076 template <int dim, typename Number>
1077 template <typename OtherNumber>
1079 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1080 {
1081  value /= s;
1082  return *this;
1083 }
1084 
1085 
1086 template <int dim, typename Number>
1089 {
1090  return -value;
1091 }
1092 
1093 
1094 template <int dim, typename Number>
1097 {
1098  Assert(dim != 0,
1099  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1100  return numbers::NumberTraits<Number>::abs(value);
1101 }
1102 
1103 
1104 template <int dim, typename Number>
1108 {
1109  // We cannot use Assert inside a CUDA kernel
1110 # ifndef __CUDA_ARCH__
1111  Assert(dim != 0,
1112  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1113 # endif
1115 }
1116 
1117 
1118 template <int dim, typename Number>
1119 template <typename OtherNumber>
1120 inline void
1121 Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1122  unsigned int & index) const
1123 {
1124  Assert(dim != 0,
1125  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1126  result[index] = value;
1127  ++index;
1128 }
1129 
1130 
1131 template <int dim, typename Number>
1132 DEAL_II_CONSTEXPR inline void
1134 {
1135  // Some auto-differentiable numbers need explicit
1136  // zero initialization.
1138 }
1139 
1140 
1141 template <int dim, typename Number>
1142 template <class Archive>
1143 inline void
1144 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1145 {
1146  ar &value;
1147 }
1148 
1149 
1150 template <int dim, typename Number>
1152 
1153 
1154 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1155 
1156 template <int rank_, int dim, typename Number>
1157 template <typename ArrayLike, std::size_t... indices>
1159  Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1160  std::index_sequence<indices...>)
1161  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
1162 {
1163  static_assert(sizeof...(indices) == dim,
1164  "dim should match the number of indices");
1165 }
1166 
1167 
1168 template <int rank_, int dim, typename Number>
1171  : Tensor(initializer, std::make_index_sequence<dim>{})
1172 {}
1173 
1174 
1175 
1176 template <int rank_, int dim, typename Number>
1177 template <typename ElementType, typename MemorySpace>
1180  const ArrayView<ElementType, MemorySpace> &initializer)
1181 {
1183 
1184  for (unsigned int i = 0; i < n_independent_components; ++i)
1185  (*this)[unrolled_to_component_indices(i)] = initializer[i];
1186 }
1187 
1188 
1189 
1190 template <int rank_, int dim, typename Number>
1191 template <typename OtherNumber>
1194  const Tensor<rank_, dim, OtherNumber> &initializer)
1195  : Tensor(initializer, std::make_index_sequence<dim>{})
1196 {}
1197 
1198 
1199 template <int rank_, int dim, typename Number>
1200 template <typename OtherNumber>
1201 constexpr DEAL_II_ALWAYS_INLINE
1203  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1204  : Tensor(initializer, std::make_index_sequence<dim>{})
1205 {}
1206 
1207 
1208 template <int rank_, int dim, typename Number>
1209 template <typename OtherNumber>
1211  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1212 {
1213  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
1214 }
1215 
1216 
1217 
1218 namespace internal
1219 {
1220  namespace TensorSubscriptor
1221  {
1222  template <typename ArrayElementType, int dim>
1224  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1225  subscript(ArrayElementType * values,
1226  const unsigned int i,
1227  std::integral_constant<int, dim>)
1228  {
1229  // We cannot use Assert in a CUDA kernel
1230 # ifndef __CUDA_ARCH__
1231  AssertIndexRange(i, dim);
1232 # endif
1233  return values[i];
1234  }
1235 
1236  // The variables within this struct will be referenced in the next function.
1237  // It is a workaround that allows returning a reference to a static variable
1238  // while allowing constexpr evaluation of the function.
1239  // It has to be defined outside the function because constexpr functions
1240  // cannot define static variables
1241  template <typename ArrayElementType>
1242  struct Uninitialized
1243  {
1244  static ArrayElementType value;
1245  };
1246 
1247  template <typename Type>
1248  Type Uninitialized<Type>::value;
1249 
1250  template <typename ArrayElementType>
1252  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1253  subscript(ArrayElementType *,
1254  const unsigned int,
1255  std::integral_constant<int, 0>)
1256  {
1257  // We cannot use Assert in a CUDA kernel
1258 # ifndef __CUDA_ARCH__
1259  Assert(
1260  false,
1261  ExcMessage(
1262  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1263 # endif
1264  return Uninitialized<ArrayElementType>::value;
1265  }
1266  } // namespace TensorSubscriptor
1267 } // namespace internal
1268 
1269 
1270 template <int rank_, int dim, typename Number>
1273  operator[](const unsigned int i)
1274 {
1275  return ::internal::TensorSubscriptor::subscript(
1276  values, i, std::integral_constant<int, dim>());
1277 }
1278 
1279 
1280 template <int rank_, int dim, typename Number>
1281 constexpr DEAL_II_ALWAYS_INLINE
1283  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1284 {
1285  AssertIndexRange(i, dim);
1286 
1287  return values[i];
1288 }
1289 
1290 
1291 template <int rank_, int dim, typename Number>
1292 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1294  operator[](const TableIndices<rank_> &indices) const
1295 {
1296  Assert(dim != 0,
1297  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1298 
1299  return TensorAccessors::extract<rank_>(*this, indices);
1300 }
1301 
1302 
1303 
1304 template <int rank_, int dim, typename Number>
1307 {
1308  Assert(dim != 0,
1309  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1310 
1311  return TensorAccessors::extract<rank_>(*this, indices);
1312 }
1313 
1314 
1315 
1316 template <int rank_, int dim, typename Number>
1317 inline Number *
1319 {
1320  return std::addressof(
1321  this->operator[](this->unrolled_to_component_indices(0)));
1322 }
1323 
1324 
1325 
1326 template <int rank_, int dim, typename Number>
1327 inline const Number *
1329 {
1330  return std::addressof(
1331  this->operator[](this->unrolled_to_component_indices(0)));
1332 }
1333 
1334 
1335 
1336 template <int rank_, int dim, typename Number>
1337 inline Number *
1339 {
1341 }
1342 
1343 
1344 
1345 template <int rank_, int dim, typename Number>
1346 inline const Number *
1348 {
1350 }
1351 
1352 
1353 
1354 template <int rank_, int dim, typename Number>
1355 template <typename OtherNumber>
1358 {
1359  // The following loop could be written more concisely using std::copy, but
1360  // that function is only constexpr from C++20 on.
1361  for (unsigned int i = 0; i < dim; ++i)
1362  values[i] = t.values[i];
1363  return *this;
1364 }
1365 
1366 
1367 template <int rank_, int dim, typename Number>
1370 {
1372  ExcMessage("Only assignment with zero is allowed"));
1373  (void)d;
1374 
1375  for (unsigned int i = 0; i < dim; ++i)
1377  return *this;
1378 }
1379 
1380 
1381 template <int rank_, int dim, typename Number>
1382 template <typename OtherNumber>
1383 DEAL_II_CONSTEXPR inline bool
1386 {
1387  for (unsigned int i = 0; i < dim; ++i)
1388  if (values[i] != p.values[i])
1389  return false;
1390  return true;
1391 }
1392 
1393 
1394 // At some places in the library, we have Point<0> for formal reasons
1395 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1396 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1397 // in the above function that the loop end check always fails, we
1398 // implement this function here
1399 template <>
1400 template <>
1401 DEAL_II_CONSTEXPR inline bool
1403 {
1404  return true;
1405 }
1406 
1407 
1408 template <int rank_, int dim, typename Number>
1409 template <typename OtherNumber>
1410 constexpr bool
1413 {
1414  return !((*this) == p);
1415 }
1416 
1417 
1418 template <int rank_, int dim, typename Number>
1419 template <typename OtherNumber>
1424 {
1425  for (unsigned int i = 0; i < dim; ++i)
1426  values[i] += p.values[i];
1427  return *this;
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>
1448  Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1449 {
1450  for (unsigned int i = 0; i < dim; ++i)
1451  values[i] *= s;
1452  return *this;
1453 }
1454 
1455 
1456 namespace internal
1457 {
1458  namespace TensorImplementation
1459  {
1460  template <int rank,
1461  int dim,
1462  typename Number,
1463  typename OtherNumber,
1464  typename std::enable_if<
1465  !std::is_integral<
1466  typename ProductType<Number, OtherNumber>::type>::value &&
1467  !std::is_same<Number, Differentiation::SD::Expression>::value,
1468  int>::type = 0>
1471  const OtherNumber &factor)
1472  {
1473  const Number inverse_factor = Number(1.) / factor;
1474  // recurse over the base objects
1475  for (unsigned int d = 0; d < dim; ++d)
1476  t[d] *= inverse_factor;
1477  }
1478 
1479 
1480  template <int rank,
1481  int dim,
1482  typename Number,
1483  typename OtherNumber,
1484  typename std::enable_if<
1485  std::is_integral<
1486  typename ProductType<Number, OtherNumber>::type>::value ||
1487  std::is_same<Number, Differentiation::SD::Expression>::value,
1488  int>::type = 0>
1491  const OtherNumber &factor)
1492  {
1493  // recurse over the base objects
1494  for (unsigned int d = 0; d < dim; ++d)
1495  t[d] /= factor;
1496  }
1497  } // namespace TensorImplementation
1498 } // namespace internal
1499 
1500 
1501 template <int rank_, int dim, typename Number>
1502 template <typename OtherNumber>
1505  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1506 {
1508  return *this;
1509 }
1510 
1511 
1512 template <int rank_, int dim, typename Number>
1516 {
1518 
1519  for (unsigned int i = 0; i < dim; ++i)
1520  tmp.values[i] = -values[i];
1521 
1522  return tmp;
1523 }
1524 
1525 
1526 template <int rank_, int dim, typename Number>
1529 {
1530  return std::sqrt(norm_square());
1531 }
1532 
1533 
1534 template <int rank_, int dim, typename Number>
1538 {
1540  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1541  for (unsigned int i = 0; i < dim; ++i)
1542  s += values[i].norm_square();
1543 
1544  return s;
1545 }
1546 
1547 
1548 template <int rank_, int dim, typename Number>
1549 template <typename OtherNumber>
1550 inline void
1551 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
1552 {
1553  AssertDimension(result.size(),
1554  (Utilities::fixed_power<rank_, unsigned int>(dim)));
1555 
1556  unsigned int index = 0;
1557  unroll_recursion(result, index);
1558 }
1559 
1560 
1561 template <int rank_, int dim, typename Number>
1562 template <typename OtherNumber>
1563 inline void
1564 Tensor<rank_, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1565  unsigned int & index) const
1566 {
1567  for (unsigned int i = 0; i < dim; ++i)
1568  values[i].unroll_recursion(result, index);
1569 }
1570 
1571 
1572 template <int rank_, int dim, typename Number>
1573 DEAL_II_CONSTEXPR inline unsigned int
1575  const TableIndices<rank_> &indices)
1576 {
1577  unsigned int index = 0;
1578  for (int r = 0; r < rank_; ++r)
1579  index = index * dim + indices[r];
1580 
1581  return index;
1582 }
1583 
1584 
1585 
1586 namespace internal
1587 {
1588  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1589  // and rank=2. Make sure we don't have compiler warnings.
1590 
1591  template <int dim>
1592  inline DEAL_II_CONSTEXPR unsigned int
1593  mod(const unsigned int x)
1594  {
1595  return x % dim;
1596  }
1597 
1598  template <>
1599  inline unsigned int
1600  mod<0>(const unsigned int x)
1601  {
1602  Assert(false, ExcInternalError());
1603  return x;
1604  }
1605 
1606  template <int dim>
1607  inline DEAL_II_CONSTEXPR unsigned int
1608  div(const unsigned int x)
1609  {
1610  return x / dim;
1611  }
1612 
1613  template <>
1614  inline unsigned int
1615  div<0>(const unsigned int x)
1616  {
1617  Assert(false, ExcInternalError());
1618  return x;
1619  }
1620 
1621 } // namespace internal
1622 
1623 
1624 
1625 template <int rank_, int dim, typename Number>
1628 {
1629  AssertIndexRange(i, n_independent_components);
1630 
1631  TableIndices<rank_> indices;
1632 
1633  unsigned int remainder = i;
1634  for (int r = rank_ - 1; r >= 0; --r)
1635  {
1636  indices[r] = internal::mod<dim>(remainder);
1637  remainder = internal::div<dim>(remainder);
1638  }
1639  Assert(remainder == 0, ExcInternalError());
1640 
1641  return indices;
1642 }
1643 
1644 
1645 template <int rank_, int dim, typename Number>
1646 DEAL_II_CONSTEXPR inline void
1648 {
1649  for (unsigned int i = 0; i < dim; ++i)
1651 }
1652 
1653 
1654 template <int rank_, int dim, typename Number>
1655 constexpr std::size_t
1657 {
1658  return sizeof(Tensor<rank_, dim, Number>);
1659 }
1660 
1661 
1662 template <int rank_, int dim, typename Number>
1663 template <class Archive>
1664 inline void
1665 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1666 {
1667  ar &values;
1668 }
1669 
1670 
1671 template <int rank_, int dim, typename Number>
1673 
1674 #endif // DOXYGEN
1675 
1676 /* ----------------- Non-member functions operating on tensors. ------------ */
1677 
1682 
1690 template <int rank_, int dim, typename Number>
1691 inline std::ostream &
1692 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1693 {
1694  for (unsigned int i = 0; i < dim; ++i)
1695  {
1696  out << p[i];
1697  if (i != dim - 1)
1698  out << ' ';
1699  }
1700 
1701  return out;
1702 }
1703 
1704 
1711 template <int dim, typename Number>
1712 inline std::ostream &
1713 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1714 {
1715  out << static_cast<const Number &>(p);
1716  return out;
1717 }
1718 
1719 
1721 
1725 
1726 
1737 template <int dim, typename Number, typename Other>
1740  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1741 {
1742  return object * static_cast<const Number &>(t);
1743 }
1744 
1745 
1746 
1757 template <int dim, typename Number, typename Other>
1760  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1761 {
1762  return static_cast<const Number &>(t) * object;
1763 }
1764 
1765 
1777 template <int dim, typename Number, typename OtherNumber>
1781  const Tensor<0, dim, OtherNumber> &src2)
1782 {
1783  return static_cast<const Number &>(src1) *
1784  static_cast<const OtherNumber &>(src2);
1785 }
1786 
1787 
1795 template <int dim, typename Number, typename OtherNumber>
1797  Tensor<0,
1798  dim,
1799  typename ProductType<Number,
1800  typename EnableIfScalar<OtherNumber>::type>::type>
1801  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
1802 {
1803  return static_cast<const Number &>(t) / factor;
1804 }
1805 
1806 
1814 template <int dim, typename Number, typename OtherNumber>
1818  const Tensor<0, dim, OtherNumber> &q)
1819 {
1820  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1821 }
1822 
1823 
1831 template <int dim, typename Number, typename OtherNumber>
1835  const Tensor<0, dim, OtherNumber> &q)
1836 {
1837  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1838 }
1839 
1840 
1853 template <int rank, int dim, typename Number, typename OtherNumber>
1855  Tensor<rank,
1856  dim,
1857  typename ProductType<Number,
1858  typename EnableIfScalar<OtherNumber>::type>::type>
1859  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1860 {
1861  // recurse over the base objects
1863  for (unsigned int d = 0; d < dim; ++d)
1864  tt[d] = t[d] * factor;
1865  return tt;
1866 }
1867 
1868 
1881 template <int rank, int dim, typename Number, typename OtherNumber>
1883  Tensor<rank,
1884  dim,
1886  OtherNumber>::type>
1887  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
1888 {
1889  // simply forward to the operator above
1890  return t * factor;
1891 }
1892 
1893 
1894 namespace internal
1895 {
1896  namespace TensorImplementation
1897  {
1898  template <int rank,
1899  int dim,
1900  typename Number,
1901  typename OtherNumber,
1902  typename std::enable_if<
1903  !std::is_integral<
1904  typename ProductType<Number, OtherNumber>::type>::value,
1905  int>::type = 0>
1909  const OtherNumber & factor)
1910  {
1912  const Number inverse_factor = Number(1.) / factor;
1913  // recurse over the base objects
1914  for (unsigned int d = 0; d < dim; ++d)
1915  tt[d] = t[d] * inverse_factor;
1916  return tt;
1917  }
1918 
1919 
1920  template <int rank,
1921  int dim,
1922  typename Number,
1923  typename OtherNumber,
1924  typename std::enable_if<
1925  std::is_integral<
1926  typename ProductType<Number, OtherNumber>::type>::value,
1927  int>::type = 0>
1931  const OtherNumber & factor)
1932  {
1934  // recurse over the base objects
1935  for (unsigned int d = 0; d < dim; ++d)
1936  tt[d] = t[d] / factor;
1937  return tt;
1938  }
1939  } // namespace TensorImplementation
1940 } // namespace internal
1941 
1942 
1952 template <int rank, int dim, typename Number, typename OtherNumber>
1954  Tensor<rank,
1955  dim,
1956  typename ProductType<Number,
1957  typename EnableIfScalar<OtherNumber>::type>::type>
1958  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1959 {
1961 }
1962 
1963 
1973 template <int rank, int dim, typename Number, typename OtherNumber>
1978 {
1980 
1981  for (unsigned int i = 0; i < dim; ++i)
1982  tmp[i] += q[i];
1983 
1984  return tmp;
1985 }
1986 
1987 
1997 template <int rank, int dim, typename Number, typename OtherNumber>
2002 {
2004 
2005  for (unsigned int i = 0; i < dim; ++i)
2006  tmp[i] -= q[i];
2007 
2008  return tmp;
2009 }
2010 
2017 template <int dim, typename Number, typename OtherNumber>
2021  const Tensor<0, dim, OtherNumber> &src2)
2022 {
2024 
2025  tmp *= src2;
2026 
2027  return tmp;
2028 }
2029 
2046 template <int rank, int dim, typename Number, typename OtherNumber>
2050  const Tensor<rank, dim, OtherNumber> &src2)
2051 {
2053 
2054  for (unsigned int i = 0; i < dim; ++i)
2055  tmp[i] = schur_product(Tensor<rank - 1, dim, Number>(src1[i]),
2057 
2058  return tmp;
2059 }
2060 
2062 
2066 
2067 
2090 template <int rank_1,
2091  int rank_2,
2092  int dim,
2093  typename Number,
2094  typename OtherNumber,
2095  typename = typename std::enable_if<rank_1 >= 1 && rank_2 >= 1>::type>
2097  typename Tensor<rank_1 + rank_2 - 2,
2098  dim,
2099  typename ProductType<Number, OtherNumber>::type>::tensor_type
2102 {
2103  typename Tensor<rank_1 + rank_2 - 2,
2104  dim,
2105  typename ProductType<Number, OtherNumber>::type>::tensor_type
2106  result{};
2107 
2108  TensorAccessors::internal::
2109  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2110  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2111  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2112 
2113  return result;
2114 }
2115 
2116 
2145 template <int index_1,
2146  int index_2,
2147  int rank_1,
2148  int rank_2,
2149  int dim,
2150  typename Number,
2151  typename OtherNumber>
2153  typename Tensor<rank_1 + rank_2 - 2,
2154  dim,
2155  typename ProductType<Number, OtherNumber>::type>::tensor_type
2158 {
2159  Assert(0 <= index_1 && index_1 < rank_1,
2160  ExcMessage(
2161  "The specified index_1 must lie within the range [0,rank_1)"));
2162  Assert(0 <= index_2 && index_2 < rank_2,
2163  ExcMessage(
2164  "The specified index_2 must lie within the range [0,rank_2)"));
2165 
2166  using namespace TensorAccessors;
2167  using namespace TensorAccessors::internal;
2168 
2169  // Reorder index_1 to the end of src1:
2171  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2172 
2173  // Reorder index_2 to the end of src2:
2175  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2176 
2177  typename Tensor<rank_1 + rank_2 - 2,
2178  dim,
2179  typename ProductType<Number, OtherNumber>::type>::tensor_type
2180  result{};
2181  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2182  return result;
2183 }
2184 
2185 
2216 template <int index_1,
2217  int index_2,
2218  int index_3,
2219  int index_4,
2220  int rank_1,
2221  int rank_2,
2222  int dim,
2223  typename Number,
2224  typename OtherNumber>
2225 DEAL_II_CONSTEXPR inline
2226  typename Tensor<rank_1 + rank_2 - 4,
2227  dim,
2228  typename ProductType<Number, OtherNumber>::type>::tensor_type
2231 {
2232  Assert(0 <= index_1 && index_1 < rank_1,
2233  ExcMessage(
2234  "The specified index_1 must lie within the range [0,rank_1)"));
2235  Assert(0 <= index_3 && index_3 < rank_1,
2236  ExcMessage(
2237  "The specified index_3 must lie within the range [0,rank_1)"));
2238  Assert(index_1 != index_3,
2239  ExcMessage("index_1 and index_3 must not be the same"));
2240  Assert(0 <= index_2 && index_2 < rank_2,
2241  ExcMessage(
2242  "The specified index_2 must lie within the range [0,rank_2)"));
2243  Assert(0 <= index_4 && index_4 < rank_2,
2244  ExcMessage(
2245  "The specified index_4 must lie within the range [0,rank_2)"));
2246  Assert(index_2 != index_4,
2247  ExcMessage("index_2 and index_4 must not be the same"));
2248 
2249  using namespace TensorAccessors;
2250  using namespace TensorAccessors::internal;
2251 
2252  // Reorder index_1 to the end of src1:
2254  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2255 
2256  // Reorder index_2 to the end of src2:
2258  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2259 
2260  // Now, reorder index_3 to the end of src1. We have to make sure to
2261  // preserve the original ordering: index_1 has been removed. If
2262  // index_3 > index_1, we have to use (index_3 - 1) instead:
2264  (index_3 < index_1 ? index_3 : index_3 - 1),
2265  rank_1,
2266  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2267  reord_3 =
2268  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2269  index_3 - 1,
2270  rank_1 > (reord_1);
2271 
2272  // Now, reorder index_4 to the end of src2. We have to make sure to
2273  // preserve the original ordering: index_2 has been removed. If
2274  // index_4 > index_2, we have to use (index_4 - 1) instead:
2275  ReorderedIndexView<
2276  (index_4 < index_2 ? index_4 : index_4 - 1),
2277  rank_2,
2278  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2279  reord_4 =
2280  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2281  index_4 - 1,
2282  rank_2 > (reord_2);
2283 
2284  typename Tensor<rank_1 + rank_2 - 4,
2285  dim,
2286  typename ProductType<Number, OtherNumber>::type>::tensor_type
2287  result{};
2288  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2289  return result;
2290 }
2291 
2292 
2305 template <int rank, int dim, typename Number, typename OtherNumber>
2309  const Tensor<rank, dim, OtherNumber> &right)
2310 {
2311  typename ProductType<Number, OtherNumber>::type result{};
2312  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2313  return result;
2314 }
2315 
2316 
2334 template <template <int, int, typename> class TensorT1,
2335  template <int, int, typename> class TensorT2,
2336  template <int, int, typename> class TensorT3,
2337  int rank_1,
2338  int rank_2,
2339  int dim,
2340  typename T1,
2341  typename T2,
2342  typename T3>
2345  contract3(const TensorT1<rank_1, dim, T1> & left,
2346  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2347  const TensorT3<rank_2, dim, T3> & right)
2348 {
2349  using return_type =
2351  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2352  middle,
2353  right);
2354 }
2355 
2356 
2367 template <int rank_1,
2368  int rank_2,
2369  int dim,
2370  typename Number,
2371  typename OtherNumber>
2376 {
2377  typename Tensor<rank_1 + rank_2,
2378  dim,
2379  typename ProductType<Number, OtherNumber>::type>::tensor_type
2380  result{};
2381  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2382  return result;
2383 }
2384 
2385 
2387 
2391 
2392 
2403 template <int dim, typename Number>
2406 {
2407  Assert(dim == 2, ExcInternalError());
2408 
2409  Tensor<1, dim, Number> result;
2410 
2411  result[0] = src[1];
2412  result[1] = -src[0];
2413 
2414  return result;
2415 }
2416 
2417 
2427 template <int dim, typename Number1, typename Number2>
2431  const Tensor<1, dim, Number2> &src2)
2432 {
2433  Assert(dim == 3, ExcInternalError());
2434 
2436 
2437  // avoid compiler warnings
2438  constexpr int s0 = 0 % dim;
2439  constexpr int s1 = 1 % dim;
2440  constexpr int s2 = 2 % dim;
2441 
2442  result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2443  result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2444  result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2445 
2446  return result;
2447 }
2448 
2449 
2451 
2455 
2456 
2462 template <int dim, typename Number>
2465 {
2466  // Compute the determinant using the Laplace expansion of the
2467  // determinant. We expand along the last row.
2468  Number det = internal::NumberType<Number>::value(0.0);
2469 
2470  for (unsigned int k = 0; k < dim; ++k)
2471  {
2472  Tensor<2, dim - 1, Number> minor;
2473  for (unsigned int i = 0; i < dim - 1; ++i)
2474  for (unsigned int j = 0; j < dim - 1; ++j)
2475  minor[i][j] = t[i][j < k ? j : j + 1];
2476 
2477  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2478 
2479  det += t[dim - 1][k] * cofactor;
2480  }
2481 
2482  return ((dim % 2 == 0) ? 1. : -1.) * det;
2483 }
2484 
2490 template <typename Number>
2491 constexpr DEAL_II_ALWAYS_INLINE Number
2493 {
2494  return t[0][0];
2495 }
2496 
2497 
2504 template <int dim, typename Number>
2507 {
2508  Number t = d[0][0];
2509  for (unsigned int i = 1; i < dim; ++i)
2510  t += d[i][i];
2511  return t;
2512 }
2513 
2514 
2523 template <int dim, typename Number>
2526 {
2527  Number return_tensor[dim][dim];
2528 
2529  // if desired, take over the
2530  // inversion of a 4x4 tensor
2531  // from the FullMatrix
2532  AssertThrow(false, ExcNotImplemented());
2533 
2534  return Tensor<2, dim, Number>(return_tensor);
2535 }
2536 
2537 
2538 #ifndef DOXYGEN
2539 
2540 template <typename Number>
2542  invert(const Tensor<2, 1, Number> &t)
2543 {
2544  Tensor<2, 1, Number> return_tensor;
2545 
2546  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2547 
2548  return return_tensor;
2549 }
2550 
2551 
2552 template <typename Number>
2554  invert(const Tensor<2, 2, Number> &t)
2555 {
2556  Tensor<2, 2, Number> return_tensor;
2557 
2558  // this is Maple output,
2559  // thus a bit unstructured
2560  const Number inv_det_t = internal::NumberType<Number>::value(
2561  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2562  return_tensor[0][0] = t[1][1];
2563  return_tensor[0][1] = -t[0][1];
2564  return_tensor[1][0] = -t[1][0];
2565  return_tensor[1][1] = t[0][0];
2566  return_tensor *= inv_det_t;
2567 
2568  return return_tensor;
2569 }
2570 
2571 
2572 template <typename Number>
2574  invert(const Tensor<2, 3, Number> &t)
2575 {
2576  Tensor<2, 3, Number> return_tensor;
2577 
2578  const Number t4 = internal::NumberType<Number>::value(t[0][0] * t[1][1]),
2579  t6 = internal::NumberType<Number>::value(t[0][0] * t[1][2]),
2580  t8 = internal::NumberType<Number>::value(t[0][1] * t[1][0]),
2581  t00 = internal::NumberType<Number>::value(t[0][2] * t[1][0]),
2582  t01 = internal::NumberType<Number>::value(t[0][1] * t[2][0]),
2583  t04 = internal::NumberType<Number>::value(t[0][2] * t[2][0]),
2585  1.0 / (t4 * t[2][2] - t6 * t[2][1] - t8 * t[2][2] +
2586  t00 * t[2][1] + t01 * t[1][2] - t04 * t[1][1]));
2587  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2588  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2589  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2590  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2591  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2592  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2593  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2594  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2595  return_tensor[1][1] =
2596  internal::NumberType<Number>::value(t[0][0] * t[2][2]) - t04;
2597  return_tensor[1][2] = t00 - t6;
2598  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2599  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2600  return_tensor[2][1] =
2601  t01 - internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2602  return_tensor[2][2] = internal::NumberType<Number>::value(t4 - t8);
2603  return_tensor *= inv_det_t;
2604 
2605  return return_tensor;
2606 }
2607 
2608 #endif /* DOXYGEN */
2609 
2610 
2616 template <int dim, typename Number>
2619 {
2621  for (unsigned int i = 0; i < dim; ++i)
2622  {
2623  tt[i][i] = t[i][i];
2624  for (unsigned int j = i + 1; j < dim; ++j)
2625  {
2626  tt[i][j] = t[j][i];
2627  tt[j][i] = t[i][j];
2628  };
2629  }
2630  return tt;
2631 }
2632 
2633 
2647 template <int dim, typename Number>
2648 constexpr Tensor<2, dim, Number>
2650 {
2651  return determinant(t) * invert(t);
2652 }
2653 
2654 
2668 template <int dim, typename Number>
2669 constexpr Tensor<2, dim, Number>
2671 {
2672  return transpose(adjugate(t));
2673 }
2674 
2675 
2739 template <int dim, typename Number>
2742 
2743 
2751 template <int dim, typename Number>
2752 inline Number
2754 {
2756  for (unsigned int j = 0; j < dim; ++j)
2757  {
2759  for (unsigned int i = 0; i < dim; ++i)
2760  sum += std::fabs(t[i][j]);
2761 
2762  if (sum > max)
2763  max = sum;
2764  }
2765 
2766  return max;
2767 }
2768 
2769 
2777 template <int dim, typename Number>
2778 inline Number
2780 {
2782  for (unsigned int i = 0; i < dim; ++i)
2783  {
2785  for (unsigned int j = 0; j < dim; ++j)
2786  sum += std::fabs(t[i][j]);
2787 
2788  if (sum > max)
2789  max = sum;
2790  }
2791 
2792  return max;
2793 }
2794 
2796 
2797 
2798 #ifndef DOXYGEN
2799 
2800 
2801 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2802 
2803 // Specialization of functions for ADOL-C number types when
2804 // the advanced branching feature is used
2805 template <int dim>
2806 inline adouble
2808 {
2810  for (unsigned int j = 0; j < dim; ++j)
2811  {
2813  for (unsigned int i = 0; i < dim; ++i)
2814  sum += std::fabs(t[i][j]);
2815 
2816  condassign(max, (sum > max), sum, max);
2817  }
2818 
2819  return max;
2820 }
2821 
2822 
2823 template <int dim>
2824 inline adouble
2826 {
2828  for (unsigned int i = 0; i < dim; ++i)
2829  {
2831  for (unsigned int j = 0; j < dim; ++j)
2832  sum += std::fabs(t[i][j]);
2833 
2834  condassign(max, (sum > max), sum, max);
2835  }
2836 
2837  return max;
2838 }
2839 
2840 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2841 
2842 
2843 #endif // DOXYGEN
2844 
2846 
2847 #endif
constexpr Tensor & operator+=(const Tensor< rank_, dim, OtherNumber > &)
constexpr Tensor()=default
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:2525
Number * begin_raw()
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1568
constexpr ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:2308
constexpr Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2670
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:2618
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:1740
constexpr Tensor operator-() const
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2405
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2753
#define AssertIndexRange(index, range)
Definition: exceptions.h:1636
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:2430
constexpr void clear()
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2779
constexpr Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:772
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:1521
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
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:1411
constexpr Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2506
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:363
#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:2229
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:2649
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:2345
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:1801
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > division_operator(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1908
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:1817
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
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:2020
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:2374
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:2156
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 Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2464