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             : /// @}