Reference documentation for deal.II version Git 4c3a1f2 2018-02-23 12:56:49 +0100
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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_tensor_h
17 #define dealii_tensor_h
18 
19 #include <deal.II/base/config.h>
20 #include <deal.II/base/exceptions.h>
21 #include <deal.II/base/numbers.h>
22 #include <deal.II/base/table_indices.h>
23 #include <deal.II/base/tensor_accessors.h>
24 #include <deal.II/base/template_constraints.h>
25 #include <deal.II/base/utilities.h>
26 
27 #include <cmath>
28 #include <ostream>
29 #include <vector>
30 
31 
32 DEAL_II_NAMESPACE_OPEN
33 
34 // Forward declarations:
35 
36 template <int dim, typename Number> class Point;
37 template <int rank_, int dim, typename Number = double> class Tensor;
38 template <typename Number> class Vector;
39 template <typename Number> class VectorizedArray;
40 
41 #ifndef DOXYGEN
42 // Overload invalid tensor types of negative rank that come up during
43 // overload resolution of operator* and related contraction variants.
44 template <int dim, typename Number>
45 class Tensor<-2, dim, Number>
46 {
47 };
48 
49 template <int dim, typename Number>
50 class Tensor<-1, dim, Number>
51 {
52 };
53 #endif /* DOXYGEN */
54 
55 
86 template <int dim, typename Number>
87 class Tensor<0,dim,Number>
88 {
89 public:
98  static const unsigned int dimension = dim;
99 
103  static const unsigned int rank = 0;
104 
108  static const unsigned int n_independent_components = 1;
109 
119 
124  typedef Number value_type;
125 
131  typedef Number array_type;
132 
138  DEAL_II_CUDA_HOST_DEV Tensor ();
139 
145  template <typename OtherNumber>
146  Tensor (const Tensor<0,dim,OtherNumber> &initializer);
147 
151  template <typename OtherNumber>
152  Tensor (const OtherNumber &initializer);
153 
157  Number *
158  begin_raw();
159 
163  const Number *
164  begin_raw() const;
165 
169  Number *
170  end_raw();
171 
176  const Number *
177  end_raw() const;
178 
188  DEAL_II_CUDA_HOST_DEV operator Number &();
189 
198  DEAL_II_CUDA_HOST_DEV operator const Number &() const;
199 
205  template <typename OtherNumber>
207 
212 
217  template <typename OtherNumber>
218  Tensor &operator = (const OtherNumber &d);
219 
223  template <typename OtherNumber>
224  bool operator == (const Tensor<0,dim,OtherNumber> &rhs) const;
225 
229  template <typename OtherNumber>
230  bool operator != (const Tensor<0,dim,OtherNumber> &rhs) const;
231 
235  template <typename OtherNumber>
237 
241  template <typename OtherNumber>
243 
249  template <typename OtherNumber>
250  DEAL_II_CUDA_HOST_DEV Tensor &operator *= (const OtherNumber &factor);
251 
255  template <typename OtherNumber>
256  Tensor &operator /= (const OtherNumber &factor);
257 
261  Tensor operator - () const;
262 
275  void clear ();
276 
282  real_type norm () const;
283 
290  DEAL_II_CUDA_HOST_DEV real_type norm_square () const;
291 
296  template <class Archive>
297  void serialize(Archive &ar, const unsigned int version);
298 
303  typedef Number tensor_type;
304 
305 private:
309  Number value;
310 
314  template <typename OtherNumber>
316  unsigned int &start_index) const;
317 
321  template <int, int, typename> friend class Tensor;
322 };
323 
324 
325 
367 template <int rank_, int dim, typename Number>
368 class Tensor
369 {
370 public:
379  static const unsigned int dimension = dim;
380 
384  static const unsigned int rank = rank_;
385 
390  static const unsigned int
392 
398  typedef typename Tensor<rank_-1,dim,Number>::tensor_type value_type;
399 
404  typedef typename Tensor<rank_-1,dim,Number>::array_type
405  array_type[(dim != 0) ? dim : 1];
406  // ... avoid a compiler warning in case of dim == 0 and ensure that the
407  // array always has positive size.
408 
414  DEAL_II_CUDA_HOST_DEV Tensor ();
415 
419  explicit Tensor (const array_type &initializer);
420 
426  template <typename OtherNumber>
427  Tensor (const Tensor<rank_,dim,OtherNumber> &initializer);
428 
432  template <typename OtherNumber>
433  Tensor (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer);
434 
438  template <typename OtherNumber>
439  operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const;
440 
446  DEAL_II_CUDA_HOST_DEV value_type &operator [] (const unsigned int i);
447 
453  DEAL_II_CUDA_HOST_DEV const value_type &operator[](const unsigned int i) const;
454 
458  const Number &operator [] (const TableIndices<rank_> &indices) const;
459 
463  Number &operator [] (const TableIndices<rank_> &indices);
464 
468  Number *
469  begin_raw();
470 
474  const Number *
475  begin_raw() const;
476 
480  Number *
481  end_raw();
482 
486  const Number *
487  end_raw() const;
488 
494  template <typename OtherNumber>
496 
503  Tensor &operator = (const Number &d);
504 
508  template <typename OtherNumber>
509  bool operator == (const Tensor<rank_,dim,OtherNumber> &) const;
510 
514  template <typename OtherNumber>
515  bool operator != (const Tensor<rank_,dim,OtherNumber> &) const;
516 
520  template <typename OtherNumber>
522 
526  template <typename OtherNumber>
528 
535  template <typename OtherNumber>
536  DEAL_II_CUDA_HOST_DEV Tensor &operator *= (const OtherNumber &factor);
537 
541  template <typename OtherNumber>
542  Tensor &operator /= (const OtherNumber &factor);
543 
547  Tensor operator - () const;
548 
561  void clear ();
562 
570 
577  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type norm_square() const;
578 
586  template <typename OtherNumber>
587  void unroll (Vector<OtherNumber> &result) const;
588 
593  static
594  unsigned int
596 
601  static
603 
608  static std::size_t memory_consumption ();
609 
614  template <class Archive>
615  void serialize(Archive &ar, const unsigned int version);
616 
622 
623 private:
627  Tensor<rank_-1, dim, Number> values[(dim != 0) ? dim : 1];
628  // ... avoid a compiler warning in case of dim == 0 and ensure that the
629  // array always has positive size.
630 
634  template <typename OtherNumber>
636  unsigned int &start_index) const;
637 
641  template <int, int, typename> friend class Tensor;
642 
647  friend class Point<dim,Number>;
648 };
649 
650 
651 namespace internal
652 {
663  template <int rank, int dim, typename T>
664  struct NumberType<Tensor<rank,dim,T> >
665  {
666  static const Tensor<rank,dim,T> &value (const Tensor<rank,dim,T> &t)
667  {
668  return t;
669  }
670 
671  static Tensor<rank,dim,T> value (const T &t)
672  {
673  Tensor<rank,dim,T> tmp;
674  tmp=t;
675  return tmp;
676  }
677  };
678 
679  template <int rank, int dim, typename T>
680  struct NumberType<Tensor<rank,dim,VectorizedArray<T> > >
681  {
682  static const Tensor<rank,dim,VectorizedArray<T> > &value (const Tensor<rank,dim,VectorizedArray<T> > &t)
683  {
684  return t;
685  }
686 
687  static Tensor<rank,dim,VectorizedArray<T> > value (const T &t)
688  {
691  return tmp;
692  }
693 
695  {
697  tmp=t;
698  return tmp;
699  }
700  };
701 }
702 
703 
704 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
705 
706 
707 template <int dim,typename Number>
708 inline
709 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::Tensor ()
710 // Some auto-differentiable numbers need explicit
711 // zero initialization.
712  : value(internal::NumberType<Number>::value(0.0))
713 {
714 }
715 
716 
717 
718 template <int dim, typename Number>
719 template <typename OtherNumber>
720 inline
721 Tensor<0,dim,Number>::Tensor (const OtherNumber &initializer)
722 {
723  value = internal::NumberType<Number>::value(initializer);
724 }
725 
726 
727 
728 template <int dim, typename Number>
729 template <typename OtherNumber>
730 inline
732 {
733  value = p.value;
734 }
735 
736 
737 
738 template <int dim, typename Number>
739 inline
740 Number *
742 {
743  return std::addressof(value);
744 }
745 
746 
747 
748 template <int dim, typename Number>
749 inline
750 const Number *
752 {
753  return std::addressof(value);
754 }
755 
756 
757 
758 template <int dim, typename Number>
759 inline
760 Number *
762 {
764 }
765 
766 
767 
768 template <int dim, typename Number>
769 inline
770 const Number *
772 {
774 }
775 
776 
777 
778 template <int dim, typename Number>
779 inline
780 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::operator Number &()
781 {
782  // We cannot use Assert inside a CUDA kernel
783 #ifndef DEAL_II_WITH_CUDA
784  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
785 #endif
786  return value;
787 }
788 
789 
790 template <int dim, typename Number>
791 inline
792 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::operator const Number &() const
793 {
794  // We cannot use Assert inside a CUDA kernel
795 #ifndef DEAL_II_WITH_CUDA
796  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
797 #endif
798  return value;
799 }
800 
801 
802 template <int dim, typename Number>
803 template <typename OtherNumber>
804 inline
806 {
808  return *this;
809 }
810 
811 
812 template <int dim, typename Number>
813 inline
815 {
816  value = p.value;
817  return *this;
818 }
819 
820 
821 template <int dim, typename Number>
822 template <typename OtherNumber>
823 inline
825 {
827  return *this;
828 }
829 
830 
831 template <int dim, typename Number>
832 template <typename OtherNumber>
833 inline
835 {
836  return (value == p.value);
837 }
838 
839 
840 template <int dim, typename Number>
841 template <typename OtherNumber>
842 inline
844 {
845  return !((*this) == p);
846 }
847 
848 
849 template <int dim, typename Number>
850 template <typename OtherNumber>
851 inline
853 {
854  value += p.value;
855  return *this;
856 }
857 
858 
859 template <int dim, typename Number>
860 template <typename OtherNumber>
861 inline
863 {
864  value -= p.value;
865  return *this;
866 }
867 
868 
869 template <int dim, typename Number>
870 template <typename OtherNumber>
871 inline
872 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator *= (const OtherNumber &s)
873 {
874  value *= s;
875  return *this;
876 }
877 
878 
879 template <int dim, typename Number>
880 template <typename OtherNumber>
881 inline
883 {
884  value /= s;
885  return *this;
886 }
887 
888 
889 template <int dim, typename Number>
890 inline
892 {
893  return -value;
894 }
895 
896 
897 template <int dim, typename Number>
898 inline
901 {
902  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
903  return numbers::NumberTraits<Number>::abs (value);
904 }
905 
906 
907 template <int dim, typename Number>
908 inline
910 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::norm_square () const
911 {
912  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
914 }
915 
916 
917 template <int dim, typename Number>
918 template <typename OtherNumber>
919 inline
920 void
922  unsigned int &index) const
923 {
924  Assert(dim != 0, ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
925  result[index] = value;
926  ++index;
927 }
928 
929 
930 template <int dim, typename Number>
931 inline
933 {
934  // Some auto-differentiable numbers need explicit
935  // zero initialization.
937 }
938 
939 
940 template <int dim, typename Number>
941 template <class Archive>
942 inline
943 void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int)
944 {
945  ar &value;
946 }
947 
948 
949 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
950 
951 
952 template <int rank_, int dim, typename Number>
953 inline
954 DEAL_II_CUDA_HOST_DEV Tensor<rank_,dim,Number>::Tensor ()
955 {
956  // All members of the c-style array values are already default initialized
957  // and thus all values are already set to zero recursively.
958 }
959 
960 
961 template <int rank_, int dim, typename Number>
962 inline
963 Tensor<rank_,dim,Number>::Tensor (const array_type &initializer)
964 {
965  for (unsigned int i=0; i<dim; ++i)
966  values[i] = Tensor<rank_-1, dim, Number>(initializer[i]);
967 }
968 
969 
970 template <int rank_, int dim, typename Number>
971 template <typename OtherNumber>
972 inline
974 {
975  for (unsigned int i=0; i!=dim; ++i)
976  values[i] = Tensor<rank_-1,dim,Number>(initializer[i]);
977 }
978 
979 
980 template <int rank_, int dim, typename Number>
981 template <typename OtherNumber>
982 inline
984 (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer)
985 {
986  for (unsigned int i=0; i<dim; ++i)
987  values[i] = initializer[i];
988 }
989 
990 
991 template <int rank_, int dim, typename Number>
992 template <typename OtherNumber>
993 inline
995 operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const
996 {
997  return Tensor<1,dim,Tensor<rank_-1,dim,Number> > (values);
998 }
999 
1000 
1001 
1002 namespace internal
1003 {
1004  namespace TensorSubscriptor
1005  {
1006  template <typename ArrayElementType, int dim>
1007  inline DEAL_II_ALWAYS_INLINE
1008  DEAL_II_CUDA_HOST_DEV
1009  ArrayElementType &
1010  subscript (ArrayElementType *values,
1011  const unsigned int i,
1012  std::integral_constant<int, dim>)
1013  {
1014  // We cannot use Assert in a CUDA kernel
1015 #ifndef DEAL_II_WITH_CUDA
1016  Assert (i<dim, ExcIndexRange(i, 0, dim));
1017 #endif
1018  return values[i];
1019  }
1020 
1021 
1022  template <typename ArrayElementType>
1023  DEAL_II_CUDA_HOST_DEV
1024  ArrayElementType &
1025  subscript (ArrayElementType *,
1026  const unsigned int,
1027  std::integral_constant<int, 0>)
1028  {
1029  // We cannot use Assert in a CUDA kernel
1030 #ifndef DEAL_II_WITH_CUDA
1031  Assert(false, ExcMessage("Cannot access elements of an object of type Tensor<rank,0,Number>."));
1032 #endif
1033  static ArrayElementType t;
1034  return t;
1035  }
1036  }
1037 }
1038 
1039 
1040 template <int rank_, int dim, typename Number>
1041 inline DEAL_II_ALWAYS_INLINE
1042 DEAL_II_CUDA_HOST_DEV
1045 {
1046  return ::internal::TensorSubscriptor::subscript(values, i, std::integral_constant<int, dim>());
1047 }
1048 
1049 
1050 template <int rank_, int dim, typename Number>
1051 inline DEAL_II_ALWAYS_INLINE
1052 DEAL_II_CUDA_HOST_DEV
1053 const typename Tensor<rank_,dim,Number>::value_type &
1054 Tensor<rank_,dim,Number>::operator[] (const unsigned int i) const
1055 {
1056  return ::internal::TensorSubscriptor::subscript(values, i, std::integral_constant<int, dim>());
1057 }
1058 
1059 
1060 template <int rank_, int dim, typename Number>
1061 inline
1062 const Number &
1064 {
1065  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1066 
1067  return TensorAccessors::extract<rank_>(*this, indices);
1068 }
1069 
1070 
1071 
1072 template <int rank_, int dim, typename Number>
1073 inline
1074 Number &
1076 {
1077  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1078 
1079  return TensorAccessors::extract<rank_>(*this, indices);
1080 }
1081 
1082 
1083 
1084 template <int rank_, int dim, typename Number>
1085 inline
1086 Number *
1088 {
1089  return std::addressof(this->operator[](this->unrolled_to_component_indices(0)));
1090 }
1091 
1092 
1093 
1094 template <int rank_, int dim, typename Number>
1095 inline
1096 const Number *
1098 {
1099  return std::addressof(this->operator[](this->unrolled_to_component_indices(0)));
1100 }
1101 
1102 
1103 
1104 template <int rank_, int dim, typename Number>
1105 inline
1106 Number *
1108 {
1109  return begin_raw()+n_independent_components;
1110 }
1111 
1112 
1113 
1114 template <int rank_, int dim, typename Number>
1115 inline
1116 const Number *
1118 {
1119  return begin_raw()+n_independent_components;
1120 }
1121 
1122 
1123 
1124 template <int rank_, int dim, typename Number>
1125 template <typename OtherNumber>
1126 inline
1129 {
1130  if (dim > 0)
1131  std::copy (&t.values[0], &t.values[0]+dim, &values[0]);
1132  return *this;
1133 }
1134 
1135 
1136 template <int rank_, int dim, typename Number>
1137 inline
1140 {
1142  ExcMessage ("Only assignment with zero is allowed"));
1143  (void) d;
1144 
1145  for (unsigned int i=0; i<dim; ++i)
1146  values[i] = internal::NumberType<Number>::value(0.0);
1147  return *this;
1148 }
1149 
1150 
1151 template <int rank_, int dim, typename Number>
1152 template <typename OtherNumber>
1153 inline
1154 bool
1156 {
1157  for (unsigned int i=0; i<dim; ++i)
1158  if (values[i] != p.values[i])
1159  return false;
1160  return true;
1161 }
1162 
1163 
1164 // At some places in the library, we have Point<0> for formal reasons
1165 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1166 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1167 // in the above function that the loop end check always fails, we
1168 // implement this function here
1169 template <>
1170 template <>
1171 inline
1173 {
1174  return true;
1175 }
1176 
1177 
1178 template <int rank_, int dim, typename Number>
1179 template <typename OtherNumber>
1180 inline
1181 bool
1183 {
1184  return !((*this) == p);
1185 }
1186 
1187 
1188 template <int rank_, int dim, typename Number>
1189 template <typename OtherNumber>
1190 inline
1193 {
1194  for (unsigned int i=0; i<dim; ++i)
1195  values[i] += p.values[i];
1196  return *this;
1197 }
1198 
1199 
1200 template <int rank_, int dim, typename Number>
1201 template <typename OtherNumber>
1202 inline
1205 {
1206  for (unsigned int i=0; i<dim; ++i)
1207  values[i] -= p.values[i];
1208  return *this;
1209 }
1210 
1211 
1212 template <int rank_, int dim, typename Number>
1213 template <typename OtherNumber>
1214 inline
1215 DEAL_II_CUDA_HOST_DEV
1217 Tensor<rank_,dim,Number>::operator *= (const OtherNumber &s)
1218 {
1219  for (unsigned int i=0; i<dim; ++i)
1220  values[i] *= s;
1221  return *this;
1222 }
1223 
1224 
1225 template <int rank_, int dim, typename Number>
1226 template <typename OtherNumber>
1227 inline
1229 Tensor<rank_,dim,Number>::operator /= (const OtherNumber &s)
1230 {
1231  for (unsigned int i=0; i<dim; ++i)
1232  values[i] /= s;
1233  return *this;
1234 }
1235 
1236 
1237 template <int rank_, int dim, typename Number>
1238 inline
1241 {
1243 
1244  for (unsigned int i=0; i<dim; ++i)
1245  tmp.values[i] = -values[i];
1246 
1247  return tmp;
1248 }
1249 
1250 
1251 template <int rank_, int dim, typename Number>
1252 inline
1255 {
1256  return std::sqrt (norm_square());
1257 }
1258 
1259 
1260 template <int rank_, int dim, typename Number>
1261 inline
1262 DEAL_II_CUDA_HOST_DEV
1265 {
1268  for (unsigned int i=0; i<dim; ++i)
1269  s += values[i].norm_square();
1270 
1271  return s;
1272 }
1273 
1274 
1275 template <int rank_, int dim, typename Number>
1276 template <typename OtherNumber>
1277 inline
1278 void
1280 {
1281  AssertDimension (result.size(),(Utilities::fixed_power<rank_, unsigned int>(dim)));
1282 
1283  unsigned int index = 0;
1284  unroll_recursion (result, index);
1285 }
1286 
1287 
1288 template <int rank_, int dim, typename Number>
1289 template <typename OtherNumber>
1290 inline
1291 void
1293  unsigned int &index) const
1294 {
1295  for (unsigned int i=0; i<dim; ++i)
1296  values[i].unroll_recursion(result, index);
1297 }
1298 
1299 
1300 template <int rank_, int dim, typename Number>
1301 inline
1302 unsigned int
1304 {
1305  unsigned int index = 0;
1306  for (int r = 0; r < rank_; ++r)
1307  index = index * dim + indices[r];
1308 
1309  return index;
1310 }
1311 
1312 
1313 template <int rank_, int dim, typename Number>
1314 inline
1317 {
1318  Assert (i < n_independent_components,
1319  ExcIndexRange (i, 0, n_independent_components));
1320 
1321  TableIndices<rank_> indices;
1322 
1323  unsigned int remainder = i;
1324  for (int r=rank_-1; r>=0; --r)
1325  {
1326  indices[r] = (remainder % dim);
1327  remainder /= dim;
1328  }
1329  Assert (remainder == 0, ExcInternalError());
1330 
1331  return indices;
1332 }
1333 
1334 
1335 template <int rank_, int dim, typename Number>
1336 inline
1338 {
1339  for (unsigned int i=0; i<dim; ++i)
1340  values[i] = internal::NumberType<Number>::value(0.0);
1341 }
1342 
1343 
1344 template <int rank_, int dim, typename Number>
1345 inline
1346 std::size_t
1348 {
1349  return sizeof(Tensor<rank_,dim,Number>);
1350 }
1351 
1352 
1353 template <int rank_, int dim, typename Number>
1354 template <class Archive>
1355 inline
1356 void
1357 Tensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
1358 {
1359  ar &values;
1360 }
1361 
1362 
1363 /* ----------------- Non-member functions operating on tensors. ------------ */
1364 
1369 
1377 template <int rank_, int dim, typename Number>
1378 inline
1379 std::ostream &operator << (std::ostream &out, const Tensor<rank_,dim,Number> &p)
1380 {
1381  for (unsigned int i = 0; i < dim; ++i)
1382  {
1383  out << p[i];
1384  if (i != dim - 1)
1385  out << ' ';
1386  }
1387 
1388  return out;
1389 }
1390 
1391 
1398 template <int dim, typename Number>
1399 inline
1400 std::ostream &operator << (std::ostream &out, const Tensor<0,dim,Number> &p)
1401 {
1402  out << static_cast<const Number &>(p);
1403  return out;
1404 }
1405 
1406 
1408 
1412 
1413 
1422 template <int dim, typename Number, typename Other>
1423 inline
1424 typename ProductType<Other, Number>::type
1425 operator * (const Other &object,
1426  const Tensor<0,dim,Number> &t)
1427 {
1428  return object * static_cast<const Number &>(t);
1429 }
1430 
1431 
1432 
1441 template <int dim, typename Number, typename Other>
1442 inline
1443 typename ProductType<Number, Other>::type
1445  const Other &object)
1446 {
1447  return static_cast<const Number &>(t) * object;
1448 }
1449 
1450 
1460 template <int dim, typename Number, typename OtherNumber>
1461 inline
1462 typename ProductType<Number, OtherNumber>::type
1464  const Tensor<0, dim, OtherNumber> &src2)
1465 {
1466  return static_cast<const Number &>(src1) *
1467  static_cast<const OtherNumber &>(src2);
1468 }
1469 
1470 
1476 template <int dim, typename Number, typename OtherNumber>
1477 inline
1480  const OtherNumber &factor)
1481 {
1482  return static_cast<const Number &>(t) / factor;
1483 }
1484 
1485 
1491 template <int dim, typename Number, typename OtherNumber>
1492 inline
1495  const Tensor<0,dim,OtherNumber> &q)
1496 {
1497  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1498 }
1499 
1500 
1506 template <int dim, typename Number, typename OtherNumber>
1507 inline
1510  const Tensor<0,dim,OtherNumber> &q)
1511 {
1512  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1513 }
1514 
1515 
1526 template <int rank, int dim,
1527  typename Number,
1528  typename OtherNumber>
1529 inline
1532  const OtherNumber &factor)
1533 {
1534  // recurse over the base objects
1536  for (unsigned int d=0; d<dim; ++d)
1537  tt[d] = t[d] * factor;
1538  return tt;
1539 }
1540 
1541 
1552 template <int rank, int dim,
1553  typename Number,
1554  typename OtherNumber>
1555 inline
1557 operator * (const Number &factor,
1559 {
1560  // simply forward to the operator above
1561  return t * factor;
1562 }
1563 
1564 
1572 template <int rank, int dim,
1573  typename Number,
1574  typename OtherNumber>
1575 inline
1578  const OtherNumber &factor)
1579 {
1580  // recurse over the base objects
1582  for (unsigned int d=0; d<dim; ++d)
1583  tt[d] = t[d] / factor;
1584  return tt;
1585 }
1586 
1587 
1595 template <int rank, int dim, typename Number, typename OtherNumber>
1596 inline
1600 {
1602 
1603  for (unsigned int i=0; i<dim; ++i)
1604  tmp[i] += q[i];
1605 
1606  return tmp;
1607 }
1608 
1609 
1617 template <int rank, int dim, typename Number, typename OtherNumber>
1618 inline
1622 {
1624 
1625  for (unsigned int i=0; i<dim; ++i)
1626  tmp[i] -= q[i];
1627 
1628  return tmp;
1629 }
1630 
1631 
1633 
1637 
1638 
1662 template <int rank_1, int rank_2, int dim,
1663  typename Number, typename OtherNumber>
1664 inline DEAL_II_ALWAYS_INLINE
1665 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1668 {
1669  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
1670 
1671  TensorAccessors::internal::ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1672  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
1673  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
1674 
1675  return result;
1676 }
1677 
1678 
1708 template <int index_1, int index_2,
1709  int rank_1, int rank_2, int dim,
1710  typename Number, typename OtherNumber>
1711 inline
1712 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1715 {
1716  Assert(0 <= index_1 && index_1 < rank_1,
1717  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1718  Assert(0 <= index_2 && index_2 < rank_2,
1719  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1720 
1721  using namespace TensorAccessors;
1722  using namespace TensorAccessors::internal;
1723 
1724  // Reorder index_1 to the end of src1:
1725  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1726  reord_01 = reordered_index_view<index_1, rank_1>(src1);
1727 
1728  // Reorder index_2 to the end of src2:
1729  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1730  reord_02 = reordered_index_view<index_2, rank_2>(src2);
1731 
1732  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1733  result;
1734  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
1735  return result;
1736 }
1737 
1738 
1770 template <int index_1, int index_2, int index_3, int index_4,
1771  int rank_1, int rank_2, int dim,
1772  typename Number, typename OtherNumber>
1773 inline
1774 typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1777 {
1778  Assert(0 <= index_1 && index_1 < rank_1,
1779  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1780  Assert(0 <= index_3 && index_3 < rank_1,
1781  ExcMessage("The specified index_3 must lie within the range [0,rank_1)"));
1782  Assert(index_1 != index_3,
1783  ExcMessage("index_1 and index_3 must not be the same"));
1784  Assert(0 <= index_2 && index_2 < rank_2,
1785  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1786  Assert(0 <= index_4 && index_4 < rank_2,
1787  ExcMessage("The specified index_4 must lie within the range [0,rank_2)"));
1788  Assert(index_2 != index_4,
1789  ExcMessage("index_2 and index_4 must not be the same"));
1790 
1791  using namespace TensorAccessors;
1792  using namespace TensorAccessors::internal;
1793 
1794  // Reorder index_1 to the end of src1:
1795  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1796  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
1797 
1798  // Reorder index_2 to the end of src2:
1799  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1800  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
1801 
1802  // Now, reorder index_3 to the end of src1. We have to make sure to
1803  // preserve the orginial ordering: index_1 has been removed. If
1804  // index_3 > index_1, we have to use (index_3 - 1) instead:
1805  ReorderedIndexView<(index_3 < index_1 ? index_3 : index_3 - 1), rank_1, ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> > >
1806  reord_3 = TensorAccessors::reordered_index_view<index_3 < index_1 ? index_3 : index_3 - 1, rank_1>(reord_1);
1807 
1808  // Now, reorder index_4 to the end of src2. We have to make sure to
1809  // preserve the orginial ordering: index_2 has been removed. If
1810  // index_4 > index_2, we have to use (index_4 - 1) instead:
1811  ReorderedIndexView<(index_4 < index_2 ? index_4 : index_4 - 1), rank_2, ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> > >
1812  reord_4 = TensorAccessors::reordered_index_view<index_4 < index_2 ? index_4 : index_4 - 1, rank_2>(reord_2);
1813 
1814  typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1815  result;
1816  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
1817  return result;
1818 }
1819 
1820 
1834 template <int rank, int dim, typename Number, typename OtherNumber>
1835 inline
1836 typename ProductType<Number, OtherNumber>::type
1838  const Tensor<rank, dim, OtherNumber> &right)
1839 {
1840  typename ProductType<Number, OtherNumber>::type result;
1841  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
1842  return result;
1843 }
1844 
1845 
1864 template <template <int, int, typename> class TensorT1,
1865  template <int, int, typename> class TensorT2,
1866  template <int, int, typename> class TensorT3,
1867  int rank_1, int rank_2, int dim,
1868  typename T1, typename T2, typename T3>
1870 contract3 (const TensorT1<rank_1, dim, T1> &left,
1871  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
1872  const TensorT3<rank_2, dim, T3> &right)
1873 {
1875  return_type;
1876  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(
1877  left, middle, right);
1878 }
1879 
1880 
1892 template <int rank_1, int rank_2, int dim,
1893  typename Number, typename OtherNumber>
1894 inline
1898 {
1900  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
1901  return result;
1902 }
1903 
1904 
1906 
1910 
1911 
1923 template <int dim, typename Number>
1924 inline
1927 {
1928  Assert (dim==2, ExcInternalError());
1929 
1930  Tensor<1, dim, Number> result;
1931 
1932  result[0] = src[1];
1933  result[1] = -src[0];
1934 
1935  return result;
1936 }
1937 
1938 
1949 template <int dim, typename Number>
1950 inline
1953  const Tensor<1,dim,Number> &src2)
1954 {
1955  Assert (dim==3, ExcInternalError());
1956 
1957  Tensor<1, dim, Number> result;
1958 
1959  result[0] = src1[1]*src2[2] - src1[2]*src2[1];
1960  result[1] = src1[2]*src2[0] - src1[0]*src2[2];
1961  result[2] = src1[0]*src2[1] - src1[1]*src2[0];
1962 
1963  return result;
1964 }
1965 
1966 
1968 
1972 
1973 
1980 template <int dim, typename Number>
1981 inline
1983 {
1984  // Compute the determinant using the Laplace expansion of the
1985  // determinant. We expand along the last row.
1986  Number det = internal::NumberType<Number>::value(0.0);
1987 
1988  for (unsigned int k=0; k<dim; ++k)
1989  {
1990  Tensor<2,dim-1,Number> minor;
1991  for (unsigned int i=0; i<dim-1; ++i)
1992  for (unsigned int j=0; j<dim-1; ++j)
1993  minor[i][j] = t[i][j<k ? j : j+1];
1994 
1995  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
1996 
1997  det += t[dim-1][k] * cofactor;
1998  }
1999 
2000  return ((dim % 2 == 0) ? 1. : -1.) * det;
2001 }
2002 
2008 template <typename Number>
2009 inline
2011 {
2012  return t[0][0];
2013 }
2014 
2015 
2023 template <int dim, typename Number>
2024 Number trace (const Tensor<2,dim,Number> &d)
2025 {
2026  Number t=d[0][0];
2027  for (unsigned int i=1; i<dim; ++i)
2028  t += d[i][i];
2029  return t;
2030 }
2031 
2032 
2042 template <int dim, typename Number>
2043 inline
2046 {
2047  Number return_tensor [dim][dim];
2048 
2049  // if desired, take over the
2050  // inversion of a 4x4 tensor
2051  // from the FullMatrix
2052  AssertThrow (false, ExcNotImplemented());
2053 
2054  return Tensor<2,dim,Number>(return_tensor);
2055 }
2056 
2057 
2058 #ifndef DOXYGEN
2059 
2060 template <typename Number>
2061 inline
2063 invert (const Tensor<2,1,Number> &t)
2064 {
2065  Number return_tensor [1][1];
2066 
2067  return_tensor[0][0] = 1.0/t[0][0];
2068 
2069  return Tensor<2,1,Number>(return_tensor);
2070 }
2071 
2072 
2073 template <typename Number>
2074 inline
2076 invert (const Tensor<2,2,Number> &t)
2077 {
2078  Tensor<2,2,Number> return_tensor;
2079 
2080  // this is Maple output,
2081  // thus a bit unstructured
2082  const Number inv_det_t = 1.0/(t[0][0]*t[1][1]-t[1][0]*t[0][1]);
2083  return_tensor[0][0] = t[1][1];
2084  return_tensor[0][1] = -t[0][1];
2085  return_tensor[1][0] = -t[1][0];
2086  return_tensor[1][1] = t[0][0];
2087  return_tensor *= inv_det_t;
2088 
2089  return return_tensor;
2090 }
2091 
2092 
2093 template <typename Number>
2094 inline
2096 invert (const Tensor<2,3,Number> &t)
2097 {
2098  Tensor<2,3,Number> return_tensor;
2099 
2100  const Number t4 = t[0][0]*t[1][1],
2101  t6 = t[0][0]*t[1][2],
2102  t8 = t[0][1]*t[1][0],
2103  t00 = t[0][2]*t[1][0],
2104  t01 = t[0][1]*t[2][0],
2105  t04 = t[0][2]*t[2][0],
2106  inv_det_t = 1.0/(t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
2107  t00*t[2][1]+t01*t[1][2]-t04*t[1][1]);
2108  return_tensor[0][0] = t[1][1]*t[2][2]-t[1][2]*t[2][1];
2109  return_tensor[0][1] = t[0][2]*t[2][1]-t[0][1]*t[2][2];
2110  return_tensor[0][2] = t[0][1]*t[1][2]-t[0][2]*t[1][1];
2111  return_tensor[1][0] = t[1][2]*t[2][0]-t[1][0]*t[2][2];
2112  return_tensor[1][1] = t[0][0]*t[2][2]-t04;
2113  return_tensor[1][2] = t00-t6;
2114  return_tensor[2][0] = t[1][0]*t[2][1]-t[1][1]*t[2][0];
2115  return_tensor[2][1] = t01-t[0][0]*t[2][1];
2116  return_tensor[2][2] = t4-t8;
2117  return_tensor *= inv_det_t;
2118 
2119  return return_tensor;
2120 }
2121 
2122 #endif /* DOXYGEN */
2123 
2124 
2131 template <int dim, typename Number>
2132 inline
2135 {
2137  for (unsigned int i=0; i<dim; ++i)
2138  {
2139  tt[i][i] = t[i][i];
2140  for (unsigned int j=i+1; j<dim; ++j)
2141  {
2142  tt[i][j] = t[j][i];
2143  tt[j][i] = t[i][j];
2144  };
2145  }
2146  return tt;
2147 }
2148 
2149 
2163 template <int dim, typename Number>
2164 inline
2167 {
2168  return determinant(t)*invert(t);
2169 }
2170 
2171 
2186 template <int dim, typename Number>
2187 inline
2190 {
2191  return transpose(adjugate(t));
2192 }
2193 
2194 
2202 template <int dim, typename Number>
2203 inline
2204 Number
2206 {
2207  Number max = internal::NumberType<Number>::value(0.0);
2208  for (unsigned int j=0; j<dim; ++j)
2209  {
2210  Number sum = internal::NumberType<Number>::value(0.0);
2211  for (unsigned int i=0; i<dim; ++i)
2212  sum += std::fabs(t[i][j]);
2213 
2214  if (sum > max)
2215  max = sum;
2216  }
2217 
2218  return max;
2219 }
2220 
2221 
2229 template <int dim, typename Number>
2230 inline
2231 Number
2233 {
2234  Number max = internal::NumberType<Number>::value(0.0);
2235  for (unsigned int i=0; i<dim; ++i)
2236  {
2237  Number sum = internal::NumberType<Number>::value(0.0);
2238  for (unsigned int j=0; j<dim; ++j)
2239  sum += std::fabs(t[i][j]);
2240 
2241  if (sum > max)
2242  max = sum;
2243  }
2244 
2245  return max;
2246 }
2247 
2249 
2250 DEAL_II_NAMESPACE_CLOSE
2251 
2252 // include deprecated non-member functions operating on Tensor
2253 #include <deal.II/base/tensor_deprecated.h>
2254 
2255 #endif
SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
DEAL_II_CUDA_HOST_DEV Tensor()
Definition: tensor.h:954
numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:118
Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2024
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1202
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
Definition: tensor.h:1303
Tensor< rank_-1, dim, Number >::array_type array_type[(dim!=0)?dim:1]
Definition: tensor.h:405
static DEAL_II_CUDA_HOST_DEV real_type abs_square(const number &x)
Definition: numbers.h:345
static const unsigned int n_independent_components
Definition: tensor.h:391
DEAL_II_CUDA_HOST_DEV Tensor & operator*=(const OtherNumber &factor)
size_type size() const
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2205
static std::size_t memory_consumption()
Definition: tensor.h:1347
Tensor & operator/=(const OtherNumber &factor)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2232
Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2166
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:1713
static const unsigned int rank
Definition: tensor.h:384
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1254
#define AssertThrow(cond, exc)
Definition: exceptions.h:410
static real_type abs(const number &x)
Definition: numbers.h:354
bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1182
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
Definition: point.h:104
Tensor< rank_-1, dim, Number > values[(dim!=0)?dim:1]
Definition: tensor.h:627
void unroll(Vector< OtherNumber > &result) const
Definition: tensor.h:1279
Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2189
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2134
static::ExceptionBase & ExcMessage(std::string arg1)
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:1870
Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:2045
Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:1926
bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1155
#define Assert(cond, exc)
Definition: exceptions.h:349
ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:1837
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: tensor.h:1292
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:621
DEAL_II_ALWAYS_INLINE internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
void serialize(Archive &ar, const unsigned int version)
Definition: tensor.h:1357
Number * end_raw()
Definition: tensor.h:1107
Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Definition: tensor.h:1316
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:1775
Tensor< 1, dim, Number > cross_product_3d(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:1952
Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1982
static const unsigned int dimension
Definition: tensor.h:379
Tensor operator-() const
Definition: tensor.h:1240
Definition: mpi.h:53
Tensor< rank_-1, dim, Number >::tensor_type value_type
Definition: tensor.h:398
Number * begin_raw()
Definition: tensor.h:1087
DEAL_II_CUDA_HOST_DEV value_type & operator[](const unsigned int i)
Definition: tensor.h:1044
static::ExceptionBase & ExcNotImplemented()
DEAL_II_CUDA_HOST_DEV numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1264
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:1896
Tensor & operator+=(const Tensor< rank_, dim, OtherNumber > &)
void clear()
Definition: tensor.h:1337
static::ExceptionBase & ExcInternalError()
Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:2010
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const