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