dune-localfunctions  2.2.0
tensor.hh
Go to the documentation of this file.
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set et ts=8 sw=2 sts=2:
00003 
00004 #ifndef DUNE_TENSOR_HH
00005 #define DUNE_TENSOR_HH
00006 
00007 #include <ostream>
00008 #include <vector>
00009 
00010 #include <dune/common/fvector.hh>
00011 #include <dune/common/misc.hh>
00012 #include <dune/common/typetraits.hh>
00013 
00014 #include <dune/localfunctions/utility/field.hh>
00015 
00016 namespace Dune 
00017 {
00018   /***********************************************
00019    * The classes here are work in progress.
00020    * Basically they provide tensor structures for
00021    * higher order derivatives of vector valued function.
00022    * Two storage structures are provided 
00023    * (either based on the components of the vector valued
00024    * functions or on the order of the derivative).
00025    * Conversions are supplied between the two storage
00026    * structures and simple operations, which make the
00027    * code difficult to use and requires rewritting...
00028    ***************************************************/
00029 
00030   // Structure for scalar tensor of order deriv
00031   template <class F,int dimD,unsigned int deriv>
00032   class LFETensor
00033   {
00034     typedef LFETensor<F,dimD,deriv> This;
00035     typedef LFETensor<F,dimD-1,deriv> BaseDim;
00036     typedef LFETensor<F,dimD,deriv-1> BaseDeriv;
00037 
00038   public:
00039     typedef F field_type;
00040     static const unsigned int size = BaseDim::size+BaseDeriv::size;
00041     typedef Dune::FieldVector<F,size> Block;
00042 
00043     template< class FF >
00044     This &operator= ( const FF &f )
00045     {
00046       block() = field_cast< F >( f );
00047       return *this;
00048     }
00049 
00050     This &operator= ( const Block &b )
00051     {
00052       block() = b;
00053       return *this;
00054     }
00055 
00056     This &operator*= ( const field_type &f )
00057     {
00058       block() *= f;
00059       return *this;
00060     }
00061 
00062     const field_type &operator[] ( const unsigned int i ) const
00063     {
00064       return block()[ i ];
00065     }
00066 
00067     field_type &operator[] ( const unsigned int i )
00068     {
00069       return block()[ i ];
00070     }
00071 
00072     Block &block() 
00073     {
00074       return block_;
00075     }
00076     const Block &block() const
00077     {
00078       return block_;
00079     }
00080     void axpy(const F& a, const This &y)
00081     {
00082       block().axpy(a,y.block());
00083     }
00084     template <class Fy>
00085     void assign(const LFETensor<Fy,dimD,deriv> &y)
00086     {
00087       field_cast(y.block(),block());
00088     }
00089     Block block_;
00090   };
00091 
00092   // ******************************************
00093   template <class F,unsigned int deriv>
00094   struct LFETensor<F,0,deriv> 
00095   {
00096     static const int size = 0;
00097   };
00098 
00099   template <class F>
00100   struct LFETensor<F,0,0> 
00101   {
00102     static const int size = 1;
00103   };
00104 
00105   template <class F,int dimD>
00106   class LFETensor<F,dimD,0> 
00107   {
00108     typedef LFETensor<F,dimD,0> This;
00109 
00110   public:
00111     typedef F field_type;
00112     static const int size = 1;
00113     typedef Dune::FieldVector<F,size> Block;
00114 
00115     template< class FF >
00116     This &operator= ( const FF &f )
00117     {
00118       block() = field_cast< F >( f );
00119       return *this;
00120     }
00121 
00122     This &operator= ( const Block &b )
00123     {
00124       block() = b;
00125       return *this;
00126     }
00127 
00128     This &operator*= ( const field_type &f )
00129     {
00130       block() *= f;
00131       return *this;
00132     }
00133 
00134     const F &operator[] ( const unsigned int i ) const
00135     {
00136       return block()[ i ];
00137     }
00138 
00139     F &operator[] ( const unsigned int i )
00140     {
00141       return block()[ i ];
00142     }
00143 
00144     void axpy(const F& a, const This &y)
00145     {
00146       block().axpy(a,y.block());
00147     }
00148     template <class Fy>
00149     void assign(const LFETensor<Fy,dimD,0> &y)
00150     {
00151       field_cast(y.block(),block());
00152     }
00153 
00154     Block &block() 
00155     {
00156       return block_;
00157     }
00158     const Block &block() const
00159     {
00160       return block_;
00161     }
00162     Block block_;
00163   };
00164   // ***********************************************************
00165   // Structure for all derivatives up to order deriv
00166   // for vector valued function
00167   enum DerivativeLayout {value,derivative};
00168   template <class F,int dimD,int dimR,unsigned int deriv,
00169             DerivativeLayout layout>
00170   struct Derivatives;
00171 
00172   // Implemnetation for valued based layout
00173   template <class F,int dimD,int dimR,unsigned int deriv>
00174   struct Derivatives<F,dimD,dimR,deriv,value>
00175   : public Derivatives<F,dimD,dimR,deriv-1,value>
00176   {
00177     typedef Derivatives<F,dimD,dimR,deriv,value> This;
00178     typedef Derivatives<F,dimD,dimR,deriv-1,value> Base;
00179     typedef LFETensor<F,dimD,deriv> ThisLFETensor;
00180 
00181     typedef F Field;
00182     typedef F field_type;
00183 
00184     static const DerivativeLayout layout = value;
00185     static const unsigned int dimDomain = dimD;
00186     static const unsigned int dimRange = dimR;
00187     // size needs to be an anonymous enum value for gcc 3.4 compatibility
00188     enum { size = Base::size+ThisLFETensor::size*dimR };
00189     typedef Dune::FieldVector<F,size> Block;
00190  
00191     This &operator=(const F& f)
00192     {
00193       block() = f;
00194       return *this;
00195     }
00196     This &operator=(const Dune::FieldVector<ThisLFETensor,dimR> &t) 
00197     {
00198       tensor_ = t;
00199       return *this;
00200     }
00201     template <unsigned int dorder>
00202     This &operator=(const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &t) 
00203     {
00204       tensor<dorder>() = t;
00205       return *this;
00206     }
00207     This &operator=(const Block &t) 
00208     {
00209       block() = t;
00210       return *this;
00211     }
00212 
00213     This &operator*= ( const field_type &f )
00214     {
00215       block() *= f;
00216       return *this;
00217     }
00218 
00219     void axpy(const F &a, const This &y)
00220     {
00221       block().axpy(a,y.block());
00222     }
00223 
00224     // assign with same layout (only diffrent Field)
00225     template <class Fy>
00226     void assign(const Derivatives<Fy,dimD,dimR,deriv,value> &y) 
00227     {
00228       field_cast(y.block(),block());
00229     }
00230     // assign with diffrent layout (same dimRange)
00231     template <class Fy>
00232     void assign(const Derivatives<Fy,dimD,dimR,deriv,derivative> &y) 
00233     {
00234       Base::assign(y);
00235       for (int rr=0;rr<dimR;++rr)
00236         tensor_[rr] = y[rr].template tensor<deriv>()[0];
00237     }
00238     // assign with rth component of function
00239     template <class Fy,int dimRy>
00240     void assign(const Derivatives<Fy,dimD,dimRy,deriv,value> &y,unsigned int r) 
00241     {
00242       assign<Fy,dimRy>(y.block(),r);
00243     }
00244     // assign with scalar functions to component r
00245     template <class Fy>
00246     void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,value> &y) 
00247     {
00248       assign(r,y.block());
00249     }
00250     template <class Fy>
00251     void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,derivative> &y) 
00252     {
00253       assign(r,y[0]);
00254     }
00255 
00256     Block &block() 
00257     {
00258       return reinterpret_cast<Block&>(*this);
00259     }
00260     const Block &block() const
00261     {
00262       return reinterpret_cast<const Block&>(*this);
00263     }
00264 
00265     template <unsigned int dorder>
00266     const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor() const
00267     {
00268       // use integral_constant<int,...> here to stay compatible with Int2Type
00269       const integral_constant<int,dorder> a = {};
00270       return tensor(a);
00271     }
00272     template <unsigned int dorder>
00273     Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor() 
00274     {
00275       // use integral_constant<int,...> here to stay compatible with Int2Type
00276       return tensor(integral_constant<int,dorder>());
00277     }
00278     template <unsigned int dorder>
00279     const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block() const
00280     {
00281       // use integral_constant<int,...> here to stay compatible with Int2Type
00282       const integral_constant<int,dorder> a = {};
00283       return reinterpret_cast<const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
00284     }
00285     template <unsigned int dorder>
00286     Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block() 
00287     {
00288       // use integral_constant<int,...> here to stay compatible with Int2Type
00289       const integral_constant<int,dorder> a = {};
00290       return reinterpret_cast<Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
00291     }
00292     ThisLFETensor &operator[](int r) {
00293       return tensor_[r];
00294     }
00295     const ThisLFETensor &operator[](int r) const {
00296       return tensor_[r];
00297     }
00298     protected:
00299     template <class Fy,int dimRy>
00300     void assign(const FieldVector<Fy,size*dimRy> &y,unsigned int r)
00301     {
00302       Base::template assign<Fy,dimRy>(reinterpret_cast<const FieldVector<Fy,Base::size*dimRy>&>(y),r);
00303       tensor_[0] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[Base::size*dimRy+r*ThisLFETensor::size]);
00304     }
00305     template <class Fy>
00306     void assign(unsigned int r,const FieldVector<Fy,size/dimR> &y) 
00307     {
00308       Base::assign(r,reinterpret_cast<const FieldVector<Fy,Base::size/dimR>&>(y));
00309       tensor_[r] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[Base::size/dimR]);
00310     }
00311     // assign with diffrent layout (same dimRange)
00312     template <class Fy,unsigned int dy>
00313     void assign(const Derivatives<Fy,dimD,dimR,dy,derivative> &y) 
00314     {
00315       Base::assign(y);
00316       for (int rr=0;rr<dimR;++rr)
00317         tensor_[rr] = y[rr].template tensor<deriv>()[0];
00318     }
00319 
00320     template <int dorder>
00321     const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
00322     tensor(const integral_constant<int,dorder> &dorderVar) const
00323     {
00324       return Base::tensor(dorderVar);
00325     }
00326     const Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
00327     tensor(const integral_constant<int,deriv> &dorderVar) const
00328     {
00329       return tensor_;
00330     }
00331     template <int dorder>
00332     Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
00333     tensor(const integral_constant<int,dorder> &dorderVar)
00334     {
00335       return Base::tensor(dorderVar);
00336     }
00337     Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
00338     tensor(const integral_constant<int,deriv> &dorderVar)
00339     {
00340       return tensor_;
00341     }
00342     Dune::FieldVector<ThisLFETensor,dimR> tensor_;
00343   };
00344 
00345   template <class F,int dimD,int dimR>
00346   struct Derivatives<F,dimD,dimR,0,value>
00347   {
00348     typedef Derivatives<F,dimD,dimR,0,value> This;
00349     typedef LFETensor<F,dimD,0> ThisLFETensor;
00350 
00351     typedef F Field;
00352     typedef F field_type;
00353 
00354     static const DerivativeLayout layout = value;
00355     static const unsigned int dimDomain = dimD;
00356     static const unsigned int dimRange = dimR;
00357     // size needs to be an anonymous enum value for gcc 3.4 compatibility
00358     enum { size = ThisLFETensor::size*dimR };
00359     typedef Dune::FieldVector<F,size> Block;
00360 
00361     template <class FF>
00362     This &operator=(const FF& f)
00363     {
00364       for (int r=0;r<dimR;++r)
00365         tensor_[r] = field_cast<F>(f);
00366       return *this;
00367     }
00368     This &operator=(const Dune::FieldVector<ThisLFETensor,dimR> &t) 
00369     {
00370       tensor_ = t;
00371       return *this;
00372     }
00373 
00374     This &operator=(const Block &t) 
00375     {
00376       block() = t;
00377       return *this;
00378     }
00379 
00380     This &operator*= ( const field_type &f )
00381     {
00382       block() *= f;
00383       return *this;
00384     }
00385 
00386     void axpy(const F &a, const This &y)
00387     {
00388       block().axpy(a,y.block());
00389     }
00390     template <class Fy>
00391     void assign(const Derivatives<Fy,dimD,dimR,0,value> &y) 
00392     {
00393       field_cast(y.block(),block());
00394     }
00395     template <class Fy>
00396     void assign(const Derivatives<Fy,dimD,dimR,0,derivative> &y) 
00397     {
00398       for (int rr=0;rr<dimR;++rr)
00399         tensor_[rr] = y[rr].template tensor<0>()[0];
00400     }
00401     template <class Fy,int dimRy>
00402     void assign(const Derivatives<Fy,dimD,dimRy,0,value> &y,unsigned int r) 
00403     {
00404       assign<Fy,dimRy>(y.block(),r);
00405     }
00406     template <class Fy>
00407     void assign(unsigned int r,const Derivatives<Fy,dimD,1,0,value> &y) 
00408     {
00409       tensor_[r].assign(y[0]);
00410     }
00411     template <class Fy>
00412     void assign(unsigned int r,const Derivatives<Fy,dimD,1,0,derivative> &y) 
00413     {
00414       tensor_[r].assign(y[0][0]);
00415     }
00416 
00417     Block &block() 
00418     {
00419       return reinterpret_cast<Block&>(*this);
00420     }
00421     const Block &block() const
00422     {
00423       return reinterpret_cast<const Block&>(*this);
00424     }
00425 
00426     ThisLFETensor &operator[](int r) {
00427       return tensor_[r];
00428     }
00429     const ThisLFETensor &operator[](int r) const {
00430       return tensor_[r];
00431     }
00432     template <int dorder>
00433     const Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &tensor() const
00434     {
00435       return tensor_;
00436     }
00437     Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &tensor() 
00438     {
00439       return tensor_;
00440     }
00441     template <unsigned int dorder>
00442     const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block() const
00443     {
00444       // use integral_constant<int,...> here to stay compatible with Int2Type
00445       const integral_constant<int,dorder> a = {};
00446       return reinterpret_cast<const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
00447     }
00448     template <unsigned int dorder>
00449     Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block() 
00450     {
00451       // use integral_constant<int,...> here to stay compatible with Int2Type
00452       const integral_constant<int,dorder> a = {};
00453       return reinterpret_cast<Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
00454     }
00455 
00456     protected:
00457     const Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &
00458     tensor(const integral_constant<int,0> &dorderVar) const
00459     {
00460       return tensor_;
00461     }
00462     Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &
00463     tensor(const integral_constant<int,0> &dorderVar)
00464     {
00465       return tensor_;
00466     }
00467     template <class Fy,unsigned int dy>
00468     void assign(const Derivatives<Fy,dimD,dimR,dy,derivative> &y) 
00469     {
00470       for (int rr=0;rr<dimR;++rr)
00471         tensor_[rr] = y[rr].template tensor<0>()[0];
00472     }
00473     template <class Fy,int dimRy>
00474     void assign(const FieldVector<Fy,size*dimRy> &y,unsigned int r) 
00475     {
00476       tensor_[0] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[r*ThisLFETensor::size]);
00477     }
00478     template <class Fy>
00479     void assign(unsigned int r,const FieldVector<Fy,size/dimR> &y) 
00480     {
00481       tensor_[r] = y;
00482     }
00483     Dune::FieldVector<ThisLFETensor,dimR> tensor_;
00484   };
00485 
00486   // Implemnetation for derivative based layout
00487   template <class F,int dimD,int dimR,unsigned int deriv>
00488   struct Derivatives<F,dimD,dimR,deriv,derivative>
00489   {
00490     typedef Derivatives<F,dimD,dimR,deriv,derivative> This;
00491     typedef Derivatives<F,dimD,1,deriv,value> ScalarDeriv;
00492 
00493     typedef F Field;
00494     typedef F field_type;
00495 
00496     static const DerivativeLayout layout = value;
00497     static const unsigned int dimDomain = dimD;
00498     static const unsigned int dimRange = dimR;
00499     // size needs to be an anonymous enum value for gcc 3.4 compatibility
00500     enum { size = ScalarDeriv::size*dimR };
00501     typedef Dune::FieldVector<F,size> Block;
00502 
00503     template <class FF>
00504     This &operator=(const FF& f)
00505     {
00506       block() = field_cast<F>(f);
00507       return *this;
00508     }
00509     This &operator=(const Block &t) 
00510     {
00511       block() = t;
00512       return *this;
00513     }
00514 
00515     This &operator*= ( const field_type &f )
00516     {
00517       block() *= f;
00518       return *this;
00519     }
00520 
00521     template <class FF>
00522     void axpy(const FF &a, const This &y)
00523     {
00524       block().axpy(field_cast<F>(a),y.block());
00525     }
00526     // assign with same layout (only diffrent Field)
00527     template <class Fy>
00528     void assign(const Derivatives<Fy,dimD,dimR,deriv,derivative> &y) 
00529     {
00530       field_cast(y.block(),block());
00531     }
00532     // assign with diffrent layout (same dimRange)
00533     template <class Fy>
00534     void assign(const Derivatives<Fy,dimD,dimR,deriv,value> &y) 
00535     {
00536       for (unsigned int rr=0;rr<dimR;++rr)
00537         deriv_[rr].assign(y,rr);
00538     }
00539     // assign with scalar functions to component r
00540     template <class Fy,DerivativeLayout layouty>
00541     void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,layouty> &y) 
00542     {
00543       deriv_[r].assign(r,y);
00544     }
00545 
00546     Block &block() 
00547     {
00548       return reinterpret_cast<Block&>(*this);
00549     }
00550     const Block &block() const
00551     {
00552       return reinterpret_cast<const Block&>(*this);
00553     }
00554 
00555     ScalarDeriv &operator[](int r) {
00556       return deriv_[r];
00557     }
00558     const ScalarDeriv &operator[](int r) const {
00559       return deriv_[r];
00560     }
00561     protected:
00562     Dune::FieldVector<ScalarDeriv,dimR> deriv_;
00563   };
00564 
00565   // ******************************************
00566   // AXPY *************************************
00567   // ******************************************
00568   template <class Vec1,class Vec2,unsigned int deriv>
00569   struct LFETensorAxpy
00570   {
00571     template <class Field>
00572     static void apply(unsigned int r,const Field &a,
00573                       const Vec1 &x, Vec2 &y) 
00574     {
00575       y.axpy(a,x);
00576     }
00577   };
00578   template <class F1,int dimD,int dimR,
00579             unsigned int d,
00580             class Vec2,
00581             unsigned int deriv>
00582   struct LFETensorAxpy<Derivatives<F1,dimD,dimR,d,value>,Vec2,deriv>
00583   {
00584     typedef Derivatives<F1,dimD,dimR,d,value> Vec1;
00585     template <class Field>
00586     static void apply(unsigned int r,const Field &a,
00587                       const Vec1 &x, Vec2 &y) 
00588     {
00589       const FieldVector<F1,Vec2::size> &xx = x.template block<deriv>();
00590       for (int i=0;i<y.size;++i)
00591         y[i] += xx[i]*a;
00592     }
00593   };
00594   template <class F1,int dimD,int dimR,
00595             unsigned int d,
00596             class Vec2,
00597             unsigned int deriv>
00598   struct LFETensorAxpy<Derivatives<F1,dimD,dimR,d,derivative>,Vec2,deriv>
00599   {
00600     typedef Derivatives<F1,dimD,dimR,d,derivative> Vec1;
00601     template <class Field>
00602     static void apply(unsigned int r,const Field &a,
00603                       const Vec1 &x, Vec2 &y) 
00604     {
00605       for (int rr=0;rr<dimR;++rr)
00606         LFETensorAxpy<Derivatives<F1,dimD,1,d,value>,
00607                    Vec2,deriv>::apply(rr,a,x[rr],y);
00608     }
00609   };
00610   template <class F1,int dimD,
00611             unsigned int d,
00612             class Vec2,
00613             unsigned int deriv>
00614   struct LFETensorAxpy<Derivatives<F1,dimD,1,d,derivative>,Vec2,deriv>
00615   {
00616     typedef Derivatives<F1,dimD,1,d,derivative> Vec1;
00617     template <class Field>
00618     static void apply(unsigned int r,const Field &a,
00619                       const Vec1 &x, Vec2 &y) 
00620     {
00621       LFETensorAxpy<Derivatives<F1,dimD,1,d,value>,
00622                  Vec2,deriv>::apply(r,a,x[0],y);
00623     }
00624   };
00625   template <class F1,int dimD,
00626             unsigned int d,
00627             class Vec2,
00628             unsigned int deriv>
00629   struct LFETensorAxpy<Derivatives<F1,dimD,1,d,value>,Vec2,deriv>
00630   {
00631     typedef Derivatives<F1,dimD,1,d,value> Vec1;
00632     template <class Field>
00633     static void apply(unsigned int r,const Field &a,
00634                       const Vec1 &x, Vec2 &y) 
00635     {
00636       typedef LFETensor<F1,dimD,deriv> LFETensorType;
00637       const unsigned int rr = r*LFETensorType::size;
00638       const FieldVector<F1,LFETensorType::size> &xx = x.template block<deriv>();
00639       for (int i=0;i<FieldVector<F1,LFETensorType::size>::dimension;++i)
00640         y[rr+i] += xx[i]*a;
00641     }
00642   };
00643 
00644   // ***********************************************
00645   // Assign ****************************************
00646   // ***********************************************
00647   template <class Vec1,class Vec2>
00648   struct DerivativeAssign
00649   {
00650     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00651     {
00652       field_cast(vec1,vec2);
00653     }
00654   };
00655   template <int dimD,int dimR,unsigned int deriv, DerivativeLayout layout,
00656            class F1,class F2>
00657   struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,layout>,
00658                           Derivatives<F2,dimD,dimR,deriv,layout> > 
00659   {
00660     typedef Derivatives<F1,dimD,dimR,deriv,layout> Vec1;
00661     typedef Derivatives<F2,dimD,dimR,deriv,layout> Vec2;
00662     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00663     {
00664       field_cast(vec1.block(),vec2.block());
00665     }
00666   };
00667   template <int dimD,int dimR,unsigned int deriv,
00668            class F1, class F2>
00669   struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,value>,
00670                           Derivatives<F2,dimD,dimR,deriv,derivative> > 
00671   {
00672     typedef Derivatives<F1,dimD,dimR,deriv,value> Vec1;
00673     typedef Derivatives<F2,dimD,dimR,deriv,derivative> Vec2;
00674     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00675     {
00676       vec2.assign(vec1);
00677     }
00678   };
00679   template <int dimD,int dimR,unsigned int deriv,
00680            class F1, class F2>
00681   struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,derivative>,
00682                           Derivatives<F2,dimD,dimR,deriv,value> > 
00683   {
00684     typedef Derivatives<F1,dimD,dimR,deriv,derivative> Vec1;
00685     typedef Derivatives<F2,dimD,dimR,deriv,value> Vec2;
00686     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00687     {
00688       vec2.assign(vec1);
00689     }
00690   };
00691   template <int dimD,int dimR,unsigned int deriv,DerivativeLayout layout,
00692            class F1, class F2>
00693   struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
00694                           Derivatives<F2,dimD,dimR,deriv,value> > 
00695   {
00696     typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
00697     typedef Derivatives<F2,dimD,dimR,deriv,value> Vec2;
00698     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00699     {
00700       vec2.assign(r,vec1);
00701     }
00702   };
00703   template <int dimD,int dimR,unsigned int deriv,DerivativeLayout layout,
00704            class F1, class F2>
00705   struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
00706                           Derivatives<F2,dimD,dimR,deriv,derivative> > 
00707   {
00708     typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
00709     typedef Derivatives<F2,dimD,dimR,deriv,derivative> Vec2;
00710     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00711     {
00712       vec2.assign(r,vec1);
00713     }
00714   };
00715   template <int dimD,unsigned int deriv,
00716            class F1, class F2>
00717   struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,value>,
00718                           Derivatives<F2,dimD,1,deriv,value> > 
00719   {
00720     typedef Derivatives<F1,dimD,1,deriv,value> Vec1;
00721     typedef Derivatives<F2,dimD,1,deriv,value> Vec2;
00722     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00723     {
00724       field_cast(vec1.block(),vec2.block());
00725     }
00726   };
00727   template <int dimD,unsigned int deriv,
00728            class F1, class F2>
00729   struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,derivative>,
00730                           Derivatives<F2,dimD,1,deriv,derivative> > 
00731   {
00732     typedef Derivatives<F1,dimD,1,deriv,derivative> Vec1;
00733     typedef Derivatives<F2,dimD,1,deriv,derivative> Vec2;
00734     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00735     {
00736       field_cast(vec1.block(),vec2.block());
00737     }
00738   };
00739   template <int dimD,unsigned int deriv,
00740            class F1, class F2>
00741   struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,derivative>,
00742                           Derivatives<F2,dimD,1,deriv,value> > 
00743   {
00744     typedef Derivatives<F1,dimD,1,deriv,derivative> Vec1;
00745     typedef Derivatives<F2,dimD,1,deriv,value> Vec2;
00746     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00747     {
00748       field_cast(vec1.block(),vec2.block());
00749     }
00750   };
00751   template <int dimD,unsigned int deriv,
00752            class F1, class F2>
00753   struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,value>,
00754                           Derivatives<F2,dimD,1,deriv,derivative> > 
00755   {
00756     typedef Derivatives<F1,dimD,1,deriv,value> Vec1;
00757     typedef Derivatives<F2,dimD,1,deriv,derivative> Vec2;
00758     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00759     {
00760       field_cast(vec1.block(),vec2.block());
00761     }
00762   };
00763   template <int dimD,unsigned int deriv,DerivativeLayout layout,
00764            class F1, class F2>
00765   struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
00766                           F2 > 
00767   {
00768     typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
00769     typedef F2 Vec2;
00770     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00771     {
00772       field_cast(vec1.block(),vec2);
00773     }
00774   };
00775   template <int dimD,int dimR,
00776             class F1,unsigned int deriv,
00777             class F2>
00778   struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,value>,FieldVector<F2,dimR> >
00779   {
00780     typedef Derivatives<F1,dimD,dimR,deriv,value> Vec1;
00781     typedef FieldVector<F2,dimR> Vec2;
00782     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00783     {
00784       field_cast(vec1.template block<0>(),vec2);
00785     }
00786   };
00787   template <int dimD,int dimR,
00788             class F1,unsigned int deriv,
00789             class F2>
00790   struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,derivative>,FieldVector<F2,dimR> >
00791   {
00792     typedef Derivatives<F1,dimD,dimR,deriv,derivative> Vec1;
00793     typedef FieldVector<F2,dimR> Vec2;
00794     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00795     {
00796       for (int rr=0;rr<dimR;++rr)
00797         field_cast(vec1[rr].template tensor<0>()[0].block(),vec2[rr]);
00798     }
00799   };
00800   template <int dimD,
00801             class F1,unsigned int deriv,
00802             class F2,int dimR>
00803   struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,value>,FieldVector<F2,dimR> >
00804   {
00805     typedef Derivatives<F1,dimD,1,deriv,value> Vec1;
00806     typedef FieldVector<F2,dimR> Vec2;
00807     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00808     {
00809       field_cast(vec1.template tensor<0>()[0].block(),vec2[r]);
00810     }
00811   };
00812   template <int dimD,
00813             class F1,unsigned int deriv,
00814             class F2,int dimR>
00815   struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,derivative>,FieldVector<F2,dimR> >
00816   {
00817     typedef Derivatives<F1,dimD,1,deriv,derivative> Vec1;
00818     typedef FieldVector<F2,dimR> Vec2;
00819     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00820     {
00821       field_cast(vec1[0].template tensor<0>()[0].block(),vec2[r]);
00822     }
00823   };
00824   template <int dimD,
00825             class F1,unsigned int deriv,
00826             class F2>
00827   struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,value>,FieldVector<F2,1> >
00828   {
00829     typedef Derivatives<F1,dimD,1,deriv,value> Vec1;
00830     typedef FieldVector<F2,1> Vec2;
00831     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00832     {
00833       field_cast(vec1.template tensor<0>()[0].block(),vec2);
00834     }
00835   };
00836   template <int dimD,
00837             class F1,unsigned int deriv,
00838             class F2>
00839   struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,derivative>,FieldVector<F2,1> >
00840   {
00841     typedef Derivatives<F1,dimD,1,deriv,derivative> Vec1;
00842     typedef FieldVector<F2,1> Vec2;
00843     static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2) 
00844     {
00845       field_cast(vec1[0].template tensor<0>()[0].block(),vec2);
00846     }
00847   };
00848 
00849   // ***********************************************
00850   // IO ********************************************
00851   // ***********************************************
00852   template <class F,int dimD,unsigned int deriv>
00853   std::ostream &operator<< ( std::ostream &out, const LFETensor< F,dimD,deriv > &tensor )
00854   {
00855     return out << tensor.block();
00856   }
00857 #if 0
00858   template <class F,int dimD,unsigned int deriv>
00859   std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,deriv > &d )
00860   {
00861     out << static_cast<const ScalarDerivatives< F,dimD,deriv-1 > &>(d);
00862     out << " , " << d.tensor() << std::endl;
00863     return out;
00864   }
00865   template <class F,int dimD>
00866   std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,0 > &d )
00867   {
00868     out << d.tensor() << std::endl;
00869     return out;
00870   }
00871 #endif
00872   template <class F,int dimD,int dimR,unsigned int deriv>
00873   std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,derivative > &d )
00874   {
00875     out << " ( ";
00876     out << d[0];
00877     for (int r=1;r<dimR;++r) 
00878     {
00879       out << " , " << d[r];
00880     }
00881     out << " ) " << std::endl;
00882     return out;
00883   }
00884   template <class F,int dimD,int dimR,unsigned int deriv>
00885   std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,value > &d )
00886   {
00887     out << static_cast<const Derivatives< F,dimD,dimR,deriv-1,value > &>(d);
00888     out << " ( ";
00889     out << d[0];
00890     for (int r=1;r<dimR;++r) 
00891     {
00892       out << " , " << d[r];
00893     }
00894     out << " ) " << std::endl;
00895     return out;
00896   }
00897   template <class F,int dimD,int dimR>
00898   std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,derivative > &d )
00899   {
00900     out << " ( ";
00901     out << d[0];
00902     for (int r=1;r<dimR;++r) 
00903     {
00904       out << " , " << d[r];
00905     }
00906     out << " ) " << std::endl;
00907     return out;
00908   }
00909   template <class F,int dimD,int dimR>
00910   std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,value > &d )
00911   {
00912     out << " ( ";
00913     out << d[0];
00914     for (int r=1;r<dimR;++r) 
00915     {
00916       out << " , " << d[r];
00917     }
00918     out << " ) " << std::endl;
00919     return out;
00920   }
00921   template <class F,int dimD,int dimR,unsigned int deriv,DerivativeLayout layout>
00922   std::ostream &operator<< ( std::ostream &out, const std::vector<Derivatives< F,dimD,dimR,deriv,layout > > &y )
00923   {
00924     out << "Number of basis functions: " << y.size() << std::endl;
00925     for (unsigned int i=0;i<y.size();++i) 
00926     {
00927       out << "Base " << i << " : " << std::endl;
00928       out << y[i];
00929       out << std::endl;
00930     }
00931     return out;
00932   }
00933 }
00934 #endif // DUNE_TENSOR_HH
00935