Reference documentation for deal.II version 8.1.0
tensor_base.h
1 // ---------------------------------------------------------------------
2 // @f$Id: tensor_base.h 30040 2013-07-18 17:06:48Z maier @f$
3 //
4 // Copyright (C) 1998 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__tensor_base_h
18 #define __deal2__tensor_base_h
19 
20 
21 // single this file out from tensor.h, since we want to derive
22 // Point<dim,Number> from Tensor<1,dim,Number>. However, the point class will
23 // not need all the tensor stuff, so we don't want the whole tensor package to
24 // be included every time we use a point.
25 
26 
27 #include <deal.II/base/config.h>
29 #include <deal.II/base/table_indices.h>
30 #include <vector>
31 
32 #include <cmath>
33 #include <ostream>
34 
36 
37 template <typename number> class Vector;
38 
39 // forward declare Point and Tensor. This is the first definition of these
40 // classes and here we set the default Number type to double (this means that
41 // this file must be included when using something like Tensor<1,dim>, and
42 // Point and Tensor must not be forward declared without the number type
43 // specified)
44 template <int dim, typename Number=double> class Point;
45 
46 // general template; specialized for rank==1; the general template is in
47 // tensor.h
48 template <int rank_, int dim, typename Number=double> class Tensor;
49 template <int dim, typename Number> class Tensor<0,dim,Number>;
50 template <int dim, typename Number> class Tensor<1,dim,Number>;
51 
52 
53 
69 template <int dim, typename Number>
70 class Tensor<0,dim,Number>
71 {
72 public:
90  static const unsigned int dimension = dim;
91 
96  static const unsigned int rank = 0;
97 
103  typedef Number value_type;
104 
120 
124  Tensor ();
125 
131  Tensor (const value_type &initializer);
132 
136  Tensor (const Tensor<0,dim,Number> &);
137 
143  operator Number () const;
144 
154  operator Number &();
155 
160 
164  Tensor<0,dim,Number> &operator = (const Number d);
165 
170  bool operator == (const Tensor<0,dim,Number> &) const;
171 
176  bool operator != (const Tensor<0,dim,Number> &) const;
177 
184 
189 
195  Tensor<0,dim,Number> &operator *= (const Number factor);
196 
200  Tensor<0,dim,Number> &operator /= (const Number factor);
201 
206  Number operator * (const Tensor<0,dim,Number> &) const;
207 
215 
224 
229 
240  real_type norm () const;
241 
254  real_type norm_square () const;
255 
274  void clear ();
275 
286  DeclException1 (ExcDimTooSmall,
287  int,
288  << "dim must be positive, but was " << arg1);
289 
294  template <class Archive>
295  void serialize(Archive &ar, const unsigned int version);
296 
297 private:
301  Number value;
302 };
303 
304 
305 
328 template <int dim,typename Number>
329 class Tensor<1,dim,Number>
330 {
331 public:
349  static const unsigned int dimension = dim;
350 
355  static const unsigned int rank = 1;
356 
361  static const unsigned int
363 
369  typedef Number value_type;
370 
386 
398  typedef Number array_type[(dim!=0) ? dim : 100000000];
399 
405  explicit Tensor (const bool initialize = true);
406 
412  Tensor (const array_type &initializer);
413 
417  Tensor (const Tensor<1,dim,Number> &);
418 
429  Number operator [] (const unsigned int index) const;
430 
441  Number &operator [] (const unsigned int index);
442 
446  Number operator [](const TableIndices<1> &indices) const;
447 
451  Number &operator [](const TableIndices<1> &indices);
452 
457 
470  Tensor<1,dim,Number> &operator = (const Number d);
471 
476  bool operator == (const Tensor<1,dim,Number> &) const;
477 
482  bool operator != (const Tensor<1,dim,Number> &) const;
483 
490 
495 
501  Tensor<1,dim,Number> &operator *= (const Number factor);
502 
506  Tensor<1,dim,Number> &operator /= (const Number factor);
507 
512  Number operator * (const Tensor<1,dim,Number> &) const;
513 
521 
530 
535 
546  real_type norm () const;
547 
560  real_type norm_square () const;
561 
580  void clear ();
581 
591  template <typename Number2>
592  void unroll (Vector<Number2> &result) const;
593 
602  static
603  unsigned int
605 
615  static
616  TableIndices<1> unrolled_to_component_indices(const unsigned int i);
617 
618 
625  static std::size_t memory_consumption ();
626 
637  DeclException1 (ExcDimTooSmall,
638  int,
639  << "dim must be positive, but was " << arg1);
640 
645  template <class Archive>
646  void serialize(Archive &ar, const unsigned int version);
647 
648 private:
658  Number values[(dim!=0) ? (dim) : (dim+1)];
659 
671  template <typename Number2>
672  void unroll_recursion (Vector<Number2> &result,
673  unsigned int &start_index) const;
674 
675 private:
689  template <int otherrank, int otherdim, typename OtherNumber> friend class ::Tensor;
690 
696  friend class Point<dim,Number>;
697 };
698 
699 
703 template <int dim,typename Number>
704 std::ostream &operator << (std::ostream &out, const Tensor<0,dim,Number> &p);
705 
710 template <int dim,typename Number>
711 std::ostream &operator << (std::ostream &out, const Tensor<1,dim,Number> &p);
712 
713 
714 #ifndef DOXYGEN
715 
716 /*---------------------------- Inline functions: Tensor<0,dim> ------------------------*/
717 
718 template <int dim,typename Number>
719 inline
721 {
722  Assert (dim>0, ExcDimTooSmall(dim));
723 
724  value = 0;
725 }
726 
727 
728 
729 template <int dim, typename Number>
730 inline
731 Tensor<0,dim,Number>::Tensor (const value_type &initializer)
732 {
733  Assert (dim>0, ExcDimTooSmall(dim));
734 
735  value = initializer;
736 }
737 
738 
739 
740 template <int dim, typename Number>
741 inline
743 {
744  Assert (dim>0, ExcDimTooSmall(dim));
745 
746  value = p.value;
747 }
748 
749 
750 
751 
752 template <int dim, typename Number>
753 inline
754 Tensor<0,dim,Number>::operator Number () const
755 {
756  return value;
757 }
758 
759 
760 
761 template <int dim, typename Number>
762 inline
764 {
765  return value;
766 }
767 
768 
769 
770 template <int dim, typename Number>
771 inline
773 {
774  value = p.value;
775  return *this;
776 }
777 
778 
779 
780 template <int dim, typename Number>
781 inline
783 {
784  value = d;
785  return *this;
786 }
787 
788 
789 
790 template <int dim, typename Number>
791 inline
793 {
794  return (value == p.value);
795 }
796 
797 
798 
799 template <int dim, typename Number>
800 inline
802 {
803  return !((*this) == p);
804 }
805 
806 
807 
808 template <int dim, typename Number>
809 inline
811 {
812  value += p.value;
813  return *this;
814 }
815 
816 
817 
818 template <int dim, typename Number>
819 inline
821 {
822  value -= p.value;
823  return *this;
824 }
825 
826 
827 
828 template <int dim, typename Number>
829 inline
831 {
832  value *= s;
833  return *this;
834 }
835 
836 
837 
838 template <int dim, typename Number>
839 inline
841 {
842  value /= s;
843  return *this;
844 }
845 
846 
847 
848 template <int dim, typename Number>
849 inline
851 {
852  return value*p.value;
853 }
854 
855 
856 
857 template <int dim, typename Number>
858 inline
860 {
861  return value+p.value;
862 }
863 
864 
865 
866 template <int dim, typename Number>
867 inline
869 {
870  return value-p.value;
871 }
872 
873 
874 
875 template <int dim, typename Number>
876 inline
878 {
879  return -value;
880 }
881 
882 
883 
884 template <int dim, typename Number>
885 inline
888 {
889  return numbers::NumberTraits<Number>::abs (value);
890 }
891 
892 
893 
894 template <int dim, typename Number>
895 inline
898 {
900 }
901 
902 
903 
904 template <int dim, typename Number>
905 inline
907 {
908  value = 0;
909 }
910 
911 
912 
913 template <int dim, typename Number>
914 template <class Archive>
915 inline
916 void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int)
917 {
918  ar &value;
919 }
920 
921 /*---------------------------- Inline functions: Tensor<1,dim,Number> ------------------------*/
922 
923 
924 template <int dim, typename Number>
925 inline
926 Tensor<1,dim,Number>::Tensor (const bool initialize)
927 {
928  if (initialize)
929  // need to create an object Number() to
930  // initialize to zero to avoid confusion with
931  // Tensor::operator=(scalar) when using
932  // something like
933  // Tensor<1,dim,Tensor<1,dim,Number> >.
934  for (unsigned int i=0; i!=dim; ++i)
935  values[i] = Number();
936 }
937 
938 
939 
940 template <int dim, typename Number>
941 inline
942 Tensor<1,dim,Number>::Tensor (const array_type &initializer)
943 {
944  Assert (dim>0, ExcDimTooSmall(dim));
945 
946  for (unsigned int i=0; i<dim; ++i)
947  values[i] = initializer[i];
948 }
949 
950 
951 
952 template <int dim, typename Number>
953 inline
955 {
956  Assert (dim>0, ExcDimTooSmall(dim));
957 
958  for (unsigned int i=0; i<dim; ++i)
959  values[i] = p.values[i];
960 }
961 
962 
963 
964 template <>
965 inline
967 {
968  // at some places in the library,
969  // we have Point<0> for formal
970  // reasons (e.g., we sometimes have
971  // Quadrature<dim-1> for faces, so
972  // we have Quadrature<0> for dim=1,
973  // and then we have Point<0>). To
974  // avoid warnings in the above
975  // function that the loop end check
976  // always fails, we implement this
977  // function here
978 }
979 
980 
981 
982 template <int dim, typename Number>
983 inline
984 Number Tensor<1,dim,Number>::operator [] (const unsigned int index) const
985 {
986  Assert (index<dim, ExcIndexRange (index, 0, dim));
987  return values[index];
988 }
989 
990 
991 
992 template <int dim, typename Number>
993 inline
994 Number &Tensor<1,dim,Number>::operator [] (const unsigned int index)
995 {
996  Assert (index<dim, ExcIndexRange (index, 0, dim));
997  return values[index];
998 }
999 
1000 template <int dim, typename Number>
1001 inline
1002 Number Tensor<1,dim,Number>::operator [] (const TableIndices<1> &indices) const
1003 {
1004  Assert (indices[0]<dim, ExcIndexRange (indices[0], 0, dim));
1005  return values[indices[0]];
1006 }
1007 
1008 template <int dim, typename Number>
1009 inline
1011 {
1012  Assert (indices[0]<dim, ExcIndexRange (indices[0], 0, dim));
1013  return values[indices[0]];
1014 }
1015 
1016 
1017 
1018 template <>
1019 inline
1021 {
1022  // at some places in the library,
1023  // we have Point<0> for formal
1024  // reasons (e.g., we sometimes have
1025  // Quadrature<dim-1> for faces, so
1026  // we have Quadrature<0> for dim=1,
1027  // and then we have Point<0>). To
1028  // avoid warnings in the above
1029  // function that the loop end check
1030  // always fails, we implement this
1031  // function here
1032  return *this;
1033 }
1034 
1035 
1036 
1037 template <int dim, typename Number>
1038 inline
1041 {
1042  for (unsigned int i=0; i<dim; ++i)
1043  values[i] = p.values[i];
1044 
1045  return *this;
1046 }
1047 
1048 
1049 
1050 template <int dim, typename Number>
1051 inline
1053 {
1054  Assert (d==Number(0), ExcMessage ("Only assignment with zero is allowed"));
1055  (void) d;
1056 
1057  for (unsigned int i=0; i<dim; ++i)
1058  values[i] = 0;
1059 
1060  return *this;
1061 }
1062 
1063 
1064 
1065 template <int dim, typename Number>
1066 inline
1068 {
1069  for (unsigned int i=0; i<dim; ++i)
1070  if (values[i] != p.values[i])
1071  return false;
1072  return true;
1073 }
1074 
1075 
1076 
1077 template <>
1078 inline
1080 {
1081  return true;
1082 }
1083 
1084 
1085 
1086 template <int dim, typename Number>
1087 inline
1089 {
1090  return !((*this) == p);
1091 }
1092 
1093 
1094 
1095 template <int dim, typename Number>
1096 inline
1098 {
1099  for (unsigned int i=0; i<dim; ++i)
1100  values[i] += p.values[i];
1101  return *this;
1102 }
1103 
1104 
1105 
1106 template <int dim, typename Number>
1107 inline
1109 {
1110  for (unsigned int i=0; i<dim; ++i)
1111  values[i] -= p.values[i];
1112  return *this;
1113 }
1114 
1115 
1116 
1117 template <int dim, typename Number>
1118 inline
1120 {
1121  for (unsigned int i=0; i<dim; ++i)
1122  values[i] *= s;
1123  return *this;
1124 }
1125 
1126 
1127 
1128 template <int dim, typename Number>
1129 inline
1131 {
1132  for (unsigned int i=0; i<dim; ++i)
1133  values[i] /= s;
1134  return *this;
1135 }
1136 
1137 
1138 
1139 template <int dim, typename Number>
1140 inline
1141 Number
1143 {
1144  // unroll by hand since this is a
1145  // frequently called function and
1146  // some compilers don't want to
1147  // always unroll the loop in the
1148  // general template
1149  switch (dim)
1150  {
1151  case 1:
1152  return (values[0] * p.values[0]);
1153  break;
1154  case 2:
1155  return (values[0] * p.values[0] +
1156  values[1] * p.values[1]);
1157  break;
1158  case 3:
1159  return (values[0] * p.values[0] +
1160  values[1] * p.values[1] +
1161  values[2] * p.values[2]);
1162  break;
1163  default:
1164  Number q = values[0] * p.values[0];
1165  for (unsigned int i=1; i<dim; ++i)
1166  q += values[i] * p.values[i];
1167  return q;
1168  }
1169 }
1170 
1171 
1172 
1173 template <int dim, typename Number>
1174 inline
1176 {
1177  return (Tensor<1,dim,Number>(*this) += p);
1178 }
1179 
1180 
1181 
1182 template <int dim, typename Number>
1183 inline
1185 {
1186  return (Tensor<1,dim,Number>(*this) -= p);
1187 }
1188 
1189 
1190 
1191 template <int dim, typename Number>
1192 inline
1194 {
1195  Tensor<1,dim,Number> result (false);
1196  for (unsigned int i=0; i<dim; ++i)
1197  result.values[i] = -values[i];
1198  return result;
1199 }
1200 
1201 
1202 
1203 template <int dim, typename Number>
1204 inline
1207 {
1208  return std::sqrt (norm_square());
1209 }
1210 
1211 
1212 
1213 template <int dim, typename Number>
1214 inline
1217 {
1219  for (unsigned int i=1; i<dim; ++i)
1221 
1222  return s;
1223 }
1224 
1225 
1226 
1227 template <int dim, typename Number>
1228 template <typename Number2>
1229 inline
1230 void
1232 {
1233  Assert (result.size()==dim,
1234  ExcDimensionMismatch(dim, result.size()));
1235 
1236  unsigned int index = 0;
1237  unroll_recursion (result,index);
1238 }
1239 
1240 
1241 
1242 template<int dim, typename Number>
1243 template <typename Number2>
1244 inline
1245 void
1247  unsigned int &index) const
1248 {
1249  for (unsigned int i=0; i<dim; ++i)
1250  result(index++) = operator[](i);
1251 }
1252 
1253 
1254 template <int dim, typename Number>
1255 inline
1256 unsigned int
1258 {
1259  return indices[0];
1260 }
1261 
1262 template <int dim, typename Number>
1263 inline
1266 {
1267  return TableIndices<1>(i);
1268 }
1269 
1270 
1271 
1272 template <int dim, typename Number>
1273 inline
1275 {
1276  for (unsigned int i=0; i<dim; ++i)
1277  values[i] = 0;
1278 }
1279 
1280 
1281 
1282 template <int dim, typename Number>
1283 inline
1284 std::size_t
1286 {
1287  return sizeof(Tensor<1,dim,Number>);
1288 }
1289 
1290 
1291 template <int dim, typename Number>
1292 template <class Archive>
1293 inline
1294 void Tensor<1,dim,Number>::serialize(Archive &ar, const unsigned int)
1295 {
1296  ar &values;
1297 }
1298 #endif // DOXYGEN
1299 
1300 
1307 template <int dim, typename Number>
1308 inline
1309 std::ostream &operator << (std::ostream &out, const Tensor<0,dim,Number> &p)
1310 {
1311  out << static_cast<Number>(p);
1312  return out;
1313 }
1314 
1315 
1316 
1323 template <int dim, typename Number>
1324 inline
1325 std::ostream &operator << (std::ostream &out, const Tensor<1,dim,Number> &p)
1326 {
1327  for (unsigned int i=0; i<dim-1; ++i)
1328  out << p[i] << ' ';
1329  out << p[dim-1];
1330 
1331  return out;
1332 }
1333 
1334 
1335 
1343 inline
1344 std::ostream &operator << (std::ostream &out, const Tensor<1,1,double> &p)
1345 {
1346  out << p[0];
1347 
1348  return out;
1349 }
1350 
1351 
1352 
1358 template <int dim, typename Number>
1359 inline
1362  const Number factor)
1363 {
1364  Tensor<1,dim,Number> tt (false);
1365  for (unsigned int d=0; d<dim; ++d)
1366  tt[d] = t[d] * factor;
1367  return tt;
1368 }
1369 
1370 
1371 
1377 template <int dim, typename Number>
1378 inline
1380 operator * (const Number factor,
1381  const Tensor<1,dim,Number> &t)
1382 {
1383  Tensor<1,dim,Number> tt (false);
1384  for (unsigned int d=0; d<dim; ++d)
1385  tt[d] = t[d] * factor;
1386  return tt;
1387 }
1388 
1389 
1390 
1396 template <int dim, typename Number>
1397 inline
1400  const Number factor)
1401 {
1402  Tensor<1,dim,Number> tt (false);
1403  for (unsigned int d=0; d<dim; ++d)
1404  tt[d] = t[d] / factor;
1405  return tt;
1406 }
1407 
1408 
1409 
1415 template <int dim>
1416 inline
1419  const double factor)
1420 {
1421  Tensor<1,dim> tt (false);
1422  for (unsigned int d=0; d<dim; ++d)
1423  tt[d] = t[d] * factor;
1424  return tt;
1425 }
1426 
1427 
1428 
1434 template <int dim>
1435 inline
1437 operator * (const double factor,
1438  const Tensor<1,dim> &t)
1439 {
1440  Tensor<1,dim> tt (false);
1441  for (unsigned int d=0; d<dim; ++d)
1442  tt[d] = t[d] * factor;
1443  return tt;
1444 }
1445 
1446 
1447 
1453 template <int dim>
1454 inline
1457  const double factor)
1458 {
1459  Tensor<1,dim> tt (false);
1460  for (unsigned int d=0; d<dim; ++d)
1461  tt[d] = t[d] / factor;
1462  return tt;
1463 }
1464 
1465 
1466 DEAL_II_NAMESPACE_CLOSE
1467 
1468 #endif
1469 
1470 
VectorizedArray< Number > operator/(const VectorizedArray< Number > &u, const VectorizedArray< Number > &v)
numbers::NumberTraits< Number >::real_type real_type
Definition: tensor_base.h:119
numbers::NumberTraits< Number >::real_type real_type
Definition: tensor_base.h:385
Tensor< rank_, dim, Number > operator-() const
VectorizedArray< Number > operator*(const VectorizedArray< Number > &u, const VectorizedArray< Number > &v)
Number array_type[(dim!=0)?dim:100000000]
Definition: tensor_base.h:398
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
static const unsigned int n_independent_components
Definition: tensor.h:85
Tensor< 1, dim, Number > operator*(const Tensor< 1, dim, Number > &t, const Number factor)
Definition: tensor_base.h:1361
::ExceptionBase & ExcMessage(std::string arg1)
real_type norm_square() const
static const unsigned int rank
Definition: tensor.h:77
static real_type abs(const number &x)
Definition: numbers.h:316
Tensor< rank_, dim, Number > & operator-=(const Tensor< rank_, dim, Number > &)
Tensor< rank_-1, dim, Number > & operator[](const unsigned int i)
Number values[(dim!=0)?(dim):(dim+1)]
Definition: tensor_base.h:658
bool operator!=(const Tensor< rank_, dim, Number > &) const
static real_type abs_square(const number &x)
Definition: numbers.h:307
#define Assert(cond, exc)
Definition: exceptions.h:299
Tensor< rank_-1, dim, Number >::array_type array_type[dim]
Definition: tensor.h:115
static std::size_t memory_consumption()
void serialize(Archive &ar, const unsigned int version)
std::size_t size() const
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Tensor< rank_, dim, Number > operator+(const Tensor< rank_, dim, Number > &) const
void unroll_recursion(Vector< Number2 > &result, unsigned int &start_index) const
DeclException1(ExcInvalidTensorIndex, int,<< "Invalid tensor index "<< arg1)
bool operator==(const Tensor< rank_, dim, Number > &) const
static const unsigned int dimension
Definition: tensor.h:71
Tensor< rank_, dim, Number > & operator/=(const Number factor)
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Tensor< rank_, dim, Number > & operator+=(const Tensor< rank_, dim, Number > &)
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
void clear()
Tensor< rank_, dim, Number > & operator*=(const Number factor)
void unroll(Vector< Number2 > &result) const
real_type norm() const
Tensor & operator=(const Tensor< rank_, dim, Number > &)