dune-localfunctions  2.9.0
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 // SPDX-FileCopyrightInfo: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
4 // SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5 
6 #ifndef DUNE_TENSOR_HH
7 #define DUNE_TENSOR_HH
8 
9 #include <ostream>
10 #include <vector>
11 
12 #include <dune/common/fvector.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
167  namespace DerivativeLayoutNS {
169  }
170  template <class F,int dimD,int dimR,unsigned int deriv,
172  struct Derivatives;
173 
174  // Implemnetation for valued based layout
175  template <class F,int dimD,int dimR,unsigned int deriv>
176  struct Derivatives<F,dimD,dimR,deriv,DerivativeLayoutNS::value>
177  : public Derivatives<F,dimD,dimR,deriv-1,DerivativeLayoutNS::value>
178  {
180  typedef Derivatives<F,dimD,dimR,deriv-1,DerivativeLayoutNS::value> Base;
182 
183  typedef F Field;
184  typedef F field_type;
185 
187  static const unsigned int dimDomain = dimD;
188  static const unsigned int dimRange = dimR;
189  constexpr static int size = Base::size+ThisLFETensor::size*dimR;
190  typedef Dune::FieldVector<F,size> Block;
191 
192  This &operator=(const F& f)
193  {
194  block() = f;
195  return *this;
196  }
197  This &operator=(const Dune::FieldVector<ThisLFETensor,dimR> &t)
198  {
199  tensor_ = t;
200  return *this;
201  }
202  template <unsigned int dorder>
203  This &operator=(const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &t)
204  {
205  tensor<dorder>() = t;
206  return *this;
207  }
208  This &operator=(const Block &t)
209  {
210  block() = t;
211  return *this;
212  }
213 
214  This &operator*= ( const field_type &f )
215  {
216  block() *= f;
217  return *this;
218  }
219 
220  void axpy(const F &a, const This &y)
221  {
222  block().axpy(a,y.block());
223  }
224 
225  // assign with same layout (only different Field)
226  template <class Fy>
228  {
229  field_cast(y.block(),block());
230  }
231  // assign with different layout (same dimRange)
232  template <class Fy>
234  {
235  Base::assign(y);
236  for (int rr=0; rr<dimR; ++rr)
237  tensor_[rr] = y[rr].template tensor<deriv>()[0];
238  }
239  // assign with rth component of function
240  template <class Fy,int dimRy>
242  {
243  assign<Fy,dimRy>(y.block(),r);
244  }
245  // assign with scalar functions to component r
246  template <class Fy>
248  {
249  assign(r,y.block());
250  }
251  template <class Fy>
253  {
254  assign(r,y[0]);
255  }
256 
258  {
259  return reinterpret_cast<Block&>(*this);
260  }
261  const Block &block() const
262  {
263  return reinterpret_cast<const Block&>(*this);
264  }
265 
266  template <unsigned int dorder>
267  const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor() const
268  {
269  // use integral_constant<int,...> here to stay compatible with Int2Type
270  const std::integral_constant<int,dorder> a = {};
271  return tensor(a);
272  }
273  template <unsigned int dorder>
274  Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor()
275  {
276  // use integral_constant<int,...> here to stay compatible with Int2Type
277  return tensor(std::integral_constant<int,dorder>());
278  }
279  template <unsigned int dorder>
280  const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block() const
281  {
282  // use integral_constant<int,...> here to stay compatible with Int2Type
283  const std::integral_constant<int,dorder> a = {};
284  return reinterpret_cast<const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
285  }
286  template <unsigned int dorder>
287  Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
288  {
289  // use integral_constant<int,...> here to stay compatible with Int2Type
290  const std::integral_constant<int,dorder> a = {};
291  return reinterpret_cast<Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
292  }
294  return tensor_[r];
295  }
296  const ThisLFETensor &operator[](int r) const {
297  return tensor_[r];
298  }
299  protected:
300  template <class Fy,int dimRy>
301  void assign(const FieldVector<Fy,size*dimRy> &y,unsigned int r)
302  {
303  Base::template assign<Fy,dimRy>(reinterpret_cast<const FieldVector<Fy,Base::size*dimRy>&>(y),r);
304  tensor_[0] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[Base::size*dimRy+r*ThisLFETensor::size]);
305  }
306  template <class Fy>
307  void assign(unsigned int r,const FieldVector<Fy,size/dimR> &y)
308  {
309  Base::assign(r,reinterpret_cast<const FieldVector<Fy,Base::size/dimR>&>(y));
310  tensor_[r] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[Base::size/dimR]);
311  }
312  // assign with different layout (same dimRange)
313  template <class Fy,unsigned int dy>
315  {
316  Base::assign(y);
317  for (int rr=0; rr<dimR; ++rr)
318  tensor_[rr] = y[rr].template tensor<deriv>()[0];
319  }
320 
321  template <int dorder>
322  const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
323  tensor(const std::integral_constant<int,dorder> &dorderVar) const
324  {
325  return Base::tensor(dorderVar);
326  }
327  const Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
328  tensor(const std::integral_constant<int,deriv> &dorderVar) const
329  {
330  return tensor_;
331  }
332  template <int dorder>
333  Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
334  tensor(const std::integral_constant<int,dorder> &dorderVar)
335  {
336  return Base::tensor(dorderVar);
337  }
338  Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
339  tensor(const std::integral_constant<int,deriv> &dorderVar)
340  {
341  return tensor_;
342  }
343  Dune::FieldVector<ThisLFETensor,dimR> tensor_;
344  };
345 
346  template <class F,int dimD,int dimR>
347  struct Derivatives<F,dimD,dimR,0,DerivativeLayoutNS::value>
348  {
351 
352  typedef F Field;
353  typedef F field_type;
354 
356  static const unsigned int dimDomain = dimD;
357  static const unsigned int dimRange = dimR;
358  constexpr static int 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>
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>
403  {
404  assign<Fy,dimRy>(y.block(),r);
405  }
406  template <class Fy>
408  {
409  tensor_[r].assign(y[0]);
410  }
411  template <class Fy>
413  {
414  tensor_[r].assign(y[0][0]);
415  }
416 
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 std::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 std::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 std::integral_constant<int,0> &dorderVar) const
459  {
460  return tensor_;
461  }
462  Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &
463  tensor(const std::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 DerivativeLayoutNS::derivative based layout
487  template <class F,int dimD,int dimR,unsigned int deriv>
488  struct Derivatives<F,dimD,dimR,deriv,DerivativeLayoutNS::derivative>
489  {
492 
493  typedef F Field;
494  typedef F field_type;
495 
497  static const unsigned int dimDomain = dimD;
498  static const unsigned int dimRange = dimR;
499  constexpr static int size = ScalarDeriv::size*dimR;
500  typedef Dune::FieldVector<F,size> Block;
501 
502  template <class FF>
503  This &operator=(const FF& f)
504  {
505  block() = field_cast<F>(f);
506  return *this;
507  }
508  This &operator=(const Block &t)
509  {
510  block() = t;
511  return *this;
512  }
513 
514  This &operator*= ( const field_type &f )
515  {
516  block() *= f;
517  return *this;
518  }
519 
520  template <class FF>
521  void axpy(const FF &a, const This &y)
522  {
523  block().axpy(field_cast<F>(a),y.block());
524  }
525  // assign with same layout (only different Field)
526  template <class Fy>
528  {
529  field_cast(y.block(),block());
530  }
531  // assign with different layout (same dimRange)
532  template <class Fy>
534  {
535  for (unsigned int rr=0; rr<dimR; ++rr)
536  deriv_[rr].assign(y,rr);
537  }
538  // assign with scalar functions to component r
539  template <class Fy,DerivativeLayoutNS::DerivativeLayout layouty>
540  void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,layouty> &y)
541  {
542  deriv_[r].assign(r,y);
543  }
544 
546  {
547  return reinterpret_cast<Block&>(*this);
548  }
549  const Block &block() const
550  {
551  return reinterpret_cast<const Block&>(*this);
552  }
553 
555  return deriv_[r];
556  }
557  const ScalarDeriv &operator[](int r) const {
558  return deriv_[r];
559  }
560  protected:
561  Dune::FieldVector<ScalarDeriv,dimR> deriv_;
562  };
563 
564  // ******************************************
565  // AXPY *************************************
566  // ******************************************
567  template <class Vec1,class Vec2,unsigned int deriv>
569  {
570  template <class Field>
571  static void apply(unsigned int r,const Field &a,
572  const Vec1 &x, Vec2 &y)
573  {
574  y.axpy(a,x);
575  }
576  };
577  template <class F1,int dimD,int dimR,
578  unsigned int d,
579  class Vec2,
580  unsigned int deriv>
581  struct LFETensorAxpy<Derivatives<F1,dimD,dimR,d,DerivativeLayoutNS::value>,Vec2,deriv>
582  {
584  template <class Field>
585  static void apply(unsigned int r,const Field &a,
586  const Vec1 &x, Vec2 &y)
587  {
588  const FieldVector<F1,Vec2::size> &xx = x.template block<deriv>();
589  for (int i=0; i<y.size; ++i)
590  y[i] += xx[i]*a;
591  }
592  };
593  template <class F1,int dimD,int dimR,
594  unsigned int d,
595  class Vec2,
596  unsigned int deriv>
597  struct LFETensorAxpy<Derivatives<F1,dimD,dimR,d,DerivativeLayoutNS::derivative>,Vec2,deriv>
598  {
600  template <class Field>
601  static void apply(unsigned int r,const Field &a,
602  const Vec1 &x, Vec2 &y)
603  {
604  for (int rr=0; rr<dimR; ++rr)
606  Vec2,deriv>::apply(rr,a,x[rr],y);
607  }
608  };
609  template <class F1,int dimD,
610  unsigned int d,
611  class Vec2,
612  unsigned int deriv>
613  struct LFETensorAxpy<Derivatives<F1,dimD,1,d,DerivativeLayoutNS::derivative>,Vec2,deriv>
614  {
616  template <class Field>
617  static void apply(unsigned int r,const Field &a,
618  const Vec1 &x, Vec2 &y)
619  {
621  Vec2,deriv>::apply(r,a,x[0],y);
622  }
623  };
624  template <class F1,int dimD,
625  unsigned int d,
626  class Vec2,
627  unsigned int deriv>
628  struct LFETensorAxpy<Derivatives<F1,dimD,1,d,DerivativeLayoutNS::value>,Vec2,deriv>
629  {
631  template <class Field>
632  static void apply(unsigned int r,const Field &a,
633  const Vec1 &x, Vec2 &y)
634  {
635  typedef LFETensor<F1,dimD,deriv> LFETensorType;
636  const unsigned int rr = r*LFETensorType::size;
637  const FieldVector<F1,LFETensorType::size> &xx = x.template block<deriv>();
638  for (int i=0; i<FieldVector<F1,LFETensorType::size>::dimension; ++i)
639  y[rr+i] += xx[i]*a;
640  }
641  };
642 
643  // ***********************************************
644  // Assign ****************************************
645  // ***********************************************
646  template <class Vec1,class Vec2>
648  {
649  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
650  {
651  field_cast(vec1,vec2);
652  }
653  };
654  template <int dimD,int dimR,unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout,
655  class F1,class F2>
656  struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,layout>,
657  Derivatives<F2,dimD,dimR,deriv,layout> >
658  {
661  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
662  {
663  field_cast(vec1.block(),vec2.block());
664  }
665  };
666  template <int dimD,int dimR,unsigned int deriv,
667  class F1, class F2>
668  struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value>,
669  Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> >
670  {
673  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
674  {
675  vec2.assign(vec1);
676  }
677  };
678  template <int dimD,int dimR,unsigned int deriv,
679  class F1, class F2>
680  struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative>,
681  Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> >
682  {
685  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
686  {
687  vec2.assign(vec1);
688  }
689  };
690  template <int dimD,int dimR,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout,
691  class F1, class F2>
692  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
693  Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> >
694  {
697  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
698  {
699  vec2.assign(r,vec1);
700  }
701  };
702  template <int dimD,int dimR,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout,
703  class F1, class F2>
704  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
705  Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> >
706  {
709  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
710  {
711  vec2.assign(r,vec1);
712  }
713  };
714  template <int dimD,unsigned int deriv,
715  class F1, class F2>
716  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,
717  Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> >
718  {
721  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
722  {
723  field_cast(vec1.block(),vec2.block());
724  }
725  };
726  template <int dimD,unsigned int deriv,
727  class F1, class F2>
728  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,
729  Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> >
730  {
733  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
734  {
735  field_cast(vec1.block(),vec2.block());
736  }
737  };
738  template <int dimD,unsigned int deriv,
739  class F1, class F2>
740  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,
741  Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> >
742  {
745  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
746  {
747  field_cast(vec1.block(),vec2.block());
748  }
749  };
750  template <int dimD,unsigned int deriv,
751  class F1, class F2>
752  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,
753  Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> >
754  {
757  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
758  {
759  field_cast(vec1.block(),vec2.block());
760  }
761  };
762  template <int dimD,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout,
763  class F1, class F2>
764  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
765  F2 >
766  {
768  typedef F2 Vec2;
769  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
770  {
771  field_cast(vec1.block(),vec2);
772  }
773  };
774  template <int dimD,int dimR,
775  class F1,unsigned int deriv,
776  class F2>
777  struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value>,FieldVector<F2,dimR> >
778  {
780  typedef FieldVector<F2,dimR> Vec2;
781  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
782  {
783  field_cast(vec1.template block<0>(),vec2);
784  }
785  };
786  template <int dimD,int dimR,
787  class F1,unsigned int deriv,
788  class F2>
789  struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative>,FieldVector<F2,dimR> >
790  {
792  typedef FieldVector<F2,dimR> Vec2;
793  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
794  {
795  for (int rr=0; rr<dimR; ++rr)
796  field_cast(vec1[rr].template tensor<0>()[0].block(),vec2[rr]);
797  }
798  };
799  template <int dimD,
800  class F1,unsigned int deriv,
801  class F2,int dimR>
802  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,FieldVector<F2,dimR> >
803  {
805  typedef FieldVector<F2,dimR> Vec2;
806  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
807  {
808  field_cast(vec1.template tensor<0>()[0].block(),vec2[r]);
809  }
810  };
811  template <int dimD,
812  class F1,unsigned int deriv,
813  class F2,int dimR>
814  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,FieldVector<F2,dimR> >
815  {
817  typedef FieldVector<F2,dimR> Vec2;
818  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
819  {
820  field_cast(vec1[0].template tensor<0>()[0].block(),vec2[r]);
821  }
822  };
823  template <int dimD,
824  class F1,unsigned int deriv,
825  class F2>
826  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,FieldVector<F2,1> >
827  {
829  typedef FieldVector<F2,1> Vec2;
830  static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
831  {
832  field_cast(vec1.template tensor<0>()[0].block(),vec2);
833  }
834  };
835  template <int dimD,
836  class F1,unsigned int deriv,
837  class F2>
838  struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,FieldVector<F2,1> >
839  {
841  typedef FieldVector<F2,1> Vec2;
842  static void apply(unsigned int /*r*/,const Vec1 &vec1,Vec2 &vec2)
843  {
844  field_cast(vec1[0].template tensor<0>()[0].block(),vec2);
845  }
846  };
847 
848  // ***********************************************
849  // IO ********************************************
850  // ***********************************************
851  template <class F,int dimD,unsigned int deriv>
852  std::ostream &operator<< ( std::ostream &out, const LFETensor< F,dimD,deriv > &tensor )
853  {
854  return out << tensor.block();
855  }
856 #if 0
857  template <class F,int dimD,unsigned int deriv>
858  std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,deriv > &d )
859  {
860  out << static_cast<const ScalarDerivatives< F,dimD,deriv-1 > &>(d);
861  out << " , " << d.tensor() << std::endl;
862  return out;
863  }
864  template <class F,int dimD>
865  std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,0 > &d )
866  {
867  out << d.tensor() << std::endl;
868  return out;
869  }
870 #endif
871  template <class F,int dimD,int dimR,unsigned int deriv>
873  {
874  out << " ( ";
875  out << d[0];
876  for (int r=1; r<dimR; ++r)
877  {
878  out << " , " << d[r];
879  }
880  out << " ) " << std::endl;
881  return out;
882  }
883  template <class F,int dimD,int dimR,unsigned int deriv>
884  std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,DerivativeLayoutNS::value > &d )
885  {
886  out << static_cast<const Derivatives< F,dimD,dimR,deriv-1,DerivativeLayoutNS::value > &>(d);
887  out << " ( ";
888  out << d[0];
889  for (int r=1; r<dimR; ++r)
890  {
891  out << " , " << d[r];
892  }
893  out << " ) " << std::endl;
894  return out;
895  }
896  template <class F,int dimD,int dimR>
898  {
899  out << " ( ";
900  out << d[0];
901  for (int r=1; r<dimR; ++r)
902  {
903  out << " , " << d[r];
904  }
905  out << " ) " << std::endl;
906  return out;
907  }
908  template <class F,int dimD,int dimR>
909  std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,DerivativeLayoutNS::value > &d )
910  {
911  out << " ( ";
912  out << d[0];
913  for (int r=1; r<dimR; ++r)
914  {
915  out << " , " << d[r];
916  }
917  out << " ) " << std::endl;
918  return out;
919  }
920  template <class F,int dimD,int dimR,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout>
921  std::ostream &operator<< ( std::ostream &out, const std::vector<Derivatives< F,dimD,dimR,deriv,layout > > &y )
922  {
923  out << "Number of basis functions: " << y.size() << std::endl;
924  for (unsigned int i=0; i<y.size(); ++i)
925  {
926  out << "Base " << i << " : " << std::endl;
927  out << y[i];
928  out << std::endl;
929  }
930  return out;
931  }
932 }
933 #endif // DUNE_TENSOR_HH
Definition: bdfmcube.hh:18
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:159
std::ostream & operator<<(std::ostream &out, const LFEMatrix< Field > &mat)
Definition: lfematrix.hh:152
DerivativeLayout
Definition: tensor.hh:168
@ derivative
Definition: tensor.hh:168
@ value
Definition: tensor.hh:168
Definition: tensor.hh:33
This & operator*=(const field_type &f)
Definition: tensor.hh:56
Dune::FieldVector< F, size > Block
Definition: tensor.hh:41
const field_type & operator[](const unsigned int i) const
Definition: tensor.hh:62
Block block_
Definition: tensor.hh:89
F field_type
Definition: tensor.hh:39
void axpy(const F &a, const This &y)
Definition: tensor.hh:80
const Block & block() const
Definition: tensor.hh:76
void assign(const LFETensor< Fy, dimD, deriv > &y)
Definition: tensor.hh:85
static const unsigned int size
Definition: tensor.hh:40
This & operator=(const FF &f)
Definition: tensor.hh:44
Block & block()
Definition: tensor.hh:72
Definition: tensor.hh:107
F field_type
Definition: tensor.hh:111
Block block_
Definition: tensor.hh:162
void assign(const LFETensor< Fy, dimD, 0 > &y)
Definition: tensor.hh:149
Block & block()
Definition: tensor.hh:154
void axpy(const F &a, const This &y)
Definition: tensor.hh:144
Dune::FieldVector< F, size > Block
Definition: tensor.hh:113
const Block & block() const
Definition: tensor.hh:158
Definition: tensor.hh:172
Derivatives< F, dimD, dimR, deriv, DerivativeLayoutNS::value > This
Definition: tensor.hh:179
Dune::FieldVector< LFETensor< F, dimD, dorder >, dimR > & tensor(const std::integral_constant< int, dorder > &dorderVar)
Definition: tensor.hh:334
Derivatives< F, dimD, dimR, deriv-1, DerivativeLayoutNS::value > Base
Definition: tensor.hh:180
const Dune::FieldVector< LFETensor< F, dimD, dorder >, dimR > & tensor() const
Definition: tensor.hh:267
void assign(const FieldVector< Fy, size *dimRy > &y, unsigned int r)
Definition: tensor.hh:301
void assign(const Derivatives< Fy, dimD, dimR, deriv, DerivativeLayoutNS::value > &y)
Definition: tensor.hh:227
void assign(const Derivatives< Fy, dimD, dimR, dy, DerivativeLayoutNS::derivative > &y)
Definition: tensor.hh:314
This & operator=(const Dune::FieldVector< ThisLFETensor, dimR > &t)
Definition: tensor.hh:197
LFETensor< F, dimD, deriv > ThisLFETensor
Definition: tensor.hh:181
void assign(unsigned int r, const Derivatives< Fy, dimD, 1, deriv, DerivativeLayoutNS::value > &y)
Definition: tensor.hh:247
This & operator=(const F &f)
Definition: tensor.hh:192
void axpy(const F &a, const This &y)
Definition: tensor.hh:220
void assign(const Derivatives< Fy, dimD, dimR, deriv, DerivativeLayoutNS::derivative > &y)
Definition: tensor.hh:233
Dune::FieldVector< LFETensor< F, dimD, dorder >, dimR > & tensor()
Definition: tensor.hh:274
This & operator=(const Block &t)
Definition: tensor.hh:208
void assign(unsigned int r, const Derivatives< Fy, dimD, 1, deriv, DerivativeLayoutNS::derivative > &y)
Definition: tensor.hh:252
const ThisLFETensor & operator[](int r) const
Definition: tensor.hh:296
Dune::FieldVector< F, LFETensor< F, dimD, dorder >::size *dimR > & block()
Definition: tensor.hh:287
ThisLFETensor & operator[](int r)
Definition: tensor.hh:293
const Dune::FieldVector< F, LFETensor< F, dimD, dorder >::size *dimR > & block() const
Definition: tensor.hh:280
Dune::FieldVector< F, size > Block
Definition: tensor.hh:190
Dune::FieldVector< ThisLFETensor, dimR > tensor_
Definition: tensor.hh:343
void assign(const Derivatives< Fy, dimD, dimRy, deriv, DerivativeLayoutNS::value > &y, unsigned int r)
Definition: tensor.hh:241
const Dune::FieldVector< LFETensor< F, dimD, deriv >, dimR > & tensor(const std::integral_constant< int, deriv > &dorderVar) const
Definition: tensor.hh:328
This & operator=(const Dune::FieldVector< LFETensor< F, dimD, dorder >, dimR > &t)
Definition: tensor.hh:203
const Dune::FieldVector< LFETensor< F, dimD, dorder >, dimR > & tensor(const std::integral_constant< int, dorder > &dorderVar) const
Definition: tensor.hh:323
void assign(unsigned int r, const FieldVector< Fy, size/dimR > &y)
Definition: tensor.hh:307
Dune::FieldVector< LFETensor< F, dimD, deriv >, dimR > & tensor(const std::integral_constant< int, deriv > &dorderVar)
Definition: tensor.hh:339
const Block & block() const
Definition: tensor.hh:261
const Dune::FieldVector< LFETensor< F, dimD, 0 >, dimR > & tensor(const std::integral_constant< int, 0 > &dorderVar) const
Definition: tensor.hh:458
Dune::FieldVector< LFETensor< F, dimD, 0 >, dimR > & tensor()
Definition: tensor.hh:437
LFETensor< F, dimD, 0 > ThisLFETensor
Definition: tensor.hh:350
void assign(const Derivatives< Fy, dimD, dimR, 0, DerivativeLayoutNS::value > &y)
Definition: tensor.hh:391
const Block & block() const
Definition: tensor.hh:421
Derivatives< F, dimD, dimR, 0, DerivativeLayoutNS::value > This
Definition: tensor.hh:349
const Dune::FieldVector< LFETensor< F, dimD, 0 >, dimR > & tensor() const
Definition: tensor.hh:433
void assign(const Derivatives< Fy, dimD, dimR, dy, DerivativeLayoutNS::derivative > &y)
Definition: tensor.hh:468
void assign(const Derivatives< Fy, dimD, dimRy, 0, DerivativeLayoutNS::value > &y, unsigned int r)
Definition: tensor.hh:402
This & operator=(const Dune::FieldVector< ThisLFETensor, dimR > &t)
Definition: tensor.hh:368
This & operator=(const Block &t)
Definition: tensor.hh:374
ThisLFETensor & operator[](int r)
Definition: tensor.hh:426
const ThisLFETensor & operator[](int r) const
Definition: tensor.hh:429
Dune::FieldVector< LFETensor< F, dimD, 0 >, dimR > & tensor(const std::integral_constant< int, 0 > &dorderVar)
Definition: tensor.hh:463
Dune::FieldVector< F, size > Block
Definition: tensor.hh:359
void assign(unsigned int r, const Derivatives< Fy, dimD, 1, 0, DerivativeLayoutNS::derivative > &y)
Definition: tensor.hh:412
void assign(const Derivatives< Fy, dimD, dimR, 0, DerivativeLayoutNS::derivative > &y)
Definition: tensor.hh:396
void assign(unsigned int r, const Derivatives< Fy, dimD, 1, 0, DerivativeLayoutNS::value > &y)
Definition: tensor.hh:407
const Dune::FieldVector< F, LFETensor< F, dimD, dorder >::size *dimR > & block() const
Definition: tensor.hh:442
void assign(const FieldVector< Fy, size *dimRy > &y, unsigned int r)
Definition: tensor.hh:474
Dune::FieldVector< ThisLFETensor, dimR > tensor_
Definition: tensor.hh:483
Dune::FieldVector< F, LFETensor< F, dimD, dorder >::size *dimR > & block()
Definition: tensor.hh:449
void assign(unsigned int r, const FieldVector< Fy, size/dimR > &y)
Definition: tensor.hh:479
This & operator=(const FF &f)
Definition: tensor.hh:362
void axpy(const F &a, const This &y)
Definition: tensor.hh:386
const ScalarDeriv & operator[](int r) const
Definition: tensor.hh:557
void assign(unsigned int r, const Derivatives< Fy, dimD, 1, deriv, layouty > &y)
Definition: tensor.hh:540
This & operator=(const Block &t)
Definition: tensor.hh:508
void assign(const Derivatives< Fy, dimD, dimR, deriv, DerivativeLayoutNS::value > &y)
Definition: tensor.hh:533
Derivatives< F, dimD, dimR, deriv, DerivativeLayoutNS::derivative > This
Definition: tensor.hh:490
Derivatives< F, dimD, 1, deriv, DerivativeLayoutNS::value > ScalarDeriv
Definition: tensor.hh:491
void axpy(const FF &a, const This &y)
Definition: tensor.hh:521
Dune::FieldVector< F, size > Block
Definition: tensor.hh:500
Dune::FieldVector< ScalarDeriv, dimR > deriv_
Definition: tensor.hh:561
void assign(const Derivatives< Fy, dimD, dimR, deriv, DerivativeLayoutNS::derivative > &y)
Definition: tensor.hh:527
Definition: tensor.hh:569
static void apply(unsigned int r, const Field &a, const Vec1 &x, Vec2 &y)
Definition: tensor.hh:571
static void apply(unsigned int r, const Field &a, const Vec1 &x, Vec2 &y)
Definition: tensor.hh:585
Derivatives< F1, dimD, dimR, d, DerivativeLayoutNS::value > Vec1
Definition: tensor.hh:583
Derivatives< F1, dimD, dimR, d, DerivativeLayoutNS::derivative > Vec1
Definition: tensor.hh:599
static void apply(unsigned int r, const Field &a, const Vec1 &x, Vec2 &y)
Definition: tensor.hh:601
Derivatives< F1, dimD, 1, d, DerivativeLayoutNS::derivative > Vec1
Definition: tensor.hh:615
static void apply(unsigned int r, const Field &a, const Vec1 &x, Vec2 &y)
Definition: tensor.hh:617
static void apply(unsigned int r, const Field &a, const Vec1 &x, Vec2 &y)
Definition: tensor.hh:632
Derivatives< F1, dimD, 1, d, DerivativeLayoutNS::value > Vec1
Definition: tensor.hh:630
Definition: tensor.hh:648
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition: tensor.hh:649
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition: tensor.hh:661
Derivatives< F2, dimD, dimR, deriv, DerivativeLayoutNS::value > Vec2
Definition: tensor.hh:696
Derivatives< F2, dimD, dimR, deriv, DerivativeLayoutNS::derivative > Vec2
Definition: tensor.hh:708
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition: tensor.hh:769
Derivatives< F1, dimD, 1, deriv, layout > Vec1
Definition: tensor.hh:767
Derivatives< F1, dimD, dimR, deriv, DerivativeLayoutNS::value > Vec1
Definition: tensor.hh:779
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition: tensor.hh:781
Derivatives< F1, dimD, dimR, deriv, DerivativeLayoutNS::derivative > Vec1
Definition: tensor.hh:791
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition: tensor.hh:793
Derivatives< F1, dimD, 1, deriv, DerivativeLayoutNS::value > Vec1
Definition: tensor.hh:804
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition: tensor.hh:806
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition: tensor.hh:818
Derivatives< F1, dimD, 1, deriv, DerivativeLayoutNS::derivative > Vec1
Definition: tensor.hh:816
Derivatives< F1, dimD, 1, deriv, DerivativeLayoutNS::value > Vec1
Definition: tensor.hh:828
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition: tensor.hh:830
Derivatives< F1, dimD, 1, deriv, DerivativeLayoutNS::derivative > Vec1
Definition: tensor.hh:840
static void apply(unsigned int, const Vec1 &vec1, Vec2 &vec2)
Definition: tensor.hh:842