Dune Core Modules (2.7.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 namespace DerivativeLayoutNS {
166 enum DerivativeLayout {value,derivative};
167 }
168 template <class F,int dimD,int dimR,unsigned int deriv,
169 DerivativeLayoutNS::DerivativeLayout layout>
170 struct Derivatives;
171
172 // Implemnetation for valued based layout
173 template <class F,int dimD,int dimR,unsigned int deriv>
174 struct Derivatives<F,dimD,dimR,deriv,DerivativeLayoutNS::value>
175 : public Derivatives<F,dimD,dimR,deriv-1,DerivativeLayoutNS::value>
176 {
177 typedef Derivatives<F,dimD,dimR,deriv,DerivativeLayoutNS::value> This;
178 typedef Derivatives<F,dimD,dimR,deriv-1,DerivativeLayoutNS::value> Base;
179 typedef LFETensor<F,dimD,deriv> ThisLFETensor;
180
181 typedef F Field;
182 typedef F field_type;
183
184 static const DerivativeLayoutNS::DerivativeLayout layout = DerivativeLayoutNS::value;
185 static const unsigned int dimDomain = dimD;
186 static const unsigned int dimRange = dimR;
187 // size needs to be an anonymous enum value for gcc 3.4 compatibility
188 enum { size = Base::size+ThisLFETensor::size*dimR };
189 typedef Dune::FieldVector<F,size> Block;
190
191 This &operator=(const F& f)
192 {
193 block() = f;
194 return *this;
195 }
196 This &operator=(const Dune::FieldVector<ThisLFETensor,dimR> &t)
197 {
198 tensor_ = t;
199 return *this;
200 }
201 template <unsigned int dorder>
202 This &operator=(const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &t)
203 {
204 tensor<dorder>() = t;
205 return *this;
206 }
207 This &operator=(const Block &t)
208 {
209 block() = t;
210 return *this;
211 }
212
213 This &operator*= ( const field_type &f )
214 {
215 block() *= f;
216 return *this;
217 }
218
219 void axpy(const F &a, const This &y)
220 {
221 block().axpy(a,y.block());
222 }
223
224 // assign with same layout (only different Field)
225 template <class Fy>
226 void assign(const Derivatives<Fy,dimD,dimR,deriv,DerivativeLayoutNS::value> &y)
227 {
228 field_cast(y.block(),block());
229 }
230 // assign with different layout (same dimRange)
231 template <class Fy>
232 void assign(const Derivatives<Fy,dimD,dimR,deriv,DerivativeLayoutNS::derivative> &y)
233 {
234 Base::assign(y);
235 for (int rr=0; rr<dimR; ++rr)
236 tensor_[rr] = y[rr].template tensor<deriv>()[0];
237 }
238 // assign with rth component of function
239 template <class Fy,int dimRy>
240 void assign(const Derivatives<Fy,dimD,dimRy,deriv,DerivativeLayoutNS::value> &y,unsigned int r)
241 {
242 assign<Fy,dimRy>(y.block(),r);
243 }
244 // assign with scalar functions to component r
245 template <class Fy>
246 void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,DerivativeLayoutNS::value> &y)
247 {
248 assign(r,y.block());
249 }
250 template <class Fy>
251 void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,DerivativeLayoutNS::derivative> &y)
252 {
253 assign(r,y[0]);
254 }
255
256 Block &block()
257 {
258 return reinterpret_cast<Block&>(*this);
259 }
260 const Block &block() const
261 {
262 return reinterpret_cast<const Block&>(*this);
263 }
264
265 template <unsigned int dorder>
266 const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor() const
267 {
268 // use integral_constant<int,...> here to stay compatible with Int2Type
269 const std::integral_constant<int,dorder> a = {};
270 return tensor(a);
271 }
272 template <unsigned int dorder>
273 Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor()
274 {
275 // use integral_constant<int,...> here to stay compatible with Int2Type
276 return tensor(std::integral_constant<int,dorder>());
277 }
278 template <unsigned int dorder>
279 const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block() const
280 {
281 // use integral_constant<int,...> here to stay compatible with Int2Type
282 const std::integral_constant<int,dorder> a = {};
283 return reinterpret_cast<const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
284 }
285 template <unsigned int dorder>
286 Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
287 {
288 // use integral_constant<int,...> here to stay compatible with Int2Type
289 const std::integral_constant<int,dorder> a = {};
290 return reinterpret_cast<Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
291 }
292 ThisLFETensor &operator[](int r) {
293 return tensor_[r];
294 }
295 const ThisLFETensor &operator[](int r) const {
296 return tensor_[r];
297 }
298 protected:
299 template <class Fy,int dimRy>
300 void assign(const FieldVector<Fy,size*dimRy> &y,unsigned int r)
301 {
302 Base::template assign<Fy,dimRy>(reinterpret_cast<const FieldVector<Fy,Base::size*dimRy>&>(y),r);
303 tensor_[0] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[Base::size*dimRy+r*ThisLFETensor::size]);
304 }
305 template <class Fy>
306 void assign(unsigned int r,const FieldVector<Fy,size/dimR> &y)
307 {
308 Base::assign(r,reinterpret_cast<const FieldVector<Fy,Base::size/dimR>&>(y));
309 tensor_[r] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[Base::size/dimR]);
310 }
311 // assign with different layout (same dimRange)
312 template <class Fy,unsigned int dy>
313 void assign(const Derivatives<Fy,dimD,dimR,dy,DerivativeLayoutNS::derivative> &y)
314 {
315 Base::assign(y);
316 for (int rr=0; rr<dimR; ++rr)
317 tensor_[rr] = y[rr].template tensor<deriv>()[0];
318 }
319
320 template <int dorder>
321 const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
322 tensor(const std::integral_constant<int,dorder> &dorderVar) const
323 {
324 return Base::tensor(dorderVar);
325 }
326 const Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
327 tensor(const std::integral_constant<int,deriv> &dorderVar) const
328 {
329 return tensor_;
330 }
331 template <int dorder>
332 Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
333 tensor(const std::integral_constant<int,dorder> &dorderVar)
334 {
335 return Base::tensor(dorderVar);
336 }
337 Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
338 tensor(const std::integral_constant<int,deriv> &dorderVar)
339 {
340 return tensor_;
341 }
342 Dune::FieldVector<ThisLFETensor,dimR> tensor_;
343 };
344
345 template <class F,int dimD,int dimR>
346 struct Derivatives<F,dimD,dimR,0,DerivativeLayoutNS::value>
347 {
348 typedef Derivatives<F,dimD,dimR,0,DerivativeLayoutNS::value> This;
349 typedef LFETensor<F,dimD,0> ThisLFETensor;
350
351 typedef F Field;
352 typedef F field_type;
353
354 static const DerivativeLayoutNS::DerivativeLayout layout = DerivativeLayoutNS::value;
355 static const unsigned int dimDomain = dimD;
356 static const unsigned int dimRange = dimR;
357 // size needs to be an anonymous enum value for gcc 3.4 compatibility
358 enum { size = ThisLFETensor::size*dimR };
359 typedef Dune::FieldVector<F,size> Block;
360
361 template <class FF>
362 This &operator=(const FF& f)
363 {
364 for (int r=0; r<dimR; ++r)
365 tensor_[r] = field_cast<F>(f);
366 return *this;
367 }
368 This &operator=(const Dune::FieldVector<ThisLFETensor,dimR> &t)
369 {
370 tensor_ = t;
371 return *this;
372 }
373
374 This &operator=(const Block &t)
375 {
376 block() = t;
377 return *this;
378 }
379
380 This &operator*= ( const field_type &f )
381 {
382 block() *= f;
383 return *this;
384 }
385
386 void axpy(const F &a, const This &y)
387 {
388 block().axpy(a,y.block());
389 }
390 template <class Fy>
391 void assign(const Derivatives<Fy,dimD,dimR,0,DerivativeLayoutNS::value> &y)
392 {
393 field_cast(y.block(),block());
394 }
395 template <class Fy>
396 void assign(const Derivatives<Fy,dimD,dimR,0,DerivativeLayoutNS::derivative> &y)
397 {
398 for (int rr=0; rr<dimR; ++rr)
399 tensor_[rr] = y[rr].template tensor<0>()[0];
400 }
401 template <class Fy,int dimRy>
402 void assign(const Derivatives<Fy,dimD,dimRy,0,DerivativeLayoutNS::value> &y,unsigned int r)
403 {
404 assign<Fy,dimRy>(y.block(),r);
405 }
406 template <class Fy>
407 void assign(unsigned int r,const Derivatives<Fy,dimD,1,0,DerivativeLayoutNS::value> &y)
408 {
409 tensor_[r].assign(y[0]);
410 }
411 template <class Fy>
412 void assign(unsigned int r,const Derivatives<Fy,dimD,1,0,DerivativeLayoutNS::derivative> &y)
413 {
414 tensor_[r].assign(y[0][0]);
415 }
416
417 Block &block()
418 {
419 return reinterpret_cast<Block&>(*this);
420 }
421 const Block &block() const
422 {
423 return reinterpret_cast<const Block&>(*this);
424 }
425
426 ThisLFETensor &operator[](int r) {
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>
468 void assign(const Derivatives<Fy,dimD,dimR,dy,DerivativeLayoutNS::derivative> &y)
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 {
490 typedef Derivatives<F,dimD,dimR,deriv,DerivativeLayoutNS::derivative> This;
491 typedef Derivatives<F,dimD,1,deriv,DerivativeLayoutNS::value> ScalarDeriv;
492
493 typedef F Field;
494 typedef F field_type;
495
496 static const DerivativeLayoutNS::DerivativeLayout layout = DerivativeLayoutNS::value;
497 static const unsigned int dimDomain = dimD;
498 static const unsigned int dimRange = dimR;
499 // size needs to be an anonymous enum value for gcc 3.4 compatibility
500 enum { size = ScalarDeriv::size*dimR };
501 typedef Dune::FieldVector<F,size> Block;
502
503 template <class FF>
504 This &operator=(const FF& f)
505 {
506 block() = field_cast<F>(f);
507 return *this;
508 }
509 This &operator=(const Block &t)
510 {
511 block() = t;
512 return *this;
513 }
514
515 This &operator*= ( const field_type &f )
516 {
517 block() *= f;
518 return *this;
519 }
520
521 template <class FF>
522 void axpy(const FF &a, const This &y)
523 {
524 block().axpy(field_cast<F>(a),y.block());
525 }
526 // assign with same layout (only different Field)
527 template <class Fy>
528 void assign(const Derivatives<Fy,dimD,dimR,deriv,DerivativeLayoutNS::derivative> &y)
529 {
530 field_cast(y.block(),block());
531 }
532 // assign with different layout (same dimRange)
533 template <class Fy>
534 void assign(const Derivatives<Fy,dimD,dimR,deriv,DerivativeLayoutNS::value> &y)
535 {
536 for (unsigned int rr=0; rr<dimR; ++rr)
537 deriv_[rr].assign(y,rr);
538 }
539 // assign with scalar functions to component r
540 template <class Fy,DerivativeLayoutNS::DerivativeLayout layouty>
541 void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,layouty> &y)
542 {
543 deriv_[r].assign(r,y);
544 }
545
546 Block &block()
547 {
548 return reinterpret_cast<Block&>(*this);
549 }
550 const Block &block() const
551 {
552 return reinterpret_cast<const Block&>(*this);
553 }
554
555 ScalarDeriv &operator[](int r) {
556 return deriv_[r];
557 }
558 const ScalarDeriv &operator[](int r) const {
559 return deriv_[r];
560 }
561 protected:
562 Dune::FieldVector<ScalarDeriv,dimR> deriv_;
563 };
564
565 // ******************************************
566 // AXPY *************************************
567 // ******************************************
568 template <class Vec1,class Vec2,unsigned int deriv>
569 struct LFETensorAxpy
570 {
571 template <class Field>
572 static void apply(unsigned int r,const Field &a,
573 const Vec1 &x, Vec2 &y)
574 {
575 y.axpy(a,x);
576 }
577 };
578 template <class F1,int dimD,int dimR,
579 unsigned int d,
580 class Vec2,
581 unsigned int deriv>
582 struct LFETensorAxpy<Derivatives<F1,dimD,dimR,d,DerivativeLayoutNS::value>,Vec2,deriv>
583 {
584 typedef Derivatives<F1,dimD,dimR,d,DerivativeLayoutNS::value> Vec1;
585 template <class Field>
586 static void apply(unsigned int r,const Field &a,
587 const Vec1 &x, Vec2 &y)
588 {
589 const FieldVector<F1,Vec2::size> &xx = x.template block<deriv>();
590 for (int i=0; i<y.size; ++i)
591 y[i] += xx[i]*a;
592 }
593 };
594 template <class F1,int dimD,int dimR,
595 unsigned int d,
596 class Vec2,
597 unsigned int deriv>
598 struct LFETensorAxpy<Derivatives<F1,dimD,dimR,d,DerivativeLayoutNS::derivative>,Vec2,deriv>
599 {
600 typedef Derivatives<F1,dimD,dimR,d,DerivativeLayoutNS::derivative> Vec1;
601 template <class Field>
602 static void apply(unsigned int r,const Field &a,
603 const Vec1 &x, Vec2 &y)
604 {
605 for (int rr=0; rr<dimR; ++rr)
606 LFETensorAxpy<Derivatives<F1,dimD,1,d,DerivativeLayoutNS::value>,
607 Vec2,deriv>::apply(rr,a,x[rr],y);
608 }
609 };
610 template <class F1,int dimD,
611 unsigned int d,
612 class Vec2,
613 unsigned int deriv>
614 struct LFETensorAxpy<Derivatives<F1,dimD,1,d,DerivativeLayoutNS::derivative>,Vec2,deriv>
615 {
616 typedef Derivatives<F1,dimD,1,d,DerivativeLayoutNS::derivative> Vec1;
617 template <class Field>
618 static void apply(unsigned int r,const Field &a,
619 const Vec1 &x, Vec2 &y)
620 {
621 LFETensorAxpy<Derivatives<F1,dimD,1,d,DerivativeLayoutNS::value>,
622 Vec2,deriv>::apply(r,a,x[0],y);
623 }
624 };
625 template <class F1,int dimD,
626 unsigned int d,
627 class Vec2,
628 unsigned int deriv>
629 struct LFETensorAxpy<Derivatives<F1,dimD,1,d,DerivativeLayoutNS::value>,Vec2,deriv>
630 {
631 typedef Derivatives<F1,dimD,1,d,DerivativeLayoutNS::value> Vec1;
632 template <class Field>
633 static void apply(unsigned int r,const Field &a,
634 const Vec1 &x, Vec2 &y)
635 {
636 typedef LFETensor<F1,dimD,deriv> LFETensorType;
637 const unsigned int rr = r*LFETensorType::size;
638 const FieldVector<F1,LFETensorType::size> &xx = x.template block<deriv>();
639 for (int i=0; i<FieldVector<F1,LFETensorType::size>::dimension; ++i)
640 y[rr+i] += xx[i]*a;
641 }
642 };
643
644 // ***********************************************
645 // Assign ****************************************
646 // ***********************************************
647 template <class Vec1,class Vec2>
648 struct DerivativeAssign
649 {
650 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
651 {
652 field_cast(vec1,vec2);
653 }
654 };
655 template <int dimD,int dimR,unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout,
656 class F1,class F2>
657 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,layout>,
658 Derivatives<F2,dimD,dimR,deriv,layout> >
659 {
660 typedef Derivatives<F1,dimD,dimR,deriv,layout> Vec1;
661 typedef Derivatives<F2,dimD,dimR,deriv,layout> Vec2;
662 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
663 {
664 field_cast(vec1.block(),vec2.block());
665 }
666 };
667 template <int dimD,int dimR,unsigned int deriv,
668 class F1, class F2>
669 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value>,
670 Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> >
671 {
672 typedef Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value> Vec1;
673 typedef Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> Vec2;
674 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
675 {
676 vec2.assign(vec1);
677 }
678 };
679 template <int dimD,int dimR,unsigned int deriv,
680 class F1, class F2>
681 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative>,
682 Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> >
683 {
684 typedef Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative> Vec1;
685 typedef Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> Vec2;
686 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
687 {
688 vec2.assign(vec1);
689 }
690 };
691 template <int dimD,int dimR,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout,
692 class F1, class F2>
693 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
694 Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> >
695 {
696 typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
697 typedef Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> Vec2;
698 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
699 {
700 vec2.assign(r,vec1);
701 }
702 };
703 template <int dimD,int dimR,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout,
704 class F1, class F2>
705 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
706 Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> >
707 {
708 typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
709 typedef Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> Vec2;
710 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
711 {
712 vec2.assign(r,vec1);
713 }
714 };
715 template <int dimD,unsigned int deriv,
716 class F1, class F2>
717 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,
718 Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> >
719 {
720 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value> Vec1;
721 typedef Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> Vec2;
722 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
723 {
724 field_cast(vec1.block(),vec2.block());
725 }
726 };
727 template <int dimD,unsigned int deriv,
728 class F1, class F2>
729 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,
730 Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> >
731 {
732 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec1;
733 typedef Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec2;
734 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
735 {
736 field_cast(vec1.block(),vec2.block());
737 }
738 };
739 template <int dimD,unsigned int deriv,
740 class F1, class F2>
741 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,
742 Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> >
743 {
744 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec1;
745 typedef Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> Vec2;
746 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
747 {
748 field_cast(vec1.block(),vec2.block());
749 }
750 };
751 template <int dimD,unsigned int deriv,
752 class F1, class F2>
753 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,
754 Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> >
755 {
756 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value> Vec1;
757 typedef Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec2;
758 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
759 {
760 field_cast(vec1.block(),vec2.block());
761 }
762 };
763 template <int dimD,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout,
764 class F1, class F2>
765 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
766 F2 >
767 {
768 typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
769 typedef F2 Vec2;
770 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
771 {
772 field_cast(vec1.block(),vec2);
773 }
774 };
775 template <int dimD,int dimR,
776 class F1,unsigned int deriv,
777 class F2>
778 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value>,FieldVector<F2,dimR> >
779 {
780 typedef Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value> Vec1;
781 typedef FieldVector<F2,dimR> Vec2;
782 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
783 {
784 field_cast(vec1.template block<0>(),vec2);
785 }
786 };
787 template <int dimD,int dimR,
788 class F1,unsigned int deriv,
789 class F2>
790 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative>,FieldVector<F2,dimR> >
791 {
792 typedef Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative> Vec1;
793 typedef FieldVector<F2,dimR> Vec2;
794 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
795 {
796 for (int rr=0; rr<dimR; ++rr)
797 field_cast(vec1[rr].template tensor<0>()[0].block(),vec2[rr]);
798 }
799 };
800 template <int dimD,
801 class F1,unsigned int deriv,
802 class F2,int dimR>
803 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,FieldVector<F2,dimR> >
804 {
805 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value> Vec1;
806 typedef FieldVector<F2,dimR> Vec2;
807 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
808 {
809 field_cast(vec1.template tensor<0>()[0].block(),vec2[r]);
810 }
811 };
812 template <int dimD,
813 class F1,unsigned int deriv,
814 class F2,int dimR>
815 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,FieldVector<F2,dimR> >
816 {
817 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec1;
818 typedef FieldVector<F2,dimR> Vec2;
819 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
820 {
821 field_cast(vec1[0].template tensor<0>()[0].block(),vec2[r]);
822 }
823 };
824 template <int dimD,
825 class F1,unsigned int deriv,
826 class F2>
827 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,FieldVector<F2,1> >
828 {
829 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value> Vec1;
830 typedef FieldVector<F2,1> Vec2;
831 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
832 {
833 field_cast(vec1.template tensor<0>()[0].block(),vec2);
834 }
835 };
836 template <int dimD,
837 class F1,unsigned int deriv,
838 class F2>
839 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,FieldVector<F2,1> >
840 {
841 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec1;
842 typedef FieldVector<F2,1> Vec2;
843 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
844 {
845 field_cast(vec1[0].template tensor<0>()[0].block(),vec2);
846 }
847 };
848
849 // ***********************************************
850 // IO ********************************************
851 // ***********************************************
852 template <class F,int dimD,unsigned int deriv>
853 std::ostream &operator<< ( std::ostream &out, const LFETensor< F,dimD,deriv > &tensor )
854 {
855 return out << tensor.block();
856 }
857#if 0
858 template <class F,int dimD,unsigned int deriv>
859 std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,deriv > &d )
860 {
861 out << static_cast<const ScalarDerivatives< F,dimD,deriv-1 > &>(d);
862 out << " , " << d.tensor() << std::endl;
863 return out;
864 }
865 template <class F,int dimD>
866 std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,0 > &d )
867 {
868 out << d.tensor() << std::endl;
869 return out;
870 }
871#endif
872 template <class F,int dimD,int dimR,unsigned int deriv>
873 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,DerivativeLayoutNS::derivative > &d )
874 {
875 out << " ( ";
876 out << d[0];
877 for (int r=1; r<dimR; ++r)
878 {
879 out << " , " << d[r];
880 }
881 out << " ) " << std::endl;
882 return out;
883 }
884 template <class F,int dimD,int dimR,unsigned int deriv>
885 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,DerivativeLayoutNS::value > &d )
886 {
887 out << static_cast<const Derivatives< F,dimD,dimR,deriv-1,DerivativeLayoutNS::value > &>(d);
888 out << " ( ";
889 out << d[0];
890 for (int r=1; r<dimR; ++r)
891 {
892 out << " , " << d[r];
893 }
894 out << " ) " << std::endl;
895 return out;
896 }
897 template <class F,int dimD,int dimR>
898 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,DerivativeLayoutNS::derivative > &d )
899 {
900 out << " ( ";
901 out << d[0];
902 for (int r=1; r<dimR; ++r)
903 {
904 out << " , " << d[r];
905 }
906 out << " ) " << std::endl;
907 return out;
908 }
909 template <class F,int dimD,int dimR>
910 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,DerivativeLayoutNS::value > &d )
911 {
912 out << " ( ";
913 out << d[0];
914 for (int r=1; r<dimR; ++r)
915 {
916 out << " , " << d[r];
917 }
918 out << " ) " << std::endl;
919 return out;
920 }
921 template <class F,int dimD,int dimR,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout>
922 std::ostream &operator<< ( std::ostream &out, const std::vector<Derivatives< F,dimD,dimR,deriv,layout > > &y )
923 {
924 out << "Number of basis functions: " << y.size() << std::endl;
925 for (unsigned int i=0; i<y.size(); ++i)
926 {
927 out << "Base " << i << " : " << std::endl;
928 out << y[i];
929 out << std::endl;
930 }
931 return out;
932 }
933}
934#endif // DUNE_TENSOR_HH
Dune::DenseVector::axpy
derived_type & axpy(const field_type &a, const DenseVector< Other > &x)
vector space axpy operation ( *this += a x )
Definition: densevector.hh:576
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:46
Dune
Dune namespace.
Definition: alignedallocator.hh:14
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:445
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)