Reference documentation for deal.II version Git e159ac89de 2020-06-06 19:38:41 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tensor.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 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 
29 
30 #ifdef DEAL_II_WITH_ADOLC
31 # include <adolc/adouble.h> // Taped double
32 #endif
33 
34 #include <cmath>
35 #include <ostream>
36 #include <utility>
37 #include <vector>
38 
39 
41 
42 // Forward declarations:
43 #ifndef DOXYGEN
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 
520  template <typename OtherNumber>
521  constexpr DEAL_II_CUDA_HOST_DEV
522  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
523 
527  template <typename OtherNumber>
528  constexpr Tensor(
529  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
530 
534  template <typename OtherNumber>
535  constexpr
536  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
537 
544  operator[](const unsigned int i);
545 
551  constexpr DEAL_II_CUDA_HOST_DEV const value_type &
552  operator[](const unsigned int i) const;
553 
557  DEAL_II_CONSTEXPR const Number &
558  operator[](const TableIndices<rank_> &indices) const;
559 
563  DEAL_II_CONSTEXPR Number &operator[](const TableIndices<rank_> &indices);
564 
568  Number *
569  begin_raw();
570 
574  const Number *
575  begin_raw() const;
576 
580  Number *
581  end_raw();
582 
586  const Number *
587  end_raw() const;
588 
596  template <typename OtherNumber>
598  operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
599 
606  DEAL_II_CONSTEXPR Tensor &
607  operator=(const Number &d);
608 
612  template <typename OtherNumber>
613  DEAL_II_CONSTEXPR bool
615 
619  template <typename OtherNumber>
620  constexpr bool
622 
628  template <typename OtherNumber>
630  operator+=(const Tensor<rank_, dim, OtherNumber> &);
631 
637  template <typename OtherNumber>
639  operator-=(const Tensor<rank_, dim, OtherNumber> &);
640 
647  template <typename OtherNumber>
649  operator*=(const OtherNumber &factor);
650 
656  template <typename OtherNumber>
658  operator/=(const OtherNumber &factor);
659 
666  operator-() const;
667 
680  DEAL_II_CONSTEXPR void
681  clear();
682 
692  norm() const;
693 
702  norm_square() const;
703 
711  template <typename OtherNumber>
712  void
713  unroll(Vector<OtherNumber> &result) const;
714 
719  static DEAL_II_CONSTEXPR unsigned int
720  component_to_unrolled_index(const TableIndices<rank_> &indices);
721 
728  unrolled_to_component_indices(const unsigned int i);
729 
734  static constexpr std::size_t
736 
741  template <class Archive>
742  void
743  serialize(Archive &ar, const unsigned int version);
744 
750 
751 private:
755  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
756  // ... avoid a compiler warning in case of dim == 0 and ensure that the
757  // array always has positive size.
758 
762  template <typename OtherNumber>
763  void
764  unroll_recursion(Vector<OtherNumber> &result,
765  unsigned int & start_index) const;
766 
773  template <typename ArrayLike, std::size_t... Indices>
774  constexpr DEAL_II_CUDA_HOST_DEV
775  Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
776 
777  // Allow an arbitrary Tensor to access the underlying values.
778  template <int, int, typename>
779  friend class Tensor;
780 
781  // Point is allowed access to the coordinates. This is supposed to improve
782  // speed.
783  friend class Point<dim, Number>;
784 };
785 
786 
787 #ifndef DOXYGEN
788 namespace internal
789 {
790  // Workaround: The following 4 overloads are necessary to be able to
791  // compile the library with Apple Clang 8 and older. We should remove
792  // these overloads again when we bump the minimal required version to
793  // something later than clang-3.6 / Apple Clang 6.3.
794  template <int rank, int dim, typename T, typename U>
795  struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
796  {
797  using type =
799  };
800 
801  template <int rank, int dim, typename T, typename U>
802  struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
803  {
804  using type =
806  };
807 
808  template <typename T, int rank, int dim, typename U>
809  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
810  {
811  using type =
813  };
814 
815  template <int rank, int dim, typename T, typename U>
816  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
817  {
818  using type =
820  };
821  // end workaround
822 
827  template <int rank, int dim, typename T>
828  struct NumberType<Tensor<rank, dim, T>>
829  {
830  static constexpr DEAL_II_ALWAYS_INLINE const Tensor<rank, dim, T> &
831  value(const Tensor<rank, dim, T> &t)
832  {
833  return t;
834  }
835 
837  value(const T &t)
838  {
840  tmp = t;
841  return tmp;
842  }
843  };
844 } // namespace internal
845 
846 
847 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
848 
849 
850 template <int dim, typename Number>
853  // Some auto-differentiable numbers need explicit
854  // zero initialization such as adtl::adouble.
855  : Tensor{0.0}
856 {}
857 
858 
859 
860 template <int dim, typename Number>
861 template <typename OtherNumber>
863  Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
864  : value(internal::NumberType<Number>::value(initializer))
865 {}
866 
867 
868 
869 template <int dim, typename Number>
870 template <typename OtherNumber>
873  : Tensor{p.value}
874 {}
875 
876 
877 
878 template <int dim, typename Number>
879 inline Number *
881 {
882  return std::addressof(value);
883 }
884 
885 
886 
887 template <int dim, typename Number>
888 inline const Number *
890 {
891  return std::addressof(value);
892 }
893 
894 
895 
896 template <int dim, typename Number>
897 inline Number *
899 {
901 }
902 
903 
904 
905 template <int dim, typename Number>
906 const Number *
908 {
910 }
911 
912 
913 
914 template <int dim, typename Number>
917 {
918  // We cannot use Assert inside a CUDA kernel
919 # ifndef __CUDA_ARCH__
920  Assert(dim != 0,
921  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
922 # endif
923  return value;
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>
941 template <typename OtherNumber>
945 {
947  return *this;
948 }
949 
950 
951 # ifdef __INTEL_COMPILER
952 template <int dim, typename Number>
956 {
957  value = p.value;
958  return *this;
959 }
960 # endif
961 
962 
963 template <int dim, typename Number>
964 template <typename OtherNumber>
967  Tensor<0, dim, Number>::operator=(const OtherNumber &d)
968 {
970  return *this;
971 }
972 
973 
974 template <int dim, typename Number>
975 template <typename OtherNumber>
976 DEAL_II_CONSTEXPR inline bool
978 {
979 # if defined(DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING)
982  ExcMessage(
983  "The Tensor equality operator for ADOL-C taped numbers has not yet "
984  "been extended to support advanced branching."));
985 # endif
986 
987  return numbers::values_are_equal(value, p.value);
988 }
989 
990 
991 template <int dim, typename Number>
992 template <typename OtherNumber>
993 constexpr bool
995 {
996  return !((*this) == p);
997 }
998 
999 
1000 template <int dim, typename Number>
1001 template <typename OtherNumber>
1005 {
1006  value += p.value;
1007  return *this;
1008 }
1009 
1010 
1011 template <int dim, typename Number>
1012 template <typename OtherNumber>
1016 {
1017  value -= p.value;
1018  return *this;
1019 }
1020 
1021 
1022 
1023 namespace internal
1024 {
1025  namespace ComplexWorkaround
1026  {
1027  template <typename Number, typename OtherNumber>
1029  multiply_assign_scalar(Number &val, const OtherNumber &s)
1030  {
1031  val *= s;
1032  }
1033 
1034 # ifdef __CUDA_ARCH__
1035  template <typename Number, typename OtherNumber>
1037  multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
1038  {
1039  printf("This function is not implemented for std::complex<Number>!\n");
1040  assert(false);
1041  }
1042 # endif
1043  } // namespace ComplexWorkaround
1044 } // namespace internal
1045 
1046 
1047 template <int dim, typename Number>
1048 template <typename OtherNumber>
1051  Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1052 {
1053  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1054  return *this;
1055 }
1056 
1057 
1058 
1059 template <int dim, typename Number>
1060 template <typename OtherNumber>
1062 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1063 {
1064  value /= s;
1065  return *this;
1066 }
1067 
1068 
1069 template <int dim, typename Number>
1072 {
1073  return -value;
1074 }
1075 
1076 
1077 template <int dim, typename Number>
1080 {
1081  Assert(dim != 0,
1082  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1084 }
1085 
1086 
1087 template <int dim, typename Number>
1091 {
1092  // We cannot use Assert inside a CUDA kernel
1093 # ifndef __CUDA_ARCH__
1094  Assert(dim != 0,
1095  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1096 # endif
1098 }
1099 
1100 
1101 template <int dim, typename Number>
1102 template <typename OtherNumber>
1103 inline void
1104 Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1105  unsigned int & index) const
1106 {
1107  Assert(dim != 0,
1108  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1109  result[index] = value;
1110  ++index;
1111 }
1112 
1113 
1114 template <int dim, typename Number>
1115 DEAL_II_CONSTEXPR inline void
1117 {
1118  // Some auto-differentiable numbers need explicit
1119  // zero initialization.
1121 }
1122 
1123 
1124 template <int dim, typename Number>
1125 template <class Archive>
1126 inline void
1127 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1128 {
1129  ar &value;
1130 }
1131 
1132 
1133 template <int dim, typename Number>
1135 
1136 
1137 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1138 
1139 template <int rank_, int dim, typename Number>
1140 template <typename ArrayLike, std::size_t... indices>
1142  Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1143  std::index_sequence<indices...>)
1144  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
1145 {
1146  static_assert(sizeof...(indices) == dim,
1147  "dim should match the number of indices");
1148 }
1149 
1150 
1151 template <int rank_, int dim, typename Number>
1154  : Tensor(initializer, std::make_index_sequence<dim>{})
1155 {}
1156 
1157 
1158 template <int rank_, int dim, typename Number>
1159 template <typename OtherNumber>
1162  const Tensor<rank_, dim, OtherNumber> &initializer)
1163  : Tensor(initializer, std::make_index_sequence<dim>{})
1164 {}
1165 
1166 
1167 template <int rank_, int dim, typename Number>
1168 template <typename OtherNumber>
1169 constexpr DEAL_II_ALWAYS_INLINE
1171  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1172  : Tensor(initializer, std::make_index_sequence<dim>{})
1173 {}
1174 
1175 
1176 template <int rank_, int dim, typename Number>
1177 template <typename OtherNumber>
1179  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1180 {
1181  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
1182 }
1183 
1184 
1185 
1186 namespace internal
1187 {
1188  namespace TensorSubscriptor
1189  {
1190  template <typename ArrayElementType, int dim>
1192  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1193  subscript(ArrayElementType * values,
1194  const unsigned int i,
1195  std::integral_constant<int, dim>)
1196  {
1197  // We cannot use Assert in a CUDA kernel
1198 # ifndef __CUDA_ARCH__
1199  AssertIndexRange(i, dim);
1200 # endif
1201  return values[i];
1202  }
1203 
1204  // The variables within this struct will be referenced in the next function.
1205  // It is a workaround that allows returning a reference to a static variable
1206  // while allowing constexpr evaluation of the function.
1207  // It has to be defined outside the function because constexpr functions
1208  // cannot define static variables
1209  template <typename ArrayElementType>
1210  struct Uninitialized
1211  {
1212  static ArrayElementType value;
1213  };
1214 
1215  template <typename Type>
1217 
1218  template <typename ArrayElementType>
1220  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1221  subscript(ArrayElementType *,
1222  const unsigned int,
1223  std::integral_constant<int, 0>)
1224  {
1225  // We cannot use Assert in a CUDA kernel
1226 # ifndef __CUDA_ARCH__
1227  Assert(
1228  false,
1229  ExcMessage(
1230  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1231 # endif
1233  }
1234  } // namespace TensorSubscriptor
1235 } // namespace internal
1236 
1237 
1238 template <int rank_, int dim, typename Number>
1241  operator[](const unsigned int i)
1242 {
1243  return ::internal::TensorSubscriptor::subscript(
1244  values, i, std::integral_constant<int, dim>());
1245 }
1246 
1247 
1248 template <int rank_, int dim, typename Number>
1249 constexpr DEAL_II_ALWAYS_INLINE
1251  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1252 {
1253  AssertIndexRange(i, dim);
1254 
1255  return values[i];
1256 }
1257 
1258 
1259 template <int rank_, int dim, typename Number>
1260 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1262  operator[](const TableIndices<rank_> &indices) const
1263 {
1264  Assert(dim != 0,
1265  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1266 
1267  return TensorAccessors::extract<rank_>(*this, indices);
1268 }
1269 
1270 
1271 
1272 template <int rank_, int dim, typename Number>
1275 {
1276  Assert(dim != 0,
1277  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1278 
1279  return TensorAccessors::extract<rank_>(*this, indices);
1280 }
1281 
1282 
1283 
1284 template <int rank_, int dim, typename Number>
1285 inline Number *
1287 {
1288  return std::addressof(
1289  this->operator[](this->unrolled_to_component_indices(0)));
1290 }
1291 
1292 
1293 
1294 template <int rank_, int dim, typename Number>
1295 inline const Number *
1297 {
1298  return std::addressof(
1299  this->operator[](this->unrolled_to_component_indices(0)));
1300 }
1301 
1302 
1303 
1304 template <int rank_, int dim, typename Number>
1305 inline Number *
1307 {
1309 }
1310 
1311 
1312 
1313 template <int rank_, int dim, typename Number>
1314 inline const Number *
1316 {
1318 }
1319 
1320 
1321 
1322 template <int rank_, int dim, typename Number>
1323 template <typename OtherNumber>
1326 {
1327  // The following loop could be written more concisely using std::copy, but
1328  // that function is only constexpr from C++20 on.
1329  for (unsigned int i = 0; i < dim; ++i)
1330  values[i] = t.values[i];
1331  return *this;
1332 }
1333 
1334 
1335 template <int rank_, int dim, typename Number>
1338 {
1340  ExcMessage("Only assignment with zero is allowed"));
1341  (void)d;
1342 
1343  for (unsigned int i = 0; i < dim; ++i)
1345  return *this;
1346 }
1347 
1348 
1349 template <int rank_, int dim, typename Number>
1350 template <typename OtherNumber>
1351 DEAL_II_CONSTEXPR inline bool
1354 {
1355  for (unsigned int i = 0; i < dim; ++i)
1356  if (values[i] != p.values[i])
1357  return false;
1358  return true;
1359 }
1360 
1361 
1362 // At some places in the library, we have Point<0> for formal reasons
1363 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1364 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1365 // in the above function that the loop end check always fails, we
1366 // implement this function here
1367 template <>
1368 template <>
1369 DEAL_II_CONSTEXPR inline bool
1371 {
1372  return true;
1373 }
1374 
1375 
1376 template <int rank_, int dim, typename Number>
1377 template <typename OtherNumber>
1378 constexpr bool
1381 {
1382  return !((*this) == p);
1383 }
1384 
1385 
1386 template <int rank_, int dim, typename Number>
1387 template <typename OtherNumber>
1392 {
1393  for (unsigned int i = 0; i < dim; ++i)
1394  values[i] += p.values[i];
1395  return *this;
1396 }
1397 
1398 
1399 template <int rank_, int dim, typename Number>
1400 template <typename OtherNumber>
1405 {
1406  for (unsigned int i = 0; i < dim; ++i)
1407  values[i] -= p.values[i];
1408  return *this;
1409 }
1410 
1411 
1412 template <int rank_, int dim, typename Number>
1413 template <typename OtherNumber>
1416  Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1417 {
1418  for (unsigned int i = 0; i < dim; ++i)
1419  values[i] *= s;
1420  return *this;
1421 }
1422 
1423 
1424 namespace internal
1425 {
1426  namespace TensorImplementation
1427  {
1428  template <int rank,
1429  int dim,
1430  typename Number,
1431  typename OtherNumber,
1432  typename std::enable_if<
1433  !std::is_integral<
1436  int>::type = 0>
1439  const OtherNumber &factor)
1440  {
1441  const Number inverse_factor = Number(1.) / factor;
1442  // recurse over the base objects
1443  for (unsigned int d = 0; d < dim; ++d)
1444  t[d] *= inverse_factor;
1445  }
1446 
1447 
1448  template <int rank,
1449  int dim,
1450  typename Number,
1451  typename OtherNumber,
1452  typename std::enable_if<
1453  std::is_integral<
1456  int>::type = 0>
1459  const OtherNumber &factor)
1460  {
1461  // recurse over the base objects
1462  for (unsigned int d = 0; d < dim; ++d)
1463  t[d] /= factor;
1464  }
1465  } // namespace TensorImplementation
1466 } // namespace internal
1467 
1468 
1469 template <int rank_, int dim, typename Number>
1470 template <typename OtherNumber>
1473  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1474 {
1476  return *this;
1477 }
1478 
1479 
1480 template <int rank_, int dim, typename Number>
1484 {
1486 
1487  for (unsigned int i = 0; i < dim; ++i)
1488  tmp.values[i] = -values[i];
1489 
1490  return tmp;
1491 }
1492 
1493 
1494 template <int rank_, int dim, typename Number>
1497 {
1498  return std::sqrt(norm_square());
1499 }
1500 
1501 
1502 template <int rank_, int dim, typename Number>
1506 {
1509  for (unsigned int i = 0; i < dim; ++i)
1510  s += values[i].norm_square();
1511 
1512  return s;
1513 }
1514 
1515 
1516 template <int rank_, int dim, typename Number>
1517 template <typename OtherNumber>
1518 inline void
1519 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
1520 {
1521  AssertDimension(result.size(),
1522  (Utilities::fixed_power<rank_, unsigned int>(dim)));
1523 
1524  unsigned int index = 0;
1525  unroll_recursion(result, index);
1526 }
1527 
1528 
1529 template <int rank_, int dim, typename Number>
1530 template <typename OtherNumber>
1531 inline void
1532 Tensor<rank_, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1533  unsigned int & index) const
1534 {
1535  for (unsigned int i = 0; i < dim; ++i)
1536  values[i].unroll_recursion(result, index);
1537 }
1538 
1539 
1540 template <int rank_, int dim, typename Number>
1541 DEAL_II_CONSTEXPR inline unsigned int
1543  const TableIndices<rank_> &indices)
1544 {
1545  unsigned int index = 0;
1546  for (int r = 0; r < rank_; ++r)
1547  index = index * dim + indices[r];
1548 
1549  return index;
1550 }
1551 
1552 
1553 
1554 namespace internal
1555 {
1556  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1557  // and rank=2. Make sure we don't have compiler warnings.
1558 
1559  template <int dim>
1560  inline DEAL_II_CONSTEXPR unsigned int
1561  mod(const unsigned int x)
1562  {
1563  return x % dim;
1564  }
1565 
1566  template <>
1567  inline unsigned int
1568  mod<0>(const unsigned int x)
1569  {
1570  Assert(false, ExcInternalError());
1571  return x;
1572  }
1573 
1574  template <int dim>
1575  inline DEAL_II_CONSTEXPR unsigned int
1576  div(const unsigned int x)
1577  {
1578  return x / dim;
1579  }
1580 
1581  template <>
1582  inline unsigned int
1583  div<0>(const unsigned int x)
1584  {
1585  Assert(false, ExcInternalError());
1586  return x;
1587  }
1588 
1589 } // namespace internal
1590 
1591 
1592 
1593 template <int rank_, int dim, typename Number>
1596 {
1597  AssertIndexRange(i, n_independent_components);
1598 
1599  TableIndices<rank_> indices;
1600 
1601  unsigned int remainder = i;
1602  for (int r = rank_ - 1; r >= 0; --r)
1603  {
1604  indices[r] = internal::mod<dim>(remainder);
1605  remainder = internal::div<dim>(remainder);
1606  }
1607  Assert(remainder == 0, ExcInternalError());
1608 
1609  return indices;
1610 }
1611 
1612 
1613 template <int rank_, int dim, typename Number>
1614 DEAL_II_CONSTEXPR inline void
1616 {
1617  for (unsigned int i = 0; i < dim; ++i)
1619 }
1620 
1621 
1622 template <int rank_, int dim, typename Number>
1623 constexpr std::size_t
1625 {
1626  return sizeof(Tensor<rank_, dim, Number>);
1627 }
1628 
1629 
1630 template <int rank_, int dim, typename Number>
1631 template <class Archive>
1632 inline void
1633 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1634 {
1635  ar &values;
1636 }
1637 
1638 
1639 template <int rank_, int dim, typename Number>
1641 
1642 #endif // DOXYGEN
1643 
1644 /* ----------------- Non-member functions operating on tensors. ------------ */
1645 
1650 
1658 template <int rank_, int dim, typename Number>
1659 inline std::ostream &
1660 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1661 {
1662  for (unsigned int i = 0; i < dim; ++i)
1663  {
1664  out << p[i];
1665  if (i != dim - 1)
1666  out << ' ';
1667  }
1668 
1669  return out;
1670 }
1671 
1672 
1679 template <int dim, typename Number>
1680 inline std::ostream &
1681 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1682 {
1683  out << static_cast<const Number &>(p);
1684  return out;
1685 }
1686 
1687 
1689 
1693 
1694 
1705 template <int dim, typename Number, typename Other>
1708  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1709 {
1710  return object * static_cast<const Number &>(t);
1711 }
1712 
1713 
1714 
1725 template <int dim, typename Number, typename Other>
1728  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1729 {
1730  return static_cast<const Number &>(t) * object;
1731 }
1732 
1733 
1745 template <int dim, typename Number, typename OtherNumber>
1749  const Tensor<0, dim, OtherNumber> &src2)
1750 {
1751  return static_cast<const Number &>(src1) *
1752  static_cast<const OtherNumber &>(src2);
1753 }
1754 
1755 
1763 template <int dim, typename Number, typename OtherNumber>
1765  Tensor<0,
1766  dim,
1767  typename ProductType<Number,
1768  typename EnableIfScalar<OtherNumber>::type>::type>
1769  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
1770 {
1771  return static_cast<const Number &>(t) / factor;
1772 }
1773 
1774 
1782 template <int dim, typename Number, typename OtherNumber>
1786  const Tensor<0, dim, OtherNumber> &q)
1787 {
1788  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1789 }
1790 
1791 
1799 template <int dim, typename Number, typename OtherNumber>
1803  const Tensor<0, dim, OtherNumber> &q)
1804 {
1805  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1806 }
1807 
1808 
1821 template <int rank, int dim, typename Number, typename OtherNumber>
1823  Tensor<rank,
1824  dim,
1825  typename ProductType<Number,
1826  typename EnableIfScalar<OtherNumber>::type>::type>
1827  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1828 {
1829  // recurse over the base objects
1831  for (unsigned int d = 0; d < dim; ++d)
1832  tt[d] = t[d] * factor;
1833  return tt;
1834 }
1835 
1836 
1849 template <int rank, int dim, typename Number, typename OtherNumber>
1851  Tensor<rank,
1852  dim,
1854  OtherNumber>::type>
1855  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
1856 {
1857  // simply forward to the operator above
1858  return t * factor;
1859 }
1860 
1861 
1862 namespace internal
1863 {
1864  namespace TensorImplementation
1865  {
1866  template <int rank,
1867  int dim,
1868  typename Number,
1869  typename OtherNumber,
1870  typename std::enable_if<
1871  !std::is_integral<
1873  int>::type = 0>
1877  const OtherNumber & factor)
1878  {
1880  const Number inverse_factor = Number(1.) / factor;
1881  // recurse over the base objects
1882  for (unsigned int d = 0; d < dim; ++d)
1883  tt[d] = t[d] * inverse_factor;
1884  return tt;
1885  }
1886 
1887 
1888  template <int rank,
1889  int dim,
1890  typename Number,
1891  typename OtherNumber,
1892  typename std::enable_if<
1893  std::is_integral<
1895  int>::type = 0>
1899  const OtherNumber & factor)
1900  {
1902  // recurse over the base objects
1903  for (unsigned int d = 0; d < dim; ++d)
1904  tt[d] = t[d] / factor;
1905  return tt;
1906  }
1907  } // namespace TensorImplementation
1908 } // namespace internal
1909 
1910 
1920 template <int rank, int dim, typename Number, typename OtherNumber>
1922  Tensor<rank,
1923  dim,
1924  typename ProductType<Number,
1925  typename EnableIfScalar<OtherNumber>::type>::type>
1926  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1927 {
1929 }
1930 
1931 
1941 template <int rank, int dim, typename Number, typename OtherNumber>
1946 {
1948 
1949  for (unsigned int i = 0; i < dim; ++i)
1950  tmp[i] += q[i];
1951 
1952  return tmp;
1953 }
1954 
1955 
1965 template <int rank, int dim, typename Number, typename OtherNumber>
1970 {
1972 
1973  for (unsigned int i = 0; i < dim; ++i)
1974  tmp[i] -= q[i];
1975 
1976  return tmp;
1977 }
1978 
1985 template <int dim, typename Number, typename OtherNumber>
1989  const Tensor<0, dim, OtherNumber> &src2)
1990 {
1992 
1993  tmp *= src2;
1994 
1995  return tmp;
1996 }
1997 
2014 template <int rank, int dim, typename Number, typename OtherNumber>
2018  const Tensor<rank, dim, OtherNumber> &src2)
2019 {
2021 
2022  for (unsigned int i = 0; i < dim; ++i)
2023  tmp[i] = schur_product(Tensor<rank - 1, dim, Number>(src1[i]),
2025 
2026  return tmp;
2027 }
2028 
2030 
2034 
2035 
2058 template <int rank_1,
2059  int rank_2,
2060  int dim,
2061  typename Number,
2062  typename OtherNumber,
2063  typename = typename std::enable_if<rank_1 >= 1 && rank_2 >= 1>::type>
2065  typename Tensor<rank_1 + rank_2 - 2,
2066  dim,
2067  typename ProductType<Number, OtherNumber>::type>::tensor_type
2070 {
2071  typename Tensor<rank_1 + rank_2 - 2,
2072  dim,
2073  typename ProductType<Number, OtherNumber>::type>::tensor_type
2074  result{};
2075 
2076  TensorAccessors::internal::
2077  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2078  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2079  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2080 
2081  return result;
2082 }
2083 
2084 
2113 template <int index_1,
2114  int index_2,
2115  int rank_1,
2116  int rank_2,
2117  int dim,
2118  typename Number,
2119  typename OtherNumber>
2121  typename Tensor<rank_1 + rank_2 - 2,
2122  dim,
2123  typename ProductType<Number, OtherNumber>::type>::tensor_type
2126 {
2127  Assert(0 <= index_1 && index_1 < rank_1,
2128  ExcMessage(
2129  "The specified index_1 must lie within the range [0,rank_1)"));
2130  Assert(0 <= index_2 && index_2 < rank_2,
2131  ExcMessage(
2132  "The specified index_2 must lie within the range [0,rank_2)"));
2133 
2134  using namespace TensorAccessors;
2135  using namespace TensorAccessors::internal;
2136 
2137  // Reorder index_1 to the end of src1:
2139  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2140 
2141  // Reorder index_2 to the end of src2:
2143  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2144 
2145  typename Tensor<rank_1 + rank_2 - 2,
2146  dim,
2147  typename ProductType<Number, OtherNumber>::type>::tensor_type
2148  result{};
2149  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2150  return result;
2151 }
2152 
2153 
2184 template <int index_1,
2185  int index_2,
2186  int index_3,
2187  int index_4,
2188  int rank_1,
2189  int rank_2,
2190  int dim,
2191  typename Number,
2192  typename OtherNumber>
2193 DEAL_II_CONSTEXPR inline
2194  typename Tensor<rank_1 + rank_2 - 4,
2195  dim,
2196  typename ProductType<Number, OtherNumber>::type>::tensor_type
2199 {
2200  Assert(0 <= index_1 && index_1 < rank_1,
2201  ExcMessage(
2202  "The specified index_1 must lie within the range [0,rank_1)"));
2203  Assert(0 <= index_3 && index_3 < rank_1,
2204  ExcMessage(
2205  "The specified index_3 must lie within the range [0,rank_1)"));
2206  Assert(index_1 != index_3,
2207  ExcMessage("index_1 and index_3 must not be the same"));
2208  Assert(0 <= index_2 && index_2 < rank_2,
2209  ExcMessage(
2210  "The specified index_2 must lie within the range [0,rank_2)"));
2211  Assert(0 <= index_4 && index_4 < rank_2,
2212  ExcMessage(
2213  "The specified index_4 must lie within the range [0,rank_2)"));
2214  Assert(index_2 != index_4,
2215  ExcMessage("index_2 and index_4 must not be the same"));
2216 
2217  using namespace TensorAccessors;
2218  using namespace TensorAccessors::internal;
2219 
2220  // Reorder index_1 to the end of src1:
2222  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2223 
2224  // Reorder index_2 to the end of src2:
2226  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2227 
2228  // Now, reorder index_3 to the end of src1. We have to make sure to
2229  // preserve the original ordering: index_1 has been removed. If
2230  // index_3 > index_1, we have to use (index_3 - 1) instead:
2232  (index_3 < index_1 ? index_3 : index_3 - 1),
2233  rank_1,
2234  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2235  reord_3 =
2236  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2237  index_3 - 1,
2238  rank_1 > (reord_1);
2239 
2240  // Now, reorder index_4 to the end of src2. We have to make sure to
2241  // preserve the original ordering: index_2 has been removed. If
2242  // index_4 > index_2, we have to use (index_4 - 1) instead:
2243  ReorderedIndexView<
2244  (index_4 < index_2 ? index_4 : index_4 - 1),
2245  rank_2,
2246  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2247  reord_4 =
2248  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2249  index_4 - 1,
2250  rank_2 > (reord_2);
2251 
2252  typename Tensor<rank_1 + rank_2 - 4,
2253  dim,
2254  typename ProductType<Number, OtherNumber>::type>::tensor_type
2255  result{};
2256  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2257  return result;
2258 }
2259 
2260 
2273 template <int rank, int dim, typename Number, typename OtherNumber>
2277  const Tensor<rank, dim, OtherNumber> &right)
2278 {
2279  typename ProductType<Number, OtherNumber>::type result{};
2280  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2281  return result;
2282 }
2283 
2284 
2302 template <template <int, int, typename> class TensorT1,
2303  template <int, int, typename> class TensorT2,
2304  template <int, int, typename> class TensorT3,
2305  int rank_1,
2306  int rank_2,
2307  int dim,
2308  typename T1,
2309  typename T2,
2310  typename T3>
2313  contract3(const TensorT1<rank_1, dim, T1> & left,
2314  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2315  const TensorT3<rank_2, dim, T3> & right)
2316 {
2317  using return_type =
2319  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2320  middle,
2321  right);
2322 }
2323 
2324 
2335 template <int rank_1,
2336  int rank_2,
2337  int dim,
2338  typename Number,
2339  typename OtherNumber>
2344 {
2345  typename Tensor<rank_1 + rank_2,
2346  dim,
2347  typename ProductType<Number, OtherNumber>::type>::tensor_type
2348  result{};
2349  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2350  return result;
2351 }
2352 
2353 
2355 
2359 
2360 
2371 template <int dim, typename Number>
2374 {
2375  Assert(dim == 2, ExcInternalError());
2376 
2377  Tensor<1, dim, Number> result;
2378 
2379  result[0] = src[1];
2380  result[1] = -src[0];
2381 
2382  return result;
2383 }
2384 
2385 
2395 template <int dim, typename Number1, typename Number2>
2399  const Tensor<1, dim, Number2> &src2)
2400 {
2401  Assert(dim == 3, ExcInternalError());
2402 
2404 
2405  // avoid compiler warnings
2406  constexpr int s0 = 0 % dim;
2407  constexpr int s1 = 1 % dim;
2408  constexpr int s2 = 2 % dim;
2409 
2410  result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2411  result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2412  result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2413 
2414  return result;
2415 }
2416 
2417 
2419 
2423 
2424 
2430 template <int dim, typename Number>
2433 {
2434  // Compute the determinant using the Laplace expansion of the
2435  // determinant. We expand along the last row.
2436  Number det = internal::NumberType<Number>::value(0.0);
2437 
2438  for (unsigned int k = 0; k < dim; ++k)
2439  {
2440  Tensor<2, dim - 1, Number> minor;
2441  for (unsigned int i = 0; i < dim - 1; ++i)
2442  for (unsigned int j = 0; j < dim - 1; ++j)
2443  minor[i][j] = t[i][j < k ? j : j + 1];
2444 
2445  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2446 
2447  det += t[dim - 1][k] * cofactor;
2448  }
2449 
2450  return ((dim % 2 == 0) ? 1. : -1.) * det;
2451 }
2452 
2458 template <typename Number>
2459 constexpr DEAL_II_ALWAYS_INLINE Number
2461 {
2462  return t[0][0];
2463 }
2464 
2465 
2472 template <int dim, typename Number>
2475 {
2476  Number t = d[0][0];
2477  for (unsigned int i = 1; i < dim; ++i)
2478  t += d[i][i];
2479  return t;
2480 }
2481 
2482 
2491 template <int dim, typename Number>
2494 {
2495  Number return_tensor[dim][dim];
2496 
2497  // if desired, take over the
2498  // inversion of a 4x4 tensor
2499  // from the FullMatrix
2500  AssertThrow(false, ExcNotImplemented());
2501 
2502  return Tensor<2, dim, Number>(return_tensor);
2503 }
2504 
2505 
2506 #ifndef DOXYGEN
2507 
2508 template <typename Number>
2510  invert(const Tensor<2, 1, Number> &t)
2511 {
2512  Tensor<2, 1, Number> return_tensor;
2513 
2514  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2515 
2516  return return_tensor;
2517 }
2518 
2519 
2520 template <typename Number>
2522  invert(const Tensor<2, 2, Number> &t)
2523 {
2524  Tensor<2, 2, Number> return_tensor;
2525 
2526  // this is Maple output,
2527  // thus a bit unstructured
2528  const Number inv_det_t = internal::NumberType<Number>::value(
2529  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2530  return_tensor[0][0] = t[1][1];
2531  return_tensor[0][1] = -t[0][1];
2532  return_tensor[1][0] = -t[1][0];
2533  return_tensor[1][1] = t[0][0];
2534  return_tensor *= inv_det_t;
2535 
2536  return return_tensor;
2537 }
2538 
2539 
2540 template <typename Number>
2542  invert(const Tensor<2, 3, Number> &t)
2543 {
2544  Tensor<2, 3, Number> return_tensor;
2545 
2546  const Number t4 = internal::NumberType<Number>::value(t[0][0] * t[1][1]),
2547  t6 = internal::NumberType<Number>::value(t[0][0] * t[1][2]),
2548  t8 = internal::NumberType<Number>::value(t[0][1] * t[1][0]),
2549  t00 = internal::NumberType<Number>::value(t[0][2] * t[1][0]),
2550  t01 = internal::NumberType<Number>::value(t[0][1] * t[2][0]),
2551  t04 = internal::NumberType<Number>::value(t[0][2] * t[2][0]),
2553  1.0 / (t4 * t[2][2] - t6 * t[2][1] - t8 * t[2][2] +
2554  t00 * t[2][1] + t01 * t[1][2] - t04 * t[1][1]));
2555  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2556  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2557  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2558  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2559  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2560  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2561  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2562  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2563  return_tensor[1][1] =
2564  internal::NumberType<Number>::value(t[0][0] * t[2][2]) - t04;
2565  return_tensor[1][2] = t00 - t6;
2566  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2567  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2568  return_tensor[2][1] =
2569  t01 - internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2570  return_tensor[2][2] = internal::NumberType<Number>::value(t4 - t8);
2571  return_tensor *= inv_det_t;
2572 
2573  return return_tensor;
2574 }
2575 
2576 #endif /* DOXYGEN */
2577 
2578 
2584 template <int dim, typename Number>
2587 {
2589  for (unsigned int i = 0; i < dim; ++i)
2590  {
2591  tt[i][i] = t[i][i];
2592  for (unsigned int j = i + 1; j < dim; ++j)
2593  {
2594  tt[i][j] = t[j][i];
2595  tt[j][i] = t[i][j];
2596  };
2597  }
2598  return tt;
2599 }
2600 
2601 
2615 template <int dim, typename Number>
2616 constexpr Tensor<2, dim, Number>
2618 {
2619  return determinant(t) * invert(t);
2620 }
2621 
2622 
2636 template <int dim, typename Number>
2637 constexpr Tensor<2, dim, Number>
2639 {
2640  return transpose(adjugate(t));
2641 }
2642 
2643 
2654 template <int dim, typename Number>
2657 {
2658  Tensor<2, dim, Number> output_tensor;
2660  LAPACKFullMatrix<Number> lapack_matrix(dim);
2661  LAPACKFullMatrix<Number> result(dim);
2662 
2663  // todo: find or add dealii functionality to copy in one step.
2664  matrix.copy_from(tensor);
2665  lapack_matrix.copy_from(matrix);
2666 
2667  // now compute the svd of the matrices
2668  lapack_matrix.compute_svd();
2669 
2670  // Use the SVD results to orthogonalize: @f$U V^T@f$
2671  lapack_matrix.get_svd_u().mmult(result, lapack_matrix.get_svd_vt());
2672 
2673  // todo: find or add dealii functionality to copy in one step.
2674  matrix = result;
2675  matrix.copy_to(output_tensor);
2676  return output_tensor;
2677 }
2678 
2679 
2687 template <int dim, typename Number>
2688 inline Number
2690 {
2692  for (unsigned int j = 0; j < dim; ++j)
2693  {
2695  for (unsigned int i = 0; i < dim; ++i)
2696  sum += std::fabs(t[i][j]);
2697 
2698  if (sum > max)
2699  max = sum;
2700  }
2701 
2702  return max;
2703 }
2704 
2705 
2713 template <int dim, typename Number>
2714 inline Number
2716 {
2718  for (unsigned int i = 0; i < dim; ++i)
2719  {
2721  for (unsigned int j = 0; j < dim; ++j)
2722  sum += std::fabs(t[i][j]);
2723 
2724  if (sum > max)
2725  max = sum;
2726  }
2727 
2728  return max;
2729 }
2730 
2732 
2733 
2734 #ifndef DOXYGEN
2735 
2736 
2737 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2738 
2739 // Specialization of functions for ADOL-C number types when
2740 // the advanced branching feature is used
2741 template <int dim>
2742 inline adouble
2744 {
2746  for (unsigned int j = 0; j < dim; ++j)
2747  {
2749  for (unsigned int i = 0; i < dim; ++i)
2750  sum += std::fabs(t[i][j]);
2751 
2752  condassign(max, (sum > max), sum, max);
2753  }
2754 
2755  return max;
2756 }
2757 
2758 
2759 template <int dim>
2760 inline adouble
2762 {
2764  for (unsigned int i = 0; i < dim; ++i)
2765  {
2767  for (unsigned int j = 0; j < dim; ++j)
2768  sum += std::fabs(t[i][j]);
2769 
2770  condassign(max, (sum > max), sum, max);
2771  }
2772 
2773  return max;
2774 }
2775 
2776 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2777 
2778 
2779 #endif // DOXYGEN
2780 
2782 
2783 #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:2493
Number * begin_raw()
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1560
constexpr ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:2276
constexpr Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2638
Contents is actually a matrix.
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:2586
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:1708
constexpr Tensor operator-() const
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2373
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2689
#define AssertIndexRange(index, range)
Definition: exceptions.h:1628
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:2398
constexpr void clear()
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2715
constexpr Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:755
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:510
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:1513
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)
const LAPACKFullMatrix< number > & get_svd_u() const
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:1403
void copy_from(const MatrixType &)
Tensor< 2, dim, Number > project_onto_orthogonal_tensors(const Tensor< 2, dim, Number > &tensor)
Definition: tensor.h:2656
constexpr Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2474
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:361
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:92
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:2197
const LAPACKFullMatrix< number > & get_svd_vt() const
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:2617
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:2313
void copy_to(Tensor< 2, dim > &T, const size_type src_r_i=0, const size_type src_r_j=dim - 1, const size_type src_c_i=0, const size_type src_c_j=dim - 1, const unsigned int dst_r=0, const unsigned int dst_c=0) const
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:1769
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > division_operator(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1876
void copy_from(const MatrixType &)
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:1785
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:360
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
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:1988
static const bool value
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:2342
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:2124
Expression operator!=(const Expression &lhs, const Expression &rhs)
T max(const T &t, const MPI_Comm &mpi_communicator)
void mmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
#define DEAL_II_CONSTEXPR
Definition: config.h:152
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:2432