Reference documentation for deal.II version Git e8a02dd 2017-08-23 11:09:52 +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 template <typename Number> class VectorizedArray;
39 
40 #ifndef DOXYGEN
41 // Overload invalid tensor types of negative rank that come up during
42 // overload resolution of operator* and related contraction variants.
43 template <int dim, typename Number>
44 class Tensor<-2, dim, Number>
45 {
46 };
47 
48 template <int dim, typename Number>
49 class Tensor<-1, dim, Number>
50 {
51 };
52 #endif /* DOXYGEN */
53 
54 
85 template <int dim, typename Number>
86 class Tensor<0,dim,Number>
87 {
88 public:
97  static const unsigned int dimension = dim;
98 
102  static const unsigned int rank = 0;
103 
107  static const unsigned int n_independent_components = 1;
108 
118 
123  typedef Number value_type;
124 
130  typedef Number array_type;
131 
137  DEAL_II_CUDA_HOST_DEV Tensor ();
138 
144  template <typename OtherNumber>
145  Tensor (const Tensor<0,dim,OtherNumber> &initializer);
146 
150  template <typename OtherNumber>
151  Tensor (const OtherNumber &initializer);
152 
162  DEAL_II_CUDA_HOST_DEV operator Number &();
163 
172  DEAL_II_CUDA_HOST_DEV operator const Number &() const;
173 
179  template <typename OtherNumber>
181 
185  template <typename OtherNumber>
186  bool operator == (const Tensor<0,dim,OtherNumber> &rhs) const;
187 
191  template <typename OtherNumber>
192  bool operator != (const Tensor<0,dim,OtherNumber> &rhs) const;
193 
197  template <typename OtherNumber>
199 
203  template <typename OtherNumber>
205 
211  template <typename OtherNumber>
212  DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number> &operator *= (const OtherNumber &factor);
213 
217  template <typename OtherNumber>
218  Tensor<0,dim,Number> &operator /= (const OtherNumber &factor);
219 
224 
237  void clear ();
238 
244  real_type norm () const;
245 
252  DEAL_II_CUDA_HOST_DEV real_type norm_square () const;
253 
258  template <class Archive>
259  void serialize(Archive &ar, const unsigned int version);
260 
265  typedef Number tensor_type;
266 
267 private:
271  Number value;
272 
276  template <typename OtherNumber>
278  unsigned int &start_index) const;
279 
283  template <int, int, typename> friend class Tensor;
284 };
285 
286 
287 
329 template <int rank_, int dim, typename Number>
330 class Tensor
331 {
332 public:
341  static const unsigned int dimension = dim;
342 
346  static const unsigned int rank = rank_;
347 
352  static const unsigned int
354 
360  typedef typename Tensor<rank_-1,dim,Number>::tensor_type value_type;
361 
366  typedef typename Tensor<rank_-1,dim,Number>::array_type
367  array_type[(dim != 0) ? dim : 1];
368  // ... avoid a compiler warning in case of dim == 0 and ensure that the
369  // array always has positive size.
370 
376  DEAL_II_CUDA_HOST_DEV Tensor ();
377 
381  explicit Tensor (const array_type &initializer);
382 
388  template <typename OtherNumber>
389  Tensor (const Tensor<rank_,dim,OtherNumber> &initializer);
390 
394  template <typename OtherNumber>
395  Tensor (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer);
396 
400  template <typename OtherNumber>
401  operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const;
402 
408  DEAL_II_CUDA_HOST_DEV value_type &operator [] (const unsigned int i);
409 
415  DEAL_II_CUDA_HOST_DEV const value_type &operator[](const unsigned int i) const;
416 
420  const Number &operator [] (const TableIndices<rank_> &indices) const;
421 
425  Number &operator [] (const TableIndices<rank_> &indices);
426 
432  template <typename OtherNumber>
434 
441  Tensor<rank_,dim,Number> &operator = (const Number &d);
442 
446  template <typename OtherNumber>
447  bool operator == (const Tensor<rank_,dim,OtherNumber> &) const;
448 
452  template <typename OtherNumber>
453  bool operator != (const Tensor<rank_,dim,OtherNumber> &) const;
454 
458  template <typename OtherNumber>
460 
464  template <typename OtherNumber>
466 
473  template <typename OtherNumber>
474  DEAL_II_CUDA_HOST_DEV Tensor<rank_,dim,Number> &operator *= (const OtherNumber &factor);
475 
479  template <typename OtherNumber>
480  Tensor<rank_,dim,Number> &operator /= (const OtherNumber &factor);
481 
486 
499  void clear ();
500 
508 
515  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type norm_square() const;
516 
524  template <typename OtherNumber>
525  void unroll (Vector<OtherNumber> &result) const;
526 
531  static
532  unsigned int
534 
539  static
541 
546  static std::size_t memory_consumption ();
547 
552  template <class Archive>
553  void serialize(Archive &ar, const unsigned int version);
554 
560 
561 private:
565  Tensor<rank_-1, dim, Number> values[(dim != 0) ? dim : 1];
566  // ... avoid a compiler warning in case of dim == 0 and ensure that the
567  // array always has positive size.
568 
572  template <typename OtherNumber>
574  unsigned int &start_index) const;
575 
579  template <int, int, typename> friend class Tensor;
580 
585  friend class Point<dim,Number>;
586 };
587 
588 
589 namespace internal
590 {
601  template <int rank, int dim, typename T>
602  struct NumberType<Tensor<rank,dim,T> >
603  {
604  static Tensor<rank,dim,T> value (const T &t)
605  {
606  Tensor<rank,dim,T> tmp;
607  tmp=t;
608  return tmp;
609  }
610  };
611 
612  template <int rank, int dim, typename T>
613  struct NumberType<Tensor<rank,dim,VectorizedArray<T> > >
614  {
615  static Tensor<rank,dim,VectorizedArray<T> > value (const T &t)
616  {
619  return tmp;
620  }
621  };
622 }
623 
624 
625 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
626 
627 
628 template <int dim,typename Number>
629 inline
630 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::Tensor ()
631 // Some auto-differentiable numbers need explicit
632 // zero initialization.
633  : value(internal::NumberType<Number>::value(0.0))
634 {
635 }
636 
637 
638 template <int dim, typename Number>
639 template <typename OtherNumber>
640 inline
641 Tensor<0,dim,Number>::Tensor (const OtherNumber &initializer)
642 {
643  value = initializer;
644 }
645 
646 
647 template <int dim, typename Number>
648 template <typename OtherNumber>
649 inline
651 {
652  value = p.value;
653 }
654 
655 
656 template <int dim, typename Number>
657 inline
658 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::operator Number &()
659 {
660  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
661  return value;
662 }
663 
664 
665 template <int dim, typename Number>
666 inline
667 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::operator const Number &() const
668 {
669  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
670  return value;
671 }
672 
673 
674 template <int dim, typename Number>
675 template <typename OtherNumber>
676 inline
678 {
679  value = p.value;
680  return *this;
681 }
682 
683 
684 template <int dim, typename Number>
685 template <typename OtherNumber>
686 inline
688 {
689  return (value == p.value);
690 }
691 
692 
693 template <int dim, typename Number>
694 template <typename OtherNumber>
695 inline
697 {
698  return !((*this) == p);
699 }
700 
701 
702 template <int dim, typename Number>
703 template <typename OtherNumber>
704 inline
706 {
707  value += p.value;
708  return *this;
709 }
710 
711 
712 template <int dim, typename Number>
713 template <typename OtherNumber>
714 inline
716 {
717  value -= p.value;
718  return *this;
719 }
720 
721 
722 template <int dim, typename Number>
723 template <typename OtherNumber>
724 inline
725 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator *= (const OtherNumber &s)
726 {
727  value *= s;
728  return *this;
729 }
730 
731 
732 template <int dim, typename Number>
733 template <typename OtherNumber>
734 inline
736 {
737  value /= s;
738  return *this;
739 }
740 
741 
742 template <int dim, typename Number>
743 inline
745 {
746  return -value;
747 }
748 
749 
750 template <int dim, typename Number>
751 inline
754 {
755  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
756  return numbers::NumberTraits<Number>::abs (value);
757 }
758 
759 
760 template <int dim, typename Number>
761 inline
763 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::norm_square () const
764 {
765  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
767 }
768 
769 
770 template <int dim, typename Number>
771 template <typename OtherNumber>
772 inline
773 void
775  unsigned int &index) const
776 {
777  Assert(dim != 0, ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
778  result[index] = value;
779  ++index;
780 }
781 
782 
783 template <int dim, typename Number>
784 inline
786 {
787  // Some auto-differentiable numbers need explicit
788  // zero initialization.
790 }
791 
792 
793 template <int dim, typename Number>
794 template <class Archive>
795 inline
796 void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int)
797 {
798  ar &value;
799 }
800 
801 
802 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
803 
804 
805 template <int rank_, int dim, typename Number>
806 inline
807 DEAL_II_CUDA_HOST_DEV Tensor<rank_,dim,Number>::Tensor ()
808 {
809  // All members of the c-style array values are already default initialized
810  // and thus all values are already set to zero recursively.
811 }
812 
813 
814 template <int rank_, int dim, typename Number>
815 inline
816 Tensor<rank_,dim,Number>::Tensor (const array_type &initializer)
817 {
818  for (unsigned int i=0; i<dim; ++i)
819  values[i] = Tensor<rank_-1, dim, Number>(initializer[i]);
820 }
821 
822 
823 template <int rank_, int dim, typename Number>
824 template <typename OtherNumber>
825 inline
827 {
828  for (unsigned int i=0; i!=dim; ++i)
829  values[i] = Tensor<rank_-1,dim,Number>(initializer[i]);
830 }
831 
832 
833 template <int rank_, int dim, typename Number>
834 template <typename OtherNumber>
835 inline
837 (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer)
838 {
839  for (unsigned int i=0; i<dim; ++i)
840  values[i] = initializer[i];
841 }
842 
843 
844 template <int rank_, int dim, typename Number>
845 template <typename OtherNumber>
846 inline
848 operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const
849 {
850  return Tensor<1,dim,Tensor<rank_-1,dim,Number> > (values);
851 }
852 
853 
854 
855 namespace internal
856 {
857  namespace TensorSubscriptor
858  {
859  template <typename ArrayElementType, int dim>
860  inline DEAL_II_ALWAYS_INLINE
861  DEAL_II_CUDA_HOST_DEV
862  ArrayElementType &
863  subscript (ArrayElementType *values,
864  const unsigned int i,
865  std::integral_constant<int, dim>)
866  {
867  Assert (i<dim, ExcIndexRange(i, 0, dim));
868  return values[i];
869  }
870 
871 
872  template <typename ArrayElementType>
873  DEAL_II_CUDA_HOST_DEV
874  ArrayElementType &
875  subscript (ArrayElementType *,
876  const unsigned int,
877  std::integral_constant<int, 0>)
878  {
879  Assert(false, ExcMessage("Cannot access elements of an object of type Tensor<rank,0,Number>."));
880  static ArrayElementType t;
881  return t;
882  }
883  }
884 }
885 
886 
887 template <int rank_, int dim, typename Number>
888 inline DEAL_II_ALWAYS_INLINE
889 DEAL_II_CUDA_HOST_DEV
892 {
893  return ::internal::TensorSubscriptor::subscript(values, i, std::integral_constant<int, dim>());
894 }
895 
896 
897 template <int rank_, int dim, typename Number>
898 inline DEAL_II_ALWAYS_INLINE
899 DEAL_II_CUDA_HOST_DEV
901 Tensor<rank_,dim,Number>::operator[] (const unsigned int i) const
902 {
903  return ::internal::TensorSubscriptor::subscript(values, i, std::integral_constant<int, dim>());
904 }
905 
906 
907 template <int rank_, int dim, typename Number>
908 inline
909 const Number &
911 {
912  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
913 
914  return TensorAccessors::extract<rank_>(*this, indices);
915 }
916 
917 
918 template <int rank_, int dim, typename Number>
919 inline
920 Number &
922 {
923  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
924 
925  return TensorAccessors::extract<rank_>(*this, indices);
926 }
927 
928 
929 template <int rank_, int dim, typename Number>
930 template <typename OtherNumber>
931 inline
934 {
935  if (dim > 0)
936  std::copy (&t.values[0], &t.values[0]+dim, &values[0]);
937  return *this;
938 }
939 
940 
941 template <int rank_, int dim, typename Number>
942 inline
945 {
947  ExcMessage ("Only assignment with zero is allowed"));
948  (void) d;
949 
950  for (unsigned int i=0; i<dim; ++i)
951  values[i] = internal::NumberType<Number>::value(0.0);
952  return *this;
953 }
954 
955 
956 template <int rank_, int dim, typename Number>
957 template <typename OtherNumber>
958 inline
959 bool
961 {
962  for (unsigned int i=0; i<dim; ++i)
963  if (values[i] != p.values[i])
964  return false;
965  return true;
966 }
967 
968 
969 // At some places in the library, we have Point<0> for formal reasons
970 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
971 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
972 // in the above function that the loop end check always fails, we
973 // implement this function here
974 template <>
975 template <>
976 inline
978 {
979  return true;
980 }
981 
982 
983 template <int rank_, int dim, typename Number>
984 template <typename OtherNumber>
985 inline
986 bool
988 {
989  return !((*this) == p);
990 }
991 
992 
993 template <int rank_, int dim, typename Number>
994 template <typename OtherNumber>
995 inline
998 {
999  for (unsigned int i=0; i<dim; ++i)
1000  values[i] += p.values[i];
1001  return *this;
1002 }
1003 
1004 
1005 template <int rank_, int dim, typename Number>
1006 template <typename OtherNumber>
1007 inline
1010 {
1011  for (unsigned int i=0; i<dim; ++i)
1012  values[i] -= p.values[i];
1013  return *this;
1014 }
1015 
1016 
1017 template <int rank_, int dim, typename Number>
1018 template <typename OtherNumber>
1019 inline
1020 DEAL_II_CUDA_HOST_DEV
1023 {
1024  for (unsigned int i=0; i<dim; ++i)
1025  values[i] *= s;
1026  return *this;
1027 }
1028 
1029 
1030 template <int rank_, int dim, typename Number>
1031 template <typename OtherNumber>
1032 inline
1035 {
1036  for (unsigned int i=0; i<dim; ++i)
1037  values[i] /= s;
1038  return *this;
1039 }
1040 
1041 
1042 template <int rank_, int dim, typename Number>
1043 inline
1046 {
1048 
1049  for (unsigned int i=0; i<dim; ++i)
1050  tmp.values[i] = -values[i];
1051 
1052  return tmp;
1053 }
1054 
1055 
1056 template <int rank_, int dim, typename Number>
1057 inline
1060 {
1061  return std::sqrt (norm_square());
1062 }
1063 
1064 
1065 template <int rank_, int dim, typename Number>
1066 inline
1067 DEAL_II_CUDA_HOST_DEV
1070 {
1073  for (unsigned int i=0; i<dim; ++i)
1074  s += values[i].norm_square();
1075 
1076  return s;
1077 }
1078 
1079 
1080 template <int rank_, int dim, typename Number>
1081 template <typename OtherNumber>
1082 inline
1083 void
1085 {
1086  AssertDimension (result.size(),(Utilities::fixed_power<rank_, unsigned int>(dim)));
1087 
1088  unsigned int index = 0;
1089  unroll_recursion (result, index);
1090 }
1091 
1092 
1093 template <int rank_, int dim, typename Number>
1094 template <typename OtherNumber>
1095 inline
1096 void
1098  unsigned int &index) const
1099 {
1100  for (unsigned int i=0; i<dim; ++i)
1101  values[i].unroll_recursion(result, index);
1102 }
1103 
1104 
1105 template <int rank_, int dim, typename Number>
1106 inline
1107 unsigned int
1109 {
1110  unsigned int index = 0;
1111  for (int r = 0; r < rank_; ++r)
1112  index = index * dim + indices[r];
1113 
1114  return index;
1115 }
1116 
1117 
1118 template <int rank_, int dim, typename Number>
1119 inline
1122 {
1123  Assert (i < n_independent_components,
1124  ExcIndexRange (i, 0, n_independent_components));
1125 
1126  TableIndices<rank_> indices;
1127 
1128  unsigned int remainder = i;
1129  for (int r=rank_-1; r>=0; --r)
1130  {
1131  indices[r] = (remainder % dim);
1132  remainder /= dim;
1133  }
1134  Assert (remainder == 0, ExcInternalError());
1135 
1136  return indices;
1137 }
1138 
1139 
1140 template <int rank_, int dim, typename Number>
1141 inline
1143 {
1144  for (unsigned int i=0; i<dim; ++i)
1145  values[i] = internal::NumberType<Number>::value(0.0);
1146 }
1147 
1148 
1149 template <int rank_, int dim, typename Number>
1150 inline
1151 std::size_t
1153 {
1154  return sizeof(Tensor<rank_,dim,Number>);
1155 }
1156 
1157 
1158 template <int rank_, int dim, typename Number>
1159 template <class Archive>
1160 inline
1161 void
1162 Tensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
1163 {
1164  ar &values;
1165 }
1166 
1167 
1168 /* ----------------- Non-member functions operating on tensors. ------------ */
1169 
1174 
1182 template <int rank_, int dim, typename Number>
1183 inline
1184 std::ostream &operator << (std::ostream &out, const Tensor<rank_,dim,Number> &p)
1185 {
1186  for (unsigned int i = 0; i < dim; ++i)
1187  {
1188  out << p[i];
1189  if (i != dim - 1)
1190  out << ' ';
1191  }
1192 
1193  return out;
1194 }
1195 
1196 
1203 template <int dim, typename Number>
1204 inline
1205 std::ostream &operator << (std::ostream &out, const Tensor<0,dim,Number> &p)
1206 {
1207  out << static_cast<const Number &>(p);
1208  return out;
1209 }
1210 
1211 
1213 
1217 
1218 
1227 template <int dim, typename Number, typename Other>
1228 inline
1230 operator * (const Other &object,
1231  const Tensor<0,dim,Number> &t)
1232 {
1233  return object * static_cast<const Number &>(t);
1234 }
1235 
1236 
1237 
1246 template <int dim, typename Number, typename Other>
1247 inline
1250  const Other &object)
1251 {
1252  return static_cast<const Number &>(t) * object;
1253 }
1254 
1255 
1265 template <int dim, typename Number, typename OtherNumber>
1266 inline
1269  const Tensor<0, dim, OtherNumber> &src2)
1270 {
1271  return static_cast<const Number &>(src1) *
1272  static_cast<const OtherNumber &>(src2);
1273 }
1274 
1275 
1281 template <int dim, typename Number, typename OtherNumber>
1282 inline
1285  const OtherNumber &factor)
1286 {
1287  return static_cast<const Number &>(t) / factor;
1288 }
1289 
1290 
1296 template <int dim, typename Number, typename OtherNumber>
1297 inline
1300  const Tensor<0,dim,OtherNumber> &q)
1301 {
1302  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1303 }
1304 
1305 
1311 template <int dim, typename Number, typename OtherNumber>
1312 inline
1315  const Tensor<0,dim,OtherNumber> &q)
1316 {
1317  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1318 }
1319 
1320 
1331 template <int rank, int dim,
1332  typename Number,
1333  typename OtherNumber>
1334 inline
1337  const OtherNumber &factor)
1338 {
1339  // recurse over the base objects
1341  for (unsigned int d=0; d<dim; ++d)
1342  tt[d] = t[d] * factor;
1343  return tt;
1344 }
1345 
1346 
1357 template <int rank, int dim,
1358  typename Number,
1359  typename OtherNumber>
1360 inline
1362 operator * (const Number &factor,
1364 {
1365  // simply forward to the operator above
1366  return t * factor;
1367 }
1368 
1369 
1377 template <int rank, int dim,
1378  typename Number,
1379  typename OtherNumber>
1380 inline
1383  const OtherNumber &factor)
1384 {
1385  // recurse over the base objects
1387  for (unsigned int d=0; d<dim; ++d)
1388  tt[d] = t[d] / factor;
1389  return tt;
1390 }
1391 
1392 
1400 template <int rank, int dim, typename Number, typename OtherNumber>
1401 inline
1405 {
1407 
1408  for (unsigned int i=0; i<dim; ++i)
1409  tmp[i] += q[i];
1410 
1411  return tmp;
1412 }
1413 
1414 
1422 template <int rank, int dim, typename Number, typename OtherNumber>
1423 inline
1427 {
1429 
1430  for (unsigned int i=0; i<dim; ++i)
1431  tmp[i] -= q[i];
1432 
1433  return tmp;
1434 }
1435 
1436 
1438 
1442 
1443 
1467 template <int rank_1, int rank_2, int dim,
1468  typename Number, typename OtherNumber>
1469 inline DEAL_II_ALWAYS_INLINE
1470 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1473 {
1474  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
1475 
1476  TensorAccessors::internal::ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1477  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
1478  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
1479 
1480  return result;
1481 }
1482 
1483 
1513 template <int index_1, int index_2,
1514  int rank_1, int rank_2, int dim,
1515  typename Number, typename OtherNumber>
1516 inline
1517 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1520 {
1521  Assert(0 <= index_1 && index_1 < rank_1,
1522  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1523  Assert(0 <= index_2 && index_2 < rank_2,
1524  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1525 
1526  using namespace TensorAccessors;
1527  using namespace TensorAccessors::internal;
1528 
1529  // Reorder index_1 to the end of src1:
1530  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1531  reord_01 = reordered_index_view<index_1, rank_1>(src1);
1532 
1533  // Reorder index_2 to the end of src2:
1534  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1535  reord_02 = reordered_index_view<index_2, rank_2>(src2);
1536 
1537  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1538  result;
1539  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
1540  return result;
1541 }
1542 
1543 
1575 template <int index_1, int index_2, int index_3, int index_4,
1576  int rank_1, int rank_2, int dim,
1577  typename Number, typename OtherNumber>
1578 inline
1579 typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1582 {
1583  Assert(0 <= index_1 && index_1 < rank_1,
1584  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1585  Assert(0 <= index_3 && index_3 < rank_1,
1586  ExcMessage("The specified index_3 must lie within the range [0,rank_1)"));
1587  Assert(index_1 != index_3,
1588  ExcMessage("index_1 and index_3 must not be the same"));
1589  Assert(0 <= index_2 && index_2 < rank_2,
1590  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1591  Assert(0 <= index_4 && index_4 < rank_2,
1592  ExcMessage("The specified index_4 must lie within the range [0,rank_2)"));
1593  Assert(index_2 != index_4,
1594  ExcMessage("index_2 and index_4 must not be the same"));
1595 
1596  using namespace TensorAccessors;
1597  using namespace TensorAccessors::internal;
1598 
1599  // Reorder index_1 to the end of src1:
1600  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1601  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
1602 
1603  // Reorder index_2 to the end of src2:
1604  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1605  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
1606 
1607  // Now, reorder index_3 to the end of src1. We have to make sure to
1608  // preserve the orginial ordering: index_1 has been removed. If
1609  // index_3 > index_1, we have to use (index_3 - 1) instead:
1610  ReorderedIndexView<(index_3 < index_1 ? index_3 : index_3 - 1), rank_1, ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> > >
1611  reord_3 = TensorAccessors::reordered_index_view<index_3 < index_1 ? index_3 : index_3 - 1, rank_1>(reord_1);
1612 
1613  // Now, reorder index_4 to the end of src2. We have to make sure to
1614  // preserve the orginial ordering: index_2 has been removed. If
1615  // index_4 > index_2, we have to use (index_4 - 1) instead:
1616  ReorderedIndexView<(index_4 < index_2 ? index_4 : index_4 - 1), rank_2, ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> > >
1617  reord_4 = TensorAccessors::reordered_index_view<index_4 < index_2 ? index_4 : index_4 - 1, rank_2>(reord_2);
1618 
1619  typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1620  result;
1621  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
1622  return result;
1623 }
1624 
1625 
1639 template <int rank, int dim, typename Number, typename OtherNumber>
1640 inline
1643  const Tensor<rank, dim, OtherNumber> &right)
1644 {
1646  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
1647  return result;
1648 }
1649 
1650 
1669 template <template <int, int, typename> class TensorT1,
1670  template <int, int, typename> class TensorT2,
1671  template <int, int, typename> class TensorT3,
1672  int rank_1, int rank_2, int dim,
1673  typename T1, typename T2, typename T3>
1675 contract3 (const TensorT1<rank_1, dim, T1> &left,
1676  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
1677  const TensorT3<rank_2, dim, T3> &right)
1678 {
1680  return_type;
1681  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(
1682  left, middle, right);
1683 }
1684 
1685 
1697 template <int rank_1, int rank_2, int dim,
1698  typename Number, typename OtherNumber>
1699 inline
1703 {
1705  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
1706  return result;
1707 }
1708 
1709 
1711 
1715 
1716 
1728 template <int dim, typename Number>
1729 inline
1732 {
1733  Assert (dim==2, ExcInternalError());
1734 
1735  Tensor<1, dim, Number> result;
1736 
1737  result[0] = src[1];
1738  result[1] = -src[0];
1739 
1740  return result;
1741 }
1742 
1743 
1754 template <int dim, typename Number>
1755 inline
1758  const Tensor<1,dim,Number> &src2)
1759 {
1760  Assert (dim==3, ExcInternalError());
1761 
1762  Tensor<1, dim, Number> result;
1763 
1764  result[0] = src1[1]*src2[2] - src1[2]*src2[1];
1765  result[1] = src1[2]*src2[0] - src1[0]*src2[2];
1766  result[2] = src1[0]*src2[1] - src1[1]*src2[0];
1767 
1768  return result;
1769 }
1770 
1771 
1773 
1777 
1778 
1785 template <int dim, typename Number>
1786 inline
1788 {
1789  // Compute the determinant using the Laplace expansion of the
1790  // determinant. We expand along the last row.
1791  Number det = internal::NumberType<Number>::value(0.0);
1792 
1793  for (unsigned int k=0; k<dim; ++k)
1794  {
1795  Tensor<2,dim-1,Number> minor;
1796  for (unsigned int i=0; i<dim-1; ++i)
1797  for (unsigned int j=0; j<dim-1; ++j)
1798  minor[i][j] = t[i][j<k ? j : j+1];
1799 
1800  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
1801 
1802  det += t[dim-1][k] * cofactor;
1803  }
1804 
1805  return ((dim % 2 == 0) ? 1. : -1.) * det;
1806 }
1807 
1813 template <typename Number>
1814 inline
1816 {
1817  return t[0][0];
1818 }
1819 
1820 
1828 template <int dim, typename Number>
1829 Number trace (const Tensor<2,dim,Number> &d)
1830 {
1831  Number t=d[0][0];
1832  for (unsigned int i=1; i<dim; ++i)
1833  t += d[i][i];
1834  return t;
1835 }
1836 
1837 
1847 template <int dim, typename Number>
1848 inline
1851 {
1852  Number return_tensor [dim][dim];
1853 
1854  // if desired, take over the
1855  // inversion of a 4x4 tensor
1856  // from the FullMatrix
1857  AssertThrow (false, ExcNotImplemented());
1858 
1859  return Tensor<2,dim,Number>(return_tensor);
1860 }
1861 
1862 
1863 #ifndef DOXYGEN
1864 
1865 template <typename Number>
1866 inline
1868 invert (const Tensor<2,1,Number> &t)
1869 {
1870  Number return_tensor [1][1];
1871 
1872  return_tensor[0][0] = 1.0/t[0][0];
1873 
1874  return Tensor<2,1,Number>(return_tensor);
1875 }
1876 
1877 
1878 template <typename Number>
1879 inline
1881 invert (const Tensor<2,2,Number> &t)
1882 {
1883  Tensor<2,2,Number> return_tensor;
1884 
1885  // this is Maple output,
1886  // thus a bit unstructured
1887  const Number inv_det_t = 1.0/(t[0][0]*t[1][1]-t[1][0]*t[0][1]);
1888  return_tensor[0][0] = t[1][1];
1889  return_tensor[0][1] = -t[0][1];
1890  return_tensor[1][0] = -t[1][0];
1891  return_tensor[1][1] = t[0][0];
1892  return_tensor *= inv_det_t;
1893 
1894  return return_tensor;
1895 }
1896 
1897 
1898 template <typename Number>
1899 inline
1901 invert (const Tensor<2,3,Number> &t)
1902 {
1903  Tensor<2,3,Number> return_tensor;
1904 
1905  const Number t4 = t[0][0]*t[1][1],
1906  t6 = t[0][0]*t[1][2],
1907  t8 = t[0][1]*t[1][0],
1908  t00 = t[0][2]*t[1][0],
1909  t01 = t[0][1]*t[2][0],
1910  t04 = t[0][2]*t[2][0],
1911  inv_det_t = 1.0/(t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
1912  t00*t[2][1]+t01*t[1][2]-t04*t[1][1]);
1913  return_tensor[0][0] = t[1][1]*t[2][2]-t[1][2]*t[2][1];
1914  return_tensor[0][1] = t[0][2]*t[2][1]-t[0][1]*t[2][2];
1915  return_tensor[0][2] = t[0][1]*t[1][2]-t[0][2]*t[1][1];
1916  return_tensor[1][0] = t[1][2]*t[2][0]-t[1][0]*t[2][2];
1917  return_tensor[1][1] = t[0][0]*t[2][2]-t04;
1918  return_tensor[1][2] = t00-t6;
1919  return_tensor[2][0] = t[1][0]*t[2][1]-t[1][1]*t[2][0];
1920  return_tensor[2][1] = t01-t[0][0]*t[2][1];
1921  return_tensor[2][2] = t4-t8;
1922  return_tensor *= inv_det_t;
1923 
1924  return return_tensor;
1925 }
1926 
1927 #endif /* DOXYGEN */
1928 
1929 
1936 template <int dim, typename Number>
1937 inline
1940 {
1942  for (unsigned int i=0; i<dim; ++i)
1943  {
1944  tt[i][i] = t[i][i];
1945  for (unsigned int j=i+1; j<dim; ++j)
1946  {
1947  tt[i][j] = t[j][i];
1948  tt[j][i] = t[i][j];
1949  };
1950  }
1951  return tt;
1952 }
1953 
1954 
1968 template <int dim, typename Number>
1969 inline
1972 {
1973  return determinant(t)*invert(t);
1974 }
1975 
1976 
1991 template <int dim, typename Number>
1992 inline
1995 {
1996  return transpose(adjugate(t));
1997 }
1998 
1999 
2007 template <int dim, typename Number>
2008 inline
2009 Number
2011 {
2012  Number max = internal::NumberType<Number>::value(0.0);
2013  for (unsigned int j=0; j<dim; ++j)
2014  {
2015  Number sum = internal::NumberType<Number>::value(0.0);
2016  for (unsigned int i=0; i<dim; ++i)
2017  sum += std::fabs(t[i][j]);
2018 
2019  if (sum > max)
2020  max = sum;
2021  }
2022 
2023  return max;
2024 }
2025 
2026 
2034 template <int dim, typename Number>
2035 inline
2036 Number
2038 {
2039  Number max = internal::NumberType<Number>::value(0.0);
2040  for (unsigned int i=0; i<dim; ++i)
2041  {
2042  Number sum = internal::NumberType<Number>::value(0.0);
2043  for (unsigned int j=0; j<dim; ++j)
2044  sum += std::fabs(t[i][j]);
2045 
2046  if (sum > max)
2047  max = sum;
2048  }
2049 
2050  return max;
2051 }
2052 
2054 
2055 DEAL_II_NAMESPACE_CLOSE
2056 
2057 // include deprecated non-member functions operating on Tensor
2058 #include <deal.II/base/tensor_deprecated.h>
2059 
2060 #endif
DEAL_II_CUDA_HOST_DEV Tensor()
Definition: tensor.h:807
numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:117
Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:1829
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1175
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
Definition: tensor.h:1108
Tensor< rank_-1, dim, Number >::array_type array_type[(dim!=0)?dim:1]
Definition: tensor.h:367
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:353
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2010
static std::size_t memory_consumption()
Definition: tensor.h:1152
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2037
Tensor< rank_, dim, Number > & operator-=(const Tensor< rank_, dim, OtherNumber > &)
Definition: tensor.h:1009
Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1971
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:1518
static const unsigned int rank
Definition: tensor.h:346
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1059
#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:987
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:565
void unroll(Vector< OtherNumber > &result) const
Definition: tensor.h:1084
Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1994
Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1939
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:1675
Tensor< rank_, dim, Number > & operator+=(const Tensor< rank_, dim, OtherNumber > &)
Definition: tensor.h:997
Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:1850
Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:1731
bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:960
#define Assert(cond, exc)
Definition: exceptions.h:337
Tensor< rank_, dim, Number > & operator/=(const OtherNumber &factor)
Definition: tensor.h:1034
ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:1642
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: tensor.h:1097
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:559
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:1162
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:1121
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:1580
Tensor< 1, dim, Number > cross_product_3d(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:1757
Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1787
static const unsigned int dimension
Definition: tensor.h:341
Tensor< rank_, dim, Number > operator-() const
Definition: tensor.h:1045
Definition: mpi.h:41
Tensor< rank_-1, dim, Number >::tensor_type value_type
Definition: tensor.h:360
DEAL_II_CUDA_HOST_DEV value_type & operator[](const unsigned int i)
Definition: tensor.h:891
static::ExceptionBase & ExcNotImplemented()
DEAL_II_CUDA_HOST_DEV numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1069
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:1701
void clear()
Definition: tensor.h:1142
static::ExceptionBase & ExcInternalError()
Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:1815
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:1022