Reference documentation for deal.II version Git 3fe3ec1 2018-08-21 21:05:51 -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  using value_type = Number;
132 
138  using array_type = Number;
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  using tensor_type = Number;
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  using value_type = typename Tensor<rank_ - 1, dim, Number>::tensor_type;
428 
433  using array_type =
434  typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1];
435 
441  DEAL_II_CUDA_HOST_DEV
442  Tensor();
443 
447  explicit Tensor(const array_type &initializer);
448 
454  template <typename OtherNumber>
455  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
456 
460  template <typename OtherNumber>
461  Tensor(
462  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
463 
467  template <typename OtherNumber>
468  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
469 
475  DEAL_II_CUDA_HOST_DEV value_type &operator[](const unsigned int i);
476 
482  DEAL_II_CUDA_HOST_DEV const value_type &
483  operator[](const unsigned int i) const;
484 
488  const Number &operator[](const TableIndices<rank_> &indices) const;
489 
493  Number &operator[](const TableIndices<rank_> &indices);
494 
498  Number *
499  begin_raw();
500 
504  const Number *
505  begin_raw() const;
506 
510  Number *
511  end_raw();
512 
516  const Number *
517  end_raw() const;
518 
524  template <typename OtherNumber>
525  Tensor &
527 
534  Tensor &
535  operator=(const Number &d);
536 
540  template <typename OtherNumber>
541  bool
543 
547  template <typename OtherNumber>
548  bool
550 
554  template <typename OtherNumber>
555  Tensor &
557 
561  template <typename OtherNumber>
562  Tensor &
564 
571  template <typename OtherNumber>
572  DEAL_II_CUDA_HOST_DEV Tensor &
573  operator*=(const OtherNumber &factor);
574 
578  template <typename OtherNumber>
579  Tensor &
580  operator/=(const OtherNumber &factor);
581 
585  Tensor
586  operator-() const;
587 
600  void
601  clear();
602 
610  norm() const;
611 
618  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
619  norm_square() const;
620 
628  template <typename OtherNumber>
629  void
630  unroll(Vector<OtherNumber> &result) const;
631 
636  static unsigned int
638 
643  static TableIndices<rank_>
644  unrolled_to_component_indices(const unsigned int i);
645 
650  static std::size_t
652 
657  template <class Archive>
658  void
659  serialize(Archive &ar, const unsigned int version);
660 
666 
667 private:
671  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
672  // ... avoid a compiler warning in case of dim == 0 and ensure that the
673  // array always has positive size.
674 
678  template <typename OtherNumber>
679  void
680  unroll_recursion(Vector<OtherNumber> &result,
681  unsigned int & start_index) const;
682 
686  template <int, int, typename>
687  friend class Tensor;
688 
693  friend class Point<dim, Number>;
694 };
695 
696 
697 namespace internal
698 {
712  template <int rank, int dim, typename T>
713  struct NumberType<Tensor<rank, dim, T>>
714  {
715  static const Tensor<rank, dim, T> &
716  value(const Tensor<rank, dim, T> &t)
717  {
718  return t;
719  }
720 
721  static Tensor<rank, dim, T>
722  value(const T &t)
723  {
725  tmp = t;
726  return tmp;
727  }
728  };
729 
730  template <int rank, int dim, typename T>
731  struct NumberType<Tensor<rank, dim, VectorizedArray<T>>>
732  {
734  value(const Tensor<rank, dim, VectorizedArray<T>> &t)
735  {
736  return t;
737  }
738 
740  value(const T &t)
741  {
744  return tmp;
745  }
746 
748  value(const VectorizedArray<T> &t)
749  {
751  tmp = t;
752  return tmp;
753  }
754  };
755 } // namespace internal
756 
757 
758 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
759 
760 
761 template <int dim, typename Number>
762 inline DEAL_II_CUDA_HOST_DEV
764  // Some auto-differentiable numbers need explicit
765  // zero initialization.
766  : value(internal::NumberType<Number>::value(0.0))
767 {}
768 
769 
770 
771 template <int dim, typename Number>
772 template <typename OtherNumber>
773 inline Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
774 {
775  value = internal::NumberType<Number>::value(initializer);
776 }
777 
778 
779 
780 template <int dim, typename Number>
781 template <typename OtherNumber>
783 {
784  value = p.value;
785 }
786 
787 
788 
789 template <int dim, typename Number>
790 inline Number *
792 {
793  return std::addressof(value);
794 }
795 
796 
797 
798 template <int dim, typename Number>
799 inline const Number *
801 {
802  return std::addressof(value);
803 }
804 
805 
806 
807 template <int dim, typename Number>
808 inline Number *
810 {
812 }
813 
814 
815 
816 template <int dim, typename Number>
817 inline const Number *
819 {
821 }
822 
823 
824 
825 template <int dim, typename Number>
826 inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::operator Number &()
827 {
828  // We cannot use Assert inside a CUDA kernel
829 #ifndef __CUDA_ARCH__
830  Assert(dim != 0,
831  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
832 #endif
833  return value;
834 }
835 
836 
837 template <int dim, typename Number>
838 inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::
839  operator const Number &() const
840 {
841  // We cannot use Assert inside a CUDA kernel
842 #ifndef __CUDA_ARCH__
843  Assert(dim != 0,
844  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
845 #endif
846  return value;
847 }
848 
849 
850 template <int dim, typename Number>
851 template <typename OtherNumber>
852 inline Tensor<0, dim, Number> &
854 {
856  return *this;
857 }
858 
859 
860 #ifdef __INTEL_COMPILER
861 template <int dim, typename Number>
862 inline Tensor<0, dim, Number> &
864 {
865  value = p.value;
866  return *this;
867 }
868 #endif
869 
870 
871 template <int dim, typename Number>
872 template <typename OtherNumber>
873 inline Tensor<0, dim, Number> &
874 Tensor<0, dim, Number>::operator=(const OtherNumber &d)
875 {
877  return *this;
878 }
879 
880 
881 template <int dim, typename Number>
882 template <typename OtherNumber>
883 inline bool
885 {
886 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
887  Assert(!(std::is_same<Number, adouble>::value ||
888  std::is_same<OtherNumber, adouble>::value),
889  ExcMessage(
890  "The Tensor equality operator for Adol-C taped numbers has not yet "
891  "been extended to support advanced branching."));
892 #endif
893 
894  return numbers::values_are_equal(value, p.value);
895 }
896 
897 
898 template <int dim, typename Number>
899 template <typename OtherNumber>
900 inline bool
902 {
903  return !((*this) == p);
904 }
905 
906 
907 template <int dim, typename Number>
908 template <typename OtherNumber>
909 inline Tensor<0, dim, Number> &
911 {
912  value += p.value;
913  return *this;
914 }
915 
916 
917 template <int dim, typename Number>
918 template <typename OtherNumber>
919 inline Tensor<0, dim, Number> &
921 {
922  value -= p.value;
923  return *this;
924 }
925 
926 
927 template <int dim, typename Number>
928 template <typename OtherNumber>
929 inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
930 Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
931 {
932  value *= s;
933  return *this;
934 }
935 
936 
937 template <int dim, typename Number>
938 template <typename OtherNumber>
939 inline Tensor<0, dim, Number> &
940 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
941 {
942  value /= s;
943  return *this;
944 }
945 
946 
947 template <int dim, typename Number>
950 {
951  return -value;
952 }
953 
954 
955 template <int dim, typename Number>
956 inline typename Tensor<0, dim, Number>::real_type
958 {
959  Assert(dim != 0,
960  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
962 }
963 
964 
965 template <int dim, typename Number>
966 inline typename Tensor<0, dim, Number>::real_type DEAL_II_CUDA_HOST_DEV
968 {
969  // We cannot use Assert inside a CUDA kernel
970 #ifndef __CUDA_ARCH__
971  Assert(dim != 0,
972  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
973 #endif
975 }
976 
977 
978 template <int dim, typename Number>
979 template <typename OtherNumber>
980 inline void
981 Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
982  unsigned int & index) const
983 {
984  Assert(dim != 0,
985  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
986  result[index] = value;
987  ++index;
988 }
989 
990 
991 template <int dim, typename Number>
992 inline void
994 {
995  // Some auto-differentiable numbers need explicit
996  // zero initialization.
998 }
999 
1000 
1001 template <int dim, typename Number>
1002 template <class Archive>
1003 inline void
1004 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1005 {
1006  ar &value;
1007 }
1008 
1009 
1010 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1011 
1012 
1013 template <int rank_, int dim, typename Number>
1014 inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1016 {
1017  // All members of the c-style array values are already default initialized
1018  // and thus all values are already set to zero recursively.
1019 }
1020 
1021 
1022 template <int rank_, int dim, typename Number>
1023 inline DEAL_II_ALWAYS_INLINE
1025 {
1026  for (unsigned int i = 0; i < dim; ++i)
1027  values[i] = Tensor<rank_ - 1, dim, Number>(initializer[i]);
1028 }
1029 
1030 
1031 template <int rank_, int dim, typename Number>
1032 template <typename OtherNumber>
1033 inline DEAL_II_ALWAYS_INLINE
1035  const Tensor<rank_, dim, OtherNumber> &initializer)
1036 {
1037  for (unsigned int i = 0; i != dim; ++i)
1038  values[i] = Tensor<rank_ - 1, dim, Number>(initializer[i]);
1039 }
1040 
1041 
1042 template <int rank_, int dim, typename Number>
1043 template <typename OtherNumber>
1044 inline DEAL_II_ALWAYS_INLINE
1046  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1047 {
1048  for (unsigned int i = 0; i < dim; ++i)
1049  values[i] = initializer[i];
1050 }
1051 
1052 
1053 template <int rank_, int dim, typename Number>
1054 template <typename OtherNumber>
1055 inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number>::
1056  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1057 {
1058  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
1059 }
1060 
1061 
1062 
1063 namespace internal
1064 {
1065  namespace TensorSubscriptor
1066  {
1067  template <typename ArrayElementType, int dim>
1068  inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV ArrayElementType &
1069  subscript(ArrayElementType * values,
1070  const unsigned int i,
1071  std::integral_constant<int, dim>)
1072  {
1073  // We cannot use Assert in a CUDA kernel
1074 #ifndef __CUDA_ARCH__
1075  Assert(i < dim, ExcIndexRange(i, 0, dim));
1076 #endif
1077  return values[i];
1078  }
1079 
1080 
1081  template <typename ArrayElementType>
1082  ArrayElementType &
1083  subscript(ArrayElementType *,
1084  const unsigned int,
1085  std::integral_constant<int, 0>)
1086  {
1087  Assert(
1088  false,
1089  ExcMessage(
1090  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1091  static ArrayElementType t;
1092  return t;
1093  }
1094  } // namespace TensorSubscriptor
1095 } // namespace internal
1096 
1097 
1098 template <int rank_, int dim, typename Number>
1099 inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1101  operator[](const unsigned int i)
1102 {
1103  return ::internal::TensorSubscriptor::subscript(
1104  values, i, std::integral_constant<int, dim>());
1105 }
1106 
1107 
1108 template <int rank_, int dim, typename Number>
1109 inline DEAL_II_ALWAYS_INLINE
1110  DEAL_II_CUDA_HOST_DEV const typename Tensor<rank_, dim, Number>::value_type &
1111  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1112 {
1113  return ::internal::TensorSubscriptor::subscript(
1114  values, i, std::integral_constant<int, dim>());
1115 }
1116 
1117 
1118 template <int rank_, int dim, typename Number>
1119 inline const Number &Tensor<rank_, dim, Number>::
1120  operator[](const TableIndices<rank_> &indices) const
1121 {
1122  Assert(dim != 0,
1123  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1124 
1125  return TensorAccessors::extract<rank_>(*this, indices);
1126 }
1127 
1128 
1129 
1130 template <int rank_, int dim, typename Number>
1131 inline Number &Tensor<rank_, dim, Number>::
1133 {
1134  Assert(dim != 0,
1135  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1136 
1137  return TensorAccessors::extract<rank_>(*this, indices);
1138 }
1139 
1140 
1141 
1142 template <int rank_, int dim, typename Number>
1143 inline Number *
1145 {
1146  return std::addressof(
1147  this->operator[](this->unrolled_to_component_indices(0)));
1148 }
1149 
1150 
1151 
1152 template <int rank_, int dim, typename Number>
1153 inline const Number *
1155 {
1156  return std::addressof(
1157  this->operator[](this->unrolled_to_component_indices(0)));
1158 }
1159 
1160 
1161 
1162 template <int rank_, int dim, typename Number>
1163 inline Number *
1165 {
1166  return begin_raw() + n_independent_components;
1167 }
1168 
1169 
1170 
1171 template <int rank_, int dim, typename Number>
1172 inline const Number *
1174 {
1175  return begin_raw() + n_independent_components;
1176 }
1177 
1178 
1179 
1180 template <int rank_, int dim, typename Number>
1181 template <typename OtherNumber>
1182 inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
1184 {
1185  if (dim > 0)
1186  std::copy(&t.values[0], &t.values[0] + dim, &values[0]);
1187  return *this;
1188 }
1189 
1190 
1191 template <int rank_, int dim, typename Number>
1192 inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
1194 {
1196  ExcMessage("Only assignment with zero is allowed"));
1197  (void)d;
1198 
1199  for (unsigned int i = 0; i < dim; ++i)
1200  values[i] = internal::NumberType<Number>::value(0.0);
1201  return *this;
1202 }
1203 
1204 
1205 template <int rank_, int dim, typename Number>
1206 template <typename OtherNumber>
1207 inline bool
1210 {
1211  for (unsigned int i = 0; i < dim; ++i)
1212  if (values[i] != p.values[i])
1213  return false;
1214  return true;
1215 }
1216 
1217 
1218 // At some places in the library, we have Point<0> for formal reasons
1219 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1220 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1221 // in the above function that the loop end check always fails, we
1222 // implement this function here
1223 template <>
1224 template <>
1225 inline bool
1227 {
1228  return true;
1229 }
1230 
1231 
1232 template <int rank_, int dim, typename Number>
1233 template <typename OtherNumber>
1234 inline bool
1237 {
1238  return !((*this) == p);
1239 }
1240 
1241 
1242 template <int rank_, int dim, typename Number>
1243 template <typename OtherNumber>
1246 {
1247  for (unsigned int i = 0; i < dim; ++i)
1248  values[i] += p.values[i];
1249  return *this;
1250 }
1251 
1252 
1253 template <int rank_, int dim, typename Number>
1254 template <typename OtherNumber>
1257 {
1258  for (unsigned int i = 0; i < dim; ++i)
1259  values[i] -= p.values[i];
1260  return *this;
1261 }
1262 
1263 
1264 template <int rank_, int dim, typename Number>
1265 template <typename OtherNumber>
1266 inline DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
1267 Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1268 {
1269  for (unsigned int i = 0; i < dim; ++i)
1270  values[i] *= s;
1271  return *this;
1272 }
1273 
1274 
1275 template <int rank_, int dim, typename Number>
1276 template <typename OtherNumber>
1278 Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1279 {
1280  for (unsigned int i = 0; i < dim; ++i)
1281  values[i] /= s;
1282  return *this;
1283 }
1284 
1285 
1286 template <int rank_, int dim, typename Number>
1289 {
1291 
1292  for (unsigned int i = 0; i < dim; ++i)
1293  tmp.values[i] = -values[i];
1294 
1295  return tmp;
1296 }
1297 
1298 
1299 template <int rank_, int dim, typename Number>
1302 {
1303  return std::sqrt(norm_square());
1304 }
1305 
1306 
1307 template <int rank_, int dim, typename Number>
1308 inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
1310 {
1312  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1313  for (unsigned int i = 0; i < dim; ++i)
1314  s += values[i].norm_square();
1315 
1316  return s;
1317 }
1318 
1319 
1320 template <int rank_, int dim, typename Number>
1321 template <typename OtherNumber>
1322 inline void
1323 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
1324 {
1325  AssertDimension(result.size(),
1326  (Utilities::fixed_power<rank_, unsigned int>(dim)));
1327 
1328  unsigned int index = 0;
1329  unroll_recursion(result, index);
1330 }
1331 
1332 
1333 template <int rank_, int dim, typename Number>
1334 template <typename OtherNumber>
1335 inline void
1337  unsigned int & index) const
1338 {
1339  for (unsigned int i = 0; i < dim; ++i)
1340  values[i].unroll_recursion(result, index);
1341 }
1342 
1343 
1344 template <int rank_, int dim, typename Number>
1345 inline unsigned int
1347  const TableIndices<rank_> &indices)
1348 {
1349  unsigned int index = 0;
1350  for (int r = 0; r < rank_; ++r)
1351  index = index * dim + indices[r];
1352 
1353  return index;
1354 }
1355 
1356 
1357 
1358 namespace internal
1359 {
1360  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1361  // and rank=2. Make sure we don't have compiler warnings.
1362 
1363  template <int dim>
1364  inline unsigned int
1365  mod(const unsigned int x)
1366  {
1367  return x % dim;
1368  }
1369 
1370  template <>
1371  inline unsigned int
1372  mod<0>(const unsigned int x)
1373  {
1374  Assert(false, ExcInternalError());
1375  return x;
1376  }
1377 
1378  template <int dim>
1379  inline unsigned int
1380  div(const unsigned int x)
1381  {
1382  return x / dim;
1383  }
1384 
1385  template <>
1386  inline unsigned int
1387  div<0>(const unsigned int x)
1388  {
1389  Assert(false, ExcInternalError());
1390  return x;
1391  }
1392 
1393 } // namespace internal
1394 
1395 
1396 
1397 template <int rank_, int dim, typename Number>
1398 inline TableIndices<rank_>
1400 {
1401  Assert(i < n_independent_components,
1402  ExcIndexRange(i, 0, n_independent_components));
1403 
1404  TableIndices<rank_> indices;
1405 
1406  unsigned int remainder = i;
1407  for (int r = rank_ - 1; r >= 0; --r)
1408  {
1409  indices[r] = internal::mod<dim>(remainder);
1410  remainder = internal::div<dim>(remainder);
1411  }
1412  Assert(remainder == 0, ExcInternalError());
1413 
1414  return indices;
1415 }
1416 
1417 
1418 template <int rank_, int dim, typename Number>
1419 inline void
1421 {
1422  for (unsigned int i = 0; i < dim; ++i)
1423  values[i] = internal::NumberType<Number>::value(0.0);
1424 }
1425 
1426 
1427 template <int rank_, int dim, typename Number>
1428 inline std::size_t
1430 {
1431  return sizeof(Tensor<rank_, dim, Number>);
1432 }
1433 
1434 
1435 template <int rank_, int dim, typename Number>
1436 template <class Archive>
1437 inline void
1438 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1439 {
1440  ar &values;
1441 }
1442 
1443 
1444 /* ----------------- Non-member functions operating on tensors. ------------ */
1445 
1450 
1458 template <int rank_, int dim, typename Number>
1459 inline std::ostream &
1460 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1461 {
1462  for (unsigned int i = 0; i < dim; ++i)
1463  {
1464  out << p[i];
1465  if (i != dim - 1)
1466  out << ' ';
1467  }
1468 
1469  return out;
1470 }
1471 
1472 
1479 template <int dim, typename Number>
1480 inline std::ostream &
1481 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1482 {
1483  out << static_cast<const Number &>(p);
1484  return out;
1485 }
1486 
1487 
1489 
1493 
1494 
1503 template <int dim, typename Number, typename Other>
1504 inline DEAL_II_ALWAYS_INLINE typename ProductType<Other, Number>::type
1505 operator*(const Other &object, const Tensor<0, dim, Number> &t)
1506 {
1507  return object * static_cast<const Number &>(t);
1508 }
1509 
1510 
1511 
1520 template <int dim, typename Number, typename Other>
1521 inline DEAL_II_ALWAYS_INLINE typename ProductType<Number, Other>::type
1522 operator*(const Tensor<0, dim, Number> &t, const Other &object)
1523 {
1524  return static_cast<const Number &>(t) * object;
1525 }
1526 
1527 
1537 template <int dim, typename Number, typename OtherNumber>
1538 inline DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
1540  const Tensor<0, dim, OtherNumber> &src2)
1541 {
1542  return static_cast<const Number &>(src1) *
1543  static_cast<const OtherNumber &>(src2);
1544 }
1545 
1546 
1552 template <int dim, typename Number, typename OtherNumber>
1553 inline DEAL_II_ALWAYS_INLINE
1554  Tensor<0,
1555  dim,
1556  typename ProductType<Number,
1557  typename EnableIfScalar<OtherNumber>::type>::type>
1558  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
1559 {
1560  return static_cast<const Number &>(t) / factor;
1561 }
1562 
1563 
1569 template <int dim, typename Number, typename OtherNumber>
1570 inline DEAL_II_ALWAYS_INLINE
1573  const Tensor<0, dim, OtherNumber> &q)
1574 {
1575  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1576 }
1577 
1578 
1584 template <int dim, typename Number, typename OtherNumber>
1585 inline DEAL_II_ALWAYS_INLINE
1588  const Tensor<0, dim, OtherNumber> &q)
1589 {
1590  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1591 }
1592 
1593 
1604 template <int rank, int dim, typename Number, typename OtherNumber>
1605 inline DEAL_II_ALWAYS_INLINE
1606  Tensor<rank,
1607  dim,
1608  typename ProductType<Number,
1609  typename EnableIfScalar<OtherNumber>::type>::type>
1610  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1611 {
1612  // recurse over the base objects
1614  for (unsigned int d = 0; d < dim; ++d)
1615  tt[d] = t[d] * factor;
1616  return tt;
1617 }
1618 
1619 
1630 template <int rank, int dim, typename Number, typename OtherNumber>
1631 inline DEAL_II_ALWAYS_INLINE
1632  Tensor<rank,
1633  dim,
1635  OtherNumber>::type>
1636  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
1637 {
1638  // simply forward to the operator above
1639  return t * factor;
1640 }
1641 
1642 
1650 template <int rank, int dim, typename Number, typename OtherNumber>
1651 inline Tensor<
1652  rank,
1653  dim,
1654  typename ProductType<Number,
1655  typename EnableIfScalar<OtherNumber>::type>::type>
1656 operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1657 {
1658  // recurse over the base objects
1660  for (unsigned int d = 0; d < dim; ++d)
1661  tt[d] = t[d] / factor;
1662  return tt;
1663 }
1664 
1665 
1673 template <int rank, int dim, typename Number, typename OtherNumber>
1674 inline DEAL_II_ALWAYS_INLINE
1678 {
1680 
1681  for (unsigned int i = 0; i < dim; ++i)
1682  tmp[i] += q[i];
1683 
1684  return tmp;
1685 }
1686 
1687 
1695 template <int rank, int dim, typename Number, typename OtherNumber>
1696 inline DEAL_II_ALWAYS_INLINE
1700 {
1702 
1703  for (unsigned int i = 0; i < dim; ++i)
1704  tmp[i] -= q[i];
1705 
1706  return tmp;
1707 }
1708 
1709 
1711 
1715 
1716 
1740 template <int rank_1,
1741  int rank_2,
1742  int dim,
1743  typename Number,
1744  typename OtherNumber>
1745 inline DEAL_II_ALWAYS_INLINE
1746  typename Tensor<rank_1 + rank_2 - 2,
1747  dim,
1748  typename ProductType<Number, OtherNumber>::type>::tensor_type
1751 {
1752  typename Tensor<rank_1 + rank_2 - 2,
1753  dim,
1754  typename ProductType<Number, OtherNumber>::type>::tensor_type
1755  result;
1756 
1757  TensorAccessors::internal::
1758  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
1759  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
1760  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
1761 
1762  return result;
1763 }
1764 
1765 
1795 template <int index_1,
1796  int index_2,
1797  int rank_1,
1798  int rank_2,
1799  int dim,
1800  typename Number,
1801  typename OtherNumber>
1802 inline DEAL_II_ALWAYS_INLINE
1803  typename Tensor<rank_1 + rank_2 - 2,
1804  dim,
1805  typename ProductType<Number, OtherNumber>::type>::tensor_type
1808 {
1809  Assert(0 <= index_1 && index_1 < rank_1,
1810  ExcMessage(
1811  "The specified index_1 must lie within the range [0,rank_1)"));
1812  Assert(0 <= index_2 && index_2 < rank_2,
1813  ExcMessage(
1814  "The specified index_2 must lie within the range [0,rank_2)"));
1815 
1816  using namespace TensorAccessors;
1817  using namespace TensorAccessors::internal;
1818 
1819  // Reorder index_1 to the end of src1:
1820  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>
1821  reord_01 = reordered_index_view<index_1, rank_1>(src1);
1822 
1823  // Reorder index_2 to the end of src2:
1824  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>
1825  reord_02 = reordered_index_view<index_2, rank_2>(src2);
1826 
1827  typename Tensor<rank_1 + rank_2 - 2,
1828  dim,
1829  typename ProductType<Number, OtherNumber>::type>::tensor_type
1830  result;
1831  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
1832  return result;
1833 }
1834 
1835 
1867 template <int index_1,
1868  int index_2,
1869  int index_3,
1870  int index_4,
1871  int rank_1,
1872  int rank_2,
1873  int dim,
1874  typename Number,
1875  typename OtherNumber>
1876 inline
1877  typename Tensor<rank_1 + rank_2 - 4,
1878  dim,
1879  typename ProductType<Number, OtherNumber>::type>::tensor_type
1882 {
1883  Assert(0 <= index_1 && index_1 < rank_1,
1884  ExcMessage(
1885  "The specified index_1 must lie within the range [0,rank_1)"));
1886  Assert(0 <= index_3 && index_3 < rank_1,
1887  ExcMessage(
1888  "The specified index_3 must lie within the range [0,rank_1)"));
1889  Assert(index_1 != index_3,
1890  ExcMessage("index_1 and index_3 must not be the same"));
1891  Assert(0 <= index_2 && index_2 < rank_2,
1892  ExcMessage(
1893  "The specified index_2 must lie within the range [0,rank_2)"));
1894  Assert(0 <= index_4 && index_4 < rank_2,
1895  ExcMessage(
1896  "The specified index_4 must lie within the range [0,rank_2)"));
1897  Assert(index_2 != index_4,
1898  ExcMessage("index_2 and index_4 must not be the same"));
1899 
1900  using namespace TensorAccessors;
1901  using namespace TensorAccessors::internal;
1902 
1903  // Reorder index_1 to the end of src1:
1904  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>
1905  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
1906 
1907  // Reorder index_2 to the end of src2:
1908  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>
1909  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
1910 
1911  // Now, reorder index_3 to the end of src1. We have to make sure to
1912  // preserve the orginial ordering: index_1 has been removed. If
1913  // index_3 > index_1, we have to use (index_3 - 1) instead:
1914  ReorderedIndexView<
1915  (index_3 < index_1 ? index_3 : index_3 - 1),
1916  rank_1,
1917  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
1918  reord_3 =
1919  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
1920  index_3 - 1,
1921  rank_1 > (reord_1);
1922 
1923  // Now, reorder index_4 to the end of src2. We have to make sure to
1924  // preserve the orginial ordering: index_2 has been removed. If
1925  // index_4 > index_2, we have to use (index_4 - 1) instead:
1926  ReorderedIndexView<
1927  (index_4 < index_2 ? index_4 : index_4 - 1),
1928  rank_2,
1929  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
1930  reord_4 =
1931  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
1932  index_4 - 1,
1933  rank_2 > (reord_2);
1934 
1935  typename Tensor<rank_1 + rank_2 - 4,
1936  dim,
1937  typename ProductType<Number, OtherNumber>::type>::tensor_type
1938  result;
1939  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
1940  return result;
1941 }
1942 
1943 
1957 template <int rank, int dim, typename Number, typename OtherNumber>
1958 inline DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
1960  const Tensor<rank, dim, OtherNumber> &right)
1961 {
1962  typename ProductType<Number, OtherNumber>::type result;
1963  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
1964  return result;
1965 }
1966 
1967 
1986 template <template <int, int, typename> class TensorT1,
1987  template <int, int, typename> class TensorT2,
1988  template <int, int, typename> class TensorT3,
1989  int rank_1,
1990  int rank_2,
1991  int dim,
1992  typename T1,
1993  typename T2,
1994  typename T3>
1996 contract3(const TensorT1<rank_1, dim, T1> & left,
1997  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
1998  const TensorT3<rank_2, dim, T3> & right)
1999 {
2000  using return_type =
2002  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2003  middle,
2004  right);
2005 }
2006 
2007 
2019 template <int rank_1,
2020  int rank_2,
2021  int dim,
2022  typename Number,
2023  typename OtherNumber>
2024 inline DEAL_II_ALWAYS_INLINE
2028 {
2029  typename Tensor<rank_1 + rank_2,
2030  dim,
2031  typename ProductType<Number, OtherNumber>::type>::tensor_type
2032  result;
2033  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2034  return result;
2035 }
2036 
2037 
2039 
2043 
2044 
2056 template <int dim, typename Number>
2057 inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
2059 {
2060  Assert(dim == 2, ExcInternalError());
2061 
2062  Tensor<1, dim, Number> result;
2063 
2064  result[0] = src[1];
2065  result[1] = -src[0];
2066 
2067  return result;
2068 }
2069 
2070 
2081 template <int dim, typename Number>
2082 inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
2084  const Tensor<1, dim, Number> &src2)
2085 {
2086  Assert(dim == 3, ExcInternalError());
2087 
2088  Tensor<1, dim, Number> result;
2089 
2090  result[0] = src1[1] * src2[2] - src1[2] * src2[1];
2091  result[1] = src1[2] * src2[0] - src1[0] * src2[2];
2092  result[2] = src1[0] * src2[1] - src1[1] * src2[0];
2093 
2094  return result;
2095 }
2096 
2097 
2099 
2103 
2104 
2111 template <int dim, typename Number>
2112 inline Number
2114 {
2115  // Compute the determinant using the Laplace expansion of the
2116  // determinant. We expand along the last row.
2117  Number det = internal::NumberType<Number>::value(0.0);
2118 
2119  for (unsigned int k = 0; k < dim; ++k)
2120  {
2121  Tensor<2, dim - 1, Number> minor;
2122  for (unsigned int i = 0; i < dim - 1; ++i)
2123  for (unsigned int j = 0; j < dim - 1; ++j)
2124  minor[i][j] = t[i][j < k ? j : j + 1];
2125 
2126  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2127 
2128  det += t[dim - 1][k] * cofactor;
2129  }
2130 
2131  return ((dim % 2 == 0) ? 1. : -1.) * det;
2132 }
2133 
2139 template <typename Number>
2140 inline Number
2142 {
2143  return t[0][0];
2144 }
2145 
2146 
2154 template <int dim, typename Number>
2155 inline DEAL_II_ALWAYS_INLINE Number
2157 {
2158  Number t = d[0][0];
2159  for (unsigned int i = 1; i < dim; ++i)
2160  t += d[i][i];
2161  return t;
2162 }
2163 
2164 
2174 template <int dim, typename Number>
2177 {
2178  Number return_tensor[dim][dim];
2179 
2180  // if desired, take over the
2181  // inversion of a 4x4 tensor
2182  // from the FullMatrix
2183  AssertThrow(false, ExcNotImplemented());
2184 
2185  return Tensor<2, dim, Number>(return_tensor);
2186 }
2187 
2188 
2189 #ifndef DOXYGEN
2190 
2191 template <typename Number>
2192 inline Tensor<2, 1, Number>
2193 invert(const Tensor<2, 1, Number> &t)
2194 {
2195  Number return_tensor[1][1];
2196 
2197  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2198 
2199  return Tensor<2, 1, Number>(return_tensor);
2200 }
2201 
2202 
2203 template <typename Number>
2204 inline Tensor<2, 2, Number>
2205 invert(const Tensor<2, 2, Number> &t)
2206 {
2207  Tensor<2, 2, Number> return_tensor;
2208 
2209  // this is Maple output,
2210  // thus a bit unstructured
2211  const Number inv_det_t = internal::NumberType<Number>::value(
2212  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2213  return_tensor[0][0] = t[1][1];
2214  return_tensor[0][1] = -t[0][1];
2215  return_tensor[1][0] = -t[1][0];
2216  return_tensor[1][1] = t[0][0];
2217  return_tensor *= inv_det_t;
2218 
2219  return return_tensor;
2220 }
2221 
2222 
2223 template <typename Number>
2224 inline Tensor<2, 3, Number>
2225 invert(const Tensor<2, 3, Number> &t)
2226 {
2227  Tensor<2, 3, Number> return_tensor;
2228 
2229  const Number t4 = internal::NumberType<Number>::value(t[0][0] * t[1][1]),
2230  t6 = internal::NumberType<Number>::value(t[0][0] * t[1][2]),
2231  t8 = internal::NumberType<Number>::value(t[0][1] * t[1][0]),
2232  t00 = internal::NumberType<Number>::value(t[0][2] * t[1][0]),
2233  t01 = internal::NumberType<Number>::value(t[0][1] * t[2][0]),
2234  t04 = internal::NumberType<Number>::value(t[0][2] * t[2][0]),
2236  1.0 / (t4 * t[2][2] - t6 * t[2][1] - t8 * t[2][2] +
2237  t00 * t[2][1] + t01 * t[1][2] - t04 * t[1][1]));
2238  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2239  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2240  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2241  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2242  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2243  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2244  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2245  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2246  return_tensor[1][1] =
2247  internal::NumberType<Number>::value(t[0][0] * t[2][2]) - t04;
2248  return_tensor[1][2] = t00 - t6;
2249  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2250  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2251  return_tensor[2][1] =
2252  t01 - internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2253  return_tensor[2][2] = internal::NumberType<Number>::value(t4 - t8);
2254  return_tensor *= inv_det_t;
2255 
2256  return return_tensor;
2257 }
2258 
2259 #endif /* DOXYGEN */
2260 
2261 
2268 template <int dim, typename Number>
2269 inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, Number>
2271 {
2273  for (unsigned int i = 0; i < dim; ++i)
2274  {
2275  tt[i][i] = t[i][i];
2276  for (unsigned int j = i + 1; j < dim; ++j)
2277  {
2278  tt[i][j] = t[j][i];
2279  tt[j][i] = t[i][j];
2280  };
2281  }
2282  return tt;
2283 }
2284 
2285 
2299 template <int dim, typename Number>
2302 {
2303  return determinant(t) * invert(t);
2304 }
2305 
2306 
2321 template <int dim, typename Number>
2324 {
2325  return transpose(adjugate(t));
2326 }
2327 
2328 
2336 template <int dim, typename Number>
2337 inline Number
2339 {
2340  Number max = internal::NumberType<Number>::value(0.0);
2341  for (unsigned int j = 0; j < dim; ++j)
2342  {
2343  Number sum = internal::NumberType<Number>::value(0.0);
2344  for (unsigned int i = 0; i < dim; ++i)
2345  sum += std::fabs(t[i][j]);
2346 
2347  if (sum > max)
2348  max = sum;
2349  }
2350 
2351  return max;
2352 }
2353 
2354 
2362 template <int dim, typename Number>
2363 inline Number
2365 {
2366  Number max = internal::NumberType<Number>::value(0.0);
2367  for (unsigned int i = 0; i < dim; ++i)
2368  {
2369  Number sum = internal::NumberType<Number>::value(0.0);
2370  for (unsigned int j = 0; j < dim; ++j)
2371  sum += std::fabs(t[i][j]);
2372 
2373  if (sum > max)
2374  max = sum;
2375  }
2376 
2377  return max;
2378 }
2379 
2381 
2382 
2383 #ifndef DOXYGEN
2384 
2385 
2386 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2387 
2388 // Specialization of functions for Adol-C number types when
2389 // the advanced branching feature is used
2390 template <int dim>
2391 inline adouble
2392 l1_norm(const Tensor<2, dim, adouble> &t)
2393 {
2394  adouble max = internal::NumberType<adouble>::value(0.0);
2395  for (unsigned int j = 0; j < dim; ++j)
2396  {
2397  adouble sum = internal::NumberType<adouble>::value(0.0);
2398  for (unsigned int i = 0; i < dim; ++i)
2399  sum += std::fabs(t[i][j]);
2400 
2401  condassign(max, (sum > max), sum, max);
2402  }
2403 
2404  return max;
2405 }
2406 
2407 
2408 template <int dim>
2409 inline adouble
2410 linfty_norm(const Tensor<2, dim, adouble> &t)
2411 {
2412  adouble max = internal::NumberType<adouble>::value(0.0);
2413  for (unsigned int i = 0; i < dim; ++i)
2414  {
2415  adouble sum = internal::NumberType<adouble>::value(0.0);
2416  for (unsigned int j = 0; j < dim; ++j)
2417  sum += std::fabs(t[i][j]);
2418 
2419  condassign(max, (sum > max), sum, max);
2420  }
2421 
2422  return max;
2423 }
2424 
2425 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2426 
2427 
2428 #endif // DOXYGEN
2429 
2430 DEAL_II_NAMESPACE_CLOSE
2431 
2432 // include deprecated non-member functions operating on Tensor
2433 #include <deal.II/base/tensor_deprecated.h>
2434 
2435 #endif
Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2156
#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:1346
Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1656
static const unsigned int n_independent_components
Definition: tensor.h:419
bool value_is_zero(const Number &value)
Definition: numbers.h:803
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2338
static std::size_t memory_consumption()
Definition: tensor.h:1429
Tensor< rank_-1, dim, Number > values[(dim!=0)?dim:1]
Definition: tensor.h:671
Tensor & operator/=(const OtherNumber &factor)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2364
Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2301
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:1806
typename Tensor< rank_-1, dim, VectorizedArray< Number > >::array_type[(dim!=0)?dim:1] array_type
Definition: tensor.h:434
static const unsigned int rank
Definition: tensor.h:413
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1301
#define AssertThrow(cond, exc)
Definition: exceptions.h:1329
static real_type abs(const number &x)
Definition: numbers.h:484
bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1236
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:787
void unroll(Vector< OtherNumber > &result) const
Definition: tensor.h:1323
Tensor & operator*=(const OtherNumber &factor)
Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2323
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:665
Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2270
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:1636
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:1996
Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:2176
static real_type abs_square(const number &x)
Definition: numbers.h:475
Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2058
bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1209
#define Assert(cond, exc)
Definition: exceptions.h:1239
ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:1959
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: tensor.h:1336
Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
Definition: tensor.h:1572
void serialize(Archive &ar, const unsigned int version)
Definition: tensor.h:1438
Tensor< 0, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const Tensor< 0, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1558
Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1610
Number * end_raw()
Definition: tensor.h:1164
typename Tensor< rank_-1, dim, VectorizedArray< Number > >::tensor_type value_type
Definition: tensor.h:427
value_type & operator[](const unsigned int i)
Definition: tensor.h:1101
Tensor()
Definition: tensor.h:1015
Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< rank, dim, Number > &p, const Tensor< rank, dim, OtherNumber > &q)
Definition: tensor.h:1676
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:1399
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:1880
Tensor< 1, dim, Number > cross_product_3d(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:2083
typename numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:125
Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2113
static const unsigned int dimension
Definition: tensor.h:408
Tensor operator-() const
Definition: tensor.h:1288
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:1698
ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1505
Number * begin_raw()
Definition: tensor.h:1144
static::ExceptionBase & ExcNotImplemented()
ProductType< Number, Other >::type operator*(const Tensor< 0, dim, Number > &t, const Other &object)
Definition: tensor.h:1522
ProductType< Number, OtherNumber >::type operator*(const Tensor< 0, dim, Number > &src1, const Tensor< 0, dim, OtherNumber > &src2)
Definition: tensor.h:1539
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:2026
numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1309
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:1587
void clear()
Definition: tensor.h:1420
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:1749
static::ExceptionBase & ExcInternalError()
Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:2141