00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__tensor_h
00013 #define __deal2__tensor_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/tensor_base.h>
00018 #include <deal.II/base/utilities.h>
00019
00020 DEAL_II_NAMESPACE_OPEN
00021 template <int rank_, int dim, typename Number> class Tensor;
00022 template <int dim, typename Number> class Tensor<1,dim,Number>;
00023
00045 template <int rank_, int dim, typename Number>
00046 class Tensor
00047 {
00048 public:
00066 static const unsigned int dimension = dim;
00067
00072 static const unsigned int rank = rank_;
00073
00078 typedef Tensor<rank_-1,dim,Number> value_type;
00079
00094 typedef typename numbers::NumberTraits<Number>::real_type real_type;
00095
00102 typedef typename Tensor<rank_-1,dim,Number>::array_type array_type[dim];
00103
00108 Tensor ();
00109
00114 Tensor (const array_type &initializer);
00115
00120 Tensor (const Tensor<1,dim,Tensor<rank_-1,dim,Number> > &tensor_in);
00121
00126 operator Tensor<1,dim,Tensor<rank_-1,dim,Number> > () const;
00127
00131 Tensor<rank_-1,dim,Number> &operator [] (const unsigned int i);
00132
00136 const Tensor<rank_-1,dim,Number> &operator [] (const unsigned int i) const;
00137
00141 Tensor & operator = (const Tensor<rank_,dim,Number> &);
00142
00155 Tensor<rank_,dim,Number> & operator = (const Number d);
00156
00160 bool operator == (const Tensor<rank_,dim,Number> &) const;
00161
00165 bool operator != (const Tensor<rank_,dim,Number> &) const;
00166
00170 Tensor<rank_,dim,Number> & operator += (const Tensor<rank_,dim,Number> &);
00171
00175 Tensor<rank_,dim,Number> & operator -= (const Tensor<rank_,dim,Number> &);
00176
00182 Tensor<rank_,dim,Number> & operator *= (const Number factor);
00183
00188 Tensor<rank_,dim,Number> & operator /= (const Number factor);
00189
00196 Tensor<rank_,dim,Number> operator + (const Tensor<rank_,dim,Number> &) const;
00197
00204 Tensor<rank_,dim,Number> operator - (const Tensor<rank_,dim,Number> &) const;
00205
00210 Tensor<rank_,dim,Number> operator - () const;
00211
00217 real_type norm () const;
00218
00231 real_type norm_square () const;
00232
00242 template <typename Number2>
00243 void unroll (Vector<Number2> & result) const;
00244
00245
00264 void clear ();
00265
00271 static std::size_t memory_consumption ();
00272
00276 DeclException1 (ExcInvalidTensorIndex,
00277 int,
00278 << "Invalid tensor index " << arg1);
00279
00284 template <class Archive>
00285 void serialize(Archive & ar, const unsigned int version);
00286
00287 private:
00292 Tensor<rank_-1,dim,Number> subtensor[dim];
00293
00297 template <typename Number2>
00298 void unroll_recursion(Vector<Number2> &result,
00299 unsigned int &start_index) const;
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311 template <int, int, typename> friend class Tensor;
00312 };
00313
00314
00315
00316
00317 #ifndef DOXYGEN
00318
00319 template <int rank_, int dim, typename Number>
00320 inline
00321 Tensor<rank_,dim,Number>::Tensor ()
00322 {
00323
00324
00325
00326 }
00327
00328
00329
00330 template <int rank_, int dim, typename Number>
00331 inline
00332 Tensor<rank_,dim,Number>::Tensor (const array_type &initializer)
00333 {
00334 for (unsigned int i=0; i<dim; ++i)
00335 subtensor[i] = Tensor<rank_-1,dim,Number>(initializer[i]);
00336 }
00337
00338
00339
00340 template <int rank_, int dim, typename Number>
00341 inline
00342 Tensor<rank_,dim,Number>::Tensor
00343 (const Tensor<1,dim,Tensor<rank_-1,dim,Number> > &tensor_in)
00344 {
00345 for (unsigned int i=0; i<dim; ++i)
00346 subtensor[i] = tensor_in[i];
00347 }
00348
00349
00350
00351 template <int rank_, int dim, typename Number>
00352 inline
00353 Tensor<rank_,dim,Number>::operator
00354 Tensor<1,dim,Tensor<rank_-1,dim,Number> > () const
00355 {
00356 return Tensor<1,dim,Tensor<rank_-1,dim,Number> > (subtensor);
00357 }
00358
00359
00360
00361 template <int rank_, int dim, typename Number>
00362 inline
00363 typename Tensor<rank_,dim,Number>::value_type&
00364 Tensor<rank_,dim,Number>::operator[] (const unsigned int i)
00365 {
00366 Assert (i<dim, ExcIndexRange(i, 0, dim));
00367 return subtensor[i];
00368 }
00369
00370
00371
00372 template <int rank_, int dim, typename Number>
00373 inline
00374 const typename Tensor<rank_,dim,Number>::value_type&
00375 Tensor<rank_,dim,Number>::operator[] (const unsigned int i) const
00376 {
00377 Assert (i<dim, ExcIndexRange(i, 0, dim));
00378
00379 return subtensor[i];
00380 }
00381
00382
00383
00384 template <int rank_, int dim, typename Number>
00385 inline
00386 Tensor<rank_,dim,Number> &
00387 Tensor<rank_,dim,Number>::operator = (const Tensor<rank_,dim,Number> &t)
00388 {
00389 for (unsigned int i=0; i<dim; ++i)
00390 subtensor[i] = t.subtensor[i];
00391 return *this;
00392 }
00393
00394
00395
00396 template <int rank_, int dim, typename Number>
00397 inline
00398 Tensor<rank_,dim,Number> &
00399 Tensor<rank_,dim,Number>::operator = (const Number d)
00400 {
00401 Assert (d==Number(0), ExcMessage ("Only assignment with zero is allowed"));
00402
00403 for (unsigned int i=0; i<dim; ++i)
00404 subtensor[i] = 0;
00405 return *this;
00406 }
00407
00408
00409
00410 template <int rank_, int dim, typename Number>
00411 inline
00412 bool
00413 Tensor<rank_,dim,Number>::operator == (const Tensor<rank_,dim,Number> &p) const
00414 {
00415 for (unsigned int i=0; i<dim; ++i)
00416 if (subtensor[i] != p.subtensor[i]) return false;
00417 return true;
00418 }
00419
00420
00421
00422 template <int rank_, int dim, typename Number>
00423 inline
00424 bool
00425 Tensor<rank_,dim,Number>::operator != (const Tensor<rank_,dim,Number> &p) const
00426 {
00427 return !((*this) == p);
00428 }
00429
00430
00431
00432 template <int rank_, int dim, typename Number>
00433 inline
00434 Tensor<rank_,dim,Number> &
00435 Tensor<rank_,dim,Number>::operator += (const Tensor<rank_,dim,Number> &p)
00436 {
00437 for (unsigned int i=0; i<dim; ++i)
00438 subtensor[i] += p.subtensor[i];
00439 return *this;
00440 }
00441
00442
00443
00444 template <int rank_, int dim, typename Number>
00445 inline
00446 Tensor<rank_,dim,Number> &
00447 Tensor<rank_,dim,Number>::operator -= (const Tensor<rank_,dim,Number> &p)
00448 {
00449 for (unsigned int i=0; i<dim; ++i)
00450 subtensor[i] -= p.subtensor[i];
00451 return *this;
00452 }
00453
00454
00455
00456 template <int rank_, int dim, typename Number>
00457 inline
00458 Tensor<rank_,dim,Number> &
00459 Tensor<rank_,dim,Number>::operator *= (const Number s)
00460 {
00461 for (unsigned int i=0; i<dim; ++i)
00462 subtensor[i] *= s;
00463 return *this;
00464 }
00465
00466
00467
00468 template <int rank_, int dim, typename Number>
00469 inline
00470 Tensor<rank_,dim,Number> &
00471 Tensor<rank_,dim,Number>::operator /= (const Number s)
00472 {
00473 for (unsigned int i=0; i<dim; ++i)
00474 subtensor[i] /= s;
00475 return *this;
00476 }
00477
00478
00479
00480 template <int rank_, int dim, typename Number>
00481 inline
00482 Tensor<rank_,dim,Number>
00483 Tensor<rank_,dim,Number>::operator + (const Tensor<rank_,dim,Number> &t) const
00484 {
00485 Tensor<rank_,dim,Number> tmp(*this);
00486
00487 for (unsigned int i=0; i<dim; ++i)
00488 tmp.subtensor[i] += t.subtensor[i];
00489
00490 return tmp;
00491 }
00492
00493
00494
00495 template <int rank_, int dim, typename Number>
00496 inline
00497 Tensor<rank_,dim,Number>
00498 Tensor<rank_,dim,Number>::operator - (const Tensor<rank_,dim,Number> &t) const
00499 {
00500 Tensor<rank_,dim,Number> tmp(*this);
00501
00502 for (unsigned int i=0; i<dim; ++i)
00503 tmp.subtensor[i] -= t.subtensor[i];
00504
00505 return tmp;
00506 }
00507
00508
00509
00510 template <int rank_, int dim, typename Number>
00511 inline
00512 Tensor<rank_,dim,Number>
00513 Tensor<rank_,dim,Number>::operator - () const
00514 {
00515 Tensor<rank_,dim,Number> tmp;
00516
00517 for (unsigned int i=0; i<dim; ++i)
00518 tmp.subtensor[i] = -subtensor[i];
00519
00520 return tmp;
00521 }
00522
00523
00524
00525 template <int rank_, int dim, typename Number>
00526 inline
00527 typename Tensor<rank_,dim,Number>::real_type
00528 Tensor<rank_,dim,Number>::norm () const
00529 {
00530 return std::sqrt (norm_square());
00531 }
00532
00533
00534
00535 template <int rank_, int dim, typename Number>
00536 inline
00537 typename Tensor<rank_,dim,Number>::real_type
00538 Tensor<rank_,dim,Number>::norm_square () const
00539 {
00540 real_type s = 0;
00541 for (unsigned int i=0; i<dim; ++i)
00542 s += subtensor[i].norm_square();
00543
00544 return s;
00545 }
00546
00547
00548
00549 template <int rank_, int dim, typename Number>
00550 template <typename Number2>
00551 inline
00552 void
00553 Tensor<rank_, dim, Number>::unroll (Vector<Number2> &result) const
00554 {
00555 AssertDimension (result.size(),(Utilities::fixed_power<rank_, unsigned int>(dim)));
00556
00557 unsigned index = 0;
00558 unroll_recursion (result, index);
00559 }
00560
00561
00562
00563 template <int rank_, int dim, typename Number>
00564 template <typename Number2>
00565 inline
00566 void
00567 Tensor<rank_, dim, Number>::unroll_recursion (Vector<Number2> &result,
00568 unsigned int &index) const
00569 {
00570 for (unsigned i=0; i<dim; ++i)
00571 {
00572 operator[](i).unroll_recursion(result, index);
00573 }
00574 }
00575
00576
00577
00578 template <int rank_, int dim, typename Number>
00579 inline
00580 void Tensor<rank_,dim,Number>::clear ()
00581 {
00582 for (unsigned int i=0; i<dim; ++i)
00583 subtensor[i].clear();
00584 }
00585
00586
00587
00588 template <int rank_, int dim, typename Number>
00589 inline
00590 std::size_t
00591 Tensor<rank_,dim,Number>::memory_consumption ()
00592 {
00593 return sizeof(Tensor<rank_,dim,Number>);
00594 }
00595
00596
00597
00598 template <int rank_, int dim, typename Number>
00599 template <class Archive>
00600 inline
00601 void
00602 Tensor<rank_,dim,Number>::serialize(Archive & ar, const unsigned int)
00603 {
00604 ar & subtensor;
00605 }
00606
00607 #endif // DOXYGEN
00608
00609
00610
00611
00619 template <int rank_, int dim, typename Number>
00620 inline
00621 std::ostream & operator << (std::ostream &out, const Tensor<rank_,dim,Number> &p)
00622 {
00623 for (unsigned int i=0; i<dim-1; ++i)
00624 out << p[i] << ' ';
00625 out << p[dim-1];
00626
00627 return out;
00628 }
00629
00630 #ifndef DOXYGEN
00631
00635 template <int rank_>
00636 inline
00637 std::ostream & operator << (std::ostream &out, const Tensor<rank_,1> &p)
00638 {
00639 out << p[0];
00640
00641 return out;
00642 }
00643
00644 #endif // DOXYGEN
00645
00646
00654 template <int dim, typename Number>
00655 inline
00656 Number contract (const Tensor<1,dim,Number> &src1,
00657 const Tensor<1,dim,Number> &src2)
00658 {
00659 Number res = 0.;
00660 for (unsigned int i=0; i<dim; ++i)
00661 res += src1[i] * src2[i];
00662
00663 return res;
00664 }
00665
00666
00683 template <int dim, typename Number>
00684 inline
00685 Number
00686 operator * (const Tensor<1,dim,Number> &src1,
00687 const Tensor<1,dim,Number> &src2)
00688 {
00689 return contract(src1, src2);
00690 }
00691
00692
00700 template <int dim, typename Number>
00701 inline
00702 void contract (Tensor<1,dim,Number> &dest,
00703 const Tensor<2,dim,Number> &src1,
00704 const Tensor<1,dim,Number> &src2)
00705 {
00706 for (unsigned int i=0; i<dim; ++i)
00707 {
00708 dest[i] = src1[i][0] * src2[0];
00709 for (unsigned int j=1; j<dim; ++j)
00710 dest[i] += src1[i][j] * src2[j];
00711 }
00712 }
00713
00714
00731 template <int dim, typename Number>
00732 Tensor<1,dim,Number>
00733 operator * (const Tensor<2,dim,Number> &src1,
00734 const Tensor<1,dim,Number> &src2)
00735 {
00736 Tensor<1,dim,Number> dest (false);
00737 for (unsigned int i=0; i<dim; ++i)
00738 {
00739 dest[i] = src1[i][0] * src2[0];
00740 for (unsigned int j=1; j<dim; ++j)
00741 dest[i] += src1[i][j] * src2[j];
00742 }
00743 return dest;
00744 }
00745
00746
00754 template <int dim, typename Number>
00755 inline
00756 void contract (Tensor<1,dim,Number> &dest,
00757 const Tensor<1,dim,Number> &src1,
00758 const Tensor<2,dim,Number> &src2)
00759 {
00760 for (unsigned int i=0; i<dim; ++i)
00761 {
00762 dest[i] = src1[0] * src2[0][i];
00763 for (unsigned int j=1; j<dim; ++j)
00764 dest[i] += src1[j] * src2[j][i];
00765 }
00766 }
00767
00768
00785 template <int dim, typename Number>
00786 inline
00787 Tensor<1,dim,Number>
00788 operator * (const Tensor<1,dim,Number> &src1,
00789 const Tensor<2,dim,Number> &src2)
00790 {
00791 Tensor<1,dim,Number> dest (false);
00792 for (unsigned int i=0; i<dim; ++i)
00793 {
00794 dest[i] = src1[0] * src2[0][i];
00795 for (unsigned int j=1; j<dim; ++j)
00796 dest[i] += src1[j] * src2[j][i];
00797 }
00798 return dest;
00799 }
00800
00801
00809 template <int dim, typename Number>
00810 inline
00811 void contract (Tensor<2,dim,Number> &dest,
00812 const Tensor<2,dim,Number> &src1,
00813 const Tensor<2,dim,Number> &src2)
00814 {
00815 for (unsigned int i=0; i<dim; ++i)
00816 for (unsigned int j=0; j<dim; ++j)
00817 {
00818 dest[i][j] = src1[i][0] * src2[0][j];
00819 for (unsigned int k=1; k<dim; ++k)
00820 dest[i][j] += src1[i][k] * src2[k][j];
00821 }
00822 }
00823
00824
00825
00842 template <int dim, typename Number>
00843 inline
00844 Tensor<2,dim,Number>
00845 operator * (const Tensor<2,dim,Number> &src1,
00846 const Tensor<2,dim,Number> &src2)
00847 {
00848 Tensor<2,dim,Number> dest;
00849 for (unsigned int i=0; i<dim; ++i)
00850 for (unsigned int j=0; j<dim; ++j)
00851 for (unsigned int k=0; k<dim; ++k)
00852 dest[i][j] += src1[i][k] * src2[k][j];
00853 return dest;
00854 }
00855
00856
00871 template <int dim, typename Number>
00872 inline
00873 void contract (Tensor<2,dim,Number> &dest,
00874 const Tensor<2,dim,Number> &src1, const unsigned int index1,
00875 const Tensor<2,dim,Number> &src2, const unsigned int index2)
00876 {
00877 dest.clear ();
00878
00879 switch (index1)
00880 {
00881 case 1:
00882 switch (index2)
00883 {
00884 case 1:
00885 for (unsigned int i=0; i<dim; ++i)
00886 for (unsigned int j=0; j<dim; ++j)
00887 for (unsigned int k=0; k<dim; ++k)
00888 dest[i][j] += src1[k][i] * src2[k][j];
00889 break;
00890 case 2:
00891 for (unsigned int i=0; i<dim; ++i)
00892 for (unsigned int j=0; j<dim; ++j)
00893 for (unsigned int k=0; k<dim; ++k)
00894 dest[i][j] += src1[k][i] * src2[j][k];
00895 break;
00896
00897 default:
00898 Assert (false,
00899 (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
00900 };
00901 break;
00902 case 2:
00903 switch (index2)
00904 {
00905 case 1:
00906 for (unsigned int i=0; i<dim; ++i)
00907 for (unsigned int j=0; j<dim; ++j)
00908 for (unsigned int k=0; k<dim; ++k)
00909 dest[i][j] += src1[i][k] * src2[k][j];
00910 break;
00911 case 2:
00912 for (unsigned int i=0; i<dim; ++i)
00913 for (unsigned int j=0; j<dim; ++j)
00914 for (unsigned int k=0; k<dim; ++k)
00915 dest[i][j] += src1[i][k] * src2[j][k];
00916 break;
00917
00918 default:
00919 Assert (false,
00920 (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
00921 };
00922 break;
00923
00924 default:
00925 Assert (false, (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index1)));
00926 };
00927 }
00928
00929
00941 template <int dim, typename Number>
00942 inline
00943 void contract (Tensor<2,dim,Number> &dest,
00944 const Tensor<3,dim,Number> &src1, const unsigned int index1,
00945 const Tensor<1,dim,Number> &src2)
00946 {
00947 dest.clear ();
00948
00949 switch (index1)
00950 {
00951 case 1:
00952 for (unsigned int i=0; i<dim; ++i)
00953 for (unsigned int j=0; j<dim; ++j)
00954 for (unsigned int k=0; k<dim; ++k)
00955 dest[i][j] += src1[k][i][j] * src2[k];
00956 break;
00957
00958 case 2:
00959 for (unsigned int i=0; i<dim; ++i)
00960 for (unsigned int j=0; j<dim; ++j)
00961 for (unsigned int k=0; k<dim; ++k)
00962 dest[i][j] += src1[i][k][j] * src2[k];
00963 break;
00964
00965 case 3:
00966 for (unsigned int i=0; i<dim; ++i)
00967 for (unsigned int j=0; j<dim; ++j)
00968 for (unsigned int k=0; k<dim; ++k)
00969 dest[i][j] += src1[i][j][k] * src2[k];
00970 break;
00971
00972 default:
00973 Assert (false,
00974 (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index1)));
00975 };
00976 }
00977
00978
00986 template <int dim, typename Number>
00987 inline
00988 void contract (Tensor<3,dim,Number> &dest,
00989 const Tensor<3,dim,Number> &src1,
00990 const Tensor<2,dim,Number> &src2)
00991 {
00992 dest.clear ();
00993 for (unsigned int i=0; i<dim; ++i)
00994 for (unsigned int j=0; j<dim; ++j)
00995 for (unsigned int k=0; k<dim; ++k)
00996 for (unsigned int l=0; l<dim; ++l)
00997 dest[i][j][k] += src1[i][j][l] * src2[l][k];
00998 }
00999
01000
01014 template <int dim, typename Number>
01015 inline
01016 void contract (Tensor<3,dim,Number> &dest,
01017 const Tensor<3,dim,Number> &src1, const unsigned int index1,
01018 const Tensor<2,dim,Number> &src2, const unsigned int index2)
01019 {
01020 dest.clear ();
01021
01022 switch (index1)
01023 {
01024 case 1:
01025 switch (index2)
01026 {
01027 case 1:
01028 for (unsigned int i=0; i<dim; ++i)
01029 for (unsigned int j=0; j<dim; ++j)
01030 for (unsigned int k=0; k<dim; ++k)
01031 for (unsigned int l=0; l<dim; ++l)
01032 dest[i][j][k] += src1[l][i][j] * src2[l][k];
01033 break;
01034 case 2:
01035 for (unsigned int i=0; i<dim; ++i)
01036 for (unsigned int j=0; j<dim; ++j)
01037 for (unsigned int k=0; k<dim; ++k)
01038 for (unsigned int l=0; l<dim; ++l)
01039 dest[i][j][k] += src1[l][i][j] * src2[k][l];
01040 break;
01041 default:
01042 Assert (false,
01043 (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
01044 }
01045
01046 break;
01047 case 2:
01048 switch (index2)
01049 {
01050 case 1:
01051 for (unsigned int i=0; i<dim; ++i)
01052 for (unsigned int j=0; j<dim; ++j)
01053 for (unsigned int k=0; k<dim; ++k)
01054 for (unsigned int l=0; l<dim; ++l)
01055 dest[i][j][k] += src1[i][l][j] * src2[l][k];
01056 break;
01057 case 2:
01058 for (unsigned int i=0; i<dim; ++i)
01059 for (unsigned int j=0; j<dim; ++j)
01060 for (unsigned int k=0; k<dim; ++k)
01061 for (unsigned int l=0; l<dim; ++l)
01062 dest[i][j][k] += src1[i][l][j] * src2[k][l];
01063 break;
01064 default:
01065 Assert (false,
01066 (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
01067 }
01068
01069 break;
01070 case 3:
01071 switch (index2)
01072 {
01073 case 1:
01074 for (unsigned int i=0; i<dim; ++i)
01075 for (unsigned int j=0; j<dim; ++j)
01076 for (unsigned int k=0; k<dim; ++k)
01077 for (unsigned int l=0; l<dim; ++l)
01078 dest[i][j][k] += src1[i][j][l] * src2[l][k];
01079 break;
01080 case 2:
01081 for (unsigned int i=0; i<dim; ++i)
01082 for (unsigned int j=0; j<dim; ++j)
01083 for (unsigned int k=0; k<dim; ++k)
01084 for (unsigned int l=0; l<dim; ++l)
01085 dest[i][j][k] += src1[i][j][l] * src2[k][l];
01086 break;
01087 default:
01088 Assert (false,
01089 (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
01090 }
01091
01092 break;
01093 default:
01094 Assert (false,
01095 (typename Tensor<3,dim,Number>::ExcInvalidTensorIndex (index1)));
01096 }
01097 }
01098
01099
01116 template <int dim, typename Number>
01117 inline
01118 Tensor<3,dim,Number>
01119 operator * (const Tensor<3,dim,Number> &src1,
01120 const Tensor<2,dim,Number> &src2)
01121 {
01122 Tensor<3,dim,Number> dest;
01123 for (unsigned int i=0; i<dim; ++i)
01124 for (unsigned int j=0; j<dim; ++j)
01125 for (unsigned int k=0; k<dim; ++k)
01126 for (unsigned int l=0; l<dim; ++l)
01127 dest[i][j][k] += src1[i][j][l] * src2[l][k];
01128 return dest;
01129 }
01130
01131
01139 template <int dim, typename Number>
01140 inline
01141 void contract (Tensor<3,dim,Number> &dest,
01142 const Tensor<2,dim,Number> &src1,
01143 const Tensor<3,dim,Number> &src2)
01144 {
01145 dest.clear ();
01146 for (unsigned int i=0; i<dim; ++i)
01147 for (unsigned int j=0; j<dim; ++j)
01148 for (unsigned int k=0; k<dim; ++k)
01149 for (unsigned int l=0; l<dim; ++l)
01150 dest[i][j][k] += src1[i][l] * src2[l][j][k];
01151 }
01152
01153
01170 template <int dim, typename Number>
01171 inline
01172 Tensor<3,dim,Number>
01173 operator * (const Tensor<2,dim,Number> &src1,
01174 const Tensor<3,dim,Number> &src2)
01175 {
01176 Tensor<3,dim,Number> dest;
01177 for (unsigned int i=0; i<dim; ++i)
01178 for (unsigned int j=0; j<dim; ++j)
01179 for (unsigned int k=0; k<dim; ++k)
01180 for (unsigned int l=0; l<dim; ++l)
01181 dest[i][j][k] += src1[i][l] * src2[l][j][k];
01182 return dest;
01183 }
01184
01185
01193 template <int dim, typename Number>
01194 inline
01195 Tensor<4,dim,Number>
01196 operator * (const Tensor<3,dim,Number> &src1,
01197 const Tensor<3,dim,Number> &src2)
01198 {
01199 Tensor<4,dim,Number> dest;
01200 for (unsigned int i=0; i<dim; ++i)
01201 for (unsigned int j=0; j<dim; ++j)
01202 for (unsigned int k=0; k<dim; ++k)
01203 for (unsigned int l=0; l<dim; ++l)
01204 for (unsigned int m=0; m<dim; ++m)
01205 dest[i][j][k][l] += src1[i][j][m] * src2[m][k][l];
01206 return dest;
01207 }
01208
01209
01218 template <int dim, typename Number>
01219 inline
01220 void double_contract (Tensor<2,dim,Number> &dest,
01221 const Tensor<4,dim,Number> &src1,
01222 const Tensor<2,dim,Number> &src2)
01223 {
01224 dest.clear ();
01225 for (unsigned int i=0; i<dim; ++i)
01226 for (unsigned int j=0; j<dim; ++j)
01227 for (unsigned int k=0; k<dim; ++k)
01228 for (unsigned int l=0; l<dim; ++l)
01229 dest[i][j] += src1[i][j][k][l] * src2[k][l];
01230 }
01231
01232
01240 template <int dim, typename Number>
01241 inline
01242 Number contract3 (const Tensor<1,dim,Number>& u,
01243 const Tensor<2,dim,Number>& A,
01244 const Tensor<1,dim,Number>& v)
01245 {
01246 Number result = 0.;
01247
01248 for (unsigned int i=0; i<dim; ++i)
01249 for (unsigned int j=0; j<dim; ++j)
01250 result += u[i] * A[i][j] * v[j];
01251 return result;
01252 }
01253
01254
01262 template <int dim, typename Number>
01263 inline
01264 Number
01265 contract3 (const Tensor<1,dim,Number> &t1,
01266 const Tensor<3,dim,Number> &t2,
01267 const Tensor<2,dim,Number> &t3)
01268 {
01269 Number s = 0;
01270 for (unsigned int i=0; i<dim; ++i)
01271 for (unsigned int j=0; j<dim; ++j)
01272 for (unsigned int k=0; k<dim; ++k)
01273 s += t1[i] * t2[i][j][k] * t3[j][k];
01274 return s;
01275 }
01276
01277
01285 template <int dim, typename Number>
01286 inline
01287 Number
01288 contract3 (const Tensor<2,dim,Number> &t1,
01289 const Tensor<3,dim,Number> &t2,
01290 const Tensor<1,dim,Number> &t3)
01291 {
01292 Number s = 0;
01293 for (unsigned int i=0; i<dim; ++i)
01294 for (unsigned int j=0; j<dim; ++j)
01295 for (unsigned int k=0; k<dim; ++k)
01296 s += t1[i][j] * t2[i][j][k] * t3[k];
01297 return s;
01298 }
01299
01300
01308 template <int dim, typename Number>
01309 inline
01310 Number
01311 contract3 (const Tensor<2,dim,Number> &t1,
01312 const Tensor<4,dim,Number> &t2,
01313 const Tensor<2,dim,Number> &t3)
01314 {
01315 Number s = 0;
01316 for (unsigned int i=0; i<dim; ++i)
01317 for (unsigned int j=0; j<dim; ++j)
01318 for (unsigned int k=0; k<dim; ++k)
01319 for (unsigned int l=0; l<dim; ++l)
01320 s += t1[i][j] * t2[i][j][k][l] * t3[k][l];
01321 return s;
01322 }
01323
01324
01332 template <int dim, typename Number>
01333 void outer_product (Tensor<2,dim,Number> &dst,
01334 const Tensor<1,dim,Number> &src1,
01335 const Tensor<1,dim,Number> &src2)
01336 {
01337 for (unsigned int i=0; i<dim; ++i)
01338 for (unsigned int j=0; j<dim; ++j)
01339 dst[i][j] = src1[i] * src2[j];
01340 }
01341
01342
01350 template <int dim, typename Number>
01351 void outer_product (Tensor<3,dim,Number> &dst,
01352 const Tensor<1,dim,Number> &src1,
01353 const Tensor<2,dim,Number> &src2)
01354 {
01355 for (unsigned int i=0; i<dim; ++i)
01356 for (unsigned int j=0; j<dim; ++j)
01357 for (unsigned int k=0; k<dim; ++k)
01358 dst[i][j][k] = src1[i] * src2[j][k];
01359 }
01360
01361
01369 template <int dim, typename Number>
01370 void outer_product (Tensor<3,dim,Number> &dst,
01371 const Tensor<2,dim,Number> &src1,
01372 const Tensor<1,dim,Number> &src2)
01373 {
01374 for (unsigned int i=0; i<dim; ++i)
01375 for (unsigned int j=0; j<dim; ++j)
01376 for (unsigned int k=0; k<dim; ++k)
01377 dst[i][j][k] = src1[i][j] * src2[k];
01378 }
01379
01380
01392 template <int dim, typename Number>
01393 void outer_product (Tensor<1,dim,Number> &dst,
01394 const Number src1,
01395 const Tensor<1,dim,Number> &src2)
01396 {
01397 for (unsigned int i=0; i<dim; ++i)
01398 dst[i] = src1 * src2[i];
01399 }
01400
01401
01402
01414 template <int dim, typename Number>
01415 void outer_product (Tensor<1,dim,Number> &dst,
01416 const Tensor<1,dim,Number> src1,
01417 const Number src2)
01418 {
01419 for (unsigned int i=0; i<dim; ++i)
01420 dst[i] = src1[i] * src2;
01421 }
01422
01423
01436 template <int dim, typename Number>
01437 inline
01438 void
01439 cross_product (Tensor<1,dim,Number> &dst,
01440 const Tensor<1,dim,Number> &src)
01441 {
01442 Assert (dim==2, ExcInternalError());
01443
01444 dst[0] = src[1];
01445 dst[1] = -src[0];
01446 }
01447
01448
01459 template <int dim, typename Number>
01460 inline
01461 void
01462 cross_product (Tensor<1,dim,Number> &dst,
01463 const Tensor<1,dim,Number> &src1,
01464 const Tensor<1,dim,Number> &src2)
01465 {
01466 Assert (dim==3, ExcInternalError());
01467
01468 dst[0] = src1[1]*src2[2] - src1[2]*src2[1];
01469 dst[1] = src1[2]*src2[0] - src1[0]*src2[2];
01470 dst[2] = src1[0]*src2[1] - src1[1]*src2[0];
01471 }
01472
01473
01484 template <int dim, typename Number>
01485 inline
01486 Number
01487 scalar_product (const Tensor<2,dim,Number> &t1,
01488 const Tensor<2,dim,Number> &t2)
01489 {
01490 Number s = 0;
01491 for (unsigned int i=0; i<dim; ++i)
01492 for (unsigned int j=0; j<dim; ++j)
01493 s += t1[i][j] * t2[i][j];
01494 return s;
01495 }
01496
01497
01506 template <int rank, typename Number>
01507 inline
01508 Number determinant (const Tensor<rank,1,Number> &t)
01509 {
01510 return determinant(t[0]);
01511 }
01512
01513
01514
01523 template <typename Number>
01524 inline
01525 Number determinant (const Tensor<1,1,Number> &t)
01526 {
01527 return t[0];
01528 }
01529
01530
01539 template <typename Number>
01540 inline
01541 Number determinant (const Tensor<2,1,Number> &t)
01542 {
01543 return t[0][0];
01544 }
01545
01546
01547
01554 template <typename Number>
01555 inline
01556 Number determinant (const Tensor<2,2,Number> &t)
01557 {
01558 return ((t[0][0] * t[1][1]) - (t[1][0] * t[0][1]));
01559 }
01560
01561
01568 template <typename Number>
01569 inline
01570 Number determinant (const Tensor<2,3,Number> &t)
01571 {
01572
01573
01574
01575
01576
01577 const Number t4 = t[0][0]*t[1][1],
01578 t6 = t[0][0]*t[1][2],
01579 t8 = t[0][1]*t[1][0],
01580 t00 = t[0][2]*t[1][0],
01581 t01 = t[0][1]*t[2][0],
01582 t04 = t[0][2]*t[2][0],
01583 det = (t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
01584 t00*t[2][1]+t01*t[1][2]-t04*t[1][1]);
01585 return det;
01586 }
01587
01588
01596 template <int dim, typename Number>
01597 inline
01598 Number determinant (const Tensor<2,dim,Number> &t)
01599 {
01600
01601
01602
01603
01604
01605
01606
01607
01608 Number det = 0;
01609
01610 for (unsigned int k=0; k<dim; ++k)
01611 {
01612 Tensor<2,dim-1,Number> minor;
01613 for (unsigned int i=0; i<dim-1; ++i)
01614 for (unsigned int j=0; j<dim-1; ++j)
01615 minor[i][j] = t[i][j<k ? j : j+1];
01616
01617 const Number cofactor = std::pow (-1., static_cast<Number>(k+1)) *
01618 determinant (minor);
01619
01620 det += t[dim-1][k] * cofactor;
01621 }
01622
01623 return std::pow (-1., static_cast<Number>(dim)) * det;
01624 }
01625
01626
01627
01635 template <int dim, typename Number>
01636 Number trace (const Tensor<2,dim,Number> &d)
01637 {
01638 Number t=0;
01639 for (unsigned int i=0; i<dim; ++i)
01640 t += d[i][i];
01641 return t;
01642 }
01643
01644
01645
01655 template <int dim, typename Number>
01656 inline
01657 Tensor<2,dim,Number>
01658 invert (const Tensor<2,dim,Number> &t)
01659 {
01660 Number return_tensor [dim][dim];
01661 switch (dim)
01662 {
01663 case 1:
01664 return_tensor[0][0] = 1.0/t[0][0];
01665 break;
01666
01667 case 2:
01668
01669
01670 {
01671 const Number det = t[0][0]*t[1][1]-t[1][0]*t[0][1];
01672 const Number t4 = 1.0/det;
01673 return_tensor[0][0] = t[1][1]*t4;
01674 return_tensor[0][1] = -t[0][1]*t4;
01675 return_tensor[1][0] = -t[1][0]*t4;
01676 return_tensor[1][1] = t[0][0]*t4;
01677 break;
01678 }
01679
01680 case 3:
01681 {
01682 const Number t4 = t[0][0]*t[1][1],
01683 t6 = t[0][0]*t[1][2],
01684 t8 = t[0][1]*t[1][0],
01685 t00 = t[0][2]*t[1][0],
01686 t01 = t[0][1]*t[2][0],
01687 t04 = t[0][2]*t[2][0],
01688 det = (t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
01689 t00*t[2][1]+t01*t[1][2]-t04*t[1][1]),
01690 t07 = 1.0/det;
01691 return_tensor[0][0] = (t[1][1]*t[2][2]-t[1][2]*t[2][1])*t07;
01692 return_tensor[0][1] = (t[0][2]*t[2][1]-t[0][1]*t[2][2])*t07;
01693 return_tensor[0][2] = (t[0][1]*t[1][2]-t[0][2]*t[1][1])*t07;
01694 return_tensor[1][0] = (t[1][2]*t[2][0]-t[1][0]*t[2][2])*t07;
01695 return_tensor[1][1] = (t[0][0]*t[2][2]-t04)*t07;
01696 return_tensor[1][2] = (t00-t6)*t07;
01697 return_tensor[2][0] = (t[1][0]*t[2][1]-t[1][1]*t[2][0])*t07;
01698 return_tensor[2][1] = (t01-t[0][0]*t[2][1])*t07;
01699 return_tensor[2][2] = (t4-t8)*t07;
01700
01701 break;
01702 }
01703
01704
01705
01706
01707 default:
01708 AssertThrow (false, ExcNotImplemented());
01709 }
01710 return Tensor<2,dim,Number>(return_tensor);
01711 }
01712
01713
01714
01725 template <int dim, typename Number>
01726 inline
01727 Tensor<2,dim,Number>
01728 transpose (const Tensor<2,dim,Number> &t)
01729 {
01730 Number tt[dim][dim];
01731 for (unsigned int i=0; i<dim; ++i)
01732 {
01733 tt[i][i] = t[i][i];
01734 for (unsigned int j=i+1; j<dim; ++j)
01735 {
01736 tt[i][j] = t[j][i];
01737 tt[j][i] = t[i][j];
01738 };
01739 }
01740 return Tensor<2,dim,Number>(tt);
01741 }
01742
01743 #ifndef DOXYGEN
01744
01752 template <typename Number>
01753 inline
01754 Tensor<2,1,Number>
01755 transpose (const Tensor<2,1,Number> &t)
01756 {
01757 return t;
01758 }
01759
01760
01761
01762
01770 template <typename Number>
01771 inline
01772 Tensor<2,2,Number>
01773 transpose (const Tensor<2,2,Number> &t)
01774 {
01775 const Number x[2][2] = {{t[0][0], t[1][0]}, {t[0][1], t[1][1]}};
01776 return Tensor<2,2,Number>(x);
01777 }
01778
01779
01780
01781
01789 template <typename Number>
01790 inline
01791 Tensor<2,3,Number>
01792 transpose (const Tensor<2,3,Number> &t)
01793 {
01794 const Number x[3][3] = {{t[0][0], t[1][0], t[2][0]},
01795 {t[0][1], t[1][1], t[2][1]},
01796 {t[0][2], t[1][2], t[2][2]}};
01797 return Tensor<2,3,Number>(x);
01798 }
01799
01800 #endif // DOXYGEN
01801
01802
01811 template <int dim, typename Number>
01812 inline
01813 double
01814 l1_norm (const Tensor<2,dim,Number> &t)
01815 {
01816 double max = 0;
01817 for (unsigned int j=0; j<dim; ++j)
01818 {
01819 double sum = 0;
01820 for (unsigned int i=0; i<dim; ++i)
01821 sum += std::fabs(t[i][j]);
01822
01823 if (sum > max)
01824 max = sum;
01825 }
01826
01827 return max;
01828 }
01829
01830
01831
01840 template <int dim, typename Number>
01841 inline
01842 double
01843 linfty_norm (const Tensor<2,dim,Number> &t)
01844 {
01845 double max = 0;
01846 for (unsigned int i=0; i<dim; ++i)
01847 {
01848 double sum = 0;
01849 for (unsigned int j=0; j<dim; ++j)
01850 sum += std::fabs(t[i][j]);
01851
01852 if (sum > max)
01853 max = sum;
01854 }
01855
01856 return max;
01857 }
01858
01859
01860
01867 template <int rank, int dim, typename Number>
01868 inline
01869 Tensor<rank,dim,Number>
01870 operator * (const Tensor<rank,dim,Number> &t,
01871 const Number factor)
01872 {
01873 Tensor<rank,dim,Number> tt = t;
01874 tt *= factor;
01875 return tt;
01876 }
01877
01878
01879
01886 template <int rank, int dim, typename Number>
01887 inline
01888 Tensor<rank,dim,Number>
01889 operator * (const Number factor,
01890 const Tensor<rank,dim,Number> &t)
01891 {
01892 Tensor<rank,dim,Number> tt = t;
01893 tt *= factor;
01894 return tt;
01895 }
01896
01897
01898
01904 template <int rank, int dim, typename Number>
01905 inline
01906 Tensor<rank,dim,Number>
01907 operator / (const Tensor<rank,dim,Number> &t,
01908 const Number factor)
01909 {
01910 Tensor<rank,dim,Number> tt = t;
01911 tt /= factor;
01912 return tt;
01913 }
01914
01915
01916
01917
01924 template <int rank, int dim>
01925 inline
01926 Tensor<rank,dim>
01927 operator * (const Tensor<rank,dim> &t,
01928 const double factor)
01929 {
01930 Tensor<rank,dim> tt = t;
01931 tt *= factor;
01932 return tt;
01933 }
01934
01935
01936
01943 template <int rank, int dim>
01944 inline
01945 Tensor<rank,dim>
01946 operator * (const double factor,
01947 const Tensor<rank,dim> &t)
01948 {
01949 Tensor<rank,dim> tt = t;
01950 tt *= factor;
01951 return tt;
01952 }
01953
01954
01955
01961 template <int rank, int dim>
01962 inline
01963 Tensor<rank,dim>
01964 operator / (const Tensor<rank,dim> &t,
01965 const double factor)
01966 {
01967 Tensor<rank,dim> tt = t;
01968 tt /= factor;
01969 return tt;
01970 }
01971
01972 DEAL_II_NAMESPACE_CLOSE
01973
01974 #endif
01975