dune-localfunctions  2.3.1
tensor.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_TENSOR_HH
5 #define DUNE_TENSOR_HH
6 
7 #include <ostream>
8 #include <vector>
9 
10 #include <dune/common/fvector.hh>
11 
13 
14 namespace Dune
15 {
16  /***********************************************
17  * The classes here are work in progress.
18  * Basically they provide tensor structures for
19  * higher order derivatives of vector valued function.
20  * Two storage structures are provided
21  * (either based on the components of the vector valued
22  * functions or on the order of the derivative).
23  * Conversions are supplied between the two storage
24  * structures and simple operations, which make the
25  * code difficult to use and requires rewritting...
26  ***************************************************/
27 
28  // Structure for scalar tensor of order deriv
29  template <class F,int dimD,unsigned int deriv>
30  class LFETensor
31  {
33  typedef LFETensor<F,dimD-1,deriv> BaseDim;
34  typedef LFETensor<F,dimD,deriv-1> BaseDeriv;
35 
36  public:
37  typedef F field_type;
38  static const unsigned int size = BaseDim::size+BaseDeriv::size;
39  typedef Dune::FieldVector<F,size> Block;
40 
41  template< class FF >
42  This &operator= ( const FF &f )
43  {
44  block() = field_cast< F >( f );
45  return *this;
46  }
47 
48  This &operator= ( const Block &b )
49  {
50  block() = b;
51  return *this;
52  }
53 
55  {
56  block() *= f;
57  return *this;
58  }
59 
60  const field_type &operator[] ( const unsigned int i ) const
61  {
62  return block()[ i ];
63  }
64 
65  field_type &operator[] ( const unsigned int i )
66  {
67  return block()[ i ];
68  }
69 
71  {
72  return block_;
73  }
74  const Block &block() const
75  {
76  return block_;
77  }
78  void axpy(const F& a, const This &y)
79  {
80  block().axpy(a,y.block());
81  }
82  template <class Fy>
84  {
85  field_cast(y.block(),block());
86  }
88  };
89 
90  // ******************************************
91  template <class F,unsigned int deriv>
92  struct LFETensor<F,0,deriv>
93  {
94  static const int size = 0;
95  };
96 
97  template <class F>
98  struct LFETensor<F,0,0>
99  {
100  static const int size = 1;
101  };
102 
103  template <class F,int dimD>
104  class LFETensor<F,dimD,0>
105  {
106  typedef LFETensor<F,dimD,0> This;
107 
108  public:
109  typedef F field_type;
110  static const int size = 1;
111  typedef Dune::FieldVector<F,size> Block;
112 
113  template< class FF >
114  This &operator= ( const FF &f )
115  {
116  block() = field_cast< F >( f );
117  return *this;
118  }
119 
120  This &operator= ( const Block &b )
121  {
122  block() = b;
123  return *this;
124  }
125 
127  {
128  block() *= f;
129  return *this;
130  }
131 
132  const F &operator[] ( const unsigned int i ) const
133  {
134  return block()[ i ];
135  }
136 
137  F &operator[] ( const unsigned int i )
138  {
139  return block()[ i ];
140  }
141 
142  void axpy(const F& a, const This &y)
143  {
144  block().axpy(a,y.block());
145  }
146  template <class Fy>
148  {
149  field_cast(y.block(),block());
150  }
151 
153  {
154  return block_;
155  }
156  const Block &block() const
157  {
158  return block_;
159  }
161  };
162  // ***********************************************************
163  // Structure for all derivatives up to order deriv
164  // for vector valued function
166  template <class F,int dimD,int dimR,unsigned int deriv,
167  DerivativeLayout layout>
168  struct Derivatives;
169 
170  // Implemnetation for valued based layout
171  template <class F,int dimD,int dimR,unsigned int deriv>
172  struct Derivatives<F,dimD,dimR,deriv,value>
173  : public Derivatives<F,dimD,dimR,deriv-1,value>
174  {
176  typedef Derivatives<F,dimD,dimR,deriv-1,value> Base;
178 
179  typedef F Field;
180  typedef F field_type;
181 
182  static const DerivativeLayout layout = value;
183  static const unsigned int dimDomain = dimD;
184  static const unsigned int dimRange = dimR;
185  // size needs to be an anonymous enum value for gcc 3.4 compatibility
186  enum { size = Base::size+ThisLFETensor::size*dimR };
187  typedef Dune::FieldVector<F,size> Block;
188 
189  This &operator=(const F& f)
190  {
191  block() = f;
192  return *this;
193  }
194  This &operator=(const Dune::FieldVector<ThisLFETensor,dimR> &t)
195  {
196  tensor_ = t;
197  return *this;
198  }
199  template <unsigned int dorder>
200  This &operator=(const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &t)
201  {
202  tensor<dorder>() = t;
203  return *this;
204  }
205  This &operator=(const Block &t)
206  {
207  block() = t;
208  return *this;
209  }
210 
211  This &operator*= ( const field_type &f )
212  {
213  block() *= f;
214  return *this;
215  }
216 
217  void axpy(const F &a, const This &y)
218  {
219  block().axpy(a,y.block());
220  }
221 
222  // assign with same layout (only diffrent Field)
223  template <class Fy>
224  void assign(const Derivatives<Fy,dimD,dimR,deriv,value> &y)
225  {
226  field_cast(y.block(),block());
227  }
228  // assign with diffrent layout (same dimRange)
229  template <class Fy>
230  void assign(const Derivatives<Fy,dimD,dimR,deriv,derivative> &y)
231  {
232  Base::assign(y);
233  for (int rr=0; rr<dimR; ++rr)
234  tensor_[rr] = y[rr].template tensor<deriv>()[0];
235  }
236  // assign with rth component of function
237  template <class Fy,int dimRy>
238  void assign(const Derivatives<Fy,dimD,dimRy,deriv,value> &y,unsigned int r)
239  {
240  assign<Fy,dimRy>(y.block(),r);
241  }
242  // assign with scalar functions to component r
243  template <class Fy>
244  void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,value> &y)
245  {
246  assign(r,y.block());
247  }
248  template <class Fy>
249  void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,derivative> &y)
250  {
251  assign(r,y[0]);
252  }
253 
254  Block &block()
255  {
256  return reinterpret_cast<Block&>(*this);
257  }
258  const Block &block() const
259  {
260  return reinterpret_cast<const Block&>(*this);
261  }
262 
263  template <unsigned int dorder>
264  const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor() const
265  {
266  // use integral_constant<int,...> here to stay compatible with Int2Type
267  const integral_constant<int,dorder> a = {};
268  return tensor(a);
269  }
270  template <unsigned int dorder>
271  Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor()
272  {
273  // use integral_constant<int,...> here to stay compatible with Int2Type
274  return tensor(integral_constant<int,dorder>());
275  }
276  template <unsigned int dorder>
277  const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block() const
278  {
279  // use integral_constant<int,...> here to stay compatible with Int2Type
280  const integral_constant<int,dorder> a = {};
281  return reinterpret_cast<const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
282  }
283  template <unsigned int dorder>
284  Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
285  {
286  // use integral_constant<int,...> here to stay compatible with Int2Type
287  const integral_constant<int,dorder> a = {};
288  return reinterpret_cast<Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
289  }
291  return tensor_[r];
292  }
293  const ThisLFETensor &operator[](int r) const {
294  return tensor_[r];
295  }
296  protected:
297  template <class Fy,int dimRy>
298  void assign(const FieldVector<Fy,size*dimRy> &y,unsigned int r)
299  {
300  Base::template assign<Fy,dimRy>(reinterpret_cast<const FieldVector<Fy,Base::size*dimRy>&>(y),r);
301  tensor_[0] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[Base::size*dimRy+r*ThisLFETensor::size]);
302  }
303  template <class Fy>
304  void assign(unsigned int r,const FieldVector<Fy,size/dimR> &y)
305  {
306  Base::assign(r,reinterpret_cast<const FieldVector<Fy,Base::size/dimR>&>(y));
307  tensor_[r] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[Base::size/dimR]);
308  }
309  // assign with diffrent layout (same dimRange)
310  template <class Fy,unsigned int dy>
311  void assign(const Derivatives<Fy,dimD,dimR,dy,derivative> &y)
312  {
313  Base::assign(y);
314  for (int rr=0; rr<dimR; ++rr)
315  tensor_[rr] = y[rr].template tensor<deriv>()[0];
316  }
317 
318  template <int dorder>
319  const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
320  tensor(const integral_constant<int,dorder> &dorderVar) const
321  {
322  return Base::tensor(dorderVar);
323  }
324  const Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
325  tensor(const integral_constant<int,deriv> &dorderVar) const
326  {
327  return tensor_;
328  }
329  template <int dorder>
330  Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
331  tensor(const integral_constant<int,dorder> &dorderVar)
332  {
333  return Base::tensor(dorderVar);
334  }
335  Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
336  tensor(const integral_constant<int,deriv> &dorderVar)
337  {
338  return tensor_;
339  }
340  Dune::FieldVector<ThisLFETensor,dimR> tensor_;
341  };
342 
343  template <class F,int dimD,int dimR>
344  struct Derivatives<F,dimD,dimR,0,value>
345  {
348 
349  typedef F Field;
350  typedef F field_type;
351 
352  static const DerivativeLayout layout = value;
353  static const unsigned int dimDomain = dimD;
354  static const unsigned int dimRange = dimR;
355  // size needs to be an anonymous enum value for gcc 3.4 compatibility
356  enum { size = ThisLFETensor::size*dimR };
357  typedef Dune::FieldVector<F,size> Block;
358 
359  template <class FF>
360  This &operator=(const FF& f)
361  {
362  for (int r=0; r<dimR; ++r)
363  tensor_[r] = field_cast<F>(f);
364  return *this;
365  }
366  This &operator=(const Dune::FieldVector<ThisLFETensor,dimR> &t)
367  {
368  tensor_ = t;
369  return *this;
370  }
371 
372  This &operator=(const Block &t)
373  {
374  block() = t;
375  return *this;
376  }
377 
378  This &operator*= ( const field_type &f )
379  {
380  block() *= f;
381  return *this;
382  }
383 
384  void axpy(const F &a, const This &y)
385  {
386  block().axpy(a,y.block());
387  }
388  template <class Fy>
389  void assign(const Derivatives<Fy,dimD,dimR,0,value> &y)
390  {
391  field_cast(y.block(),block());
392  }
393  template <class Fy>
394  void assign(const Derivatives<Fy,dimD,dimR,0,derivative> &y)
395  {
396  for (int rr=0; rr<dimR; ++rr)
397  tensor_[rr] = y[rr].template tensor<0>()[0];
398  }
399  template <class Fy,int dimRy>
400  void assign(const Derivatives<Fy,dimD,dimRy,0,value> &y,unsigned int r)
401  {
402  assign<Fy,dimRy>(y.block(),r);
403  }
404  template <class Fy>
405  void assign(unsigned int r,const Derivatives<Fy,dimD,1,0,value> &y)
406  {
407  tensor_[r].assign(y[0]);
408  }
409  template <class Fy>
410  void assign(unsigned int r,const Derivatives<Fy,dimD,1,0,derivative> &y)
411  {
412  tensor_[r].assign(y[0][0]);
413  }
414 
415  Block &block()
416  {
417  return reinterpret_cast<Block&>(*this);
418  }
419  const Block &block() const
420  {
421  return reinterpret_cast<const Block&>(*this);
422  }
423 
425  return tensor_[r];
426  }
427  const ThisLFETensor &operator[](int r) const {
428  return tensor_[r];
429  }
430  template <int dorder>
431  const Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &tensor() const
432  {
433  return tensor_;
434  }
435  Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &tensor()
436  {
437  return tensor_;
438  }
439  template <unsigned int dorder>
440  const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block() const
441  {
442  // use integral_constant<int,...> here to stay compatible with Int2Type
443  const integral_constant<int,dorder> a = {};
444  return reinterpret_cast<const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
445  }
446  template <unsigned int dorder>
447  Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
448  {
449  // use integral_constant<int,...> here to stay compatible with Int2Type
450  const integral_constant<int,dorder> a = {};
451  return reinterpret_cast<Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
452  }
453 
454  protected:
455  const Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &
456  tensor(const integral_constant<int,0> &dorderVar) const
457  {
458  return tensor_;
459  }
460  Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &
461  tensor(const integral_constant<int,0> &dorderVar)
462  {
463  return tensor_;
464  }
465  template <class Fy,unsigned int dy>
466  void assign(const Derivatives<Fy,dimD,dimR,dy,derivative> &y)
467  {
468  for (int rr=0; rr<dimR; ++rr)
469  tensor_[rr] = y[rr].template tensor<0>()[0];
470  }
471  template <class Fy,int dimRy>
472  void assign(const FieldVector<Fy,size*dimRy> &y,unsigned int r)
473  {
474  tensor_[0] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[r*ThisLFETensor::size]);
475  }
476  template <class Fy>
477  void assign(unsigned int r,const FieldVector<Fy,size/dimR> &y)
478  {
479  tensor_[r] = y;
480  }
481  Dune::FieldVector<ThisLFETensor,dimR> tensor_;
482  };
483 
484  // Implemnetation for derivative based layout
485  template <class F,int dimD,int dimR,unsigned int deriv>
486  struct Derivatives<F,dimD,dimR,deriv,derivative>
487  {
489  typedef Derivatives<F,dimD,1,deriv,value> ScalarDeriv;
490 
491  typedef F Field;
492  typedef F field_type;
493 
494  static const DerivativeLayout layout = value;
495  static const unsigned int dimDomain = dimD;
496  static const unsigned int dimRange = dimR;
497  // size needs to be an anonymous enum value for gcc 3.4 compatibility
498  enum { size = ScalarDeriv::size*dimR };
499  typedef Dune::FieldVector<F,size> Block;
500 
501  template <class FF>
502  This &operator=(const FF& f)
503  {
504  block() = field_cast<F>(f);
505  return *this;
506  }
507  This &operator=(const Block &t)
508  {
509  block() = t;
510  return *this;
511  }
512 
513  This &operator*= ( const field_type &f )
514  {
515  block() *= f;
516  return *this;
517  }
518 
519  template <class FF>
520  void axpy(const FF &a, const This &y)
521  {
522  block().axpy(field_cast<F>(a),y.block());
523  }
524  // assign with same layout (only diffrent Field)
525  template <class Fy>
526  void assign(const Derivatives<Fy,dimD,dimR,deriv,derivative> &y)
527  {
528  field_cast(y.block(),block());
529  }
530  // assign with diffrent layout (same dimRange)
531  template <class Fy>
532  void assign(const Derivatives<Fy,dimD,dimR,deriv,value> &y)
533  {
534  for (unsigned int rr=0; rr<dimR; ++rr)
535  deriv_[rr].assign(y,rr);
536  }
537  // assign with scalar functions to component r
538  template <class Fy,DerivativeLayout layouty>
539  void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,layouty> &y)
540  {
541  deriv_[r].assign(r,y);
542  }
543 
544  Block &block()
545  {
546  return reinterpret_cast<Block&>(*this);
547  }
548  const Block &block() const
549  {
550  return reinterpret_cast<const Block&>(*this);
551  }
552 
554  return deriv_[r];
555  }
556  const ScalarDeriv &operator[](int r) const {
557  return deriv_[r];
558  }
559  protected:
560  Dune::FieldVector<ScalarDeriv,dimR> deriv_;
561  };
562 
563  // ******************************************
564  // AXPY *************************************
565  // ******************************************
566  template <class Vec1,class Vec2,unsigned int deriv>
568  {
569  template <class Field>
570  static void apply(unsigned int r,const Field &a,
571  const Vec1 &x, Vec2 &y)
572  {
573  y.axpy(a,x);
574  }
575  };
576  template <class F1,int dimD,int dimR,
577  unsigned int d,
578  class Vec2,
579  unsigned int deriv>
580  struct LFETensorAxpy<Derivatives<F1,dimD,dimR,d,value>,Vec2,deriv>
581  {
582  typedef Derivatives<F1,dimD,dimR,d,value> Vec1;
583  template <class Field>
584  static void apply(unsigned int r,const Field &a,
585  const Vec1 &x, Vec2 &y)
586  {
587  const FieldVector<F1,Vec2::size> &xx = x.template block<deriv>();
588  for (int i=0; i<y.size; ++i)
589  y[i] += xx[i]*a;
590  }
591  };
592  template <class F1,int dimD,int dimR,
593  unsigned int d,
594  class Vec2,
595  unsigned int deriv>
596  struct LFETensorAxpy<Derivatives<F1,dimD,dimR,d,derivative>,Vec2,deriv>
597  {
598  typedef Derivatives<F1,dimD,dimR,d,derivative> Vec1;
599  template <class Field>
600  static void apply(unsigned int r,const Field &a,
601  const Vec1 &x, Vec2 &y)
602  {
603  for (int rr=0; rr<dimR; ++rr)
604  LFETensorAxpy<Derivatives<F1,dimD,1,d,value>,
605  Vec2,deriv>::apply(rr,a,x[rr],y);
606  }
607  };
608  template <class F1,int dimD,
609  unsigned int d,
610  class Vec2,
611  unsigned int deriv>
612  struct LFETensorAxpy<Derivatives<F1,dimD,1,d,derivative>,Vec2,deriv>
613  {
614  typedef Derivatives<F1,dimD,1,d,derivative> Vec1;
615  template <class Field>
616  static void apply(unsigned int r,const Field &a,
617  const Vec1 &x, Vec2 &y)
618  {
620  Vec2,deriv>::apply(r,a,x[0],y);
621  }
622  };
623  template <class F1,int dimD,
624  unsigned int d,
625  class Vec2,
626  unsigned int deriv>
627  struct LFETensorAxpy<Derivatives<F1,dimD,1,d,value>,Vec2,deriv>
628  {
629  typedef Derivatives<F1,dimD,1,d,value> Vec1;
630  template <class Field>
631  static void apply(unsigned int r,const Field &a,
632  const Vec1 &x, Vec2 &y)
633  {
634  typedef LFETensor<F1,dimD,deriv> LFETensorType;
635  const unsigned int rr = r*LFETensorType::size;
636  const FieldVector<F1,LFETensorType::size> &xx = x.template block<deriv>();
637  for (int i=0; i<FieldVector<F1,LFETensorType::size>::dimension; ++i)
638  y[rr+i] += xx[i]*a;
639  }
640  };
641 
642  // ***********************************************
643  // Assign ****************************************
644  // ***********************************************
645  template <class Vec1,class Vec2>
647  {
648  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
649  {
650  field_cast(vec1,vec2);
651  }
652  };
653  template <int dimD,int dimR,unsigned int deriv, DerivativeLayout layout,
654  class F1,class F2>
655  struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,layout>,
656  Derivatives<F2,dimD,dimR,deriv,layout> >
657  {
658  typedef Derivatives<F1,dimD,dimR,deriv,layout> Vec1;
659  typedef Derivatives<F2,dimD,dimR,deriv,layout> Vec2;
660  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
661  {
662  field_cast(vec1.block(),vec2.block());
663  }
664  };
665  template <int dimD,int dimR,unsigned int deriv,
666  class F1, class F2>
667  struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,value>,
668  Derivatives<F2,dimD,dimR,deriv,derivative> >
669  {
670  typedef Derivatives<F1,dimD,dimR,deriv,value> Vec1;
671  typedef Derivatives<F2,dimD,dimR,deriv,derivative> Vec2;
672  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
673  {
674  vec2.assign(vec1);
675  }
676  };
677  template <int dimD,int dimR,unsigned int deriv,
678  class F1, class F2>
679  struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,derivative>,
680  Derivatives<F2,dimD,dimR,deriv,value> >
681  {
682  typedef Derivatives<F1,dimD,dimR,deriv,derivative> Vec1;
683  typedef Derivatives<F2,dimD,dimR,deriv,value> Vec2;
684  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
685  {
686  vec2.assign(vec1);
687  }
688  };
689  template <int dimD,int dimR,unsigned int deriv,DerivativeLayout layout,
690  class F1, class F2>
691  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
692  Derivatives<F2,dimD,dimR,deriv,value> >
693  {
694  typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
695  typedef Derivatives<F2,dimD,dimR,deriv,value> Vec2;
696  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
697  {
698  vec2.assign(r,vec1);
699  }
700  };
701  template <int dimD,int dimR,unsigned int deriv,DerivativeLayout layout,
702  class F1, class F2>
703  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
704  Derivatives<F2,dimD,dimR,deriv,derivative> >
705  {
706  typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
707  typedef Derivatives<F2,dimD,dimR,deriv,derivative> Vec2;
708  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
709  {
710  vec2.assign(r,vec1);
711  }
712  };
713  template <int dimD,unsigned int deriv,
714  class F1, class F2>
715  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,value>,
716  Derivatives<F2,dimD,1,deriv,value> >
717  {
718  typedef Derivatives<F1,dimD,1,deriv,value> Vec1;
719  typedef Derivatives<F2,dimD,1,deriv,value> Vec2;
720  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
721  {
722  field_cast(vec1.block(),vec2.block());
723  }
724  };
725  template <int dimD,unsigned int deriv,
726  class F1, class F2>
727  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,derivative>,
728  Derivatives<F2,dimD,1,deriv,derivative> >
729  {
730  typedef Derivatives<F1,dimD,1,deriv,derivative> Vec1;
731  typedef Derivatives<F2,dimD,1,deriv,derivative> Vec2;
732  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
733  {
734  field_cast(vec1.block(),vec2.block());
735  }
736  };
737  template <int dimD,unsigned int deriv,
738  class F1, class F2>
739  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,derivative>,
740  Derivatives<F2,dimD,1,deriv,value> >
741  {
742  typedef Derivatives<F1,dimD,1,deriv,derivative> Vec1;
743  typedef Derivatives<F2,dimD,1,deriv,value> Vec2;
744  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
745  {
746  field_cast(vec1.block(),vec2.block());
747  }
748  };
749  template <int dimD,unsigned int deriv,
750  class F1, class F2>
751  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,value>,
752  Derivatives<F2,dimD,1,deriv,derivative> >
753  {
754  typedef Derivatives<F1,dimD,1,deriv,value> Vec1;
755  typedef Derivatives<F2,dimD,1,deriv,derivative> Vec2;
756  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
757  {
758  field_cast(vec1.block(),vec2.block());
759  }
760  };
761  template <int dimD,unsigned int deriv,DerivativeLayout layout,
762  class F1, class F2>
763  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
764  F2 >
765  {
766  typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
767  typedef F2 Vec2;
768  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
769  {
770  field_cast(vec1.block(),vec2);
771  }
772  };
773  template <int dimD,int dimR,
774  class F1,unsigned int deriv,
775  class F2>
776  struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,value>,FieldVector<F2,dimR> >
777  {
778  typedef Derivatives<F1,dimD,dimR,deriv,value> Vec1;
779  typedef FieldVector<F2,dimR> Vec2;
780  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
781  {
782  field_cast(vec1.template block<0>(),vec2);
783  }
784  };
785  template <int dimD,int dimR,
786  class F1,unsigned int deriv,
787  class F2>
788  struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,derivative>,FieldVector<F2,dimR> >
789  {
790  typedef Derivatives<F1,dimD,dimR,deriv,derivative> Vec1;
791  typedef FieldVector<F2,dimR> Vec2;
792  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
793  {
794  for (int rr=0; rr<dimR; ++rr)
795  field_cast(vec1[rr].template tensor<0>()[0].block(),vec2[rr]);
796  }
797  };
798  template <int dimD,
799  class F1,unsigned int deriv,
800  class F2,int dimR>
801  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,value>,FieldVector<F2,dimR> >
802  {
803  typedef Derivatives<F1,dimD,1,deriv,value> Vec1;
804  typedef FieldVector<F2,dimR> Vec2;
805  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
806  {
807  field_cast(vec1.template tensor<0>()[0].block(),vec2[r]);
808  }
809  };
810  template <int dimD,
811  class F1,unsigned int deriv,
812  class F2,int dimR>
813  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,derivative>,FieldVector<F2,dimR> >
814  {
815  typedef Derivatives<F1,dimD,1,deriv,derivative> Vec1;
816  typedef FieldVector<F2,dimR> Vec2;
817  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
818  {
819  field_cast(vec1[0].template tensor<0>()[0].block(),vec2[r]);
820  }
821  };
822  template <int dimD,
823  class F1,unsigned int deriv,
824  class F2>
825  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,value>,FieldVector<F2,1> >
826  {
827  typedef Derivatives<F1,dimD,1,deriv,value> Vec1;
828  typedef FieldVector<F2,1> Vec2;
829  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
830  {
831  field_cast(vec1.template tensor<0>()[0].block(),vec2);
832  }
833  };
834  template <int dimD,
835  class F1,unsigned int deriv,
836  class F2>
837  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,derivative>,FieldVector<F2,1> >
838  {
839  typedef Derivatives<F1,dimD,1,deriv,derivative> Vec1;
840  typedef FieldVector<F2,1> Vec2;
841  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
842  {
843  field_cast(vec1[0].template tensor<0>()[0].block(),vec2);
844  }
845  };
846 
847  // ***********************************************
848  // IO ********************************************
849  // ***********************************************
850  template <class F,int dimD,unsigned int deriv>
851  std::ostream &operator<< ( std::ostream &out, const LFETensor< F,dimD,deriv > &tensor )
852  {
853  return out << tensor.block();
854  }
855 #if 0
856  template <class F,int dimD,unsigned int deriv>
857  std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,deriv > &d )
858  {
859  out << static_cast<const ScalarDerivatives< F,dimD,deriv-1 > &>(d);
860  out << " , " << d.tensor() << std::endl;
861  return out;
862  }
863  template <class F,int dimD>
864  std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,0 > &d )
865  {
866  out << d.tensor() << std::endl;
867  return out;
868  }
869 #endif
870  template <class F,int dimD,int dimR,unsigned int deriv>
871  std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,derivative > &d )
872  {
873  out << " ( ";
874  out << d[0];
875  for (int r=1; r<dimR; ++r)
876  {
877  out << " , " << d[r];
878  }
879  out << " ) " << std::endl;
880  return out;
881  }
882  template <class F,int dimD,int dimR,unsigned int deriv>
883  std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,value > &d )
884  {
885  out << static_cast<const Derivatives< F,dimD,dimR,deriv-1,value > &>(d);
886  out << " ( ";
887  out << d[0];
888  for (int r=1; r<dimR; ++r)
889  {
890  out << " , " << d[r];
891  }
892  out << " ) " << std::endl;
893  return out;
894  }
895  template <class F,int dimD,int dimR>
896  std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,derivative > &d )
897  {
898  out << " ( ";
899  out << d[0];
900  for (int r=1; r<dimR; ++r)
901  {
902  out << " , " << d[r];
903  }
904  out << " ) " << std::endl;
905  return out;
906  }
907  template <class F,int dimD,int dimR>
908  std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,value > &d )
909  {
910  out << " ( ";
911  out << d[0];
912  for (int r=1; r<dimR; ++r)
913  {
914  out << " , " << d[r];
915  }
916  out << " ) " << std::endl;
917  return out;
918  }
919  template <class F,int dimD,int dimR,unsigned int deriv,DerivativeLayout layout>
920  std::ostream &operator<< ( std::ostream &out, const std::vector<Derivatives< F,dimD,dimR,deriv,layout > > &y )
921  {
922  out << "Number of basis functions: " << y.size() << std::endl;
923  for (unsigned int i=0; i<y.size(); ++i)
924  {
925  out << "Base " << i << " : " << std::endl;
926  out << y[i];
927  out << std::endl;
928  }
929  return out;
930  }
931 }
932 #endif // DUNE_TENSOR_HH