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 DataType, typename Index>
25 1 : void dot_product(
26 : const gsl::not_null<Scalar<DataType>*> dot_product,
27 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
28 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_b) {
29 : get(*dot_product) = get<0>(vector_a) * get<0>(vector_b);
30 : for (size_t d = 1; d < Index::dim; ++d) {
31 : get(*dot_product) += vector_a.get(d) * vector_b.get(d);
32 : }
33 : }
34 :
35 : template <typename DataType, typename Index>
36 1 : Scalar<DataType> dot_product(
37 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
38 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_b) {
39 : Scalar<DataType> dot_product(get_size(get<0>(vector_a)));
40 : ::dot_product(make_not_null(&dot_product), vector_a, vector_b);
41 : return dot_product;
42 : }
43 : /// @}
44 :
45 : /// @{
46 : /*!
47 : * \ingroup TensorGroup
48 : * \brief Compute the dot product of a vector and a one form
49 : *
50 : * \details
51 : * Returns \f$A^a B_b \delta_{a}^b\f$ for input vector \f$A^a\f$ and
52 : * input one form \f$B_b\f$
53 : * or \f$A_a B^b \delta^a_b\f$ for input one form \f$A_a\f$ and
54 : * input vector \f$B^b\f$.
55 : */
56 : template <typename DataType, typename Index>
57 1 : void dot_product(
58 : // NOLINTNEXTLINE(readability-avoid-const-params-in-decls) false positive
59 : const gsl::not_null<Scalar<DataType>*> dot_product,
60 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
61 : const Tensor<DataType, Symmetry<1>, index_list<change_index_up_lo<Index>>>&
62 : vector_b) {
63 : get(*dot_product) = get<0>(vector_a) * get<0>(vector_b);
64 : for (size_t d = 1; d < Index::dim; ++d) {
65 : get(*dot_product) += vector_a.get(d) * vector_b.get(d);
66 : }
67 : }
68 :
69 : template <typename DataType, typename Index>
70 1 : Scalar<DataType> dot_product(
71 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
72 : const Tensor<DataType, Symmetry<1>, index_list<change_index_up_lo<Index>>>&
73 : vector_b) {
74 : Scalar<DataType> dot_product(get_size(get<0>(vector_a)));
75 : ::dot_product(make_not_null(&dot_product), vector_a, vector_b);
76 : return dot_product;
77 : }
78 : /// @}
79 :
80 : /// @{
81 : /*!
82 : * \ingroup TensorGroup
83 : * \brief Compute the dot_product of two vectors or one forms
84 : *
85 : * \details
86 : * Returns \f$g_{ab} A^a B^b\f$, where \f$g_{ab}\f$ is the metric,
87 : * \f$A^a\f$ is vector_a, and \f$B^b\f$ is vector_b.
88 : * Or, returns \f$g^{ab} A_a B_b\f$ when given one forms \f$A_a\f$
89 : * and \f$B_b\f$ with an inverse metric \f$g^{ab}\f$.
90 : */
91 : template <typename DataType, typename Index>
92 1 : void dot_product(
93 : const gsl::not_null<Scalar<DataType>*> dot_product,
94 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
95 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_b,
96 : const Tensor<DataType, Symmetry<1, 1>,
97 : index_list<change_index_up_lo<Index>,
98 : change_index_up_lo<Index>>>& metric) {
99 : if constexpr (Index::dim == 1) {
100 : get(*dot_product) = get<0>(vector_a) * get<0>(vector_b) * get<0, 0>(metric);
101 : } else {
102 : if (&vector_a == &vector_b) {
103 : get(*dot_product) =
104 : get<0>(vector_a) * get<1>(vector_a) * get<0, 1>(metric);
105 : for (size_t j = 2; j < Index::dim; ++j) {
106 : get(*dot_product) +=
107 : vector_a.get(0) * vector_a.get(j) * metric.get(0, j);
108 : }
109 : for (size_t i = 1; i < Index::dim; ++i) {
110 : for (size_t j = i + 1; j < Index::dim; ++j) {
111 : get(*dot_product) +=
112 : vector_a.get(i) * vector_a.get(j) * metric.get(i, j);
113 : }
114 : }
115 : get(*dot_product) *= 2.0;
116 :
117 : for (size_t i = 0; i < Index::dim; ++i) {
118 : get(*dot_product) += square(vector_a.get(i)) * metric.get(i, i);
119 : }
120 : } else {
121 : get(*dot_product) =
122 : get<0>(vector_a) * get<0>(vector_b) * get<0, 0>(metric);
123 : for (size_t b = 1; b < Index::dim; ++b) {
124 : get(*dot_product) +=
125 : get<0>(vector_a) * vector_b.get(b) * metric.get(0, b);
126 : }
127 :
128 : for (size_t a = 1; a < Index::dim; ++a) {
129 : for (size_t b = 0; b < Index::dim; ++b) {
130 : get(*dot_product) +=
131 : vector_a.get(a) * vector_b.get(b) * metric.get(a, b);
132 : }
133 : }
134 : }
135 : }
136 : }
137 :
138 : template <typename DataType, typename Index>
139 1 : Scalar<DataType> dot_product(
140 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
141 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_b,
142 : const Tensor<DataType, Symmetry<1, 1>,
143 : index_list<change_index_up_lo<Index>,
144 : change_index_up_lo<Index>>>& metric) {
145 : Scalar<DataType> dot_product(get_size(get<0>(vector_a)));
146 : ::dot_product(make_not_null(&dot_product), vector_a, vector_b, metric);
147 : return dot_product;
148 : }
149 : /// @}
|
