Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines function for taking the determinant of a rank-2 tensor
6 :
7 : #pragma once
8 :
9 : #include "DataStructures/Tensor/Tensor.hpp"
10 : #include "Utilities/Gsl.hpp"
11 :
12 : namespace detail {
13 : template <typename Symm, typename Index, typename = std::nullptr_t>
14 : struct DeterminantImpl;
15 :
16 : template <typename Symm, typename Index>
17 : struct DeterminantImpl<Symm, Index, Requires<Index::dim == 1>> {
18 : template <typename T>
19 : static typename T::type apply(const T& tensor) {
20 : return get<0, 0>(tensor);
21 : }
22 : };
23 :
24 : template <typename Symm, typename Index>
25 : struct DeterminantImpl<Symm, Index, Requires<Index::dim == 2>> {
26 : template <typename T>
27 : static typename T::type apply(const T& tensor) {
28 : const auto& t00 = get<0, 0>(tensor);
29 : const auto& t01 = get<0, 1>(tensor);
30 : const auto& t10 = get<1, 0>(tensor);
31 : const auto& t11 = get<1, 1>(tensor);
32 : return t00 * t11 - t01 * t10;
33 : }
34 : };
35 :
36 : template <typename Index>
37 : struct DeterminantImpl<Symmetry<2, 1>, Index, Requires<Index::dim == 3>> {
38 : template <typename T>
39 : static typename T::type apply(const T& tensor) {
40 : const auto& t00 = get<0, 0>(tensor);
41 : const auto& t01 = get<0, 1>(tensor);
42 : const auto& t02 = get<0, 2>(tensor);
43 : const auto& t10 = get<1, 0>(tensor);
44 : const auto& t11 = get<1, 1>(tensor);
45 : const auto& t12 = get<1, 2>(tensor);
46 : const auto& t20 = get<2, 0>(tensor);
47 : const auto& t21 = get<2, 1>(tensor);
48 : const auto& t22 = get<2, 2>(tensor);
49 : return t00 * (t11 * t22 - t12 * t21) - t01 * (t10 * t22 - t12 * t20) +
50 : t02 * (t10 * t21 - t11 * t20);
51 : }
52 : };
53 :
54 : template <typename Index>
55 : struct DeterminantImpl<Symmetry<1, 1>, Index, Requires<Index::dim == 3>> {
56 : template <typename T>
57 : static typename T::type apply(const T& tensor) {
58 : const auto& t00 = get<0, 0>(tensor);
59 : const auto& t01 = get<0, 1>(tensor);
60 : const auto& t02 = get<0, 2>(tensor);
61 : const auto& t11 = get<1, 1>(tensor);
62 : const auto& t12 = get<1, 2>(tensor);
63 : const auto& t22 = get<2, 2>(tensor);
64 : return t00 * (t11 * t22 - t12 * t12) - t01 * (t01 * t22 - t12 * t02) +
65 : t02 * (t01 * t12 - t11 * t02);
66 : }
67 : };
68 :
69 : template <typename Index>
70 : struct DeterminantImpl<Symmetry<2, 1>, Index, Requires<Index::dim == 4>> {
71 : template <typename T>
72 : static typename T::type apply(const T& tensor) {
73 : const auto& t00 = get<0, 0>(tensor);
74 : const auto& t01 = get<0, 1>(tensor);
75 : const auto& t02 = get<0, 2>(tensor);
76 : const auto& t03 = get<0, 3>(tensor);
77 : const auto& t10 = get<1, 0>(tensor);
78 : const auto& t11 = get<1, 1>(tensor);
79 : const auto& t12 = get<1, 2>(tensor);
80 : const auto& t13 = get<1, 3>(tensor);
81 : const auto& t20 = get<2, 0>(tensor);
82 : const auto& t21 = get<2, 1>(tensor);
83 : const auto& t22 = get<2, 2>(tensor);
84 : const auto& t23 = get<2, 3>(tensor);
85 : const auto& t30 = get<3, 0>(tensor);
86 : const auto& t31 = get<3, 1>(tensor);
87 : const auto& t32 = get<3, 2>(tensor);
88 : const auto& t33 = get<3, 3>(tensor);
89 : const auto minor1 = t22 * t33 - t23 * t32;
90 : const auto minor2 = t21 * t33 - t23 * t31;
91 : const auto minor3 = t20 * t33 - t23 * t30;
92 : const auto minor4 = t21 * t32 - t22 * t31;
93 : const auto minor5 = t20 * t32 - t22 * t30;
94 : const auto minor6 = t20 * t31 - t21 * t30;
95 : return t00 * (t11 * minor1 - t12 * minor2 + t13 * minor4) -
96 : t01 * (t10 * minor1 - t12 * minor3 + t13 * minor5) +
97 : t02 * (t10 * minor2 - t11 * minor3 + t13 * minor6) -
98 : t03 * (t10 * minor4 - t11 * minor5 + t12 * minor6);
99 : }
100 : };
101 :
102 : template <typename Index>
103 : struct DeterminantImpl<Symmetry<1, 1>, Index, Requires<Index::dim == 4>> {
104 : template <typename T>
105 : static typename T::type apply(const T& tensor) {
106 : const auto& t00 = get<0, 0>(tensor);
107 : const auto& t01 = get<0, 1>(tensor);
108 : const auto& t02 = get<0, 2>(tensor);
109 : const auto& t03 = get<0, 3>(tensor);
110 : const auto& t11 = get<1, 1>(tensor);
111 : const auto& t12 = get<1, 2>(tensor);
112 : const auto& t13 = get<1, 3>(tensor);
113 : const auto& t22 = get<2, 2>(tensor);
114 : const auto& t23 = get<2, 3>(tensor);
115 : const auto& t33 = get<3, 3>(tensor);
116 : const auto minor1 = t22 * t33 - t23 * t23;
117 : const auto minor2 = t12 * t33 - t23 * t13;
118 : const auto minor3 = t02 * t33 - t23 * t03;
119 : const auto minor4 = t12 * t23 - t22 * t13;
120 : const auto minor5 = t02 * t23 - t22 * t03;
121 : const auto minor6 = t02 * t13 - t12 * t03;
122 : return t00 * (t11 * minor1 - t12 * minor2 + t13 * minor4) -
123 : t01 * (t01 * minor1 - t12 * minor3 + t13 * minor5) +
124 : t02 * (t01 * minor2 - t11 * minor3 + t13 * minor6) -
125 : t03 * (t01 * minor4 - t11 * minor5 + t12 * minor6);
126 : }
127 : };
128 : } // namespace detail
129 :
130 : /// @{
131 : /*!
132 : * \ingroup TensorGroup
133 : * \brief Computes the determinant of a rank-2 Tensor `tensor`.
134 : *
135 : * \requires That `tensor` be a rank-2 Tensor, with both indices sharing the
136 : * same dimension and type.
137 : */
138 : template <typename T, typename Symm, typename Index0, typename Index1>
139 1 : void determinant(const gsl::not_null<Scalar<T>*> det_tensor,
140 : const Tensor<T, Symm, index_list<Index0, Index1>>& tensor) {
141 : static_assert(Index0::dim == Index1::dim,
142 : "Cannot take the determinant of a Tensor whose Indices are not "
143 : "of the same dimensionality.");
144 : static_assert(Index0::index_type == Index1::index_type,
145 : "Taking the determinant of a mixed Spatial and Spacetime index "
146 : "Tensor is not allowed since it's not clear what that means.");
147 : get(*det_tensor) = detail::DeterminantImpl<Symm, Index0>::apply(tensor);
148 : }
149 :
150 : template <typename T, typename Symm, typename Index0, typename Index1>
151 1 : Scalar<T> determinant(
152 : const Tensor<T, Symm, index_list<Index0, Index1>>& tensor) {
153 : Scalar<T> result{};
154 : determinant(make_not_null(&result), tensor);
155 : return result;
156 : }
157 : /// @}
|
