Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines functions `cross_product` for flat and curved-space cross products.
6 :
7 : #pragma once
8 :
9 : #include <cstddef>
10 :
11 : #include "DataStructures/LeviCivitaIterator.hpp"
12 : #include "DataStructures/Tensor/IndexType.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "Utilities/MakeWithValue.hpp"
15 :
16 : /*!
17 : * \ingroup TensorGroup
18 : * \brief Compute the Euclidean cross product of two vectors or one forms
19 : *
20 : * \details
21 : * Returns \f$A^j B^k \epsilon_{ljk} \delta^{il}\f$ for input vectors \f$A^j\f$
22 : * and \f$B^k\f$ or \f$A_j B_k \epsilon^{ljk} \delta_{il}\f$ for input one
23 : * forms \f$A_j\f$ and \f$B_k\f$.
24 : */
25 : template <typename DataType, typename Index>
26 1 : Tensor<DataType, Symmetry<1>, index_list<Index>> cross_product(
27 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
28 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_b) {
29 : static_assert(Index::dim == 3, "cross_product vectors must have dimension 3");
30 : static_assert(Index::index_type == IndexType::Spatial,
31 : "cross_product vectors must be spatial");
32 :
33 : auto cross_product =
34 : make_with_value<Tensor<DataType, Symmetry<1>, index_list<Index>>>(
35 : vector_a, 0.0);
36 : for (LeviCivitaIterator<3> it; it; ++it) {
37 : cross_product.get(it[0]) +=
38 : it.sign() * vector_a.get(it[1]) * vector_b.get(it[2]);
39 : }
40 : return cross_product;
41 : }
42 :
43 : /*!
44 : * \ingroup TensorGroup
45 : * \brief Compute the Euclidean cross product of a vector and a one form
46 : *
47 : * \details
48 : * Returns \f$A^j B_l \delta^{lk} \epsilon_{ijk}\f$ for input vector \f$A^j\f$
49 : * and input one form \f$B_l\f$
50 : * or \f$A_j B^l \delta_{lk} \epsilon^{ijk}\f$ for input one form \f$A_j\f$ and
51 : * input vector \f$B^l\f$. Note that this function returns a vector if
52 : * `vector_b` is a vector and a one form if `vector_b` is a one form.
53 : */
54 : template <typename DataType, typename Index>
55 : Tensor<DataType, Symmetry<1>, index_list<change_index_up_lo<Index>>>
56 1 : cross_product(const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
57 : const Tensor<DataType, Symmetry<1>,
58 : index_list<change_index_up_lo<Index>>>& vector_b) {
59 : static_assert(Index::dim == 3, "cross_product vectors must have dimension 3");
60 : static_assert(Index::index_type == IndexType::Spatial,
61 : "cross_product vectors must be spatial");
62 :
63 : auto cross_product = make_with_value<
64 : Tensor<DataType, Symmetry<1>, index_list<change_index_up_lo<Index>>>>(
65 : vector_b, 0.0);
66 : for (LeviCivitaIterator<3> it; it; ++it) {
67 : cross_product.get(it[0]) +=
68 : it.sign() * vector_a.get(it[1]) * vector_b.get(it[2]);
69 : }
70 : return cross_product;
71 : }
72 :
73 : /*!
74 : * \ingroup TensorGroup
75 : * \brief Compute the cross product of two vectors or one forms
76 : *
77 : * \details
78 : * Returns \f$\sqrt{g} g^{li} A^j B^k \epsilon_{ljk}\f$, where
79 : * \f$A^j\f$ and \f$B^k\f$ are vectors and \f$g^{li}\f$ and \f$g\f$
80 : * are the inverse and determinant, respectively, of
81 : * the spatial metric (computed via `determinant_and_inverse`). In this case,
82 : * the arguments `vector_a` and `vector_b` should be vectors,
83 : * the argument `metric_or_inverse_metric` should be the inverse spatial
84 : * metric \f$g^{ij}\f$, and the argument `metric_determinant` should be the
85 : * determinant of the spatial metric \f$\det(g_{ij})\f$.
86 : * Or, returns \f$\sqrt{g}^{-1} g_{li} A_j B_k \epsilon^{ljk}\f$, where
87 : * \f$A_j\f$ and \f$B_k\f$ are one forms and \f$g_{li}\f$ and \f$g\f$
88 : * are the spatial metric and its determinant. In this case,
89 : * the arguments `vector_a` and `vector_b` should be one forms,
90 : * the argument `metric_or_inverse_metric` should be the spatial metric
91 : * \f$g_{ij}\f$, and the argument `metric_determinant` should be the
92 : * determinant of the spatial metric \f$\det(g_{ij})\f$.
93 : */
94 : template <typename DataType, typename Index>
95 1 : Tensor<DataType, Symmetry<1>, index_list<Index>> cross_product(
96 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
97 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_b,
98 : const Tensor<DataType, Symmetry<1, 1>, index_list<Index, Index>>&
99 : metric_or_inverse_metric,
100 : const Scalar<DataType>& metric_determinant) {
101 : static_assert(Index::dim == 3, "cross_product vectors must have dimension 3");
102 : static_assert(Index::index_type == IndexType::Spatial,
103 : "cross_product vectors must be spatial");
104 :
105 : auto cross_product =
106 : make_with_value<Tensor<DataType, Symmetry<1>, index_list<Index>>>(
107 : vector_a, 0.);
108 : for (size_t i = 0; i < Index::dim; ++i) {
109 : for (LeviCivitaIterator<3> it; it; ++it) {
110 : cross_product.get(i) += it.sign() * vector_a.get(it[1]) *
111 : vector_b.get(it[2]) *
112 : metric_or_inverse_metric.get(it[0], i);
113 : }
114 : if (Index::ul == UpLo::Up) {
115 : cross_product.get(i) *= sqrt(get(metric_determinant));
116 : } else {
117 : cross_product.get(i) /= sqrt(get(metric_determinant));
118 : }
119 : }
120 : return cross_product;
121 : }
122 :
123 : /*!
124 : * \ingroup TensorGroup
125 : * \brief Compute the cross product of a vector and a one form
126 : *
127 : * \details
128 : * Returns \f$\sqrt{g} A^j B_l g^{lk} \epsilon_{ijk}\f$ for input vector
129 : * \f$A^j\f$ and input one form \f$B_l\f$. In this case,
130 : * the argument `vector_a` should be a vector, `vector_b` should be a one form,
131 : * `metric_or_inverse_metric` should be the inverse spatial
132 : * metric \f$g^{ij}\f$, and `metric_determinant` should be the
133 : * determinant of the spatial metric \f$\det(g_{ij})\f$.
134 : * Or, returns \f$\sqrt{g}^{-1} A_j B^l g_{lk}
135 : * \epsilon^{ijk}\f$ for input one form \f$A_j\f$ and input vector \f$B^l\f$.
136 : * In this case,
137 : * the argument `vector_a` should be a one form, `vector_b` should be a vector,
138 : * `metric_or_inverse_metric` should be the spatial metric
139 : * \f$g_{ij}\f$, and `metric_determinant` should be the
140 : * determinant of the spatial metric \f$\det(g_{ij})\f$.
141 : * Note that this function returns a vector if `vector_b` is a vector and a
142 : * one form if `vector_b` is a one form.
143 : */
144 : template <typename DataType, typename Index>
145 : Tensor<DataType, Symmetry<1>, index_list<change_index_up_lo<Index>>>
146 1 : cross_product(const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
147 : const Tensor<DataType, Symmetry<1>,
148 : index_list<change_index_up_lo<Index>>>& vector_b,
149 : const Tensor<DataType, Symmetry<1, 1>, index_list<Index, Index>>&
150 : metric_or_inverse_metric,
151 : const Scalar<DataType>& metric_determinant) {
152 : static_assert(Index::dim == 3, "cross_product vectors must have dimension 3");
153 : static_assert(Index::index_type == IndexType::Spatial,
154 : "cross_product vectors must be spatial");
155 :
156 : auto cross_product = make_with_value<
157 : Tensor<DataType, Symmetry<1>, index_list<change_index_up_lo<Index>>>>(
158 : vector_b, 0.0);
159 : for (LeviCivitaIterator<3> it; it; ++it) {
160 : if (Index::ul == UpLo::Up) {
161 : for (size_t l = 0; l < Index::dim; ++l) {
162 : cross_product.get(it[0]) += it.sign() * vector_a.get(it[1]) *
163 : vector_b.get(l) *
164 : metric_or_inverse_metric.get(it[2], l);
165 : }
166 : } else {
167 : for (size_t l = 0; l < Index::dim; ++l) {
168 : cross_product.get(it[0]) += it.sign() * vector_a.get(l) *
169 : vector_b.get(it[2]) *
170 : metric_or_inverse_metric.get(it[1], l);
171 : }
172 : }
173 : }
174 :
175 : for (size_t i = 0; i < Index::dim; ++i) {
176 : if (Index::ul == UpLo::Up) {
177 : cross_product.get(i) *= sqrt(get(metric_determinant));
178 : } else {
179 : cross_product.get(i) /= sqrt(get(metric_determinant));
180 : }
181 : }
182 : return cross_product;
183 : }
184 :
185 : /*!
186 : * \brief Compute the cross product of three spacetime vectors or one forms
187 : *
188 : * \details
189 : * Returns $\sqrt{g} g^{ed} A^a B^b C^c \epsilon_{abcd}$, where $\{A,B,C\}$ are
190 : * the input vectors and $g^{ab}$ and $g$ are the inverse spacetime metric and
191 : * the spacetime metric determinant. Or, returns $-\sqrt{g}^{-1} g_{ed} A_a B_b
192 : * C_c \epsilon^{abcd}$, where $\{A,B,C\}$ are one forms and $g_{ab}$ and $g$
193 : * are the spacetime metric and the spacetime metric determinant.
194 : * See \cite Misner1973 page 202 for details.
195 : *
196 : * The output buffer must be sized correctly on entry.
197 : */
198 : template <typename DataType, typename Index>
199 1 : void cross_product(
200 : const gsl::not_null<Tensor<DataType, Symmetry<1>, index_list<Index>>*>
201 : result,
202 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
203 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_b,
204 : const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_c,
205 : const Tensor<DataType, Symmetry<1, 1>, index_list<Index, Index>>&
206 : metric_or_inverse_metric,
207 : const Scalar<DataType>& metric_determinant) {
208 : static_assert(Index::dim == 4,
209 : "cross_product vectors must have spacetime dimension 4");
210 : static_assert(Index::index_type == IndexType::Spacetime,
211 : "cross product of 3 vectors must be spacetime");
212 : for (size_t i = 0; i < Index::dim; ++i) {
213 : result->get(i) = 0.0;
214 : for (LeviCivitaIterator<4> it; it; ++it) {
215 : result->get(i) += it.sign() * vector_a.get(it[0]) * vector_b.get(it[1]) *
216 : vector_c.get(it[2]) *
217 : metric_or_inverse_metric.get(it[3], i);
218 : }
219 : if (Index::ul == UpLo::Up) {
220 : result->get(i) *= sqrt(abs(get(metric_determinant)));
221 : } else {
222 : result->get(i) /= -sqrt(abs(get(metric_determinant)));
223 : }
224 : }
225 : }
|
