dune-localfunctions 2.10
Loading...
Searching...
No Matches
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 © 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
16namespace 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 rewriting...
28 ***************************************************/
29
30 // Structure for scalar tensor of order deriv
31 template <class F,int dimD,unsigned int deriv>
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,int dimD,unsigned int deriv>
94 struct FieldTraits<LFETensor<F,dimD,deriv>>
95 {
96 using field_type = F;
97 using real_type = typename FieldTraits<field_type>::real_type;
98 };
99
100 // ******************************************
101 template <class F,unsigned int deriv>
102 struct LFETensor<F,0,deriv>
103 {
104 static const int size = 0;
105 };
106
107 template <class F>
108 struct LFETensor<F,0,0>
109 {
110 static const int size = 1;
111 };
112
113 template <class F,int dimD>
114 class LFETensor<F,dimD,0>
115 {
117
118 public:
119 typedef F field_type;
120 static const int size = 1;
121 typedef Dune::FieldVector<F,size> Block;
122
123 template< class FF >
124 This &operator= ( const FF &f )
125 {
126 block() = field_cast< F >( f );
127 return *this;
128 }
129
130 This &operator= ( const Block &b )
131 {
132 block() = b;
133 return *this;
134 }
135
137 {
138 block() *= f;
139 return *this;
140 }
141
142 const F &operator[] ( const unsigned int i ) const
143 {
144 return block()[ i ];
145 }
146
147 F &operator[] ( const unsigned int i )
148 {
149 return block()[ i ];
150 }
151
152 void axpy(const F& a, const This &y)
153 {
154 block().axpy(a,y.block());
155 }
156 template <class Fy>
158 {
159 field_cast(y.block(),block());
160 }
161
163 {
164 return block_;
165 }
166 const Block &block() const
167 {
168 return block_;
169 }
171 };
172 // ***********************************************************
173 // Structure for all derivatives up to order deriv
174 // for vector valued function
175 namespace DerivativeLayoutNS {
177 }
178 template <class F,int dimD,int dimR,unsigned int deriv,
181
182 template <class F,int dimD,int dimR,unsigned int deriv,
184 struct FieldTraits<Derivatives<F,dimD,dimR,deriv,layout>>
185 {
186 using field_type = F;
187 using real_type = typename FieldTraits<field_type>::real_type;
188 };
189
190 // Implemnetation for valued based layout
191 template <class F,int dimD,int dimR,unsigned int deriv>
192 struct Derivatives<F,dimD,dimR,deriv,DerivativeLayoutNS::value>
193 : public Derivatives<F,dimD,dimR,deriv-1,DerivativeLayoutNS::value>
194 {
196 typedef Derivatives<F,dimD,dimR,deriv-1,DerivativeLayoutNS::value> Base;
198
199 typedef F Field;
200 typedef F field_type;
201
203 static const unsigned int dimDomain = dimD;
204 static const unsigned int dimRange = dimR;
205 constexpr static int size = Base::size+ThisLFETensor::size*dimR;
206 typedef Dune::FieldVector<F,size> Block;
207
208 This &operator=(const F& f)
209 {
210 block() = f;
211 return *this;
212 }
213 This &operator=(const Dune::FieldVector<ThisLFETensor,dimR> &t)
214 {
215 tensor_ = t;
216 return *this;
217 }
218 template <unsigned int dorder>
219 This &operator=(const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &t)
220 {
221 tensor<dorder>() = t;
222 return *this;
223 }
225 {
226 block() = t;
227 return *this;
228 }
229
230 This &operator*= ( const field_type &f )
231 {
232 block() *= f;
233 return *this;
234 }
235
236 void axpy(const F &a, const This &y)
237 {
238 block().axpy(a,y.block());
239 }
240
241 // assign with same layout (only different Field)
242 template <class Fy>
244 {
245 field_cast(y.block(),block());
246 }
247 // assign with different layout (same dimRange)
248 template <class Fy>
250 {
251 Base::assign(y);
252 for (int rr=0; rr<dimR; ++rr)
253 tensor_[rr] = y[rr].template tensor<deriv>()[0];
254 }
255 // assign with rth component of function
256 template <class Fy,int dimRy>
258 {
259 assign<Fy,dimRy>(y.block(),r);
260 }
261 // assign with scalar functions to component r
262 template <class Fy>
264 {
265 assign(r,y.block());
266 }
267 template <class Fy>
269 {
270 assign(r,y[0]);
271 }
272
274 {
275 return reinterpret_cast<Block&>(*this);
276 }
277 const Block &block() const
278 {
279 return reinterpret_cast<const Block&>(*this);
280 }
281
282 template <unsigned int dorder>
283 const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor() const
284 {
285 // use integral_constant<int,...> here to stay compatible with Int2Type
286 const std::integral_constant<int,dorder> a = {};
287 return tensor(a);
288 }
289 template <unsigned int dorder>
290 Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor()
291 {
292 // use integral_constant<int,...> here to stay compatible with Int2Type
293 return tensor(std::integral_constant<int,dorder>());
294 }
295 template <unsigned int dorder>
296 const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block() const
297 {
298 // use integral_constant<int,...> here to stay compatible with Int2Type
299 const std::integral_constant<int,dorder> a = {};
300 return reinterpret_cast<const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
301 }
302 template <unsigned int dorder>
303 Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
304 {
305 // use integral_constant<int,...> here to stay compatible with Int2Type
306 const std::integral_constant<int,dorder> a = {};
307 return reinterpret_cast<Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
308 }
310 return tensor_[r];
311 }
312 const ThisLFETensor &operator[](int r) const {
313 return tensor_[r];
314 }
315 protected:
316 template <class Fy,int dimRy>
317 void assign(const FieldVector<Fy,size*dimRy> &y,unsigned int r)
318 {
319 Base::template assign<Fy,dimRy>(reinterpret_cast<const FieldVector<Fy,Base::size*dimRy>&>(y),r);
320 tensor_[0] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[Base::size*dimRy+r*ThisLFETensor::size]);
321 }
322 template <class Fy>
323 void assign(unsigned int r,const FieldVector<Fy,size/dimR> &y)
324 {
325 Base::assign(r,reinterpret_cast<const FieldVector<Fy,Base::size/dimR>&>(y));
326 tensor_[r] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[Base::size/dimR]);
327 }
328 // assign with different layout (same dimRange)
329 template <class Fy,unsigned int dy>
331 {
332 Base::assign(y);
333 for (int rr=0; rr<dimR; ++rr)
334 tensor_[rr] = y[rr].template tensor<deriv>()[0];
335 }
336
337 template <int dorder>
338 const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
339 tensor(const std::integral_constant<int,dorder> &dorderVar) const
340 {
341 return Base::tensor(dorderVar);
342 }
343 const Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
344 tensor(const std::integral_constant<int,deriv> &dorderVar) const
345 {
346 return tensor_;
347 }
348 template <int dorder>
349 Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
350 tensor(const std::integral_constant<int,dorder> &dorderVar)
351 {
352 return Base::tensor(dorderVar);
353 }
354 Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
355 tensor(const std::integral_constant<int,deriv> &dorderVar)
356 {
357 return tensor_;
358 }
359 Dune::FieldVector<ThisLFETensor,dimR> tensor_;
360 };
361
362 template <class F,int dimD,int dimR>
363 struct Derivatives<F,dimD,dimR,0,DerivativeLayoutNS::value>
364 {
367
368 typedef F Field;
369 typedef F field_type;
370
372 static const unsigned int dimDomain = dimD;
373 static const unsigned int dimRange = dimR;
374 constexpr static int size = ThisLFETensor::size*dimR;
375 typedef Dune::FieldVector<F,size> Block;
376
377 template <class FF>
378 This &operator=(const FF& f)
379 {
380 for (int r=0; r<dimR; ++r)
381 tensor_[r] = field_cast<F>(f);
382 return *this;
383 }
384 This &operator=(const Dune::FieldVector<ThisLFETensor,dimR> &t)
385 {
386 tensor_ = t;
387 return *this;
388 }
389
391 {
392 block() = t;
393 return *this;
394 }
395
396 This &operator*= ( const field_type &f )
397 {
398 block() *= f;
399 return *this;
400 }
401
402 void axpy(const F &a, const This &y)
403 {
404 block().axpy(a,y.block());
405 }
406 template <class Fy>
408 {
409 field_cast(y.block(),block());
410 }
411 template <class Fy>
413 {
414 for (int rr=0; rr<dimR; ++rr)
415 tensor_[rr] = y[rr].template tensor<0>()[0];
416 }
417 template <class Fy,int dimRy>
419 {
420 assign<Fy,dimRy>(y.block(),r);
421 }
422 template <class Fy>
424 {
425 tensor_[r].assign(y[0]);
426 }
427 template <class Fy>
429 {
430 tensor_[r].assign(y[0][0]);
431 }
432
434 {
435 return reinterpret_cast<Block&>(*this);
436 }
437 const Block &block() const
438 {
439 return reinterpret_cast<const Block&>(*this);
440 }
441
443 return tensor_[r];
444 }
445 const ThisLFETensor &operator[](int r) const {
446 return tensor_[r];
447 }
448 template <int dorder>
449 const Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &tensor() const
450 {
451 return tensor_;
452 }
453 Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &tensor()
454 {
455 return tensor_;
456 }
457 template <unsigned int dorder>
458 const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block() const
459 {
460 // use integral_constant<int,...> here to stay compatible with Int2Type
461 const std::integral_constant<int,dorder> a = {};
462 return reinterpret_cast<const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
463 }
464 template <unsigned int dorder>
465 Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
466 {
467 // use integral_constant<int,...> here to stay compatible with Int2Type
468 const std::integral_constant<int,dorder> a = {};
469 return reinterpret_cast<Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
470 }
471
472 protected:
473 const Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &
474 tensor(const std::integral_constant<int,0> &dorderVar) const
475 {
476 return tensor_;
477 }
478 Dune::FieldVector<LFETensor<F,dimD,0>,dimR> &
479 tensor(const std::integral_constant<int,0> &dorderVar)
480 {
481 return tensor_;
482 }
483 template <class Fy,unsigned int dy>
485 {
486 for (int rr=0; rr<dimR; ++rr)
487 tensor_[rr] = y[rr].template tensor<0>()[0];
488 }
489 template <class Fy,int dimRy>
490 void assign(const FieldVector<Fy,size*dimRy> &y,unsigned int r)
491 {
492 tensor_[0] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[r*ThisLFETensor::size]);
493 }
494 template <class Fy>
495 void assign(unsigned int r,const FieldVector<Fy,size/dimR> &y)
496 {
497 tensor_[r] = y;
498 }
499 Dune::FieldVector<ThisLFETensor,dimR> tensor_;
500 };
501
502 // Implemnetation for DerivativeLayoutNS::derivative based layout
503 template <class F,int dimD,int dimR,unsigned int deriv>
504 struct Derivatives<F,dimD,dimR,deriv,DerivativeLayoutNS::derivative>
505 {
508
509 typedef F Field;
510 typedef F field_type;
511
513 static const unsigned int dimDomain = dimD;
514 static const unsigned int dimRange = dimR;
515 constexpr static int size = ScalarDeriv::size*dimR;
516 typedef Dune::FieldVector<F,size> Block;
517
518 template <class FF>
519 This &operator=(const FF& f)
520 {
521 block() = field_cast<F>(f);
522 return *this;
523 }
525 {
526 block() = t;
527 return *this;
528 }
529
530 This &operator*= ( const field_type &f )
531 {
532 block() *= f;
533 return *this;
534 }
535
536 template <class FF>
537 void axpy(const FF &a, const This &y)
538 {
539 block().axpy(field_cast<F>(a),y.block());
540 }
541 // assign with same layout (only different Field)
542 template <class Fy>
544 {
545 field_cast(y.block(),block());
546 }
547 // assign with different layout (same dimRange)
548 template <class Fy>
550 {
551 for (unsigned int rr=0; rr<dimR; ++rr)
552 deriv_[rr].assign(y,rr);
553 }
554 // assign with scalar functions to component r
555 template <class Fy,DerivativeLayoutNS::DerivativeLayout layouty>
556 void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,layouty> &y)
557 {
558 deriv_[r].assign(r,y);
559 }
560
562 {
563 return reinterpret_cast<Block&>(*this);
564 }
565 const Block &block() const
566 {
567 return reinterpret_cast<const Block&>(*this);
568 }
569
571 return deriv_[r];
572 }
573 const ScalarDeriv &operator[](int r) const {
574 return deriv_[r];
575 }
576 protected:
577 Dune::FieldVector<ScalarDeriv,dimR> deriv_;
578 };
579
580 // ******************************************
581 // AXPY *************************************
582 // ******************************************
583 template <class Vec1,class Vec2,unsigned int deriv>
585 {
586 template <class Field>
587 static void apply(unsigned int r,const Field &a,
588 const Vec1 &x, Vec2 &y)
589 {
590 y.axpy(a,x);
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::value>,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 const FieldVector<F1,Vec2::size> &xx = x.template block<deriv>();
605 for (int i=0; i<y.size; ++i)
606 y[i] += xx[i]*a;
607 }
608 };
609 template <class F1,int dimD,int dimR,
610 unsigned int d,
611 class Vec2,
612 unsigned int deriv>
613 struct LFETensorAxpy<Derivatives<F1,dimD,dimR,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 {
620 for (int rr=0; rr<dimR; ++rr)
622 Vec2,deriv>::apply(rr,a,x[rr],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,DerivativeLayoutNS::derivative>,Vec2,deriv>
630 {
632 template <class Field>
633 static void apply(unsigned int r,const Field &a,
634 const Vec1 &x, Vec2 &y)
635 {
637 Vec2,deriv>::apply(r,a,x[0],y);
638 }
639 };
640 template <class F1,int dimD,
641 unsigned int d,
642 class Vec2,
643 unsigned int deriv>
644 struct LFETensorAxpy<Derivatives<F1,dimD,1,d,DerivativeLayoutNS::value>,Vec2,deriv>
645 {
647 template <class Field>
648 static void apply(unsigned int r,const Field &a,
649 const Vec1 &x, Vec2 &y)
650 {
651 typedef LFETensor<F1,dimD,deriv> LFETensorType;
652 const unsigned int rr = r*LFETensorType::size;
653 const FieldVector<F1,LFETensorType::size> &xx = x.template block<deriv>();
654 for (int i=0; i<FieldVector<F1,LFETensorType::size>::dimension; ++i)
655 y[rr+i] += xx[i]*a;
656 }
657 };
658
659 // ***********************************************
660 // Assign ****************************************
661 // ***********************************************
662 template <class Vec1,class Vec2>
664 {
665 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
666 {
667 field_cast(vec1,vec2);
668 }
669 };
670 template <int dimD,int dimR,unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout,
671 class F1,class F2>
672 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,layout>,
673 Derivatives<F2,dimD,dimR,deriv,layout> >
674 {
677 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
678 {
679 field_cast(vec1.block(),vec2.block());
680 }
681 };
682 template <int dimD,int dimR,unsigned int deriv,
683 class F1, class F2>
684 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value>,
685 Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> >
686 {
689 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
690 {
691 vec2.assign(vec1);
692 }
693 };
694 template <int dimD,int dimR,unsigned int deriv,
695 class F1, class F2>
696 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative>,
697 Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> >
698 {
701 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
702 {
703 vec2.assign(vec1);
704 }
705 };
706 template <int dimD,int dimR,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout,
707 class F1, class F2>
708 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
709 Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> >
710 {
713 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
714 {
715 vec2.assign(r,vec1);
716 }
717 };
718 template <int dimD,int dimR,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout,
719 class F1, class F2>
720 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
721 Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> >
722 {
725 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
726 {
727 vec2.assign(r,vec1);
728 }
729 };
730 template <int dimD,unsigned int deriv,
731 class F1, class F2>
732 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,
733 Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> >
734 {
737 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
738 {
739 field_cast(vec1.block(),vec2.block());
740 }
741 };
742 template <int dimD,unsigned int deriv,
743 class F1, class F2>
744 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,
745 Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> >
746 {
749 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
750 {
751 field_cast(vec1.block(),vec2.block());
752 }
753 };
754 template <int dimD,unsigned int deriv,
755 class F1, class F2>
756 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,
757 Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> >
758 {
761 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
762 {
763 field_cast(vec1.block(),vec2.block());
764 }
765 };
766 template <int dimD,unsigned int deriv,
767 class F1, class F2>
768 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,
769 Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> >
770 {
773 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
774 {
775 field_cast(vec1.block(),vec2.block());
776 }
777 };
778 template <int dimD,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout,
779 class F1, class F2>
780 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
781 F2 >
782 {
784 typedef F2 Vec2;
785 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
786 {
787 field_cast(vec1.block(),vec2);
788 }
789 };
790 template <int dimD,int dimR,
791 class F1,unsigned int deriv,
792 class F2>
793 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value>,FieldVector<F2,dimR> >
794 {
796 typedef FieldVector<F2,dimR> Vec2;
797 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
798 {
799 field_cast(vec1.template block<0>(),vec2);
800 }
801 };
802 template <int dimD,int dimR,
803 class F1,unsigned int deriv,
804 class F2>
805 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative>,FieldVector<F2,dimR> >
806 {
808 typedef FieldVector<F2,dimR> Vec2;
809 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
810 {
811 for (int rr=0; rr<dimR; ++rr)
812 field_cast(vec1[rr].template tensor<0>()[0].block(),vec2[rr]);
813 }
814 };
815 template <int dimD,
816 class F1,unsigned int deriv,
817 class F2,int dimR>
818 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,FieldVector<F2,dimR> >
819 {
821 typedef FieldVector<F2,dimR> Vec2;
822 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
823 {
824 field_cast(vec1.template tensor<0>()[0].block(),vec2[r]);
825 }
826 };
827 template <int dimD,
828 class F1,unsigned int deriv,
829 class F2,int dimR>
830 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,FieldVector<F2,dimR> >
831 {
833 typedef FieldVector<F2,dimR> Vec2;
834 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
835 {
836 field_cast(vec1[0].template tensor<0>()[0].block(),vec2[r]);
837 }
838 };
839 template <int dimD,
840 class F1,unsigned int deriv,
841 class F2>
842 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,FieldVector<F2,1> >
843 {
845 typedef FieldVector<F2,1> Vec2;
846 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
847 {
848 field_cast(vec1.template tensor<0>()[0].block(),vec2);
849 }
850 };
851 template <int dimD,
852 class F1,unsigned int deriv,
853 class F2>
854 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,FieldVector<F2,1> >
855 {
857 typedef FieldVector<F2,1> Vec2;
858 static void apply(unsigned int /*r*/,const Vec1 &vec1,Vec2 &vec2)
859 {
860 field_cast(vec1[0].template tensor<0>()[0].block(),vec2);
861 }
862 };
863
864 // ***********************************************
865 // IO ********************************************
866 // ***********************************************
867 template <class F,int dimD,unsigned int deriv>
868 std::ostream &operator<< ( std::ostream &out, const LFETensor< F,dimD,deriv > &tensor )
869 {
870 return out << tensor.block();
871 }
872#if 0
873 template <class F,int dimD,unsigned int deriv>
874 std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,deriv > &d )
875 {
876 out << static_cast<const ScalarDerivatives< F,dimD,deriv-1 > &>(d);
877 out << " , " << d.tensor() << std::endl;
878 return out;
879 }
880 template <class F,int dimD>
881 std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,0 > &d )
882 {
883 out << d.tensor() << std::endl;
884 return out;
885 }
886#endif
887 template <class F,int dimD,int dimR,unsigned int deriv>
889 {
890 out << " ( ";
891 out << d[0];
892 for (int r=1; r<dimR; ++r)
893 {
894 out << " , " << d[r];
895 }
896 out << " ) " << std::endl;
897 return out;
898 }
899 template <class F,int dimD,int dimR,unsigned int deriv>
901 {
902 out << static_cast<const Derivatives< F,dimD,dimR,deriv-1,DerivativeLayoutNS::value > &>(d);
903 out << " ( ";
904 out << d[0];
905 for (int r=1; r<dimR; ++r)
906 {
907 out << " , " << d[r];
908 }
909 out << " ) " << std::endl;
910 return out;
911 }
912 template <class F,int dimD,int dimR>
914 {
915 out << " ( ";
916 out << d[0];
917 for (int r=1; r<dimR; ++r)
918 {
919 out << " , " << d[r];
920 }
921 out << " ) " << std::endl;
922 return out;
923 }
924 template <class F,int dimD,int dimR>
925 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,DerivativeLayoutNS::value > &d )
926 {
927 out << " ( ";
928 out << d[0];
929 for (int r=1; r<dimR; ++r)
930 {
931 out << " , " << d[r];
932 }
933 out << " ) " << std::endl;
934 return out;
935 }
936 template <class F,int dimD,int dimR,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout>
937 std::ostream &operator<< ( std::ostream &out, const std::vector<Derivatives< F,dimD,dimR,deriv,layout > > &y )
938 {
939 out << "Number of basis functions: " << y.size() << std::endl;
940 for (unsigned int i=0; i<y.size(); ++i)
941 {
942 out << "Base " << i << " : " << std::endl;
943 out << y[i];
944 out << std::endl;
945 }
946 return out;
947 }
948}
949#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:151
DerivativeLayout
Definition tensor.hh:176
@ derivative
Definition tensor.hh:176
@ value
Definition tensor.hh:176
Definition tensor.hh:33
const Block & block() const
Definition tensor.hh:76
This & operator*=(const field_type &f)
Definition tensor.hh:56
Dune::FieldVector< F, size > Block
Definition tensor.hh:41
This & operator=(const FF &f)
Definition tensor.hh:44
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
Block & block()
Definition tensor.hh:72
void assign(const LFETensor< Fy, dimD, deriv > &y)
Definition tensor.hh:85
static const unsigned int size
Definition tensor.hh:40
const field_type & operator[](const unsigned int i) const
Definition tensor.hh:62
typename FieldTraits< field_type >::real_type real_type
Definition tensor.hh:97
Definition tensor.hh:115
Block & block()
Definition tensor.hh:162
F field_type
Definition tensor.hh:119
Block block_
Definition tensor.hh:170
void assign(const LFETensor< Fy, dimD, 0 > &y)
Definition tensor.hh:157
const Block & block() const
Definition tensor.hh:166
void axpy(const F &a, const This &y)
Definition tensor.hh:152
Dune::FieldVector< F, size > Block
Definition tensor.hh:121
Definition tensor.hh:180
typename FieldTraits< field_type >::real_type real_type
Definition tensor.hh:187
Derivatives< F, dimD, dimR, deriv, DerivativeLayoutNS::value > This
Definition tensor.hh:195
This & operator=(const Dune::FieldVector< LFETensor< F, dimD, dorder >, dimR > &t)
Definition tensor.hh:219
Derivatives< F, dimD, dimR, deriv-1, DerivativeLayoutNS::value > Base
Definition tensor.hh:196
const Dune::FieldVector< LFETensor< F, dimD, dorder >, dimR > & tensor(const std::integral_constant< int, dorder > &dorderVar) const
Definition tensor.hh:339
void assign(const FieldVector< Fy, size *dimRy > &y, unsigned int r)
Definition tensor.hh:317
void assign(const Derivatives< Fy, dimD, dimR, deriv, DerivativeLayoutNS::value > &y)
Definition tensor.hh:243
void assign(const Derivatives< Fy, dimD, dimR, dy, DerivativeLayoutNS::derivative > &y)
Definition tensor.hh:330
const ThisLFETensor & operator[](int r) const
Definition tensor.hh:312
Dune::FieldVector< F, LFETensor< F, dimD, dorder >::size *dimR > & block()
Definition tensor.hh:303
LFETensor< F, dimD, deriv > ThisLFETensor
Definition tensor.hh:197
void assign(unsigned int r, const Derivatives< Fy, dimD, 1, deriv, DerivativeLayoutNS::value > &y)
Definition tensor.hh:263
const Dune::FieldVector< F, LFETensor< F, dimD, dorder >::size *dimR > & block() const
Definition tensor.hh:296
Dune::FieldVector< LFETensor< F, dimD, dorder >, dimR > & tensor(const std::integral_constant< int, dorder > &dorderVar)
Definition tensor.hh:350
const Block & block() const
Definition tensor.hh:277
void axpy(const F &a, const This &y)
Definition tensor.hh:236
void assign(const Derivatives< Fy, dimD, dimR, deriv, DerivativeLayoutNS::derivative > &y)
Definition tensor.hh:249
This & operator=(const Dune::FieldVector< ThisLFETensor, dimR > &t)
Definition tensor.hh:213
void assign(unsigned int r, const Derivatives< Fy, dimD, 1, deriv, DerivativeLayoutNS::derivative > &y)
Definition tensor.hh:268
This & operator=(const Block &t)
Definition tensor.hh:224
Dune::FieldVector< F, size > Block
Definition tensor.hh:206
Dune::FieldVector< ThisLFETensor, dimR > tensor_
Definition tensor.hh:359
Dune::FieldVector< LFETensor< F, dimD, deriv >, dimR > & tensor(const std::integral_constant< int, deriv > &dorderVar)
Definition tensor.hh:355
Dune::FieldVector< LFETensor< F, dimD, dorder >, dimR > & tensor()
Definition tensor.hh:290
void assign(const Derivatives< Fy, dimD, dimRy, deriv, DerivativeLayoutNS::value > &y, unsigned int r)
Definition tensor.hh:257
const Dune::FieldVector< LFETensor< F, dimD, dorder >, dimR > & tensor() const
Definition tensor.hh:283
const Dune::FieldVector< LFETensor< F, dimD, deriv >, dimR > & tensor(const std::integral_constant< int, deriv > &dorderVar) const
Definition tensor.hh:344
ThisLFETensor & operator[](int r)
Definition tensor.hh:309
void assign(unsigned int r, const FieldVector< Fy, size/dimR > &y)
Definition tensor.hh:323
ThisLFETensor & operator[](int r)
Definition tensor.hh:442
LFETensor< F, dimD, 0 > ThisLFETensor
Definition tensor.hh:366
Dune::FieldVector< F, LFETensor< F, dimD, dorder >::size *dimR > & block()
Definition tensor.hh:465
Dune::FieldVector< LFETensor< F, dimD, 0 >, dimR > & tensor(const std::integral_constant< int, 0 > &dorderVar)
Definition tensor.hh:479
void assign(const Derivatives< Fy, dimD, dimR, 0, DerivativeLayoutNS::value > &y)
Definition tensor.hh:407
Derivatives< F, dimD, dimR, 0, DerivativeLayoutNS::value > This
Definition tensor.hh:365
void assign(const Derivatives< Fy, dimD, dimR, dy, DerivativeLayoutNS::derivative > &y)
Definition tensor.hh:484
void assign(const Derivatives< Fy, dimD, dimRy, 0, DerivativeLayoutNS::value > &y, unsigned int r)
Definition tensor.hh:418
const Dune::FieldVector< LFETensor< F, dimD, 0 >, dimR > & tensor() const
Definition tensor.hh:449
Dune::FieldVector< LFETensor< F, dimD, 0 >, dimR > & tensor()
Definition tensor.hh:453
This & operator=(const Block &t)
Definition tensor.hh:390
const ThisLFETensor & operator[](int r) const
Definition tensor.hh:445
const Dune::FieldVector< LFETensor< F, dimD, 0 >, dimR > & tensor(const std::integral_constant< int, 0 > &dorderVar) const
Definition tensor.hh:474
Dune::FieldVector< F, size > Block
Definition tensor.hh:375
const Block & block() const
Definition tensor.hh:437
const Dune::FieldVector< F, LFETensor< F, dimD, dorder >::size *dimR > & block() const
Definition tensor.hh:458
void assign(unsigned int r, const Derivatives< Fy, dimD, 1, 0, DerivativeLayoutNS::derivative > &y)
Definition tensor.hh:428
void assign(const Derivatives< Fy, dimD, dimR, 0, DerivativeLayoutNS::derivative > &y)
Definition tensor.hh:412
void assign(unsigned int r, const Derivatives< Fy, dimD, 1, 0, DerivativeLayoutNS::value > &y)
Definition tensor.hh:423
This & operator=(const Dune::FieldVector< ThisLFETensor, dimR > &t)
Definition tensor.hh:384
void assign(const FieldVector< Fy, size *dimRy > &y, unsigned int r)
Definition tensor.hh:490
Dune::FieldVector< ThisLFETensor, dimR > tensor_
Definition tensor.hh:499
void assign(unsigned int r, const FieldVector< Fy, size/dimR > &y)
Definition tensor.hh:495
This & operator=(const FF &f)
Definition tensor.hh:378
void axpy(const F &a, const This &y)
Definition tensor.hh:402
void assign(unsigned int r, const Derivatives< Fy, dimD, 1, deriv, layouty > &y)
Definition tensor.hh:556
void assign(const Derivatives< Fy, dimD, dimR, deriv, DerivativeLayoutNS::value > &y)
Definition tensor.hh:549
Derivatives< F, dimD, dimR, deriv, DerivativeLayoutNS::derivative > This
Definition tensor.hh:506
Derivatives< F, dimD, 1, deriv, DerivativeLayoutNS::value > ScalarDeriv
Definition tensor.hh:507
const ScalarDeriv & operator[](int r) const
Definition tensor.hh:573
void axpy(const FF &a, const This &y)
Definition tensor.hh:537
This & operator=(const Block &t)
Definition tensor.hh:524
Dune::FieldVector< F, size > Block
Definition tensor.hh:516
Dune::FieldVector< ScalarDeriv, dimR > deriv_
Definition tensor.hh:577
void assign(const Derivatives< Fy, dimD, dimR, deriv, DerivativeLayoutNS::derivative > &y)
Definition tensor.hh:543
Definition tensor.hh:585
static void apply(unsigned int r, const Field &a, const Vec1 &x, Vec2 &y)
Definition tensor.hh:587
static void apply(unsigned int r, const Field &a, const Vec1 &x, Vec2 &y)
Definition tensor.hh:601
Derivatives< F1, dimD, dimR, d, DerivativeLayoutNS::value > Vec1
Definition tensor.hh:599
Derivatives< F1, dimD, dimR, 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
Derivatives< F1, dimD, 1, d, DerivativeLayoutNS::derivative > Vec1
Definition tensor.hh:631
static void apply(unsigned int r, const Field &a, const Vec1 &x, Vec2 &y)
Definition tensor.hh:633
static void apply(unsigned int r, const Field &a, const Vec1 &x, Vec2 &y)
Definition tensor.hh:648
Derivatives< F1, dimD, 1, d, DerivativeLayoutNS::value > Vec1
Definition tensor.hh:646
Definition tensor.hh:664
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition tensor.hh:665
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition tensor.hh:677
Derivatives< F2, dimD, dimR, deriv, DerivativeLayoutNS::value > Vec2
Definition tensor.hh:712
Derivatives< F2, dimD, dimR, deriv, DerivativeLayoutNS::derivative > Vec2
Definition tensor.hh:724
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition tensor.hh:785
Derivatives< F1, dimD, 1, deriv, layout > Vec1
Definition tensor.hh:783
Derivatives< F1, dimD, dimR, deriv, DerivativeLayoutNS::value > Vec1
Definition tensor.hh:795
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition tensor.hh:797
Derivatives< F1, dimD, dimR, deriv, DerivativeLayoutNS::derivative > Vec1
Definition tensor.hh:807
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition tensor.hh:809
Derivatives< F1, dimD, 1, deriv, DerivativeLayoutNS::value > Vec1
Definition tensor.hh:820
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition tensor.hh:822
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition tensor.hh:834
Derivatives< F1, dimD, 1, deriv, DerivativeLayoutNS::derivative > Vec1
Definition tensor.hh:832
Derivatives< F1, dimD, 1, deriv, DerivativeLayoutNS::value > Vec1
Definition tensor.hh:844
static void apply(unsigned int r, const Vec1 &vec1, Vec2 &vec2)
Definition tensor.hh:846
Derivatives< F1, dimD, 1, deriv, DerivativeLayoutNS::derivative > Vec1
Definition tensor.hh:856
static void apply(unsigned int, const Vec1 &vec1, Vec2 &vec2)
Definition tensor.hh:858