Reference documentation for deal.II version Git f70953c 2018-04-22 22:20:09 +0200
tensor.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2018 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 
157  Number *
158  begin_raw();
159 
163  const Number *
164  begin_raw() const;
165 
169  Number *
170  end_raw();
171 
176  const Number *
177  end_raw() const;
178 
188  DEAL_II_CUDA_HOST_DEV operator Number &();
189 
198  DEAL_II_CUDA_HOST_DEV operator const Number &() const;
199 
205  template <typename OtherNumber>
207 
212 
217  template <typename OtherNumber>
218  Tensor &operator = (const OtherNumber &d);
219 
223  template <typename OtherNumber>
224  bool operator == (const Tensor<0,dim,OtherNumber> &rhs) const;
225 
229  template <typename OtherNumber>
230  bool operator != (const Tensor<0,dim,OtherNumber> &rhs) const;
231 
235  template <typename OtherNumber>
237 
241  template <typename OtherNumber>
243 
249  template <typename OtherNumber>
250  DEAL_II_CUDA_HOST_DEV Tensor &operator *= (const OtherNumber &factor);
251 
255  template <typename OtherNumber>
256  Tensor &operator /= (const OtherNumber &factor);
257 
261  Tensor operator - () const;
262 
275  void clear ();
276 
282  real_type norm () const;
283 
290  DEAL_II_CUDA_HOST_DEV real_type norm_square () const;
291 
296  template <class Archive>
297  void serialize(Archive &ar, const unsigned int version);
298 
303  typedef Number tensor_type;
304 
305 private:
309  Number value;
310 
314  template <typename OtherNumber>
316  unsigned int &start_index) const;
317 
321  template <int, int, typename> friend class Tensor;
322 };
323 
324 
325 
367 template <int rank_, int dim, typename Number>
368 class Tensor
369 {
370 public:
379  static const unsigned int dimension = dim;
380 
384  static const unsigned int rank = rank_;
385 
390  static const unsigned int
392 
398  typedef typename Tensor<rank_-1,dim,Number>::tensor_type value_type;
399 
404  typedef typename Tensor<rank_-1,dim,Number>::array_type
405  array_type[(dim != 0) ? dim : 1];
406  // ... avoid a compiler warning in case of dim == 0 and ensure that the
407  // array always has positive size.
408 
414  DEAL_II_CUDA_HOST_DEV Tensor ();
415 
419  explicit Tensor (const array_type &initializer);
420 
426  template <typename OtherNumber>
427  Tensor (const Tensor<rank_,dim,OtherNumber> &initializer);
428 
432  template <typename OtherNumber>
433  Tensor (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer);
434 
438  template <typename OtherNumber>
439  operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const;
440 
446  DEAL_II_CUDA_HOST_DEV value_type &operator [] (const unsigned int i);
447 
453  DEAL_II_CUDA_HOST_DEV const value_type &operator[](const unsigned int i) const;
454 
458  const Number &operator [] (const TableIndices<rank_> &indices) const;
459 
463  Number &operator [] (const TableIndices<rank_> &indices);
464 
468  Number *
469  begin_raw();
470 
474  const Number *
475  begin_raw() const;
476 
480  Number *
481  end_raw();
482 
486  const Number *
487  end_raw() const;
488 
494  template <typename OtherNumber>
496 
503  Tensor &operator = (const Number &d);
504 
508  template <typename OtherNumber>
509  bool operator == (const Tensor<rank_,dim,OtherNumber> &) const;
510 
514  template <typename OtherNumber>
515  bool operator != (const Tensor<rank_,dim,OtherNumber> &) const;
516 
520  template <typename OtherNumber>
522 
526  template <typename OtherNumber>
528 
535  template <typename OtherNumber>
536  DEAL_II_CUDA_HOST_DEV Tensor &operator *= (const OtherNumber &factor);
537 
541  template <typename OtherNumber>
542  Tensor &operator /= (const OtherNumber &factor);
543 
547  Tensor operator - () const;
548 
561  void clear ();
562 
570 
577  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type norm_square() const;
578 
586  template <typename OtherNumber>
587  void unroll (Vector<OtherNumber> &result) const;
588 
593  static
594  unsigned int
596 
601  static
603 
608  static std::size_t memory_consumption ();
609 
614  template <class Archive>
615  void serialize(Archive &ar, const unsigned int version);
616 
622 
623 private:
627  Tensor<rank_-1, dim, Number> values[(dim != 0) ? dim : 1];
628  // ... avoid a compiler warning in case of dim == 0 and ensure that the
629  // array always has positive size.
630 
634  template <typename OtherNumber>
636  unsigned int &start_index) const;
637 
641  template <int, int, typename> friend class Tensor;
642 
647  friend class Point<dim,Number>;
648 };
649 
650 
651 namespace internal
652 {
663  template <int rank, int dim, typename T>
664  struct NumberType<Tensor<rank,dim,T> >
665  {
666  static const Tensor<rank,dim,T> &value (const Tensor<rank,dim,T> &t)
667  {
668  return t;
669  }
670 
671  static Tensor<rank,dim,T> value (const T &t)
672  {
673  Tensor<rank,dim,T> tmp;
674  tmp=t;
675  return tmp;
676  }
677  };
678 
679  template <int rank, int dim, typename T>
680  struct NumberType<Tensor<rank,dim,VectorizedArray<T> > >
681  {
682  static const Tensor<rank,dim,VectorizedArray<T> > &value (const Tensor<rank,dim,VectorizedArray<T> > &t)
683  {
684  return t;
685  }
686 
687  static Tensor<rank,dim,VectorizedArray<T> > value (const T &t)
688  {
691  return tmp;
692  }
693 
695  {
697  tmp=t;
698  return tmp;
699  }
700  };
701 }
702 
703 
704 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
705 
706 
707 template <int dim,typename Number>
708 inline
709 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::Tensor ()
710 // Some auto-differentiable numbers need explicit
711 // zero initialization.
712  : value(internal::NumberType<Number>::value(0.0))
713 {
714 }
715 
716 
717 
718 template <int dim, typename Number>
719 template <typename OtherNumber>
720 inline
721 Tensor<0,dim,Number>::Tensor (const OtherNumber &initializer)
722 {
723  value = internal::NumberType<Number>::value(initializer);
724 }
725 
726 
727 
728 template <int dim, typename Number>
729 template <typename OtherNumber>
730 inline
732 {
733  value = p.value;
734 }
735 
736 
737 
738 template <int dim, typename Number>
739 inline
740 Number *
742 {
743  return std::addressof(value);
744 }
745 
746 
747 
748 template <int dim, typename Number>
749 inline
750 const Number *
752 {
753  return std::addressof(value);
754 }
755 
756 
757 
758 template <int dim, typename Number>
759 inline
760 Number *
762 {
764 }
765 
766 
767 
768 template <int dim, typename Number>
769 inline
770 const Number *
772 {
774 }
775 
776 
777 
778 template <int dim, typename Number>
779 inline
780 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::operator Number &()
781 {
782  // We cannot use Assert inside a CUDA kernel
783 #ifndef DEAL_II_WITH_CUDA
784  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
785 #endif
786  return value;
787 }
788 
789 
790 template <int dim, typename Number>
791 inline
792 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::operator const Number &() const
793 {
794  // We cannot use Assert inside a CUDA kernel
795 #ifndef DEAL_II_WITH_CUDA
796  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
797 #endif
798  return value;
799 }
800 
801 
802 template <int dim, typename Number>
803 template <typename OtherNumber>
804 inline
806 {
808  return *this;
809 }
810 
811 
812 template <int dim, typename Number>
813 inline
815 {
816  value = p.value;
817  return *this;
818 }
819 
820 
821 template <int dim, typename Number>
822 template <typename OtherNumber>
823 inline
825 {
827  return *this;
828 }
829 
830 
831 template <int dim, typename Number>
832 template <typename OtherNumber>
833 inline
835 {
836  return (value == p.value);
837 }
838 
839 
840 template <int dim, typename Number>
841 template <typename OtherNumber>
842 inline
844 {
845  return !((*this) == p);
846 }
847 
848 
849 template <int dim, typename Number>
850 template <typename OtherNumber>
851 inline
853 {
854  value += p.value;
855  return *this;
856 }
857 
858 
859 template <int dim, typename Number>
860 template <typename OtherNumber>
861 inline
863 {
864  value -= p.value;
865  return *this;
866 }
867 
868 
869 template <int dim, typename Number>
870 template <typename OtherNumber>
871 inline
872 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator *= (const OtherNumber &s)
873 {
874  value *= s;
875  return *this;
876 }
877 
878 
879 template <int dim, typename Number>
880 template <typename OtherNumber>
881 inline
883 {
884  value /= s;
885  return *this;
886 }
887 
888 
889 template <int dim, typename Number>
890 inline
892 {
893  return -value;
894 }
895 
896 
897 template <int dim, typename Number>
898 inline
901 {
902  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
903  return numbers::NumberTraits<Number>::abs (value);
904 }
905 
906 
907 template <int dim, typename Number>
908 inline
910 DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::norm_square () const
911 {
912  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
914 }
915 
916 
917 template <int dim, typename Number>
918 template <typename OtherNumber>
919 inline
920 void
922  unsigned int &index) const
923 {
924  Assert(dim != 0, ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
925  result[index] = value;
926  ++index;
927 }
928 
929 
930 template <int dim, typename Number>
931 inline
933 {
934  // Some auto-differentiable numbers need explicit
935  // zero initialization.
937 }
938 
939 
940 template <int dim, typename Number>
941 template <class Archive>
942 inline
943 void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int)
944 {
945  ar &value;
946 }
947 
948 
949 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
950 
951 
952 template <int rank_, int dim, typename Number>
953 inline DEAL_II_ALWAYS_INLINE
954 DEAL_II_CUDA_HOST_DEV Tensor<rank_,dim,Number>::Tensor ()
955 {
956  // All members of the c-style array values are already default initialized
957  // and thus all values are already set to zero recursively.
958 }
959 
960 
961 template <int rank_, int dim, typename Number>
962 inline DEAL_II_ALWAYS_INLINE
963 Tensor<rank_,dim,Number>::Tensor (const array_type &initializer)
964 {
965  for (unsigned int i=0; i<dim; ++i)
966  values[i] = Tensor<rank_-1, dim, Number>(initializer[i]);
967 }
968 
969 
970 template <int rank_, int dim, typename Number>
971 template <typename OtherNumber>
972 inline DEAL_II_ALWAYS_INLINE
974 {
975  for (unsigned int i=0; i!=dim; ++i)
976  values[i] = Tensor<rank_-1,dim,Number>(initializer[i]);
977 }
978 
979 
980 template <int rank_, int dim, typename Number>
981 template <typename OtherNumber>
982 inline DEAL_II_ALWAYS_INLINE
984 (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer)
985 {
986  for (unsigned int i=0; i<dim; ++i)
987  values[i] = initializer[i];
988 }
989 
990 
991 template <int rank_, int dim, typename Number>
992 template <typename OtherNumber>
993 inline DEAL_II_ALWAYS_INLINE
995 operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const
996 {
997  return Tensor<1,dim,Tensor<rank_-1,dim,Number> > (values);
998 }
999 
1000 
1001 
1002 namespace internal
1003 {
1004  namespace TensorSubscriptor
1005  {
1006  template <typename ArrayElementType, int dim>
1007  inline DEAL_II_ALWAYS_INLINE
1008  DEAL_II_CUDA_HOST_DEV
1009  ArrayElementType &
1010  subscript (ArrayElementType *values,
1011  const unsigned int i,
1012  std::integral_constant<int, dim>)
1013  {
1014  // We cannot use Assert in a CUDA kernel
1015 #ifndef DEAL_II_WITH_CUDA
1016  Assert (i<dim, ExcIndexRange(i, 0, dim));
1017 #endif
1018  return values[i];
1019  }
1020 
1021 
1022  template <typename ArrayElementType>
1023  DEAL_II_CUDA_HOST_DEV
1024  ArrayElementType &
1025  subscript (ArrayElementType *,
1026  const unsigned int,
1027  std::integral_constant<int, 0>)
1028  {
1029  // We cannot use Assert in a CUDA kernel
1030 #ifndef DEAL_II_WITH_CUDA
1031  Assert(false, ExcMessage("Cannot access elements of an object of type Tensor<rank,0,Number>."));
1032 #endif
1033  static ArrayElementType t;
1034  return t;
1035  }
1036  }
1037 }
1038 
1039 
1040 template <int rank_, int dim, typename Number>
1041 inline DEAL_II_ALWAYS_INLINE
1042 DEAL_II_CUDA_HOST_DEV
1045 {
1046  return ::internal::TensorSubscriptor::subscript(values, i, std::integral_constant<int, dim>());
1047 }
1048 
1049 
1050 template <int rank_, int dim, typename Number>
1051 inline DEAL_II_ALWAYS_INLINE
1052 DEAL_II_CUDA_HOST_DEV
1053 const typename Tensor<rank_,dim,Number>::value_type &
1054 Tensor<rank_,dim,Number>::operator[] (const unsigned int i) const
1055 {
1056  return ::internal::TensorSubscriptor::subscript(values, i, std::integral_constant<int, dim>());
1057 }
1058 
1059 
1060 template <int rank_, int dim, typename Number>
1061 inline
1062 const Number &
1064 {
1065  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1066 
1067  return TensorAccessors::extract<rank_>(*this, indices);
1068 }
1069 
1070 
1071 
1072 template <int rank_, int dim, typename Number>
1073 inline
1074 Number &
1076 {
1077  Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1078 
1079  return TensorAccessors::extract<rank_>(*this, indices);
1080 }
1081 
1082 
1083 
1084 template <int rank_, int dim, typename Number>
1085 inline
1086 Number *
1088 {
1089  return std::addressof(this->operator[](this->unrolled_to_component_indices(0)));
1090 }
1091 
1092 
1093 
1094 template <int rank_, int dim, typename Number>
1095 inline
1096 const Number *
1098 {
1099  return std::addressof(this->operator[](this->unrolled_to_component_indices(0)));
1100 }
1101 
1102 
1103 
1104 template <int rank_, int dim, typename Number>
1105 inline
1106 Number *
1108 {
1109  return begin_raw()+n_independent_components;
1110 }
1111 
1112 
1113 
1114 template <int rank_, int dim, typename Number>
1115 inline
1116 const Number *
1118 {
1119  return begin_raw()+n_independent_components;
1120 }
1121 
1122 
1123 
1124 template <int rank_, int dim, typename Number>
1125 template <typename OtherNumber>
1126 inline DEAL_II_ALWAYS_INLINE
1129 {
1130  if (dim > 0)
1131  std::copy (&t.values[0], &t.values[0]+dim, &values[0]);
1132  return *this;
1133 }
1134 
1135 
1136 template <int rank_, int dim, typename Number>
1137 inline DEAL_II_ALWAYS_INLINE
1140 {
1142  ExcMessage ("Only assignment with zero is allowed"));
1143  (void) d;
1144 
1145  for (unsigned int i=0; i<dim; ++i)
1146  values[i] = internal::NumberType<Number>::value(0.0);
1147  return *this;
1148 }
1149 
1150 
1151 template <int rank_, int dim, typename Number>
1152 template <typename OtherNumber>
1153 inline
1154 bool
1156 {
1157  for (unsigned int i=0; i<dim; ++i)
1158  if (values[i] != p.values[i])
1159  return false;
1160  return true;
1161 }
1162 
1163 
1164 // At some places in the library, we have Point<0> for formal reasons
1165 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1166 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1167 // in the above function that the loop end check always fails, we
1168 // implement this function here
1169 template <>
1170 template <>
1171 inline
1173 {
1174  return true;
1175 }
1176 
1177 
1178 template <int rank_, int dim, typename Number>
1179 template <typename OtherNumber>
1180 inline
1181 bool
1183 {
1184  return !((*this) == p);
1185 }
1186 
1187 
1188 template <int rank_, int dim, typename Number>
1189 template <typename OtherNumber>
1190 inline
1193 {
1194  for (unsigned int i=0; i<dim; ++i)
1195  values[i] += p.values[i];
1196  return *this;
1197 }
1198 
1199 
1200 template <int rank_, int dim, typename Number>
1201 template <typename OtherNumber>
1202 inline
1205 {
1206  for (unsigned int i=0; i<dim; ++i)
1207  values[i] -= p.values[i];
1208  return *this;
1209 }
1210 
1211 
1212 template <int rank_, int dim, typename Number>
1213 template <typename OtherNumber>
1214 inline
1215 DEAL_II_CUDA_HOST_DEV
1217 Tensor<rank_,dim,Number>::operator *= (const OtherNumber &s)
1218 {
1219  for (unsigned int i=0; i<dim; ++i)
1220  values[i] *= s;
1221  return *this;
1222 }
1223 
1224 
1225 template <int rank_, int dim, typename Number>
1226 template <typename OtherNumber>
1227 inline
1229 Tensor<rank_,dim,Number>::operator /= (const OtherNumber &s)
1230 {
1231  for (unsigned int i=0; i<dim; ++i)
1232  values[i] /= s;
1233  return *this;
1234 }
1235 
1236 
1237 template <int rank_, int dim, typename Number>
1238 inline
1241 {
1243 
1244  for (unsigned int i=0; i<dim; ++i)
1245  tmp.values[i] = -values[i];
1246 
1247  return tmp;
1248 }
1249 
1250 
1251 template <int rank_, int dim, typename Number>
1252 inline
1255 {
1256  return std::sqrt (norm_square());
1257 }
1258 
1259 
1260 template <int rank_, int dim, typename Number>
1261 inline
1262 DEAL_II_CUDA_HOST_DEV
1265 {
1268  for (unsigned int i=0; i<dim; ++i)
1269  s += values[i].norm_square();
1270 
1271  return s;
1272 }
1273 
1274 
1275 template <int rank_, int dim, typename Number>
1276 template <typename OtherNumber>
1277 inline
1278 void
1280 {
1281  AssertDimension (result.size(),(Utilities::fixed_power<rank_, unsigned int>(dim)));
1282 
1283  unsigned int index = 0;
1284  unroll_recursion (result, index);
1285 }
1286 
1287 
1288 template <int rank_, int dim, typename Number>
1289 template <typename OtherNumber>
1290 inline
1291 void
1293  unsigned int &index) const
1294 {
1295  for (unsigned int i=0; i<dim; ++i)
1296  values[i].unroll_recursion(result, index);
1297 }
1298 
1299 
1300 template <int rank_, int dim, typename Number>
1301 inline
1302 unsigned int
1304 {
1305  unsigned int index = 0;
1306  for (int r = 0; r < rank_; ++r)
1307  index = index * dim + indices[r];
1308 
1309  return index;
1310 }
1311 
1312 
1313 template <int rank_, int dim, typename Number>
1314 inline
1317 {
1318  Assert (i < n_independent_components,
1319  ExcIndexRange (i, 0, n_independent_components));
1320 
1321  TableIndices<rank_> indices;
1322 
1323  unsigned int remainder = i;
1324  for (int r=rank_-1; r>=0; --r)
1325  {
1326  indices[r] = (remainder % dim);
1327  remainder /= dim;
1328  }
1329  Assert (remainder == 0, ExcInternalError());
1330 
1331  return indices;
1332 }
1333 
1334 
1335 template <int rank_, int dim, typename Number>
1336 inline
1338 {
1339  for (unsigned int i=0; i<dim; ++i)
1340  values[i] = internal::NumberType<Number>::value(0.0);
1341 }
1342 
1343 
1344 template <int rank_, int dim, typename Number>
1345 inline
1346 std::size_t
1348 {
1349  return sizeof(Tensor<rank_,dim,Number>);
1350 }
1351 
1352 
1353 template <int rank_, int dim, typename Number>
1354 template <class Archive>
1355 inline
1356 void
1357 Tensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
1358 {
1359  ar &values;
1360 }
1361 
1362 
1363 /* ----------------- Non-member functions operating on tensors. ------------ */
1364 
1369 
1377 template <int rank_, int dim, typename Number>
1378 inline
1379 std::ostream &operator << (std::ostream &out, const Tensor<rank_,dim,Number> &p)
1380 {
1381  for (unsigned int i = 0; i < dim; ++i)
1382  {
1383  out << p[i];
1384  if (i != dim - 1)
1385  out << ' ';
1386  }
1387 
1388  return out;
1389 }
1390 
1391 
1398 template <int dim, typename Number>
1399 inline
1400 std::ostream &operator << (std::ostream &out, const Tensor<0,dim,Number> &p)
1401 {
1402  out << static_cast<const Number &>(p);
1403  return out;
1404 }
1405 
1406 
1408 
1412 
1413 
1422 template <int dim, typename Number, typename Other>
1423 inline DEAL_II_ALWAYS_INLINE
1424 typename ProductType<Other, Number>::type
1425 operator * (const Other &object,
1426  const Tensor<0,dim,Number> &t)
1427 {
1428  return object * static_cast<const Number &>(t);
1429 }
1430 
1431 
1432 
1441 template <int dim, typename Number, typename Other>
1442 inline DEAL_II_ALWAYS_INLINE
1443 typename ProductType<Number, Other>::type
1445  const Other &object)
1446 {
1447  return static_cast<const Number &>(t) * object;
1448 }
1449 
1450 
1460 template <int dim, typename Number, typename OtherNumber>
1461 inline DEAL_II_ALWAYS_INLINE
1462 typename ProductType<Number, OtherNumber>::type
1464  const Tensor<0, dim, OtherNumber> &src2)
1465 {
1466  return static_cast<const Number &>(src1) *
1467  static_cast<const OtherNumber &>(src2);
1468 }
1469 
1470 
1476 template <int dim, typename Number, typename OtherNumber>
1477 inline DEAL_II_ALWAYS_INLINE
1480  const OtherNumber &factor)
1481 {
1482  return static_cast<const Number &>(t) / factor;
1483 }
1484 
1485 
1491 template <int dim, typename Number, typename OtherNumber>
1492 inline DEAL_II_ALWAYS_INLINE
1495  const Tensor<0,dim,OtherNumber> &q)
1496 {
1497  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1498 }
1499 
1500 
1506 template <int dim, typename Number, typename OtherNumber>
1507 inline DEAL_II_ALWAYS_INLINE
1510  const Tensor<0,dim,OtherNumber> &q)
1511 {
1512  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1513 }
1514 
1515 
1526 template <int rank, int dim,
1527  typename Number,
1528  typename OtherNumber>
1529 inline DEAL_II_ALWAYS_INLINE
1532  const OtherNumber &factor)
1533 {
1534  // recurse over the base objects
1536  for (unsigned int d=0; d<dim; ++d)
1537  tt[d] = t[d] * factor;
1538  return tt;
1539 }
1540 
1541 
1552 template <int rank, int dim,
1553  typename Number,
1554  typename OtherNumber>
1555 inline DEAL_II_ALWAYS_INLINE
1557 operator * (const Number &factor,
1559 {
1560  // simply forward to the operator above
1561  return t * factor;
1562 }
1563 
1564 
1572 template <int rank, int dim,
1573  typename Number,
1574  typename OtherNumber>
1575 inline
1578  const OtherNumber &factor)
1579 {
1580  // recurse over the base objects
1582  for (unsigned int d=0; d<dim; ++d)
1583  tt[d] = t[d] / factor;
1584  return tt;
1585 }
1586 
1587 
1595 template <int rank, int dim, typename Number, typename OtherNumber>
1596 inline DEAL_II_ALWAYS_INLINE
1600 {
1602 
1603  for (unsigned int i=0; i<dim; ++i)
1604  tmp[i] += q[i];
1605 
1606  return tmp;
1607 }
1608 
1609 
1617 template <int rank, int dim, typename Number, typename OtherNumber>
1618 inline DEAL_II_ALWAYS_INLINE
1622 {
1624 
1625  for (unsigned int i=0; i<dim; ++i)
1626  tmp[i] -= q[i];
1627 
1628  return tmp;
1629 }
1630 
1631 
1633 
1637 
1638 
1662 template <int rank_1, int rank_2, int dim,
1663  typename Number, typename OtherNumber>
1664 inline DEAL_II_ALWAYS_INLINE
1665 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1668 {
1669  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
1670 
1671  TensorAccessors::internal::ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1672  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
1673  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
1674 
1675  return result;
1676 }
1677 
1678 
1708 template <int index_1, int index_2,
1709  int rank_1, int rank_2, int dim,
1710  typename Number, typename OtherNumber>
1711 inline DEAL_II_ALWAYS_INLINE
1712 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1715 {
1716  Assert(0 <= index_1 && index_1 < rank_1,
1717  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1718  Assert(0 <= index_2 && index_2 < rank_2,
1719  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1720 
1721  using namespace TensorAccessors;
1722  using namespace TensorAccessors::internal;
1723 
1724  // Reorder index_1 to the end of src1:
1725  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1726  reord_01 = reordered_index_view<index_1, rank_1>(src1);
1727 
1728  // Reorder index_2 to the end of src2:
1729  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1730  reord_02 = reordered_index_view<index_2, rank_2>(src2);
1731 
1732  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1733  result;
1734  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
1735  return result;
1736 }
1737 
1738 
1770 template <int index_1, int index_2, int index_3, int index_4,
1771  int rank_1, int rank_2, int dim,
1772  typename Number, typename OtherNumber>
1773 inline
1774 typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1777 {
1778  Assert(0 <= index_1 && index_1 < rank_1,
1779  ExcMessage("The specified index_1 must lie within the range [0,rank_1)"));
1780  Assert(0 <= index_3 && index_3 < rank_1,
1781  ExcMessage("The specified index_3 must lie within the range [0,rank_1)"));
1782  Assert(index_1 != index_3,
1783  ExcMessage("index_1 and index_3 must not be the same"));
1784  Assert(0 <= index_2 && index_2 < rank_2,
1785  ExcMessage("The specified index_2 must lie within the range [0,rank_2)"));
1786  Assert(0 <= index_4 && index_4 < rank_2,
1787  ExcMessage("The specified index_4 must lie within the range [0,rank_2)"));
1788  Assert(index_2 != index_4,
1789  ExcMessage("index_2 and index_4 must not be the same"));
1790 
1791  using namespace TensorAccessors;
1792  using namespace TensorAccessors::internal;
1793 
1794  // Reorder index_1 to the end of src1:
1795  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> >
1796  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
1797 
1798  // Reorder index_2 to the end of src2:
1799  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> >
1800  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
1801 
1802  // Now, reorder index_3 to the end of src1. We have to make sure to
1803  // preserve the orginial ordering: index_1 has been removed. If
1804  // index_3 > index_1, we have to use (index_3 - 1) instead:
1805  ReorderedIndexView<(index_3 < index_1 ? index_3 : index_3 - 1), rank_1, ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number> > >
1806  reord_3 = TensorAccessors::reordered_index_view<index_3 < index_1 ? index_3 : index_3 - 1, rank_1>(reord_1);
1807 
1808  // Now, reorder index_4 to the end of src2. We have to make sure to
1809  // preserve the orginial ordering: index_2 has been removed. If
1810  // index_4 > index_2, we have to use (index_4 - 1) instead:
1811  ReorderedIndexView<(index_4 < index_2 ? index_4 : index_4 - 1), rank_2, ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber> > >
1812  reord_4 = TensorAccessors::reordered_index_view<index_4 < index_2 ? index_4 : index_4 - 1, rank_2>(reord_2);
1813 
1814  typename Tensor<rank_1 + rank_2 - 4, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
1815  result;
1816  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
1817  return result;
1818 }
1819 
1820 
1834 template <int rank, int dim, typename Number, typename OtherNumber>
1835 inline DEAL_II_ALWAYS_INLINE
1836 typename ProductType<Number, OtherNumber>::type
1838  const Tensor<rank, dim, OtherNumber> &right)
1839 {
1840  typename ProductType<Number, OtherNumber>::type result;
1841  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
1842  return result;
1843 }
1844 
1845 
1864 template <template <int, int, typename> class TensorT1,
1865  template <int, int, typename> class TensorT2,
1866  template <int, int, typename> class TensorT3,
1867  int rank_1, int rank_2, int dim,
1868  typename T1, typename T2, typename T3>
1870 contract3 (const TensorT1<rank_1, dim, T1> &left,
1871  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
1872  const TensorT3<rank_2, dim, T3> &right)
1873 {
1875  return_type;
1876  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(
1877  left, middle, right);
1878 }
1879 
1880 
1892 template <int rank_1, int rank_2, int dim,
1893  typename Number, typename OtherNumber>
1894 inline DEAL_II_ALWAYS_INLINE
1898 {
1900  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
1901  return result;
1902 }
1903 
1904 
1906 
1910 
1911 
1923 template <int dim, typename Number>
1924 inline DEAL_II_ALWAYS_INLINE
1927 {
1928  Assert (dim==2, ExcInternalError());
1929 
1930  Tensor<1, dim, Number> result;
1931 
1932  result[0] = src[1];
1933  result[1] = -src[0];
1934 
1935  return result;
1936 }
1937 
1938 
1949 template <int dim, typename Number>
1950 inline DEAL_II_ALWAYS_INLINE
1953  const Tensor<1,dim,Number> &src2)
1954 {
1955  Assert (dim==3, ExcInternalError());
1956 
1957  Tensor<1, dim, Number> result;
1958 
1959  result[0] = src1[1]*src2[2] - src1[2]*src2[1];
1960  result[1] = src1[2]*src2[0] - src1[0]*src2[2];
1961  result[2] = src1[0]*src2[1] - src1[1]*src2[0];
1962 
1963  return result;
1964 }
1965 
1966 
1968 
1972 
1973 
1980 template <int dim, typename Number>
1981 inline
1983 {
1984  // Compute the determinant using the Laplace expansion of the
1985  // determinant. We expand along the last row.
1986  Number det = internal::NumberType<Number>::value(0.0);
1987 
1988  for (unsigned int k=0; k<dim; ++k)
1989  {
1990  Tensor<2,dim-1,Number> minor;
1991  for (unsigned int i=0; i<dim-1; ++i)
1992  for (unsigned int j=0; j<dim-1; ++j)
1993  minor[i][j] = t[i][j<k ? j : j+1];
1994 
1995  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
1996 
1997  det += t[dim-1][k] * cofactor;
1998  }
1999 
2000  return ((dim % 2 == 0) ? 1. : -1.) * det;
2001 }
2002 
2008 template <typename Number>
2009 inline
2011 {
2012  return t[0][0];
2013 }
2014 
2015 
2023 template <int dim, typename Number>
2024 inline DEAL_II_ALWAYS_INLINE
2025 Number trace (const Tensor<2,dim,Number> &d)
2026 {
2027  Number t=d[0][0];
2028  for (unsigned int i=1; i<dim; ++i)
2029  t += d[i][i];
2030  return t;
2031 }
2032 
2033 
2043 template <int dim, typename Number>
2044 inline
2047 {
2048  Number return_tensor [dim][dim];
2049 
2050  // if desired, take over the
2051  // inversion of a 4x4 tensor
2052  // from the FullMatrix
2053  AssertThrow (false, ExcNotImplemented());
2054 
2055  return Tensor<2,dim,Number>(return_tensor);
2056 }
2057 
2058 
2059 #ifndef DOXYGEN
2060 
2061 template <typename Number>
2062 inline
2064 invert (const Tensor<2,1,Number> &t)
2065 {
2066  Number return_tensor [1][1];
2067 
2068  return_tensor[0][0] = internal::NumberType<Number>::value(1.0/t[0][0]);
2069 
2070  return Tensor<2,1,Number>(return_tensor);
2071 }
2072 
2073 
2074 template <typename Number>
2075 inline
2077 invert (const Tensor<2,2,Number> &t)
2078 {
2079  Tensor<2,2,Number> return_tensor;
2080 
2081  // this is Maple output,
2082  // thus a bit unstructured
2083  const Number inv_det_t = internal::NumberType<Number>::value(1.0/(t[0][0]*t[1][1]-t[1][0]*t[0][1]));
2084  return_tensor[0][0] = t[1][1];
2085  return_tensor[0][1] = -t[0][1];
2086  return_tensor[1][0] = -t[1][0];
2087  return_tensor[1][1] = t[0][0];
2088  return_tensor *= inv_det_t;
2089 
2090  return return_tensor;
2091 }
2092 
2093 
2094 template <typename Number>
2095 inline
2097 invert (const Tensor<2,3,Number> &t)
2098 {
2099  Tensor<2,3,Number> return_tensor;
2100 
2101  const Number t4 = internal::NumberType<Number>::value(t[0][0]*t[1][1]),
2102  t6 = internal::NumberType<Number>::value(t[0][0]*t[1][2]),
2103  t8 = internal::NumberType<Number>::value(t[0][1]*t[1][0]),
2104  t00 = internal::NumberType<Number>::value(t[0][2]*t[1][0]),
2105  t01 = internal::NumberType<Number>::value(t[0][1]*t[2][0]),
2106  t04 = internal::NumberType<Number>::value(t[0][2]*t[2][0]),
2108  1.0/(t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
2109  t00*t[2][1]+t01*t[1][2]-t04*t[1][1]));
2110  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1]*t[2][2])-internal::NumberType<Number>::value(t[1][2]*t[2][1]);
2111  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2]*t[2][1])-internal::NumberType<Number>::value(t[0][1]*t[2][2]);
2112  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1]*t[1][2])-internal::NumberType<Number>::value(t[0][2]*t[1][1]);
2113  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2]*t[2][0])-internal::NumberType<Number>::value(t[1][0]*t[2][2]);
2114  return_tensor[1][1] = internal::NumberType<Number>::value(t[0][0]*t[2][2])-t04;
2115  return_tensor[1][2] = t00-t6;
2116  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0]*t[2][1])-internal::NumberType<Number>::value(t[1][1]*t[2][0]);
2117  return_tensor[2][1] = t01-internal::NumberType<Number>::value(t[0][0]*t[2][1]);
2118  return_tensor[2][2] = internal::NumberType<Number>::value(t4-t8);
2119  return_tensor *= inv_det_t;
2120 
2121  return return_tensor;
2122 }
2123 
2124 #endif /* DOXYGEN */
2125 
2126 
2133 template <int dim, typename Number>
2134 inline DEAL_II_ALWAYS_INLINE
2137 {
2139  for (unsigned int i=0; i<dim; ++i)
2140  {
2141  tt[i][i] = t[i][i];
2142  for (unsigned int j=i+1; j<dim; ++j)
2143  {
2144  tt[i][j] = t[j][i];
2145  tt[j][i] = t[i][j];
2146  };
2147  }
2148  return tt;
2149 }
2150 
2151 
2165 template <int dim, typename Number>
2166 inline
2169 {
2170  return determinant(t)*invert(t);
2171 }
2172 
2173 
2188 template <int dim, typename Number>
2189 inline
2192 {
2193  return transpose(adjugate(t));
2194 }
2195 
2196 
2204 template <int dim, typename Number>
2205 inline
2206 Number
2208 {
2209  Number max = internal::NumberType<Number>::value(0.0);
2210  for (unsigned int j=0; j<dim; ++j)
2211  {
2212  Number sum = internal::NumberType<Number>::value(0.0);
2213  for (unsigned int i=0; i<dim; ++i)
2214  sum += std::fabs(t[i][j]);
2215 
2216  if (sum > max)
2217  max = sum;
2218  }
2219 
2220  return max;
2221 }
2222 
2223 
2231 template <int dim, typename Number>
2232 inline
2233 Number
2235 {
2236  Number max = internal::NumberType<Number>::value(0.0);
2237  for (unsigned int i=0; i<dim; ++i)
2238  {
2239  Number sum = internal::NumberType<Number>::value(0.0);
2240  for (unsigned int j=0; j<dim; ++j)
2241  sum += std::fabs(t[i][j]);
2242 
2243  if (sum > max)
2244  max = sum;
2245  }
2246 
2247  return max;
2248 }
2249 
2251 
2252 DEAL_II_NAMESPACE_CLOSE
2253 
2254 // include deprecated non-member functions operating on Tensor
2255 #include <deal.II/base/tensor_deprecated.h>
2256 
2257 #endif
SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:118
Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2025
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1217
Number determinant(const SymmetricTensor< 2, dim, Number > &)
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
Definition: tensor.h:1303
Tensor< rank_-1, dim, Number >::array_type array_type[(dim!=0)?dim:1]
Definition: tensor.h:405
static const unsigned int n_independent_components
Definition: tensor.h:391
size_type size() const
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2207
static std::size_t memory_consumption()
Definition: tensor.h:1347
Tensor & operator/=(const OtherNumber &factor)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2234
Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2168
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:1713
static const unsigned int rank
Definition: tensor.h:384
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1254
#define AssertThrow(cond, exc)
Definition: exceptions.h:1190
static real_type abs(const number &x)
Definition: numbers.h:365
bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1182
internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
Definition: point.h:104
Tensor< rank_-1, dim, Number > values[(dim!=0)?dim:1]
Definition: tensor.h:627
void unroll(Vector< OtherNumber > &result) const
Definition: tensor.h:1279
Tensor & operator*=(const OtherNumber &factor)
Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2191
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2136
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:1870
Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:2046
static real_type abs_square(const number &x)
Definition: numbers.h:356
Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:1926
bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1155
#define Assert(cond, exc)
Definition: exceptions.h:1113
ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:1837
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: tensor.h:1292
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:621
void serialize(Archive &ar, const unsigned int version)
Definition: tensor.h:1357
Number * end_raw()
Definition: tensor.h:1107
value_type & operator[](const unsigned int i)
Definition: tensor.h:1044
Tensor()
Definition: tensor.h:954
Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Definition: tensor.h:1316
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:1775
Tensor< 1, dim, Number > cross_product_3d(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:1952
Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1982
static const unsigned int dimension
Definition: tensor.h:379
Tensor operator-() const
Definition: tensor.h:1240
Definition: mpi.h:53
Tensor< rank_-1, dim, Number >::tensor_type value_type
Definition: tensor.h:398
Number * begin_raw()
Definition: tensor.h:1087
static::ExceptionBase & ExcNotImplemented()
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:1896
numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1264
Tensor & operator+=(const Tensor< rank_, dim, OtherNumber > &)
void clear()
Definition: tensor.h:1337
static::ExceptionBase & ExcInternalError()
Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:2010
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const