Dune Core Modules (2.9.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// SPDX-FileCopyrightInfo: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5
6#ifndef DUNE_TENSOR_HH
7#define DUNE_TENSOR_HH
8
9#include <ostream>
10#include <vector>
11
12#include <dune/common/fvector.hh>
13
14#include <dune/localfunctions/utility/field.hh>
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 rewritting...
28 ***************************************************/
29
30 // Structure for scalar tensor of order deriv
31 template <class F,int dimD,unsigned int deriv>
32 class LFETensor
33 {
34 typedef LFETensor<F,dimD,deriv> This;
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
56 This &operator*= ( const field_type &f )
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
72 Block &block()
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>
85 void assign(const LFETensor<Fy,dimD,deriv> &y)
86 {
87 field_cast(y.block(),block());
88 }
89 Block block_;
90 };
91
92 // ******************************************
93 template <class F,unsigned int deriv>
94 struct LFETensor<F,0,deriv>
95 {
96 static const int size = 0;
97 };
98
99 template <class F>
100 struct LFETensor<F,0,0>
101 {
102 static const int size = 1;
103 };
104
105 template <class F,int dimD>
106 class LFETensor<F,dimD,0>
107 {
108 typedef LFETensor<F,dimD,0> This;
109
110 public:
111 typedef F field_type;
112 static const int size = 1;
113 typedef Dune::FieldVector<F,size> Block;
114
115 template< class FF >
116 This &operator= ( const FF &f )
117 {
118 block() = field_cast< F >( f );
119 return *this;
120 }
121
122 This &operator= ( const Block &b )
123 {
124 block() = b;
125 return *this;
126 }
127
128 This &operator*= ( const field_type &f )
129 {
130 block() *= f;
131 return *this;
132 }
133
134 const F &operator[] ( const unsigned int i ) const
135 {
136 return block()[ i ];
137 }
138
139 F &operator[] ( const unsigned int i )
140 {
141 return block()[ i ];
142 }
143
144 void axpy(const F& a, const This &y)
145 {
146 block().axpy(a,y.block());
147 }
148 template <class Fy>
149 void assign(const LFETensor<Fy,dimD,0> &y)
150 {
151 field_cast(y.block(),block());
152 }
153
154 Block &block()
155 {
156 return block_;
157 }
158 const Block &block() const
159 {
160 return block_;
161 }
162 Block block_;
163 };
164 // ***********************************************************
165 // Structure for all derivatives up to order deriv
166 // for vector valued function
167 namespace DerivativeLayoutNS {
168 enum DerivativeLayout {value,derivative};
169 }
170 template <class F,int dimD,int dimR,unsigned int deriv,
171 DerivativeLayoutNS::DerivativeLayout layout>
172 struct Derivatives;
173
174 // Implemnetation for valued based layout
175 template <class F,int dimD,int dimR,unsigned int deriv>
176 struct Derivatives<F,dimD,dimR,deriv,DerivativeLayoutNS::value>
177 : public Derivatives<F,dimD,dimR,deriv-1,DerivativeLayoutNS::value>
178 {
179 typedef Derivatives<F,dimD,dimR,deriv,DerivativeLayoutNS::value> This;
180 typedef Derivatives<F,dimD,dimR,deriv-1,DerivativeLayoutNS::value> Base;
181 typedef LFETensor<F,dimD,deriv> ThisLFETensor;
182
183 typedef F Field;
184 typedef F field_type;
185
186 static const DerivativeLayoutNS::DerivativeLayout layout = DerivativeLayoutNS::value;
187 static const unsigned int dimDomain = dimD;
188 static const unsigned int dimRange = dimR;
189 constexpr static int size = Base::size+ThisLFETensor::size*dimR;
190 typedef Dune::FieldVector<F,size> Block;
191
192 This &operator=(const F& f)
193 {
194 block() = f;
195 return *this;
196 }
197 This &operator=(const Dune::FieldVector<ThisLFETensor,dimR> &t)
198 {
199 tensor_ = t;
200 return *this;
201 }
202 template <unsigned int dorder>
203 This &operator=(const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &t)
204 {
205 tensor<dorder>() = t;
206 return *this;
207 }
208 This &operator=(const Block &t)
209 {
210 block() = t;
211 return *this;
212 }
213
214 This &operator*= ( const field_type &f )
215 {
216 block() *= f;
217 return *this;
218 }
219
220 void axpy(const F &a, const This &y)
221 {
222 block().axpy(a,y.block());
223 }
224
225 // assign with same layout (only different Field)
226 template <class Fy>
227 void assign(const Derivatives<Fy,dimD,dimR,deriv,DerivativeLayoutNS::value> &y)
228 {
229 field_cast(y.block(),block());
230 }
231 // assign with different layout (same dimRange)
232 template <class Fy>
233 void assign(const Derivatives<Fy,dimD,dimR,deriv,DerivativeLayoutNS::derivative> &y)
234 {
235 Base::assign(y);
236 for (int rr=0; rr<dimR; ++rr)
237 tensor_[rr] = y[rr].template tensor<deriv>()[0];
238 }
239 // assign with rth component of function
240 template <class Fy,int dimRy>
241 void assign(const Derivatives<Fy,dimD,dimRy,deriv,DerivativeLayoutNS::value> &y,unsigned int r)
242 {
243 assign<Fy,dimRy>(y.block(),r);
244 }
245 // assign with scalar functions to component r
246 template <class Fy>
247 void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,DerivativeLayoutNS::value> &y)
248 {
249 assign(r,y.block());
250 }
251 template <class Fy>
252 void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,DerivativeLayoutNS::derivative> &y)
253 {
254 assign(r,y[0]);
255 }
256
257 Block &block()
258 {
259 return reinterpret_cast<Block&>(*this);
260 }
261 const Block &block() const
262 {
263 return reinterpret_cast<const Block&>(*this);
264 }
265
266 template <unsigned int dorder>
267 const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor() const
268 {
269 // use integral_constant<int,...> here to stay compatible with Int2Type
270 const std::integral_constant<int,dorder> a = {};
271 return tensor(a);
272 }
273 template <unsigned int dorder>
274 Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &tensor()
275 {
276 // use integral_constant<int,...> here to stay compatible with Int2Type
277 return tensor(std::integral_constant<int,dorder>());
278 }
279 template <unsigned int dorder>
280 const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block() const
281 {
282 // use integral_constant<int,...> here to stay compatible with Int2Type
283 const std::integral_constant<int,dorder> a = {};
284 return reinterpret_cast<const Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
285 }
286 template <unsigned int dorder>
287 Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR> &block()
288 {
289 // use integral_constant<int,...> here to stay compatible with Int2Type
290 const std::integral_constant<int,dorder> a = {};
291 return reinterpret_cast<Dune::FieldVector<F,LFETensor<F,dimD,dorder>::size*dimR>&>(tensor(a));
292 }
293 ThisLFETensor &operator[](int r) {
294 return tensor_[r];
295 }
296 const ThisLFETensor &operator[](int r) const {
297 return tensor_[r];
298 }
299 protected:
300 template <class Fy,int dimRy>
301 void assign(const FieldVector<Fy,size*dimRy> &y,unsigned int r)
302 {
303 Base::template assign<Fy,dimRy>(reinterpret_cast<const FieldVector<Fy,Base::size*dimRy>&>(y),r);
304 tensor_[0] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[Base::size*dimRy+r*ThisLFETensor::size]);
305 }
306 template <class Fy>
307 void assign(unsigned int r,const FieldVector<Fy,size/dimR> &y)
308 {
309 Base::assign(r,reinterpret_cast<const FieldVector<Fy,Base::size/dimR>&>(y));
310 tensor_[r] = reinterpret_cast<const FieldVector<Fy,ThisLFETensor::size>&>(y[Base::size/dimR]);
311 }
312 // assign with different layout (same dimRange)
313 template <class Fy,unsigned int dy>
314 void assign(const Derivatives<Fy,dimD,dimR,dy,DerivativeLayoutNS::derivative> &y)
315 {
316 Base::assign(y);
317 for (int rr=0; rr<dimR; ++rr)
318 tensor_[rr] = y[rr].template tensor<deriv>()[0];
319 }
320
321 template <int dorder>
322 const Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
323 tensor(const std::integral_constant<int,dorder> &dorderVar) const
324 {
325 return Base::tensor(dorderVar);
326 }
327 const Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
328 tensor(const std::integral_constant<int,deriv> &dorderVar) const
329 {
330 return tensor_;
331 }
332 template <int dorder>
333 Dune::FieldVector<LFETensor<F,dimD,dorder>,dimR> &
334 tensor(const std::integral_constant<int,dorder> &dorderVar)
335 {
336 return Base::tensor(dorderVar);
337 }
338 Dune::FieldVector<LFETensor<F,dimD,deriv>,dimR> &
339 tensor(const std::integral_constant<int,deriv> &dorderVar)
340 {
341 return tensor_;
342 }
343 Dune::FieldVector<ThisLFETensor,dimR> tensor_;
344 };
345
346 template <class F,int dimD,int dimR>
347 struct Derivatives<F,dimD,dimR,0,DerivativeLayoutNS::value>
348 {
349 typedef Derivatives<F,dimD,dimR,0,DerivativeLayoutNS::value> This;
350 typedef LFETensor<F,dimD,0> ThisLFETensor;
351
352 typedef F Field;
353 typedef F field_type;
354
355 static const DerivativeLayoutNS::DerivativeLayout layout = DerivativeLayoutNS::value;
356 static const unsigned int dimDomain = dimD;
357 static const unsigned int dimRange = dimR;
358 constexpr static int size = ThisLFETensor::size*dimR;
359 typedef Dune::FieldVector<F,size> Block;
360
361 template <class FF>
362 This &operator=(const FF& f)
363 {
364 for (int r=0; r<dimR; ++r)
365 tensor_[r] = field_cast<F>(f);
366 return *this;
367 }
368 This &operator=(const Dune::FieldVector<ThisLFETensor,dimR> &t)
369 {
370 tensor_ = t;
371 return *this;
372 }
373
374 This &operator=(const Block &t)
375 {
376 block() = t;
377 return *this;
378 }
379
380 This &operator*= ( const field_type &f )
381 {
382 block() *= f;
383 return *this;
384 }
385
386 void axpy(const F &a, const This &y)
387 {
388 block().axpy(a,y.block());
389 }
390 template <class Fy>
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 constexpr static int size = ScalarDeriv::size*dimR;
500 typedef Dune::FieldVector<F,size> Block;
501
502 template <class FF>
503 This &operator=(const FF& f)
504 {
505 block() = field_cast<F>(f);
506 return *this;
507 }
508 This &operator=(const Block &t)
509 {
510 block() = t;
511 return *this;
512 }
513
514 This &operator*= ( const field_type &f )
515 {
516 block() *= f;
517 return *this;
518 }
519
520 template <class FF>
521 void axpy(const FF &a, const This &y)
522 {
523 block().axpy(field_cast<F>(a),y.block());
524 }
525 // assign with same layout (only different Field)
526 template <class Fy>
527 void assign(const Derivatives<Fy,dimD,dimR,deriv,DerivativeLayoutNS::derivative> &y)
528 {
529 field_cast(y.block(),block());
530 }
531 // assign with different layout (same dimRange)
532 template <class Fy>
533 void assign(const Derivatives<Fy,dimD,dimR,deriv,DerivativeLayoutNS::value> &y)
534 {
535 for (unsigned int rr=0; rr<dimR; ++rr)
536 deriv_[rr].assign(y,rr);
537 }
538 // assign with scalar functions to component r
539 template <class Fy,DerivativeLayoutNS::DerivativeLayout layouty>
540 void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,layouty> &y)
541 {
542 deriv_[r].assign(r,y);
543 }
544
545 Block &block()
546 {
547 return reinterpret_cast<Block&>(*this);
548 }
549 const Block &block() const
550 {
551 return reinterpret_cast<const Block&>(*this);
552 }
553
554 ScalarDeriv &operator[](int r) {
555 return deriv_[r];
556 }
557 const ScalarDeriv &operator[](int r) const {
558 return deriv_[r];
559 }
560 protected:
561 Dune::FieldVector<ScalarDeriv,dimR> deriv_;
562 };
563
564 // ******************************************
565 // AXPY *************************************
566 // ******************************************
567 template <class Vec1,class Vec2,unsigned int deriv>
568 struct LFETensorAxpy
569 {
570 template <class Field>
571 static void apply(unsigned int r,const Field &a,
572 const Vec1 &x, Vec2 &y)
573 {
574 y.axpy(a,x);
575 }
576 };
577 template <class F1,int dimD,int dimR,
578 unsigned int d,
579 class Vec2,
580 unsigned int deriv>
581 struct LFETensorAxpy<Derivatives<F1,dimD,dimR,d,DerivativeLayoutNS::value>,Vec2,deriv>
582 {
583 typedef Derivatives<F1,dimD,dimR,d,DerivativeLayoutNS::value> Vec1;
584 template <class Field>
585 static void apply(unsigned int r,const Field &a,
586 const Vec1 &x, Vec2 &y)
587 {
588 const FieldVector<F1,Vec2::size> &xx = x.template block<deriv>();
589 for (int i=0; i<y.size; ++i)
590 y[i] += xx[i]*a;
591 }
592 };
593 template <class F1,int dimD,int dimR,
594 unsigned int d,
595 class Vec2,
596 unsigned int deriv>
597 struct LFETensorAxpy<Derivatives<F1,dimD,dimR,d,DerivativeLayoutNS::derivative>,Vec2,deriv>
598 {
599 typedef Derivatives<F1,dimD,dimR,d,DerivativeLayoutNS::derivative> Vec1;
600 template <class Field>
601 static void apply(unsigned int r,const Field &a,
602 const Vec1 &x, Vec2 &y)
603 {
604 for (int rr=0; rr<dimR; ++rr)
605 LFETensorAxpy<Derivatives<F1,dimD,1,d,DerivativeLayoutNS::value>,
606 Vec2,deriv>::apply(rr,a,x[rr],y);
607 }
608 };
609 template <class F1,int dimD,
610 unsigned int d,
611 class Vec2,
612 unsigned int deriv>
613 struct LFETensorAxpy<Derivatives<F1,dimD,1,d,DerivativeLayoutNS::derivative>,Vec2,deriv>
614 {
615 typedef Derivatives<F1,dimD,1,d,DerivativeLayoutNS::derivative> Vec1;
616 template <class Field>
617 static void apply(unsigned int r,const Field &a,
618 const Vec1 &x, Vec2 &y)
619 {
620 LFETensorAxpy<Derivatives<F1,dimD,1,d,DerivativeLayoutNS::value>,
621 Vec2,deriv>::apply(r,a,x[0],y);
622 }
623 };
624 template <class F1,int dimD,
625 unsigned int d,
626 class Vec2,
627 unsigned int deriv>
628 struct LFETensorAxpy<Derivatives<F1,dimD,1,d,DerivativeLayoutNS::value>,Vec2,deriv>
629 {
630 typedef Derivatives<F1,dimD,1,d,DerivativeLayoutNS::value> Vec1;
631 template <class Field>
632 static void apply(unsigned int r,const Field &a,
633 const Vec1 &x, Vec2 &y)
634 {
635 typedef LFETensor<F1,dimD,deriv> LFETensorType;
636 const unsigned int rr = r*LFETensorType::size;
637 const FieldVector<F1,LFETensorType::size> &xx = x.template block<deriv>();
638 for (int i=0; i<FieldVector<F1,LFETensorType::size>::dimension; ++i)
639 y[rr+i] += xx[i]*a;
640 }
641 };
642
643 // ***********************************************
644 // Assign ****************************************
645 // ***********************************************
646 template <class Vec1,class Vec2>
647 struct DerivativeAssign
648 {
649 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
650 {
651 field_cast(vec1,vec2);
652 }
653 };
654 template <int dimD,int dimR,unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout,
655 class F1,class F2>
656 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,layout>,
657 Derivatives<F2,dimD,dimR,deriv,layout> >
658 {
659 typedef Derivatives<F1,dimD,dimR,deriv,layout> Vec1;
660 typedef Derivatives<F2,dimD,dimR,deriv,layout> Vec2;
661 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
662 {
663 field_cast(vec1.block(),vec2.block());
664 }
665 };
666 template <int dimD,int dimR,unsigned int deriv,
667 class F1, class F2>
668 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value>,
669 Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> >
670 {
671 typedef Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value> Vec1;
672 typedef Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> Vec2;
673 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
674 {
675 vec2.assign(vec1);
676 }
677 };
678 template <int dimD,int dimR,unsigned int deriv,
679 class F1, class F2>
680 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative>,
681 Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> >
682 {
683 typedef Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative> Vec1;
684 typedef Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> Vec2;
685 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
686 {
687 vec2.assign(vec1);
688 }
689 };
690 template <int dimD,int dimR,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout,
691 class F1, class F2>
692 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
693 Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> >
694 {
695 typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
696 typedef Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> Vec2;
697 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
698 {
699 vec2.assign(r,vec1);
700 }
701 };
702 template <int dimD,int dimR,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout,
703 class F1, class F2>
704 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
705 Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> >
706 {
707 typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
708 typedef Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> Vec2;
709 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
710 {
711 vec2.assign(r,vec1);
712 }
713 };
714 template <int dimD,unsigned int deriv,
715 class F1, class F2>
716 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,
717 Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> >
718 {
719 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value> Vec1;
720 typedef Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> Vec2;
721 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
722 {
723 field_cast(vec1.block(),vec2.block());
724 }
725 };
726 template <int dimD,unsigned int deriv,
727 class F1, class F2>
728 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,
729 Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> >
730 {
731 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec1;
732 typedef Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec2;
733 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
734 {
735 field_cast(vec1.block(),vec2.block());
736 }
737 };
738 template <int dimD,unsigned int deriv,
739 class F1, class F2>
740 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,
741 Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> >
742 {
743 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec1;
744 typedef Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> Vec2;
745 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
746 {
747 field_cast(vec1.block(),vec2.block());
748 }
749 };
750 template <int dimD,unsigned int deriv,
751 class F1, class F2>
752 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,
753 Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> >
754 {
755 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value> Vec1;
756 typedef Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec2;
757 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
758 {
759 field_cast(vec1.block(),vec2.block());
760 }
761 };
762 template <int dimD,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout,
763 class F1, class F2>
764 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,layout>,
765 F2 >
766 {
767 typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
768 typedef F2 Vec2;
769 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
770 {
771 field_cast(vec1.block(),vec2);
772 }
773 };
774 template <int dimD,int dimR,
775 class F1,unsigned int deriv,
776 class F2>
777 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value>,FieldVector<F2,dimR> >
778 {
779 typedef Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value> Vec1;
780 typedef FieldVector<F2,dimR> Vec2;
781 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
782 {
783 field_cast(vec1.template block<0>(),vec2);
784 }
785 };
786 template <int dimD,int dimR,
787 class F1,unsigned int deriv,
788 class F2>
789 struct DerivativeAssign<Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative>,FieldVector<F2,dimR> >
790 {
791 typedef Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative> Vec1;
792 typedef FieldVector<F2,dimR> Vec2;
793 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
794 {
795 for (int rr=0; rr<dimR; ++rr)
796 field_cast(vec1[rr].template tensor<0>()[0].block(),vec2[rr]);
797 }
798 };
799 template <int dimD,
800 class F1,unsigned int deriv,
801 class F2,int dimR>
802 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,FieldVector<F2,dimR> >
803 {
804 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value> Vec1;
805 typedef FieldVector<F2,dimR> Vec2;
806 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
807 {
808 field_cast(vec1.template tensor<0>()[0].block(),vec2[r]);
809 }
810 };
811 template <int dimD,
812 class F1,unsigned int deriv,
813 class F2,int dimR>
814 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,FieldVector<F2,dimR> >
815 {
816 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec1;
817 typedef FieldVector<F2,dimR> Vec2;
818 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
819 {
820 field_cast(vec1[0].template tensor<0>()[0].block(),vec2[r]);
821 }
822 };
823 template <int dimD,
824 class F1,unsigned int deriv,
825 class F2>
826 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value>,FieldVector<F2,1> >
827 {
828 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value> Vec1;
829 typedef FieldVector<F2,1> Vec2;
830 static void apply(unsigned int r,const Vec1 &vec1,Vec2 &vec2)
831 {
832 field_cast(vec1.template tensor<0>()[0].block(),vec2);
833 }
834 };
835 template <int dimD,
836 class F1,unsigned int deriv,
837 class F2>
838 struct DerivativeAssign<Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative>,FieldVector<F2,1> >
839 {
840 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec1;
841 typedef FieldVector<F2,1> Vec2;
842 static void apply(unsigned int /*r*/,const Vec1 &vec1,Vec2 &vec2)
843 {
844 field_cast(vec1[0].template tensor<0>()[0].block(),vec2);
845 }
846 };
847
848 // ***********************************************
849 // IO ********************************************
850 // ***********************************************
851 template <class F,int dimD,unsigned int deriv>
852 std::ostream &operator<< ( std::ostream &out, const LFETensor< F,dimD,deriv > &tensor )
853 {
854 return out << tensor.block();
855 }
856#if 0
857 template <class F,int dimD,unsigned int deriv>
858 std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,deriv > &d )
859 {
860 out << static_cast<const ScalarDerivatives< F,dimD,deriv-1 > &>(d);
861 out << " , " << d.tensor() << std::endl;
862 return out;
863 }
864 template <class F,int dimD>
865 std::ostream &operator<< ( std::ostream &out, const ScalarDerivatives< F,dimD,0 > &d )
866 {
867 out << d.tensor() << std::endl;
868 return out;
869 }
870#endif
871 template <class F,int dimD,int dimR,unsigned int deriv>
872 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,DerivativeLayoutNS::derivative > &d )
873 {
874 out << " ( ";
875 out << d[0];
876 for (int r=1; r<dimR; ++r)
877 {
878 out << " , " << d[r];
879 }
880 out << " ) " << std::endl;
881 return out;
882 }
883 template <class F,int dimD,int dimR,unsigned int deriv>
884 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,DerivativeLayoutNS::value > &d )
885 {
886 out << static_cast<const Derivatives< F,dimD,dimR,deriv-1,DerivativeLayoutNS::value > &>(d);
887 out << " ( ";
888 out << d[0];
889 for (int r=1; r<dimR; ++r)
890 {
891 out << " , " << d[r];
892 }
893 out << " ) " << std::endl;
894 return out;
895 }
896 template <class F,int dimD,int dimR>
897 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,DerivativeLayoutNS::derivative > &d )
898 {
899 out << " ( ";
900 out << d[0];
901 for (int r=1; r<dimR; ++r)
902 {
903 out << " , " << d[r];
904 }
905 out << " ) " << std::endl;
906 return out;
907 }
908 template <class F,int dimD,int dimR>
909 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,DerivativeLayoutNS::value > &d )
910 {
911 out << " ( ";
912 out << d[0];
913 for (int r=1; r<dimR; ++r)
914 {
915 out << " , " << d[r];
916 }
917 out << " ) " << std::endl;
918 return out;
919 }
920 template <class F,int dimD,int dimR,unsigned int deriv,DerivativeLayoutNS::DerivativeLayout layout>
921 std::ostream &operator<< ( std::ostream &out, const std::vector<Derivatives< F,dimD,dimR,deriv,layout > > &y )
922 {
923 out << "Number of basis functions: " << y.size() << std::endl;
924 for (unsigned int i=0; i<y.size(); ++i)
925 {
926 out << "Base " << i << " : " << std::endl;
927 out << y[i];
928 out << std::endl;
929 }
930 return out;
931 }
932}
933#endif // DUNE_TENSOR_HH
Dune::DenseVector::axpy
derived_type & axpy(const field_type &a, const DenseVector< Other > &x)
vector space axpy operation ( *this += a x )
Definition: densevector.hh:575
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune::Functions::derivative
TrigonometricFunction< K, -cosFactor, sinFactor > derivative(const TrigonometricFunction< K, sinFactor, cosFactor > &f)
Obtain derivative of TrigonometricFunction function.
Definition: trigonometricfunction.hh:39
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::field_cast
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:159
Dune::assign
void assign(T &dst, const T &src, bool mask)
masked Simd assignment (scalar version)
Definition: simd.hh:447
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 21, 23:30, 2024)