00001
00002
00003
00004
00005
00006
00007
00008
00009
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
00407
00408
00409
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
00526
00527
00528
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
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
01191
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 }
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
01514
01515
01516
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
01536
01537
01538
01539
01540
01541
01542
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
01558
01559
01560 if (indices[0] == indices[1])
01561 return data[indices[0]];
01562
01563
01564
01565
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
01573
01574
01575 if (indices[0] == indices[1])
01576 return data[indices[0]];
01577
01578
01579
01580
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
01615
01616
01617 if (indices[0] == indices[1])
01618 return data[indices[0]];
01619
01620
01621
01622
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
01630
01631
01632 if (indices[0] == indices[1])
01633 return data[indices[0]];
01634
01635
01636
01637
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
01672
01673
01674
01675
01676
01677
01678
01679
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
01701
01702
01703
01704
01705
01706
01707
01708
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
01776
01777
01778
01779
01780
01781
01782
01783
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
01805
01806
01807
01808
01809
01810
01811
01812
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 }
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 }
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
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
02162
02163
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
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
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
02400
02401
02402
02403
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
02468 for (unsigned int i=0; i<dim; ++i)
02469 tmp.data[i][i] = 1;
02470
02471
02472
02473
02474
02475
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
02544
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558
02559
02560
02561
02562
02563
02564
02565
02566
02567
02568
02569
02570
02571
02572
02573
02574
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
02596
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
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
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
03077
03078
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
03105
03106
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