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