00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__tensor_base_h
00013 #define __deal2__tensor_base_h
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include <deal.II/base/config.h>
00023 #include <deal.II/base/exceptions.h>
00024 #include <vector>
00025
00026 #include <cmath>
00027
00028 DEAL_II_NAMESPACE_OPEN
00029
00030
00031 #ifdef HAVE_STD_OSTREAM_HEADER
00032 # include <ostream>
00033 #else
00034 # include <iostream>
00035 #endif
00036
00037 template <typename number> class Vector;
00038
00039
00040
00041
00042
00043
00044 template <int dim, typename Number=double> class Point;
00045
00046
00047
00048 template <int rank_, int dim, typename Number=double> class Tensor;
00049 template <int dim, typename Number> class Tensor<0,dim,Number>;
00050 template <int dim, typename Number> class Tensor<1,dim,Number>;
00051
00052
00053
00069 template <int dim, typename Number>
00070 class Tensor<0,dim,Number>
00071 {
00072 public:
00090 static const unsigned int dimension = dim;
00091
00096 static const unsigned int rank = 0;
00097
00103 typedef Number value_type;
00104
00119 typedef typename numbers::NumberTraits<Number>::real_type real_type;
00120
00124 Tensor ();
00125
00131 Tensor (const value_type &initializer);
00132
00136 Tensor (const Tensor<0,dim,Number> &);
00137
00143 operator Number () const;
00144
00154 operator Number& ();
00155
00159 Tensor<0,dim,Number> & operator = (const Tensor<0,dim,Number> &);
00160
00164 Tensor<0,dim,Number> & operator = (const Number d);
00165
00170 bool operator == (const Tensor<0,dim,Number> &) const;
00171
00176 bool operator != (const Tensor<0,dim,Number> &) const;
00177
00183 Tensor<0,dim,Number> & operator += (const Tensor<0,dim,Number> &);
00184
00188 Tensor<0,dim,Number> & operator -= (const Tensor<0,dim,Number> &);
00189
00195 Tensor<0,dim,Number> & operator *= (const Number factor);
00196
00200 Tensor<0,dim,Number> & operator /= (const Number factor);
00201
00206 Number operator * (const Tensor<0,dim,Number> &) const;
00207
00214 Tensor<0,dim,Number> operator + (const Tensor<0,dim,Number> &) const;
00215
00223 Tensor<0,dim,Number> operator - (const Tensor<0,dim,Number> &) const;
00224
00228 Tensor<0,dim,Number> operator - () const;
00229
00240 real_type norm () const;
00241
00254 real_type norm_square () const;
00255
00274 void clear ();
00275
00286 DeclException1 (ExcDimTooSmall,
00287 int,
00288 << "dim must be positive, but was " << arg1);
00289
00294 template <class Archive>
00295 void serialize(Archive & ar, const unsigned int version);
00296
00297 private:
00301 Number value;
00302 };
00303
00304
00305
00328 template <int dim,typename Number>
00329 class Tensor<1,dim,Number>
00330 {
00331 public:
00349 static const unsigned int dimension = dim;
00350
00355 static const unsigned int rank = 1;
00356
00362 typedef Number value_type;
00363
00378 typedef typename numbers::NumberTraits<Number>::real_type real_type;
00379
00391 typedef Number array_type[(dim!=0) ? dim : 100000000];
00392
00398 explicit Tensor (const bool initialize = true);
00399
00405 Tensor (const array_type &initializer);
00406
00410 Tensor (const Tensor<1,dim,Number> &);
00411
00422 Number operator [] (const unsigned int index) const;
00423
00434 Number & operator [] (const unsigned int index);
00435
00439 Tensor<1,dim,Number> & operator = (const Tensor<1,dim,Number> &);
00440
00453 Tensor<1,dim,Number> & operator = (const Number d);
00454
00459 bool operator == (const Tensor<1,dim,Number> &) const;
00460
00465 bool operator != (const Tensor<1,dim,Number> &) const;
00466
00472 Tensor<1,dim,Number> & operator += (const Tensor<1,dim,Number> &);
00473
00477 Tensor<1,dim,Number> & operator -= (const Tensor<1,dim,Number> &);
00478
00484 Tensor<1,dim,Number> & operator *= (const Number factor);
00485
00489 Tensor<1,dim,Number> & operator /= (const Number factor);
00490
00495 Number operator * (const Tensor<1,dim,Number> &) const;
00496
00503 Tensor<1,dim,Number> operator + (const Tensor<1,dim,Number> &) const;
00504
00512 Tensor<1,dim,Number> operator - (const Tensor<1,dim,Number> &) const;
00513
00517 Tensor<1,dim,Number> operator - () const;
00518
00529 real_type norm () const;
00530
00543 real_type norm_square () const;
00544
00563 void clear ();
00564
00574 template <typename Number2>
00575 void unroll (Vector<Number2> &result) const;
00576
00583 static std::size_t memory_consumption ();
00584
00595 DeclException1 (ExcDimTooSmall,
00596 int,
00597 << "dim must be positive, but was " << arg1);
00598
00603 template <class Archive>
00604 void serialize(Archive & ar, const unsigned int version);
00605
00606 private:
00616 Number values[(dim!=0) ? (dim) : 1];
00617
00618 #ifdef DEAL_II_TEMPLATE_SPEC_ACCESS_WORKAROUND
00619 public:
00620 #endif
00621
00632 template <typename Number2>
00633 void unroll_recursion (Vector<Number2> &result,
00634 unsigned int &start_index) const;
00635
00636 private:
00650 template <int otherrank, int otherdim, typename OtherNumber> friend class ::Tensor;
00651
00657 friend class Point<dim,Number>;
00658 };
00659
00660
00664 template <int dim,typename Number>
00665 std::ostream & operator << (std::ostream &out, const Tensor<0,dim,Number> &p);
00666
00671 template <int dim,typename Number>
00672 std::ostream & operator << (std::ostream &out, const Tensor<1,dim,Number> &p);
00673
00674
00675 #ifndef DOXYGEN
00676
00677
00678
00679 template <int dim,typename Number>
00680 inline
00681 Tensor<0,dim,Number>::Tensor ()
00682 {
00683 Assert (dim>0, ExcDimTooSmall(dim));
00684
00685 value = 0;
00686 }
00687
00688
00689
00690 template <int dim, typename Number>
00691 inline
00692 Tensor<0,dim,Number>::Tensor (const value_type &initializer)
00693 {
00694 Assert (dim>0, ExcDimTooSmall(dim));
00695
00696 value = initializer;
00697 }
00698
00699
00700
00701 template <int dim, typename Number>
00702 inline
00703 Tensor<0,dim,Number>::Tensor (const Tensor<0,dim,Number> &p)
00704 {
00705 Assert (dim>0, ExcDimTooSmall(dim));
00706
00707 value = p.value;
00708 }
00709
00710
00711
00712
00713 template <int dim, typename Number>
00714 inline
00715 Tensor<0,dim,Number>::operator Number () const
00716 {
00717 return value;
00718 }
00719
00720
00721
00722 template <int dim, typename Number>
00723 inline
00724 Tensor<0,dim,Number>::operator Number & ()
00725 {
00726 return value;
00727 }
00728
00729
00730
00731 template <int dim, typename Number>
00732 inline
00733 Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator = (const Tensor<0,dim,Number> &p)
00734 {
00735 value = p.value;
00736 return *this;
00737 }
00738
00739
00740
00741 template <int dim, typename Number>
00742 inline
00743 Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator = (const Number d)
00744 {
00745 value = d;
00746 return *this;
00747 }
00748
00749
00750
00751 template <int dim, typename Number>
00752 inline
00753 bool Tensor<0,dim,Number>::operator == (const Tensor<0,dim,Number> &p) const
00754 {
00755 return (value == p.value);
00756 }
00757
00758
00759
00760 template <int dim, typename Number>
00761 inline
00762 bool Tensor<0,dim,Number>::operator != (const Tensor<0,dim,Number> &p) const
00763 {
00764 return !((*this) == p);
00765 }
00766
00767
00768
00769 template <int dim, typename Number>
00770 inline
00771 Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator += (const Tensor<0,dim,Number> &p)
00772 {
00773 value += p.value;
00774 return *this;
00775 }
00776
00777
00778
00779 template <int dim, typename Number>
00780 inline
00781 Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator -= (const Tensor<0,dim,Number> &p)
00782 {
00783 value -= p.value;
00784 return *this;
00785 }
00786
00787
00788
00789 template <int dim, typename Number>
00790 inline
00791 Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator *= (const Number s)
00792 {
00793 value *= s;
00794 return *this;
00795 }
00796
00797
00798
00799 template <int dim, typename Number>
00800 inline
00801 Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator /= (const Number s)
00802 {
00803 value /= s;
00804 return *this;
00805 }
00806
00807
00808
00809 template <int dim, typename Number>
00810 inline
00811 Number Tensor<0,dim,Number>::operator * (const Tensor<0,dim,Number> &p) const
00812 {
00813 return value*p.value;
00814 }
00815
00816
00817
00818 template <int dim, typename Number>
00819 inline
00820 Tensor<0,dim,Number> Tensor<0,dim,Number>::operator + (const Tensor<0,dim,Number> &p) const
00821 {
00822 return value+p.value;
00823 }
00824
00825
00826
00827 template <int dim, typename Number>
00828 inline
00829 Tensor<0,dim,Number> Tensor<0,dim,Number>::operator - (const Tensor<0,dim,Number> &p) const
00830 {
00831 return value-p.value;
00832 }
00833
00834
00835
00836 template <int dim, typename Number>
00837 inline
00838 Tensor<0,dim,Number> Tensor<0,dim,Number>::operator - () const
00839 {
00840 return -value;
00841 }
00842
00843
00844
00845 template <int dim, typename Number>
00846 inline
00847 typename Tensor<0,dim,Number>::real_type
00848 Tensor<0,dim,Number>::norm () const
00849 {
00850 return numbers::NumberTraits<Number>::abs (value);
00851 }
00852
00853
00854
00855 template <int dim, typename Number>
00856 inline
00857 typename Tensor<0,dim,Number>::real_type
00858 Tensor<0,dim,Number>::norm_square () const
00859 {
00860 return numbers::NumberTraits<Number>::abs_square (value);
00861 }
00862
00863
00864
00865 template <int dim, typename Number>
00866 inline
00867 void Tensor<0,dim,Number>::clear ()
00868 {
00869 value = 0;
00870 }
00871
00872
00873
00874 template <int dim, typename Number>
00875 template <class Archive>
00876 inline
00877 void Tensor<0,dim,Number>::serialize(Archive & ar, const unsigned int)
00878 {
00879 ar & value;
00880 }
00881
00882
00883
00884
00885 template <int dim, typename Number>
00886 inline
00887 Tensor<1,dim,Number>::Tensor (const bool initialize)
00888 {
00889 if (initialize)
00890
00891
00892
00893
00894
00895 for (unsigned int i=0; i!=dim; ++i)
00896 values[i] = Number();
00897 }
00898
00899
00900
00901 template <int dim, typename Number>
00902 inline
00903 Tensor<1,dim,Number>::Tensor (const array_type &initializer)
00904 {
00905 Assert (dim>0, ExcDimTooSmall(dim));
00906
00907 for (unsigned int i=0; i<dim; ++i)
00908 values[i] = initializer[i];
00909 }
00910
00911
00912
00913 template <int dim, typename Number>
00914 inline
00915 Tensor<1,dim,Number>::Tensor (const Tensor<1,dim,Number> &p)
00916 {
00917 Assert (dim>0, ExcDimTooSmall(dim));
00918
00919 for (unsigned int i=0; i<dim; ++i)
00920 values[i] = p.values[i];
00921 }
00922
00923
00924
00925 template <>
00926 inline
00927 Tensor<1,0,double>::Tensor (const Tensor<1,0,double> &)
00928 {
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939 }
00940
00941
00942
00943 template <int dim, typename Number>
00944 inline
00945 Number Tensor<1,dim,Number>::operator [] (const unsigned int index) const
00946 {
00947 Assert (index<dim, ExcIndexRange (index, 0, dim));
00948 return values[index];
00949 }
00950
00951
00952
00953 template <int dim, typename Number>
00954 inline
00955 Number & Tensor<1,dim,Number>::operator [] (const unsigned int index)
00956 {
00957 Assert (index<dim, ExcIndexRange (index, 0, dim));
00958 return values[index];
00959 }
00960
00961
00962
00963 template <>
00964 inline
00965 Tensor<1,0,double> & Tensor<1,0,double>::operator = (const Tensor<1,0,double> &)
00966 {
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977 return *this;
00978 }
00979
00980
00981
00982 template <int dim, typename Number>
00983 inline
00984 Tensor<1,dim,Number> &
00985 Tensor<1,dim,Number>::operator = (const Tensor<1,dim,Number> &p)
00986 {
00987
00988
00989
00990
00991
00992 switch (dim)
00993 {
00994 case 1:
00995 values[0] = p.values[0];
00996 break;
00997 case 2:
00998 values[0] = p.values[0];
00999 values[1] = p.values[1];
01000 break;
01001 case 3:
01002 values[0] = p.values[0];
01003 values[1] = p.values[1];
01004 values[2] = p.values[2];
01005 break;
01006 default:
01007 for (unsigned int i=0; i<dim; ++i)
01008 values[i] = p.values[i];
01009 }
01010 return *this;
01011 }
01012
01013
01014
01015 template <int dim, typename Number>
01016 inline
01017 Tensor<1,dim,Number> & Tensor<1,dim,Number>::operator = (const Number d)
01018 {
01019 Assert (d==Number(0), ExcMessage ("Only assignment with zero is allowed"));
01020
01021 for (unsigned int i=0; i<dim; ++i)
01022 values[i] = 0;
01023
01024 return *this;
01025 }
01026
01027
01028
01029 template <int dim, typename Number>
01030 inline
01031 bool Tensor<1,dim,Number>::operator == (const Tensor<1,dim,Number> &p) const
01032 {
01033 for (unsigned int i=0; i<dim; ++i)
01034 if (values[i] != p.values[i])
01035 return false;
01036 return true;
01037 }
01038
01039
01040
01041 template <>
01042 inline
01043 bool Tensor<1,0,double>::operator == (const Tensor<1,0,double> &) const
01044 {
01045 return true;
01046 }
01047
01048
01049
01050 template <int dim, typename Number>
01051 inline
01052 bool Tensor<1,dim,Number>::operator != (const Tensor<1,dim,Number> &p) const
01053 {
01054 return !((*this) == p);
01055 }
01056
01057
01058
01059 template <int dim, typename Number>
01060 inline
01061 Tensor<1,dim,Number> & Tensor<1,dim,Number>::operator += (const Tensor<1,dim,Number> &p)
01062 {
01063 for (unsigned int i=0; i<dim; ++i)
01064 values[i] += p.values[i];
01065 return *this;
01066 }
01067
01068
01069
01070 template <int dim, typename Number>
01071 inline
01072 Tensor<1,dim,Number> & Tensor<1,dim,Number>::operator -= (const Tensor<1,dim,Number> &p)
01073 {
01074 for (unsigned int i=0; i<dim; ++i)
01075 values[i] -= p.values[i];
01076 return *this;
01077 }
01078
01079
01080
01081 template <int dim, typename Number>
01082 inline
01083 Tensor<1,dim,Number> & Tensor<1,dim,Number>::operator *= (const Number s)
01084 {
01085 for (unsigned int i=0; i<dim; ++i)
01086 values[i] *= s;
01087 return *this;
01088 }
01089
01090
01091
01092 template <int dim, typename Number>
01093 inline
01094 Tensor<1,dim,Number> & Tensor<1,dim,Number>::operator /= (const Number s)
01095 {
01096 for (unsigned int i=0; i<dim; ++i)
01097 values[i] /= s;
01098 return *this;
01099 }
01100
01101
01102
01103 template <int dim, typename Number>
01104 inline
01105 Number
01106 Tensor<1,dim,Number>::operator * (const Tensor<1,dim,Number> &p) const
01107 {
01108
01109
01110
01111
01112
01113 switch (dim)
01114 {
01115 case 1:
01116 return (values[0] * p.values[0]);
01117 break;
01118 case 2:
01119 return (values[0] * p.values[0] +
01120 values[1] * p.values[1]);
01121 break;
01122 case 3:
01123 return (values[0] * p.values[0] +
01124 values[1] * p.values[1] +
01125 values[2] * p.values[2]);
01126 break;
01127 default:
01128 Number q = values[0] * p.values[0];
01129 for (unsigned int i=1; i<dim; ++i)
01130 q += values[i] * p.values[i];
01131 return q;
01132 }
01133 }
01134
01135
01136
01137 template <int dim, typename Number>
01138 inline
01139 Tensor<1,dim,Number> Tensor<1,dim,Number>::operator + (const Tensor<1,dim,Number> &p) const
01140 {
01141 return (Tensor<1,dim,Number>(*this) += p);
01142 }
01143
01144
01145
01146 template <int dim, typename Number>
01147 inline
01148 Tensor<1,dim,Number> Tensor<1,dim,Number>::operator - (const Tensor<1,dim,Number> &p) const
01149 {
01150 return (Tensor<1,dim,Number>(*this) -= p);
01151 }
01152
01153
01154
01155 template <int dim, typename Number>
01156 inline
01157 Tensor<1,dim,Number> Tensor<1,dim,Number>::operator - () const
01158 {
01159 Tensor<1,dim,Number> result (false);
01160 for (unsigned int i=0; i<dim; ++i)
01161 result.values[i] = -values[i];
01162 return result;
01163 }
01164
01165
01166
01167 template <int dim, typename Number>
01168 inline
01169 typename Tensor<1,dim,Number>::real_type
01170 Tensor<1,dim,Number>::norm () const
01171 {
01172 return std::sqrt (norm_square());
01173 }
01174
01175
01176
01177 template <int dim, typename Number>
01178 inline
01179 typename Tensor<1,dim,Number>::real_type
01180 Tensor<1,dim,Number>::norm_square () const
01181 {
01182 real_type s = numbers::NumberTraits<Number>::abs_square(values[0]);
01183 for (unsigned int i=1; i<dim; ++i)
01184 s += numbers::NumberTraits<Number>::abs_square(values[i]);
01185
01186 return s;
01187 }
01188
01189
01190
01191 template <int dim, typename Number>
01192 template <typename Number2>
01193 inline
01194 void
01195 Tensor<1,dim,Number>::unroll (Vector<Number2> &result) const
01196 {
01197 Assert (result.size()==dim,
01198 ExcDimensionMismatch(dim, result.size()));
01199
01200 unsigned index = 0;
01201 unroll_recursion (result,index);
01202 }
01203
01204
01205
01206 template<int dim, typename Number>
01207 template <typename Number2>
01208 inline
01209 void
01210 Tensor<1,dim,Number>::unroll_recursion (Vector<Number2> &result,
01211 unsigned int &index) const
01212 {
01213 for (unsigned i=0; i<dim; ++i)
01214 result(index++) = operator[](i);
01215 }
01216
01217
01218
01219 template <int dim, typename Number>
01220 inline
01221 void Tensor<1,dim,Number>::clear ()
01222 {
01223 for (unsigned int i=0; i<dim; ++i)
01224 values[i] = 0;
01225 }
01226
01227
01228
01229 template <int dim, typename Number>
01230 inline
01231 std::size_t
01232 Tensor<1,dim,Number>::memory_consumption ()
01233 {
01234 return sizeof(Tensor<1,dim,Number>);
01235 }
01236
01237
01238 template <int dim, typename Number>
01239 template <class Archive>
01240 inline
01241 void Tensor<1,dim,Number>::serialize(Archive & ar, const unsigned int)
01242 {
01243 ar & values;
01244 }
01245 #endif // DOXYGEN
01246
01247
01254 template <int dim, typename Number>
01255 inline
01256 std::ostream & operator << (std::ostream &out, const Tensor<0,dim,Number> &p)
01257 {
01258 out << static_cast<Number>(p);
01259 return out;
01260 }
01261
01262
01263
01270 template <int dim, typename Number>
01271 inline
01272 std::ostream & operator << (std::ostream &out, const Tensor<1,dim,Number> &p)
01273 {
01274 for (unsigned int i=0; i<dim-1; ++i)
01275 out << p[i] << ' ';
01276 out << p[dim-1];
01277
01278 return out;
01279 }
01280
01281
01282
01290 inline
01291 std::ostream & operator << (std::ostream &out, const Tensor<1,1,double> &p)
01292 {
01293 out << p[0];
01294
01295 return out;
01296 }
01297
01298
01299
01305 template <int dim, typename Number>
01306 inline
01307 Tensor<1,dim,Number>
01308 operator * (const Tensor<1,dim,Number> &t,
01309 const Number factor)
01310 {
01311 Tensor<1,dim,Number> tt (false);
01312 for (unsigned int d=0; d<dim; ++d)
01313 tt[d] = t[d] * factor;
01314 return tt;
01315 }
01316
01317
01318
01324 template <int dim, typename Number>
01325 inline
01326 Tensor<1,dim,Number>
01327 operator * (const Number factor,
01328 const Tensor<1,dim,Number> &t)
01329 {
01330 Tensor<1,dim,Number> tt (false);
01331 for (unsigned int d=0; d<dim; ++d)
01332 tt[d] = t[d] * factor;
01333 return tt;
01334 }
01335
01336
01337
01343 template <int dim, typename Number>
01344 inline
01345 Tensor<1,dim,Number>
01346 operator / (const Tensor<1,dim,Number> &t,
01347 const Number factor)
01348 {
01349 Tensor<1,dim,Number> tt (false);
01350 for (unsigned int d=0; d<dim; ++d)
01351 tt[d] = t[d] / factor;
01352 return tt;
01353 }
01354
01355
01356
01362 template <int dim>
01363 inline
01364 Tensor<1,dim>
01365 operator * (const Tensor<1,dim> &t,
01366 const double factor)
01367 {
01368 Tensor<1,dim> tt (false);
01369 for (unsigned int d=0; d<dim; ++d)
01370 tt[d] = t[d] * factor;
01371 return tt;
01372 }
01373
01374
01375
01381 template <int dim>
01382 inline
01383 Tensor<1,dim>
01384 operator * (const double factor,
01385 const Tensor<1,dim> &t)
01386 {
01387 Tensor<1,dim> tt (false);
01388 for (unsigned int d=0; d<dim; ++d)
01389 tt[d] = t[d] * factor;
01390 return tt;
01391 }
01392
01393
01394
01400 template <int dim>
01401 inline
01402 Tensor<1,dim>
01403 operator / (const Tensor<1,dim> &t,
01404 const double factor)
01405 {
01406 Tensor<1,dim> tt (false);
01407 for (unsigned int d=0; d<dim; ++d)
01408 tt[d] = t[d] / factor;
01409 return tt;
01410 }
01411
01412
01413 DEAL_II_NAMESPACE_CLOSE
01414
01415 #endif
01416
01417
01418