Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines functions euclidean dot_product and dot_product with a metric
       6             : 
       7             : #pragma once
       8             : 
       9             : #include <cstddef>
      10             : 
      11             : #include "DataStructures/Tensor/Tensor.hpp"
      12             : #include "Utilities/ContainerHelpers.hpp"
      13             : #include "Utilities/MakeWithValue.hpp"
      14             : 
      15             : /// @{
      16             : /*!
      17             :  * \ingroup TensorGroup
      18             :  * \brief Compute the Euclidean dot product of two vectors or one forms
      19             :  *
      20             :  * \details
      21             :  * Returns \f$A^a B^b \delta_{ab}\f$ for input vectors \f$A^a\f$ and \f$B^b\f$
      22             :  * or \f$A_a B_b \delta^{ab}\f$ for input one forms \f$A_a\f$ and \f$B_b\f$.
      23             :  */
      24             : template <typename DataTypeLhs, typename DataTypeRhs, typename Index,
      25             :           typename DataTypeResult = decltype(blaze::evaluate(DataTypeLhs() *
      26             :                                                              DataTypeRhs()))>
      27           1 : void dot_product(
      28             :     const gsl::not_null<Scalar<DataTypeResult>*> dot_product,
      29             :     const Tensor<DataTypeLhs, Symmetry<1>, index_list<Index>>& vector_a,
      30             :     const Tensor<DataTypeRhs, Symmetry<1>, index_list<Index>>& vector_b) {
      31             :   get(*dot_product) = get<0>(vector_a) * get<0>(vector_b);
      32             :   for (size_t d = 1; d < Index::dim; ++d) {
      33             :     get(*dot_product) += vector_a.get(d) * vector_b.get(d);
      34             :   }
      35             : }
      36             : 
      37             : template <typename DataTypeLhs, typename DataTypeRhs, typename Index,
      38             :           typename DataTypeResult = decltype(blaze::evaluate(DataTypeLhs() *
      39             :                                                              DataTypeRhs()))>
      40           1 : Scalar<DataTypeResult> dot_product(
      41             :     const Tensor<DataTypeLhs, Symmetry<1>, index_list<Index>>& vector_a,
      42             :     const Tensor<DataTypeRhs, Symmetry<1>, index_list<Index>>& vector_b) {
      43             :   Scalar<DataTypeResult> dot_product(get_size(get<0>(vector_a)));
      44             :   ::dot_product(make_not_null(&dot_product), vector_a, vector_b);
      45             :   return dot_product;
      46             : }
      47             : /// @}
      48             : 
      49             : /// @{
      50             : /*!
      51             :  * \ingroup TensorGroup
      52             :  * \brief Compute the dot product of a vector and a one form
      53             :  *
      54             :  * \details
      55             :  * Returns \f$A^a B_b \delta_{a}^b\f$ for input vector \f$A^a\f$ and
      56             :  * input one form \f$B_b\f$
      57             :  * or \f$A_a B^b \delta^a_b\f$ for input one form \f$A_a\f$ and
      58             :  * input vector \f$B^b\f$.
      59             :  */
      60             : template <typename DataTypeLhs, typename DataTypeRhs, typename Index,
      61             :           typename DataTypeResult = decltype(blaze::evaluate(DataTypeLhs() *
      62             :                                                              DataTypeRhs()))>
      63           1 : void dot_product(
      64             :     // NOLINTNEXTLINE(readability-avoid-const-params-in-decls) false positive
      65             :     const gsl::not_null<Scalar<DataTypeResult>*> dot_product,
      66             :     const Tensor<DataTypeLhs, Symmetry<1>, index_list<Index>>& vector_a,
      67             :     const Tensor<DataTypeRhs, Symmetry<1>,
      68             :                  index_list<change_index_up_lo<Index>>>& vector_b) {
      69             :   get(*dot_product) = get<0>(vector_a) * get<0>(vector_b);
      70             :   for (size_t d = 1; d < Index::dim; ++d) {
      71             :     get(*dot_product) += vector_a.get(d) * vector_b.get(d);
      72             :   }
      73             : }
      74             : 
      75             : template <typename DataTypeLhs, typename DataTypeRhs, typename Index,
      76             :           typename DataTypeResult = decltype(blaze::evaluate(DataTypeLhs() *
      77             :                                                              DataTypeRhs()))>
      78           1 : Scalar<DataTypeResult> dot_product(
      79             :     const Tensor<DataTypeLhs, Symmetry<1>, index_list<Index>>& vector_a,
      80             :     const Tensor<DataTypeRhs, Symmetry<1>,
      81             :                  index_list<change_index_up_lo<Index>>>& vector_b) {
      82             :   Scalar<DataTypeResult> dot_product(get_size(get<0>(vector_a)));
      83             :   ::dot_product(make_not_null(&dot_product), vector_a, vector_b);
      84             :   return dot_product;
      85             : }
      86             : /// @}
      87             : 
      88             : /// @{
      89             : /*!
      90             :  * \ingroup TensorGroup
      91             :  * \brief Compute the dot_product of two vectors or one forms
      92             :  *
      93             :  * \details
      94             :  * Returns \f$g_{ab} A^a B^b\f$, where \f$g_{ab}\f$ is the metric,
      95             :  * \f$A^a\f$ is vector_a, and \f$B^b\f$ is vector_b.
      96             :  * Or, returns \f$g^{ab} A_a B_b\f$ when given one forms \f$A_a\f$
      97             :  * and \f$B_b\f$ with an inverse metric \f$g^{ab}\f$.
      98             :  */
      99             : template <typename DataTypeLhs, typename DataTypeRhs, typename DataTypeMetric,
     100             :           typename Index,
     101             :           typename DataTypeResult = decltype(blaze::evaluate(
     102             :               DataTypeLhs() * DataTypeRhs() * DataTypeMetric()))>
     103           1 : void dot_product(
     104             :     const gsl::not_null<Scalar<DataTypeResult>*> dot_product,
     105             :     const Tensor<DataTypeLhs, Symmetry<1>, index_list<Index>>& vector_a,
     106             :     const Tensor<DataTypeRhs, Symmetry<1>, index_list<Index>>& vector_b,
     107             :     const Tensor<DataTypeMetric, Symmetry<1, 1>,
     108             :                  index_list<change_index_up_lo<Index>,
     109             :                             change_index_up_lo<Index>>>& metric) {
     110             :   if constexpr (Index::dim == 1) {
     111             :     get(*dot_product) = get<0>(vector_a) * get<0>(vector_b) * get<0, 0>(metric);
     112             :   } else {
     113             :     if (&vector_a == &vector_b) {
     114             :       get(*dot_product) =
     115             :           get<0>(vector_a) * get<1>(vector_a) * get<0, 1>(metric);
     116             :       for (size_t j = 2; j < Index::dim; ++j) {
     117             :         get(*dot_product) +=
     118             :             vector_a.get(0) * vector_a.get(j) * metric.get(0, j);
     119             :       }
     120             :       for (size_t i = 1; i < Index::dim; ++i) {
     121             :         for (size_t j = i + 1; j < Index::dim; ++j) {
     122             :           get(*dot_product) +=
     123             :               vector_a.get(i) * vector_a.get(j) * metric.get(i, j);
     124             :         }
     125             :       }
     126             :       get(*dot_product) *= 2.0;
     127             : 
     128             :       for (size_t i = 0; i < Index::dim; ++i) {
     129             :         get(*dot_product) += square(vector_a.get(i)) * metric.get(i, i);
     130             :       }
     131             :     } else {
     132             :       get(*dot_product) =
     133             :           get<0>(vector_a) * get<0>(vector_b) * get<0, 0>(metric);
     134             :       for (size_t b = 1; b < Index::dim; ++b) {
     135             :         get(*dot_product) +=
     136             :             get<0>(vector_a) * vector_b.get(b) * metric.get(0, b);
     137             :       }
     138             : 
     139             :       for (size_t a = 1; a < Index::dim; ++a) {
     140             :         for (size_t b = 0; b < Index::dim; ++b) {
     141             :           get(*dot_product) +=
     142             :               vector_a.get(a) * vector_b.get(b) * metric.get(a, b);
     143             :         }
     144             :       }
     145             :     }
     146             :   }
     147             : }
     148             : 
     149             : template <typename DataTypeLhs, typename DataTypeRhs, typename DataTypeMetric,
     150             :           typename Index,
     151             :           typename DataTypeResult = decltype(blaze::evaluate(
     152             :               DataTypeLhs() * DataTypeRhs() * DataTypeMetric()))>
     153           1 : Scalar<DataTypeResult> dot_product(
     154             :     const Tensor<DataTypeLhs, Symmetry<1>, index_list<Index>>& vector_a,
     155             :     const Tensor<DataTypeRhs, Symmetry<1>, index_list<Index>>& vector_b,
     156             :     const Tensor<DataTypeMetric, Symmetry<1, 1>,
     157             :                  index_list<change_index_up_lo<Index>,
     158             :                             change_index_up_lo<Index>>>& metric) {
     159             :   Scalar<DataTypeResult> dot_product(get_size(get<0>(vector_a)));
     160             :   ::dot_product(make_not_null(&dot_product), vector_a, vector_b, metric);
     161             :   return dot_product;
     162             : }
     163             : /// @}