CrossProduct.hpp
Go to the documentation of this file.
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"
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 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) noexcept {
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>>>
57  const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
58  const Tensor<DataType, Symmetry<1>, index_list<change_index_up_lo<Index>>>&
59  vector_b) noexcept {
60  static_assert(Index::dim == 3, "cross_product vectors must have dimension 3");
61  static_assert(Index::index_type == IndexType::Spatial,
62  "cross_product vectors must be spatial");
63 
65  Tensor<DataType, Symmetry<1>, index_list<change_index_up_lo<Index>>>>(
66  vector_b, 0.0);
67  for (LeviCivitaIterator<3> it; it; ++it) {
68  cross_product.get(it[0]) +=
69  it.sign() * vector_a.get(it[1]) * vector_b.get(it[2]);
70  }
71  return cross_product;
72 }
73 
74 /*!
75  * \ingroup TensorGroup
76  * \brief Compute the cross product of two vectors or one forms
77  *
78  * \details
79  * Returns \f$\sqrt{g} g^{li} A^j B^k \epsilon_{ljk}\f$, where
80  * \f$A^j\f$ and \f$B^k\f$ are vectors and \f$g^{li}\f$ and \f$g\f$
81  * are the inverse and determinant, respectively, of
82  * the spatial metric (computed via `determinant_and_inverse`). In this case,
83  * the arguments `vector_a` and `vector_b` should be vectors,
84  * the argument `metric_or_inverse_metric` should be the inverse spatial
85  * metric \f$g^{ij}\f$, and the argument `metric_determinant` should be the
86  * determinant of the spatial metric \f$\det(g_{ij})\f$.
87  * Or, returns \f$\sqrt{g}^{-1} g_{li} A_j B_k \epsilon^{ljk}\f$, where
88  * \f$A_j\f$ and \f$B_k\f$ are one forms and \f$g_{li}\f$ and \f$g\f$
89  * are the spatial metric and its determinant. In this case,
90  * the arguments `vector_a` and `vector_b` should be one forms,
91  * the argument `metric_or_inverse_metric` should be the spatial metric
92  * \f$g_{ij}\f$, and the argument `metric_determinant` should be the
93  * determinant of the spatial metric \f$\det(g_{ij})\f$.
94  */
95 template <typename DataType, typename Index>
96 Tensor<DataType, Symmetry<1>, index_list<Index>> cross_product(
97  const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
98  const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_b,
99  const Tensor<DataType, Symmetry<1, 1>, index_list<Index, Index>>&
100  metric_or_inverse_metric,
101  const Scalar<DataType>& metric_determinant) noexcept {
102  static_assert(Index::dim == 3, "cross_product vectors must have dimension 3");
103  static_assert(Index::index_type == IndexType::Spatial,
104  "cross_product vectors must be spatial");
105 
106  auto cross_product =
107  make_with_value<Tensor<DataType, Symmetry<1>, index_list<Index>>>(
108  vector_a, 0.);
109  for (size_t i = 0; i < Index::dim; ++i) {
110  for (LeviCivitaIterator<3> it; it; ++it) {
111  cross_product.get(i) += it.sign() * vector_a.get(it[1]) *
112  vector_b.get(it[2]) *
113  metric_or_inverse_metric.get(it[0], i);
114  }
115  if (Index::ul == UpLo::Up) {
116  cross_product.get(i) *= sqrt(get(metric_determinant));
117  } else {
118  cross_product.get(i) /= sqrt(get(metric_determinant));
119  }
120  }
121  return cross_product;
122 }
123 
124 /*!
125  * \ingroup TensorGroup
126  * \brief Compute the cross product of a vector and a one form
127  *
128  * \details
129  * Returns \f$\sqrt{g} A^j B_l g^{lk} \epsilon_{ijk}\f$ for input vector
130  * \f$A^j\f$ and input one form \f$B_l\f$. In this case,
131  * the argument `vector_a` should be a vector, `vector_b` should be a one form,
132  * `metric_or_inverse_metric` should be the inverse spatial
133  * metric \f$g^{ij}\f$, and `metric_determinant` should be the
134  * determinant of the spatial metric \f$\det(g_{ij})\f$.
135  * Or, returns \f$\sqrt{g}^{-1} A_j B^l g_{lk}
136  * \epsilon^{ijk}\f$ for input one form \f$A_j\f$ and input vector \f$B^l\f$.
137  * In this case,
138  * the argument `vector_a` should be a one form, `vector_b` should be a vector,
139  * `metric_or_inverse_metric` should be the spatial metric
140  * \f$g_{ij}\f$, and `metric_determinant` should be the
141  * determinant of the spatial metric \f$\det(g_{ij})\f$.
142  * Note that this function returns a vector if `vector_b` is a vector and a
143  * one form if `vector_b` is a one form.
144  */
145 template <typename DataType, typename Index>
146 Tensor<DataType, Symmetry<1>, index_list<change_index_up_lo<Index>>>
147 cross_product(const Tensor<DataType, Symmetry<1>, index_list<Index>>& vector_a,
148  const Tensor<DataType, Symmetry<1>,
149  index_list<change_index_up_lo<Index>>>& vector_b,
150  const Tensor<DataType, Symmetry<1, 1>, index_list<Index, Index>>&
151  metric_or_inverse_metric,
152  const Scalar<DataType>& metric_determinant) noexcept {
153  static_assert(Index::dim == 3, "cross_product vectors must have dimension 3");
154  static_assert(Index::index_type == IndexType::Spatial,
155  "cross_product vectors must be spatial");
156 
158  Tensor<DataType, Symmetry<1>, index_list<change_index_up_lo<Index>>>>(
159  vector_b, 0.0);
160  for (LeviCivitaIterator<3> it; it; ++it) {
161  if (Index::ul == UpLo::Up) {
162  for (size_t l = 0; l < Index::dim; ++l) {
163  cross_product.get(it[0]) += it.sign() * vector_a.get(it[1]) *
164  vector_b.get(l) *
165  metric_or_inverse_metric.get(it[2], l);
166  }
167  } else {
168  for (size_t l = 0; l < Index::dim; ++l) {
169  cross_product.get(it[0]) += it.sign() * vector_a.get(l) *
170  vector_b.get(it[2]) *
171  metric_or_inverse_metric.get(it[1], l);
172  }
173  }
174  }
175 
176  for (size_t i = 0; i < Index::dim; ++i) {
177  if (Index::ul == UpLo::Up) {
178  cross_product.get(i) *= sqrt(get(metric_determinant));
179  } else {
180  cross_product.get(i) /= sqrt(get(metric_determinant));
181  }
182  }
183  return cross_product;
184 }
The TensorIndexType is purely spatial.
Contravariant, or Upper index.
Iterate over all nonzero index permutations for a Levi-Civita symbol.
Definition: LeviCivitaIterator.hpp:33
typename detail::SymmetryImpl< std::make_index_sequence< sizeof...(T)>, tmpl::integral_list< std::int32_t, T... > >::type Symmetry
Computes the canonical symmetry from the integers T
Definition: Symmetry.hpp:81
std::remove_const_t< R > make_with_value(const T &input, const ValueType &value) noexcept
Given an object of type T, create an object of type R whose elements are initialized to value...
Definition: MakeWithValue.hpp:42
Defines classes for Tensor.
Tensor_detail::TensorIndexType< Index::index_type==IndexType::Spatial ? Index::value :Index::value - 1, Index::ul==UpLo::Up ? UpLo::Lo :UpLo::Up, typename Index::Frame, Index::index_type > change_index_up_lo
Change the TensorIndexType to be covariant if it&#39;s contravariant and vice-versa.
Definition: IndexType.hpp:233
Tensor< DataType, Symmetry< 1 >, index_list< Index > > cross_product(const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector_a, const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector_b) noexcept
Compute the Euclidean cross product of two vectors or one forms.
Definition: CrossProduct.hpp:26
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
Defines make_with_value.
Defines classes representing tensor indices.