deal.II version GIT relicensing-2848-g5241f990fb 2025-03-16 19:30:00+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
23
24#include <cstddef>
25
26
28
72{
73 // forward declarations
74 namespace internal
75 {
76 template <int index, int rank, typename T>
78 template <int position, int rank>
79 struct ExtractHelper;
80 template <int no_contr, int rank_1, int rank_2, int dim>
81 class Contract;
82 template <int rank_1, int rank_2, int dim>
83 class Contract3;
84 } // namespace internal
85
86
103 template <typename T>
105 {
106 using value_type = typename T::value_type;
107 };
108
109 template <typename T>
110 struct ValueType<const T>
111 {
112 using value_type = const typename T::value_type;
113 };
114
115 template <typename T, std::size_t N>
116 struct ValueType<T[N]>
117 {
118 using value_type = T;
119 };
120
121 template <typename T, std::size_t N>
122 struct ValueType<const T[N]>
123 {
124 using value_type = const T;
125 };
126
127
135 template <int deref_steps, typename T>
137 {
139 typename ReturnType<deref_steps - 1,
141 };
142
143 template <typename T>
144 struct ReturnType<0, T>
145 {
146 using value_type = T;
147 };
148
149
187 template <int index, int rank, typename T>
190 {
191 static_assert(0 <= index && index < rank,
192 "The specified index must lie within the range [0,rank)");
193
195 }
196
197
219 template <int rank, typename T, typename ArrayType>
221 extract(T &t, const ArrayType &indices)
222 {
223 return internal::ExtractHelper<0, rank>::template extract<T, ArrayType>(
224 t, indices);
225 }
226
227
266 template <int no_contr,
267 int rank_1,
268 int rank_2,
269 int dim,
270 typename T1,
271 typename T2,
272 typename T3>
273 constexpr inline DEAL_II_ALWAYS_INLINE void
274 contract(T1 &result, const T2 &left, const T3 &right)
275 {
276 static_assert(rank_1 >= no_contr,
277 "The rank of the left tensor must be "
278 "equal or greater than the number of "
279 "contractions");
280 static_assert(rank_2 >= no_contr,
281 "The rank of the right tensor must be "
282 "equal or greater than the number of "
283 "contractions");
284
286 template contract<T1, T2, T3>(result, left, right);
287 }
288
289
318 template <int rank_1,
319 int rank_2,
320 int dim,
321 typename T1,
322 typename T2,
323 typename T3,
324 typename T4>
325 constexpr T1
326 contract3(const T2 &left, const T3 &middle, const T4 &right)
327 {
329 template contract3<T1, T2, T3, T4>(left, middle, right);
330 }
331
332
333 namespace internal
334 {
335 // -------------------------------------------------------------------------
336 // Forward declarations and type traits
337 // -------------------------------------------------------------------------
338
339 template <int rank, typename S>
340 class StoreIndex;
341 template <typename T>
342 class Identity;
343 template <int no_contr, int dim>
344 class Contract2;
345
355 template <typename T>
357 {
358 using type = T &;
359 };
360
361 template <int rank, typename S>
362 struct ReferenceType<StoreIndex<rank, S>>
363 {
365 };
366
367 template <int index, int rank, typename T>
372
373
374 // TODO: Is there a possibility to just have the following block of
375 // explanation on an internal page in doxygen? If, yes. Doxygen
376 // wizards, your call!
377
378 // -------------------------------------------------------------------------
379 // Implementation of helper classes for reordered_index_view
380 // -------------------------------------------------------------------------
381
382 // OK. This is utterly brutal template magic. Therefore, we will not
383 // comment on the individual internal helper classes, because this is
384 // of not much value, but explain the general recursion procedure.
385 //
386 // (In order of appearance)
387 //
388 // Our task is to reorder access to a tensor object where a specified
389 // index is moved to the end. Thus we want to construct an object
390 // <code>reordered</code> out of a <code>tensor</code> where the
391 // following access patterns are equivalent:
392 // @code
393 // tensor [i_0]...[i_index-1][i_index][i_index+1]...[i_n]
394 // reordered [i_0]...[i_index_1][i_index+1]...[i_n][i_index]
395 // @endcode
396 //
397 // The first task is to get rid of the application of
398 // [i_0]...[i_index-1]. This is a classical recursion pattern - relay
399 // the task from <index, rank> to <index-1, rank-1> by accessing the
400 // subtensor object:
401
402 template <int index, int rank, typename T>
404 {
405 public:
407 : t_(t)
408 {}
409
411 rank - 1,
412 typename ValueType<T>::value_type>;
413
414 // Recurse by applying index j directly:
416 operator[](unsigned int j) const
417 {
418 return value_type(t_[j]);
419 }
420
421 private:
423 };
424
425 // At some point we hit the condition index == 0 and rank > 1, i.e.,
426 // the first index should be reordered to the end.
427 //
428 // At this point we cannot be lazy any more and have to start storing
429 // indices because we get them in the wrong order. The user supplies
430 // [i_0][i_1]...[i_{rank - 1}]
431 // but we have to call the subtensor object with
432 // [i_{rank - 1}[i_0][i_1]...[i_{rank-2}]
433 //
434 // So give up and relay the task to the StoreIndex class:
435
436 template <int rank, typename T>
437 class ReorderedIndexView<0, rank, T>
438 {
439 public:
441 : t_(t)
442 {}
443
445
447 operator[](unsigned int j) const
448 {
449 return value_type(Identity<T>(t_), j);
450 }
451
452 private:
454 };
455
456 // Sometimes, we're lucky and don't have to do anything. In this case
457 // just return the original tensor.
458
459 template <typename T>
460 class ReorderedIndexView<0, 1, T>
461 {
462 public:
464 : t_(t)
465 {}
466
469
471 operator[](unsigned int j) const
472 {
473 return t_[j];
474 }
475
476 private:
478 };
479
480 // Here, Identity is a helper class to ground the recursion in
481 // StoreIndex. Its implementation is easy - we haven't stored any
482 // indices yet. So, we just provide a function apply that returns the
483 // application of an index j to the stored tensor t_:
484
485 template <typename T>
487 {
488 public:
489 constexpr Identity(typename ReferenceType<T>::type t)
490 : t_(t)
491 {}
492
494
496 apply(unsigned int j) const
497 {
498 return t_[j];
499 }
500
501 private:
503 };
504
505 // StoreIndex is a class that stores an index recursively with every
506 // invocation of operator[](unsigned int j): We do this by recursively
507 // creating a new StoreIndex class of lower rank that stores the
508 // supplied index j and holds a copy of the current class (with all
509 // other stored indices). Again, we provide an apply member function
510 // that knows how to apply an index on the highest rank and all
511 // subsequently stored indices:
512
513 template <int rank, typename S>
515 {
516 public:
517 constexpr StoreIndex(S s, int i)
518 : s_(s)
519 , i_(i)
520 {}
521
523
525 operator[](unsigned int j) const
526 {
527 return value_type(*this, j);
528 }
529
532
533 constexpr typename ReferenceType<return_type>::type
534 apply(unsigned int j) const
535 {
536 return s_.apply(j)[i_];
537 }
538
539 private:
540 const S s_;
541 const int i_;
542 };
543
544 // We have to store indices until we hit rank == 1. Then, upon the next
545 // invocation of operator[](unsigned int j) we have all necessary
546 // information available to return the actual object.
547
548 template <typename S>
549 class StoreIndex<1, S>
550 {
551 public:
552 constexpr StoreIndex(S s, int i)
553 : s_(s)
554 , i_(i)
555 {}
556
560
562 operator[](unsigned int j) const
563 {
564 return s_.apply(j)[i_];
565 }
566
567 private:
568 const S s_;
569 const int i_;
570 };
571
572
573 // -------------------------------------------------------------------------
574 // Implementation of helper classes for extract
575 // -------------------------------------------------------------------------
576
577 // Straightforward recursion implemented by specializing ExtractHelper
578 // for position == rank. We use the type trait ReturnType<rank, T> to
579 // have an idea what the final type will be.
580 template <int position, int rank>
582 {
583 template <typename T, typename ArrayType>
584 constexpr static typename ReturnType<rank - position, T>::value_type &
585 extract(T &t, const ArrayType &indices)
586 {
589 ArrayType>(t[indices[position]], indices);
590 }
591 };
592
593 // For position == rank there is nothing to extract, just return the
594 // object.
595 template <int rank>
596 struct ExtractHelper<rank, rank>
597 {
598 template <typename T, typename ArrayType>
599 constexpr static T &
600 extract(T &t, const ArrayType &)
601 {
602 return t;
603 }
604 };
605
606
607 // -------------------------------------------------------------------------
608 // Implementation of helper classes for contract
609 // -------------------------------------------------------------------------
610
611 // Straightforward recursive pattern:
612 //
613 // As long as rank_1 > no_contr, assign indices from the left tensor to
614 // result. This builds up the first part of the nested outer loops:
615 //
616 // for(unsigned int i_0; i_0 < dim; ++i_0)
617 // ...
618 // for(i_; i_ < dim; ++i_)
619 // [...]
620 // result[i_0]..[i_] ... left[i_0]..[i_] ...
621
622 template <int no_contr, int rank_1, int rank_2, int dim>
624 {
625 public:
626 template <typename T1, typename T2, typename T3>
627 constexpr inline DEAL_II_ALWAYS_INLINE static void
628 contract(T1 &result, const T2 &left, const T3 &right)
629 {
630 for (unsigned int i = 0; i < dim; ++i)
632 left[i],
633 right);
634 }
635 };
636
637 // If rank_1 == no_contr leave out the remaining no_contr indices for
638 // the contraction and assign indices from the right tensor to the
639 // result. This builds up the second part of the nested loops:
640 //
641 // for(unsigned int i_0 = 0; i_0 < dim; ++i_0)
642 // ...
643 // for(unsigned int i_ = 0; i_ < dim; ++i_)
644 // for(unsigned int j_0 = 0; j_0 < dim; ++j_0)
645 // ...
646 // for(unsigned int j_ = 0; j_ < dim; ++j_)
647 // [...]
648 // result[i_0]..[i_][j_0]..[j_] ... left[i_0]..[i_] ...
649 // right[j_0]..[j_]
650 //
651
652 template <int no_contr, int rank_2, int dim>
653 class Contract<no_contr, no_contr, rank_2, dim>
654 {
655 public:
656 template <typename T1, typename T2, typename T3>
657 constexpr inline DEAL_II_ALWAYS_INLINE static void
658 contract(T1 &result, const T2 &left, const T3 &right)
659 {
660 for (unsigned int i = 0; i < dim; ++i)
662 left,
663 right[i]);
664 }
665 };
666
667 // If rank_1 == rank_2 == no_contr we have built up all of the outer
668 // loop. Now, it is time to do the actual contraction:
669 //
670 // [...]
671 // {
672 // result[i_0]..[i_][j_0]..[j_] = 0.;
673 // for(unsigned int k_0 = 0; k_0 < dim; ++k_0)
674 // ...
675 // for(unsigned int k_ = 0; k_ < dim; ++k_)
676 // result[i_0]..[i_][j_0]..[j_] += left[i_0]..[i_][k_0]..[k_] *
677 // right[j_0]..[j_][k_0]..[k_];
678 // }
679 //
680 // Relay this summation to another helper class.
681
682 template <int no_contr, int dim>
683 class Contract<no_contr, no_contr, no_contr, dim>
684 {
685 public:
686 template <typename T1, typename T2, typename T3>
687 constexpr inline DEAL_II_ALWAYS_INLINE static void
688 contract(T1 &result, const T2 &left, const T3 &right)
689 {
690 result = Contract2<no_contr, dim>::template contract2<T1>(left, right);
691 }
692 };
693
694 // Straightforward recursion:
695 //
696 // Contract leftmost index and recurse one down.
697
698 template <int no_contr, int dim>
700 {
701 public:
702 template <typename T1, typename T2, typename T3>
703 constexpr inline DEAL_II_ALWAYS_INLINE static T1
704 contract2(const T2 &left, const T3 &right)
705 {
706 // Some auto-differentiable numbers need explicit
707 // zero initialization.
708 if (dim == 0)
709 {
710 T1 result = ::internal::NumberType<T1>::value(0.0);
711 return result;
712 }
713 else
714 {
715 T1 result =
716 Contract2<no_contr - 1, dim>::template contract2<T1>(left[0],
717 right[0]);
718 for (unsigned int i = 1; i < dim; ++i)
719 result +=
720 Contract2<no_contr - 1, dim>::template contract2<T1>(left[i],
721 right[i]);
722 return result;
723 }
724 }
725 };
726
727 // A contraction of two objects of order 0 is just a scalar
728 // multiplication:
729
730 template <int dim>
731 class Contract2<0, dim>
732 {
733 public:
734 template <typename T1, typename T2, typename T3>
735 constexpr DEAL_II_ALWAYS_INLINE static T1
736 contract2(const T2 &left, const T3 &right)
737 {
738 return left * right;
739 }
740 };
741
742
743 // -------------------------------------------------------------------------
744 // Implementation of helper classes for contract3
745 // -------------------------------------------------------------------------
746
747 // Fully contract three tensorial objects
748 //
749 // As long as rank_1 > 0, recurse over left and middle:
750 //
751 // for(unsigned int i_0; i_0 < dim; ++i_0)
752 // ...
753 // for(i_; i_ < dim; ++i_)
754 // [...]
755 // left[i_0]..[i_] ... middle[i_0]..[i_] ... right
756
757 template <int rank_1, int rank_2, int dim>
759 {
760 public:
761 template <typename T1, typename T2, typename T3, typename T4>
762 constexpr static inline T1
763 contract3(const T2 &left, const T3 &middle, const T4 &right)
764 {
765 // Some auto-differentiable numbers need explicit
766 // zero initialization.
767 T1 result = ::internal::NumberType<T1>::value(0.0);
768 for (unsigned int i = 0; i < dim; ++i)
769 result += Contract3<rank_1 - 1, rank_2, dim>::template contract3<T1>(
770 left[i], middle[i], right);
771 return result;
772 }
773 };
774
775 // If rank_1 ==0, continue to recurse over middle and right:
776 //
777 // for(unsigned int i_0; i_0 < dim; ++i_0)
778 // ...
779 // for(i_; i_ < dim; ++i_)
780 // for(unsigned int j_0; j_0 < dim; ++j_0)
781 // ...
782 // for(j_; j_ < dim; ++j_)
783 // [...]
784 // left[i_0]..[i_] ... middle[i_0]..[i_][j_0]..[j_] ...
785 // right[j_0]..[j_]
786
787 template <int rank_2, int dim>
788 class Contract3<0, rank_2, dim>
789 {
790 public:
791 template <typename T1, typename T2, typename T3, typename T4>
792 constexpr static inline T1
793 contract3(const T2 &left, const T3 &middle, const T4 &right)
794 {
795 // Some auto-differentiable numbers need explicit
796 // zero initialization.
797 T1 result = ::internal::NumberType<T1>::value(0.0);
798 for (unsigned int i = 0; i < dim; ++i)
799 result +=
801 middle[i],
802 right[i]);
803 return result;
804 }
805 };
806
807 // Contraction of three tensorial objects of rank 0 is just a scalar
808 // multiplication.
809
810 template <int dim>
811 class Contract3<0, 0, dim>
812 {
813 public:
814 template <typename T1, typename T2, typename T3, typename T4>
815 constexpr static T1
816 contract3(const T2 &left, const T3 &middle, const T4 &right)
817 {
818 return left * middle * right;
819 }
820 };
821
822 // -------------------------------------------------------------------------
823
824 } /* namespace internal */
825} /* namespace TensorAccessors */
826
828
829#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:164
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
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 const T & value(const T &t)
Definition numbers.h:682