Reference documentation for deal.II version 8.1.0
tensor.h
1 // ---------------------------------------------------------------------
2 // @f$Id: tensor.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_h
18 #define __deal2__tensor_h
19 
20 
21 #include <deal.II/base/config.h>
22 #include <deal.II/base/tensor_base.h>
23 #include <deal.II/base/utilities.h>
24 
26 template <int rank_, int dim, typename Number> class Tensor;
27 template <int dim, typename Number> class Tensor<1,dim,Number>;
28 
50 template <int rank_, int dim, typename Number>
51 class Tensor
52 {
53 public:
71  static const unsigned int dimension = dim;
72 
77  static const unsigned int rank = rank_;
78 
84  static const unsigned int
86 
91  typedef Tensor<rank_-1,dim,Number> value_type;
92 
108 
115  typedef typename Tensor<rank_-1,dim,Number>::array_type array_type[dim];
116 
121  Tensor ();
122 
127  Tensor (const array_type &initializer);
128 
133  Tensor (const Tensor<1,dim,Tensor<rank_-1,dim,Number> > &tensor_in);
134 
139  operator Tensor<1,dim,Tensor<rank_-1,dim,Number> > () const;
140 
144  Tensor<rank_-1,dim,Number> &operator [] (const unsigned int i);
145 
149  const Tensor<rank_-1,dim,Number> &operator [] (const unsigned int i) const;
150 
154  Number operator [](const TableIndices<rank_> &indices) const;
155 
159  Number &operator [](const TableIndices<rank_> &indices);
160 
165 
178  Tensor<rank_,dim,Number> &operator = (const Number d);
179 
183  bool operator == (const Tensor<rank_,dim,Number> &) const;
184 
188  bool operator != (const Tensor<rank_,dim,Number> &) const;
189 
194 
199 
205  Tensor<rank_,dim,Number> &operator *= (const Number factor);
206 
211  Tensor<rank_,dim,Number> &operator /= (const Number factor);
212 
220 
228 
234 
240  real_type norm () const;
241 
254  real_type norm_square () const;
255 
265  template <typename Number2>
266  void unroll (Vector<Number2> &result) const;
267 
273  static
274  unsigned int
276 
282  static
284 
285 
286 
305  void clear ();
306 
312  static std::size_t memory_consumption ();
313 
317  DeclException1 (ExcInvalidTensorIndex,
318  int,
319  << "Invalid tensor index " << arg1);
320 
325  template <class Archive>
326  void serialize(Archive &ar, const unsigned int version);
327 
328 private:
333  Tensor<rank_-1,dim,Number> subtensor[dim];
334 
338  template <typename Number2>
339  void unroll_recursion(Vector<Number2> &result,
340  unsigned int &start_index) const;
341 
342  // make the following class a
343  // friend to this class. in principle,
344  // it would suffice if otherrank==rank+1,
345  // but then the compiler complains
346  // that this be an explicit specialization
347  // which is not what we want
348  //
349  // also, it would be sufficient to make
350  // the function unroll_loops a friend,
351  // but that seems to be impossible as well.
352  template <int, int, typename> friend class Tensor;
353 };
354 
355 
356 /*--------------------------- Inline functions -----------------------------*/
357 
358 #ifndef DOXYGEN
359 
360 template <int rank_, int dim, typename Number>
361 inline
363 {
364 // default constructor. not specifying an initializer list calls
365 // the default constructor of the subobjects, which initialize them
366 // selves. therefore, the tensor is set to zero this way
367 }
368 
369 
370 
371 template <int rank_, int dim, typename Number>
372 inline
373 Tensor<rank_,dim,Number>::Tensor (const array_type &initializer)
374 {
375  for (unsigned int i=0; i<dim; ++i)
376  subtensor[i] = Tensor<rank_-1,dim,Number>(initializer[i]);
377 }
378 
379 
380 
381 template <int rank_, int dim, typename Number>
382 inline
384 (const Tensor<1,dim,Tensor<rank_-1,dim,Number> > &tensor_in)
385 {
386  for (unsigned int i=0; i<dim; ++i)
387  subtensor[i] = tensor_in[i];
388 }
389 
390 
391 
392 template <int rank_, int dim, typename Number>
393 inline
395 Tensor<1,dim,Tensor<rank_-1,dim,Number> > () const
396 {
397  return Tensor<1,dim,Tensor<rank_-1,dim,Number> > (subtensor);
398 }
399 
400 
401 
402 template <int rank_, int dim, typename Number>
403 inline
405 Tensor<rank_,dim,Number>::operator[] (const unsigned int i)
406 {
407  Assert (i<dim, ExcIndexRange(i, 0, dim));
408  return subtensor[i];
409 }
410 
411 
412 
413 template <int rank_, int dim, typename Number>
414 inline
416 Tensor<rank_,dim,Number>::operator[] (const unsigned int i) const
417 {
418  Assert (i<dim, ExcIndexRange(i, 0, dim));
419 
420  return subtensor[i];
421 }
422 
423 template <int rank_, int dim, typename Number>
424 inline
425 Number
427 {
428  const unsigned int inner_ind = indices[0];
429  Assert (inner_ind<dim, ExcIndexRange(inner_ind, 0, dim));
430 
431  TableIndices<rank_-1> indices1;
432  for (unsigned int i = 0; i < rank_-1; i++)
433  indices1[i] = indices[i+1];
434  return (subtensor[inner_ind])[indices1];
435 }
436 
437 template <int rank_, int dim, typename Number>
438 inline
439 Number &
441 {
442  const unsigned int inner_ind = indices[0];
443  Assert (inner_ind<dim, ExcIndexRange(inner_ind, 0, dim));
444 
445  TableIndices<rank_-1> indices1;
446  for (unsigned int i = 0; i < rank_-1; i++)
447  indices1[i] = indices[i+1];
448  return (subtensor[inner_ind])[indices1];
449 }
450 
451 template <int rank_, int dim, typename Number>
452 inline
455 {
456  for (unsigned int i=0; i<dim; ++i)
457  subtensor[i] = t.subtensor[i];
458  return *this;
459 }
460 
461 
462 
463 template <int rank_, int dim, typename Number>
464 inline
467 {
468  Assert (d==Number(0), ExcMessage ("Only assignment with zero is allowed"));
469  (void) d;
470 
471  for (unsigned int i=0; i<dim; ++i)
472  subtensor[i] = 0;
473  return *this;
474 }
475 
476 
477 
478 template <int rank_, int dim, typename Number>
479 inline
480 bool
482 {
483  for (unsigned int i=0; i<dim; ++i)
484  if (subtensor[i] != p.subtensor[i]) return false;
485  return true;
486 }
487 
488 
489 
490 template <int rank_, int dim, typename Number>
491 inline
492 bool
494 {
495  return !((*this) == p);
496 }
497 
498 
499 
500 template <int rank_, int dim, typename Number>
501 inline
504 {
505  for (unsigned int i=0; i<dim; ++i)
506  subtensor[i] += p.subtensor[i];
507  return *this;
508 }
509 
510 
511 
512 template <int rank_, int dim, typename Number>
513 inline
516 {
517  for (unsigned int i=0; i<dim; ++i)
518  subtensor[i] -= p.subtensor[i];
519  return *this;
520 }
521 
522 
523 
524 template <int rank_, int dim, typename Number>
525 inline
528 {
529  for (unsigned int i=0; i<dim; ++i)
530  subtensor[i] *= s;
531  return *this;
532 }
533 
534 
535 
536 template <int rank_, int dim, typename Number>
537 inline
540 {
541  for (unsigned int i=0; i<dim; ++i)
542  subtensor[i] /= s;
543  return *this;
544 }
545 
546 
547 
548 template <int rank_, int dim, typename Number>
549 inline
552 {
553  Tensor<rank_,dim,Number> tmp(*this);
554 
555  for (unsigned int i=0; i<dim; ++i)
556  tmp.subtensor[i] += t.subtensor[i];
557 
558  return tmp;
559 }
560 
561 
562 
563 template <int rank_, int dim, typename Number>
564 inline
567 {
568  Tensor<rank_,dim,Number> tmp(*this);
569 
570  for (unsigned int i=0; i<dim; ++i)
571  tmp.subtensor[i] -= t.subtensor[i];
572 
573  return tmp;
574 }
575 
576 
577 
578 template <int rank_, int dim, typename Number>
579 inline
582 {
584 
585  for (unsigned int i=0; i<dim; ++i)
586  tmp.subtensor[i] = -subtensor[i];
587 
588  return tmp;
589 }
590 
591 
592 
593 template <int rank_, int dim, typename Number>
594 inline
597 {
598  return std::sqrt (norm_square());
599 }
600 
601 
602 
603 template <int rank_, int dim, typename Number>
604 inline
607 {
608  real_type s = 0;
609  for (unsigned int i=0; i<dim; ++i)
610  s += subtensor[i].norm_square();
611 
612  return s;
613 }
614 
615 
616 
617 template <int rank_, int dim, typename Number>
618 template <typename Number2>
619 inline
620 void
622 {
623  AssertDimension (result.size(),(Utilities::fixed_power<rank_, unsigned int>(dim)));
624 
625  unsigned int index = 0;
626  unroll_recursion (result, index);
627 }
628 
629 
630 
631 template <int rank_, int dim, typename Number>
632 template <typename Number2>
633 inline
634 void
636  unsigned int &index) const
637 {
638  for (unsigned int i=0; i<dim; ++i)
639  {
640  operator[](i).unroll_recursion(result, index);
641  }
642 }
643 
644 template <int rank_, int dim, typename Number>
645 inline
646 unsigned int
648 {
649  TableIndices<rank_-1> indices1;
650  for (unsigned int i = 0; i < rank_-1; i++)
651  indices1[i] = indices[i];
652 
653  Assert (indices[rank_-1] < dim,
654  ExcIndexRange (indices[rank_-1], 0, dim));
655  return ( Tensor<rank_-1,dim,Number>::component_to_unrolled_index(indices1) * dim + indices[rank_-1]);
656 }
657 
658 template <int rank_, int dim, typename Number>
659 inline
662 {
665 
666  TableIndices<rank_> indices;
667 
668  unsigned int remainder = i;
669  for (int r=rank_-1; r>=0; --r)
670  {
671  indices[r] = (remainder % dim);
672  remainder /= dim;
673  }
674  Assert (remainder == 0, ExcInternalError());
675 
676  return indices;
677 }
678 
679 
680 
681 template <int rank_, int dim, typename Number>
682 inline
684 {
685  for (unsigned int i=0; i<dim; ++i)
686  subtensor[i].clear();
687 }
688 
689 
690 
691 template <int rank_, int dim, typename Number>
692 inline
693 std::size_t
695 {
696  return sizeof(Tensor<rank_,dim,Number>);
697 }
698 
699 
700 
701 template <int rank_, int dim, typename Number>
702 template <class Archive>
703 inline
704 void
705 Tensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
706 {
707  ar &subtensor;
708 }
709 
710 #endif // DOXYGEN
711 /* ----------------- Non-member functions operating on tensors. ------------ */
712 
713 
714 
722 template <int rank_, int dim, typename Number>
723 inline
724 std::ostream &operator << (std::ostream &out, const Tensor<rank_,dim,Number> &p)
725 {
726  for (unsigned int i=0; i<dim-1; ++i)
727  out << p[i] << ' ';
728  out << p[dim-1];
729 
730  return out;
731 }
732 
733 #ifndef DOXYGEN
734 
738 template <int rank_>
739 inline
740 std::ostream &operator << (std::ostream &out, const Tensor<rank_,1> &p)
741 {
742  out << p[0];
743 
744  return out;
745 }
746 
747 #endif // DOXYGEN
748 
749 
757 template <int dim, typename Number>
758 inline
759 Number contract (const Tensor<1,dim,Number> &src1,
760  const Tensor<1,dim,Number> &src2)
761 {
762  Number res = 0.;
763  for (unsigned int i=0; i<dim; ++i)
764  res += src1[i] * src2[i];
765 
766  return res;
767 }
768 
769 
786 template <int dim, typename Number>
787 inline
788 Number
790  const Tensor<1,dim,Number> &src2)
791 {
792  return contract(src1, src2);
793 }
794 
795 
803 template <int dim, typename Number>
804 inline
806  const Tensor<2, dim, Number> &src2)
807 {
808  Number res = 0.;
809  for (unsigned int i=0; i<dim; ++i)
810  res += contract(src1[i],src2[i]);
811 
812  return res;
813 }
814 
815 
823 template <int dim, typename Number>
824 inline
826  const Tensor<2,dim,Number> &src1,
827  const Tensor<1,dim,Number> &src2)
828 {
829  for (unsigned int i=0; i<dim; ++i)
830  {
831  dest[i] = src1[i][0] * src2[0];
832  for (unsigned int j=1; j<dim; ++j)
833  dest[i] += src1[i][j] * src2[j];
834  }
835 }
836 
837 
854 template <int dim, typename Number>
857  const Tensor<1,dim,Number> &src2)
858 {
859  Tensor<1,dim,Number> dest (false);
860  for (unsigned int i=0; i<dim; ++i)
861  {
862  dest[i] = src1[i][0] * src2[0];
863  for (unsigned int j=1; j<dim; ++j)
864  dest[i] += src1[i][j] * src2[j];
865  }
866  return dest;
867 }
868 
869 
877 template <int dim, typename Number>
878 inline
880  const Tensor<1,dim,Number> &src1,
881  const Tensor<2,dim,Number> &src2)
882 {
883  for (unsigned int i=0; i<dim; ++i)
884  {
885  dest[i] = src1[0] * src2[0][i];
886  for (unsigned int j=1; j<dim; ++j)
887  dest[i] += src1[j] * src2[j][i];
888  }
889 }
890 
891 
908 template <int dim, typename Number>
909 inline
912  const Tensor<2,dim,Number> &src2)
913 {
914  Tensor<1,dim,Number> dest (false);
915  for (unsigned int i=0; i<dim; ++i)
916  {
917  dest[i] = src1[0] * src2[0][i];
918  for (unsigned int j=1; j<dim; ++j)
919  dest[i] += src1[j] * src2[j][i];
920  }
921  return dest;
922 }
923 
924 
932 template <int dim, typename Number>
933 inline
935  const Tensor<2,dim,Number> &src1,
936  const Tensor<2,dim,Number> &src2)
937 {
938  for (unsigned int i=0; i<dim; ++i)
939  for (unsigned int j=0; j<dim; ++j)
940  {
941  dest[i][j] = src1[i][0] * src2[0][j];
942  for (unsigned int k=1; k<dim; ++k)
943  dest[i][j] += src1[i][k] * src2[k][j];
944  }
945 }
946 
947 
948 
965 template <int dim, typename Number>
966 inline
969  const Tensor<2,dim,Number> &src2)
970 {
972  for (unsigned int i=0; i<dim; ++i)
973  for (unsigned int j=0; j<dim; ++j)
974  for (unsigned int k=0; k<dim; ++k)
975  dest[i][j] += src1[i][k] * src2[k][j];
976  return dest;
977 }
978 
979 
994 template <int dim, typename Number>
995 inline
997  const Tensor<2,dim,Number> &src1, const unsigned int index1,
998  const Tensor<2,dim,Number> &src2, const unsigned int index2)
999 {
1000  dest.clear ();
1001 
1002  switch (index1)
1003  {
1004  case 1:
1005  switch (index2)
1006  {
1007  case 1:
1008  for (unsigned int i=0; i<dim; ++i)
1009  for (unsigned int j=0; j<dim; ++j)
1010  for (unsigned int k=0; k<dim; ++k)
1011  dest[i][j] += src1[k][i] * src2[k][j];
1012  break;
1013  case 2:
1014  for (unsigned int i=0; i<dim; ++i)
1015  for (unsigned int j=0; j<dim; ++j)
1016  for (unsigned int k=0; k<dim; ++k)
1017  dest[i][j] += src1[k][i] * src2[j][k];
1018  break;
1019 
1020  default:
1021  Assert (false,
1022  (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
1023  };
1024  break;
1025  case 2:
1026  switch (index2)
1027  {
1028  case 1:
1029  for (unsigned int i=0; i<dim; ++i)
1030  for (unsigned int j=0; j<dim; ++j)
1031  for (unsigned int k=0; k<dim; ++k)
1032  dest[i][j] += src1[i][k] * src2[k][j];
1033  break;
1034  case 2:
1035  for (unsigned int i=0; i<dim; ++i)
1036  for (unsigned int j=0; j<dim; ++j)
1037  for (unsigned int k=0; k<dim; ++k)
1038  dest[i][j] += src1[i][k] * src2[j][k];
1039  break;
1040 
1041  default:
1042  Assert (false,
1043  (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
1044  };
1045  break;
1046 
1047  default:
1048  Assert (false, (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index1)));
1049  };
1050 }
1051 
1052 
1064 template <int dim, typename Number>
1065 inline
1067  const Tensor<3,dim,Number> &src1, const unsigned int index1,
1068  const Tensor<1,dim,Number> &src2)
1069 {
1070  dest.clear ();
1071 
1072  switch (index1)
1073  {
1074  case 1:
1075  for (unsigned int i=0; i<dim; ++i)
1076  for (unsigned int j=0; j<dim; ++j)
1077  for (unsigned int k=0; k<dim; ++k)
1078  dest[i][j] += src1[k][i][j] * src2[k];
1079  break;
1080 
1081  case 2:
1082  for (unsigned int i=0; i<dim; ++i)
1083  for (unsigned int j=0; j<dim; ++j)
1084  for (unsigned int k=0; k<dim; ++k)
1085  dest[i][j] += src1[i][k][j] * src2[k];
1086  break;
1087 
1088  case 3:
1089  for (unsigned int i=0; i<dim; ++i)
1090  for (unsigned int j=0; j<dim; ++j)
1091  for (unsigned int k=0; k<dim; ++k)
1092  dest[i][j] += src1[i][j][k] * src2[k];
1093  break;
1094 
1095  default:
1096  Assert (false,
1097  (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index1)));
1098  };
1099 }
1100 
1101 
1109 template <int dim, typename Number>
1110 inline
1112  const Tensor<3,dim,Number> &src1,
1113  const Tensor<2,dim,Number> &src2)
1114 {
1115  dest.clear ();
1116  for (unsigned int i=0; i<dim; ++i)
1117  for (unsigned int j=0; j<dim; ++j)
1118  for (unsigned int k=0; k<dim; ++k)
1119  for (unsigned int l=0; l<dim; ++l)
1120  dest[i][j][k] += src1[i][j][l] * src2[l][k];
1121 }
1122 
1123 
1137 template <int dim, typename Number>
1138 inline
1140  const Tensor<3,dim,Number> &src1, const unsigned int index1,
1141  const Tensor<2,dim,Number> &src2, const unsigned int index2)
1142 {
1143  dest.clear ();
1144 
1145  switch (index1)
1146  {
1147  case 1:
1148  switch (index2)
1149  {
1150  case 1:
1151  for (unsigned int i=0; i<dim; ++i)
1152  for (unsigned int j=0; j<dim; ++j)
1153  for (unsigned int k=0; k<dim; ++k)
1154  for (unsigned int l=0; l<dim; ++l)
1155  dest[i][j][k] += src1[l][i][j] * src2[l][k];
1156  break;
1157  case 2:
1158  for (unsigned int i=0; i<dim; ++i)
1159  for (unsigned int j=0; j<dim; ++j)
1160  for (unsigned int k=0; k<dim; ++k)
1161  for (unsigned int l=0; l<dim; ++l)
1162  dest[i][j][k] += src1[l][i][j] * src2[k][l];
1163  break;
1164  default:
1165  Assert (false,
1166  (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
1167  }
1168 
1169  break;
1170  case 2:
1171  switch (index2)
1172  {
1173  case 1:
1174  for (unsigned int i=0; i<dim; ++i)
1175  for (unsigned int j=0; j<dim; ++j)
1176  for (unsigned int k=0; k<dim; ++k)
1177  for (unsigned int l=0; l<dim; ++l)
1178  dest[i][j][k] += src1[i][l][j] * src2[l][k];
1179  break;
1180  case 2:
1181  for (unsigned int i=0; i<dim; ++i)
1182  for (unsigned int j=0; j<dim; ++j)
1183  for (unsigned int k=0; k<dim; ++k)
1184  for (unsigned int l=0; l<dim; ++l)
1185  dest[i][j][k] += src1[i][l][j] * src2[k][l];
1186  break;
1187  default:
1188  Assert (false,
1189  (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
1190  }
1191 
1192  break;
1193  case 3:
1194  switch (index2)
1195  {
1196  case 1:
1197  for (unsigned int i=0; i<dim; ++i)
1198  for (unsigned int j=0; j<dim; ++j)
1199  for (unsigned int k=0; k<dim; ++k)
1200  for (unsigned int l=0; l<dim; ++l)
1201  dest[i][j][k] += src1[i][j][l] * src2[l][k];
1202  break;
1203  case 2:
1204  for (unsigned int i=0; i<dim; ++i)
1205  for (unsigned int j=0; j<dim; ++j)
1206  for (unsigned int k=0; k<dim; ++k)
1207  for (unsigned int l=0; l<dim; ++l)
1208  dest[i][j][k] += src1[i][j][l] * src2[k][l];
1209  break;
1210  default:
1211  Assert (false,
1212  (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
1213  }
1214 
1215  break;
1216  default:
1217  Assert (false,
1218  (typename Tensor<3,dim,Number>::ExcInvalidTensorIndex (index1)));
1219  }
1220 }
1221 
1222 
1239 template <int dim, typename Number>
1240 inline
1243  const Tensor<2,dim,Number> &src2)
1244 {
1245  Tensor<3,dim,Number> dest;
1246  for (unsigned int i=0; i<dim; ++i)
1247  for (unsigned int j=0; j<dim; ++j)
1248  for (unsigned int k=0; k<dim; ++k)
1249  for (unsigned int l=0; l<dim; ++l)
1250  dest[i][j][k] += src1[i][j][l] * src2[l][k];
1251  return dest;
1252 }
1253 
1254 
1262 template <int dim, typename Number>
1263 inline
1265  const Tensor<2,dim,Number> &src1,
1266  const Tensor<3,dim,Number> &src2)
1267 {
1268  dest.clear ();
1269  for (unsigned int i=0; i<dim; ++i)
1270  for (unsigned int j=0; j<dim; ++j)
1271  for (unsigned int k=0; k<dim; ++k)
1272  for (unsigned int l=0; l<dim; ++l)
1273  dest[i][j][k] += src1[i][l] * src2[l][j][k];
1274 }
1275 
1276 
1293 template <int dim, typename Number>
1294 inline
1297  const Tensor<3,dim,Number> &src2)
1298 {
1299  Tensor<3,dim,Number> dest;
1300  for (unsigned int i=0; i<dim; ++i)
1301  for (unsigned int j=0; j<dim; ++j)
1302  for (unsigned int k=0; k<dim; ++k)
1303  for (unsigned int l=0; l<dim; ++l)
1304  dest[i][j][k] += src1[i][l] * src2[l][j][k];
1305  return dest;
1306 }
1307 
1308 
1316 template <int dim, typename Number>
1317 inline
1320  const Tensor<3,dim,Number> &src2)
1321 {
1322  Tensor<4,dim,Number> dest;
1323  for (unsigned int i=0; i<dim; ++i)
1324  for (unsigned int j=0; j<dim; ++j)
1325  for (unsigned int k=0; k<dim; ++k)
1326  for (unsigned int l=0; l<dim; ++l)
1327  for (unsigned int m=0; m<dim; ++m)
1328  dest[i][j][k][l] += src1[i][j][m] * src2[m][k][l];
1329  return dest;
1330 }
1331 
1332 
1341 template <int dim, typename Number>
1342 inline
1344  const Tensor<4,dim,Number> &src1,
1345  const Tensor<2,dim,Number> &src2)
1346 {
1347  dest.clear ();
1348  for (unsigned int i=0; i<dim; ++i)
1349  for (unsigned int j=0; j<dim; ++j)
1350  for (unsigned int k=0; k<dim; ++k)
1351  for (unsigned int l=0; l<dim; ++l)
1352  dest[i][j] += src1[i][j][k][l] * src2[k][l];
1353 }
1354 
1355 
1363 template <int dim, typename Number>
1364 inline
1366  const Tensor<2,dim,Number> &A,
1367  const Tensor<1,dim,Number> &v)
1368 {
1369  Number result = 0.;
1370 
1371  for (unsigned int i=0; i<dim; ++i)
1372  for (unsigned int j=0; j<dim; ++j)
1373  result += u[i] * A[i][j] * v[j];
1374  return result;
1375 }
1376 
1377 
1385 template <int dim, typename Number>
1386 inline
1387 Number
1389  const Tensor<3,dim,Number> &t2,
1390  const Tensor<2,dim,Number> &t3)
1391 {
1392  Number s = 0;
1393  for (unsigned int i=0; i<dim; ++i)
1394  for (unsigned int j=0; j<dim; ++j)
1395  for (unsigned int k=0; k<dim; ++k)
1396  s += t1[i] * t2[i][j][k] * t3[j][k];
1397  return s;
1398 }
1399 
1400 
1408 template <int dim, typename Number>
1409 inline
1410 Number
1412  const Tensor<3,dim,Number> &t2,
1413  const Tensor<1,dim,Number> &t3)
1414 {
1415  Number s = 0;
1416  for (unsigned int i=0; i<dim; ++i)
1417  for (unsigned int j=0; j<dim; ++j)
1418  for (unsigned int k=0; k<dim; ++k)
1419  s += t1[i][j] * t2[i][j][k] * t3[k];
1420  return s;
1421 }
1422 
1423 
1431 template <int dim, typename Number>
1432 inline
1433 Number
1435  const Tensor<4,dim,Number> &t2,
1436  const Tensor<2,dim,Number> &t3)
1437 {
1438  Number s = 0;
1439  for (unsigned int i=0; i<dim; ++i)
1440  for (unsigned int j=0; j<dim; ++j)
1441  for (unsigned int k=0; k<dim; ++k)
1442  for (unsigned int l=0; l<dim; ++l)
1443  s += t1[i][j] * t2[i][j][k][l] * t3[k][l];
1444  return s;
1445 }
1446 
1447 
1455 template <int dim, typename Number>
1457  const Tensor<1,dim,Number> &src1,
1458  const Tensor<1,dim,Number> &src2)
1459 {
1460  for (unsigned int i=0; i<dim; ++i)
1461  for (unsigned int j=0; j<dim; ++j)
1462  dst[i][j] = src1[i] * src2[j];
1463 }
1464 
1465 
1473 template <int dim, typename Number>
1475  const Tensor<1,dim,Number> &src1,
1476  const Tensor<2,dim,Number> &src2)
1477 {
1478  for (unsigned int i=0; i<dim; ++i)
1479  for (unsigned int j=0; j<dim; ++j)
1480  for (unsigned int k=0; k<dim; ++k)
1481  dst[i][j][k] = src1[i] * src2[j][k];
1482 }
1483 
1484 
1492 template <int dim, typename Number>
1494  const Tensor<2,dim,Number> &src1,
1495  const Tensor<1,dim,Number> &src2)
1496 {
1497  for (unsigned int i=0; i<dim; ++i)
1498  for (unsigned int j=0; j<dim; ++j)
1499  for (unsigned int k=0; k<dim; ++k)
1500  dst[i][j][k] = src1[i][j] * src2[k];
1501 }
1502 
1503 
1515 template <int dim, typename Number>
1517  const Number src1,
1518  const Tensor<1,dim,Number> &src2)
1519 {
1520  for (unsigned int i=0; i<dim; ++i)
1521  dst[i] = src1 * src2[i];
1522 }
1523 
1524 
1525 
1537 template <int dim, typename Number>
1539  const Tensor<1,dim,Number> src1,
1540  const Number src2)
1541 {
1542  for (unsigned int i=0; i<dim; ++i)
1543  dst[i] = src1[i] * src2;
1544 }
1545 
1546 
1559 template <int dim, typename Number>
1560 inline
1561 void
1563  const Tensor<1,dim,Number> &src)
1564 {
1565  Assert (dim==2, ExcInternalError());
1566 
1567  dst[0] = src[1];
1568  dst[1] = -src[0];
1569 }
1570 
1571 
1582 template <int dim, typename Number>
1583 inline
1584 void
1586  const Tensor<1,dim,Number> &src1,
1587  const Tensor<1,dim,Number> &src2)
1588 {
1589  Assert (dim==3, ExcInternalError());
1590 
1591  dst[0] = src1[1]*src2[2] - src1[2]*src2[1];
1592  dst[1] = src1[2]*src2[0] - src1[0]*src2[2];
1593  dst[2] = src1[0]*src2[1] - src1[1]*src2[0];
1594 }
1595 
1596 
1607 template <int dim, typename Number>
1608 inline
1609 Number
1611  const Tensor<2,dim,Number> &t2)
1612 {
1613  Number s = 0;
1614  for (unsigned int i=0; i<dim; ++i)
1615  for (unsigned int j=0; j<dim; ++j)
1616  s += t1[i][j] * t2[i][j];
1617  return s;
1618 }
1619 
1620 
1629 template <int rank, typename Number>
1630 inline
1632 {
1633  return determinant(t[0]);
1634 }
1635 
1636 
1637 
1646 template <typename Number>
1647 inline
1649 {
1650  return t[0];
1651 }
1652 
1653 
1662 template <typename Number>
1663 inline
1665 {
1666  return t[0][0];
1667 }
1668 
1669 
1670 
1677 template <typename Number>
1678 inline
1680 {
1681  return ((t[0][0] * t[1][1]) - (t[1][0] * t[0][1]));
1682 }
1683 
1684 
1691 template <typename Number>
1692 inline
1694 {
1695  // use exactly the same expression with the
1696  // same order of operations as for the inverse
1697  // to let the compiler use common
1698  // subexpression elimination when using
1699  // determinant and inverse in nearby code
1700  const Number t4 = t[0][0]*t[1][1],
1701  t6 = t[0][0]*t[1][2],
1702  t8 = t[0][1]*t[1][0],
1703  t00 = t[0][2]*t[1][0],
1704  t01 = t[0][1]*t[2][0],
1705  t04 = t[0][2]*t[2][0],
1706  det = (t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
1707  t00*t[2][1]+t01*t[1][2]-t04*t[1][1]);
1708  return det;
1709 }
1710 
1711 
1719 template <int dim, typename Number>
1720 inline
1722 {
1723  // compute the determinant using the
1724  // Laplace expansion of the
1725  // determinant. this may not be the most
1726  // efficient algorithm, but it does for
1727  // small n.
1728  //
1729  // for some algorithmic simplicity, we use
1730  // the expansion along the last row
1731  Number det = 0;
1732 
1733  for (unsigned int k=0; k<dim; ++k)
1734  {
1735  Tensor<2,dim-1,Number> minor;
1736  for (unsigned int i=0; i<dim-1; ++i)
1737  for (unsigned int j=0; j<dim-1; ++j)
1738  minor[i][j] = t[i][j<k ? j : j+1];
1739 
1740  const Number cofactor = std::pow (-1., static_cast<Number>(k+1)) *
1741  determinant (minor);
1742 
1743  det += t[dim-1][k] * cofactor;
1744  }
1745 
1746  return std::pow (-1., static_cast<Number>(dim)) * det;
1747 }
1748 
1749 
1750 
1758 template <int dim, typename Number>
1759 Number trace (const Tensor<2,dim,Number> &d)
1760 {
1761  Number t=0;
1762  for (unsigned int i=0; i<dim; ++i)
1763  t += d[i][i];
1764  return t;
1765 }
1766 
1767 
1768 
1778 template <int dim, typename Number>
1779 inline
1782 {
1783  Number return_tensor [dim][dim];
1784  switch (dim)
1785  {
1786  case 1:
1787  return_tensor[0][0] = 1.0/t[0][0];
1788  break;
1789 
1790  case 2:
1791  // this is Maple output,
1792  // thus a bit unstructured
1793  {
1794  const Number det = t[0][0]*t[1][1]-t[1][0]*t[0][1];
1795  const Number t4 = 1.0/det;
1796  return_tensor[0][0] = t[1][1]*t4;
1797  return_tensor[0][1] = -t[0][1]*t4;
1798  return_tensor[1][0] = -t[1][0]*t4;
1799  return_tensor[1][1] = t[0][0]*t4;
1800  break;
1801  }
1802 
1803  case 3:
1804  {
1805  const Number t4 = t[0][0]*t[1][1],
1806  t6 = t[0][0]*t[1][2],
1807  t8 = t[0][1]*t[1][0],
1808  t00 = t[0][2]*t[1][0],
1809  t01 = t[0][1]*t[2][0],
1810  t04 = t[0][2]*t[2][0],
1811  det = (t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
1812  t00*t[2][1]+t01*t[1][2]-t04*t[1][1]),
1813  t07 = 1.0/det;
1814  return_tensor[0][0] = (t[1][1]*t[2][2]-t[1][2]*t[2][1])*t07;
1815  return_tensor[0][1] = (t[0][2]*t[2][1]-t[0][1]*t[2][2])*t07;
1816  return_tensor[0][2] = (t[0][1]*t[1][2]-t[0][2]*t[1][1])*t07;
1817  return_tensor[1][0] = (t[1][2]*t[2][0]-t[1][0]*t[2][2])*t07;
1818  return_tensor[1][1] = (t[0][0]*t[2][2]-t04)*t07;
1819  return_tensor[1][2] = (t00-t6)*t07;
1820  return_tensor[2][0] = (t[1][0]*t[2][1]-t[1][1]*t[2][0])*t07;
1821  return_tensor[2][1] = (t01-t[0][0]*t[2][1])*t07;
1822  return_tensor[2][2] = (t4-t8)*t07;
1823 
1824  break;
1825  }
1826 
1827  // if desired, take over the
1828  // inversion of a 4x4 tensor
1829  // from the FullMatrix
1830  default:
1831  AssertThrow (false, ExcNotImplemented());
1832  }
1833  return Tensor<2,dim,Number>(return_tensor);
1834 }
1835 
1836 
1837 
1848 template <int dim, typename Number>
1849 inline
1852 {
1853  Number tt[dim][dim];
1854  for (unsigned int i=0; i<dim; ++i)
1855  {
1856  tt[i][i] = t[i][i];
1857  for (unsigned int j=i+1; j<dim; ++j)
1858  {
1859  tt[i][j] = t[j][i];
1860  tt[j][i] = t[i][j];
1861  };
1862  }
1863  return Tensor<2,dim,Number>(tt);
1864 }
1865 
1866 #ifndef DOXYGEN
1867 
1875 template <typename Number>
1876 inline
1878 transpose (const Tensor<2,1,Number> &t)
1879 {
1880  return t;
1881 }
1882 
1883 
1884 
1885 
1893 template <typename Number>
1894 inline
1896 transpose (const Tensor<2,2,Number> &t)
1897 {
1898  const Number x[2][2] = {{t[0][0], t[1][0]}, {t[0][1], t[1][1]}};
1899  return Tensor<2,2,Number>(x);
1900 }
1901 
1902 
1903 
1904 
1912 template <typename Number>
1913 inline
1915 transpose (const Tensor<2,3,Number> &t)
1916 {
1917  const Number x[3][3] = {{t[0][0], t[1][0], t[2][0]},
1918  {t[0][1], t[1][1], t[2][1]},
1919  {t[0][2], t[1][2], t[2][2]}
1920  };
1921  return Tensor<2,3,Number>(x);
1922 }
1923 
1924 #endif // DOXYGEN
1925 
1926 
1935 template <int dim, typename Number>
1936 inline
1937 double
1939 {
1940  double max = 0;
1941  for (unsigned int j=0; j<dim; ++j)
1942  {
1943  double sum = 0;
1944  for (unsigned int i=0; i<dim; ++i)
1945  sum += std::fabs(t[i][j]);
1946 
1947  if (sum > max)
1948  max = sum;
1949  }
1950 
1951  return max;
1952 }
1953 
1954 
1955 
1964 template <int dim, typename Number>
1965 inline
1966 double
1968 {
1969  double max = 0;
1970  for (unsigned int i=0; i<dim; ++i)
1971  {
1972  double sum = 0;
1973  for (unsigned int j=0; j<dim; ++j)
1974  sum += std::fabs(t[i][j]);
1975 
1976  if (sum > max)
1977  max = sum;
1978  }
1979 
1980  return max;
1981 }
1982 
1983 
1984 
1991 template <int rank, int dim, typename Number>
1992 inline
1995  const Number factor)
1996 {
1997  Tensor<rank,dim,Number> tt = t;
1998  tt *= factor;
1999  return tt;
2000 }
2001 
2002 
2003 
2010 template <int rank, int dim, typename Number>
2011 inline
2013 operator * (const Number factor,
2014  const Tensor<rank,dim,Number> &t)
2015 {
2016  Tensor<rank,dim,Number> tt = t;
2017  tt *= factor;
2018  return tt;
2019 }
2020 
2021 
2022 
2028 template <int rank, int dim, typename Number>
2029 inline
2032  const Number factor)
2033 {
2034  Tensor<rank,dim,Number> tt = t;
2035  tt /= factor;
2036  return tt;
2037 }
2038 
2039 
2040 
2041 
2048 template <int rank, int dim>
2049 inline
2052  const double factor)
2053 {
2054  Tensor<rank,dim> tt = t;
2055  tt *= factor;
2056  return tt;
2057 }
2058 
2059 
2060 
2067 template <int rank, int dim>
2068 inline
2070 operator * (const double factor,
2071  const Tensor<rank,dim> &t)
2072 {
2073  Tensor<rank,dim> tt = t;
2074  tt *= factor;
2075  return tt;
2076 }
2077 
2078 
2079 
2085 template <int rank, int dim>
2086 inline
2089  const double factor)
2090 {
2091  Tensor<rank,dim> tt = t;
2092  tt /= factor;
2093  return tt;
2094 }
2095 
2096 DEAL_II_NAMESPACE_CLOSE
2097 
2098 #endif
Number determinant(const Tensor< 1, 1, Number > &t)
Definition: tensor.h:1648
void contract(Tensor< 1, dim, Number > &dest, const Tensor< 2, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:825
Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:1759
Number scalar_product(const Tensor< 2, dim, Number > &t1, const Tensor< 2, dim, Number > &t2)
Definition: tensor.h:1610
void contract(Tensor< 2, dim, Number > &dest, const Tensor< 2, dim, Number > &src1, const Tensor< 2, dim, Number > &src2)
Definition: tensor.h:934
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
Tensor< rank_, dim, Number > operator-() const
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
static const unsigned int n_independent_components
Definition: tensor.h:85
Number contract3(const Tensor< 2, dim, Number > &t1, const Tensor< 4, dim, Number > &t2, const Tensor< 2, dim, Number > &t3)
Definition: tensor.h:1434
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
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
Tensor< 1, dim, Number > operator/(const Tensor< 1, dim, Number > &t, const Number factor)
Definition: tensor_base.h:1399
Tensor< rank_, dim, Number > & operator-=(const Tensor< rank_, dim, Number > &)
Tensor< rank_-1, dim, Number > & operator[](const unsigned int i)
Number determinant(const Tensor< 2, 2, Number > &t)
Definition: tensor.h:1679
void double_contract(Tensor< 2, dim, Number > &dest, const Tensor< 4, dim, Number > &src1, const Tensor< 2, dim, Number > &src2)
Definition: tensor.h:1343
void cross_product(Tensor< 1, dim, Number > &dst, const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:1585
void cross_product(Tensor< 1, dim, Number > &dst, const Tensor< 1, dim, Number > &src)
Definition: tensor.h:1562
bool operator!=(const Tensor< rank_, dim, Number > &) const
void contract(Tensor< 3, dim, Number > &dest, const Tensor< 3, dim, Number > &src1, const Tensor< 2, dim, Number > &src2)
Definition: tensor.h:1111
Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1851
Tensor< rank_-1, dim, Number > subtensor[dim]
Definition: tensor.h:333
void outer_product(Tensor< 3, dim, Number > &dst, const Tensor< 2, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:1493
numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:107
#define Assert(cond, exc)
Definition: exceptions.h:299
Tensor< rank_-1, dim, Number >::array_type array_type[dim]
Definition: tensor.h:115
void outer_product(Tensor< 1, dim, Number > &dst, const Tensor< 1, dim, Number > src1, const Number src2)
Definition: tensor.h:1538
Number contract3(const Tensor< 1, dim, Number > &t1, const Tensor< 3, dim, Number > &t2, const Tensor< 2, dim, Number > &t3)
Definition: tensor.h:1388
static std::size_t memory_consumption()
void serialize(Archive &ar, const unsigned int version)
double linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1967
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
void contract(Tensor< 3, dim, Number > &dest, const Tensor< 3, dim, Number > &src1, const unsigned int index1, const Tensor< 2, dim, Number > &src2, const unsigned int index2)
Definition: tensor.h:1139
double l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1938
DeclException1(ExcInvalidTensorIndex, int,<< "Invalid tensor index "<< arg1)
Number contract3(const Tensor< 2, dim, Number > &t1, const Tensor< 3, dim, Number > &t2, const Tensor< 1, dim, Number > &t3)
Definition: tensor.h:1411
bool operator==(const Tensor< rank_, dim, Number > &) const
void contract(Tensor< 3, dim, Number > &dest, const Tensor< 2, dim, Number > &src1, const Tensor< 3, dim, Number > &src2)
Definition: tensor.h:1264
void contract(Tensor< 2, dim, Number > &dest, const Tensor< 3, dim, Number > &src1, const unsigned int index1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:1066
Number contract3(const Tensor< 1, dim, Number > &u, const Tensor< 2, dim, Number > &A, const Tensor< 1, dim, Number > &v)
Definition: tensor.h:1365
Number double_contract(const Tensor< 2, dim, Number > &src1, const Tensor< 2, dim, Number > &src2)
Definition: tensor.h:805
Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1721
static const unsigned int dimension
Definition: tensor.h:71
void contract(Tensor< 1, dim, Number > &dest, const Tensor< 1, dim, Number > &src1, const Tensor< 2, dim, Number > &src2)
Definition: tensor.h:879
void contract(Tensor< 2, dim, Number > &dest, const Tensor< 2, dim, Number > &src1, const unsigned int index1, const Tensor< 2, dim, Number > &src2, const unsigned int index2)
Definition: tensor.h:996
Tensor< rank_, dim, Number > & operator/=(const Number factor)
Number determinant(const Tensor< rank, 1, Number > &t)
Definition: tensor.h:1631
void outer_product(Tensor< 2, dim, Number > &dst, const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:1456
DerivativeForm< 1, spacedim, dim > transpose(const DerivativeForm< 1, dim, spacedim > &DF)
::ExceptionBase & ExcNotImplemented()
void outer_product(Tensor< 1, dim, Number > &dst, const Number src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:1516
void outer_product(Tensor< 3, dim, Number > &dst, const Tensor< 1, dim, Number > &src1, const Tensor< 2, dim, Number > &src2)
Definition: tensor.h:1474
Number contract(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:759
Tensor< rank_-1, dim, Number > value_type
Definition: tensor.h:91
Tensor< rank_, dim, Number > & operator+=(const Tensor< rank_, dim, Number > &)
::ExceptionBase & ExcInternalError()
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
void clear()
Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1781
Tensor< rank_, dim, Number > & operator*=(const Number factor)
void unroll(Vector< Number2 > &result) const
VectorizedArray< Number > max(const ::VectorizedArray< Number > &x, const ::VectorizedArray< Number > &y)
real_type norm() const
Tensor & operator=(const Tensor< rank_, dim, Number > &)
Number determinant(const Tensor< 2, 3, Number > &t)
Definition: tensor.h:1693
Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:1664