Reference documentation for deal.II version Git 65001fc 2018-06-18 17:06:20 -0400
tensor.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2018 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 
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/numbers.h>
23 #include <deal.II/base/table_indices.h>
24 #include <deal.II/base/template_constraints.h>
25 #include <deal.II/base/tensor_accessors.h>
26 #include <deal.II/base/utilities.h>
27 
28 #ifdef DEAL_II_WITH_ADOLC
29 # include <adolc/adouble.h> // Taped double
30 #endif
31 
32 #include <cmath>
33 #include <ostream>
34 #include <vector>
35 
36 
37 DEAL_II_NAMESPACE_OPEN
38 
39 // Forward declarations:
40 
41 template <int dim, typename Number>
42 class Point;
43 template <int rank_, int dim, typename Number = double>
44 class Tensor;
45 template <typename Number>
46 class Vector;
47 template <typename Number>
48 class VectorizedArray;
49 
50 #ifndef DOXYGEN
51 // Overload invalid tensor types of negative rank that come up during
52 // overload resolution of operator* and related contraction variants.
53 template <int dim, typename Number>
54 class Tensor<-2, dim, Number>
55 {};
56 
57 template <int dim, typename Number>
58 class Tensor<-1, dim, Number>
59 {};
60 #endif /* DOXYGEN */
61 
62 
93 template <int dim, typename Number>
94 class Tensor<0, dim, Number>
95 {
96 public:
105  static const unsigned int dimension = dim;
106 
110  static const unsigned int rank = 0;
111 
115  static const unsigned int n_independent_components = 1;
116 
126 
131  typedef Number value_type;
132 
138  typedef Number array_type;
139 
145  DEAL_II_CUDA_HOST_DEV
146  Tensor();
147 
153  template <typename OtherNumber>
154  Tensor(const Tensor<0, dim, OtherNumber> &initializer);
155 
159  template <typename OtherNumber>
160  Tensor(const OtherNumber &initializer);
161 
165  Number *
166  begin_raw();
167 
171  const Number *
172  begin_raw() const;
173 
177  Number *
178  end_raw();
179 
184  const Number *
185  end_raw() const;
186 
196  DEAL_II_CUDA_HOST_DEV operator Number &();
197 
206  DEAL_II_CUDA_HOST_DEV operator const Number &() const;
207 
213  template <typename OtherNumber>
214  Tensor &
216 
217 #ifdef __INTEL_COMPILER
218 
224  Tensor &
225  operator=(const Tensor<0, dim, Number> &rhs);
226 #endif
227 
232  template <typename OtherNumber>
233  Tensor &
234  operator=(const OtherNumber &d);
235 
239  template <typename OtherNumber>
240  bool
241  operator==(const Tensor<0, dim, OtherNumber> &rhs) const;
242 
246  template <typename OtherNumber>
247  bool
248  operator!=(const Tensor<0, dim, OtherNumber> &rhs) const;
249 
253  template <typename OtherNumber>
254  Tensor &
256 
260  template <typename OtherNumber>
261  Tensor &
263 
269  template <typename OtherNumber>
270  DEAL_II_CUDA_HOST_DEV Tensor &
271  operator*=(const OtherNumber &factor);
272 
276  template <typename OtherNumber>
277  Tensor &
278  operator/=(const OtherNumber &factor);
279 
283  Tensor
284  operator-() const;
285 
298  void
299  clear();
300 
306  real_type
307  norm() const;
308 
315  DEAL_II_CUDA_HOST_DEV real_type
316  norm_square() const;
317 
322  template <class Archive>
323  void
324  serialize(Archive &ar, const unsigned int version);
325 
330  typedef Number tensor_type;
331 
332 private:
336  Number value;
337 
341  template <typename OtherNumber>
342  void
343  unroll_recursion(Vector<OtherNumber> &result,
344  unsigned int & start_index) const;
345 
349  template <int, int, typename>
350  friend class Tensor;
351 };
352 
353 
354 
396 template <int rank_, int dim, typename Number>
397 class Tensor
398 {
399 public:
408  static const unsigned int dimension = dim;
409 
413  static const unsigned int rank = rank_;
414 
419  static const unsigned int n_independent_components =
420  Tensor<rank_ - 1, dim>::n_independent_components * dim;
421 
427  typedef typename Tensor<rank_ - 1, dim, Number>::tensor_type value_type;
428 
433  typedef typename Tensor<rank_ - 1, dim, Number>::array_type
434  array_type[(dim != 0) ? dim : 1];
435  // ... avoid a compiler warning in case of dim == 0 and ensure that the
436  // array always has positive size.
437 
443  DEAL_II_CUDA_HOST_DEV
444  Tensor();
445 
449  explicit Tensor(const array_type &initializer);
450 
456  template <typename OtherNumber>
457  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
458 
462  template <typename OtherNumber>
463  Tensor(
464  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
465 
469  template <typename OtherNumber>
470  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
471 
477  DEAL_II_CUDA_HOST_DEV value_type &operator[](const unsigned int i);
478 
484  DEAL_II_CUDA_HOST_DEV const value_type &
485  operator[](const unsigned int i) const;
486 
490  const Number &operator[](const TableIndices<rank_> &indices) const;
491 
495  Number &operator[](const TableIndices<rank_> &indices);
496 
500  Number *
501  begin_raw();
502 
506  const Number *
507  begin_raw() const;
508 
512  Number *
513  end_raw();
514 
518  const Number *
519  end_raw() const;
520 
526  template <typename OtherNumber>
527  Tensor &
529 
536  Tensor &
537  operator=(const Number &d);
538 
542  template <typename OtherNumber>
543  bool
545 
549  template <typename OtherNumber>
550  bool
552 
556  template <typename OtherNumber>
557  Tensor &
559 
563  template <typename OtherNumber>
564  Tensor &
566 
573  template <typename OtherNumber>
574  DEAL_II_CUDA_HOST_DEV Tensor &
575  operator*=(const OtherNumber &factor);
576 
580  template <typename OtherNumber>
581  Tensor &
582  operator/=(const OtherNumber &factor);
583 
587  Tensor
588  operator-() const;
589 
602  void
603  clear();
604 
612  norm() const;
613 
620  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
621  norm_square() const;
622 
630  template <typename OtherNumber>
631  void
632  unroll(Vector<OtherNumber> &result) const;
633 
638  static unsigned int
640 
645  static TableIndices<rank_>
646  unrolled_to_component_indices(const unsigned int i);
647 
652  static std::size_t
654 
659  template <class Archive>
660  void
661  serialize(Archive &ar, const unsigned int version);
662 
668 
669 private:
673  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
674  // ... avoid a compiler warning in case of dim == 0 and ensure that the
675  // array always has positive size.
676 
680  template <typename OtherNumber>
681  void
682  unroll_recursion(Vector<OtherNumber> &result,
683  unsigned int & start_index) const;
684 
688  template <int, int, typename>
689  friend class Tensor;
690 
695  friend class Point<dim, Number>;
696 };
697 
698 
699 namespace internal
700 {
714  template <int rank, int dim, typename T>
715  struct NumberType<Tensor<rank, dim, T>>
716  {
717  static const Tensor<rank, dim, T> &
718  value(const Tensor<rank, dim, T> &t)
719  {
720  return t;
721  }
722 
723  static Tensor<rank, dim, T>
724  value(const T &t)
725  {
727  tmp = t;
728  return tmp;
729  }
730  };
731 
732  template <int rank, int dim, typename T>
733  struct NumberType<Tensor<rank, dim, VectorizedArray<T>>>
734  {
736  value(const Tensor<rank, dim, VectorizedArray<T>> &t)
737  {
738  return t;
739  }
740 
742  value(const T &t)
743  {
746  return tmp;
747  }
748 
750  value(const VectorizedArray<T> &t)
751  {
753  tmp = t;
754  return tmp;
755  }
756  };
757 } // namespace internal
758 
759 
760 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
761 
762 
763 template <int dim, typename Number>
764 inline DEAL_II_CUDA_HOST_DEV
766  // Some auto-differentiable numbers need explicit
767  // zero initialization.
768  : value(internal::NumberType<Number>::value(0.0))
769 {}
770 
771 
772 
773 template <int dim, typename Number>
774 template <typename OtherNumber>
775 inline Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
776 {
777  value = internal::NumberType<Number>::value(initializer);
778 }
779 
780 
781 
782 template <int dim, typename Number>
783 template <typename OtherNumber>
785 {
786  value = p.value;
787 }
788 
789 
790 
791 template <int dim, typename Number>
792 inline Number *
794 {
795  return std::addressof(value);
796 }
797 
798 
799 
800 template <int dim, typename Number>
801 inline const Number *
803 {
804  return std::addressof(value);
805 }
806 
807 
808 
809 template <int dim, typename Number>
810 inline Number *
812 {
814 }
815 
816 
817 
818 template <int dim, typename Number>
819 inline const Number *
821 {
823 }
824 
825 
826 
827 template <int dim, typename Number>
828 inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::operator Number &()
829 {
830  // We cannot use Assert inside a CUDA kernel
831 #ifndef __CUDA_ARCH__
832  Assert(dim != 0,
833  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
834 #endif
835  return value;
836 }
837 
838 
839 template <int dim, typename Number>
840 inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::
841  operator const Number &() const
842 {
843  // We cannot use Assert inside a CUDA kernel
844 #ifndef __CUDA_ARCH__
845  Assert(dim != 0,
846  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
847 #endif
848  return value;
849 }
850 
851 
852 template <int dim, typename Number>
853 template <typename OtherNumber>
854 inline Tensor<0, dim, Number> &
856 {
858  return *this;
859 }
860 
861 
862 #ifdef __INTEL_COMPILER
863 template <int dim, typename Number>
864 inline Tensor<0, dim, Number> &
866 {
867  value = p.value;
868  return *this;
869 }
870 #endif
871 
872 
873 template <int dim, typename Number>
874 template <typename OtherNumber>
875 inline Tensor<0, dim, Number> &
876 Tensor<0, dim, Number>::operator=(const OtherNumber &d)
877 {
879  return *this;
880 }
881 
882 
883 template <int dim, typename Number>
884 template <typename OtherNumber>
885 inline bool
887 {
888 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
889  Assert(!(std::is_same<Number, adouble>::value ||
890  std::is_same<OtherNumber, adouble>::value),
891  ExcMessage(
892  "The Tensor equality operator for Adol-C taped numbers has not yet "
893  "been extended to support advanced branching."));
894 #endif
895 
896  return numbers::values_are_equal(value, p.value);
897 }
898 
899 
900 template <int dim, typename Number>
901 template <typename OtherNumber>
902 inline bool
904 {
905  return !((*this) == p);
906 }
907 
908 
909 template <int dim, typename Number>
910 template <typename OtherNumber>
911 inline Tensor<0, dim, Number> &
913 {
914  value += p.value;
915  return *this;
916 }
917 
918 
919 template <int dim, typename Number>
920 template <typename OtherNumber>
921 inline Tensor<0, dim, Number> &
923 {
924  value -= p.value;
925  return *this;
926 }
927 
928 
929 template <int dim, typename Number>
930 template <typename OtherNumber>
931 inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
932 Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
933 {
934  value *= s;
935  return *this;
936 }
937 
938 
939 template <int dim, typename Number>
940 template <typename OtherNumber>
941 inline Tensor<0, dim, Number> &
942 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
943 {
944  value /= s;
945  return *this;
946 }
947 
948 
949 template <int dim, typename Number>
952 {
953  return -value;
954 }
955 
956 
957 template <int dim, typename Number>
958 inline typename Tensor<0, dim, Number>::real_type
960 {
961  Assert(dim != 0,
962  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
964 }
965 
966 
967 template <int dim, typename Number>
968 inline typename Tensor<0, dim, Number>::real_type DEAL_II_CUDA_HOST_DEV
970 {
971  // We cannot use Assert inside a CUDA kernel
972 #ifndef __CUDA_ARCH__
973  Assert(dim != 0,
974  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
975 #endif
977 }
978 
979 
980 template <int dim, typename Number>
981 template <typename OtherNumber>
982 inline void
983 Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
984  unsigned int & index) const
985 {
986  Assert(dim != 0,
987  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
988  result[index] = value;
989  ++index;
990 }
991 
992 
993 template <int dim, typename Number>
994 inline void
996 {
997  // Some auto-differentiable numbers need explicit
998  // zero initialization.
1000 }
1001 
1002 
1003 template <int dim, typename Number>
1004 template <class Archive>
1005 inline void
1006 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1007 {
1008  ar &value;
1009 }
1010 
1011 
1012 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1013 
1014 
1015 template <int rank_, int dim, typename Number>
1016 inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1018 {
1019  // All members of the c-style array values are already default initialized
1020  // and thus all values are already set to zero recursively.
1021 }
1022 
1023 
1024 template <int rank_, int dim, typename Number>
1025 inline DEAL_II_ALWAYS_INLINE
1026 Tensor<rank_, dim, Number>::Tensor(const array_type &initializer)
1027 {
1028  for (unsigned int i = 0; i < dim; ++i)
1029  values[i] = Tensor<rank_ - 1, dim, Number>(initializer[i]);
1030 }
1031 
1032 
1033 template <int rank_, int dim, typename Number>
1034 template <typename OtherNumber>
1035 inline DEAL_II_ALWAYS_INLINE
1037  const Tensor<rank_, dim, OtherNumber> &initializer)
1038 {
1039  for (unsigned int i = 0; i != dim; ++i)
1040  values[i] = Tensor<rank_ - 1, dim, Number>(initializer[i]);
1041 }
1042 
1043 
1044 template <int rank_, int dim, typename Number>
1045 template <typename OtherNumber>
1046 inline DEAL_II_ALWAYS_INLINE
1048  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1049 {
1050  for (unsigned int i = 0; i < dim; ++i)
1051  values[i] = initializer[i];
1052 }
1053 
1054 
1055 template <int rank_, int dim, typename Number>
1056 template <typename OtherNumber>
1057 inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number>::
1058  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1059 {
1060  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
1061 }
1062 
1063 
1064 
1065 namespace internal
1066 {
1067  namespace TensorSubscriptor
1068  {
1069  template <typename ArrayElementType, int dim>
1070  inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV ArrayElementType &
1071  subscript(ArrayElementType * values,
1072  const unsigned int i,
1073  std::integral_constant<int, dim>)
1074  {
1075  // We cannot use Assert in a CUDA kernel
1076 #ifndef __CUDA_ARCH__
1077  Assert(i < dim, ExcIndexRange(i, 0, dim));
1078 #endif
1079  return values[i];
1080  }
1081 
1082 
1083  template <typename ArrayElementType>
1084  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1085  subscript(ArrayElementType *,
1086  const unsigned int,
1087  std::integral_constant<int, 0>)
1088  {
1089  // We cannot use Assert in a CUDA kernel
1090 #ifndef __CUDA_ARCH__
1091  Assert(
1092  false,
1093  ExcMessage(
1094  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1095 #endif
1096  static ArrayElementType t;
1097  return t;
1098  }
1099  } // namespace TensorSubscriptor
1100 } // namespace internal
1101 
1102 
1103 template <int rank_, int dim, typename Number>
1104 inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1106  operator[](const unsigned int i)
1107 {
1108  return ::internal::TensorSubscriptor::subscript(
1109  values, i, std::integral_constant<int, dim>());
1110 }
1111 
1112 
1113 template <int rank_, int dim, typename Number>
1114 inline DEAL_II_ALWAYS_INLINE
1115  DEAL_II_CUDA_HOST_DEV const typename Tensor<rank_, dim, Number>::value_type &
1116  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1117 {
1118  return ::internal::TensorSubscriptor::subscript(
1119  values, i, std::integral_constant<int, dim>());
1120 }
1121 
1122 
1123 template <int rank_, int dim, typename Number>
1124 inline const Number &Tensor<rank_, dim, Number>::
1125  operator[](const TableIndices<rank_> &indices) const
1126 {
1127  Assert(dim != 0,
1128  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1129 
1130  return TensorAccessors::extract<rank_>(*this, indices);
1131 }
1132 
1133 
1134 
1135 template <int rank_, int dim, typename Number>
1136 inline Number &Tensor<rank_, dim, Number>::
1138 {
1139  Assert(dim != 0,
1140  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1141 
1142  return TensorAccessors::extract<rank_>(*this, indices);
1143 }
1144 
1145 
1146 
1147 template <int rank_, int dim, typename Number>
1148 inline Number *
1150 {
1151  return std::addressof(
1152  this->operator[](this->unrolled_to_component_indices(0)));
1153 }
1154 
1155 
1156 
1157 template <int rank_, int dim, typename Number>
1158 inline const Number *
1160 {
1161  return std::addressof(
1162  this->operator[](this->unrolled_to_component_indices(0)));
1163 }
1164 
1165 
1166 
1167 template <int rank_, int dim, typename Number>
1168 inline Number *
1170 {
1171  return begin_raw() + n_independent_components;
1172 }
1173 
1174 
1175 
1176 template <int rank_, int dim, typename Number>
1177 inline const Number *
1179 {
1180  return begin_raw() + n_independent_components;
1181 }
1182 
1183 
1184 
1185 template <int rank_, int dim, typename Number>
1186 template <typename OtherNumber>
1187 inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
1189 {
1190  if (dim > 0)
1191  std::copy(&t.values[0], &t.values[0] + dim, &values[0]);
1192  return *this;
1193 }
1194 
1195 
1196 template <int rank_, int dim, typename Number>
1197 inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
1199 {
1201  ExcMessage("Only assignment with zero is allowed"));
1202  (void)d;
1203 
1204  for (unsigned int i = 0; i < dim; ++i)
1205  values[i] = internal::NumberType<Number>::value(0.0);
1206  return *this;
1207 }
1208 
1209 
1210 template <int rank_, int dim, typename Number>
1211 template <typename OtherNumber>
1212 inline bool
1215 {
1216  for (unsigned int i = 0; i < dim; ++i)
1217  if (values[i] != p.values[i])
1218  return false;
1219  return true;
1220 }
1221 
1222 
1223 // At some places in the library, we have Point<0> for formal reasons
1224 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1225 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1226 // in the above function that the loop end check always fails, we
1227 // implement this function here
1228 template <>
1229 template <>
1230 inline bool
1232 {
1233  return true;
1234 }
1235 
1236 
1237 template <int rank_, int dim, typename Number>
1238 template <typename OtherNumber>
1239 inline bool
1242 {
1243  return !((*this) == p);
1244 }
1245 
1246 
1247 template <int rank_, int dim, typename Number>
1248 template <typename OtherNumber>
1251 {
1252  for (unsigned int i = 0; i < dim; ++i)
1253  values[i] += p.values[i];
1254  return *this;
1255 }
1256 
1257 
1258 template <int rank_, int dim, typename Number>
1259 template <typename OtherNumber>
1262 {
1263  for (unsigned int i = 0; i < dim; ++i)
1264  values[i] -= p.values[i];
1265  return *this;
1266 }
1267 
1268 
1269 template <int rank_, int dim, typename Number>
1270 template <typename OtherNumber>
1271 inline DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
1272 Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1273 {
1274  for (unsigned int i = 0; i < dim; ++i)
1275  values[i] *= s;
1276  return *this;
1277 }
1278 
1279 
1280 template <int rank_, int dim, typename Number>
1281 template <typename OtherNumber>
1283 Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1284 {
1285  for (unsigned int i = 0; i < dim; ++i)
1286  values[i] /= s;
1287  return *this;
1288 }
1289 
1290 
1291 template <int rank_, int dim, typename Number>
1294 {
1296 
1297  for (unsigned int i = 0; i < dim; ++i)
1298  tmp.values[i] = -values[i];
1299 
1300  return tmp;
1301 }
1302 
1303 
1304 template <int rank_, int dim, typename Number>
1307 {
1308  return std::sqrt(norm_square());
1309 }
1310 
1311 
1312 template <int rank_, int dim, typename Number>
1313 inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
1315 {
1317  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1318  for (unsigned int i = 0; i < dim; ++i)
1319  s += values[i].norm_square();
1320 
1321  return s;
1322 }
1323 
1324 
1325 template <int rank_, int dim, typename Number>
1326 template <typename OtherNumber>
1327 inline void
1328 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
1329 {
1330  AssertDimension(result.size(),
1331  (Utilities::fixed_power<rank_, unsigned int>(dim)));
1332 
1333  unsigned int index = 0;
1334  unroll_recursion(result, index);
1335 }
1336 
1337 
1338 template <int rank_, int dim, typename Number>
1339 template <typename OtherNumber>
1340 inline void
1342  unsigned int & index) const
1343 {
1344  for (unsigned int i = 0; i < dim; ++i)
1345  values[i].unroll_recursion(result, index);
1346 }
1347 
1348 
1349 template <int rank_, int dim, typename Number>
1350 inline unsigned int
1352  const TableIndices<rank_> &indices)
1353 {
1354  unsigned int index = 0;
1355  for (int r = 0; r < rank_; ++r)
1356  index = index * dim + indices[r];
1357 
1358  return index;
1359 }
1360 
1361 
1362 template <int rank_, int dim, typename Number>
1363 inline TableIndices<rank_>
1365 {
1366  Assert(i < n_independent_components,
1367  ExcIndexRange(i, 0, n_independent_components));
1368 
1369  TableIndices<rank_> indices;
1370 
1371  unsigned int remainder = i;
1372  for (int r = rank_ - 1; r >= 0; --r)
1373  {
1374  indices[r] = (remainder % dim);
1375  remainder /= dim;
1376  }
1377  Assert(remainder == 0, ExcInternalError());
1378 
1379  return indices;
1380 }
1381 
1382 
1383 template <int rank_, int dim, typename Number>
1384 inline void
1386 {
1387  for (unsigned int i = 0; i < dim; ++i)
1388  values[i] = internal::NumberType<Number>::value(0.0);
1389 }
1390 
1391 
1392 template <int rank_, int dim, typename Number>
1393 inline std::size_t
1395 {
1396  return sizeof(Tensor<rank_, dim, Number>);
1397 }
1398 
1399 
1400 template <int rank_, int dim, typename Number>
1401 template <class Archive>
1402 inline void
1403 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1404 {
1405  ar &values;
1406 }
1407 
1408 
1409 /* ----------------- Non-member functions operating on tensors. ------------ */
1410 
1415 
1423 template <int rank_, int dim, typename Number>
1424 inline std::ostream &
1425 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1426 {
1427  for (unsigned int i = 0; i < dim; ++i)
1428  {
1429  out << p[i];
1430  if (i != dim - 1)
1431  out << ' ';
1432  }
1433 
1434  return out;
1435 }
1436 
1437 
1444 template <int dim, typename Number>
1445 inline std::ostream &
1446 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1447 {
1448  out << static_cast<const Number &>(p);
1449  return out;
1450 }
1451 
1452 
1454 
1458 
1459 
1468 template <int dim, typename Number, typename Other>
1469 inline DEAL_II_ALWAYS_INLINE typename ProductType<Other, Number>::type
1470 operator*(const Other &object, const Tensor<0, dim, Number> &t)
1471 {
1472  return object * static_cast<const Number &>(t);
1473 }
1474 
1475 
1476 
1485 template <int dim, typename Number, typename Other>
1486 inline DEAL_II_ALWAYS_INLINE typename ProductType<Number, Other>::type
1487 operator*(const Tensor<0, dim, Number> &t, const Other &object)
1488 {
1489  return static_cast<const Number &>(t) * object;
1490 }
1491 
1492 
1502 template <int dim, typename Number, typename OtherNumber>
1503 inline DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
1505  const Tensor<0, dim, OtherNumber> &src2)
1506 {
1507  return static_cast<const Number &>(src1) *
1508  static_cast<const OtherNumber &>(src2);
1509 }
1510 
1511 
1517 template <int dim, typename Number, typename OtherNumber>
1518 inline DEAL_II_ALWAYS_INLINE
1519  Tensor<0,
1520  dim,
1521  typename ProductType<Number,
1522  typename EnableIfScalar<OtherNumber>::type>::type>
1523  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
1524 {
1525  return static_cast<const Number &>(t) / factor;
1526 }
1527 
1528 
1534 template <int dim, typename Number, typename OtherNumber>
1535 inline DEAL_II_ALWAYS_INLINE
1538  const Tensor<0, dim, OtherNumber> &q)
1539 {
1540  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1541 }
1542 
1543 
1549 template <int dim, typename Number, typename OtherNumber>
1550 inline DEAL_II_ALWAYS_INLINE
1553  const Tensor<0, dim, OtherNumber> &q)
1554 {
1555  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1556 }
1557 
1558 
1569 template <int rank, int dim, typename Number, typename OtherNumber>
1570 inline DEAL_II_ALWAYS_INLINE
1571  Tensor<rank,
1572  dim,
1573  typename ProductType<Number,
1574  typename EnableIfScalar<OtherNumber>::type>::type>
1575  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1576 {
1577  // recurse over the base objects
1579  for (unsigned int d = 0; d < dim; ++d)
1580  tt[d] = t[d] * factor;
1581  return tt;
1582 }
1583 
1584 
1595 template <int rank, int dim, typename Number, typename OtherNumber>
1596 inline DEAL_II_ALWAYS_INLINE
1597  Tensor<rank,
1598  dim,
1600  OtherNumber>::type>
1601  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
1602 {
1603  // simply forward to the operator above
1604  return t * factor;
1605 }
1606 
1607 
1615 template <int rank, int dim, typename Number, typename OtherNumber>
1616 inline Tensor<
1617  rank,
1618  dim,
1619  typename ProductType<Number,
1620  typename EnableIfScalar<OtherNumber>::type>::type>
1621 operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1622 {
1623  // recurse over the base objects
1625  for (unsigned int d = 0; d < dim; ++d)
1626  tt[d] = t[d] / factor;
1627  return tt;
1628 }
1629 
1630 
1638 template <int rank, int dim, typename Number, typename OtherNumber>
1639 inline DEAL_II_ALWAYS_INLINE
1643 {
1645 
1646  for (unsigned int i = 0; i < dim; ++i)
1647  tmp[i] += q[i];
1648 
1649  return tmp;
1650 }
1651 
1652 
1660 template <int rank, int dim, typename Number, typename OtherNumber>
1661 inline DEAL_II_ALWAYS_INLINE
1665 {
1667 
1668  for (unsigned int i = 0; i < dim; ++i)
1669  tmp[i] -= q[i];
1670 
1671  return tmp;
1672 }
1673 
1674 
1676 
1680 
1681 
1705 template <int rank_1,
1706  int rank_2,
1707  int dim,
1708  typename Number,
1709  typename OtherNumber>
1710 inline DEAL_II_ALWAYS_INLINE
1711  typename Tensor<rank_1 + rank_2 - 2,
1712  dim,
1713  typename ProductType<Number, OtherNumber>::type>::tensor_type
1716 {
1717  typename Tensor<rank_1 + rank_2 - 2,
1718  dim,
1719  typename ProductType<Number, OtherNumber>::type>::tensor_type
1720  result;
1721 
1722  TensorAccessors::internal::
1723  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
1724  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
1725  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
1726 
1727  return result;
1728 }
1729 
1730 
1760 template <int index_1,
1761  int index_2,
1762  int rank_1,
1763  int rank_2,
1764  int dim,
1765  typename Number,
1766  typename OtherNumber>
1767 inline DEAL_II_ALWAYS_INLINE
1768  typename Tensor<rank_1 + rank_2 - 2,
1769  dim,
1770  typename ProductType<Number, OtherNumber>::type>::tensor_type
1773 {
1774  Assert(0 <= index_1 && index_1 < rank_1,
1775  ExcMessage(
1776  "The specified index_1 must lie within the range [0,rank_1)"));
1777  Assert(0 <= index_2 && index_2 < rank_2,
1778  ExcMessage(
1779  "The specified index_2 must lie within the range [0,rank_2)"));
1780 
1781  using namespace TensorAccessors;
1782  using namespace TensorAccessors::internal;
1783 
1784  // Reorder index_1 to the end of src1:
1785  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>
1786  reord_01 = reordered_index_view<index_1, rank_1>(src1);
1787 
1788  // Reorder index_2 to the end of src2:
1789  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>
1790  reord_02 = reordered_index_view<index_2, rank_2>(src2);
1791 
1792  typename Tensor<rank_1 + rank_2 - 2,
1793  dim,
1794  typename ProductType<Number, OtherNumber>::type>::tensor_type
1795  result;
1796  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
1797  return result;
1798 }
1799 
1800 
1832 template <int index_1,
1833  int index_2,
1834  int index_3,
1835  int index_4,
1836  int rank_1,
1837  int rank_2,
1838  int dim,
1839  typename Number,
1840  typename OtherNumber>
1841 inline
1842  typename Tensor<rank_1 + rank_2 - 4,
1843  dim,
1844  typename ProductType<Number, OtherNumber>::type>::tensor_type
1847 {
1848  Assert(0 <= index_1 && index_1 < rank_1,
1849  ExcMessage(
1850  "The specified index_1 must lie within the range [0,rank_1)"));
1851  Assert(0 <= index_3 && index_3 < rank_1,
1852  ExcMessage(
1853  "The specified index_3 must lie within the range [0,rank_1)"));
1854  Assert(index_1 != index_3,
1855  ExcMessage("index_1 and index_3 must not be the same"));
1856  Assert(0 <= index_2 && index_2 < rank_2,
1857  ExcMessage(
1858  "The specified index_2 must lie within the range [0,rank_2)"));
1859  Assert(0 <= index_4 && index_4 < rank_2,
1860  ExcMessage(
1861  "The specified index_4 must lie within the range [0,rank_2)"));
1862  Assert(index_2 != index_4,
1863  ExcMessage("index_2 and index_4 must not be the same"));
1864 
1865  using namespace TensorAccessors;
1866  using namespace TensorAccessors::internal;
1867 
1868  // Reorder index_1 to the end of src1:
1869  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>
1870  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
1871 
1872  // Reorder index_2 to the end of src2:
1873  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>
1874  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
1875 
1876  // Now, reorder index_3 to the end of src1. We have to make sure to
1877  // preserve the orginial ordering: index_1 has been removed. If
1878  // index_3 > index_1, we have to use (index_3 - 1) instead:
1879  ReorderedIndexView<
1880  (index_3 < index_1 ? index_3 : index_3 - 1),
1881  rank_1,
1882  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
1883  reord_3 =
1884  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
1885  index_3 - 1,
1886  rank_1 > (reord_1);
1887 
1888  // Now, reorder index_4 to the end of src2. We have to make sure to
1889  // preserve the orginial ordering: index_2 has been removed. If
1890  // index_4 > index_2, we have to use (index_4 - 1) instead:
1891  ReorderedIndexView<
1892  (index_4 < index_2 ? index_4 : index_4 - 1),
1893  rank_2,
1894  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
1895  reord_4 =
1896  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
1897  index_4 - 1,
1898  rank_2 > (reord_2);
1899 
1900  typename Tensor<rank_1 + rank_2 - 4,
1901  dim,
1902  typename ProductType<Number, OtherNumber>::type>::tensor_type
1903  result;
1904  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
1905  return result;
1906 }
1907 
1908 
1922 template <int rank, int dim, typename Number, typename OtherNumber>
1923 inline DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
1925  const Tensor<rank, dim, OtherNumber> &right)
1926 {
1927  typename ProductType<Number, OtherNumber>::type result;
1928  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
1929  return result;
1930 }
1931 
1932 
1951 template <template <int, int, typename> class TensorT1,
1952  template <int, int, typename> class TensorT2,
1953  template <int, int, typename> class TensorT3,
1954  int rank_1,
1955  int rank_2,
1956  int dim,
1957  typename T1,
1958  typename T2,
1959  typename T3>
1961 contract3(const TensorT1<rank_1, dim, T1> & left,
1962  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
1963  const TensorT3<rank_2, dim, T3> & right)
1964 {
1966  return_type;
1967  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
1968  middle,
1969  right);
1970 }
1971 
1972 
1984 template <int rank_1,
1985  int rank_2,
1986  int dim,
1987  typename Number,
1988  typename OtherNumber>
1989 inline DEAL_II_ALWAYS_INLINE
1993 {
1994  typename Tensor<rank_1 + rank_2,
1995  dim,
1996  typename ProductType<Number, OtherNumber>::type>::tensor_type
1997  result;
1998  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
1999  return result;
2000 }
2001 
2002 
2004 
2008 
2009 
2021 template <int dim, typename Number>
2022 inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
2024 {
2025  Assert(dim == 2, ExcInternalError());
2026 
2027  Tensor<1, dim, Number> result;
2028 
2029  result[0] = src[1];
2030  result[1] = -src[0];
2031 
2032  return result;
2033 }
2034 
2035 
2046 template <int dim, typename Number>
2047 inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
2049  const Tensor<1, dim, Number> &src2)
2050 {
2051  Assert(dim == 3, ExcInternalError());
2052 
2053  Tensor<1, dim, Number> result;
2054 
2055  result[0] = src1[1] * src2[2] - src1[2] * src2[1];
2056  result[1] = src1[2] * src2[0] - src1[0] * src2[2];
2057  result[2] = src1[0] * src2[1] - src1[1] * src2[0];
2058 
2059  return result;
2060 }
2061 
2062 
2064 
2068 
2069 
2076 template <int dim, typename Number>
2077 inline Number
2079 {
2080  // Compute the determinant using the Laplace expansion of the
2081  // determinant. We expand along the last row.
2082  Number det = internal::NumberType<Number>::value(0.0);
2083 
2084  for (unsigned int k = 0; k < dim; ++k)
2085  {
2086  Tensor<2, dim - 1, Number> minor;
2087  for (unsigned int i = 0; i < dim - 1; ++i)
2088  for (unsigned int j = 0; j < dim - 1; ++j)
2089  minor[i][j] = t[i][j < k ? j : j + 1];
2090 
2091  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2092 
2093  det += t[dim - 1][k] * cofactor;
2094  }
2095 
2096  return ((dim % 2 == 0) ? 1. : -1.) * det;
2097 }
2098 
2104 template <typename Number>
2105 inline Number
2107 {
2108  return t[0][0];
2109 }
2110 
2111 
2119 template <int dim, typename Number>
2120 inline DEAL_II_ALWAYS_INLINE Number
2122 {
2123  Number t = d[0][0];
2124  for (unsigned int i = 1; i < dim; ++i)
2125  t += d[i][i];
2126  return t;
2127 }
2128 
2129 
2139 template <int dim, typename Number>
2142 {
2143  Number return_tensor[dim][dim];
2144 
2145  // if desired, take over the
2146  // inversion of a 4x4 tensor
2147  // from the FullMatrix
2148  AssertThrow(false, ExcNotImplemented());
2149 
2150  return Tensor<2, dim, Number>(return_tensor);
2151 }
2152 
2153 
2154 #ifndef DOXYGEN
2155 
2156 template <typename Number>
2157 inline Tensor<2, 1, Number>
2158 invert(const Tensor<2, 1, Number> &t)
2159 {
2160  Number return_tensor[1][1];
2161 
2162  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2163 
2164  return Tensor<2, 1, Number>(return_tensor);
2165 }
2166 
2167 
2168 template <typename Number>
2169 inline Tensor<2, 2, Number>
2170 invert(const Tensor<2, 2, Number> &t)
2171 {
2172  Tensor<2, 2, Number> return_tensor;
2173 
2174  // this is Maple output,
2175  // thus a bit unstructured
2176  const Number inv_det_t = internal::NumberType<Number>::value(
2177  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2178  return_tensor[0][0] = t[1][1];
2179  return_tensor[0][1] = -t[0][1];
2180  return_tensor[1][0] = -t[1][0];
2181  return_tensor[1][1] = t[0][0];
2182  return_tensor *= inv_det_t;
2183 
2184  return return_tensor;
2185 }
2186 
2187 
2188 template <typename Number>
2189 inline Tensor<2, 3, Number>
2190 invert(const Tensor<2, 3, Number> &t)
2191 {
2192  Tensor<2, 3, Number> return_tensor;
2193 
2194  const Number t4 = internal::NumberType<Number>::value(t[0][0] * t[1][1]),
2195  t6 = internal::NumberType<Number>::value(t[0][0] * t[1][2]),
2196  t8 = internal::NumberType<Number>::value(t[0][1] * t[1][0]),
2197  t00 = internal::NumberType<Number>::value(t[0][2] * t[1][0]),
2198  t01 = internal::NumberType<Number>::value(t[0][1] * t[2][0]),
2199  t04 = internal::NumberType<Number>::value(t[0][2] * t[2][0]),
2201  1.0 / (t4 * t[2][2] - t6 * t[2][1] - t8 * t[2][2] +
2202  t00 * t[2][1] + t01 * t[1][2] - t04 * t[1][1]));
2203  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2204  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2205  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2206  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2207  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2208  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2209  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2210  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2211  return_tensor[1][1] =
2212  internal::NumberType<Number>::value(t[0][0] * t[2][2]) - t04;
2213  return_tensor[1][2] = t00 - t6;
2214  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2215  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2216  return_tensor[2][1] =
2217  t01 - internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2218  return_tensor[2][2] = internal::NumberType<Number>::value(t4 - t8);
2219  return_tensor *= inv_det_t;
2220 
2221  return return_tensor;
2222 }
2223 
2224 #endif /* DOXYGEN */
2225 
2226 
2233 template <int dim, typename Number>
2234 inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, Number>
2236 {
2238  for (unsigned int i = 0; i < dim; ++i)
2239  {
2240  tt[i][i] = t[i][i];
2241  for (unsigned int j = i + 1; j < dim; ++j)
2242  {
2243  tt[i][j] = t[j][i];
2244  tt[j][i] = t[i][j];
2245  };
2246  }
2247  return tt;
2248 }
2249 
2250 
2264 template <int dim, typename Number>
2267 {
2268  return determinant(t) * invert(t);
2269 }
2270 
2271 
2286 template <int dim, typename Number>
2289 {
2290  return transpose(adjugate(t));
2291 }
2292 
2293 
2301 template <int dim, typename Number>
2302 inline Number
2304 {
2305  Number max = internal::NumberType<Number>::value(0.0);
2306  for (unsigned int j = 0; j < dim; ++j)
2307  {
2308  Number sum = internal::NumberType<Number>::value(0.0);
2309  for (unsigned int i = 0; i < dim; ++i)
2310  sum += std::fabs(t[i][j]);
2311 
2312  if (sum > max)
2313  max = sum;
2314  }
2315 
2316  return max;
2317 }
2318 
2319 
2327 template <int dim, typename Number>
2328 inline Number
2330 {
2331  Number max = internal::NumberType<Number>::value(0.0);
2332  for (unsigned int i = 0; i < dim; ++i)
2333  {
2334  Number sum = internal::NumberType<Number>::value(0.0);
2335  for (unsigned int j = 0; j < dim; ++j)
2336  sum += std::fabs(t[i][j]);
2337 
2338  if (sum > max)
2339  max = sum;
2340  }
2341 
2342  return max;
2343 }
2344 
2346 
2347 
2348 #ifndef DOXYGEN
2349 
2350 
2351 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2352 
2353 // Specialization of functions for Adol-C number types when
2354 // the advanced branching feature is used
2355 template <int dim>
2356 inline adouble
2357 l1_norm(const Tensor<2, dim, adouble> &t)
2358 {
2359  adouble max = internal::NumberType<adouble>::value(0.0);
2360  for (unsigned int j = 0; j < dim; ++j)
2361  {
2362  adouble sum = internal::NumberType<adouble>::value(0.0);
2363  for (unsigned int i = 0; i < dim; ++i)
2364  sum += std::fabs(t[i][j]);
2365 
2366  condassign(max, (sum > max), sum, max);
2367  }
2368 
2369  return max;
2370 }
2371 
2372 
2373 template <int dim>
2374 inline adouble
2375 linfty_norm(const Tensor<2, dim, adouble> &t)
2376 {
2377  adouble max = internal::NumberType<adouble>::value(0.0);
2378  for (unsigned int i = 0; i < dim; ++i)
2379  {
2380  adouble sum = internal::NumberType<adouble>::value(0.0);
2381  for (unsigned int j = 0; j < dim; ++j)
2382  sum += std::fabs(t[i][j]);
2383 
2384  condassign(max, (sum > max), sum, max);
2385  }
2386 
2387  return max;
2388 }
2389 
2390 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2391 
2392 
2393 #endif // DOXYGEN
2394 
2395 DEAL_II_NAMESPACE_CLOSE
2396 
2397 // include deprecated non-member functions operating on Tensor
2398 #include <deal.II/base/tensor_deprecated.h>
2399 
2400 #endif
numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:125
Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2121
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1366
Number determinant(const SymmetricTensor< 2, dim, Number > &)
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
Definition: tensor.h:1351
Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1621
static const unsigned int n_independent_components
Definition: tensor.h:419
bool value_is_zero(const Number &value)
Definition: numbers.h:802
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2303
static std::size_t memory_consumption()
Definition: tensor.h:1394
Tensor< rank_-1, dim, Number > values[(dim!=0)?dim:1]
Definition: tensor.h:673
Tensor & operator/=(const OtherNumber &factor)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2329
Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2266
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:1771
static const unsigned int rank
Definition: tensor.h:413
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1306
#define AssertThrow(cond, exc)
Definition: exceptions.h:1329
static real_type abs(const number &x)
Definition: numbers.h:483
bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1241
internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
Definition: point.h:106
bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:786
void unroll(Vector< OtherNumber > &result) const
Definition: tensor.h:1328
LinearAlgebra::distributed::Vector< Number > Vector
Tensor & operator*=(const OtherNumber &factor)
Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2288
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2235
Tensor< rank_-1, dim, Number >::tensor_type value_type
Definition: tensor.h:427
static::ExceptionBase & ExcMessage(std::string arg1)
Tensor< rank, dim, typename ProductType< typename EnableIfScalar< Number >::type, OtherNumber >::type > operator*(const Number &factor, const Tensor< rank, dim, OtherNumber > &t)
Definition: tensor.h:1601
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:1961
Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:2141
static real_type abs_square(const number &x)
Definition: numbers.h:474
Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2023
bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1214
#define Assert(cond, exc)
Definition: exceptions.h:1227
ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:1924
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: tensor.h:1341
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:667
Tensor< rank_-1, dim, Number >::array_type array_type[(dim!=0)?dim:1]
Definition: tensor.h:434
Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
Definition: tensor.h:1537
void serialize(Archive &ar, const unsigned int version)
Definition: tensor.h:1403
Tensor< 0, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const Tensor< 0, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1523
Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1575
Number * end_raw()
Definition: tensor.h:1169
value_type & operator[](const unsigned int i)
Definition: tensor.h:1106
Tensor()
Definition: tensor.h:1017
Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< rank, dim, Number > &p, const Tensor< rank, dim, OtherNumber > &q)
Definition: tensor.h:1641
Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Definition: tensor.h:1364
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:1845
Tensor< 1, dim, Number > cross_product_3d(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:2048
Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2078
static const unsigned int dimension
Definition: tensor.h:408
Tensor operator-() const
Definition: tensor.h:1293
Definition: mpi.h:55
Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< rank, dim, Number > &p, const Tensor< rank, dim, OtherNumber > &q)
Definition: tensor.h:1663
ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1470
Number * begin_raw()
Definition: tensor.h:1149
static::ExceptionBase & ExcNotImplemented()
ProductType< Number, Other >::type operator*(const Tensor< 0, dim, Number > &t, const Other &object)
Definition: tensor.h:1487
ProductType< Number, OtherNumber >::type operator*(const Tensor< 0, dim, Number > &src1, const Tensor< 0, dim, OtherNumber > &src2)
Definition: tensor.h:1504
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:1991
numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1314
Tensor & operator+=(const Tensor< rank_, dim, OtherNumber > &)
Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
Definition: tensor.h:1552
void clear()
Definition: tensor.h:1385
Tensor< rank_1+rank_2-2, dim, typename ProductType< Number, OtherNumber >::type >::tensor_type operator*(const Tensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
Definition: tensor.h:1714
static::ExceptionBase & ExcInternalError()
Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:2106