Reference documentation for deal.II version Git 1ab60d1 2017-10-16 22:02:15 +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/numbers.h>
22 #include <deal.II/base/table_indices.h>
23 #include <deal.II/base/tensor_accessors.h>
24 #include <deal.II/base/template_constraints.h>
25 #include <deal.II/base/utilities.h>
26 
27 #include <cmath>
28 #include <ostream>
29 #include <vector>
30 
31 
32 DEAL_II_NAMESPACE_OPEN
33 
34 // Forward declarations:
35 
36 template <int dim, typename Number> class Point;
37 template <int rank_, int dim, typename Number = double> class Tensor;
38 template <typename Number> class Vector;
39 template <typename Number> class VectorizedArray;
40 
41 #ifndef DOXYGEN
42 // Overload invalid tensor types of negative rank that come up during
43 // overload resolution of operator* and related contraction variants.
44 template <int dim, typename Number>
45 class Tensor<-2, dim, Number>
46 {
47 };
48 
49 template <int dim, typename Number>
50 class Tensor<-1, dim, Number>
51 {
52 };
53 #endif /* DOXYGEN */
54 
55 
86 template <int dim, typename Number>
87 class Tensor<0,dim,Number>
88 {
89 public:
98  static const unsigned int dimension = dim;
99 
103  static const unsigned int rank = 0;
104 
108  static const unsigned int n_independent_components = 1;
109 
119 
124  typedef Number value_type;
125 
131  typedef Number array_type;
132 
138  DEAL_II_CUDA_HOST_DEV Tensor ();
139 
145  template <typename OtherNumber>
146  Tensor (const Tensor<0,dim,OtherNumber> &initializer);
147 
151  template <typename OtherNumber>
152  Tensor (const OtherNumber &initializer);
153 
163  DEAL_II_CUDA_HOST_DEV operator Number &();
164 
173  DEAL_II_CUDA_HOST_DEV operator const Number &() const;
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 
212  template <typename OtherNumber>
213  DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number> &operator *= (const OtherNumber &factor);
214 
218  template <typename OtherNumber>
219  Tensor<0,dim,Number> &operator /= (const OtherNumber &factor);
220 
225 
238  void clear ();
239 
245  real_type norm () const;
246 
253  DEAL_II_CUDA_HOST_DEV real_type norm_square () const;
254 
259  template <class Archive>
260  void serialize(Archive &ar, const unsigned int version);
261 
266  typedef Number tensor_type;
267 
268 private:
272  Number value;
273 
277  template <typename OtherNumber>
279  unsigned int &start_index) const;
280 
284  template <int, int, typename> friend class Tensor;
285 };
286 
287 
288 
330 template <int rank_, int dim, typename Number>
331 class Tensor
332 {
333 public:
342  static const unsigned int dimension = dim;
343 
347  static const unsigned int rank = rank_;
348 
353  static const unsigned int
355 
361  typedef typename Tensor<rank_-1,dim,Number>::tensor_type value_type;
362 
367  typedef typename Tensor<rank_-1,dim,Number>::array_type
368  array_type[(dim != 0) ? dim : 1];
369  // ... avoid a compiler warning in case of dim == 0 and ensure that the
370  // array always has positive size.
371 
377  DEAL_II_CUDA_HOST_DEV Tensor ();
378 
382  explicit Tensor (const array_type &initializer);
383 
389  template <typename OtherNumber>
390  Tensor (const Tensor<rank_,dim,OtherNumber> &initializer);
391 
395  template <typename OtherNumber>
396  Tensor (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer);
397 
401  template <typename OtherNumber>
402  operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const;
403 
409  DEAL_II_CUDA_HOST_DEV value_type &operator [] (const unsigned int i);
410 
416  DEAL_II_CUDA_HOST_DEV const value_type &operator[](const unsigned int i) const;
417 
421  const Number &operator [] (const TableIndices<rank_> &indices) const;
422 
426  Number &operator [] (const TableIndices<rank_> &indices);
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 
474  template <typename OtherNumber>
475  DEAL_II_CUDA_HOST_DEV Tensor<rank_,dim,Number> &operator *= (const OtherNumber &factor);
476 
480  template <typename OtherNumber>
481  Tensor<rank_,dim,Number> &operator /= (const OtherNumber &factor);
482 
487 
500  void clear ();
501 
509 
516  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type norm_square() const;
517 
525  template <typename OtherNumber>
526  void unroll (Vector<OtherNumber> &result) const;
527 
532  static
533  unsigned int
535 
540  static
542 
547  static std::size_t memory_consumption ();
548 
553  template <class Archive>
554  void serialize(Archive &ar, const unsigned int version);
555 
561 
562 private:
566  Tensor<rank_-1, dim, Number> values[(dim != 0) ? dim : 1];
567  // ... avoid a compiler warning in case of dim == 0 and ensure that the
568  // array always has positive size.
569 
573  template <typename OtherNumber>
575  unsigned int &start_index) const;
576 
580  template <int, int, typename> friend class Tensor;
581 
586  friend class Point<dim,Number>;
587 };
588 
589 
590 namespace internal
591 {
602  template <int rank, int dim, typename T>
603  struct NumberType<Tensor<rank,dim,T> >
604  {
605  static const Tensor<rank,dim,T> &value (const Tensor<rank,dim,T> &t)
606  {
607  return t;
608  }
609 
610  static Tensor<rank,dim,T> value (const T &t)
611  {
612  Tensor<rank,dim,T> tmp;
613  tmp=t;
614  return tmp;
615  }
616  };
617 
618  template <int rank, int dim, typename T>
619  struct NumberType<Tensor<rank,dim,VectorizedArray<T> > >
620  {
621  static const Tensor<rank,dim,VectorizedArray<T> > &value (const Tensor<rank,dim,VectorizedArray<T> > &t)
622  {
623  return t;
624  }
625 
626  static Tensor<rank,dim,VectorizedArray<T> > value (const T &t)
627  {
630  return tmp;
631  }
632 
634  {
636  tmp=t;
637  return tmp;
638  }
639  };
640 }
641 
642 
643 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
644 
645 
646 template <int dim,typename Number>
647 inline
648 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::Tensor ()
649 // Some auto-differentiable numbers need explicit
650 // zero initialization.
651  : value(internal::NumberType<Number>::value(0.0))
652 {
653 }
654 
655 
656 template <int dim, typename Number>
657 template <typename OtherNumber>
658 inline
659 Tensor<0,dim,Number>::Tensor (const OtherNumber &initializer)
660 {
661  value = internal::NumberType<Number>::value(initializer);
662 }
663 
664 
665 template <int dim, typename Number>
666 template <typename OtherNumber>
667 inline
669 {
670  value = p.value;
671 }
672 
673 
674 template <int dim, typename Number>
675 inline
676 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::operator Number &()
677 {
678  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
679  return value;
680 }
681 
682 
683 template <int dim, typename Number>
684 inline
685 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::operator const Number &() const
686 {
687  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
688  return value;
689 }
690 
691 
692 template <int dim, typename Number>
693 template <typename OtherNumber>
694 inline
696 {
697  value = p.value;
698  return *this;
699 }
700 
701 
702 template <int dim, typename Number>
703 template <typename OtherNumber>
704 inline
706 {
707  return (value == p.value);
708 }
709 
710 
711 template <int dim, typename Number>
712 template <typename OtherNumber>
713 inline
715 {
716  return !((*this) == p);
717 }
718 
719 
720 template <int dim, typename Number>
721 template <typename OtherNumber>
722 inline
724 {
725  value += p.value;
726  return *this;
727 }
728 
729 
730 template <int dim, typename Number>
731 template <typename OtherNumber>
732 inline
734 {
735  value -= p.value;
736  return *this;
737 }
738 
739 
740 template <int dim, typename Number>
741 template <typename OtherNumber>
742 inline
743 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator *= (const OtherNumber &s)
744 {
745  value *= s;
746  return *this;
747 }
748 
749 
750 template <int dim, typename Number>
751 template <typename OtherNumber>
752 inline
754 {
755  value /= s;
756  return *this;
757 }
758 
759 
760 template <int dim, typename Number>
761 inline
763 {
764  return -value;
765 }
766 
767 
768 template <int dim, typename Number>
769 inline
772 {
773  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
774  return numbers::NumberTraits<Number>::abs (value);
775 }
776 
777 
778 template <int dim, typename Number>
779 inline
781 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::norm_square () const
782 {
783  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
785 }
786 
787 
788 template <int dim, typename Number>
789 template <typename OtherNumber>
790 inline
791 void
793  unsigned int &index) const
794 {
795  Assert(dim != 0, ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
796  result[index] = value;
797  ++index;
798 }
799 
800 
801 template <int dim, typename Number>
802 inline
804 {
805  // Some auto-differentiable numbers need explicit
806  // zero initialization.
808 }
809 
810 
811 template <int dim, typename Number>
812 template <class Archive>
813 inline
814 void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int)
815 {
816  ar &value;
817 }
818 
819 
820 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
821 
822 
823 template <int rank_, int dim, typename Number>
824 inline
825 DEAL_II_CUDA_HOST_DEV Tensor<rank_,dim,Number>::Tensor ()
826 {
827  // All members of the c-style array values are already default initialized
828  // and thus all values are already set to zero recursively.
829 }
830 
831 
832 template <int rank_, int dim, typename Number>
833 inline
834 Tensor<rank_,dim,Number>::Tensor (const array_type &initializer)
835 {
836  for (unsigned int i=0; i<dim; ++i)
837  values[i] = Tensor<rank_-1, dim, Number>(initializer[i]);
838 }
839 
840 
841 template <int rank_, int dim, typename Number>
842 template <typename OtherNumber>
843 inline
845 {
846  for (unsigned int i=0; i!=dim; ++i)
847  values[i] = Tensor<rank_-1,dim,Number>(initializer[i]);
848 }
849 
850 
851 template <int rank_, int dim, typename Number>
852 template <typename OtherNumber>
853 inline
855 (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer)
856 {
857  for (unsigned int i=0; i<dim; ++i)
858  values[i] = initializer[i];
859 }
860 
861 
862 template <int rank_, int dim, typename Number>
863 template <typename OtherNumber>
864 inline
866 operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const
867 {
868  return Tensor<1,dim,Tensor<rank_-1,dim,Number> > (values);
869 }
870 
871 
872 
873 namespace internal
874 {
875  namespace TensorSubscriptor
876  {
877  template <typename ArrayElementType, int dim>
878  inline DEAL_II_ALWAYS_INLINE
879  DEAL_II_CUDA_HOST_DEV
880  ArrayElementType &
881  subscript (ArrayElementType *values,
882  const unsigned int i,
883  std::integral_constant<int, dim>)
884  {
885  Assert (i<dim, ExcIndexRange(i, 0, dim));
886  return values[i];
887  }
888 
889 
890  template <typename ArrayElementType>
891  DEAL_II_CUDA_HOST_DEV
892  ArrayElementType &
893  subscript (ArrayElementType *,
894  const unsigned int,
895  std::integral_constant<int, 0>)
896  {
897  Assert(false, ExcMessage("Cannot access elements of an object of type Tensor<rank,0,Number>."));
898  static ArrayElementType t;
899  return t;
900  }
901  }
902 }
903 
904 
905 template <int rank_, int dim, typename Number>
906 inline DEAL_II_ALWAYS_INLINE
907 DEAL_II_CUDA_HOST_DEV
910 {
911  return ::internal::TensorSubscriptor::subscript(values, i, std::integral_constant<int, dim>());
912 }
913 
914 
915 template <int rank_, int dim, typename Number>
916 inline DEAL_II_ALWAYS_INLINE
917 DEAL_II_CUDA_HOST_DEV
919 Tensor<rank_,dim,Number>::operator[] (const unsigned int i) const
920 {
921  return ::internal::TensorSubscriptor::subscript(values, i, std::integral_constant<int, dim>());
922 }
923 
924 
925 template <int rank_, int dim, typename Number>
926 inline
927 const Number &
929 {
930  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
931 
932  return TensorAccessors::extract<rank_>(*this, indices);
933 }
934 
935 
936 template <int rank_, int dim, typename Number>
937 inline
938 Number &
940 {
941  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
942 
943  return TensorAccessors::extract<rank_>(*this, indices);
944 }
945 
946 
947 template <int rank_, int dim, typename Number>
948 template <typename OtherNumber>
949 inline
952 {
953  if (dim > 0)
954  std::copy (&t.values[0], &t.values[0]+dim, &values[0]);
955  return *this;
956 }
957 
958 
959 template <int rank_, int dim, typename Number>
960 inline
963 {
965  ExcMessage ("Only assignment with zero is allowed"));
966  (void) d;
967 
968  for (unsigned int i=0; i<dim; ++i)
969  values[i] = internal::NumberType<Number>::value(0.0);
970  return *this;
971 }
972 
973 
974 template <int rank_, int dim, typename Number>
975 template <typename OtherNumber>
976 inline
977 bool
979 {
980  for (unsigned int i=0; i<dim; ++i)
981  if (values[i] != p.values[i])
982  return false;
983  return true;
984 }
985 
986 
987 // At some places in the library, we have Point<0> for formal reasons
988 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
989 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
990 // in the above function that the loop end check always fails, we
991 // implement this function here
992 template <>
993 template <>
994 inline
996 {
997  return true;
998 }
999 
1000 
1001 template <int rank_, int dim, typename Number>
1002 template <typename OtherNumber>
1003 inline
1004 bool
1006 {
1007  return !((*this) == p);
1008 }
1009 
1010 
1011 template <int rank_, int dim, typename Number>
1012 template <typename OtherNumber>
1013 inline
1016 {
1017  for (unsigned int i=0; i<dim; ++i)
1018  values[i] += p.values[i];
1019  return *this;
1020 }
1021 
1022 
1023 template <int rank_, int dim, typename Number>
1024 template <typename OtherNumber>
1025 inline
1028 {
1029  for (unsigned int i=0; i<dim; ++i)
1030  values[i] -= p.values[i];
1031  return *this;
1032 }
1033 
1034 
1035 template <int rank_, int dim, typename Number>
1036 template <typename OtherNumber>
1037 inline
1038 DEAL_II_CUDA_HOST_DEV
1041 {
1042  for (unsigned int i=0; i<dim; ++i)
1043  values[i] *= s;
1044  return *this;
1045 }
1046 
1047 
1048 template <int rank_, int dim, typename Number>
1049 template <typename OtherNumber>
1050 inline
1053 {
1054  for (unsigned int i=0; i<dim; ++i)
1055  values[i] /= s;
1056  return *this;
1057 }
1058 
1059 
1060 template <int rank_, int dim, typename Number>
1061 inline
1064 {
1066 
1067  for (unsigned int i=0; i<dim; ++i)
1068  tmp.values[i] = -values[i];
1069 
1070  return tmp;
1071 }
1072 
1073 
1074 template <int rank_, int dim, typename Number>
1075 inline
1078 {
1079  return std::sqrt (norm_square());
1080 }
1081 
1082 
1083 template <int rank_, int dim, typename Number>
1084 inline
1085 DEAL_II_CUDA_HOST_DEV
1088 {
1091  for (unsigned int i=0; i<dim; ++i)
1092  s += values[i].norm_square();
1093 
1094  return s;
1095 }
1096 
1097 
1098 template <int rank_, int dim, typename Number>
1099 template <typename OtherNumber>
1100 inline
1101 void
1103 {
1104  AssertDimension (result.size(),(Utilities::fixed_power<rank_, unsigned int>(dim)));
1105 
1106  unsigned int index = 0;
1107  unroll_recursion (result, index);
1108 }
1109 
1110 
1111 template <int rank_, int dim, typename Number>
1112 template <typename OtherNumber>
1113 inline
1114 void
1116  unsigned int &index) const
1117 {
1118  for (unsigned int i=0; i<dim; ++i)
1119  values[i].unroll_recursion(result, index);
1120 }
1121 
1122 
1123 template <int rank_, int dim, typename Number>
1124 inline
1125 unsigned int
1127 {
1128  unsigned int index = 0;
1129  for (int r = 0; r < rank_; ++r)
1130  index = index * dim + indices[r];
1131 
1132  return index;
1133 }
1134 
1135 
1136 template <int rank_, int dim, typename Number>
1137 inline
1140 {
1141  Assert (i < n_independent_components,
1142  ExcIndexRange (i, 0, n_independent_components));
1143 
1144  TableIndices<rank_> indices;
1145 
1146  unsigned int remainder = i;
1147  for (int r=rank_-1; r>=0; --r)
1148  {
1149  indices[r] = (remainder % dim);
1150  remainder /= dim;
1151  }
1152  Assert (remainder == 0, ExcInternalError());
1153 
1154  return indices;
1155 }
1156 
1157 
1158 template <int rank_, int dim, typename Number>
1159 inline
1161 {
1162  for (unsigned int i=0; i<dim; ++i)
1163  values[i] = internal::NumberType<Number>::value(0.0);
1164 }
1165 
1166 
1167 template <int rank_, int dim, typename Number>
1168 inline
1169 std::size_t
1171 {
1172  return sizeof(Tensor<rank_,dim,Number>);
1173 }
1174 
1175 
1176 template <int rank_, int dim, typename Number>
1177 template <class Archive>
1178 inline
1179 void
1180 Tensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
1181 {
1182  ar &values;
1183 }
1184 
1185 
1186 /* ----------------- Non-member functions operating on tensors. ------------ */
1187 
1192 
1200 template <int rank_, int dim, typename Number>
1201 inline
1202 std::ostream &operator << (std::ostream &out, const Tensor<rank_,dim,Number> &p)
1203 {
1204  for (unsigned int i = 0; i < dim; ++i)
1205  {
1206  out << p[i];
1207  if (i != dim - 1)
1208  out << ' ';
1209  }
1210 
1211  return out;
1212 }
1213 
1214 
1221 template <int dim, typename Number>
1222 inline
1223 std::ostream &operator << (std::ostream &out, const Tensor<0,dim,Number> &p)
1224 {
1225  out << static_cast<const Number &>(p);
1226  return out;
1227 }
1228 
1229 
1231 
1235 
1236 
1245 template <int dim, typename Number, typename Other>
1246 inline
1247 typename ProductType<Other, Number>::type
1248 operator * (const Other &object,
1249  const Tensor<0,dim,Number> &t)
1250 {
1251  return object * static_cast<const Number &>(t);
1252 }
1253 
1254 
1255 
1264 template <int dim, typename Number, typename Other>
1265 inline
1266 typename ProductType<Number, Other>::type
1268  const Other &object)
1269 {
1270  return static_cast<const Number &>(t) * object;
1271 }
1272 
1273 
1283 template <int dim, typename Number, typename OtherNumber>
1284 inline
1285 typename ProductType<Number, OtherNumber>::type
1287  const Tensor<0, dim, OtherNumber> &src2)
1288 {
1289  return static_cast<const Number &>(src1) *
1290  static_cast<const OtherNumber &>(src2);
1291 }
1292 
1293 
1299 template <int dim, typename Number, typename OtherNumber>
1300 inline
1303  const OtherNumber &factor)
1304 {
1305  return static_cast<const Number &>(t) / factor;
1306 }
1307 
1308 
1314 template <int dim, typename Number, typename OtherNumber>
1315 inline
1318  const Tensor<0,dim,OtherNumber> &q)
1319 {
1320  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1321 }
1322 
1323 
1329 template <int dim, typename Number, typename OtherNumber>
1330 inline
1333  const Tensor<0,dim,OtherNumber> &q)
1334 {
1335  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1336 }
1337 
1338 
1349 template <int rank, int dim,
1350  typename Number,
1351  typename OtherNumber>
1352 inline
1355  const OtherNumber &factor)
1356 {
1357  // recurse over the base objects
1359  for (unsigned int d=0; d<dim; ++d)
1360  tt[d] = t[d] * factor;
1361  return tt;
1362 }
1363 
1364 
1375 template <int rank, int dim,
1376  typename Number,
1377  typename OtherNumber>
1378 inline
1380 operator * (const Number &factor,
1382 {
1383  // simply forward to the operator above
1384  return t * factor;
1385 }
1386 
1387 
1395 template <int rank, int dim,
1396  typename Number,
1397  typename OtherNumber>
1398 inline
1401  const OtherNumber &factor)
1402 {
1403  // recurse over the base objects
1405  for (unsigned int d=0; d<dim; ++d)
1406  tt[d] = t[d] / factor;
1407  return tt;
1408 }
1409 
1410 
1418 template <int rank, int dim, typename Number, typename OtherNumber>
1419 inline
1423 {
1425 
1426  for (unsigned int i=0; i<dim; ++i)
1427  tmp[i] += q[i];
1428 
1429  return tmp;
1430 }
1431 
1432 
1440 template <int rank, int dim, typename Number, typename OtherNumber>
1441 inline
1445 {
1447 
1448  for (unsigned int i=0; i<dim; ++i)
1449  tmp[i] -= q[i];
1450 
1451  return tmp;
1452 }
1453 
1454 
1456 
1460 
1461 
1485 template <int rank_1, int rank_2, int dim,
1486  typename Number, typename OtherNumber>
1487 inline DEAL_II_ALWAYS_INLINE
1488 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1491 {
1492  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
1493 
1494  TensorAccessors::internal::ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1495  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
1496  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
1497 
1498  return result;
1499 }
1500 
1501 
1531 template <int index_1, int index_2,
1532  int rank_1, int rank_2, int dim,
1533  typename Number, typename OtherNumber>
1534 inline
1535 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1538 {
1539  Assert(0 <= index_1 && index_1 < rank_1,
1540  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1541  Assert(0 <= index_2 && index_2 < rank_2,
1542  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1543 
1544  using namespace TensorAccessors;
1545  using namespace TensorAccessors::internal;
1546 
1547  // Reorder index_1 to the end of src1:
1548  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1549  reord_01 = reordered_index_view<index_1, rank_1>(src1);
1550 
1551  // Reorder index_2 to the end of src2:
1552  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1553  reord_02 = reordered_index_view<index_2, rank_2>(src2);
1554 
1555  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1556  result;
1557  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
1558  return result;
1559 }
1560 
1561 
1593 template <int index_1, int index_2, int index_3, int index_4,
1594  int rank_1, int rank_2, int dim,
1595  typename Number, typename OtherNumber>
1596 inline
1597 typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1600 {
1601  Assert(0 <= index_1 && index_1 < rank_1,
1602  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1603  Assert(0 <= index_3 && index_3 < rank_1,
1604  ExcMessage("The specified index_3 must lie within the range [0,rank_1)"));
1605  Assert(index_1 != index_3,
1606  ExcMessage("index_1 and index_3 must not be the same"));
1607  Assert(0 <= index_2 && index_2 < rank_2,
1608  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1609  Assert(0 <= index_4 && index_4 < rank_2,
1610  ExcMessage("The specified index_4 must lie within the range [0,rank_2)"));
1611  Assert(index_2 != index_4,
1612  ExcMessage("index_2 and index_4 must not be the same"));
1613 
1614  using namespace TensorAccessors;
1615  using namespace TensorAccessors::internal;
1616 
1617  // Reorder index_1 to the end of src1:
1618  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1619  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
1620 
1621  // Reorder index_2 to the end of src2:
1622  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1623  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
1624 
1625  // Now, reorder index_3 to the end of src1. We have to make sure to
1626  // preserve the orginial ordering: index_1 has been removed. If
1627  // index_3 > index_1, we have to use (index_3 - 1) instead:
1628  ReorderedIndexView<(index_3 < index_1 ? index_3 : index_3 - 1), rank_1, ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> > >
1629  reord_3 = TensorAccessors::reordered_index_view<index_3 < index_1 ? index_3 : index_3 - 1, rank_1>(reord_1);
1630 
1631  // Now, reorder index_4 to the end of src2. We have to make sure to
1632  // preserve the orginial ordering: index_2 has been removed. If
1633  // index_4 > index_2, we have to use (index_4 - 1) instead:
1634  ReorderedIndexView<(index_4 < index_2 ? index_4 : index_4 - 1), rank_2, ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> > >
1635  reord_4 = TensorAccessors::reordered_index_view<index_4 < index_2 ? index_4 : index_4 - 1, rank_2>(reord_2);
1636 
1637  typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1638  result;
1639  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
1640  return result;
1641 }
1642 
1643 
1657 template <int rank, int dim, typename Number, typename OtherNumber>
1658 inline
1659 typename ProductType<Number, OtherNumber>::type
1661  const Tensor<rank, dim, OtherNumber> &right)
1662 {
1663  typename ProductType<Number, OtherNumber>::type result;
1664  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
1665  return result;
1666 }
1667 
1668 
1687 template <template <int, int, typename> class TensorT1,
1688  template <int, int, typename> class TensorT2,
1689  template <int, int, typename> class TensorT3,
1690  int rank_1, int rank_2, int dim,
1691  typename T1, typename T2, typename T3>
1693 contract3 (const TensorT1<rank_1, dim, T1> &left,
1694  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
1695  const TensorT3<rank_2, dim, T3> &right)
1696 {
1698  return_type;
1699  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(
1700  left, middle, right);
1701 }
1702 
1703 
1715 template <int rank_1, int rank_2, int dim,
1716  typename Number, typename OtherNumber>
1717 inline
1721 {
1723  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
1724  return result;
1725 }
1726 
1727 
1729 
1733 
1734 
1746 template <int dim, typename Number>
1747 inline
1750 {
1751  Assert (dim==2, ExcInternalError());
1752 
1753  Tensor<1, dim, Number> result;
1754 
1755  result[0] = src[1];
1756  result[1] = -src[0];
1757 
1758  return result;
1759 }
1760 
1761 
1772 template <int dim, typename Number>
1773 inline
1776  const Tensor<1,dim,Number> &src2)
1777 {
1778  Assert (dim==3, ExcInternalError());
1779 
1780  Tensor<1, dim, Number> result;
1781 
1782  result[0] = src1[1]*src2[2] - src1[2]*src2[1];
1783  result[1] = src1[2]*src2[0] - src1[0]*src2[2];
1784  result[2] = src1[0]*src2[1] - src1[1]*src2[0];
1785 
1786  return result;
1787 }
1788 
1789 
1791 
1795 
1796 
1803 template <int dim, typename Number>
1804 inline
1806 {
1807  // Compute the determinant using the Laplace expansion of the
1808  // determinant. We expand along the last row.
1809  Number det = internal::NumberType<Number>::value(0.0);
1810 
1811  for (unsigned int k=0; k<dim; ++k)
1812  {
1813  Tensor<2,dim-1,Number> minor;
1814  for (unsigned int i=0; i<dim-1; ++i)
1815  for (unsigned int j=0; j<dim-1; ++j)
1816  minor[i][j] = t[i][j<k ? j : j+1];
1817 
1818  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
1819 
1820  det += t[dim-1][k] * cofactor;
1821  }
1822 
1823  return ((dim % 2 == 0) ? 1. : -1.) * det;
1824 }
1825 
1831 template <typename Number>
1832 inline
1834 {
1835  return t[0][0];
1836 }
1837 
1838 
1846 template <int dim, typename Number>
1847 Number trace (const Tensor<2,dim,Number> &d)
1848 {
1849  Number t=d[0][0];
1850  for (unsigned int i=1; i<dim; ++i)
1851  t += d[i][i];
1852  return t;
1853 }
1854 
1855 
1865 template <int dim, typename Number>
1866 inline
1869 {
1870  Number return_tensor [dim][dim];
1871 
1872  // if desired, take over the
1873  // inversion of a 4x4 tensor
1874  // from the FullMatrix
1875  AssertThrow (false, ExcNotImplemented());
1876 
1877  return Tensor<2,dim,Number>(return_tensor);
1878 }
1879 
1880 
1881 #ifndef DOXYGEN
1882 
1883 template <typename Number>
1884 inline
1886 invert (const Tensor<2,1,Number> &t)
1887 {
1888  Number return_tensor [1][1];
1889 
1890  return_tensor[0][0] = 1.0/t[0][0];
1891 
1892  return Tensor<2,1,Number>(return_tensor);
1893 }
1894 
1895 
1896 template <typename Number>
1897 inline
1899 invert (const Tensor<2,2,Number> &t)
1900 {
1901  Tensor<2,2,Number> return_tensor;
1902 
1903  // this is Maple output,
1904  // thus a bit unstructured
1905  const Number inv_det_t = 1.0/(t[0][0]*t[1][1]-t[1][0]*t[0][1]);
1906  return_tensor[0][0] = t[1][1];
1907  return_tensor[0][1] = -t[0][1];
1908  return_tensor[1][0] = -t[1][0];
1909  return_tensor[1][1] = t[0][0];
1910  return_tensor *= inv_det_t;
1911 
1912  return return_tensor;
1913 }
1914 
1915 
1916 template <typename Number>
1917 inline
1919 invert (const Tensor<2,3,Number> &t)
1920 {
1921  Tensor<2,3,Number> return_tensor;
1922 
1923  const Number t4 = t[0][0]*t[1][1],
1924  t6 = t[0][0]*t[1][2],
1925  t8 = t[0][1]*t[1][0],
1926  t00 = t[0][2]*t[1][0],
1927  t01 = t[0][1]*t[2][0],
1928  t04 = t[0][2]*t[2][0],
1929  inv_det_t = 1.0/(t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
1930  t00*t[2][1]+t01*t[1][2]-t04*t[1][1]);
1931  return_tensor[0][0] = t[1][1]*t[2][2]-t[1][2]*t[2][1];
1932  return_tensor[0][1] = t[0][2]*t[2][1]-t[0][1]*t[2][2];
1933  return_tensor[0][2] = t[0][1]*t[1][2]-t[0][2]*t[1][1];
1934  return_tensor[1][0] = t[1][2]*t[2][0]-t[1][0]*t[2][2];
1935  return_tensor[1][1] = t[0][0]*t[2][2]-t04;
1936  return_tensor[1][2] = t00-t6;
1937  return_tensor[2][0] = t[1][0]*t[2][1]-t[1][1]*t[2][0];
1938  return_tensor[2][1] = t01-t[0][0]*t[2][1];
1939  return_tensor[2][2] = t4-t8;
1940  return_tensor *= inv_det_t;
1941 
1942  return return_tensor;
1943 }
1944 
1945 #endif /* DOXYGEN */
1946 
1947 
1954 template <int dim, typename Number>
1955 inline
1958 {
1960  for (unsigned int i=0; i<dim; ++i)
1961  {
1962  tt[i][i] = t[i][i];
1963  for (unsigned int j=i+1; j<dim; ++j)
1964  {
1965  tt[i][j] = t[j][i];
1966  tt[j][i] = t[i][j];
1967  };
1968  }
1969  return tt;
1970 }
1971 
1972 
1986 template <int dim, typename Number>
1987 inline
1990 {
1991  return determinant(t)*invert(t);
1992 }
1993 
1994 
2009 template <int dim, typename Number>
2010 inline
2013 {
2014  return transpose(adjugate(t));
2015 }
2016 
2017 
2025 template <int dim, typename Number>
2026 inline
2027 Number
2029 {
2030  Number max = internal::NumberType<Number>::value(0.0);
2031  for (unsigned int j=0; j<dim; ++j)
2032  {
2033  Number sum = internal::NumberType<Number>::value(0.0);
2034  for (unsigned int i=0; i<dim; ++i)
2035  sum += std::fabs(t[i][j]);
2036 
2037  if (sum > max)
2038  max = sum;
2039  }
2040 
2041  return max;
2042 }
2043 
2044 
2052 template <int dim, typename Number>
2053 inline
2054 Number
2056 {
2057  Number max = internal::NumberType<Number>::value(0.0);
2058  for (unsigned int i=0; i<dim; ++i)
2059  {
2060  Number sum = internal::NumberType<Number>::value(0.0);
2061  for (unsigned int j=0; j<dim; ++j)
2062  sum += std::fabs(t[i][j]);
2063 
2064  if (sum > max)
2065  max = sum;
2066  }
2067 
2068  return max;
2069 }
2070 
2072 
2073 DEAL_II_NAMESPACE_CLOSE
2074 
2075 // include deprecated non-member functions operating on Tensor
2076 #include <deal.II/base/tensor_deprecated.h>
2077 
2078 #endif
DEAL_II_CUDA_HOST_DEV Tensor()
Definition: tensor.h:825
numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:118
Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:1847
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1183
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
Definition: tensor.h:1126
Tensor< rank_-1, dim, Number >::array_type array_type[(dim!=0)?dim:1]
Definition: tensor.h:368
static DEAL_II_CUDA_HOST_DEV real_type abs_square(const number &x)
Definition: numbers.h:342
static const unsigned int n_independent_components
Definition: tensor.h:354
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2028
static std::size_t memory_consumption()
Definition: tensor.h:1170
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2055
Tensor< rank_, dim, Number > & operator-=(const Tensor< rank_, dim, OtherNumber > &)
Definition: tensor.h:1027
Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1989
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:1536
static const unsigned int rank
Definition: tensor.h:347
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1077
#define AssertThrow(cond, exc)
Definition: exceptions.h:398
static real_type abs(const number &x)
Definition: numbers.h:351
bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1005
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:566
void unroll(Vector< OtherNumber > &result) const
Definition: tensor.h:1102
Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2012
Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1957
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:1693
Tensor< rank_, dim, Number > & operator+=(const Tensor< rank_, dim, OtherNumber > &)
Definition: tensor.h:1015
Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:1868
Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:1749
bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:978
#define Assert(cond, exc)
Definition: exceptions.h:337
Tensor< rank_, dim, Number > & operator/=(const OtherNumber &factor)
Definition: tensor.h:1052
ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:1660
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: tensor.h:1115
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:560
DEAL_II_ALWAYS_INLINE internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
void serialize(Archive &ar, const unsigned int version)
Definition: tensor.h:1180
SymmetricTensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const SymmetricTensor< rank, dim, Number > &t, const OtherNumber &factor)
std::size_t size() const
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Definition: tensor.h:1139
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:1598
Tensor< 1, dim, Number > cross_product_3d(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:1775
Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1805
static const unsigned int dimension
Definition: tensor.h:342
Tensor< rank_, dim, Number > operator-() const
Definition: tensor.h:1063
Definition: mpi.h:52
Tensor< rank_-1, dim, Number >::tensor_type value_type
Definition: tensor.h:361
DEAL_II_CUDA_HOST_DEV value_type & operator[](const unsigned int i)
Definition: tensor.h:909
static::ExceptionBase & ExcNotImplemented()
DEAL_II_CUDA_HOST_DEV numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1087
SymmetricTensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank, dim, Number > &left, const SymmetricTensor< rank, dim, OtherNumber > &right)
SymmetricTensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank, dim, Number > &left, const SymmetricTensor< rank, dim, OtherNumber > &right)
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:1719
void clear()
Definition: tensor.h:1160
static::ExceptionBase & ExcInternalError()
Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:1833
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const
DEAL_II_CUDA_HOST_DEV Tensor< rank_, dim, Number > & operator*=(const OtherNumber &factor)
Definition: tensor.h:1040