Reference documentation for deal.II version Git d1151b6 2017-02-19 18:14:48 -0700
tensor.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2017 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 template <typename Number> class Vector;
37 
38 #ifndef DOXYGEN
39 // Overload invalid tensor types of negative rank that come up during
40 // overload resolution of operator* and related contraction variants.
41 template <int dim, typename Number>
42 class Tensor<-2, dim, Number>
43 {
44 };
45 
46 template <int dim, typename Number>
47 class Tensor<-1, dim, Number>
48 {
49 };
50 #endif /* DOXYGEN */
51 
52 
83 template <int dim, typename Number>
84 class Tensor<0,dim,Number>
85 {
86 public:
95  static const unsigned int dimension = dim;
96 
100  static const unsigned int rank = 0;
101 
105  static const unsigned int n_independent_components = 1;
106 
116 
121  typedef Number value_type;
122 
128  typedef Number array_type;
129 
133  Tensor ();
134 
140  template <typename OtherNumber>
141  Tensor (const Tensor<0,dim,OtherNumber> &initializer);
142 
146  template <typename OtherNumber>
147  Tensor (const OtherNumber initializer);
148 
156  operator Number &();
157 
164  operator const Number &() const;
165 
171  template <typename OtherNumber>
173 
177  template<typename OtherNumber>
178  bool operator == (const Tensor<0,dim,OtherNumber> &rhs) const;
179 
183  template<typename OtherNumber>
184  bool operator != (const Tensor<0,dim,OtherNumber> &rhs) const;
185 
189  template<typename OtherNumber>
191 
195  template<typename OtherNumber>
197 
201  template<typename OtherNumber>
202  Tensor<0,dim,Number> &operator *= (const OtherNumber factor);
203 
207  template<typename OtherNumber>
208  Tensor<0,dim,Number> &operator /= (const OtherNumber factor);
209 
214 
227  void clear ();
228 
234  real_type norm () const;
235 
240  real_type norm_square () const;
241 
246  template <class Archive>
247  void serialize(Archive &ar, const unsigned int version);
248 
253  typedef Number tensor_type;
254 
255 private:
259  Number value;
260 
264  template <typename OtherNumber>
266  unsigned int &start_index) const;
267 
271  template <int, int, typename> friend class Tensor;
272 };
273 
274 
275 
317 template <int rank_, int dim, typename Number>
318 class Tensor
319 {
320 public:
329  static const unsigned int dimension = dim;
330 
334  static const unsigned int rank = rank_;
335 
340  static const unsigned int
342 
348  typedef typename Tensor<rank_-1,dim,Number>::tensor_type value_type;
349 
354  typedef typename Tensor<rank_-1,dim,Number>::array_type
355  array_type[(dim != 0) ? dim : 1];
356  // ... avoid a compiler warning in case of dim == 0 and ensure that the
357  // array always has positive size.
358 
362  Tensor ();
363 
367  Tensor (const array_type &initializer);
368 
374  template <typename OtherNumber>
375  Tensor (const Tensor<rank_,dim,OtherNumber> &initializer);
376 
380  template <typename OtherNumber>
381  Tensor (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer);
382 
386  template <typename OtherNumber>
387  operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const;
388 
392  value_type &operator [] (const unsigned int i);
393 
397  const value_type &operator[](const unsigned int i) const;
398 
402  const Number &operator [] (const TableIndices<rank_> &indices) const;
403 
407  Number &operator [] (const TableIndices<rank_> &indices);
408 
414  template <typename OtherNumber>
416 
423  Tensor<rank_,dim,Number> &operator = (const Number d);
424 
428  template <typename OtherNumber>
429  bool operator == (const Tensor<rank_,dim,OtherNumber> &) const;
430 
434  template <typename OtherNumber>
435  bool operator != (const Tensor<rank_,dim,OtherNumber> &) const;
436 
440  template <typename OtherNumber>
442 
446  template <typename OtherNumber>
448 
453  template <typename OtherNumber>
454  Tensor<rank_,dim,Number> &operator *= (const OtherNumber factor);
455 
459  template <typename OtherNumber>
460  Tensor<rank_,dim,Number> &operator /= (const OtherNumber factor);
461 
466 
479  void clear ();
480 
488 
494 
502  template <typename OtherNumber>
503  void unroll (Vector<OtherNumber> &result) const;
504 
509  static
510  unsigned int
512 
517  static
519 
524  static std::size_t memory_consumption ();
525 
530  template <class Archive>
531  void serialize(Archive &ar, const unsigned int version);
532 
538 
539 private:
543  Tensor<rank_-1, dim, Number> values[(dim != 0) ? dim : 1];
544  // ... avoid a compiler warning in case of dim == 0 and ensure that the
545  // array always has positive size.
546 
550  template <typename OtherNumber>
552  unsigned int &start_index) const;
553 
557  template <int, int, typename> friend class Tensor;
558 
563  friend class Point<dim,Number>;
564 };
565 
566 
567 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
568 
569 
570 template <int dim,typename Number>
571 inline
573  : value()
574 {
575 }
576 
577 
578 template <int dim, typename Number>
579 template <typename OtherNumber>
580 inline
581 Tensor<0,dim,Number>::Tensor (const OtherNumber initializer)
582 {
583  value = initializer;
584 }
585 
586 
587 template <int dim, typename Number>
588 template <typename OtherNumber>
589 inline
591 {
592  value = p.value;
593 }
594 
595 
596 template <int dim, typename Number>
597 inline
599 {
600  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
601  return value;
602 }
603 
604 
605 template <int dim, typename Number>
606 inline
607 Tensor<0,dim,Number>::operator const Number &() const
608 {
609  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
610  return value;
611 }
612 
613 
614 template <int dim, typename Number>
615 template <typename OtherNumber>
616 inline
618 {
619  value = p.value;
620  return *this;
621 }
622 
623 
624 template <int dim, typename Number>
625 template <typename OtherNumber>
626 inline
628 {
629  return (value == p.value);
630 }
631 
632 
633 template <int dim, typename Number>
634 template <typename OtherNumber>
635 inline
637 {
638  return !((*this) == p);
639 }
640 
641 
642 template <int dim, typename Number>
643 template <typename OtherNumber>
644 inline
646 {
647  value += p.value;
648  return *this;
649 }
650 
651 
652 template <int dim, typename Number>
653 template <typename OtherNumber>
654 inline
656 {
657  value -= p.value;
658  return *this;
659 }
660 
661 
662 template <int dim, typename Number>
663 template <typename OtherNumber>
664 inline
666 {
667  value *= s;
668  return *this;
669 }
670 
671 
672 template <int dim, typename Number>
673 template <typename OtherNumber>
674 inline
676 {
677  value /= s;
678  return *this;
679 }
680 
681 
682 template <int dim, typename Number>
683 inline
685 {
686  return -value;
687 }
688 
689 
690 template <int dim, typename Number>
691 inline
694 {
695  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
696  return numbers::NumberTraits<Number>::abs (value);
697 }
698 
699 
700 template <int dim, typename Number>
701 inline
704 {
705  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
707 }
708 
709 
710 template <int dim, typename Number>
711 template <typename OtherNumber>
712 inline
713 void
715  unsigned int &index) const
716 {
717  Assert(dim != 0, ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
718  result[index] = value;
719  ++index;
720 }
721 
722 
723 template <int dim, typename Number>
724 inline
726 {
727  value = value_type();
728 }
729 
730 
731 template <int dim, typename Number>
732 template <class Archive>
733 inline
734 void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int)
735 {
736  ar &value;
737 }
738 
739 
740 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
741 
742 
743 template <int rank_, int dim, typename Number>
744 inline
746 {
747  // All members of the c-style array values are already default initialized
748  // and thus all values are already set to zero recursively.
749 }
750 
751 
752 template <int rank_, int dim, typename Number>
753 inline
755 {
756  for (unsigned int i=0; i<dim; ++i)
757  values[i] = initializer[i];
758 }
759 
760 
761 template <int rank_, int dim, typename Number>
762 template <typename OtherNumber>
763 inline
765 {
766  for (unsigned int i=0; i!=dim; ++i)
767  values[i] = initializer[i];
768 }
769 
770 
771 template <int rank_, int dim, typename Number>
772 template <typename OtherNumber>
773 inline
775 (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer)
776 {
777  for (unsigned int i=0; i<dim; ++i)
778  values[i] = initializer[i];
779 }
780 
781 
782 template <int rank_, int dim, typename Number>
783 template <typename OtherNumber>
784 inline
786 operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const
787 {
788  return Tensor<1,dim,Tensor<rank_-1,dim,Number> > (values);
789 }
790 
791 
792 
793 namespace internal
794 {
795  namespace TensorSubscriptor
796  {
797  template <typename ArrayElementType, int dim>
798  inline DEAL_II_ALWAYS_INLINE
799  ArrayElementType &
800  subscript (ArrayElementType *values,
801  const unsigned int i,
803  {
804  Assert (i<dim, ExcIndexRange(i, 0, dim));
805  return values[i];
806  }
807 
808 
809  template <typename ArrayElementType>
810  ArrayElementType &
811  subscript (ArrayElementType *,
812  const unsigned int,
814  {
815  Assert(false, ExcMessage("Cannot access elements of an object of type Tensor<rank,0,Number>."));
816  static ArrayElementType t;
817  return t;
818  }
819  }
820 }
821 
822 
823 template <int rank_, int dim, typename Number>
824 inline DEAL_II_ALWAYS_INLINE
827 {
828  return ::internal::TensorSubscriptor::subscript(values, i, ::internal::int2type<dim>());
829 }
830 
831 
832 template <int rank_, int dim, typename Number>
833 inline DEAL_II_ALWAYS_INLINE
835 Tensor<rank_,dim,Number>::operator[] (const unsigned int i) const
836 {
837  return ::internal::TensorSubscriptor::subscript(values, i, ::internal::int2type<dim>());
838 }
839 
840 
841 template <int rank_, int dim, typename Number>
842 inline
843 const Number &
845 {
846  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
847 
848  return TensorAccessors::extract<rank_>(*this, indices);
849 }
850 
851 
852 template <int rank_, int dim, typename Number>
853 inline
854 Number &
856 {
857  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
858 
859  return TensorAccessors::extract<rank_>(*this, indices);
860 }
861 
862 
863 template <int rank_, int dim, typename Number>
864 template <typename OtherNumber>
865 inline
868 {
869  if (dim > 0)
870  std::copy (&t.values[0], &t.values[0]+dim, &values[0]);
871  return *this;
872 }
873 
874 
875 template <int rank_, int dim, typename Number>
876 inline
879 {
880  Assert (d == Number(), ExcMessage ("Only assignment with zero is allowed"));
881  (void) d;
882 
883  for (unsigned int i=0; i<dim; ++i)
884  values[i] = Number();
885  return *this;
886 }
887 
888 
889 template <int rank_, int dim, typename Number>
890 template <typename OtherNumber>
891 inline
892 bool
894 {
895  for (unsigned int i=0; i<dim; ++i)
896  if (values[i] != p.values[i])
897  return false;
898  return true;
899 }
900 
901 
902 // At some places in the library, we have Point<0> for formal reasons
903 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
904 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
905 // in the above function that the loop end check always fails, we
906 // implement this function here
907 template <>
908 template <>
909 inline
911 {
912  return true;
913 }
914 
915 
916 template <int rank_, int dim, typename Number>
917 template <typename OtherNumber>
918 inline
919 bool
921 {
922  return !((*this) == p);
923 }
924 
925 
926 template <int rank_, int dim, typename Number>
927 template <typename OtherNumber>
928 inline
931 {
932  for (unsigned int i=0; i<dim; ++i)
933  values[i] += p.values[i];
934  return *this;
935 }
936 
937 
938 template <int rank_, int dim, typename Number>
939 template <typename OtherNumber>
940 inline
943 {
944  for (unsigned int i=0; i<dim; ++i)
945  values[i] -= p.values[i];
946  return *this;
947 }
948 
949 
950 template <int rank_, int dim, typename Number>
951 template <typename OtherNumber>
952 inline
955 {
956  for (unsigned int i=0; i<dim; ++i)
957  values[i] *= s;
958  return *this;
959 }
960 
961 
962 template <int rank_, int dim, typename Number>
963 template <typename OtherNumber>
964 inline
967 {
968  for (unsigned int i=0; i<dim; ++i)
969  values[i] /= s;
970  return *this;
971 }
972 
973 
974 template <int rank_, int dim, typename Number>
975 inline
978 {
980 
981  for (unsigned int i=0; i<dim; ++i)
982  tmp.values[i] = -values[i];
983 
984  return tmp;
985 }
986 
987 
988 template <int rank_, int dim, typename Number>
989 inline
992 {
993  return std::sqrt (norm_square());
994 }
995 
996 
997 template <int rank_, int dim, typename Number>
998 inline
1001 {
1003  for (unsigned int i=0; i<dim; ++i)
1004  s += values[i].norm_square();
1005 
1006  return s;
1007 }
1008 
1009 
1010 template <int rank_, int dim, typename Number>
1011 template <typename OtherNumber>
1012 inline
1013 void
1015 {
1016  AssertDimension (result.size(),(Utilities::fixed_power<rank_, unsigned int>(dim)));
1017 
1018  unsigned int index = 0;
1019  unroll_recursion (result, index);
1020 }
1021 
1022 
1023 template <int rank_, int dim, typename Number>
1024 template <typename OtherNumber>
1025 inline
1026 void
1028  unsigned int &index) const
1029 {
1030  for (unsigned int i=0; i<dim; ++i)
1031  values[i].unroll_recursion(result, index);
1032 }
1033 
1034 
1035 template <int rank_, int dim, typename Number>
1036 inline
1037 unsigned int
1039 {
1040  unsigned int index = 0;
1041  for (int r = 0; r < rank_; ++r)
1042  index = index * dim + indices[r];
1043 
1044  return index;
1045 }
1046 
1047 
1048 template <int rank_, int dim, typename Number>
1049 inline
1052 {
1053  Assert (i < n_independent_components,
1054  ExcIndexRange (i, 0, n_independent_components));
1055 
1056  TableIndices<rank_> indices;
1057 
1058  unsigned int remainder = i;
1059  for (int r=rank_-1; r>=0; --r)
1060  {
1061  indices[r] = (remainder % dim);
1062  remainder /= dim;
1063  }
1064  Assert (remainder == 0, ExcInternalError());
1065 
1066  return indices;
1067 }
1068 
1069 
1070 template <int rank_, int dim, typename Number>
1071 inline
1073 {
1074  for (unsigned int i=0; i<dim; ++i)
1075  values[i] = value_type();
1076 }
1077 
1078 
1079 template <int rank_, int dim, typename Number>
1080 inline
1081 std::size_t
1083 {
1084  return sizeof(Tensor<rank_,dim,Number>);
1085 }
1086 
1087 
1088 template <int rank_, int dim, typename Number>
1089 template <class Archive>
1090 inline
1091 void
1092 Tensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
1093 {
1094  ar &values;
1095 }
1096 
1097 
1098 /* ----------------- Non-member functions operating on tensors. ------------ */
1099 
1104 
1112 template <int rank_, int dim, typename Number>
1113 inline
1114 std::ostream &operator << (std::ostream &out, const Tensor<rank_,dim,Number> &p)
1115 {
1116  for (unsigned int i = 0; i < dim; ++i)
1117  {
1118  out << p[i];
1119  if (i != dim - 1)
1120  out << ' ';
1121  }
1122 
1123  return out;
1124 }
1125 
1126 
1133 template <int dim, typename Number>
1134 inline
1135 std::ostream &operator << (std::ostream &out, const Tensor<0,dim,Number> &p)
1136 {
1137  out << static_cast<const Number &>(p);
1138  return out;
1139 }
1140 
1141 
1143 
1147 
1148 
1149 #ifndef DEAL_II_WITH_CXX11
1150 template <typename T, typename U, int rank, int dim>
1151 struct ProductType<T,Tensor<rank,dim,U> >
1152 {
1154 };
1155 
1156 template <typename T, typename U, int rank, int dim>
1157 struct ProductType<Tensor<rank,dim,T>,U>
1158 {
1160 };
1161 #endif
1162 
1163 
1164 
1173 template <int dim, typename Number, typename Other>
1174 inline
1176 operator * (const Other object,
1177  const Tensor<0,dim,Number> &t)
1178 {
1179  return object * static_cast<const Number &>(t);
1180 }
1181 
1182 
1191 template <int dim, typename Number, typename Other>
1192 inline
1195  const Other object)
1196 {
1197  return static_cast<const Number &>(t) * object;
1198 }
1199 
1200 
1210 template <int dim, typename Number, typename OtherNumber>
1211 inline
1214  const Tensor<0, dim, OtherNumber> &src2)
1215 {
1216  return static_cast<const Number &>(src1) *
1217  static_cast<const OtherNumber &>(src2);
1218 }
1219 
1220 
1226 template <int dim, typename Number, typename OtherNumber>
1227 inline
1230  const OtherNumber factor)
1231 {
1232  return static_cast<Number>(t) / factor;
1233 }
1234 
1235 
1241 template <int dim, typename Number, typename OtherNumber>
1242 inline
1245 {
1246  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1247 }
1248 
1249 
1255 template <int dim, typename Number, typename OtherNumber>
1256 inline
1259 {
1260  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1261 }
1262 
1263 
1274 template <int rank, int dim,
1275  typename Number,
1276  typename OtherNumber>
1277 inline
1280  const OtherNumber factor)
1281 {
1282  // recurse over the base objects
1284  for (unsigned int d=0; d<dim; ++d)
1285  tt[d] = t[d] * factor;
1286  return tt;
1287 }
1288 
1289 
1300 template <int rank, int dim,
1301  typename Number,
1302  typename OtherNumber>
1303 inline
1305 operator * (const Number factor,
1307 {
1308  // simply forward to the operator above
1309  return t * factor;
1310 }
1311 
1312 
1320 template <int rank, int dim,
1321  typename Number,
1322  typename OtherNumber>
1323 inline
1326  const OtherNumber factor)
1327 {
1328  // recurse over the base objects
1330  for (unsigned int d=0; d<dim; ++d)
1331  tt[d] = t[d] / factor;
1332  return tt;
1333 }
1334 
1335 
1343 template <int rank, int dim, typename Number, typename OtherNumber>
1344 inline
1347 {
1349 
1350  for (unsigned int i=0; i<dim; ++i)
1351  tmp[i] += q[i];
1352 
1353  return tmp;
1354 }
1355 
1356 
1364 template <int rank, int dim, typename Number, typename OtherNumber>
1365 inline
1368 {
1370 
1371  for (unsigned int i=0; i<dim; ++i)
1372  tmp[i] -= q[i];
1373 
1374  return tmp;
1375 }
1376 
1377 
1379 
1383 
1384 
1408 template <int rank_1, int rank_2, int dim,
1409  typename Number, typename OtherNumber>
1410 inline DEAL_II_ALWAYS_INLINE
1411 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1414 {
1415  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
1416 
1417  TensorAccessors::internal::ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1418  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
1419  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
1420 
1421  return result;
1422 }
1423 
1424 
1454 template <int index_1, int index_2,
1455  int rank_1, int rank_2, int dim,
1456  typename Number, typename OtherNumber>
1457 inline
1458 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1461 {
1462  Assert(0 <= index_1 && index_1 < rank_1,
1463  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1464  Assert(0 <= index_2 && index_2 < rank_2,
1465  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1466 
1467  using namespace TensorAccessors;
1468  using namespace TensorAccessors::internal;
1469 
1470  // Reorder index_1 to the end of src1:
1471  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1472  reord_01 = reordered_index_view<index_1, rank_1>(src1);
1473 
1474  // Reorder index_2 to the end of src2:
1475  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1476  reord_02 = reordered_index_view<index_2, rank_2>(src2);
1477 
1478  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1479  result;
1480  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
1481  return result;
1482 }
1483 
1484 
1516 template <int index_1, int index_2, int index_3, int index_4,
1517  int rank_1, int rank_2, int dim,
1518  typename Number, typename OtherNumber>
1519 inline
1520 typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1523 {
1524  Assert(0 <= index_1 && index_1 < rank_1,
1525  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1526  Assert(0 <= index_3 && index_3 < rank_1,
1527  ExcMessage("The specified index_3 must lie within the range [0,rank_1)"));
1528  Assert(index_1 != index_3,
1529  ExcMessage("index_1 and index_3 must not be the same"));
1530  Assert(0 <= index_2 && index_2 < rank_2,
1531  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1532  Assert(0 <= index_4 && index_4 < rank_2,
1533  ExcMessage("The specified index_4 must lie within the range [0,rank_2)"));
1534  Assert(index_2 != index_4,
1535  ExcMessage("index_2 and index_4 must not be the same"));
1536 
1537  using namespace TensorAccessors;
1538  using namespace TensorAccessors::internal;
1539 
1540  // Reorder index_1 to the end of src1:
1541  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1542  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
1543 
1544  // Reorder index_2 to the end of src2:
1545  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1546  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
1547 
1548  // Now, reorder index_3 to the end of src1. We have to make sure to
1549  // preserve the orginial ordering: index_1 has been removed. If
1550  // index_3 > index_1, we have to use (index_3 - 1) instead:
1551  ReorderedIndexView<(index_3 < index_1 ? index_3 : index_3 - 1), rank_1, ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> > >
1552  reord_3 = TensorAccessors::reordered_index_view<index_3 < index_1 ? index_3 : index_3 - 1, rank_1>(reord_1);
1553 
1554  // Now, reorder index_4 to the end of src2. We have to make sure to
1555  // preserve the orginial ordering: index_2 has been removed. If
1556  // index_4 > index_2, we have to use (index_4 - 1) instead:
1557  ReorderedIndexView<(index_4 < index_2 ? index_4 : index_4 - 1), rank_2, ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> > >
1558  reord_4 = TensorAccessors::reordered_index_view<index_4 < index_2 ? index_4 : index_4 - 1, rank_2>(reord_2);
1559 
1560  typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1561  result;
1562  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
1563  return result;
1564 }
1565 
1566 
1580 template <int rank, int dim, typename Number, typename OtherNumber>
1581 inline
1584  const Tensor<rank, dim, OtherNumber> &right)
1585 {
1587  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
1588  return result;
1589 }
1590 
1591 
1610 template <template<int, int, typename> class TensorT1,
1611  template<int, int, typename> class TensorT2,
1612  template<int, int, typename> class TensorT3,
1613  int rank_1, int rank_2, int dim,
1614  typename T1, typename T2, typename T3>
1616 contract3 (const TensorT1<rank_1, dim, T1> &left,
1617  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
1618  const TensorT3<rank_2, dim, T3> &right)
1619 {
1621  return_type;
1622  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(
1623  left, middle, right);
1624 }
1625 
1626 
1638 template <int rank_1, int rank_2, int dim,
1639  typename Number, typename OtherNumber>
1640 inline
1644 {
1646  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
1647  return result;
1648 }
1649 
1650 
1652 
1656 
1657 
1669 template <int dim, typename Number>
1670 inline
1673 {
1674  Assert (dim==2, ExcInternalError());
1675 
1676  Tensor<1, dim, Number> result;
1677 
1678  result[0] = src[1];
1679  result[1] = -src[0];
1680 
1681  return result;
1682 }
1683 
1684 
1695 template <int dim, typename Number>
1696 inline
1699  const Tensor<1,dim,Number> &src2)
1700 {
1701  Assert (dim==3, ExcInternalError());
1702 
1703  Tensor<1, dim, Number> result;
1704 
1705  result[0] = src1[1]*src2[2] - src1[2]*src2[1];
1706  result[1] = src1[2]*src2[0] - src1[0]*src2[2];
1707  result[2] = src1[0]*src2[1] - src1[1]*src2[0];
1708 
1709  return result;
1710 }
1711 
1712 
1714 
1718 
1719 
1726 template <int dim, typename Number>
1727 inline
1729 {
1730  // Compute the determinant using the Laplace expansion of the
1731  // determinant. We expand along the last row.
1732  Number det = Number();
1733 
1734  for (unsigned int k=0; k<dim; ++k)
1735  {
1736  Tensor<2,dim-1,Number> minor;
1737  for (unsigned int i=0; i<dim-1; ++i)
1738  for (unsigned int j=0; j<dim-1; ++j)
1739  minor[i][j] = t[i][j<k ? j : j+1];
1740 
1741  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
1742 
1743  det += t[dim-1][k] * cofactor;
1744  }
1745 
1746  return ((dim % 2 == 0) ? 1. : -1.) * det;
1747 }
1748 
1754 template <typename Number>
1755 inline
1757 {
1758  return t[0][0];
1759 }
1760 
1761 
1769 template <int dim, typename Number>
1770 Number trace (const Tensor<2,dim,Number> &d)
1771 {
1772  Number t=d[0][0];
1773  for (unsigned int i=1; i<dim; ++i)
1774  t += d[i][i];
1775  return t;
1776 }
1777 
1778 
1788 template <int dim, typename Number>
1789 inline
1792 {
1793  Number return_tensor [dim][dim];
1794 
1795  // if desired, take over the
1796  // inversion of a 4x4 tensor
1797  // from the FullMatrix
1798  AssertThrow (false, ExcNotImplemented());
1799 
1800  return Tensor<2,dim,Number>(return_tensor);
1801 }
1802 
1803 
1804 #ifndef DOXYGEN
1805 
1806 template <typename Number>
1807 inline
1809 invert (const Tensor<2,1,Number> &t)
1810 {
1811  Number return_tensor [1][1];
1812 
1813  return_tensor[0][0] = 1.0/t[0][0];
1814 
1815  return Tensor<2,1,Number>(return_tensor);
1816 }
1817 
1818 
1819 template <typename Number>
1820 inline
1822 invert (const Tensor<2,2,Number> &t)
1823 {
1824  Tensor<2,2,Number> return_tensor;
1825 
1826  // this is Maple output,
1827  // thus a bit unstructured
1828  const Number inv_det_t = 1.0/(t[0][0]*t[1][1]-t[1][0]*t[0][1]);
1829  return_tensor[0][0] = t[1][1];
1830  return_tensor[0][1] = -t[0][1];
1831  return_tensor[1][0] = -t[1][0];
1832  return_tensor[1][1] = t[0][0];
1833  return_tensor *= inv_det_t;
1834 
1835  return return_tensor;
1836 }
1837 
1838 
1839 template <typename Number>
1840 inline
1842 invert (const Tensor<2,3,Number> &t)
1843 {
1844  Tensor<2,3,Number> return_tensor;
1845 
1846  const Number t4 = t[0][0]*t[1][1],
1847  t6 = t[0][0]*t[1][2],
1848  t8 = t[0][1]*t[1][0],
1849  t00 = t[0][2]*t[1][0],
1850  t01 = t[0][1]*t[2][0],
1851  t04 = t[0][2]*t[2][0],
1852  inv_det_t = 1.0/(t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
1853  t00*t[2][1]+t01*t[1][2]-t04*t[1][1]);
1854  return_tensor[0][0] = t[1][1]*t[2][2]-t[1][2]*t[2][1];
1855  return_tensor[0][1] = t[0][2]*t[2][1]-t[0][1]*t[2][2];
1856  return_tensor[0][2] = t[0][1]*t[1][2]-t[0][2]*t[1][1];
1857  return_tensor[1][0] = t[1][2]*t[2][0]-t[1][0]*t[2][2];
1858  return_tensor[1][1] = t[0][0]*t[2][2]-t04;
1859  return_tensor[1][2] = t00-t6;
1860  return_tensor[2][0] = t[1][0]*t[2][1]-t[1][1]*t[2][0];
1861  return_tensor[2][1] = t01-t[0][0]*t[2][1];
1862  return_tensor[2][2] = t4-t8;
1863  return_tensor *= inv_det_t;
1864 
1865  return return_tensor;
1866 }
1867 
1868 #endif /* DOXYGEN */
1869 
1870 
1877 template <int dim, typename Number>
1878 inline
1881 {
1883  for (unsigned int i=0; i<dim; ++i)
1884  {
1885  tt[i][i] = t[i][i];
1886  for (unsigned int j=i+1; j<dim; ++j)
1887  {
1888  tt[i][j] = t[j][i];
1889  tt[j][i] = t[i][j];
1890  };
1891  }
1892  return tt;
1893 }
1894 
1895 
1909 template <int dim, typename Number>
1910 inline
1913 {
1914  return determinant(t)*invert(t);
1915 }
1916 
1917 
1932 template <int dim, typename Number>
1933 inline
1936 {
1937  return transpose(adjugate(t));
1938 }
1939 
1940 
1948 template <int dim, typename Number>
1949 inline
1950 double
1952 {
1953  double max = 0;
1954  for (unsigned int j=0; j<dim; ++j)
1955  {
1956  double sum = 0;
1957  for (unsigned int i=0; i<dim; ++i)
1958  sum += std::fabs(t[i][j]);
1959 
1960  if (sum > max)
1961  max = sum;
1962  }
1963 
1964  return max;
1965 }
1966 
1967 
1975 template <int dim, typename Number>
1976 inline
1977 double
1979 {
1980  double max = 0;
1981  for (unsigned int i=0; i<dim; ++i)
1982  {
1983  double sum = 0;
1984  for (unsigned int j=0; j<dim; ++j)
1985  sum += std::fabs(t[i][j]);
1986 
1987  if (sum > max)
1988  max = sum;
1989  }
1990 
1991  return max;
1992 }
1993 
1995 
1996 DEAL_II_NAMESPACE_CLOSE
1997 
1998 // include deprecated non-member functions operating on Tensor
1999 #include <deal.II/base/tensor_deprecated.h>
2000 
2001 #endif
numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:115
Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:1770
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1146
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
Definition: tensor.h:1038
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:355
SymmetricTensor< rank, dim, Number > operator/(const SymmetricTensor< rank, dim, Number > &t, const Number factor)
static const unsigned int n_independent_components
Definition: tensor.h:341
static std::size_t memory_consumption()
Definition: tensor.h:1082
Tensor< rank_, dim, Number > & operator-=(const Tensor< rank_, dim, OtherNumber > &)
Definition: tensor.h:942
Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1912
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:1459
static const unsigned int rank
Definition: tensor.h:334
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:991
#define AssertThrow(cond, exc)
Definition: exceptions.h:369
static real_type abs(const number &x)
Definition: numbers.h:322
bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:920
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
Definition: point.h:89
Tensor< rank_-1, dim, Number > values[(dim!=0)?dim:1]
Definition: tensor.h:543
void unroll(Vector< OtherNumber > &result) const
Definition: tensor.h:1014
Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1935
Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1880
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:1616
Tensor< rank_, dim, Number > & operator+=(const Tensor< rank_, dim, OtherNumber > &)
Definition: tensor.h:930
Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:1791
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:1672
bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:893
#define Assert(cond, exc)
Definition: exceptions.h:313
ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:1583
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: tensor.h:1027
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:537
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:1092
double linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1978
std::size_t size() const
Tensor()
Definition: tensor.h:745
double l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1951
Tensor< rank_, dim, Number > & operator/=(const OtherNumber factor)
Definition: tensor.h:966
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Definition: tensor.h:1051
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:1521
value_type & operator[](const unsigned int i)
Definition: tensor.h:826
Tensor< 1, dim, Number > cross_product_3d(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:1698
Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1728
static const unsigned int dimension
Definition: tensor.h:329
Tensor< rank_, dim, Number > operator-() const
Definition: tensor.h:977
Definition: mpi.h:41
Tensor< rank_-1, dim, Number >::tensor_type value_type
Definition: tensor.h:348
static::ExceptionBase & ExcNotImplemented()
Tensor< rank_, dim, Number > & operator*=(const OtherNumber factor)
Definition: tensor.h:954
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:1642
numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1000
void clear()
Definition: tensor.h:1072
static::ExceptionBase & ExcInternalError()
Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:1756
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const