Reference documentation for deal.II version 9.2.0
\(\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>
27 #include <deal.II/base/utilities.h>
28 
30 
31 #ifdef DEAL_II_WITH_ADOLC
32 # include <adolc/adouble.h> // Taped double
33 #endif
34 
35 #include <cmath>
36 #include <ostream>
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 
92 template <int dim, typename Number>
93 class Tensor<0, dim, Number>
94 {
95 public:
96  static_assert(dim >= 0,
97  "Tensors must have a dimension greater than or equal to one.");
98 
107  static constexpr unsigned int dimension = dim;
108 
112  static constexpr unsigned int rank = 0;
113 
117  static constexpr unsigned int n_independent_components = 1;
118 
128 
133  using value_type = Number;
134 
140  using array_type = Number;
141 
147  constexpr DEAL_II_CUDA_HOST_DEV
148  Tensor();
149 
157  template <typename OtherNumber>
158  constexpr DEAL_II_CUDA_HOST_DEV
159  Tensor(const Tensor<0, dim, OtherNumber> &initializer);
160 
166  template <typename OtherNumber>
167  constexpr DEAL_II_CUDA_HOST_DEV
168  Tensor(const OtherNumber &initializer);
169 
173  Number *
174  begin_raw();
175 
179  const Number *
180  begin_raw() const;
181 
185  Number *
186  end_raw();
187 
192  const Number *
193  end_raw() const;
194 
204  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV operator Number &();
205 
214  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV operator const Number &() const;
215 
223  template <typename OtherNumber>
225  operator=(const Tensor<0, dim, OtherNumber> &rhs);
226 
227 #ifdef __INTEL_COMPILER
228 
237  operator=(const Tensor<0, dim, Number> &rhs);
238 #endif
239 
246  template <typename OtherNumber>
248  operator=(const OtherNumber &d);
249 
253  template <typename OtherNumber>
254  DEAL_II_CONSTEXPR bool
255  operator==(const Tensor<0, dim, OtherNumber> &rhs) const;
256 
260  template <typename OtherNumber>
261  constexpr bool
262  operator!=(const Tensor<0, dim, OtherNumber> &rhs) const;
263 
269  template <typename OtherNumber>
271  operator+=(const Tensor<0, dim, OtherNumber> &rhs);
272 
278  template <typename OtherNumber>
280  operator-=(const Tensor<0, dim, OtherNumber> &rhs);
281 
287  template <typename OtherNumber>
289  operator*=(const OtherNumber &factor);
290 
296  template <typename OtherNumber>
298  operator/=(const OtherNumber &factor);
299 
305  constexpr DEAL_II_CUDA_HOST_DEV Tensor
306  operator-() const;
307 
320  DEAL_II_CONSTEXPR void
321  clear();
322 
328  real_type
329  norm() const;
330 
338  norm_square() const;
339 
344  template <class Archive>
345  void
346  serialize(Archive &ar, const unsigned int version);
347 
352  using tensor_type = Number;
353 
354 private:
358  Number value;
359 
363  template <typename OtherNumber>
364  void
365  unroll_recursion(Vector<OtherNumber> &result,
366  unsigned int & start_index) const;
367 
368  // Allow an arbitrary Tensor to access the underlying values.
369  template <int, int, typename>
370  friend class Tensor;
371 };
372 
373 
374 
449 template <int rank_, int dim, typename Number>
450 class Tensor
451 {
452 public:
453  static_assert(rank_ >= 0,
454  "Tensors must have a rank greater than or equal to one.");
455  static_assert(dim >= 0,
456  "Tensors must have a dimension greater than or equal to one.");
465  static constexpr unsigned int dimension = dim;
466 
470  static constexpr unsigned int rank = rank_;
471 
476  static constexpr unsigned int n_independent_components =
477  Tensor<rank_ - 1, dim>::n_independent_components * dim;
478 
484  using value_type = typename Tensor<rank_ - 1, dim, Number>::tensor_type;
485 
490  using array_type =
491  typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1];
492 
499  Tensor()
500 #ifdef DEAL_II_MSVC
501  : values{}
502  {}
503 #else
504  = default;
505 #endif
506 
512  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
513  const array_type &initializer);
514 
522  template <typename OtherNumber>
523  constexpr DEAL_II_CUDA_HOST_DEV
524  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
525 
529  template <typename OtherNumber>
530  constexpr Tensor(
531  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
532 
536  template <typename OtherNumber>
537  constexpr
538  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
539 
546  operator[](const unsigned int i);
547 
553  constexpr DEAL_II_CUDA_HOST_DEV const value_type &
554  operator[](const unsigned int i) const;
555 
559  DEAL_II_CONSTEXPR const Number &
560  operator[](const TableIndices<rank_> &indices) const;
561 
565  DEAL_II_CONSTEXPR Number &operator[](const TableIndices<rank_> &indices);
566 
570  Number *
571  begin_raw();
572 
576  const Number *
577  begin_raw() const;
578 
582  Number *
583  end_raw();
584 
588  const Number *
589  end_raw() const;
590 
598  template <typename OtherNumber>
600  operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
601 
608  DEAL_II_CONSTEXPR Tensor &
609  operator=(const Number &d);
610 
614  template <typename OtherNumber>
615  DEAL_II_CONSTEXPR bool
617 
621  template <typename OtherNumber>
622  constexpr bool
624 
630  template <typename OtherNumber>
632  operator+=(const Tensor<rank_, dim, OtherNumber> &);
633 
639  template <typename OtherNumber>
641  operator-=(const Tensor<rank_, dim, OtherNumber> &);
642 
649  template <typename OtherNumber>
651  operator*=(const OtherNumber &factor);
652 
658  template <typename OtherNumber>
660  operator/=(const OtherNumber &factor);
661 
668  operator-() const;
669 
682  DEAL_II_CONSTEXPR void
683  clear();
684 
694  norm() const;
695 
704  norm_square() const;
705 
713  template <typename OtherNumber>
714  void
715  unroll(Vector<OtherNumber> &result) const;
716 
721  static DEAL_II_CONSTEXPR unsigned int
722  component_to_unrolled_index(const TableIndices<rank_> &indices);
723 
730  unrolled_to_component_indices(const unsigned int i);
731 
736  static constexpr std::size_t
738 
743  template <class Archive>
744  void
745  serialize(Archive &ar, const unsigned int version);
746 
752 
753 private:
757  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
758  // ... avoid a compiler warning in case of dim == 0 and ensure that the
759  // array always has positive size.
760 
764  template <typename OtherNumber>
765  void
766  unroll_recursion(Vector<OtherNumber> &result,
767  unsigned int & start_index) const;
768 
775  template <typename ArrayLike, std::size_t... Indices>
776  constexpr DEAL_II_CUDA_HOST_DEV
777  Tensor(const ArrayLike &initializer, std_cxx14::index_sequence<Indices...>);
778 
779  // Allow an arbitrary Tensor to access the underlying values.
780  template <int, int, typename>
781  friend class Tensor;
782 
783  // Point is allowed access to the coordinates. This is supposed to improve
784  // speed.
785  friend class Point<dim, Number>;
786 };
787 
788 
789 #ifndef DOXYGEN
790 namespace internal
791 {
792  // Workaround: The following 4 overloads are necessary to be able to
793  // compile the library with Apple Clang 8 and older. We should remove
794  // these overloads again when we bump the minimal required version to
795  // something later than clang-3.6 / Apple Clang 6.3.
796  // - Jean-Paul Pelteret, Matthias Maier, Daniel Arndt 2020
797  template <int rank, int dim, typename T, typename U>
798  struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
799  {
800  using type =
802  };
803 
804  template <int rank, int dim, typename T, typename U>
805  struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
806  {
807  using type =
809  };
810 
811  template <typename T, int rank, int dim, typename U>
812  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
813  {
814  using type =
816  };
817 
818  template <int rank, int dim, typename T, typename U>
819  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
820  {
821  using type =
823  };
824  // end workaround
825 
830  template <int rank, int dim, typename T>
831  struct NumberType<Tensor<rank, dim, T>>
832  {
833  static constexpr DEAL_II_ALWAYS_INLINE const Tensor<rank, dim, T> &
834  value(const Tensor<rank, dim, T> &t)
835  {
836  return t;
837  }
838 
840  value(const T &t)
841  {
843  tmp = t;
844  return tmp;
845  }
846  };
847 } // namespace internal
848 
849 
850 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
851 
852 
853 template <int dim, typename Number>
856  // Some auto-differentiable numbers need explicit
857  // zero initialization such as adtl::adouble.
858  : Tensor{0.0}
859 {}
860 
861 
862 
863 template <int dim, typename Number>
864 template <typename OtherNumber>
866  Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
867  : value(internal::NumberType<Number>::value(initializer))
868 {}
869 
870 
871 
872 template <int dim, typename Number>
873 template <typename OtherNumber>
876  : Tensor{p.value}
877 {}
878 
879 
880 
881 template <int dim, typename Number>
882 inline Number *
884 {
885  return std::addressof(value);
886 }
887 
888 
889 
890 template <int dim, typename Number>
891 inline const Number *
893 {
894  return std::addressof(value);
895 }
896 
897 
898 
899 template <int dim, typename Number>
900 inline Number *
902 {
904 }
905 
906 
907 
908 template <int dim, typename Number>
909 const Number *
911 {
913 }
914 
915 
916 
917 template <int dim, typename Number>
920 {
921  // We cannot use Assert inside a CUDA kernel
922 # ifndef __CUDA_ARCH__
923  Assert(dim != 0,
924  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
925 # endif
926  return value;
927 }
928 
929 
930 template <int dim, typename Number>
933 {
934  // We cannot use Assert inside a CUDA kernel
935 # ifndef __CUDA_ARCH__
936  Assert(dim != 0,
937  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
938 # endif
939  return value;
940 }
941 
942 
943 template <int dim, typename Number>
944 template <typename OtherNumber>
948 {
950  return *this;
951 }
952 
953 
954 # ifdef __INTEL_COMPILER
955 template <int dim, typename Number>
959 {
960  value = p.value;
961  return *this;
962 }
963 # endif
964 
965 
966 template <int dim, typename Number>
967 template <typename OtherNumber>
970  Tensor<0, dim, Number>::operator=(const OtherNumber &d)
971 {
973  return *this;
974 }
975 
976 
977 template <int dim, typename Number>
978 template <typename OtherNumber>
979 DEAL_II_CONSTEXPR inline bool
981 {
982 # if defined(DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING)
985  ExcMessage(
986  "The Tensor equality operator for ADOL-C taped numbers has not yet "
987  "been extended to support advanced branching."));
988 # endif
989 
990  return numbers::values_are_equal(value, p.value);
991 }
992 
993 
994 template <int dim, typename Number>
995 template <typename OtherNumber>
996 constexpr bool
998 {
999  return !((*this) == p);
1000 }
1001 
1002 
1003 template <int dim, typename Number>
1004 template <typename OtherNumber>
1008 {
1009  value += p.value;
1010  return *this;
1011 }
1012 
1013 
1014 template <int dim, typename Number>
1015 template <typename OtherNumber>
1019 {
1020  value -= p.value;
1021  return *this;
1022 }
1023 
1024 
1025 
1026 namespace internal
1027 {
1028  namespace ComplexWorkaround
1029  {
1030  template <typename Number, typename OtherNumber>
1032  multiply_assign_scalar(Number &val, const OtherNumber &s)
1033  {
1034  val *= s;
1035  }
1036 
1037 # ifdef __CUDA_ARCH__
1038  template <typename Number, typename OtherNumber>
1040  multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
1041  {
1042  printf("This function is not implemented for std::complex<Number>!\n");
1043  assert(false);
1044  }
1045 # endif
1046  } // namespace ComplexWorkaround
1047 } // namespace internal
1048 
1049 
1050 template <int dim, typename Number>
1051 template <typename OtherNumber>
1054  Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1055 {
1056  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1057  return *this;
1058 }
1059 
1060 
1061 
1062 template <int dim, typename Number>
1063 template <typename OtherNumber>
1065 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1066 {
1067  value /= s;
1068  return *this;
1069 }
1070 
1071 
1072 template <int dim, typename Number>
1075 {
1076  return -value;
1077 }
1078 
1079 
1080 template <int dim, typename Number>
1083 {
1084  Assert(dim != 0,
1085  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1087 }
1088 
1089 
1090 template <int dim, typename Number>
1094 {
1095  // We cannot use Assert inside a CUDA kernel
1096 # ifndef __CUDA_ARCH__
1097  Assert(dim != 0,
1098  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1099 # endif
1101 }
1102 
1103 
1104 template <int dim, typename Number>
1105 template <typename OtherNumber>
1106 inline void
1107 Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1108  unsigned int & index) const
1109 {
1110  Assert(dim != 0,
1111  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1112  result[index] = value;
1113  ++index;
1114 }
1115 
1116 
1117 template <int dim, typename Number>
1118 DEAL_II_CONSTEXPR inline void
1120 {
1121  // Some auto-differentiable numbers need explicit
1122  // zero initialization.
1124 }
1125 
1126 
1127 template <int dim, typename Number>
1128 template <class Archive>
1129 inline void
1130 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1131 {
1132  ar &value;
1133 }
1134 
1135 
1136 template <int dim, typename Number>
1138 
1139 
1140 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1141 
1142 template <int rank_, int dim, typename Number>
1143 template <typename ArrayLike, std::size_t... indices>
1145  Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1146  std_cxx14::index_sequence<indices...>)
1147  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
1148 {
1149  static_assert(sizeof...(indices) == dim,
1150  "dim should match the number of indices");
1151 }
1152 
1153 
1154 template <int rank_, int dim, typename Number>
1157  : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
1158 {}
1159 
1160 
1161 template <int rank_, int dim, typename Number>
1162 template <typename OtherNumber>
1165  const Tensor<rank_, dim, OtherNumber> &initializer)
1166  : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
1167 {}
1168 
1169 
1170 template <int rank_, int dim, typename Number>
1171 template <typename OtherNumber>
1172 constexpr DEAL_II_ALWAYS_INLINE
1174  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1175  : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
1176 {}
1177 
1178 
1179 template <int rank_, int dim, typename Number>
1180 template <typename OtherNumber>
1182  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1183 {
1184  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
1185 }
1186 
1187 
1188 
1189 namespace internal
1190 {
1191  namespace TensorSubscriptor
1192  {
1193  template <typename ArrayElementType, int dim>
1195  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1196  subscript(ArrayElementType * values,
1197  const unsigned int i,
1198  std::integral_constant<int, dim>)
1199  {
1200  // We cannot use Assert in a CUDA kernel
1201 # ifndef __CUDA_ARCH__
1202  AssertIndexRange(i, dim);
1203 # endif
1204  return values[i];
1205  }
1206 
1207  // The variables within this struct will be referenced in the next function.
1208  // It is a workaround that allows returning a reference to a static variable
1209  // while allowing constexpr evaluation of the function.
1210  // It has to be defined outside the function because constexpr functions
1211  // cannot define static variables
1212  template <typename ArrayElementType>
1213  struct Uninitialized
1214  {
1215  static ArrayElementType value;
1216  };
1217 
1218  template <typename Type>
1220 
1221  template <typename ArrayElementType>
1223  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1224  subscript(ArrayElementType *,
1225  const unsigned int,
1226  std::integral_constant<int, 0>)
1227  {
1228  // We cannot use Assert in a CUDA kernel
1229 # ifndef __CUDA_ARCH__
1230  Assert(
1231  false,
1232  ExcMessage(
1233  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1234 # endif
1236  }
1237  } // namespace TensorSubscriptor
1238 } // namespace internal
1239 
1240 
1241 template <int rank_, int dim, typename Number>
1244  operator[](const unsigned int i)
1245 {
1246  return ::internal::TensorSubscriptor::subscript(
1247  values, i, std::integral_constant<int, dim>());
1248 }
1249 
1250 
1251 template <int rank_, int dim, typename Number>
1252 constexpr DEAL_II_ALWAYS_INLINE
1254  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1255 {
1256  return values[i];
1257 }
1258 
1259 
1260 template <int rank_, int dim, typename Number>
1261 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1263  operator[](const TableIndices<rank_> &indices) const
1264 {
1265  Assert(dim != 0,
1266  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1267 
1268  return TensorAccessors::extract<rank_>(*this, indices);
1269 }
1270 
1271 
1272 
1273 template <int rank_, int dim, typename Number>
1276 {
1277  Assert(dim != 0,
1278  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1279 
1280  return TensorAccessors::extract<rank_>(*this, indices);
1281 }
1282 
1283 
1284 
1285 template <int rank_, int dim, typename Number>
1286 inline Number *
1288 {
1289  return std::addressof(
1290  this->operator[](this->unrolled_to_component_indices(0)));
1291 }
1292 
1293 
1294 
1295 template <int rank_, int dim, typename Number>
1296 inline const Number *
1298 {
1299  return std::addressof(
1300  this->operator[](this->unrolled_to_component_indices(0)));
1301 }
1302 
1303 
1304 
1305 template <int rank_, int dim, typename Number>
1306 inline Number *
1308 {
1310 }
1311 
1312 
1313 
1314 template <int rank_, int dim, typename Number>
1315 inline const Number *
1317 {
1319 }
1320 
1321 
1322 
1323 template <int rank_, int dim, typename Number>
1324 template <typename OtherNumber>
1327 {
1328  // The following loop could be written more concisely using std::copy, but
1329  // that function is only constexpr from C++20 on.
1330  for (unsigned int i = 0; i < dim; ++i)
1331  values[i] = t.values[i];
1332  return *this;
1333 }
1334 
1335 
1336 template <int rank_, int dim, typename Number>
1339 {
1341  ExcMessage("Only assignment with zero is allowed"));
1342  (void)d;
1343 
1344  for (unsigned int i = 0; i < dim; ++i)
1346  return *this;
1347 }
1348 
1349 
1350 template <int rank_, int dim, typename Number>
1351 template <typename OtherNumber>
1352 DEAL_II_CONSTEXPR inline bool
1355 {
1356  for (unsigned int i = 0; i < dim; ++i)
1357  if (values[i] != p.values[i])
1358  return false;
1359  return true;
1360 }
1361 
1362 
1363 // At some places in the library, we have Point<0> for formal reasons
1364 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1365 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1366 // in the above function that the loop end check always fails, we
1367 // implement this function here
1368 template <>
1369 template <>
1370 DEAL_II_CONSTEXPR inline bool
1372 {
1373  return true;
1374 }
1375 
1376 
1377 template <int rank_, int dim, typename Number>
1378 template <typename OtherNumber>
1379 constexpr bool
1382 {
1383  return !((*this) == p);
1384 }
1385 
1386 
1387 template <int rank_, int dim, typename Number>
1388 template <typename OtherNumber>
1393 {
1394  for (unsigned int i = 0; i < dim; ++i)
1395  values[i] += p.values[i];
1396  return *this;
1397 }
1398 
1399 
1400 template <int rank_, int dim, typename Number>
1401 template <typename OtherNumber>
1406 {
1407  for (unsigned int i = 0; i < dim; ++i)
1408  values[i] -= p.values[i];
1409  return *this;
1410 }
1411 
1412 
1413 template <int rank_, int dim, typename Number>
1414 template <typename OtherNumber>
1417  Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1418 {
1419  for (unsigned int i = 0; i < dim; ++i)
1420  values[i] *= s;
1421  return *this;
1422 }
1423 
1424 
1425 namespace internal
1426 {
1427  namespace TensorImplementation
1428  {
1429  template <int rank,
1430  int dim,
1431  typename Number,
1432  typename OtherNumber,
1433  typename std::enable_if<
1434  !std::is_integral<
1437  int>::type = 0>
1440  const OtherNumber &factor)
1441  {
1442  const Number inverse_factor = Number(1.) / factor;
1443  // recurse over the base objects
1444  for (unsigned int d = 0; d < dim; ++d)
1445  t[d] *= inverse_factor;
1446  }
1447 
1448 
1449  template <int rank,
1450  int dim,
1451  typename Number,
1452  typename OtherNumber,
1453  typename std::enable_if<
1454  std::is_integral<
1457  int>::type = 0>
1460  const OtherNumber &factor)
1461  {
1462  // recurse over the base objects
1463  for (unsigned int d = 0; d < dim; ++d)
1464  t[d] /= factor;
1465  }
1466  } // namespace TensorImplementation
1467 } // namespace internal
1468 
1469 
1470 template <int rank_, int dim, typename Number>
1471 template <typename OtherNumber>
1474  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1475 {
1477  return *this;
1478 }
1479 
1480 
1481 template <int rank_, int dim, typename Number>
1485 {
1487 
1488  for (unsigned int i = 0; i < dim; ++i)
1489  tmp.values[i] = -values[i];
1490 
1491  return tmp;
1492 }
1493 
1494 
1495 template <int rank_, int dim, typename Number>
1498 {
1499  return std::sqrt(norm_square());
1500 }
1501 
1502 
1503 template <int rank_, int dim, typename Number>
1507 {
1510  for (unsigned int i = 0; i < dim; ++i)
1511  s += values[i].norm_square();
1512 
1513  return s;
1514 }
1515 
1516 
1517 template <int rank_, int dim, typename Number>
1518 template <typename OtherNumber>
1519 inline void
1520 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
1521 {
1522  AssertDimension(result.size(),
1523  (Utilities::fixed_power<rank_, unsigned int>(dim)));
1524 
1525  unsigned int index = 0;
1526  unroll_recursion(result, index);
1527 }
1528 
1529 
1530 template <int rank_, int dim, typename Number>
1531 template <typename OtherNumber>
1532 inline void
1533 Tensor<rank_, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1534  unsigned int & index) const
1535 {
1536  for (unsigned int i = 0; i < dim; ++i)
1537  values[i].unroll_recursion(result, index);
1538 }
1539 
1540 
1541 template <int rank_, int dim, typename Number>
1542 DEAL_II_CONSTEXPR inline unsigned int
1544  const TableIndices<rank_> &indices)
1545 {
1546  unsigned int index = 0;
1547  for (int r = 0; r < rank_; ++r)
1548  index = index * dim + indices[r];
1549 
1550  return index;
1551 }
1552 
1553 
1554 
1555 namespace internal
1556 {
1557  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1558  // and rank=2. Make sure we don't have compiler warnings.
1559 
1560  template <int dim>
1561  inline DEAL_II_CONSTEXPR unsigned int
1562  mod(const unsigned int x)
1563  {
1564  return x % dim;
1565  }
1566 
1567  template <>
1568  inline unsigned int
1569  mod<0>(const unsigned int x)
1570  {
1571  Assert(false, ExcInternalError());
1572  return x;
1573  }
1574 
1575  template <int dim>
1576  inline DEAL_II_CONSTEXPR unsigned int
1577  div(const unsigned int x)
1578  {
1579  return x / dim;
1580  }
1581 
1582  template <>
1583  inline unsigned int
1584  div<0>(const unsigned int x)
1585  {
1586  Assert(false, ExcInternalError());
1587  return x;
1588  }
1589 
1590 } // namespace internal
1591 
1592 
1593 
1594 template <int rank_, int dim, typename Number>
1597 {
1598  AssertIndexRange(i, n_independent_components);
1599 
1600  TableIndices<rank_> indices;
1601 
1602  unsigned int remainder = i;
1603  for (int r = rank_ - 1; r >= 0; --r)
1604  {
1605  indices[r] = internal::mod<dim>(remainder);
1606  remainder = internal::div<dim>(remainder);
1607  }
1608  Assert(remainder == 0, ExcInternalError());
1609 
1610  return indices;
1611 }
1612 
1613 
1614 template <int rank_, int dim, typename Number>
1615 DEAL_II_CONSTEXPR inline void
1617 {
1618  for (unsigned int i = 0; i < dim; ++i)
1620 }
1621 
1622 
1623 template <int rank_, int dim, typename Number>
1624 constexpr std::size_t
1626 {
1627  return sizeof(Tensor<rank_, dim, Number>);
1628 }
1629 
1630 
1631 template <int rank_, int dim, typename Number>
1632 template <class Archive>
1633 inline void
1634 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1635 {
1636  ar &values;
1637 }
1638 
1639 
1640 template <int rank_, int dim, typename Number>
1642 
1643 #endif // DOXYGEN
1644 
1645 /* ----------------- Non-member functions operating on tensors. ------------ */
1646 
1651 
1659 template <int rank_, int dim, typename Number>
1660 inline std::ostream &
1661 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1662 {
1663  for (unsigned int i = 0; i < dim; ++i)
1664  {
1665  out << p[i];
1666  if (i != dim - 1)
1667  out << ' ';
1668  }
1669 
1670  return out;
1671 }
1672 
1673 
1680 template <int dim, typename Number>
1681 inline std::ostream &
1682 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1683 {
1684  out << static_cast<const Number &>(p);
1685  return out;
1686 }
1687 
1688 
1690 
1694 
1695 
1706 template <int dim, typename Number, typename Other>
1709  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1710 {
1711  return object * static_cast<const Number &>(t);
1712 }
1713 
1714 
1715 
1726 template <int dim, typename Number, typename Other>
1729  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1730 {
1731  return static_cast<const Number &>(t) * object;
1732 }
1733 
1734 
1746 template <int dim, typename Number, typename OtherNumber>
1750  const Tensor<0, dim, OtherNumber> &src2)
1751 {
1752  return static_cast<const Number &>(src1) *
1753  static_cast<const OtherNumber &>(src2);
1754 }
1755 
1756 
1764 template <int dim, typename Number, typename OtherNumber>
1766  Tensor<0,
1767  dim,
1768  typename ProductType<Number,
1769  typename EnableIfScalar<OtherNumber>::type>::type>
1770  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
1771 {
1772  return static_cast<const Number &>(t) / factor;
1773 }
1774 
1775 
1783 template <int dim, typename Number, typename OtherNumber>
1787  const Tensor<0, dim, OtherNumber> &q)
1788 {
1789  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1790 }
1791 
1792 
1800 template <int dim, typename Number, typename OtherNumber>
1804  const Tensor<0, dim, OtherNumber> &q)
1805 {
1806  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1807 }
1808 
1809 
1822 template <int rank, int dim, typename Number, typename OtherNumber>
1824  Tensor<rank,
1825  dim,
1826  typename ProductType<Number,
1827  typename EnableIfScalar<OtherNumber>::type>::type>
1828  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1829 {
1830  // recurse over the base objects
1832  for (unsigned int d = 0; d < dim; ++d)
1833  tt[d] = t[d] * factor;
1834  return tt;
1835 }
1836 
1837 
1850 template <int rank, int dim, typename Number, typename OtherNumber>
1852  Tensor<rank,
1853  dim,
1855  OtherNumber>::type>
1856  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
1857 {
1858  // simply forward to the operator above
1859  return t * factor;
1860 }
1861 
1862 
1863 namespace internal
1864 {
1865  namespace TensorImplementation
1866  {
1867  template <int rank,
1868  int dim,
1869  typename Number,
1870  typename OtherNumber,
1871  typename std::enable_if<
1872  !std::is_integral<
1874  int>::type = 0>
1878  const OtherNumber & factor)
1879  {
1881  const Number inverse_factor = Number(1.) / factor;
1882  // recurse over the base objects
1883  for (unsigned int d = 0; d < dim; ++d)
1884  tt[d] = t[d] * inverse_factor;
1885  return tt;
1886  }
1887 
1888 
1889  template <int rank,
1890  int dim,
1891  typename Number,
1892  typename OtherNumber,
1893  typename std::enable_if<
1894  std::is_integral<
1896  int>::type = 0>
1900  const OtherNumber & factor)
1901  {
1903  // recurse over the base objects
1904  for (unsigned int d = 0; d < dim; ++d)
1905  tt[d] = t[d] / factor;
1906  return tt;
1907  }
1908  } // namespace TensorImplementation
1909 } // namespace internal
1910 
1911 
1921 template <int rank, int dim, typename Number, typename OtherNumber>
1923  Tensor<rank,
1924  dim,
1925  typename ProductType<Number,
1926  typename EnableIfScalar<OtherNumber>::type>::type>
1927  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1928 {
1930 }
1931 
1932 
1942 template <int rank, int dim, typename Number, typename OtherNumber>
1947 {
1949 
1950  for (unsigned int i = 0; i < dim; ++i)
1951  tmp[i] += q[i];
1952 
1953  return tmp;
1954 }
1955 
1956 
1966 template <int rank, int dim, typename Number, typename OtherNumber>
1971 {
1973 
1974  for (unsigned int i = 0; i < dim; ++i)
1975  tmp[i] -= q[i];
1976 
1977  return tmp;
1978 }
1979 
1986 template <int dim, typename Number, typename OtherNumber>
1990  const Tensor<0, dim, OtherNumber> &src2)
1991 {
1993 
1994  tmp *= src2;
1995 
1996  return tmp;
1997 }
1998 
2015 template <int rank, int dim, typename Number, typename OtherNumber>
2019  const Tensor<rank, dim, OtherNumber> &src2)
2020 {
2022 
2023  for (unsigned int i = 0; i < dim; ++i)
2024  tmp[i] = schur_product(Tensor<rank - 1, dim, Number>(src1[i]),
2026 
2027  return tmp;
2028 }
2029 
2031 
2035 
2036 
2060 template <int rank_1,
2061  int rank_2,
2062  int dim,
2063  typename Number,
2064  typename OtherNumber,
2065  typename = typename std::enable_if<rank_1 >= 1 && rank_2 >= 1>::type>
2067  typename Tensor<rank_1 + rank_2 - 2,
2068  dim,
2069  typename ProductType<Number, OtherNumber>::type>::tensor_type
2072 {
2073  typename Tensor<rank_1 + rank_2 - 2,
2074  dim,
2075  typename ProductType<Number, OtherNumber>::type>::tensor_type
2076  result{};
2077 
2078  TensorAccessors::internal::
2079  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2080  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2081  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2082 
2083  return result;
2084 }
2085 
2086 
2116 template <int index_1,
2117  int index_2,
2118  int rank_1,
2119  int rank_2,
2120  int dim,
2121  typename Number,
2122  typename OtherNumber>
2124  typename Tensor<rank_1 + rank_2 - 2,
2125  dim,
2126  typename ProductType<Number, OtherNumber>::type>::tensor_type
2129 {
2130  Assert(0 <= index_1 && index_1 < rank_1,
2131  ExcMessage(
2132  "The specified index_1 must lie within the range [0,rank_1)"));
2133  Assert(0 <= index_2 && index_2 < rank_2,
2134  ExcMessage(
2135  "The specified index_2 must lie within the range [0,rank_2)"));
2136 
2137  using namespace TensorAccessors;
2138  using namespace TensorAccessors::internal;
2139 
2140  // Reorder index_1 to the end of src1:
2142  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2143 
2144  // Reorder index_2 to the end of src2:
2146  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2147 
2148  typename Tensor<rank_1 + rank_2 - 2,
2149  dim,
2150  typename ProductType<Number, OtherNumber>::type>::tensor_type
2151  result{};
2152  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2153  return result;
2154 }
2155 
2156 
2188 template <int index_1,
2189  int index_2,
2190  int index_3,
2191  int index_4,
2192  int rank_1,
2193  int rank_2,
2194  int dim,
2195  typename Number,
2196  typename OtherNumber>
2197 DEAL_II_CONSTEXPR inline
2198  typename Tensor<rank_1 + rank_2 - 4,
2199  dim,
2200  typename ProductType<Number, OtherNumber>::type>::tensor_type
2203 {
2204  Assert(0 <= index_1 && index_1 < rank_1,
2205  ExcMessage(
2206  "The specified index_1 must lie within the range [0,rank_1)"));
2207  Assert(0 <= index_3 && index_3 < rank_1,
2208  ExcMessage(
2209  "The specified index_3 must lie within the range [0,rank_1)"));
2210  Assert(index_1 != index_3,
2211  ExcMessage("index_1 and index_3 must not be the same"));
2212  Assert(0 <= index_2 && index_2 < rank_2,
2213  ExcMessage(
2214  "The specified index_2 must lie within the range [0,rank_2)"));
2215  Assert(0 <= index_4 && index_4 < rank_2,
2216  ExcMessage(
2217  "The specified index_4 must lie within the range [0,rank_2)"));
2218  Assert(index_2 != index_4,
2219  ExcMessage("index_2 and index_4 must not be the same"));
2220 
2221  using namespace TensorAccessors;
2222  using namespace TensorAccessors::internal;
2223 
2224  // Reorder index_1 to the end of src1:
2226  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2227 
2228  // Reorder index_2 to the end of src2:
2230  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2231 
2232  // Now, reorder index_3 to the end of src1. We have to make sure to
2233  // preserve the original ordering: index_1 has been removed. If
2234  // index_3 > index_1, we have to use (index_3 - 1) instead:
2236  (index_3 < index_1 ? index_3 : index_3 - 1),
2237  rank_1,
2238  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2239  reord_3 =
2240  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2241  index_3 - 1,
2242  rank_1 > (reord_1);
2243 
2244  // Now, reorder index_4 to the end of src2. We have to make sure to
2245  // preserve the original ordering: index_2 has been removed. If
2246  // index_4 > index_2, we have to use (index_4 - 1) instead:
2247  ReorderedIndexView<
2248  (index_4 < index_2 ? index_4 : index_4 - 1),
2249  rank_2,
2250  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2251  reord_4 =
2252  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2253  index_4 - 1,
2254  rank_2 > (reord_2);
2255 
2256  typename Tensor<rank_1 + rank_2 - 4,
2257  dim,
2258  typename ProductType<Number, OtherNumber>::type>::tensor_type
2259  result{};
2260  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2261  return result;
2262 }
2263 
2264 
2278 template <int rank, int dim, typename Number, typename OtherNumber>
2282  const Tensor<rank, dim, OtherNumber> &right)
2283 {
2284  typename ProductType<Number, OtherNumber>::type result{};
2285  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2286  return result;
2287 }
2288 
2289 
2308 template <template <int, int, typename> class TensorT1,
2309  template <int, int, typename> class TensorT2,
2310  template <int, int, typename> class TensorT3,
2311  int rank_1,
2312  int rank_2,
2313  int dim,
2314  typename T1,
2315  typename T2,
2316  typename T3>
2319  contract3(const TensorT1<rank_1, dim, T1> & left,
2320  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2321  const TensorT3<rank_2, dim, T3> & right)
2322 {
2323  using return_type =
2325  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2326  middle,
2327  right);
2328 }
2329 
2330 
2342 template <int rank_1,
2343  int rank_2,
2344  int dim,
2345  typename Number,
2346  typename OtherNumber>
2351 {
2352  typename Tensor<rank_1 + rank_2,
2353  dim,
2354  typename ProductType<Number, OtherNumber>::type>::tensor_type
2355  result{};
2356  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2357  return result;
2358 }
2359 
2360 
2362 
2366 
2367 
2379 template <int dim, typename Number>
2382 {
2383  Assert(dim == 2, ExcInternalError());
2384 
2385  Tensor<1, dim, Number> result;
2386 
2387  result[0] = src[1];
2388  result[1] = -src[0];
2389 
2390  return result;
2391 }
2392 
2393 
2404 template <int dim, typename Number1, typename Number2>
2408  const Tensor<1, dim, Number2> &src2)
2409 {
2410  Assert(dim == 3, ExcInternalError());
2411 
2413 
2414  // avoid compiler warnings
2415  constexpr int s0 = 0 % dim;
2416  constexpr int s1 = 1 % dim;
2417  constexpr int s2 = 2 % dim;
2418 
2419  result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2420  result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2421  result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2422 
2423  return result;
2424 }
2425 
2426 
2428 
2432 
2433 
2440 template <int dim, typename Number>
2443 {
2444  // Compute the determinant using the Laplace expansion of the
2445  // determinant. We expand along the last row.
2446  Number det = internal::NumberType<Number>::value(0.0);
2447 
2448  for (unsigned int k = 0; k < dim; ++k)
2449  {
2450  Tensor<2, dim - 1, Number> minor;
2451  for (unsigned int i = 0; i < dim - 1; ++i)
2452  for (unsigned int j = 0; j < dim - 1; ++j)
2453  minor[i][j] = t[i][j < k ? j : j + 1];
2454 
2455  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2456 
2457  det += t[dim - 1][k] * cofactor;
2458  }
2459 
2460  return ((dim % 2 == 0) ? 1. : -1.) * det;
2461 }
2462 
2468 template <typename Number>
2469 constexpr DEAL_II_ALWAYS_INLINE Number
2471 {
2472  return t[0][0];
2473 }
2474 
2475 
2483 template <int dim, typename Number>
2486 {
2487  Number t = d[0][0];
2488  for (unsigned int i = 1; i < dim; ++i)
2489  t += d[i][i];
2490  return t;
2491 }
2492 
2493 
2503 template <int dim, typename Number>
2506 {
2507  Number return_tensor[dim][dim];
2508 
2509  // if desired, take over the
2510  // inversion of a 4x4 tensor
2511  // from the FullMatrix
2512  AssertThrow(false, ExcNotImplemented());
2513 
2514  return Tensor<2, dim, Number>(return_tensor);
2515 }
2516 
2517 
2518 #ifndef DOXYGEN
2519 
2520 template <typename Number>
2522  invert(const Tensor<2, 1, Number> &t)
2523 {
2524  Tensor<2, 1, Number> return_tensor;
2525 
2526  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2527 
2528  return return_tensor;
2529 }
2530 
2531 
2532 template <typename Number>
2534  invert(const Tensor<2, 2, Number> &t)
2535 {
2536  Tensor<2, 2, Number> return_tensor;
2537 
2538  // this is Maple output,
2539  // thus a bit unstructured
2540  const Number inv_det_t = internal::NumberType<Number>::value(
2541  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2542  return_tensor[0][0] = t[1][1];
2543  return_tensor[0][1] = -t[0][1];
2544  return_tensor[1][0] = -t[1][0];
2545  return_tensor[1][1] = t[0][0];
2546  return_tensor *= inv_det_t;
2547 
2548  return return_tensor;
2549 }
2550 
2551 
2552 template <typename Number>
2554  invert(const Tensor<2, 3, Number> &t)
2555 {
2556  Tensor<2, 3, Number> return_tensor;
2557 
2558  const Number t4 = internal::NumberType<Number>::value(t[0][0] * t[1][1]),
2559  t6 = internal::NumberType<Number>::value(t[0][0] * t[1][2]),
2560  t8 = internal::NumberType<Number>::value(t[0][1] * t[1][0]),
2561  t00 = internal::NumberType<Number>::value(t[0][2] * t[1][0]),
2562  t01 = internal::NumberType<Number>::value(t[0][1] * t[2][0]),
2563  t04 = internal::NumberType<Number>::value(t[0][2] * t[2][0]),
2565  1.0 / (t4 * t[2][2] - t6 * t[2][1] - t8 * t[2][2] +
2566  t00 * t[2][1] + t01 * t[1][2] - t04 * t[1][1]));
2567  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2568  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2569  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2570  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2571  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2572  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2573  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2574  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2575  return_tensor[1][1] =
2576  internal::NumberType<Number>::value(t[0][0] * t[2][2]) - t04;
2577  return_tensor[1][2] = t00 - t6;
2578  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2579  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2580  return_tensor[2][1] =
2581  t01 - internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2582  return_tensor[2][2] = internal::NumberType<Number>::value(t4 - t8);
2583  return_tensor *= inv_det_t;
2584 
2585  return return_tensor;
2586 }
2587 
2588 #endif /* DOXYGEN */
2589 
2590 
2597 template <int dim, typename Number>
2600 {
2602  for (unsigned int i = 0; i < dim; ++i)
2603  {
2604  tt[i][i] = t[i][i];
2605  for (unsigned int j = i + 1; j < dim; ++j)
2606  {
2607  tt[i][j] = t[j][i];
2608  tt[j][i] = t[i][j];
2609  };
2610  }
2611  return tt;
2612 }
2613 
2614 
2629 template <int dim, typename Number>
2630 constexpr Tensor<2, dim, Number>
2632 {
2633  return determinant(t) * invert(t);
2634 }
2635 
2636 
2651 template <int dim, typename Number>
2652 constexpr Tensor<2, dim, Number>
2654 {
2655  return transpose(adjugate(t));
2656 }
2657 
2658 
2669 template <int dim, typename Number>
2672 {
2673  Tensor<2, dim, Number> output_tensor;
2675  LAPACKFullMatrix<Number> lapack_matrix(dim);
2676  LAPACKFullMatrix<Number> result(dim);
2677 
2678  // todo: find or add dealii functionality to copy in one step.
2679  matrix.copy_from(tensor);
2680  lapack_matrix.copy_from(matrix);
2681 
2682  // now compute the svd of the matrices
2683  lapack_matrix.compute_svd();
2684 
2685  // Use the SVD results to orthogonalize: @f$U V^T@f$
2686  lapack_matrix.get_svd_u().mmult(result, lapack_matrix.get_svd_vt());
2687 
2688  // todo: find or add dealii functionality to copy in one step.
2689  matrix = result;
2690  matrix.copy_to(output_tensor);
2691  return output_tensor;
2692 }
2693 
2694 
2703 template <int dim, typename Number>
2704 inline Number
2706 {
2708  for (unsigned int j = 0; j < dim; ++j)
2709  {
2711  for (unsigned int i = 0; i < dim; ++i)
2712  sum += std::fabs(t[i][j]);
2713 
2714  if (sum > max)
2715  max = sum;
2716  }
2717 
2718  return max;
2719 }
2720 
2721 
2730 template <int dim, typename Number>
2731 inline Number
2733 {
2735  for (unsigned int i = 0; i < dim; ++i)
2736  {
2738  for (unsigned int j = 0; j < dim; ++j)
2739  sum += std::fabs(t[i][j]);
2740 
2741  if (sum > max)
2742  max = sum;
2743  }
2744 
2745  return max;
2746 }
2747 
2749 
2750 
2751 #ifndef DOXYGEN
2752 
2753 
2754 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2755 
2756 // Specialization of functions for ADOL-C number types when
2757 // the advanced branching feature is used
2758 template <int dim>
2759 inline adouble
2761 {
2763  for (unsigned int j = 0; j < dim; ++j)
2764  {
2766  for (unsigned int i = 0; i < dim; ++i)
2767  sum += std::fabs(t[i][j]);
2768 
2769  condassign(max, (sum > max), sum, max);
2770  }
2771 
2772  return max;
2773 }
2774 
2775 
2776 template <int dim>
2777 inline adouble
2779 {
2781  for (unsigned int i = 0; i < dim; ++i)
2782  {
2784  for (unsigned int j = 0; j < dim; ++j)
2785  sum += std::fabs(t[i][j]);
2786 
2787  condassign(max, (sum > max), sum, max);
2788  }
2789 
2790  return max;
2791 }
2792 
2793 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2794 
2795 
2796 #endif // DOXYGEN
2797 
2799 
2800 #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:2505
Number * begin_raw()
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
constexpr ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:2281
constexpr Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2653
Contents is actually a matrix.
static constexpr const T & value(const T &t)
Definition: numbers.h:703
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:925
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:2599
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:587
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1709
constexpr Tensor operator-() const
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2381
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2705
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
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:2407
constexpr void clear()
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2732
constexpr Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:757
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
STL namespace.
typename Tensor< rank_ - 1, dim, Number >::array_type[(dim !=0) ? dim :1] array_type
Definition: tensor.h:491
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
static real_type abs(const number &x)
Definition: numbers.h:609
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: point.h:111
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
Definition: c++.h:129
#define Assert(cond, exc)
Definition: exceptions.h:1419
void copy_from(const MatrixType &)
Tensor< 2, dim, Number > project_onto_orthogonal_tensors(const Tensor< 2, dim, Number > &tensor)
Definition: tensor.h:2671
constexpr Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2485
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:359
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:99
typename Tensor< rank_ - 1, dim, Number >::tensor_type value_type
Definition: tensor.h:484
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:2201
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:470
constexpr Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2631
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:2319
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:127
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:1770
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > division_operator(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1877
void copy_from(const MatrixType &)
Definition: tensor.h:450
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:1786
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:941
#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:476
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:1989
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:2349
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:2127
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:102
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:2442