DUNE-FEM (unstable)

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 © 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 rewriting...
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,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 {
116 typedef LFETensor<F,dimD,0> This;
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
136 This &operator*= ( const field_type &f )
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>
157 void assign(const LFETensor<Fy,dimD,0> &y)
158 {
159 field_cast(y.block(),block());
160 }
161
162 Block &block()
163 {
164 return block_;
165 }
166 const Block &block() const
167 {
168 return block_;
169 }
170 Block block_;
171 };
172 // ***********************************************************
173 // Structure for all derivatives up to order deriv
174 // for vector valued function
175 namespace DerivativeLayoutNS {
176 enum DerivativeLayout {value,derivative};
177 }
178 template <class F,int dimD,int dimR,unsigned int deriv,
179 DerivativeLayoutNS::DerivativeLayout layout>
180 struct Derivatives;
181
182 template <class F,int dimD,int dimR,unsigned int deriv,
183 DerivativeLayoutNS::DerivativeLayout layout>
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 {
195 typedef Derivatives<F,dimD,dimR,deriv,DerivativeLayoutNS::value> This;
196 typedef Derivatives<F,dimD,dimR,deriv-1,DerivativeLayoutNS::value> Base;
197 typedef LFETensor<F,dimD,deriv> ThisLFETensor;
198
199 typedef F Field;
200 typedef F field_type;
201
202 static const DerivativeLayoutNS::DerivativeLayout layout = DerivativeLayoutNS::value;
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 }
224 This &operator=(const Block &t)
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>
243 void assign(const Derivatives<Fy,dimD,dimR,deriv,DerivativeLayoutNS::value> &y)
244 {
245 field_cast(y.block(),block());
246 }
247 // assign with different layout (same dimRange)
248 template <class Fy>
249 void assign(const Derivatives<Fy,dimD,dimR,deriv,DerivativeLayoutNS::derivative> &y)
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>
257 void assign(const Derivatives<Fy,dimD,dimRy,deriv,DerivativeLayoutNS::value> &y,unsigned int r)
258 {
259 assign<Fy,dimRy>(y.block(),r);
260 }
261 // assign with scalar functions to component r
262 template <class Fy>
263 void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,DerivativeLayoutNS::value> &y)
264 {
265 assign(r,y.block());
266 }
267 template <class Fy>
268 void assign(unsigned int r,const Derivatives<Fy,dimD,1,deriv,DerivativeLayoutNS::derivative> &y)
269 {
270 assign(r,y[0]);
271 }
272
273 Block &block()
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 }
309 ThisLFETensor &operator[](int r) {
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>
330 void assign(const Derivatives<Fy,dimD,dimR,dy,DerivativeLayoutNS::derivative> &y)
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 {
365 typedef Derivatives<F,dimD,dimR,0,DerivativeLayoutNS::value> This;
366 typedef LFETensor<F,dimD,0> ThisLFETensor;
367
368 typedef F Field;
369 typedef F field_type;
370
371 static const DerivativeLayoutNS::DerivativeLayout layout = DerivativeLayoutNS::value;
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
390 This &operator=(const Block &t)
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>
407 void assign(const Derivatives<Fy,dimD,dimR,0,DerivativeLayoutNS::value> &y)
408 {
409 field_cast(y.block(),block());
410 }
411 template <class Fy>
412 void assign(const Derivatives<Fy,dimD,dimR,0,DerivativeLayoutNS::derivative> &y)
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>
418 void assign(const Derivatives<Fy,dimD,dimRy,0,DerivativeLayoutNS::value> &y,unsigned int r)
419 {
420 assign<Fy,dimRy>(y.block(),r);
421 }
422 template <class Fy>
423 void assign(unsigned int r,const Derivatives<Fy,dimD,1,0,DerivativeLayoutNS::value> &y)
424 {
425 tensor_[r].assign(y[0]);
426 }
427 template <class Fy>
428 void assign(unsigned int r,const Derivatives<Fy,dimD,1,0,DerivativeLayoutNS::derivative> &y)
429 {
430 tensor_[r].assign(y[0][0]);
431 }
432
433 Block &block()
434 {
435 return reinterpret_cast<Block&>(*this);
436 }
437 const Block &block() const
438 {
439 return reinterpret_cast<const Block&>(*this);
440 }
441
442 ThisLFETensor &operator[](int r) {
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>
484 void assign(const Derivatives<Fy,dimD,dimR,dy,DerivativeLayoutNS::derivative> &y)
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 {
506 typedef Derivatives<F,dimD,dimR,deriv,DerivativeLayoutNS::derivative> This;
507 typedef Derivatives<F,dimD,1,deriv,DerivativeLayoutNS::value> ScalarDeriv;
508
509 typedef F Field;
510 typedef F field_type;
511
512 static const DerivativeLayoutNS::DerivativeLayout layout = DerivativeLayoutNS::value;
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 }
524 This &operator=(const Block &t)
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>
543 void assign(const Derivatives<Fy,dimD,dimR,deriv,DerivativeLayoutNS::derivative> &y)
544 {
545 field_cast(y.block(),block());
546 }
547 // assign with different layout (same dimRange)
548 template <class Fy>
549 void assign(const Derivatives<Fy,dimD,dimR,deriv,DerivativeLayoutNS::value> &y)
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
561 Block &block()
562 {
563 return reinterpret_cast<Block&>(*this);
564 }
565 const Block &block() const
566 {
567 return reinterpret_cast<const Block&>(*this);
568 }
569
570 ScalarDeriv &operator[](int r) {
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>
584 struct LFETensorAxpy
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 {
599 typedef Derivatives<F1,dimD,dimR,d,DerivativeLayoutNS::value> Vec1;
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 {
615 typedef Derivatives<F1,dimD,dimR,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 for (int rr=0; rr<dimR; ++rr)
621 LFETensorAxpy<Derivatives<F1,dimD,1,d,DerivativeLayoutNS::value>,
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 {
631 typedef Derivatives<F1,dimD,1,d,DerivativeLayoutNS::derivative> Vec1;
632 template <class Field>
633 static void apply(unsigned int r,const Field &a,
634 const Vec1 &x, Vec2 &y)
635 {
636 LFETensorAxpy<Derivatives<F1,dimD,1,d,DerivativeLayoutNS::value>,
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 {
646 typedef Derivatives<F1,dimD,1,d,DerivativeLayoutNS::value> Vec1;
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>
663 struct DerivativeAssign
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 {
675 typedef Derivatives<F1,dimD,dimR,deriv,layout> Vec1;
676 typedef Derivatives<F2,dimD,dimR,deriv,layout> Vec2;
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 {
687 typedef Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value> Vec1;
688 typedef Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> Vec2;
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 {
699 typedef Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative> Vec1;
700 typedef Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> Vec2;
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 {
711 typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
712 typedef Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::value> Vec2;
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 {
723 typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
724 typedef Derivatives<F2,dimD,dimR,deriv,DerivativeLayoutNS::derivative> Vec2;
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 {
735 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value> Vec1;
736 typedef Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> Vec2;
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 {
747 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec1;
748 typedef Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec2;
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 {
759 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec1;
760 typedef Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::value> Vec2;
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 {
771 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value> Vec1;
772 typedef Derivatives<F2,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec2;
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 {
783 typedef Derivatives<F1,dimD,1,deriv,layout> Vec1;
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 {
795 typedef Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::value> Vec1;
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 {
807 typedef Derivatives<F1,dimD,dimR,deriv,DerivativeLayoutNS::derivative> Vec1;
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 {
820 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value> Vec1;
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 {
832 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec1;
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 {
844 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::value> Vec1;
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 {
856 typedef Derivatives<F1,dimD,1,deriv,DerivativeLayoutNS::derivative> Vec1;
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>
888 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,DerivativeLayoutNS::derivative > &d )
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>
900 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,deriv,DerivativeLayoutNS::value > &d )
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>
913 std::ostream &operator<< ( std::ostream &out, const Derivatives< F,dimD,dimR,0,DerivativeLayoutNS::derivative > &d )
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
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
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::size
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
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 (Nov 12, 23:30, 2024)