include/deal.II/base/tensor.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: tensor.h 25492 2012-05-05 15:11:20Z kormann @f$
00003 //
00004 //    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
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                                      // make the following class a
00302                                      // friend to this class. in principle,
00303                                      // it would suffice if otherrank==rank+1,
00304                                      // but then the compiler complains
00305                                      // that this be an explicit specialization
00306                                      // which is not what we want
00307                                      //
00308                                      // also, it would be sufficient to make
00309                                      // the function unroll_loops a friend,
00310                                      // but that seems to be impossible as well.
00311     template <int, int, typename> friend class Tensor;
00312 };
00313 
00314 
00315 /*--------------------------- Inline functions -----------------------------*/
00316 
00317 #ifndef DOXYGEN
00318 
00319 template <int rank_, int dim, typename Number>
00320 inline
00321 Tensor<rank_,dim,Number>::Tensor ()
00322 {
00323 // default constructor. not specifying an initializer list calls
00324 // the default constructor of the subobjects, which initialize them
00325 // selves. therefore, the tensor is set to zero this way
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 /* ----------------- Non-member functions operating on tensors. ------------ */
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                                 // use exactly the same expression with the
01573                                 // same order of operations as for the inverse
01574                                 // to let the compiler use common
01575                                 // subexpression elimination when using
01576                                 // determinant and inverse in nearby code
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                                    // compute the determinant using the
01601                                    // Laplace expansion of the
01602                                    // determinant. this may not be the most
01603                                    // efficient algorithm, but it does for
01604                                    // small n.
01605                                    //
01606                                    // for some algorithmic simplicity, we use
01607                                    // the expansion along the last row
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                                              // this is Maple output,
01669                                              // thus a bit unstructured
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                                         // if desired, take over the
01705                                         // inversion of a 4x4 tensor
01706                                         // from the FullMatrix
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 12:03:21 by doxygen 1.7.3