include/deal.II/base/symmetric_tensor.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: symmetric_tensor.h 25493 2012-05-05 18:36:25Z kormann @f$
00003 //
00004 //    Copyright (C) 2005, 2006, 2008, 2009, 2010, 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__symmetric_tensor_h
00013 #define __deal2__symmetric_tensor_h
00014 
00015 
00016 #include <deal.II/base/tensor.h>
00017 #include <deal.II/base/table_indices.h>
00018 
00019 DEAL_II_NAMESPACE_OPEN
00020 
00021 template <int rank, int dim, typename Number=double> class SymmetricTensor;
00022 
00023 template <int dim, typename Number> SymmetricTensor<2,dim,Number>
00024 unit_symmetric_tensor ();
00025 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
00026 deviator_tensor ();
00027 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
00028 identity_tensor ();
00029 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
00030 invert (const SymmetricTensor<4,dim,Number> &);
00031 template <int dim2, typename Number> Number
00032 trace (const SymmetricTensor<2,dim2,Number> &);
00033 
00034 template <int dim, typename Number> SymmetricTensor<2,dim,Number>
00035 deviator (const SymmetricTensor<2,dim,Number> &);
00036 template <int dim, typename Number> Number
00037 determinant (const SymmetricTensor<2,dim,Number> &);
00038 
00039 
00040 
00041 namespace internal
00042 {
00048   namespace SymmetricTensorAccessors
00049   {
00061     inline
00062     TableIndices<2> merge (const TableIndices<2> &previous_indices,
00063                            const unsigned int     new_index,
00064                            const unsigned int     position)
00065     {
00066       Assert (position < 2, ExcIndexRange (position, 0, 2));
00067 
00068       if (position == 0)
00069         return TableIndices<2>(new_index);
00070       else
00071         return TableIndices<2>(previous_indices[0], new_index);
00072     }
00073 
00074 
00075 
00087     inline
00088     TableIndices<4> merge (const TableIndices<4> &previous_indices,
00089                            const unsigned int     new_index,
00090                            const unsigned int     position)
00091     {
00092       Assert (position < 4, ExcIndexRange (position, 0, 4));
00093 
00094       switch (position)
00095         {
00096           case 0:
00097                 return TableIndices<4>(new_index);
00098           case 1:
00099                 return TableIndices<4>(previous_indices[0],
00100                                        new_index);
00101           case 2:
00102                 return TableIndices<4>(previous_indices[0],
00103                                        previous_indices[1],
00104                                        new_index);
00105           case 3:
00106                 return TableIndices<4>(previous_indices[0],
00107                                        previous_indices[1],
00108                                        previous_indices[2],
00109                                        new_index);
00110         }
00111       Assert (false, ExcInternalError());
00112       return TableIndices<4>();
00113     }
00114 
00115 
00130     template <int rank1, int rank2, int dim, typename Number>
00131     struct double_contraction_result
00132     {
00133       typedef ::SymmetricTensor<rank1+rank2-4,dim,Number> type;
00134     };
00135 
00136 
00151     template <int dim, typename Number>
00152     struct double_contraction_result<2,2,dim,Number>
00153     {
00154         typedef Number type;
00155     };
00156 
00157 
00158 
00178     template <int rank, int dim, typename Number>
00179     struct StorageType;
00180 
00185     template <int dim, typename Number>
00186     struct StorageType<2,dim,Number>
00187     {
00193         static const unsigned int
00194         n_independent_components = (dim*dim + dim)/2;
00195 
00200         typedef Tensor<1,n_independent_components,Number> base_tensor_type;
00201     };
00202 
00203 
00204 
00209     template <int dim, typename Number>
00210     struct StorageType<4,dim,Number>
00211     {
00219         static const unsigned int
00220         n_rank2_components = (dim*dim + dim)/2;
00221 
00226         static const unsigned int
00227         n_independent_components = (n_rank2_components *
00228                                     StorageType<2,dim,Number>::n_independent_components);
00229 
00239         typedef Tensor<2,n_rank2_components,Number> base_tensor_type;
00240     };
00241 
00242 
00243 
00250     template <int rank, int dim, bool constness, typename Number>
00251     struct AccessorTypes;
00252 
00261     template <int rank, int dim, typename Number>
00262     struct AccessorTypes<rank,dim,true,Number>
00263     {
00264         typedef const ::SymmetricTensor<rank,dim,Number> tensor_type;
00265 
00266         typedef Number reference;
00267     };
00268 
00278     template <int rank, int dim, typename Number>
00279     struct AccessorTypes<rank,dim,false,Number>
00280     {
00281         typedef ::SymmetricTensor<rank,dim,Number> tensor_type;
00282 
00283         typedef Number &reference;
00284     };
00285 
00286 
00322     template <int rank, int dim, bool constness, int P, typename Number>
00323     class Accessor
00324     {
00325       public:
00330         typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
00331         typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
00332 
00333       private:
00374         Accessor (tensor_type              &tensor,
00375                   const TableIndices<rank> &previous_indices);
00376 
00382         Accessor ();
00383 
00389         Accessor (const Accessor &a);
00390 
00391       public:
00392 
00396         Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i);
00397 
00398       private:
00403         tensor_type             &tensor;
00404         const TableIndices<rank> previous_indices;
00405 
00406                                          // declare some other classes
00407                                          // as friends. make sure to
00408                                          // work around bugs in some
00409                                          // compilers
00410 #ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG
00411         template <int,int,typename> friend class SymmetricTensor;
00412         template <int,int,bool,int,typename>
00413         friend class Accessor;
00414 #  ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
00415         friend class ::SymmetricTensor<rank,dim,Number>;
00416         friend class Accessor<rank,dim,constness,P+1,Number>;
00417 #  endif
00418 #else
00419         friend class SymmetricTensor<rank,dim,Number>;
00420         friend class Accessor<rank,dim,constness,P+1,Number>;
00421 #endif
00422     };
00423 
00424 
00425 
00436     template <int rank, int dim, bool constness, typename Number>
00437     class Accessor<rank,dim,constness,1,Number>
00438     {
00439       public:
00444         typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
00445         typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
00446 
00447       private:
00493         Accessor (tensor_type              &tensor,
00494                   const TableIndices<rank> &previous_indices);
00495 
00501         Accessor ();
00502 
00508         Accessor (const Accessor &a);
00509 
00510       public:
00511 
00515         reference operator [] (const unsigned int);
00516 
00517       private:
00522         tensor_type             &tensor;
00523         const TableIndices<rank> previous_indices;
00524 
00525                                          // declare some other classes
00526                                          // as friends. make sure to
00527                                          // work around bugs in some
00528                                          // compilers
00529 #ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG
00530         template <int,int,typename> friend class SymmetricTensor;
00531         template <int,int,bool,int,typename>
00532         friend class SymmetricTensorAccessors::Accessor;
00533 #  ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
00534         friend class ::SymmetricTensor<rank,dim,Number>;
00535         friend class SymmetricTensorAccessors::Accessor<rank,dim,constness,2,Number>;
00536 #  endif
00537 #else
00538         friend class SymmetricTensor<rank,dim,Number>;
00539         friend class Accessor<rank,dim,constness,2,Number>;
00540 #endif
00541     };
00542   }
00543 }
00544 
00545 
00546 
00613 template <int rank, int dim, typename Number>
00614 class SymmetricTensor
00615 {
00616   public:
00634     static const unsigned int dimension = dim;
00635 
00644     static const unsigned int n_independent_components
00645     = internal::SymmetricTensorAccessors::StorageType<rank,dim,Number>::
00646       n_independent_components;
00647 
00652     SymmetricTensor ();
00653 
00672     SymmetricTensor (const Tensor<2,dim,Number> &t);
00673 
00701     SymmetricTensor (const Number (&array) [internal::SymmetricTensorAccessors::StorageType<rank,dim,Number>::n_independent_components]);
00702 
00706     SymmetricTensor & operator = (const SymmetricTensor &);
00707 
00720     SymmetricTensor & operator = (const Number d);
00721 
00728     operator Tensor<rank,dim,Number> () const;
00729 
00733     bool operator == (const SymmetricTensor &) const;
00734 
00738     bool operator != (const SymmetricTensor &) const;
00739 
00743     SymmetricTensor & operator += (const SymmetricTensor &);
00744 
00748     SymmetricTensor & operator -= (const SymmetricTensor &);
00749 
00755     SymmetricTensor & operator *= (const Number factor);
00756 
00761     SymmetricTensor & operator /= (const Number factor);
00762 
00769     SymmetricTensor   operator + (const SymmetricTensor &s) const;
00770 
00777     SymmetricTensor   operator - (const SymmetricTensor &s) const;
00778 
00783     SymmetricTensor   operator - () const;
00784 
00836   typename internal::SymmetricTensorAccessors::double_contraction_result<rank,2,dim,Number>::type
00837   operator * (const SymmetricTensor<2,dim,Number> &s) const;
00838 
00845     typename internal::SymmetricTensorAccessors::double_contraction_result<rank,4,dim,Number>::type
00846     operator * (const SymmetricTensor<4,dim,Number> &s) const;
00847 
00852     Number & operator() (const TableIndices<rank> &indices);
00853 
00867     Number operator() (const TableIndices<rank> &indices) const;
00868 
00874     internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
00875     operator [] (const unsigned int row) const;
00876 
00882     internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
00883     operator [] (const unsigned int row);
00884 
00890     Number
00891     operator [] (const TableIndices<rank> &indices) const;
00892 
00898     Number &
00899     operator [] (const TableIndices<rank> &indices);
00900 
00910     Number
00911     access_raw_entry (const unsigned int unrolled_index) const;
00912 
00922     Number &
00923     access_raw_entry (const unsigned int unrolled_index);
00924 
00940     Number norm () const;
00941 
00954     static
00955     unsigned int
00956     component_to_unrolled_index (const TableIndices<rank> &indices);
00957 
00968     static
00969     TableIndices<rank>
00970     unrolled_to_component_indices (const unsigned int i);
00971 
00990     void clear ();
00991 
00998     static std::size_t memory_consumption ();
00999 
01004     template <class Archive>
01005     void serialize(Archive & ar, const unsigned int version);
01006 
01007   private:
01012     typedef
01013     internal::SymmetricTensorAccessors::StorageType<rank,dim,Number>
01014     base_tensor_descriptor;
01015 
01020     typedef typename base_tensor_descriptor::base_tensor_type base_tensor_type;
01021 
01026     base_tensor_type data;
01027 
01031     template <int, int, typename> friend class SymmetricTensor;
01032 
01036     template <int dim2, typename Number2>
01037     friend Number2 trace (const SymmetricTensor<2,dim2,Number2> &d);
01038 
01039     template <int dim2, typename Number2>
01040     friend Number2 determinant (const SymmetricTensor<2,dim2,Number2> &t);
01041 
01042     template <int dim2, typename Number2>
01043     friend SymmetricTensor<2,dim2,Number2>
01044     deviator (const SymmetricTensor<2,dim2,Number2> &t);
01045 
01046     template <int dim2, typename Number2>
01047     friend SymmetricTensor<2,dim2,Number2> unit_symmetric_tensor ();
01048 
01049     template <int dim2, typename Number2>
01050     friend SymmetricTensor<4,dim2,Number2> deviator_tensor ();
01051 
01052     template <int dim2, typename Number2>
01053     friend SymmetricTensor<4,dim2,Number2> identity_tensor ();
01054 
01055     template <int dim2, typename Number2>
01056     friend SymmetricTensor<4,dim2,Number2> invert (const SymmetricTensor<4,dim2,Number2> &);
01057 };
01058 
01059 
01060 
01061 // ------------------------- inline functions ------------------------
01062 
01063 #ifndef DOXYGEN
01064 
01065 namespace internal
01066 {
01067   namespace SymmetricTensorAccessors
01068   {
01069     template <int rank, int dim, bool constness, int P, typename Number>
01070     Accessor<rank,dim,constness,P,Number>::
01071     Accessor (tensor_type              &tensor,
01072               const TableIndices<rank> &previous_indices)
01073                     :
01074                     tensor (tensor),
01075                     previous_indices (previous_indices)
01076     {}
01077 
01078 
01079 
01080     template <int rank, int dim, bool constness, int P, typename Number>
01081     Accessor<rank,dim,constness,P-1,Number>
01082     Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i)
01083     {
01084       return Accessor<rank,dim,constness,P-1,Number> (tensor,
01085                                                       merge (previous_indices, i, rank-P));
01086     }
01087 
01088 
01089 
01090     template <int rank, int dim, bool constness, typename Number>
01091     Accessor<rank,dim,constness,1,Number>::
01092     Accessor (tensor_type              &tensor,
01093               const TableIndices<rank> &previous_indices)
01094                     :
01095                     tensor (tensor),
01096                     previous_indices (previous_indices)
01097     {}
01098 
01099 
01100 
01101     template <int rank, int dim, bool constness, typename Number>
01102     typename Accessor<rank,dim,constness,1,Number>::reference
01103     Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i)
01104     {
01105       return tensor(merge (previous_indices, i, rank-1));
01106     }
01107 
01108 
01109   }
01110 }
01111 
01112 
01113 
01114 template <int rank, int dim, typename Number>
01115 inline
01116 SymmetricTensor<rank,dim,Number>::SymmetricTensor ()
01117 {}
01118 
01119 
01120 
01121 template <int rank, int dim, typename Number>
01122 inline
01123 SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Tensor<2,dim,Number> &t)
01124 {
01125   Assert (rank == 2, ExcNotImplemented());
01126   switch (dim)
01127     {
01128     case 2:
01129       Assert (t[0][1] == t[1][0], ExcInternalError());
01130 
01131       data[0] = t[0][0];
01132       data[1] = t[1][1];
01133       data[2] = t[0][1];
01134 
01135       break;
01136     case 3:
01137       Assert (t[0][1] == t[1][0], ExcInternalError());
01138       Assert (t[0][2] == t[2][0], ExcInternalError());
01139       Assert (t[1][2] == t[2][1], ExcInternalError());
01140 
01141       data[0] = t[0][0];
01142       data[1] = t[1][1];
01143       data[2] = t[2][2];
01144       data[3] = t[0][1];
01145       data[4] = t[0][2];
01146       data[5] = t[1][2];
01147 
01148       break;
01149     default:
01150       Assert (false, ExcNotImplemented());
01151     }
01152 }
01153 
01154 
01155 
01156 template <int rank, int dim, typename Number>
01157 inline
01158 SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Number (&array) [internal::SymmetricTensorAccessors::StorageType<rank,dim,Number>::n_independent_components])
01159                 :
01160                 data (array)
01161 {}
01162 
01163 
01164 
01165 template <int rank, int dim, typename Number>
01166 inline
01167 SymmetricTensor<rank,dim,Number> &
01168 SymmetricTensor<rank,dim,Number>::operator = (const SymmetricTensor<rank,dim,Number> &t)
01169 {
01170   data = t.data;
01171   return *this;
01172 }
01173 
01174 
01175 
01176 template <int rank, int dim, typename Number>
01177 inline
01178 SymmetricTensor<rank,dim,Number> &
01179 SymmetricTensor<rank,dim,Number>::operator = (const Number d)
01180 {
01181   Assert (d==0, ExcMessage ("Only assignment with zero is allowed"));
01182 
01183   data = 0;
01184 
01185   return *this;
01186 }
01187 
01188 
01189 
01190                                 // helper function to convert symmetric tensor
01191                                 // to generic tensor
01192 namespace internal
01193 {
01194   template <typename Number>
01195   inline
01196   Tensor<2,1,Number>
01197   conversion (const Tensor<1,1,Number> &data)
01198   {
01199     const Number t[1][1] = {{data[0]}};
01200     return Tensor<2,1,Number>(t);
01201   }
01202 
01203   template <typename Number>
01204   inline
01205   Tensor<2,2,Number>
01206   conversion (const Tensor<1,3,Number> &data)
01207   {
01208     const Number t[2][2] = {{data[0], data[2]},
01209                             {data[2], data[1]}};
01210     return Tensor<2,2,Number>(t);
01211   }
01212 
01213   template <typename Number>
01214   inline
01215   Tensor<2,3,Number>
01216   conversion (const Tensor<1,6,Number> &data)
01217   {
01218     const Number t[3][3] = {{data[0], data[3], data[4]},
01219                             {data[3], data[1], data[5]},
01220                             {data[4], data[5], data[2]}};
01221     return Tensor<2,3,Number>(t);
01222   }
01223 }
01224 
01225 
01226 
01227 template <int rank, int dim, typename Number>
01228 inline
01229 SymmetricTensor<rank,dim,Number>::
01230 operator Tensor<rank,dim,Number> () const
01231 {
01232   Assert (rank == 2, ExcNotImplemented());
01233   return internal::conversion(data);
01234 }
01235 
01236 
01237 
01238 template <int rank, int dim, typename Number>
01239 inline
01240 bool
01241 SymmetricTensor<rank,dim,Number>::operator ==
01242 (const SymmetricTensor<rank,dim,Number> &t) const
01243 {
01244   return data == t.data;
01245 }
01246 
01247 
01248 
01249 template <int rank, int dim, typename Number>
01250 inline
01251 bool
01252 SymmetricTensor<rank,dim,Number>::operator !=
01253 (const SymmetricTensor<rank,dim,Number> &t) const
01254 {
01255   return data != t.data;
01256 }
01257 
01258 
01259 
01260 template <int rank, int dim, typename Number>
01261 inline
01262 SymmetricTensor<rank,dim,Number> &
01263 SymmetricTensor<rank,dim,Number>::operator +=
01264 (const SymmetricTensor<rank,dim,Number> &t)
01265 {
01266   data += t.data;
01267   return *this;
01268 }
01269 
01270 
01271 
01272 template <int rank, int dim, typename Number>
01273 inline
01274 SymmetricTensor<rank,dim,Number> &
01275 SymmetricTensor<rank,dim,Number>::operator -=
01276 (const SymmetricTensor<rank,dim,Number> &t)
01277 {
01278   data -= t.data;
01279   return *this;
01280 }
01281 
01282 
01283 
01284 template <int rank, int dim, typename Number>
01285 inline
01286 SymmetricTensor<rank,dim,Number> &
01287 SymmetricTensor<rank,dim,Number>::operator *= (const Number d)
01288 {
01289   data *= d;
01290   return *this;
01291 }
01292 
01293 
01294 
01295 template <int rank, int dim, typename Number>
01296 inline
01297 SymmetricTensor<rank,dim,Number> &
01298 SymmetricTensor<rank,dim,Number>::operator /= (const Number d)
01299 {
01300   data /= d;
01301   return *this;
01302 }
01303 
01304 
01305 
01306 template <int rank, int dim, typename Number>
01307 inline
01308 SymmetricTensor<rank,dim,Number>
01309 SymmetricTensor<rank,dim,Number>::operator + (const SymmetricTensor &t) const
01310 {
01311   SymmetricTensor tmp = *this;
01312   tmp.data += t.data;
01313   return tmp;
01314 }
01315 
01316 
01317 
01318 template <int rank, int dim, typename Number>
01319 inline
01320 SymmetricTensor<rank,dim,Number>
01321 SymmetricTensor<rank,dim,Number>::operator - (const SymmetricTensor &t) const
01322 {
01323   SymmetricTensor tmp = *this;
01324   tmp.data -= t.data;
01325   return tmp;
01326 }
01327 
01328 
01329 
01330 template <int rank, int dim, typename Number>
01331 inline
01332 SymmetricTensor<rank,dim,Number>
01333 SymmetricTensor<rank,dim,Number>::operator - () const
01334 {
01335   SymmetricTensor tmp = *this;
01336   tmp.data = -tmp.data;
01337   return tmp;
01338 }
01339 
01340 
01341 
01342 template <int rank, int dim, typename Number>
01343 inline
01344 void
01345 SymmetricTensor<rank,dim,Number>::clear ()
01346 {
01347   data.clear ();
01348 }
01349 
01350 
01351 
01352 template <int rank, int dim, typename Number>
01353 inline
01354 std::size_t
01355 SymmetricTensor<rank,dim,Number>::memory_consumption ()
01356 {
01357   return
01358     internal::SymmetricTensorAccessors::StorageType<rank,dim,Number>::memory_consumption ();
01359 }
01360 
01361 
01362 
01363 namespace internal
01364 {
01365 
01366   template <int dim, typename Number>
01367   inline
01368   typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number>::type
01369   perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
01370                               const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &sdata)
01371   {
01372     switch (dim)
01373       {
01374       case 1:
01375         return data[0] * sdata[0];
01376       case 2:
01377         return (data[0] * sdata[0] +
01378                 data[1] * sdata[1] +
01379                 2*data[2] * sdata[2]);
01380       case 3:
01381         return (data[0] * sdata[0] +
01382                 data[1] * sdata[1] +
01383                 data[2] * sdata[2] +
01384                 2*data[3] * sdata[3] +
01385                 2*data[4] * sdata[4] +
01386                 2*data[5] * sdata[5]);
01387       default:
01388         Number sum = 0;
01389         for (unsigned int d=0; d<dim; ++d)
01390           sum += data[d] * sdata[d];
01391         for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
01392           sum += Number(2.) * data[d] * sdata[d];
01393         return sum;
01394       }
01395   }
01396 
01397 
01398 
01399   template <int dim, typename Number>
01400   inline
01401   typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number>::type
01402   perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
01403                               const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &sdata)
01404   {
01405     Number tmp [SymmetricTensorAccessors::StorageType<2,dim,Number>::n_independent_components];
01406     switch (dim)
01407       {
01408       case 1:
01409         tmp[0] = data[0][0] * sdata[0];
01410         break;
01411       case 2:
01412         for (unsigned int i=0; i<3; ++i)
01413           tmp[i] = (data[i][0] * sdata[0] +
01414                     data[i][1] * sdata[1] +
01415                     2 * data[i][2] * sdata[2]);
01416         break;
01417       case 3:
01418         for (unsigned int i=0; i<6; ++i)
01419           tmp[i] = (data[i][0] * sdata[0] +
01420                     data[i][1] * sdata[1] +
01421                     data[i][2] * sdata[2] +
01422                     2 * data[i][3] * sdata[3] +
01423                     2 * data[i][4] * sdata[4] +
01424                     2 * data[i][5] * sdata[5]);
01425         break;
01426       default:
01427         Assert (false, ExcNotImplemented());
01428       }
01429     return SymmetricTensor<2,dim,Number>(tmp);
01430   }
01431 
01432 
01433 
01434   template <int dim, typename Number>
01435   inline
01436   typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type
01437   perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
01438                               const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &sdata)
01439   {
01440     typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type tmp;
01441     switch (dim)
01442       {
01443       case 1:
01444         tmp[0] = data[0] * sdata[0][0];
01445         break;
01446       case 2:
01447         for (unsigned int i=0; i<3; ++i)
01448           tmp[i] = (data[0] * sdata[0][i] +
01449                     data[1] * sdata[1][i] +
01450                     2 * data[2] * sdata[2][i]);
01451         break;
01452       case 3:
01453         for (unsigned int i=0; i<6; ++i)
01454           tmp[i] = (data[0] * sdata[0][i] +
01455                     data[1] * sdata[1][i] +
01456                     data[2] * sdata[2][i] +
01457                     2 * data[3] * sdata[3][i] +
01458                     2 * data[4] * sdata[4][i] +
01459                     2 * data[5] * sdata[5][i]);
01460         break;
01461       default:
01462         Assert (false, ExcNotImplemented());
01463       }
01464     return tmp;
01465   }
01466 
01467 
01468 
01469   template <int dim, typename Number>
01470   inline
01471   typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type
01472   perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
01473                               const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &sdata)
01474   {
01475     typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type tmp;
01476     switch (dim)
01477       {
01478       case 1:
01479         tmp[0][0] = data[0][0] * sdata[0][0];
01480         break;
01481       case 2:
01482         for (unsigned int i=0; i<3; ++i)
01483           for (unsigned int j=0; j<3; ++j)
01484             tmp[i][j] = (data[i][0] * sdata[0][j] +
01485                          data[i][1] * sdata[1][j] +
01486                          2*data[i][2] * sdata[2][j]);
01487         break;
01488       case 3:
01489         for (unsigned int i=0; i<6; ++i)
01490           for (unsigned int j=0; j<6; ++j)
01491             tmp[i][j] = (data[i][0] * sdata[0][j] +
01492                          data[i][1] * sdata[1][j] +
01493                          data[i][2] * sdata[2][j] +
01494                          2*data[i][3] * sdata[3][j] +
01495                          2*data[i][4] * sdata[4][j] +
01496                          2*data[i][5] * sdata[5][j]);
01497         break;
01498       default:
01499         Assert (false, ExcNotImplemented());
01500       }
01501     return tmp;
01502   }
01503 
01504 } // end of namespace internal
01505 
01506 
01507 
01508 template <int rank, int dim, typename Number>
01509 inline
01510 typename internal::SymmetricTensorAccessors::double_contraction_result<rank,2,dim,Number>::type
01511 SymmetricTensor<rank,dim,Number>::operator * (const SymmetricTensor<2,dim,Number> &s) const
01512 {
01513                                 // need to have two different function calls
01514                                 // because a scalar and rank-2 tensor are not
01515                                 // the same data type (see internal function
01516                                 // above)
01517   return internal::perform_double_contraction<dim,Number> (data, s.data);
01518 }
01519 
01520 
01521 
01522 template <int rank, int dim, typename Number>
01523 inline
01524 typename internal::SymmetricTensorAccessors::double_contraction_result<rank,4,dim,Number>::type
01525 SymmetricTensor<rank,dim,Number>::operator * (const SymmetricTensor<4,dim,Number> &s) const
01526 {
01527   typename internal::SymmetricTensorAccessors::
01528     double_contraction_result<rank,4,dim,Number>::type tmp;
01529   tmp.data = internal::perform_double_contraction<dim,Number> (data,s.data);
01530   return tmp;
01531 }
01532 
01533 
01534 
01535                                 // internal namespace to switch between the
01536                                 // access of different tensors. There used to
01537                                 // be explicit instantiations before for
01538                                 // different ranks and dimensions, but since
01539                                 // we now allow for templates on the data
01540                                 // type, and since we cannot partially
01541                                 // specialize the implementation, this got
01542                                 // into a separate namespace
01543 namespace internal
01544 {
01545   template <int dim, typename Number>
01546   inline
01547   Number &
01548   symmetric_tensor_access (const TableIndices<2> &indices,
01549                            typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data)
01550   {
01551     switch (dim)
01552       {
01553       case 1:
01554         return data[0];
01555 
01556       case 2:
01557                                    // first treat the main diagonal
01558                                    // elements, which are stored
01559                                    // consecutively at the beginning
01560         if (indices[0] == indices[1])
01561           return data[indices[0]];
01562 
01563                                    // the rest is messier and requires a few
01564                                    // switches. at least for the 2x2 case it
01565                                    // is reasonably simple
01566         Assert (((indices[0]==1) && (indices[1]==0)) ||
01567                 ((indices[0]==0) && (indices[1]==1)),
01568                 ExcInternalError());
01569         return data[2];
01570 
01571       case 3:
01572                                    // first treat the main diagonal
01573                                    // elements, which are stored
01574                                    // consecutively at the beginning
01575         if (indices[0] == indices[1])
01576           return data[indices[0]];
01577 
01578                                    // the rest is messier and requires a few
01579                                    // switches, but simpler if we just sort
01580                                    // our indices
01581         {
01582           TableIndices<2> sorted_indices (indices);
01583           sorted_indices.sort ();
01584 
01585           if ((sorted_indices[0]==0) && (sorted_indices[1]==1))
01586             return data[3];
01587           else if ((sorted_indices[0]==0) && (sorted_indices[1]==2))
01588             return data[4];
01589           else if ((sorted_indices[0]==1) && (sorted_indices[1]==2))
01590             return data[5];
01591           else
01592             Assert (false, ExcInternalError());
01593         }
01594       }
01595 
01596     static Number dummy_but_referenceable = Number();
01597     return dummy_but_referenceable;
01598   }
01599 
01600 
01601 
01602   template <int dim, typename Number>
01603   inline
01604   Number
01605   symmetric_tensor_access (const TableIndices<2> &indices,
01606                            const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data)
01607   {
01608     switch (dim)
01609       {
01610       case 1:
01611         return data[0];
01612 
01613       case 2:
01614                                    // first treat the main diagonal
01615                                    // elements, which are stored
01616                                    // consecutively at the beginning
01617         if (indices[0] == indices[1])
01618           return data[indices[0]];
01619 
01620                                    // the rest is messier and requires a few
01621                                    // switches. at least for the 2x2 case it
01622                                    // is reasonably simple
01623         Assert (((indices[0]==1) && (indices[1]==0)) ||
01624                 ((indices[0]==0) && (indices[1]==1)),
01625                 ExcInternalError());
01626         return data[2];
01627 
01628       case 3:
01629                                    // first treat the main diagonal
01630                                    // elements, which are stored
01631                                    // consecutively at the beginning
01632         if (indices[0] == indices[1])
01633           return data[indices[0]];
01634 
01635                                    // the rest is messier and requires a few
01636                                    // switches, but simpler if we just sort
01637                                    // our indices
01638         {
01639           TableIndices<2> sorted_indices (indices);
01640           sorted_indices.sort ();
01641 
01642           if ((sorted_indices[0]==0) && (sorted_indices[1]==1))
01643             return data[3];
01644           else if ((sorted_indices[0]==0) && (sorted_indices[1]==2))
01645             return data[4];
01646           else if ((sorted_indices[0]==1) && (sorted_indices[1]==2))
01647             return data[5];
01648           else
01649             Assert (false, ExcInternalError());
01650         }
01651       }
01652 
01653     static Number dummy_but_referenceable = 0;
01654     return dummy_but_referenceable;
01655   }
01656 
01657 
01658 
01659   template <int dim, typename Number>
01660   inline
01661   Number &
01662   symmetric_tensor_access (const TableIndices<4> &indices,
01663                            typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data)
01664   {
01665     switch (dim)
01666       {
01667         case 1:
01668               return data[0][0];
01669 
01670         case 2:
01671                                                // each entry of the tensor can be
01672                                                // thought of as an entry in a
01673                                                // matrix that maps the rolled-out
01674                                                // rank-2 tensors into rolled-out
01675                                                // rank-2 tensors. this is the
01676                                                // format in which we store rank-4
01677                                                // tensors. determine which
01678                                                // position the present entry is
01679                                                // stored in
01680         {
01681           unsigned int base_index[2] ;
01682           if ((indices[0] == 0) && (indices[1] == 0))
01683             base_index[0] = 0;
01684           else if ((indices[0] == 1) && (indices[1] == 1))
01685             base_index[0] = 1;
01686           else
01687             base_index[0] = 2;
01688 
01689           if ((indices[2] == 0) && (indices[3] == 0))
01690             base_index[1] = 0;
01691           else if ((indices[2] == 1) && (indices[3] == 1))
01692             base_index[1] = 1;
01693           else
01694             base_index[1] = 2;
01695 
01696           return data[base_index[0]][base_index[1]];
01697         }
01698 
01699         case 3:
01700                                                // each entry of the tensor can be
01701                                                // thought of as an entry in a
01702                                                // matrix that maps the rolled-out
01703                                                // rank-2 tensors into rolled-out
01704                                                // rank-2 tensors. this is the
01705                                                // format in which we store rank-4
01706                                                // tensors. determine which
01707                                                // position the present entry is
01708                                                // stored in
01709         {
01710           unsigned int base_index[2] ;
01711           if ((indices[0] == 0) && (indices[1] == 0))
01712             base_index[0] = 0;
01713           else if ((indices[0] == 1) && (indices[1] == 1))
01714             base_index[0] = 1;
01715           else if ((indices[0] == 2) && (indices[1] == 2))
01716             base_index[0] = 2;
01717           else if (((indices[0] == 0) && (indices[1] == 1)) ||
01718                    ((indices[0] == 1) && (indices[1] == 0)))
01719             base_index[0] = 3;
01720           else if (((indices[0] == 0) && (indices[1] == 2)) ||
01721                    ((indices[0] == 2) && (indices[1] == 0)))
01722             base_index[0] = 4;
01723           else
01724             {
01725               Assert (((indices[0] == 1) && (indices[1] == 2)) ||
01726                       ((indices[0] == 2) && (indices[1] == 1)),
01727                       ExcInternalError());
01728               base_index[0] = 5;
01729             }
01730 
01731           if ((indices[2] == 0) && (indices[3] == 0))
01732             base_index[1] = 0;
01733           else if ((indices[2] == 1) && (indices[3] == 1))
01734             base_index[1] = 1;
01735           else if ((indices[2] == 2) && (indices[3] == 2))
01736             base_index[1] = 2;
01737           else if (((indices[2] == 0) && (indices[3] == 1)) ||
01738                    ((indices[2] == 1) && (indices[3] == 0)))
01739             base_index[1] = 3;
01740           else if (((indices[2] == 0) && (indices[3] == 2)) ||
01741                    ((indices[2] == 2) && (indices[3] == 0)))
01742             base_index[1] = 4;
01743           else
01744             {
01745               Assert (((indices[2] == 1) && (indices[3] == 2)) ||
01746                       ((indices[2] == 2) && (indices[3] == 1)),
01747                       ExcInternalError());
01748               base_index[1] = 5;
01749             }
01750 
01751           return data[base_index[0]][base_index[1]];
01752         }
01753 
01754         default:
01755               Assert (false, ExcNotImplemented());
01756       }
01757 
01758     static Number dummy;
01759     return dummy;
01760   }
01761 
01762 
01763   template <int dim, typename Number>
01764   inline
01765   Number
01766   symmetric_tensor_access (const TableIndices<4> &indices,
01767                            const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data)
01768   {
01769     switch (dim)
01770       {
01771       case 1:
01772         return data[0][0];
01773 
01774       case 2:
01775                                    // each entry of the tensor can be
01776                                    // thought of as an entry in a
01777                                    // matrix that maps the rolled-out
01778                                    // rank-2 tensors into rolled-out
01779                                    // rank-2 tensors. this is the
01780                                    // format in which we store rank-4
01781                                    // tensors. determine which
01782                                    // position the present entry is
01783                                    // stored in
01784         {
01785           unsigned int base_index[2] ;
01786           if ((indices[0] == 0) && (indices[1] == 0))
01787             base_index[0] = 0;
01788           else if ((indices[0] == 1) && (indices[1] == 1))
01789             base_index[0] = 1;
01790           else
01791             base_index[0] = 2;
01792 
01793           if ((indices[2] == 0) && (indices[3] == 0))
01794             base_index[1] = 0;
01795           else if ((indices[2] == 1) && (indices[3] == 1))
01796             base_index[1] = 1;
01797           else
01798             base_index[1] = 2;
01799 
01800           return data[base_index[0]][base_index[1]];
01801         }
01802 
01803       case 3:
01804                                    // each entry of the tensor can be
01805                                    // thought of as an entry in a
01806                                    // matrix that maps the rolled-out
01807                                    // rank-2 tensors into rolled-out
01808                                    // rank-2 tensors. this is the
01809                                    // format in which we store rank-4
01810                                    // tensors. determine which
01811                                    // position the present entry is
01812                                    // stored in
01813         {
01814           unsigned int base_index[2] ;
01815           if ((indices[0] == 0) && (indices[1] == 0))
01816             base_index[0] = 0;
01817           else if ((indices[0] == 1) && (indices[1] == 1))
01818             base_index[0] = 1;
01819           else if ((indices[0] == 2) && (indices[1] == 2))
01820             base_index[0] = 2;
01821           else if (((indices[0] == 0) && (indices[1] == 1)) ||
01822                    ((indices[0] == 1) && (indices[1] == 0)))
01823             base_index[0] = 3;
01824           else if (((indices[0] == 0) && (indices[1] == 2)) ||
01825                    ((indices[0] == 2) && (indices[1] == 0)))
01826             base_index[0] = 4;
01827           else
01828             {
01829               Assert (((indices[0] == 1) && (indices[1] == 2)) ||
01830                       ((indices[0] == 2) && (indices[1] == 1)),
01831                       ExcInternalError());
01832               base_index[0] = 5;
01833             }
01834 
01835           if ((indices[2] == 0) && (indices[3] == 0))
01836             base_index[1] = 0;
01837           else if ((indices[2] == 1) && (indices[3] == 1))
01838             base_index[1] = 1;
01839           else if ((indices[2] == 2) && (indices[3] == 2))
01840             base_index[1] = 2;
01841           else if (((indices[2] == 0) && (indices[3] == 1)) ||
01842                    ((indices[2] == 1) && (indices[3] == 0)))
01843             base_index[1] = 3;
01844           else if (((indices[2] == 0) && (indices[3] == 2)) ||
01845                    ((indices[2] == 2) && (indices[3] == 0)))
01846             base_index[1] = 4;
01847           else
01848             {
01849               Assert (((indices[2] == 1) && (indices[3] == 2)) ||
01850                       ((indices[2] == 2) && (indices[3] == 1)),
01851                       ExcInternalError());
01852               base_index[1] = 5;
01853             }
01854 
01855           return data[base_index[0]][base_index[1]];
01856         }
01857 
01858         default:
01859               Assert (false, ExcNotImplemented());
01860       }
01861 
01862     static Number dummy;
01863     return dummy;
01864   }
01865 
01866 } // end of namespace internal
01867 
01868 
01869 
01870 template <int rank, int dim, typename Number>
01871 inline
01872 Number &
01873 SymmetricTensor<rank,dim,Number>::operator () (const TableIndices<rank> &indices)
01874 {
01875   for (unsigned int r=0; r<rank; ++r)
01876     Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
01877   return internal::symmetric_tensor_access<dim,Number> (indices, data);
01878 }
01879 
01880 
01881 
01882 template <int rank, int dim, typename Number>
01883 inline
01884 Number
01885 SymmetricTensor<rank,dim,Number>::operator ()
01886   (const TableIndices<rank> &indices) const
01887 {
01888   for (unsigned int r=0; r<rank; ++r)
01889     Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
01890   return internal::symmetric_tensor_access<dim,Number> (indices, data);
01891 }
01892 
01893 
01894 
01895 template <int rank, int dim, typename Number>
01896 internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
01897 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row) const
01898 {
01899   return
01900     internal::SymmetricTensorAccessors::
01901     Accessor<rank,dim,true,rank-1,Number> (*this, TableIndices<rank> (row));
01902 }
01903 
01904 
01905 
01906 template <int rank, int dim, typename Number>
01907 internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
01908 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row)
01909 {
01910   return
01911     internal::SymmetricTensorAccessors::
01912     Accessor<rank,dim,false,rank-1,Number> (*this, TableIndices<rank> (row));
01913 }
01914 
01915 
01916 
01917 template <int rank, int dim, typename Number>
01918 inline
01919 Number
01920 SymmetricTensor<rank,dim,Number>::operator [] (const TableIndices<rank> &indices) const
01921 {
01922   return data[component_to_unrolled_index(indices)];
01923 }
01924 
01925 
01926 
01927 template <int rank, int dim, typename Number>
01928 inline
01929 Number &
01930 SymmetricTensor<rank,dim,Number>::operator [] (const TableIndices<rank> &indices)
01931 {
01932   return data[component_to_unrolled_index(indices)];
01933 }
01934 
01935 
01936 
01937 template <int rank, int dim, typename Number>
01938 inline
01939 Number
01940 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index) const
01941 {
01942   AssertIndexRange (index, data.dimension);
01943   return data[index];
01944 }
01945 
01946 
01947 
01948 template <int rank, int dim, typename Number>
01949 inline
01950 Number &
01951 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index)
01952 {
01953   AssertIndexRange (index, data.dimension);
01954   return data[index];
01955 }
01956 
01957 
01958 
01959 namespace internal
01960 {
01961   template <int dim, typename Number>
01962   inline
01963   Number
01964   compute_norm (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data)
01965   {
01966     Number return_value;
01967     switch (dim)
01968       {
01969       case 1:
01970         return_value = std::fabs(data[0]);
01971         break;
01972       case 2:
01973         return_value = std::sqrt(data[0]*data[0] + data[1]*data[1] +
01974                                  2*data[2]*data[2]);
01975         break;
01976       case 3:
01977         return_value =  std::sqrt(data[0]*data[0] + data[1]*data[1] +
01978                                   data[2]*data[2] + 2*data[3]*data[3] +
01979                                   2*data[4]*data[4] + 2*data[5]*data[5]);
01980         break;
01981       default:
01982         return_value = 0;
01983         for (unsigned int d=0; d<dim; ++d)
01984           return_value += data[d] * data[d];
01985         for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
01986           return_value += 2 * data[d] * data[d];
01987         return_value = std::sqrt(return_value);
01988       }
01989     return return_value;
01990   }
01991 
01992 
01993 
01994   template <int dim, typename Number>
01995   inline
01996   Number
01997   compute_norm (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data)
01998   {
01999     Number return_value;
02000     const unsigned int n_independent_components = data.dimension;
02001 
02002     switch (dim)
02003       {
02004       case 1:
02005         return_value = std::fabs (data[0][0]);
02006         break;
02007       default:
02008         return_value = 0;
02009         for (unsigned int i=0; i<dim; ++i)
02010           for (unsigned int j=0; j<dim; ++j)
02011             return_value += data[i][j] * data[i][j];
02012         for (unsigned int i=0; i<dim; ++i)
02013           for (unsigned int j=dim; j<n_independent_components; ++j)
02014             return_value += 2 * data[i][j] * data[i][j];
02015         for (unsigned int i=dim; i<n_independent_components; ++i)
02016           for (unsigned int j=0; j<dim; ++j)
02017             return_value += 2 * data[i][j] * data[i][j];
02018         for (unsigned int i=dim; i<n_independent_components; ++i)
02019           for (unsigned int j=dim; j<n_independent_components; ++j)
02020             return_value += 4 * data[i][j] * data[i][j];
02021         return_value = std::sqrt(return_value);
02022       }
02023 
02024     return return_value;
02025   }
02026 
02027 } // end of namespace internal
02028 
02029 
02030 
02031 template <int rank, int dim, typename Number>
02032 inline
02033 Number
02034 SymmetricTensor<rank,dim,Number>::norm () const
02035 {
02036   return internal::compute_norm<dim,Number> (data);
02037 }
02038 
02039 
02040 
02041 template <int rank, int dim, typename Number>
02042 inline
02043 unsigned int
02044 SymmetricTensor<rank,dim,Number>::component_to_unrolled_index
02045   (const TableIndices<rank> &indices)
02046 {
02047   Assert (rank == 2, ExcNotImplemented());
02048   Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
02049   Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
02050 
02051   switch(dim)
02052     {
02053     case 1:
02054       return 0;
02055     case 2:
02056       {
02057         static const unsigned int table[2][2] = {{0, 2},
02058                                                  {2, 1}};
02059         return table[indices[0]][indices[1]];
02060       }
02061     case 3:
02062       {
02063         static const unsigned int table[3][3] = {{0, 3, 4},
02064                                                  {3, 1, 5},
02065                                                  {4, 5, 2}};
02066         return table[indices[0]][indices[1]];
02067       }
02068     case 4:
02069       {
02070         static const unsigned int table[4][4] = {{0, 4, 5, 6},
02071                                                  {4, 1, 7, 8},
02072                                                  {5, 7, 2, 9},
02073                                                  {6, 8, 9, 3}};
02074         return table[indices[0]][indices[1]];
02075       }
02076     default:
02077       Assert (false, ExcNotImplemented());
02078       return 0;
02079     }
02080 }
02081 
02082 
02083 
02084 template <int rank, int dim, typename Number>
02085 inline
02086 TableIndices<rank>
02087 SymmetricTensor<rank,dim,Number>::unrolled_to_component_indices
02088   (const unsigned int i)
02089 {
02090   Assert (rank == 2, ExcNotImplemented());
02091   Assert (i < n_independent_components, ExcIndexRange(i, 0, n_independent_components));
02092   switch (dim)
02093     {
02094     case 1:
02095       return TableIndices<2>(0,0);
02096     case 2:
02097       {
02098         static const TableIndices<2> table[3] =
02099           { TableIndices<2> (0,0),
02100             TableIndices<2> (1,1),
02101             TableIndices<2> (0,1) };
02102         return table[i];
02103       }
02104     case 3:
02105       {
02106         static const TableIndices<2> table[6] =
02107           { TableIndices<2> (0,0),
02108             TableIndices<2> (1,1),
02109             TableIndices<2> (2,2),
02110             TableIndices<2> (0,1),
02111             TableIndices<2> (0,2),
02112             TableIndices<2> (1,2) };
02113         return table[i];
02114       }
02115     default:
02116       Assert (false, ExcNotImplemented());
02117       return TableIndices<2>(0,0);
02118     }
02119 }
02120 
02121 
02122 
02123 template <int rank, int dim, typename Number>
02124 template <class Archive>
02125 inline
02126 void
02127 SymmetricTensor<rank,dim,Number>::serialize(Archive & ar, const unsigned int)
02128 {
02129   ar & data;
02130 }
02131 
02132 
02133 #endif // DOXYGEN
02134 
02135 /* ----------------- Non-member functions operating on tensors. ------------ */
02136 
02150 template <int dim, typename Number>
02151 inline
02152 Number determinant (const SymmetricTensor<2,dim,Number> &t)
02153 {
02154   switch (dim)
02155     {
02156     case 1:
02157       return t.data[0];
02158     case 2:
02159       return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
02160     case 3:
02161                                    // in analogy to general tensors, but
02162                                    // there's something to be simplified for
02163                                    // the present case
02164       return ( t.data[0]*t.data[1]*t.data[2]
02165                -t.data[0]*t.data[5]*t.data[5]
02166                -t.data[1]*t.data[4]*t.data[4]
02167                -t.data[2]*t.data[3]*t.data[3]
02168                +2*t.data[3]*t.data[4]*t.data[5] );
02169     default:
02170       Assert (false, ExcNotImplemented());
02171       return 0;
02172     }
02173 }
02174 
02175 
02176 
02186 template <int dim, typename Number>
02187 inline
02188 double third_invariant (const SymmetricTensor<2,dim,Number> &t)
02189 {
02190   return determinant (t);
02191 }
02192 
02193 
02194 
02203 template <int dim, typename Number>
02204 Number trace (const SymmetricTensor<2,dim,Number> &d)
02205 {
02206   Number t=0;
02207   for (unsigned int i=0; i<dim; ++i)
02208     t += d.data[i];
02209   return t;
02210 }
02211 
02212 
02222 template <int dim, typename Number>
02223 inline
02224 Number first_invariant (const SymmetricTensor<2,dim,Number> &t)
02225 {
02226   return trace (t);
02227 }
02228 
02229 
02237 template <typename Number>
02238 inline
02239 Number second_invariant (const SymmetricTensor<2,1,Number> &)
02240 {
02241   return 0;
02242 }
02243 
02244 
02245 
02253 template <typename Number>
02254 inline
02255 Number second_invariant (const SymmetricTensor<2,2,Number> &t)
02256 {
02257   return t[0][0]*t[1][1] - t[0][1]*t[0][1];
02258 }
02259 
02260 
02261 
02269 template <typename Number>
02270 inline
02271 Number second_invariant (const SymmetricTensor<2,3,Number> &t)
02272 {
02273   return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
02274           - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
02275 }
02276 
02277 
02278 
02279 
02289 template <int rank, int dim, typename Number>
02290 inline
02291 SymmetricTensor<rank,dim,Number>
02292 transpose (const SymmetricTensor<rank,dim,Number> &t)
02293 {
02294   return t;
02295 }
02296 
02297 
02298 
02308 template <int dim, typename Number>
02309 inline
02310 SymmetricTensor<2,dim,Number>
02311 deviator (const SymmetricTensor<2,dim,Number> &t)
02312 {
02313   SymmetricTensor<2,dim,Number> tmp = t;
02314 
02315                                    // subtract scaled trace from the diagonal
02316   const Number tr = trace(t) / dim;
02317   for (unsigned int i=0; i<dim; ++i)
02318     tmp.data[i] -= tr;
02319 
02320   return tmp;
02321 }
02322 
02323 
02324 
02331 template <int dim, typename Number>
02332 inline
02333 SymmetricTensor<2,dim,Number>
02334 unit_symmetric_tensor ()
02335 {
02336   SymmetricTensor<2,dim,Number> tmp;
02337   switch (dim)
02338     {
02339     case 1:
02340       tmp.data[0] = 1;
02341       break;
02342     case 2:
02343       tmp.data[0] = tmp.data[1] = 1;
02344       break;
02345     case 3:
02346       tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
02347       break;
02348     default:
02349       for (unsigned int d=0; d<dim; ++d)
02350         tmp.data[d] = 1;
02351     }
02352   return tmp;
02353 }
02354 
02355 
02356 
02363 template <int dim>
02364 inline
02365 SymmetricTensor<2,dim>
02366 unit_symmetric_tensor ()
02367 {
02368   return unit_symmetric_tensor<dim,double>();
02369 }
02370 
02371 
02372 
02387 template <int dim, typename Number>
02388 inline
02389 SymmetricTensor<4,dim,Number>
02390 deviator_tensor ()
02391 {
02392   SymmetricTensor<4,dim,Number> tmp;
02393 
02394                                    // fill the elements treating the diagonal
02395   for (unsigned int i=0; i<dim; ++i)
02396     for (unsigned int j=0; j<dim; ++j)
02397       tmp.data[i][j] = (i==j ? 1 : 0) - 1./dim;
02398 
02399                                    // then fill the ones that copy over the
02400                                    // non-diagonal elements. note that during
02401                                    // the double-contraction, we handle the
02402                                    // off-diagonal elements twice, so simply
02403                                    // copying requires a weight of 1/2
02404   for (unsigned int i=dim;
02405        i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
02406        ++i)
02407     tmp.data[i][i] = 0.5;
02408 
02409   return tmp;
02410 }
02411 
02412 
02413 
02428 template <int dim>
02429 inline
02430 SymmetricTensor<4,dim>
02431 deviator_tensor ()
02432 {
02433   return deviator_tensor<dim,double>();
02434 }
02435 
02436 
02437 
02460 template <int dim, typename Number>
02461 inline
02462 SymmetricTensor<4,dim,Number>
02463 identity_tensor ()
02464 {
02465   SymmetricTensor<4,dim,Number> tmp;
02466 
02467                                    // fill the elements treating the diagonal
02468   for (unsigned int i=0; i<dim; ++i)
02469     tmp.data[i][i] = 1;
02470 
02471                                    // then fill the ones that copy over the
02472                                    // non-diagonal elements. note that during
02473                                    // the double-contraction, we handle the
02474                                    // off-diagonal elements twice, so simply
02475                                    // copying requires a weight of 1/2
02476   for (unsigned int i=dim;
02477        i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
02478        ++i)
02479     tmp.data[i][i] = 0.5;
02480 
02481   return tmp;
02482 }
02483 
02484 
02485 
02507 template <int dim>
02508 inline
02509 SymmetricTensor<4,dim>
02510 identity_tensor ()
02511 {
02512   return identity_tensor<dim,double>();
02513 }
02514 
02515 
02516 
02530 template <int dim, typename Number>
02531 inline
02532 SymmetricTensor<4,dim,Number>
02533 invert (const SymmetricTensor<4,dim,Number> &t)
02534 {
02535   SymmetricTensor<4,dim,Number> tmp;
02536   switch (dim)
02537     {
02538     case 1:
02539       tmp.data[0][0] = 1./t.data[0][0];
02540       break;
02541     case 2:
02542 
02543                                    // inverting this tensor is a little more
02544                                    // complicated than necessary, since we
02545                                    // store the data of 't' as a 3x3 matrix
02546                                    // t.data, but the product between a rank-4
02547                                    // and a rank-2 tensor is really not the
02548                                    // product between this matrix and the
02549                                    // 3-vector of a rhs, but rather
02550                                    //
02551                                    // B.vec = t.data * mult * A.vec
02552                                    //
02553                                    // where mult is a 3x3 matrix with
02554                                    // entries [[1,0,0],[0,1,0],[0,0,2]] to
02555                                    // capture the fact that we need to add up
02556                                    // both the c_ij12*a_12 and the c_ij21*a_21
02557                                    // terms
02558                                    //
02559                                    // in addition, in this scheme, the
02560                                    // identity tensor has the matrix
02561                                    // representation mult^-1.
02562                                    //
02563                                    // the inverse of 't' therefore has the
02564                                    // matrix representation
02565                                    //
02566                                    // inv.data = mult^-1 * t.data^-1 * mult^-1
02567                                    //
02568                                    // in order to compute it, let's first
02569                                    // compute the inverse of t.data and put it
02570                                    // into tmp.data; at the end of the
02571                                    // function we then scale the last row and
02572                                    // column of the inverse by 1/2,
02573                                    // corresponding to the left and right
02574                                    // multiplication with mult^-1
02575       {
02576         const Number t4 = t.data[0][0]*t.data[1][1],
02577                      t6 = t.data[0][0]*t.data[1][2],
02578                      t8 = t.data[0][1]*t.data[1][0],
02579                      t00 = t.data[0][2]*t.data[1][0],
02580                      t01 = t.data[0][1]*t.data[2][0],
02581                      t04 = t.data[0][2]*t.data[2][0],
02582                      t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-
02583                                 t8*t.data[2][2]+t00*t.data[2][1]+
02584                                 t01*t.data[1][2]-t04*t.data[1][1]);
02585         tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
02586         tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
02587         tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
02588         tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
02589         tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
02590         tmp.data[1][2] = -(t6-t00)*t07;
02591         tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
02592         tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
02593         tmp.data[2][2] = (t4-t8)*t07;
02594 
02595                                    // scale last row and column as mentioned
02596                                    // above
02597         tmp.data[2][0] /= 2;
02598         tmp.data[2][1] /= 2;
02599         tmp.data[0][2] /= 2;
02600         tmp.data[1][2] /= 2;
02601         tmp.data[2][2] /= 4;
02602       }
02603       break;
02604     default:
02605       Assert (false, ExcNotImplemented());
02606     }
02607   return tmp;
02608 }
02609 
02610 
02611 
02625 template <>
02626 SymmetricTensor<4,3,double>
02627 invert (const SymmetricTensor<4,3,double> &t);
02628 // this function is implemented in the .cc file for double data types
02629 
02630 
02631 
02646 template <int dim, typename Number>
02647 inline
02648 SymmetricTensor<4,dim,Number>
02649 outer_product (const SymmetricTensor<2,dim,Number> &t1,
02650                const SymmetricTensor<2,dim,Number> &t2)
02651 {
02652   SymmetricTensor<4,dim,Number> tmp;
02653 
02654                                    // fill only the elements really needed
02655   for (unsigned int i=0; i<dim; ++i)
02656     for (unsigned int j=i; j<dim; ++j)
02657       for (unsigned int k=0; k<dim; ++k)
02658         for (unsigned int l=k; l<dim; ++l)
02659           tmp[i][j][k][l] = t1[i][j] * t2[k][l];
02660 
02661   return tmp;
02662 }
02663 
02664 
02665 
02674 template <typename Number>
02675 inline
02676 SymmetricTensor<2,1,Number>
02677 symmetrize (const Tensor<2,1,Number> &t)
02678 {
02679   const Number array[1]
02680     = { t[0][0] };
02681   return SymmetricTensor<2,1,Number>(array);
02682 }
02683 
02684 
02685 
02694 template <typename Number>
02695 inline
02696 SymmetricTensor<2,2,Number>
02697 symmetrize (const Tensor<2,2,Number> &t)
02698 {
02699   const Number array[3]
02700     = { t[0][0], t[1][1], (t[0][1] + t[1][0])/2 };
02701   return SymmetricTensor<2,2,Number>(array);
02702 }
02703 
02704 
02705 
02714 template <typename Number>
02715 inline
02716 SymmetricTensor<2,3,Number>
02717 symmetrize (const Tensor<2,3,Number> &t)
02718 {
02719   const Number array[6]
02720     = { t[0][0], t[1][1], t[2][2],
02721         (t[0][1] + t[1][0])/2,
02722         (t[0][2] + t[2][0])/2,
02723         (t[1][2] + t[2][1])/2 };
02724   return SymmetricTensor<2,3,Number>(array);
02725 }
02726 
02727 
02728 
02735 template <int rank, int dim, typename Number>
02736 inline
02737 SymmetricTensor<rank,dim,Number>
02738 operator * (const SymmetricTensor<rank,dim,Number> &t,
02739             const Number                            factor)
02740 {
02741   SymmetricTensor<rank,dim,Number> tt = t;
02742   tt *= factor;
02743   return tt;
02744 }
02745 
02746 
02747 
02754 template <int rank, int dim, typename Number>
02755 inline
02756 SymmetricTensor<rank,dim,Number>
02757 operator * (const Number                            factor,
02758             const SymmetricTensor<rank,dim,Number> &t)
02759 {
02760   SymmetricTensor<rank,dim,Number> tt = t;
02761   tt *= factor;
02762   return tt;
02763 }
02764 
02765 
02766 
02772 template <int rank, int dim, typename Number>
02773 inline
02774 SymmetricTensor<rank,dim,Number>
02775 operator / (const SymmetricTensor<rank,dim,Number> &t,
02776             const Number                            factor)
02777 {
02778   SymmetricTensor<rank,dim,Number> tt = t;
02779   tt /= factor;
02780   return tt;
02781 }
02782 
02783 
02784 
02791 template <int rank, int dim>
02792 inline
02793 SymmetricTensor<rank,dim>
02794 operator * (const SymmetricTensor<rank,dim> &t,
02795             const double                     factor)
02796 {
02797   SymmetricTensor<rank,dim> tt = t;
02798   tt *= factor;
02799   return tt;
02800 }
02801 
02802 
02803 
02810 template <int rank, int dim>
02811 inline
02812 SymmetricTensor<rank,dim>
02813 operator * (const double                     factor,
02814             const SymmetricTensor<rank,dim> &t)
02815 {
02816   SymmetricTensor<rank,dim> tt = t;
02817   tt *= factor;
02818   return tt;
02819 }
02820 
02821 
02822 
02828 template <int rank, int dim>
02829 inline
02830 SymmetricTensor<rank,dim>
02831 operator / (const SymmetricTensor<rank,dim> &t,
02832             const double                     factor)
02833 {
02834   SymmetricTensor<rank,dim> tt = t;
02835   tt /= factor;
02836   return tt;
02837 }
02838 
02839 
02840 
02841 
02857 template <typename Number>
02858 inline
02859 void
02860 double_contract (SymmetricTensor<2,1,Number> &tmp,
02861                  const SymmetricTensor<4,1,Number> &t,
02862                  const SymmetricTensor<2,1,Number> &s)
02863 {
02864   tmp[0][0] = t[0][0][0][0] * s[0][0];
02865 }
02866 
02867 
02868 
02884 template <typename Number>
02885 inline
02886 void
02887 double_contract (SymmetricTensor<2,1,Number> &tmp,
02888                  const SymmetricTensor<2,1,Number> &s,
02889                  const SymmetricTensor<4,1,Number> &t)
02890 {
02891   tmp[0][0] = t[0][0][0][0] * s[0][0];
02892 }
02893 
02894 
02895 
02910 template <typename Number>
02911 inline
02912 void
02913 double_contract (SymmetricTensor<2,2,Number> &tmp,
02914                  const SymmetricTensor<4,2,Number> &t,
02915                  const SymmetricTensor<2,2,Number> &s)
02916 {
02917   const unsigned int dim = 2;
02918 
02919   for (unsigned int i=0; i<dim; ++i)
02920     for (unsigned int j=i; j<dim; ++j)
02921       tmp[i][j] = t[i][j][0][0] * s[0][0] +
02922                   t[i][j][1][1] * s[1][1] +
02923                   2 * t[i][j][0][1] * s[0][1];
02924 }
02925 
02926 
02927 
02943 template <typename Number>
02944 inline
02945 void
02946 double_contract (SymmetricTensor<2,2,Number> &tmp,
02947                  const SymmetricTensor<2,2,Number> &s,
02948                  const SymmetricTensor<4,2,Number> &t)
02949 {
02950   const unsigned int dim = 2;
02951 
02952   for (unsigned int i=0; i<dim; ++i)
02953     for (unsigned int j=i; j<dim; ++j)
02954       tmp[i][j] = s[0][0] * t[0][0][i][j] * +
02955                   s[1][1] * t[1][1][i][j] +
02956                   2 * s[0][1] * t[0][1][i][j];
02957 }
02958 
02959 
02960 
02976 template <typename Number>
02977 inline
02978 void
02979 double_contract (SymmetricTensor<2,3,Number> &tmp,
02980                  const SymmetricTensor<4,3,Number> &t,
02981                  const SymmetricTensor<2,3,Number> &s)
02982 {
02983   const unsigned int dim = 3;
02984 
02985   for (unsigned int i=0; i<dim; ++i)
02986     for (unsigned int j=i; j<dim; ++j)
02987       tmp[i][j] = t[i][j][0][0] * s[0][0] +
02988                   t[i][j][1][1] * s[1][1] +
02989                   t[i][j][2][2] * s[2][2] +
02990                   2 * t[i][j][0][1] * s[0][1] +
02991                   2 * t[i][j][0][2] * s[0][2] +
02992                   2 * t[i][j][1][2] * s[1][2];
02993 }
02994 
02995 
02996 
03012 template <typename Number>
03013 inline
03014 void
03015 double_contract (SymmetricTensor<2,3,Number> &tmp,
03016                  const SymmetricTensor<2,3,Number> &s,
03017                  const SymmetricTensor<4,3,Number> &t)
03018 {
03019   const unsigned int dim = 3;
03020 
03021   for (unsigned int i=0; i<dim; ++i)
03022     for (unsigned int j=i; j<dim; ++j)
03023       tmp[i][j] = s[0][0] * t[0][0][i][j] +
03024                   s[1][1] * t[1][1][i][j] +
03025                   s[2][2] * t[2][2][i][j] +
03026                   2 * s[0][1] * t[0][1][i][j] +
03027                   2 * s[0][2] * t[0][2][i][j] +
03028                   2 * s[1][2] * t[1][2][i][j];
03029 }
03030 
03031 
03032 
03049 template <int dim, typename Number>
03050 Tensor<1,dim,Number>
03051 operator * (const SymmetricTensor<2,dim,Number> &src1,
03052             const Tensor<1,dim,Number> &src2)
03053 {
03054   Tensor<1,dim,Number> dest;
03055   for (unsigned int i=0; i<dim; ++i)
03056     for (unsigned int j=0; j<dim; ++j)
03057       dest[i] += src1[i][j] * src2[j];
03058   return dest;
03059 }
03060 
03061 
03071 template <int dim, typename Number>
03072 inline
03073 std::ostream & operator << (std::ostream &out,
03074                             const SymmetricTensor<2,dim,Number> &t)
03075 {
03076                                    //make out lives a bit simpler by outputing
03077                                    //the tensor through the operator for the
03078                                    //general Tensor class
03079   Tensor<2,dim,Number> tt;
03080 
03081   for (unsigned int i=0; i<dim; ++i)
03082     for (unsigned int j=0; j<dim; ++j)
03083       tt[i][j] = t[i][j];
03084 
03085   return out << tt;
03086 }
03087 
03088 
03089 
03099 template <int dim, typename Number>
03100 inline
03101 std::ostream & operator << (std::ostream &out,
03102                             const SymmetricTensor<4,dim,Number> &t)
03103 {
03104                                    //make out lives a bit simpler by outputing
03105                                    //the tensor through the operator for the
03106                                    //general Tensor class
03107   Tensor<4,dim,Number> tt;
03108 
03109   for (unsigned int i=0; i<dim; ++i)
03110     for (unsigned int j=0; j<dim; ++j)
03111       for (unsigned int k=0; k<dim; ++k)
03112         for (unsigned int l=0; l<dim; ++l)
03113           tt[i][j][k][l] = t[i][j][k][l];
03114 
03115   return out << tt;
03116 }
03117 
03118 
03119 DEAL_II_NAMESPACE_CLOSE
03120 
03121 #endif
03122 
 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