Reference documentation for deal.II version Git 0ec6fa4 2017-05-22 12:08:17 +0200
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 
31 DEAL_II_NAMESPACE_OPEN
32 
33 // Forward declarations:
34 
35 template <int dim, typename Number> class Point;
36 template <int rank_, int dim, typename Number = double> class Tensor;
37 template <typename Number> class Vector;
38 
39 #ifndef DOXYGEN
40 // Overload invalid tensor types of negative rank that come up during
41 // overload resolution of operator* and related contraction variants.
42 template <int dim, typename Number>
43 class Tensor<-2, dim, Number>
44 {
45 };
46 
47 template <int dim, typename Number>
48 class Tensor<-1, dim, Number>
49 {
50 };
51 #endif /* DOXYGEN */
52 
53 
84 template <int dim, typename Number>
85 class Tensor<0,dim,Number>
86 {
87 public:
96  static const unsigned int dimension = dim;
97 
101  static const unsigned int rank = 0;
102 
106  static const unsigned int n_independent_components = 1;
107 
117 
122  typedef Number value_type;
123 
129  typedef Number array_type;
130 
136  DEAL_II_CUDA_HOST_DEV Tensor ();
137 
143  template <typename OtherNumber>
144  Tensor (const Tensor<0,dim,OtherNumber> &initializer);
145 
149  template <typename OtherNumber>
150  Tensor (const OtherNumber initializer);
151 
161  DEAL_II_CUDA_HOST_DEV operator Number &();
162 
171  DEAL_II_CUDA_HOST_DEV operator const Number &() const;
172 
178  template <typename OtherNumber>
180 
184  template<typename OtherNumber>
185  bool operator == (const Tensor<0,dim,OtherNumber> &rhs) const;
186 
190  template<typename OtherNumber>
191  bool operator != (const Tensor<0,dim,OtherNumber> &rhs) const;
192 
196  template<typename OtherNumber>
198 
202  template<typename OtherNumber>
204 
210  template<typename OtherNumber>
211  DEAL_II_CUDA_HOST_DEV 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 
251  DEAL_II_CUDA_HOST_DEV real_type norm_square () const;
252 
257  template <class Archive>
258  void serialize(Archive &ar, const unsigned int version);
259 
264  typedef Number tensor_type;
265 
266 private:
270  Number value;
271 
275  template <typename OtherNumber>
277  unsigned int &start_index) const;
278 
282  template <int, int, typename> friend class Tensor;
283 };
284 
285 
286 
328 template <int rank_, int dim, typename Number>
329 class Tensor
330 {
331 public:
340  static const unsigned int dimension = dim;
341 
345  static const unsigned int rank = rank_;
346 
351  static const unsigned int
353 
359  typedef typename Tensor<rank_-1,dim,Number>::tensor_type value_type;
360 
365  typedef typename Tensor<rank_-1,dim,Number>::array_type
366  array_type[(dim != 0) ? dim : 1];
367  // ... avoid a compiler warning in case of dim == 0 and ensure that the
368  // array always has positive size.
369 
375  DEAL_II_CUDA_HOST_DEV Tensor ();
376 
380  Tensor (const array_type &initializer);
381 
387  template <typename OtherNumber>
388  Tensor (const Tensor<rank_,dim,OtherNumber> &initializer);
389 
393  template <typename OtherNumber>
394  Tensor (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer);
395 
399  template <typename OtherNumber>
400  operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const;
401 
407  DEAL_II_CUDA_HOST_DEV value_type &operator [] (const unsigned int i);
408 
414  DEAL_II_CUDA_HOST_DEV const value_type &operator[](const unsigned int i) const;
415 
419  const Number &operator [] (const TableIndices<rank_> &indices) const;
420 
424  Number &operator [] (const TableIndices<rank_> &indices);
425 
431  template <typename OtherNumber>
433 
440  Tensor<rank_,dim,Number> &operator = (const Number d);
441 
445  template <typename OtherNumber>
446  bool operator == (const Tensor<rank_,dim,OtherNumber> &) const;
447 
451  template <typename OtherNumber>
452  bool operator != (const Tensor<rank_,dim,OtherNumber> &) const;
453 
457  template <typename OtherNumber>
459 
463  template <typename OtherNumber>
465 
472  template <typename OtherNumber>
473  DEAL_II_CUDA_HOST_DEV 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 
514  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type norm_square() const;
515 
523  template <typename OtherNumber>
524  void unroll (Vector<OtherNumber> &result) const;
525 
530  static
531  unsigned int
533 
538  static
540 
545  static std::size_t memory_consumption ();
546 
551  template <class Archive>
552  void serialize(Archive &ar, const unsigned int version);
553 
559 
560 private:
564  Tensor<rank_-1, dim, Number> values[(dim != 0) ? dim : 1];
565  // ... avoid a compiler warning in case of dim == 0 and ensure that the
566  // array always has positive size.
567 
571  template <typename OtherNumber>
573  unsigned int &start_index) const;
574 
578  template <int, int, typename> friend class Tensor;
579 
584  friend class Point<dim,Number>;
585 };
586 
587 
588 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
589 
590 
591 template <int dim,typename Number>
592 inline
593 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::Tensor ()
594  : value()
595 {
596 }
597 
598 
599 template <int dim, typename Number>
600 template <typename OtherNumber>
601 inline
602 Tensor<0,dim,Number>::Tensor (const OtherNumber initializer)
603 {
604  value = initializer;
605 }
606 
607 
608 template <int dim, typename Number>
609 template <typename OtherNumber>
610 inline
612 {
613  value = p.value;
614 }
615 
616 
617 template <int dim, typename Number>
618 inline
619 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::operator Number &()
620 {
621  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
622  return value;
623 }
624 
625 
626 template <int dim, typename Number>
627 inline
628 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::operator const Number &() const
629 {
630  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
631  return value;
632 }
633 
634 
635 template <int dim, typename Number>
636 template <typename OtherNumber>
637 inline
639 {
640  value = p.value;
641  return *this;
642 }
643 
644 
645 template <int dim, typename Number>
646 template <typename OtherNumber>
647 inline
649 {
650  return (value == p.value);
651 }
652 
653 
654 template <int dim, typename Number>
655 template <typename OtherNumber>
656 inline
658 {
659  return !((*this) == p);
660 }
661 
662 
663 template <int dim, typename Number>
664 template <typename OtherNumber>
665 inline
667 {
668  value += p.value;
669  return *this;
670 }
671 
672 
673 template <int dim, typename Number>
674 template <typename OtherNumber>
675 inline
677 {
678  value -= p.value;
679  return *this;
680 }
681 
682 
683 template <int dim, typename Number>
684 template <typename OtherNumber>
685 inline
686 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator *= (const OtherNumber s)
687 {
688  value *= s;
689  return *this;
690 }
691 
692 
693 template <int dim, typename Number>
694 template <typename OtherNumber>
695 inline
697 {
698  value /= s;
699  return *this;
700 }
701 
702 
703 template <int dim, typename Number>
704 inline
706 {
707  return -value;
708 }
709 
710 
711 template <int dim, typename Number>
712 inline
715 {
716  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
717  return numbers::NumberTraits<Number>::abs (value);
718 }
719 
720 
721 template <int dim, typename Number>
722 inline
724 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::norm_square () const
725 {
726  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
728 }
729 
730 
731 template <int dim, typename Number>
732 template <typename OtherNumber>
733 inline
734 void
736  unsigned int &index) const
737 {
738  Assert(dim != 0, ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
739  result[index] = value;
740  ++index;
741 }
742 
743 
744 template <int dim, typename Number>
745 inline
747 {
748  value = value_type();
749 }
750 
751 
752 template <int dim, typename Number>
753 template <class Archive>
754 inline
755 void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int)
756 {
757  ar &value;
758 }
759 
760 
761 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
762 
763 
764 template <int rank_, int dim, typename Number>
765 inline
766 DEAL_II_CUDA_HOST_DEV Tensor<rank_,dim,Number>::Tensor ()
767 {
768  // All members of the c-style array values are already default initialized
769  // and thus all values are already set to zero recursively.
770 }
771 
772 
773 template <int rank_, int dim, typename Number>
774 inline
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 {
787  for (unsigned int i=0; i!=dim; ++i)
788  values[i] = initializer[i];
789 }
790 
791 
792 template <int rank_, int dim, typename Number>
793 template <typename OtherNumber>
794 inline
796 (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer)
797 {
798  for (unsigned int i=0; i<dim; ++i)
799  values[i] = initializer[i];
800 }
801 
802 
803 template <int rank_, int dim, typename Number>
804 template <typename OtherNumber>
805 inline
807 operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const
808 {
809  return Tensor<1,dim,Tensor<rank_-1,dim,Number> > (values);
810 }
811 
812 
813 
814 namespace internal
815 {
816  namespace TensorSubscriptor
817  {
818  template <typename ArrayElementType, int dim>
819  inline DEAL_II_ALWAYS_INLINE
820  DEAL_II_CUDA_HOST_DEV
821  ArrayElementType &
822  subscript (ArrayElementType *values,
823  const unsigned int i,
825  {
826  Assert (i<dim, ExcIndexRange(i, 0, dim));
827  return values[i];
828  }
829 
830 
831  template <typename ArrayElementType>
832  DEAL_II_CUDA_HOST_DEV
833  ArrayElementType &
834  subscript (ArrayElementType *,
835  const unsigned int,
837  {
838  Assert(false, ExcMessage("Cannot access elements of an object of type Tensor<rank,0,Number>."));
839  static ArrayElementType t;
840  return t;
841  }
842  }
843 }
844 
845 
846 template <int rank_, int dim, typename Number>
847 inline DEAL_II_ALWAYS_INLINE
848 DEAL_II_CUDA_HOST_DEV
851 {
852  return ::internal::TensorSubscriptor::subscript(values, i, ::internal::int2type<dim>());
853 }
854 
855 
856 template <int rank_, int dim, typename Number>
857 inline DEAL_II_ALWAYS_INLINE
858 DEAL_II_CUDA_HOST_DEV
860 Tensor<rank_,dim,Number>::operator[] (const unsigned int i) const
861 {
862  return ::internal::TensorSubscriptor::subscript(values, i, ::internal::int2type<dim>());
863 }
864 
865 
866 template <int rank_, int dim, typename Number>
867 inline
868 const Number &
870 {
871  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
872 
873  return TensorAccessors::extract<rank_>(*this, indices);
874 }
875 
876 
877 template <int rank_, int dim, typename Number>
878 inline
879 Number &
881 {
882  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
883 
884  return TensorAccessors::extract<rank_>(*this, indices);
885 }
886 
887 
888 template <int rank_, int dim, typename Number>
889 template <typename OtherNumber>
890 inline
893 {
894  if (dim > 0)
895  std::copy (&t.values[0], &t.values[0]+dim, &values[0]);
896  return *this;
897 }
898 
899 
900 template <int rank_, int dim, typename Number>
901 inline
904 {
905  Assert (d == Number(), ExcMessage ("Only assignment with zero is allowed"));
906  (void) d;
907 
908  for (unsigned int i=0; i<dim; ++i)
909  values[i] = Number();
910  return *this;
911 }
912 
913 
914 template <int rank_, int dim, typename Number>
915 template <typename OtherNumber>
916 inline
917 bool
919 {
920  for (unsigned int i=0; i<dim; ++i)
921  if (values[i] != p.values[i])
922  return false;
923  return true;
924 }
925 
926 
927 // At some places in the library, we have Point<0> for formal reasons
928 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
929 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
930 // in the above function that the loop end check always fails, we
931 // implement this function here
932 template <>
933 template <>
934 inline
936 {
937  return true;
938 }
939 
940 
941 template <int rank_, int dim, typename Number>
942 template <typename OtherNumber>
943 inline
944 bool
946 {
947  return !((*this) == p);
948 }
949 
950 
951 template <int rank_, int dim, typename Number>
952 template <typename OtherNumber>
953 inline
956 {
957  for (unsigned int i=0; i<dim; ++i)
958  values[i] += p.values[i];
959  return *this;
960 }
961 
962 
963 template <int rank_, int dim, typename Number>
964 template <typename OtherNumber>
965 inline
968 {
969  for (unsigned int i=0; i<dim; ++i)
970  values[i] -= p.values[i];
971  return *this;
972 }
973 
974 
975 template <int rank_, int dim, typename Number>
976 template <typename OtherNumber>
977 inline
978 DEAL_II_CUDA_HOST_DEV
981 {
982  for (unsigned int i=0; i<dim; ++i)
983  values[i] *= s;
984  return *this;
985 }
986 
987 
988 template <int rank_, int dim, typename Number>
989 template <typename OtherNumber>
990 inline
993 {
994  for (unsigned int i=0; i<dim; ++i)
995  values[i] /= s;
996  return *this;
997 }
998 
999 
1000 template <int rank_, int dim, typename Number>
1001 inline
1004 {
1006 
1007  for (unsigned int i=0; i<dim; ++i)
1008  tmp.values[i] = -values[i];
1009 
1010  return tmp;
1011 }
1012 
1013 
1014 template <int rank_, int dim, typename Number>
1015 inline
1018 {
1019  return std::sqrt (norm_square());
1020 }
1021 
1022 
1023 template <int rank_, int dim, typename Number>
1024 inline
1025 DEAL_II_CUDA_HOST_DEV
1028 {
1030  for (unsigned int i=0; i<dim; ++i)
1031  s += values[i].norm_square();
1032 
1033  return s;
1034 }
1035 
1036 
1037 template <int rank_, int dim, typename Number>
1038 template <typename OtherNumber>
1039 inline
1040 void
1042 {
1043  AssertDimension (result.size(),(Utilities::fixed_power<rank_, unsigned int>(dim)));
1044 
1045  unsigned int index = 0;
1046  unroll_recursion (result, index);
1047 }
1048 
1049 
1050 template <int rank_, int dim, typename Number>
1051 template <typename OtherNumber>
1052 inline
1053 void
1055  unsigned int &index) const
1056 {
1057  for (unsigned int i=0; i<dim; ++i)
1058  values[i].unroll_recursion(result, index);
1059 }
1060 
1061 
1062 template <int rank_, int dim, typename Number>
1063 inline
1064 unsigned int
1066 {
1067  unsigned int index = 0;
1068  for (int r = 0; r < rank_; ++r)
1069  index = index * dim + indices[r];
1070 
1071  return index;
1072 }
1073 
1074 
1075 template <int rank_, int dim, typename Number>
1076 inline
1079 {
1080  Assert (i < n_independent_components,
1081  ExcIndexRange (i, 0, n_independent_components));
1082 
1083  TableIndices<rank_> indices;
1084 
1085  unsigned int remainder = i;
1086  for (int r=rank_-1; r>=0; --r)
1087  {
1088  indices[r] = (remainder % dim);
1089  remainder /= dim;
1090  }
1091  Assert (remainder == 0, ExcInternalError());
1092 
1093  return indices;
1094 }
1095 
1096 
1097 template <int rank_, int dim, typename Number>
1098 inline
1100 {
1101  for (unsigned int i=0; i<dim; ++i)
1102  values[i] = value_type();
1103 }
1104 
1105 
1106 template <int rank_, int dim, typename Number>
1107 inline
1108 std::size_t
1110 {
1111  return sizeof(Tensor<rank_,dim,Number>);
1112 }
1113 
1114 
1115 template <int rank_, int dim, typename Number>
1116 template <class Archive>
1117 inline
1118 void
1119 Tensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
1120 {
1121  ar &values;
1122 }
1123 
1124 
1125 /* ----------------- Non-member functions operating on tensors. ------------ */
1126 
1131 
1139 template <int rank_, int dim, typename Number>
1140 inline
1141 std::ostream &operator << (std::ostream &out, const Tensor<rank_,dim,Number> &p)
1142 {
1143  for (unsigned int i = 0; i < dim; ++i)
1144  {
1145  out << p[i];
1146  if (i != dim - 1)
1147  out << ' ';
1148  }
1149 
1150  return out;
1151 }
1152 
1153 
1160 template <int dim, typename Number>
1161 inline
1162 std::ostream &operator << (std::ostream &out, const Tensor<0,dim,Number> &p)
1163 {
1164  out << static_cast<const Number &>(p);
1165  return out;
1166 }
1167 
1168 
1170 
1174 
1175 
1184 template <int dim, typename Number, typename Other>
1185 inline
1187 operator * (const Other object,
1188  const Tensor<0,dim,Number> &t)
1189 {
1190  return object * static_cast<const Number &>(t);
1191 }
1192 
1193 
1194 
1203 template <int dim, typename Number, typename Other>
1204 inline
1207  const Other object)
1208 {
1209  return static_cast<const Number &>(t) * object;
1210 }
1211 
1212 
1222 template <int dim, typename Number, typename OtherNumber>
1223 inline
1226  const Tensor<0, dim, OtherNumber> &src2)
1227 {
1228  return static_cast<const Number &>(src1) *
1229  static_cast<const OtherNumber &>(src2);
1230 }
1231 
1232 
1238 template <int dim, typename Number, typename OtherNumber>
1239 inline
1242  const OtherNumber factor)
1243 {
1244  return static_cast<Number>(t) / factor;
1245 }
1246 
1247 
1253 template <int dim, typename Number, typename OtherNumber>
1254 inline
1257 {
1258  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1259 }
1260 
1261 
1267 template <int dim, typename Number, typename OtherNumber>
1268 inline
1271 {
1272  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1273 }
1274 
1275 
1286 template <int rank, int dim,
1287  typename Number,
1288  typename OtherNumber>
1289 inline
1292  const OtherNumber factor)
1293 {
1294  // recurse over the base objects
1296  for (unsigned int d=0; d<dim; ++d)
1297  tt[d] = t[d] * factor;
1298  return tt;
1299 }
1300 
1301 
1312 template <int rank, int dim,
1313  typename Number,
1314  typename OtherNumber>
1315 inline
1317 operator * (const Number factor,
1319 {
1320  // simply forward to the operator above
1321  return t * factor;
1322 }
1323 
1324 
1332 template <int rank, int dim,
1333  typename Number,
1334  typename OtherNumber>
1335 inline
1338  const OtherNumber factor)
1339 {
1340  // recurse over the base objects
1342  for (unsigned int d=0; d<dim; ++d)
1343  tt[d] = t[d] / factor;
1344  return tt;
1345 }
1346 
1347 
1355 template <int rank, int dim, typename Number, typename OtherNumber>
1356 inline
1359 {
1361 
1362  for (unsigned int i=0; i<dim; ++i)
1363  tmp[i] += q[i];
1364 
1365  return tmp;
1366 }
1367 
1368 
1376 template <int rank, int dim, typename Number, typename OtherNumber>
1377 inline
1380 {
1382 
1383  for (unsigned int i=0; i<dim; ++i)
1384  tmp[i] -= q[i];
1385 
1386  return tmp;
1387 }
1388 
1389 
1391 
1395 
1396 
1420 template <int rank_1, int rank_2, int dim,
1421  typename Number, typename OtherNumber>
1422 inline DEAL_II_ALWAYS_INLINE
1423 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1426 {
1427  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
1428 
1429  TensorAccessors::internal::ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1430  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
1431  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
1432 
1433  return result;
1434 }
1435 
1436 
1466 template <int index_1, int index_2,
1467  int rank_1, int rank_2, int dim,
1468  typename Number, typename OtherNumber>
1469 inline
1470 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1473 {
1474  Assert(0 <= index_1 && index_1 < rank_1,
1475  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1476  Assert(0 <= index_2 && index_2 < rank_2,
1477  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1478 
1479  using namespace TensorAccessors;
1480  using namespace TensorAccessors::internal;
1481 
1482  // Reorder index_1 to the end of src1:
1483  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1484  reord_01 = reordered_index_view<index_1, rank_1>(src1);
1485 
1486  // Reorder index_2 to the end of src2:
1487  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1488  reord_02 = reordered_index_view<index_2, rank_2>(src2);
1489 
1490  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1491  result;
1492  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
1493  return result;
1494 }
1495 
1496 
1528 template <int index_1, int index_2, int index_3, int index_4,
1529  int rank_1, int rank_2, int dim,
1530  typename Number, typename OtherNumber>
1531 inline
1532 typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1535 {
1536  Assert(0 <= index_1 && index_1 < rank_1,
1537  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1538  Assert(0 <= index_3 && index_3 < rank_1,
1539  ExcMessage("The specified index_3 must lie within the range [0,rank_1)"));
1540  Assert(index_1 != index_3,
1541  ExcMessage("index_1 and index_3 must not be the same"));
1542  Assert(0 <= index_2 && index_2 < rank_2,
1543  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1544  Assert(0 <= index_4 && index_4 < rank_2,
1545  ExcMessage("The specified index_4 must lie within the range [0,rank_2)"));
1546  Assert(index_2 != index_4,
1547  ExcMessage("index_2 and index_4 must not be the same"));
1548 
1549  using namespace TensorAccessors;
1550  using namespace TensorAccessors::internal;
1551 
1552  // Reorder index_1 to the end of src1:
1553  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1554  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
1555 
1556  // Reorder index_2 to the end of src2:
1557  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1558  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
1559 
1560  // Now, reorder index_3 to the end of src1. We have to make sure to
1561  // preserve the orginial ordering: index_1 has been removed. If
1562  // index_3 > index_1, we have to use (index_3 - 1) instead:
1563  ReorderedIndexView<(index_3 < index_1 ? index_3 : index_3 - 1), rank_1, ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> > >
1564  reord_3 = TensorAccessors::reordered_index_view<index_3 < index_1 ? index_3 : index_3 - 1, rank_1>(reord_1);
1565 
1566  // Now, reorder index_4 to the end of src2. We have to make sure to
1567  // preserve the orginial ordering: index_2 has been removed. If
1568  // index_4 > index_2, we have to use (index_4 - 1) instead:
1569  ReorderedIndexView<(index_4 < index_2 ? index_4 : index_4 - 1), rank_2, ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> > >
1570  reord_4 = TensorAccessors::reordered_index_view<index_4 < index_2 ? index_4 : index_4 - 1, rank_2>(reord_2);
1571 
1572  typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1573  result;
1574  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
1575  return result;
1576 }
1577 
1578 
1592 template <int rank, int dim, typename Number, typename OtherNumber>
1593 inline
1596  const Tensor<rank, dim, OtherNumber> &right)
1597 {
1599  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
1600  return result;
1601 }
1602 
1603 
1622 template <template<int, int, typename> class TensorT1,
1623  template<int, int, typename> class TensorT2,
1624  template<int, int, typename> class TensorT3,
1625  int rank_1, int rank_2, int dim,
1626  typename T1, typename T2, typename T3>
1628 contract3 (const TensorT1<rank_1, dim, T1> &left,
1629  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
1630  const TensorT3<rank_2, dim, T3> &right)
1631 {
1633  return_type;
1634  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(
1635  left, middle, right);
1636 }
1637 
1638 
1650 template <int rank_1, int rank_2, int dim,
1651  typename Number, typename OtherNumber>
1652 inline
1656 {
1658  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
1659  return result;
1660 }
1661 
1662 
1664 
1668 
1669 
1681 template <int dim, typename Number>
1682 inline
1685 {
1686  Assert (dim==2, ExcInternalError());
1687 
1688  Tensor<1, dim, Number> result;
1689 
1690  result[0] = src[1];
1691  result[1] = -src[0];
1692 
1693  return result;
1694 }
1695 
1696 
1707 template <int dim, typename Number>
1708 inline
1711  const Tensor<1,dim,Number> &src2)
1712 {
1713  Assert (dim==3, ExcInternalError());
1714 
1715  Tensor<1, dim, Number> result;
1716 
1717  result[0] = src1[1]*src2[2] - src1[2]*src2[1];
1718  result[1] = src1[2]*src2[0] - src1[0]*src2[2];
1719  result[2] = src1[0]*src2[1] - src1[1]*src2[0];
1720 
1721  return result;
1722 }
1723 
1724 
1726 
1730 
1731 
1738 template <int dim, typename Number>
1739 inline
1741 {
1742  // Compute the determinant using the Laplace expansion of the
1743  // determinant. We expand along the last row.
1744  Number det = Number();
1745 
1746  for (unsigned int k=0; k<dim; ++k)
1747  {
1748  Tensor<2,dim-1,Number> minor;
1749  for (unsigned int i=0; i<dim-1; ++i)
1750  for (unsigned int j=0; j<dim-1; ++j)
1751  minor[i][j] = t[i][j<k ? j : j+1];
1752 
1753  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
1754 
1755  det += t[dim-1][k] * cofactor;
1756  }
1757 
1758  return ((dim % 2 == 0) ? 1. : -1.) * det;
1759 }
1760 
1766 template <typename Number>
1767 inline
1769 {
1770  return t[0][0];
1771 }
1772 
1773 
1781 template <int dim, typename Number>
1782 Number trace (const Tensor<2,dim,Number> &d)
1783 {
1784  Number t=d[0][0];
1785  for (unsigned int i=1; i<dim; ++i)
1786  t += d[i][i];
1787  return t;
1788 }
1789 
1790 
1800 template <int dim, typename Number>
1801 inline
1804 {
1805  Number return_tensor [dim][dim];
1806 
1807  // if desired, take over the
1808  // inversion of a 4x4 tensor
1809  // from the FullMatrix
1810  AssertThrow (false, ExcNotImplemented());
1811 
1812  return Tensor<2,dim,Number>(return_tensor);
1813 }
1814 
1815 
1816 #ifndef DOXYGEN
1817 
1818 template <typename Number>
1819 inline
1821 invert (const Tensor<2,1,Number> &t)
1822 {
1823  Number return_tensor [1][1];
1824 
1825  return_tensor[0][0] = 1.0/t[0][0];
1826 
1827  return Tensor<2,1,Number>(return_tensor);
1828 }
1829 
1830 
1831 template <typename Number>
1832 inline
1834 invert (const Tensor<2,2,Number> &t)
1835 {
1836  Tensor<2,2,Number> return_tensor;
1837 
1838  // this is Maple output,
1839  // thus a bit unstructured
1840  const Number inv_det_t = 1.0/(t[0][0]*t[1][1]-t[1][0]*t[0][1]);
1841  return_tensor[0][0] = t[1][1];
1842  return_tensor[0][1] = -t[0][1];
1843  return_tensor[1][0] = -t[1][0];
1844  return_tensor[1][1] = t[0][0];
1845  return_tensor *= inv_det_t;
1846 
1847  return return_tensor;
1848 }
1849 
1850 
1851 template <typename Number>
1852 inline
1854 invert (const Tensor<2,3,Number> &t)
1855 {
1856  Tensor<2,3,Number> return_tensor;
1857 
1858  const Number t4 = t[0][0]*t[1][1],
1859  t6 = t[0][0]*t[1][2],
1860  t8 = t[0][1]*t[1][0],
1861  t00 = t[0][2]*t[1][0],
1862  t01 = t[0][1]*t[2][0],
1863  t04 = t[0][2]*t[2][0],
1864  inv_det_t = 1.0/(t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
1865  t00*t[2][1]+t01*t[1][2]-t04*t[1][1]);
1866  return_tensor[0][0] = t[1][1]*t[2][2]-t[1][2]*t[2][1];
1867  return_tensor[0][1] = t[0][2]*t[2][1]-t[0][1]*t[2][2];
1868  return_tensor[0][2] = t[0][1]*t[1][2]-t[0][2]*t[1][1];
1869  return_tensor[1][0] = t[1][2]*t[2][0]-t[1][0]*t[2][2];
1870  return_tensor[1][1] = t[0][0]*t[2][2]-t04;
1871  return_tensor[1][2] = t00-t6;
1872  return_tensor[2][0] = t[1][0]*t[2][1]-t[1][1]*t[2][0];
1873  return_tensor[2][1] = t01-t[0][0]*t[2][1];
1874  return_tensor[2][2] = t4-t8;
1875  return_tensor *= inv_det_t;
1876 
1877  return return_tensor;
1878 }
1879 
1880 #endif /* DOXYGEN */
1881 
1882 
1889 template <int dim, typename Number>
1890 inline
1893 {
1895  for (unsigned int i=0; i<dim; ++i)
1896  {
1897  tt[i][i] = t[i][i];
1898  for (unsigned int j=i+1; j<dim; ++j)
1899  {
1900  tt[i][j] = t[j][i];
1901  tt[j][i] = t[i][j];
1902  };
1903  }
1904  return tt;
1905 }
1906 
1907 
1921 template <int dim, typename Number>
1922 inline
1925 {
1926  return determinant(t)*invert(t);
1927 }
1928 
1929 
1944 template <int dim, typename Number>
1945 inline
1948 {
1949  return transpose(adjugate(t));
1950 }
1951 
1952 
1960 template <int dim, typename Number>
1961 inline
1962 double
1964 {
1965  double max = 0;
1966  for (unsigned int j=0; j<dim; ++j)
1967  {
1968  double sum = 0;
1969  for (unsigned int i=0; i<dim; ++i)
1970  sum += std::fabs(t[i][j]);
1971 
1972  if (sum > max)
1973  max = sum;
1974  }
1975 
1976  return max;
1977 }
1978 
1979 
1987 template <int dim, typename Number>
1988 inline
1989 double
1991 {
1992  double max = 0;
1993  for (unsigned int i=0; i<dim; ++i)
1994  {
1995  double sum = 0;
1996  for (unsigned int j=0; j<dim; ++j)
1997  sum += std::fabs(t[i][j]);
1998 
1999  if (sum > max)
2000  max = sum;
2001  }
2002 
2003  return max;
2004 }
2005 
2007 
2008 DEAL_II_NAMESPACE_CLOSE
2009 
2010 // include deprecated non-member functions operating on Tensor
2011 #include <deal.II/base/tensor_deprecated.h>
2012 
2013 #endif
DEAL_II_CUDA_HOST_DEV Tensor()
Definition: tensor.h:766
numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:116
Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:1782
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1146
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
Definition: tensor.h:1065
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:366
static DEAL_II_CUDA_HOST_DEV real_type abs_square(const number &x)
Definition: numbers.h:356
SymmetricTensor< rank, dim, Number > operator/(const SymmetricTensor< rank, dim, Number > &t, const Number factor)
static const unsigned int n_independent_components
Definition: tensor.h:352
static std::size_t memory_consumption()
Definition: tensor.h:1109
Tensor< rank_, dim, Number > & operator-=(const Tensor< rank_, dim, OtherNumber > &)
Definition: tensor.h:967
Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1924
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:1471
static const unsigned int rank
Definition: tensor.h:345
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1017
#define AssertThrow(cond, exc)
Definition: exceptions.h:369
static real_type abs(const number &x)
Definition: numbers.h:365
bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:945
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:564
void unroll(Vector< OtherNumber > &result) const
Definition: tensor.h:1041
Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1947
Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1892
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:1628
Tensor< rank_, dim, Number > & operator+=(const Tensor< rank_, dim, OtherNumber > &)
Definition: tensor.h:955
Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:1803
Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:1684
bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:918
#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:1595
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: tensor.h:1054
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:558
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:1119
double linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1990
std::size_t size() const
DEAL_II_CUDA_HOST_DEV Tensor< rank_, dim, Number > & operator*=(const OtherNumber factor)
Definition: tensor.h:980
double l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1963
Tensor< rank_, dim, Number > & operator/=(const OtherNumber factor)
Definition: tensor.h:992
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Definition: tensor.h:1078
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:1533
Tensor< 1, dim, Number > cross_product_3d(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:1710
Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1740
static const unsigned int dimension
Definition: tensor.h:340
Tensor< rank_, dim, Number > operator-() const
Definition: tensor.h:1003
Definition: mpi.h:41
Tensor< rank_-1, dim, Number >::tensor_type value_type
Definition: tensor.h:359
DEAL_II_CUDA_HOST_DEV value_type & operator[](const unsigned int i)
Definition: tensor.h:850
static::ExceptionBase & ExcNotImplemented()
DEAL_II_CUDA_HOST_DEV numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1027
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:1654
void clear()
Definition: tensor.h:1099
static::ExceptionBase & ExcInternalError()
Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:1768
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const