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