Reference documentation for deal.II version 8.4.1
tensor.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2016 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/table_indices.h>
22 #include <deal.II/base/tensor_accessors.h>
23 #include <deal.II/base/template_constraints.h>
24 #include <deal.II/base/utilities.h>
25 
26 #include <cmath>
27 #include <ostream>
28 #include <vector>
29 
30 DEAL_II_NAMESPACE_OPEN
31 
32 // Forward declarations:
33 
34 template <int dim, typename Number> class Point;
35 template <int rank_, int dim, typename Number = double> class Tensor;
36 
37 #ifndef DOXYGEN
38 // Overload invalid tensor types of negative rank that come up during
39 // overload resolution of operator* and related contraction variants.
40 template <int dim, typename Number>
41 class Tensor<-2, dim, Number>
42 {
43 };
44 
45 template <int dim, typename Number>
46 class Tensor<-1, dim, Number>
47 {
48 };
49 #endif /* DOXYGEN */
50 
51 
82 template <int dim, typename Number>
83 class Tensor<0,dim,Number>
84 {
85 public:
94  static const unsigned int dimension = dim;
95 
99  static const unsigned int rank = 0;
100 
104  static const unsigned int n_independent_components = 1;
105 
115 
120  typedef Number value_type;
121 
127  typedef Number array_type;
128 
132  Tensor ();
133 
137  Tensor (const Tensor<0,dim,Number> &initializer);
138 
144  template <typename OtherNumber>
145  Tensor (const Tensor<0,dim,OtherNumber> &initializer);
146 
150  template <typename OtherNumber>
151  Tensor (const OtherNumber initializer);
152 
160  operator Number &();
161 
168  operator const Number &() const;
169 
174 
180  template <typename OtherNumber>
182 
186  template<typename OtherNumber>
187  bool operator == (const Tensor<0,dim,OtherNumber> &rhs) const;
188 
192  template<typename OtherNumber>
193  bool operator != (const Tensor<0,dim,OtherNumber> &rhs) const;
194 
198  template<typename OtherNumber>
200 
204  template<typename OtherNumber>
206 
210  template<typename OtherNumber>
211  Tensor<0,dim,Number> &operator *= (const OtherNumber factor);
212 
216  template<typename OtherNumber>
217  Tensor<0,dim,Number> &operator /= (const OtherNumber factor);
218 
223 
236  void clear ();
237 
243  real_type norm () const;
244 
249  real_type norm_square () const;
250 
255  template <class Archive>
256  void serialize(Archive &ar, const unsigned int version);
257 
262  typedef Number tensor_type;
263 
264 private:
268  Number value;
269 
273  template <typename OtherNumber>
275  unsigned int &start_index) const;
276 
280  template <int, int, typename> friend class Tensor;
281 };
282 
283 
284 
326 template <int rank_, int dim, typename Number>
327 class Tensor
328 {
329 public:
338  static const unsigned int dimension = dim;
339 
343  static const unsigned int rank = rank_;
344 
349  static const unsigned int
351 
357  typedef typename Tensor<rank_-1,dim,Number>::tensor_type value_type;
358 
363  typedef typename Tensor<rank_-1,dim,Number>::array_type
364  array_type[(dim != 0) ? dim : 1];
365  // ... avoid a compiler warning in case of dim == 0 and ensure that the
366  // array always has positive size.
367 
371  Tensor ();
372 
376  Tensor (const Tensor<rank_,dim,Number> &initializer);
377 
381  Tensor (const array_type &initializer);
382 
388  template <typename OtherNumber>
389  Tensor (const Tensor<rank_,dim,OtherNumber> &initializer);
390 
394  template <typename OtherNumber>
395  Tensor (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer);
396 
400  template <typename OtherNumber>
401  operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const;
402 
406  value_type &operator [] (const unsigned int i);
407 
411  const value_type &operator[](const unsigned int i) const;
412 
416  const Number &operator [] (const TableIndices<rank_> &indices) const;
417 
421  Number &operator [] (const TableIndices<rank_> &indices);
422 
427 
433  template <typename OtherNumber>
435 
442  Tensor<rank_,dim,Number> &operator = (const Number d);
443 
447  template <typename OtherNumber>
448  bool operator == (const Tensor<rank_,dim,OtherNumber> &) const;
449 
453  template <typename OtherNumber>
454  bool operator != (const Tensor<rank_,dim,OtherNumber> &) const;
455 
459  template <typename OtherNumber>
461 
465  template <typename OtherNumber>
467 
472  template <typename OtherNumber>
473  Tensor<rank_,dim,Number> &operator *= (const OtherNumber factor);
474 
478  template <typename OtherNumber>
479  Tensor<rank_,dim,Number> &operator /= (const OtherNumber factor);
480 
485 
498  void clear ();
499 
507 
513 
521  template <typename OtherNumber>
522  void unroll (Vector<OtherNumber> &result) const;
523 
528  static
529  unsigned int
531 
536  static
538 
543  static std::size_t memory_consumption ();
544 
549  template <class Archive>
550  void serialize(Archive &ar, const unsigned int version);
551 
557 
558 private:
562  Tensor<rank_-1, dim, Number> values[(dim != 0) ? dim : 1];
563  // ... avoid a compiler warning in case of dim == 0 and ensure that the
564  // array always has positive size.
565 
569  template <typename OtherNumber>
571  unsigned int &start_index) const;
572 
576  template <int, int, typename> friend class Tensor;
577 
582  friend class Point<dim,Number>;
583 };
584 
585 
586 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
587 
588 
589 template <int dim,typename Number>
590 inline
592  : value()
593 {
594 }
595 
596 
597 template <int dim, typename Number>
598 inline
600 {
601  value = p.value;
602 }
603 
604 
605 template <int dim, typename Number>
606 template <typename OtherNumber>
607 inline
608 Tensor<0,dim,Number>::Tensor (const OtherNumber initializer)
609 {
610  value = initializer;
611 }
612 
613 
614 template <int dim, typename Number>
615 template <typename OtherNumber>
616 inline
618 {
619  value = p.value;
620 }
621 
622 
623 template <int dim, typename Number>
624 inline
626 {
627  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
628  return value;
629 }
630 
631 
632 template <int dim, typename Number>
633 inline
634 Tensor<0,dim,Number>::operator const Number &() const
635 {
636  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
637  return value;
638 }
639 
640 
641 template <int dim, typename Number>
642 inline
644 {
645  value = p.value;
646  return *this;
647 }
648 
649 
650 template <int dim, typename Number>
651 template <typename OtherNumber>
652 inline
654 {
655  value = p.value;
656  return *this;
657 }
658 
659 
660 template <int dim, typename Number>
661 template <typename OtherNumber>
662 inline
664 {
665  return (value == p.value);
666 }
667 
668 
669 template <int dim, typename Number>
670 template <typename OtherNumber>
671 inline
673 {
674  return !((*this) == p);
675 }
676 
677 
678 template <int dim, typename Number>
679 template <typename OtherNumber>
680 inline
682 {
683  value += p.value;
684  return *this;
685 }
686 
687 
688 template <int dim, typename Number>
689 template <typename OtherNumber>
690 inline
692 {
693  value -= p.value;
694  return *this;
695 }
696 
697 
698 template <int dim, typename Number>
699 template <typename OtherNumber>
700 inline
702 {
703  value *= s;
704  return *this;
705 }
706 
707 
708 template <int dim, typename Number>
709 template <typename OtherNumber>
710 inline
712 {
713  value /= s;
714  return *this;
715 }
716 
717 
718 template <int dim, typename Number>
719 inline
721 {
722  return -value;
723 }
724 
725 
726 template <int dim, typename Number>
727 inline
730 {
731  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
732  return numbers::NumberTraits<Number>::abs (value);
733 }
734 
735 
736 template <int dim, typename Number>
737 inline
740 {
741  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
743 }
744 
745 
746 template <int dim, typename Number>
747 template <typename OtherNumber>
748 inline
749 void
751  unsigned int &index) const
752 {
753  Assert(dim != 0, ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
754  result[index] = value;
755  ++index;
756 }
757 
758 
759 template <int dim, typename Number>
760 inline
762 {
763  value = value_type();
764 }
765 
766 
767 template <int dim, typename Number>
768 template <class Archive>
769 inline
770 void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int)
771 {
772  ar &value;
773 }
774 
775 
776 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
777 
778 
779 template <int rank_, int dim, typename Number>
780 inline
782 {
783  // All members of the c-style array values are already default initialized
784  // and thus all values are already set to zero recursively.
785 }
786 
787 
788 template <int rank_, int dim, typename Number>
789 inline
791 {
792  if (dim > 0)
793  std::copy (&initializer[0], &initializer[0]+dim, &values[0]);
794 }
795 
796 
797 template <int rank_, int dim, typename Number>
798 inline
800 {
801  for (unsigned int i=0; i<dim; ++i)
802  values[i] = initializer[i];
803 }
804 
805 
806 template <int rank_, int dim, typename Number>
807 template <typename OtherNumber>
808 inline
810 {
811  for (unsigned int i=0; i!=dim; ++i)
812  values[i] = initializer[i];
813 }
814 
815 
816 template <int rank_, int dim, typename Number>
817 template <typename OtherNumber>
818 inline
820 (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer)
821 {
822  for (unsigned int i=0; i<dim; ++i)
823  values[i] = initializer[i];
824 }
825 
826 
827 template <int rank_, int dim, typename Number>
828 template <typename OtherNumber>
829 inline
831 operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const
832 {
833  return Tensor<1,dim,Tensor<rank_-1,dim,Number> > (values);
834 }
835 
836 
837 
838 namespace internal
839 {
840  namespace TensorSubscriptor
841  {
842  template <typename ArrayElementType, int dim>
843  inline DEAL_II_ALWAYS_INLINE
844  ArrayElementType &
845  subscript (ArrayElementType *values,
846  const unsigned int i,
848  {
849  Assert (i<dim, ExcIndexRange(i, 0, dim));
850  return values[i];
851  }
852 
853 
854  template <typename ArrayElementType>
855  ArrayElementType &
856  subscript (ArrayElementType *,
857  const unsigned int,
859  {
860  Assert(false, ExcMessage("Cannot access elements of an object of type Tensor<rank,0,Number>."));
861  static ArrayElementType t;
862  return t;
863  }
864  }
865 }
866 
867 
868 template <int rank_, int dim, typename Number>
869 inline DEAL_II_ALWAYS_INLINE
872 {
873  return ::internal::TensorSubscriptor::subscript(values, i, ::internal::int2type<dim>());
874 }
875 
876 
877 template <int rank_, int dim, typename Number>
878 inline DEAL_II_ALWAYS_INLINE
880 Tensor<rank_,dim,Number>::operator[] (const unsigned int i) const
881 {
882  return ::internal::TensorSubscriptor::subscript(values, i, ::internal::int2type<dim>());
883 }
884 
885 
886 template <int rank_, int dim, typename Number>
887 inline
888 const Number &
890 {
891  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
892 
893  return TensorAccessors::extract<rank_>(*this, indices);
894 }
895 
896 
897 template <int rank_, int dim, typename Number>
898 inline
899 Number &
901 {
902  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
903 
904  return TensorAccessors::extract<rank_>(*this, indices);
905 }
906 
907 
908 template <int rank_, int dim, typename Number>
909 inline
912 {
913  if (dim > 0)
914  std::copy (&t.values[0], &t.values[0]+dim, &values[0]);
915  return *this;
916 }
917 
918 
919 template <int rank_, int dim, typename Number>
920 template <typename OtherNumber>
921 inline
924 {
925  if (dim > 0)
926  std::copy (&t.values[0], &t.values[0]+dim, &values[0]);
927  return *this;
928 }
929 
930 
931 template <int rank_, int dim, typename Number>
932 inline
935 {
936  Assert (d == Number(), ExcMessage ("Only assignment with zero is allowed"));
937  (void) d;
938 
939  for (unsigned int i=0; i<dim; ++i)
940  values[i] = Number();
941  return *this;
942 }
943 
944 
945 template <int rank_, int dim, typename Number>
946 template <typename OtherNumber>
947 inline
948 bool
950 {
951  for (unsigned int i=0; i<dim; ++i)
952  if (values[i] != p.values[i])
953  return false;
954  return true;
955 }
956 
957 
958 // At some places in the library, we have Point<0> for formal reasons
959 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
960 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
961 // in the above function that the loop end check always fails, we
962 // implement this function here
963 template <>
964 template <>
965 inline
967 {
968  return true;
969 }
970 
971 
972 template <int rank_, int dim, typename Number>
973 template <typename OtherNumber>
974 inline
975 bool
977 {
978  return !((*this) == p);
979 }
980 
981 
982 template <int rank_, int dim, typename Number>
983 template <typename OtherNumber>
984 inline
987 {
988  for (unsigned int i=0; i<dim; ++i)
989  values[i] += p.values[i];
990  return *this;
991 }
992 
993 
994 template <int rank_, int dim, typename Number>
995 template <typename OtherNumber>
996 inline
999 {
1000  for (unsigned int i=0; i<dim; ++i)
1001  values[i] -= p.values[i];
1002  return *this;
1003 }
1004 
1005 
1006 template <int rank_, int dim, typename Number>
1007 template <typename OtherNumber>
1008 inline
1011 {
1012  for (unsigned int i=0; i<dim; ++i)
1013  values[i] *= s;
1014  return *this;
1015 }
1016 
1017 
1018 template <int rank_, int dim, typename Number>
1019 template <typename OtherNumber>
1020 inline
1023 {
1024  for (unsigned int i=0; i<dim; ++i)
1025  values[i] /= s;
1026  return *this;
1027 }
1028 
1029 
1030 template <int rank_, int dim, typename Number>
1031 inline
1034 {
1036 
1037  for (unsigned int i=0; i<dim; ++i)
1038  tmp.values[i] = -values[i];
1039 
1040  return tmp;
1041 }
1042 
1043 
1044 template <int rank_, int dim, typename Number>
1045 inline
1048 {
1049  return std::sqrt (norm_square());
1050 }
1051 
1052 
1053 template <int rank_, int dim, typename Number>
1054 inline
1057 {
1059  for (unsigned int i=0; i<dim; ++i)
1060  s += values[i].norm_square();
1061 
1062  return s;
1063 }
1064 
1065 
1066 template <int rank_, int dim, typename Number>
1067 template <typename OtherNumber>
1068 inline
1069 void
1071 {
1072  AssertDimension (result.size(),(Utilities::fixed_power<rank_, unsigned int>(dim)));
1073 
1074  unsigned int index = 0;
1075  unroll_recursion (result, index);
1076 }
1077 
1078 
1079 template <int rank_, int dim, typename Number>
1080 template <typename OtherNumber>
1081 inline
1082 void
1084  unsigned int &index) const
1085 {
1086  for (unsigned int i=0; i<dim; ++i)
1087  values[i].unroll_recursion(result, index);
1088 }
1089 
1090 
1091 template <int rank_, int dim, typename Number>
1092 inline
1093 unsigned int
1095 {
1096  unsigned int index = 0;
1097  for (int r = 0; r < rank_; ++r)
1098  index = index * dim + indices[r];
1099 
1100  return index;
1101 }
1102 
1103 
1104 template <int rank_, int dim, typename Number>
1105 inline
1108 {
1109  Assert (i < n_independent_components,
1110  ExcIndexRange (i, 0, n_independent_components));
1111 
1112  TableIndices<rank_> indices;
1113 
1114  unsigned int remainder = i;
1115  for (int r=rank_-1; r>=0; --r)
1116  {
1117  indices[r] = (remainder % dim);
1118  remainder /= dim;
1119  }
1120  Assert (remainder == 0, ExcInternalError());
1121 
1122  return indices;
1123 }
1124 
1125 
1126 template <int rank_, int dim, typename Number>
1127 inline
1129 {
1130  for (unsigned int i=0; i<dim; ++i)
1131  values[i] = value_type();
1132 }
1133 
1134 
1135 template <int rank_, int dim, typename Number>
1136 inline
1137 std::size_t
1139 {
1140  return sizeof(Tensor<rank_,dim,Number>);
1141 }
1142 
1143 
1144 template <int rank_, int dim, typename Number>
1145 template <class Archive>
1146 inline
1147 void
1148 Tensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
1149 {
1150  ar &values;
1151 }
1152 
1153 
1154 /* ----------------- Non-member functions operating on tensors. ------------ */
1155 
1160 
1168 template <int rank_, int dim, typename Number>
1169 inline
1170 std::ostream &operator << (std::ostream &out, const Tensor<rank_,dim,Number> &p)
1171 {
1172  for (unsigned int i = 0; i < dim; ++i)
1173  {
1174  out << p[i];
1175  if (i != dim - 1)
1176  out << ' ';
1177  }
1178 
1179  return out;
1180 }
1181 
1182 
1189 template <int dim, typename Number>
1190 inline
1191 std::ostream &operator << (std::ostream &out, const Tensor<0,dim,Number> &p)
1192 {
1193  out << static_cast<const Number &>(p);
1194  return out;
1195 }
1196 
1197 
1199 
1203 
1204 
1205 #ifndef DEAL_II_WITH_CXX11
1206 template <typename T, typename U, int rank, int dim>
1207 struct ProductType<T,Tensor<rank,dim,U> >
1208 {
1210 };
1211 
1212 template <typename T, typename U, int rank, int dim>
1213 struct ProductType<Tensor<rank,dim,T>,U>
1214 {
1216 };
1217 #endif
1218 
1219 
1220 
1229 template <int dim, typename Number, typename Other>
1230 inline
1232 operator * (const Other object,
1233  const Tensor<0,dim,Number> &t)
1234 {
1235  return object * static_cast<const Number &>(t);
1236 }
1237 
1238 
1247 template <int dim, typename Number, typename Other>
1248 inline
1251  const Other object)
1252 {
1253  return static_cast<const Number &>(t) * object;
1254 }
1255 
1256 
1266 template <int dim, typename Number, typename OtherNumber>
1267 inline
1270  const Tensor<0, dim, OtherNumber> &src2)
1271 {
1272  return static_cast<const Number &>(src1) *
1273  static_cast<const OtherNumber &>(src2);
1274 }
1275 
1276 
1282 template <int dim, typename Number, typename OtherNumber>
1283 inline
1286  const OtherNumber factor)
1287 {
1288  return static_cast<Number>(t) / factor;
1289 }
1290 
1291 
1297 template <int dim, typename Number, typename OtherNumber>
1298 inline
1301 {
1302  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1303 }
1304 
1305 
1311 template <int dim, typename Number, typename OtherNumber>
1312 inline
1315 {
1316  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1317 }
1318 
1319 
1330 template <int rank, int dim,
1331  typename Number,
1332  typename OtherNumber>
1333 inline
1336  const OtherNumber factor)
1337 {
1338  // recurse over the base objects
1340  for (unsigned int d=0; d<dim; ++d)
1341  tt[d] = t[d] * factor;
1342  return tt;
1343 }
1344 
1345 
1356 template <int rank, int dim,
1357  typename Number,
1358  typename OtherNumber>
1359 inline
1361 operator * (const Number factor,
1363 {
1364  // simply forward to the operator above
1365  return t * factor;
1366 }
1367 
1368 
1376 template <int rank, int dim,
1377  typename Number,
1378  typename OtherNumber>
1379 inline
1382  const OtherNumber factor)
1383 {
1384  // recurse over the base objects
1386  for (unsigned int d=0; d<dim; ++d)
1387  tt[d] = t[d] / factor;
1388  return tt;
1389 }
1390 
1391 
1399 template <int rank, int dim, typename Number, typename OtherNumber>
1400 inline
1403 {
1405 
1406  for (unsigned int i=0; i<dim; ++i)
1407  tmp[i] += q[i];
1408 
1409  return tmp;
1410 }
1411 
1412 
1420 template <int rank, int dim, typename Number, typename OtherNumber>
1421 inline
1424 {
1426 
1427  for (unsigned int i=0; i<dim; ++i)
1428  tmp[i] -= q[i];
1429 
1430  return tmp;
1431 }
1432 
1433 
1435 
1439 
1440 
1464 template <int rank_1, int rank_2, int dim,
1465  typename Number, typename OtherNumber>
1466 inline DEAL_II_ALWAYS_INLINE
1467 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1470 {
1471  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
1472 
1473  TensorAccessors::internal::ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1474  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
1475  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
1476 
1477  return result;
1478 }
1479 
1480 
1510 template <int index_1, int index_2,
1511  int rank_1, int rank_2, int dim,
1512  typename Number, typename OtherNumber>
1513 inline
1514 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1517 {
1518  Assert(0 <= index_1 && index_1 < rank_1,
1519  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1520  Assert(0 <= index_2 && index_2 < rank_2,
1521  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1522 
1523  using namespace TensorAccessors;
1524  using namespace TensorAccessors::internal;
1525 
1526  // Reorder index_1 to the end of src1:
1527  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1528  reord_01 = reordered_index_view<index_1, rank_1>(src1);
1529 
1530  // Reorder index_2 to the end of src2:
1531  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1532  reord_02 = reordered_index_view<index_2, rank_2>(src2);
1533 
1534  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1535  result;
1536  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
1537  return result;
1538 }
1539 
1540 
1572 template <int index_1, int index_2, int index_3, int index_4,
1573  int rank_1, int rank_2, int dim,
1574  typename Number, typename OtherNumber>
1575 inline
1576 typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1579 {
1580  Assert(0 <= index_1 && index_1 < rank_1,
1581  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1582  Assert(0 <= index_3 && index_3 < rank_1,
1583  ExcMessage("The specified index_3 must lie within the range [0,rank_1)"));
1584  Assert(index_1 != index_3,
1585  ExcMessage("index_1 and index_3 must not be the same"));
1586  Assert(0 <= index_2 && index_2 < rank_2,
1587  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1588  Assert(0 <= index_4 && index_4 < rank_2,
1589  ExcMessage("The specified index_4 must lie within the range [0,rank_2)"));
1590  Assert(index_2 != index_4,
1591  ExcMessage("index_2 and index_4 must not be the same"));
1592 
1593  using namespace TensorAccessors;
1594  using namespace TensorAccessors::internal;
1595 
1596  // Reorder index_1 to the end of src1:
1597  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1598  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
1599 
1600  // Reorder index_2 to the end of src2:
1601  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1602  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
1603 
1604  // Now, reorder index_3 to the end of src1. We have to make sure to
1605  // preserve the orginial ordering: index_1 has been removed. If
1606  // index_3 > index_1, we have to use (index_3 - 1) instead:
1607  ReorderedIndexView<(index_3 < index_1 ? index_3 : index_3 - 1), rank_1, ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> > >
1608  reord_3 = TensorAccessors::reordered_index_view<index_3 < index_1 ? index_3 : index_3 - 1, rank_1>(reord_1);
1609 
1610  // Now, reorder index_4 to the end of src2. We have to make sure to
1611  // preserve the orginial ordering: index_2 has been removed. If
1612  // index_4 > index_2, we have to use (index_4 - 1) instead:
1613  ReorderedIndexView<(index_4 < index_2 ? index_4 : index_4 - 1), rank_2, ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> > >
1614  reord_4 = TensorAccessors::reordered_index_view<index_4 < index_2 ? index_4 : index_4 - 1, rank_2>(reord_2);
1615 
1616  typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1617  result;
1618  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
1619  return result;
1620 }
1621 
1622 
1636 template <int rank, int dim, typename Number, typename OtherNumber>
1637 inline
1640  const Tensor<rank, dim, OtherNumber> &right)
1641 {
1643  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
1644  return result;
1645 }
1646 
1647 
1663 template <int rank_1, int rank_2, int dim,
1664  typename T1, typename T2, typename T3>
1667  const Tensor<rank_1 + rank_2, dim, T2> &middle,
1668  const Tensor<rank_2, dim, T3> &right)
1669 {
1671  return_type;
1672  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(
1673  left, middle, right);
1674 }
1675 
1676 
1688 template <int rank_1, int rank_2, int dim,
1689  typename Number, typename OtherNumber>
1690 inline
1694 {
1696  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
1697  return result;
1698 }
1699 
1700 
1702 
1706 
1707 
1719 template <int dim, typename Number>
1720 inline
1723 {
1724  Assert (dim==2, ExcInternalError());
1725 
1726  Tensor<1, dim, Number> result;
1727 
1728  result[0] = src[1];
1729  result[1] = -src[0];
1730 
1731  return result;
1732 }
1733 
1734 
1745 template <int dim, typename Number>
1746 inline
1749  const Tensor<1,dim,Number> &src2)
1750 {
1751  Assert (dim==3, ExcInternalError());
1752 
1753  Tensor<1, dim, Number> result;
1754 
1755  result[0] = src1[1]*src2[2] - src1[2]*src2[1];
1756  result[1] = src1[2]*src2[0] - src1[0]*src2[2];
1757  result[2] = src1[0]*src2[1] - src1[1]*src2[0];
1758 
1759  return result;
1760 }
1761 
1762 
1764 
1768 
1769 
1776 template <int dim, typename Number>
1777 inline
1779 {
1780  // Compute the determinant using the Laplace expansion of the
1781  // determinant. We expand along the last row.
1782  Number det = Number();
1783 
1784  for (unsigned int k=0; k<dim; ++k)
1785  {
1786  Tensor<2,dim-1,Number> minor;
1787  for (unsigned int i=0; i<dim-1; ++i)
1788  for (unsigned int j=0; j<dim-1; ++j)
1789  minor[i][j] = t[i][j<k ? j : j+1];
1790 
1791  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
1792 
1793  det += t[dim-1][k] * cofactor;
1794  }
1795 
1796  return ((dim % 2 == 0) ? 1. : -1.) * det;
1797 }
1798 
1804 template <typename Number>
1805 inline
1807 {
1808  return t[0][0];
1809 }
1810 
1811 
1819 template <int dim, typename Number>
1820 Number trace (const Tensor<2,dim,Number> &d)
1821 {
1822  Number t=d[0][0];
1823  for (unsigned int i=1; i<dim; ++i)
1824  t += d[i][i];
1825  return t;
1826 }
1827 
1828 
1838 template <int dim, typename Number>
1839 inline
1842 {
1843  Number return_tensor [dim][dim];
1844  switch (dim)
1845  {
1846  case 1:
1847  return_tensor[0][0] = 1.0/t[0][0];
1848  break;
1849 
1850  case 2:
1851  // this is Maple output,
1852  // thus a bit unstructured
1853  {
1854  const Number det = t[0][0]*t[1][1]-t[1][0]*t[0][1];
1855  const Number t4 = 1.0/det;
1856  return_tensor[0][0] = t[1][1]*t4;
1857  return_tensor[0][1] = -t[0][1]*t4;
1858  return_tensor[1][0] = -t[1][0]*t4;
1859  return_tensor[1][1] = t[0][0]*t4;
1860  break;
1861  }
1862 
1863  case 3:
1864  {
1865  const Number t4 = t[0][0]*t[1][1],
1866  t6 = t[0][0]*t[1][2],
1867  t8 = t[0][1]*t[1][0],
1868  t00 = t[0][2]*t[1][0],
1869  t01 = t[0][1]*t[2][0],
1870  t04 = t[0][2]*t[2][0],
1871  det = (t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
1872  t00*t[2][1]+t01*t[1][2]-t04*t[1][1]),
1873  t07 = 1.0/det;
1874  return_tensor[0][0] = (t[1][1]*t[2][2]-t[1][2]*t[2][1])*t07;
1875  return_tensor[0][1] = (t[0][2]*t[2][1]-t[0][1]*t[2][2])*t07;
1876  return_tensor[0][2] = (t[0][1]*t[1][2]-t[0][2]*t[1][1])*t07;
1877  return_tensor[1][0] = (t[1][2]*t[2][0]-t[1][0]*t[2][2])*t07;
1878  return_tensor[1][1] = (t[0][0]*t[2][2]-t04)*t07;
1879  return_tensor[1][2] = (t00-t6)*t07;
1880  return_tensor[2][0] = (t[1][0]*t[2][1]-t[1][1]*t[2][0])*t07;
1881  return_tensor[2][1] = (t01-t[0][0]*t[2][1])*t07;
1882  return_tensor[2][2] = (t4-t8)*t07;
1883 
1884  break;
1885  }
1886 
1887  // if desired, take over the
1888  // inversion of a 4x4 tensor
1889  // from the FullMatrix
1890  default:
1891  AssertThrow (false, ExcNotImplemented());
1892  }
1893  return Tensor<2,dim,Number>(return_tensor);
1894 }
1895 
1896 
1903 template <int dim, typename Number>
1904 inline
1907 {
1909  for (unsigned int i=0; i<dim; ++i)
1910  {
1911  tt[i][i] = t[i][i];
1912  for (unsigned int j=i+1; j<dim; ++j)
1913  {
1914  tt[i][j] = t[j][i];
1915  tt[j][i] = t[i][j];
1916  };
1917  }
1918  return tt;
1919 }
1920 
1921 
1929 template <int dim, typename Number>
1930 inline
1931 double
1933 {
1934  double max = 0;
1935  for (unsigned int j=0; j<dim; ++j)
1936  {
1937  double sum = 0;
1938  for (unsigned int i=0; i<dim; ++i)
1939  sum += std::fabs(t[i][j]);
1940 
1941  if (sum > max)
1942  max = sum;
1943  }
1944 
1945  return max;
1946 }
1947 
1948 
1956 template <int dim, typename Number>
1957 inline
1958 double
1960 {
1961  double max = 0;
1962  for (unsigned int i=0; i<dim; ++i)
1963  {
1964  double sum = 0;
1965  for (unsigned int j=0; j<dim; ++j)
1966  sum += std::fabs(t[i][j]);
1967 
1968  if (sum > max)
1969  max = sum;
1970  }
1971 
1972  return max;
1973 }
1974 
1976 
1977 DEAL_II_NAMESPACE_CLOSE
1978 
1979 // include deprecated non-member functions operating on Tensor
1980 #include <deal.II/base/tensor_deprecated.h>
1981 
1982 #endif
1983 
numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:114
Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:1820
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1052
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
Definition: tensor.h:1094
Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Tensor< rank_-1, dim, Number >::array_type array_type[(dim!=0)?dim:1]
Definition: tensor.h:364
SymmetricTensor< rank, dim, Number > operator/(const SymmetricTensor< rank, dim, Number > &t, const Number factor)
static const unsigned int n_independent_components
Definition: tensor.h:350
::ExceptionBase & ExcMessage(std::string arg1)
static std::size_t memory_consumption()
Definition: tensor.h:1138
Tensor< rank_, dim, Number > & operator-=(const Tensor< rank_, dim, OtherNumber > &)
Definition: tensor.h:998
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:1515
static const unsigned int rank
Definition: tensor.h:343
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1047
#define AssertThrow(cond, exc)
Definition: exceptions.h:358
static real_type abs(const number &x)
Definition: numbers.h:322
bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:976
Definition: point.h:89
Tensor< rank_-1, dim, Number > values[(dim!=0)?dim:1]
Definition: tensor.h:562
void unroll(Vector< OtherNumber > &result) const
Definition: tensor.h:1070
Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1906
Tensor< rank_, dim, Number > & operator+=(const Tensor< rank_, dim, OtherNumber > &)
Definition: tensor.h:986
static real_type abs_square(const number &x)
Definition: numbers.h:313
Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:1722
bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:949
#define Assert(cond, exc)
Definition: exceptions.h:294
ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:1639
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: tensor.h:1083
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:556
DEAL_II_ALWAYS_INLINE internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
void serialize(Archive &ar, const unsigned int version)
Definition: tensor.h:1148
double linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1959
std::size_t size() const
Tensor()
Definition: tensor.h:781
double l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1932
Tensor< rank_, dim, Number > & operator/=(const OtherNumber factor)
Definition: tensor.h:1022
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Definition: tensor.h:1107
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:1577
value_type & operator[](const unsigned int i)
Definition: tensor.h:871
Tensor< 1, dim, Number > cross_product_3d(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:1748
ProductType< T1, typename ProductType< T2, T3 >::type >::type contract3(const Tensor< rank_1, dim, T1 > &left, const Tensor< rank_1+rank_2, dim, T2 > &middle, const Tensor< rank_2, dim, T3 > &right)
Definition: tensor.h:1666
Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1778
static const unsigned int dimension
Definition: tensor.h:338
Tensor< rank_, dim, Number > operator-() const
Definition: tensor.h:1033
Definition: mpi.h:48
Tensor< rank_-1, dim, Number >::tensor_type value_type
Definition: tensor.h:357
Tensor & operator=(const Tensor< rank_, dim, Number > &rhs)
Definition: tensor.h:911
Tensor< rank_, dim, Number > & operator*=(const OtherNumber factor)
Definition: tensor.h:1010
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:1692
numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1056
void clear()
Definition: tensor.h:1128
Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1841
Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:1806
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const