Reference documentation for deal.II version GIT relicensing-245-g36f19064f7 2024-03-29 07:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
tensor_accessors.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2015 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_tensor_accessors_h
16#define dealii_tensor_accessors_h
17
18#include <deal.II/base/config.h>
19
22
23
25
69{
70 // forward declarations
71 namespace internal
72 {
73 template <int index, int rank, typename T>
75 template <int position, int rank>
76 struct ExtractHelper;
77 template <int no_contr, int rank_1, int rank_2, int dim>
78 class Contract;
79 template <int rank_1, int rank_2, int dim>
80 class Contract3;
81 } // namespace internal
82
83
100 template <typename T>
102 {
103 using value_type = typename T::value_type;
104 };
105
106 template <typename T>
107 struct ValueType<const T>
108 {
109 using value_type = const typename T::value_type;
110 };
111
112 template <typename T, std::size_t N>
113 struct ValueType<T[N]>
114 {
115 using value_type = T;
116 };
117
118 template <typename T, std::size_t N>
119 struct ValueType<const T[N]>
120 {
121 using value_type = const T;
122 };
123
124
132 template <int deref_steps, typename T>
134 {
136 typename ReturnType<deref_steps - 1,
138 };
139
140 template <typename T>
141 struct ReturnType<0, T>
142 {
143 using value_type = T;
144 };
145
146
184 template <int index, int rank, typename T>
187 {
188 static_assert(0 <= index && index < rank,
189 "The specified index must lie within the range [0,rank)");
190
192 }
193
194
216 template <int rank, typename T, typename ArrayType>
218 extract(T &t, const ArrayType &indices)
219 {
220 return internal::ExtractHelper<0, rank>::template extract<T, ArrayType>(
221 t, indices);
222 }
223
224
263 template <int no_contr,
264 int rank_1,
265 int rank_2,
266 int dim,
267 typename T1,
268 typename T2,
269 typename T3>
270 constexpr inline DEAL_II_ALWAYS_INLINE void
271 contract(T1 &result, const T2 &left, const T3 &right)
272 {
273 static_assert(rank_1 >= no_contr,
274 "The rank of the left tensor must be "
275 "equal or greater than the number of "
276 "contractions");
277 static_assert(rank_2 >= no_contr,
278 "The rank of the right tensor must be "
279 "equal or greater than the number of "
280 "contractions");
281
283 template contract<T1, T2, T3>(result, left, right);
284 }
285
286
315 template <int rank_1,
316 int rank_2,
317 int dim,
318 typename T1,
319 typename T2,
320 typename T3,
321 typename T4>
322 constexpr T1
323 contract3(const T2 &left, const T3 &middle, const T4 &right)
324 {
326 template contract3<T1, T2, T3, T4>(left, middle, right);
327 }
328
329
330 namespace internal
331 {
332 // -------------------------------------------------------------------------
333 // Forward declarations and type traits
334 // -------------------------------------------------------------------------
335
336 template <int rank, typename S>
337 class StoreIndex;
338 template <typename T>
339 class Identity;
340 template <int no_contr, int dim>
341 class Contract2;
342
352 template <typename T>
354 {
355 using type = T &;
356 };
357
358 template <int rank, typename S>
359 struct ReferenceType<StoreIndex<rank, S>>
360 {
362 };
363
364 template <int index, int rank, typename T>
369
370
371 // TODO: Is there a possibility to just have the following block of
372 // explanation on an internal page in doxygen? If, yes. Doxygen
373 // wizards, your call!
374
375 // -------------------------------------------------------------------------
376 // Implementation of helper classes for reordered_index_view
377 // -------------------------------------------------------------------------
378
379 // OK. This is utterly brutal template magic. Therefore, we will not
380 // comment on the individual internal helper classes, because this is
381 // of not much value, but explain the general recursion procedure.
382 //
383 // (In order of appearance)
384 //
385 // Our task is to reorder access to a tensor object where a specified
386 // index is moved to the end. Thus we want to construct an object
387 // <code>reordered</code> out of a <code>tensor</code> where the
388 // following access patterns are equivalent:
389 // @code
390 // tensor [i_0]...[i_index-1][i_index][i_index+1]...[i_n]
391 // reordered [i_0]...[i_index_1][i_index+1]...[i_n][i_index]
392 // @endcode
393 //
394 // The first task is to get rid of the application of
395 // [i_0]...[i_index-1]. This is a classical recursion pattern - relay
396 // the task from <index, rank> to <index-1, rank-1> by accessing the
397 // subtensor object:
398
399 template <int index, int rank, typename T>
401 {
402 public:
404 : t_(t)
405 {}
406
408 rank - 1,
409 typename ValueType<T>::value_type>;
410
411 // Recurse by applying index j directly:
413 operator[](unsigned int j) const
414 {
415 return value_type(t_[j]);
416 }
417
418 private:
420 };
421
422 // At some point we hit the condition index == 0 and rank > 1, i.e.,
423 // the first index should be reordered to the end.
424 //
425 // At this point we cannot be lazy any more and have to start storing
426 // indices because we get them in the wrong order. The user supplies
427 // [i_0][i_1]...[i_{rank - 1}]
428 // but we have to call the subtensor object with
429 // [i_{rank - 1}[i_0][i_1]...[i_{rank-2}]
430 //
431 // So give up and relay the task to the StoreIndex class:
432
433 template <int rank, typename T>
434 class ReorderedIndexView<0, rank, T>
435 {
436 public:
438 : t_(t)
439 {}
440
442
444 operator[](unsigned int j) const
445 {
446 return value_type(Identity<T>(t_), j);
447 }
448
449 private:
451 };
452
453 // Sometimes, we're lucky and don't have to do anything. In this case
454 // just return the original tensor.
455
456 template <typename T>
457 class ReorderedIndexView<0, 1, T>
458 {
459 public:
461 : t_(t)
462 {}
463
466
468 operator[](unsigned int j) const
469 {
470 return t_[j];
471 }
472
473 private:
475 };
476
477 // Here, Identity is a helper class to ground the recursion in
478 // StoreIndex. Its implementation is easy - we haven't stored any
479 // indices yet. So, we just provide a function apply that returns the
480 // application of an index j to the stored tensor t_:
481
482 template <typename T>
484 {
485 public:
486 constexpr Identity(typename ReferenceType<T>::type t)
487 : t_(t)
488 {}
489
491
493 apply(unsigned int j) const
494 {
495 return t_[j];
496 }
497
498 private:
500 };
501
502 // StoreIndex is a class that stores an index recursively with every
503 // invocation of operator[](unsigned int j): We do this by recursively
504 // creating a new StoreIndex class of lower rank that stores the
505 // supplied index j and holds a copy of the current class (with all
506 // other stored indices). Again, we provide an apply member function
507 // that knows how to apply an index on the highest rank and all
508 // subsequently stored indices:
509
510 template <int rank, typename S>
512 {
513 public:
514 constexpr StoreIndex(S s, int i)
515 : s_(s)
516 , i_(i)
517 {}
518
520
522 operator[](unsigned int j) const
523 {
524 return value_type(*this, j);
525 }
526
529
530 constexpr typename ReferenceType<return_type>::type
531 apply(unsigned int j) const
532 {
533 return s_.apply(j)[i_];
534 }
535
536 private:
537 const S s_;
538 const int i_;
539 };
540
541 // We have to store indices until we hit rank == 1. Then, upon the next
542 // invocation of operator[](unsigned int j) we have all necessary
543 // information available to return the actual object.
544
545 template <typename S>
546 class StoreIndex<1, S>
547 {
548 public:
549 constexpr StoreIndex(S s, int i)
550 : s_(s)
551 , i_(i)
552 {}
553
557
559 operator[](unsigned int j) const
560 {
561 return s_.apply(j)[i_];
562 }
563
564 private:
565 const S s_;
566 const int i_;
567 };
568
569
570 // -------------------------------------------------------------------------
571 // Implementation of helper classes for extract
572 // -------------------------------------------------------------------------
573
574 // Straightforward recursion implemented by specializing ExtractHelper
575 // for position == rank. We use the type trait ReturnType<rank, T> to
576 // have an idea what the final type will be.
577 template <int position, int rank>
579 {
580 template <typename T, typename ArrayType>
581 constexpr static typename ReturnType<rank - position, T>::value_type &
582 extract(T &t, const ArrayType &indices)
583 {
586 ArrayType>(t[indices[position]], indices);
587 }
588 };
589
590 // For position == rank there is nothing to extract, just return the
591 // object.
592 template <int rank>
593 struct ExtractHelper<rank, rank>
594 {
595 template <typename T, typename ArrayType>
596 constexpr static T &
597 extract(T &t, const ArrayType &)
598 {
599 return t;
600 }
601 };
602
603
604 // -------------------------------------------------------------------------
605 // Implementation of helper classes for contract
606 // -------------------------------------------------------------------------
607
608 // Straightforward recursive pattern:
609 //
610 // As long as rank_1 > no_contr, assign indices from the left tensor to
611 // result. This builds up the first part of the nested outer loops:
612 //
613 // for(unsigned int i_0; i_0 < dim; ++i_0)
614 // ...
615 // for(i_; i_ < dim; ++i_)
616 // [...]
617 // result[i_0]..[i_] ... left[i_0]..[i_] ...
618
619 template <int no_contr, int rank_1, int rank_2, int dim>
621 {
622 public:
623 template <typename T1, typename T2, typename T3>
624 constexpr inline DEAL_II_ALWAYS_INLINE static void
625 contract(T1 &result, const T2 &left, const T3 &right)
626 {
627 for (unsigned int i = 0; i < dim; ++i)
629 left[i],
630 right);
631 }
632 };
633
634 // If rank_1 == no_contr leave out the remaining no_contr indices for
635 // the contraction and assign indices from the right tensor to the
636 // result. This builds up the second part of the nested loops:
637 //
638 // for(unsigned int i_0 = 0; i_0 < dim; ++i_0)
639 // ...
640 // for(unsigned int i_ = 0; i_ < dim; ++i_)
641 // for(unsigned int j_0 = 0; j_0 < dim; ++j_0)
642 // ...
643 // for(unsigned int j_ = 0; j_ < dim; ++j_)
644 // [...]
645 // result[i_0]..[i_][j_0]..[j_] ... left[i_0]..[i_] ...
646 // right[j_0]..[j_]
647 //
648
649 template <int no_contr, int rank_2, int dim>
650 class Contract<no_contr, no_contr, rank_2, dim>
651 {
652 public:
653 template <typename T1, typename T2, typename T3>
654 constexpr inline DEAL_II_ALWAYS_INLINE static void
655 contract(T1 &result, const T2 &left, const T3 &right)
656 {
657 for (unsigned int i = 0; i < dim; ++i)
659 left,
660 right[i]);
661 }
662 };
663
664 // If rank_1 == rank_2 == no_contr we have built up all of the outer
665 // loop. Now, it is time to do the actual contraction:
666 //
667 // [...]
668 // {
669 // result[i_0]..[i_][j_0]..[j_] = 0.;
670 // for(unsigned int k_0 = 0; k_0 < dim; ++k_0)
671 // ...
672 // for(unsigned int k_ = 0; k_ < dim; ++k_)
673 // result[i_0]..[i_][j_0]..[j_] += left[i_0]..[i_][k_0]..[k_] *
674 // right[j_0]..[j_][k_0]..[k_];
675 // }
676 //
677 // Relay this summation to another helper class.
678
679 template <int no_contr, int dim>
680 class Contract<no_contr, no_contr, no_contr, dim>
681 {
682 public:
683 template <typename T1, typename T2, typename T3>
684 constexpr inline DEAL_II_ALWAYS_INLINE static void
685 contract(T1 &result, const T2 &left, const T3 &right)
686 {
687 result = Contract2<no_contr, dim>::template contract2<T1>(left, right);
688 }
689 };
690
691 // Straightforward recursion:
692 //
693 // Contract leftmost index and recurse one down.
694
695 template <int no_contr, int dim>
697 {
698 public:
699 template <typename T1, typename T2, typename T3>
700 constexpr inline DEAL_II_ALWAYS_INLINE static T1
701 contract2(const T2 &left, const T3 &right)
702 {
703 // Some auto-differentiable numbers need explicit
704 // zero initialization.
705 if (dim == 0)
706 {
707 T1 result = ::internal::NumberType<T1>::value(0.0);
708 return result;
709 }
710 else
711 {
712 T1 result =
713 Contract2<no_contr - 1, dim>::template contract2<T1>(left[0],
714 right[0]);
715 for (unsigned int i = 1; i < dim; ++i)
716 result +=
717 Contract2<no_contr - 1, dim>::template contract2<T1>(left[i],
718 right[i]);
719 return result;
720 }
721 }
722 };
723
724 // A contraction of two objects of order 0 is just a scalar
725 // multiplication:
726
727 template <int dim>
728 class Contract2<0, dim>
729 {
730 public:
731 template <typename T1, typename T2, typename T3>
732 constexpr DEAL_II_ALWAYS_INLINE static T1
733 contract2(const T2 &left, const T3 &right)
734 {
735 return left * right;
736 }
737 };
738
739
740 // -------------------------------------------------------------------------
741 // Implementation of helper classes for contract3
742 // -------------------------------------------------------------------------
743
744 // Fully contract three tensorial objects
745 //
746 // As long as rank_1 > 0, recurse over left and middle:
747 //
748 // for(unsigned int i_0; i_0 < dim; ++i_0)
749 // ...
750 // for(i_; i_ < dim; ++i_)
751 // [...]
752 // left[i_0]..[i_] ... middle[i_0]..[i_] ... right
753
754 template <int rank_1, int rank_2, int dim>
756 {
757 public:
758 template <typename T1, typename T2, typename T3, typename T4>
759 constexpr static inline T1
760 contract3(const T2 &left, const T3 &middle, const T4 &right)
761 {
762 // Some auto-differentiable numbers need explicit
763 // zero initialization.
764 T1 result = ::internal::NumberType<T1>::value(0.0);
765 for (unsigned int i = 0; i < dim; ++i)
766 result += Contract3<rank_1 - 1, rank_2, dim>::template contract3<T1>(
767 left[i], middle[i], right);
768 return result;
769 }
770 };
771
772 // If rank_1 ==0, continue to recurse over middle and right:
773 //
774 // for(unsigned int i_0; i_0 < dim; ++i_0)
775 // ...
776 // for(i_; i_ < dim; ++i_)
777 // for(unsigned int j_0; j_0 < dim; ++j_0)
778 // ...
779 // for(j_; j_ < dim; ++j_)
780 // [...]
781 // left[i_0]..[i_] ... middle[i_0]..[i_][j_0]..[j_] ...
782 // right[j_0]..[j_]
783
784 template <int rank_2, int dim>
785 class Contract3<0, rank_2, dim>
786 {
787 public:
788 template <typename T1, typename T2, typename T3, typename T4>
789 constexpr static inline T1
790 contract3(const T2 &left, const T3 &middle, const T4 &right)
791 {
792 // Some auto-differentiable numbers need explicit
793 // zero initialization.
794 T1 result = ::internal::NumberType<T1>::value(0.0);
795 for (unsigned int i = 0; i < dim; ++i)
796 result +=
798 middle[i],
799 right[i]);
800 return result;
801 }
802 };
803
804 // Contraction of three tensorial objects of rank 0 is just a scalar
805 // multiplication.
806
807 template <int dim>
808 class Contract3<0, 0, dim>
809 {
810 public:
811 template <typename T1, typename T2, typename T3, typename T4>
812 constexpr static T1
813 contract3(const T2 &left, const T3 &middle, const T4 &right)
814 {
815 return left * middle * right;
816 }
817 };
818
819 // -------------------------------------------------------------------------
820
821 } /* namespace internal */
822} /* namespace TensorAccessors */
823
825
826#endif /* dealii_tensor_accessors_h */
static constexpr T1 contract2(const T2 &left, const T3 &right)
static constexpr T1 contract2(const T2 &left, const T3 &right)
static constexpr T1 contract3(const T2 &left, const T3 &middle, const T4 &right)
static constexpr T1 contract3(const T2 &left, const T3 &middle, const T4 &right)
static constexpr T1 contract3(const T2 &left, const T3 &middle, const T4 &right)
static constexpr void contract(T1 &result, const T2 &left, const T3 &right)
static constexpr void contract(T1 &result, const T2 &left, const T3 &right)
static constexpr void contract(T1 &result, const T2 &left, const T3 &right)
constexpr ReferenceType< return_type >::type apply(unsigned int j) const
constexpr Identity(typename ReferenceType< T >::type t)
typename ValueType< T >::value_type return_type
typename ReferenceType< typename ValueType< T >::value_type >::type value_type
constexpr value_type operator[](unsigned int j) const
constexpr ReorderedIndexView(typename ReferenceType< T >::type t)
constexpr value_type operator[](unsigned int j) const
constexpr ReorderedIndexView(typename ReferenceType< T >::type t)
ReorderedIndexView< index - 1, rank - 1, typename ValueType< T >::value_type > value_type
constexpr value_type operator[](unsigned int j) const
constexpr ReorderedIndexView(typename ReferenceType< T >::type t)
constexpr return_type & operator[](unsigned int j) const
typename ValueType< typename S::return_type >::value_type return_type
typename ValueType< typename S::return_type >::value_type return_type
constexpr value_type operator[](unsigned int j) const
StoreIndex< rank - 1, StoreIndex< rank, S > > value_type
constexpr ReferenceType< return_type >::type apply(unsigned int j) const
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
constexpr internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
constexpr void contract(T1 &result, const T2 &left, const T3 &right)
constexpr T1 contract3(const T2 &left, const T3 &middle, const T4 &right)
typename ReturnType< deref_steps - 1, typename ValueType< T >::value_type >::value_type value_type
const typename T::value_type value_type
typename T::value_type value_type
static constexpr T & extract(T &t, const ArrayType &)
static constexpr ReturnType< rank-position, T >::value_type & extract(T &t, const ArrayType &indices)
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition numbers.h:702