Dune Core Modules (2.6.0)

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