diff --git a/math/fastor/include/algebra/math/impl/fastor_vector.hpp b/math/fastor/include/algebra/math/impl/fastor_vector.hpp index b50cbe38..7cc68cec 100644 --- a/math/fastor/include/algebra/math/impl/fastor_vector.hpp +++ b/math/fastor/include/algebra/math/impl/fastor_vector.hpp @@ -8,6 +8,7 @@ #pragma once // Project include(s). +#include "algebra/concepts.hpp" #include "algebra/qualifiers.hpp" // Fastor include(s). @@ -21,30 +22,56 @@ namespace algebra::fastor::math { +// Note that for all `Fastor::AbstractTensors`, the `N` template parameters +// refers to the number of dimensions, not the number of elements. + /// This method retrieves phi from a vector @param v -template -requires(N >= 2) ALGEBRA_HOST_DEVICE - constexpr auto phi(const Fastor::Tensor &v) { +template +ALGEBRA_HOST_DEVICE constexpr auto phi( + const Fastor::AbstractTensor &a) { + // we first force evaluation of whatever was passed in. + auto v = Fastor::evaluate(a); + // using `auto` relieves us from having to extract the exact dimension of the + // vector from somewhere. For all intents and purposes, we can consider the + // type to be `Fastor::Tensor`. + + // we use the cmath version of `atan2` because Fastor's `atan2` works on + // `Fastor::Tensor`s element-wise, which we don't want. return algebra::math::atan2(v[1], v[0]); } /// This method retrieves theta from a vector, vector base with rows >= 3 /// /// @param v the input vector -template -requires(N >= 3) ALGEBRA_HOST constexpr scalar_t - theta(const Fastor::Tensor &v) noexcept { - +template +ALGEBRA_HOST constexpr auto theta( + const Fastor::AbstractTensor &a) noexcept { + // we first force evaluation of whatever was passed in. + auto v = Fastor::evaluate(a); + // using `auto` relieves us from having to extract the exact dimension of the + // vector from somewhere. For all intents and purposes, we can consider the + // type to be `Fastor::Tensor`. + + // we use the cmath version of `atan2` because Fastor's `atan2` works on + // `Fastor::Tensor`s element-wise, which we don't want. return algebra::math::atan2(Fastor::norm(v(Fastor::fseq<0, 2>())), v[2]); } /// This method retrieves the perpendicular magnitude of a vector with rows >= 2 /// /// @param v the input vector -template -requires(N >= 2) ALGEBRA_HOST constexpr scalar_t - perp(const Fastor::Tensor &v) noexcept { - +template +ALGEBRA_HOST constexpr auto perp( + const Fastor::AbstractTensor &a) noexcept { + + // we first force evaluation of whatever was passed in. + // using `auto` relieves us from having to extract the exact dimension of the + // vector from somewhere. For all intents and purposes, we can consider the + // type to be `Fastor::Tensor`. + auto v = Fastor::evaluate(a); + + // we use the cmath version of `sqrt` because Fastor's `sqrt` works on + // `Fastor::Tensor`s element-wise, which we don't want. return algebra::math::sqrt( Fastor::inner(v(Fastor::fseq<0, 2>()), v(Fastor::fseq<0, 2>()))); } @@ -52,8 +79,8 @@ requires(N >= 2) ALGEBRA_HOST constexpr scalar_t /// This method retrieves the norm of a vector, no dimension restriction /// /// @param v the input vector -template -ALGEBRA_HOST constexpr scalar_t norm(const Fastor::Tensor &v) { +template +ALGEBRA_HOST constexpr auto norm(const Fastor::AbstractTensor &v) { return Fastor::norm(v); } @@ -62,21 +89,24 @@ ALGEBRA_HOST constexpr scalar_t norm(const Fastor::Tensor &v) { /// rows >= 3 /// /// @param v the input vector -template -requires(N >= 3) ALGEBRA_HOST constexpr scalar_t - eta(const Fastor::Tensor &v) noexcept { +template +ALGEBRA_HOST constexpr auto eta( + const Fastor::AbstractTensor &a) noexcept { + auto v = Fastor::evaluate(a); + // we use the cmath version of `atanh` because Fastor's `atanh` works on + // `Fastor::Tensor`s element-wise, which we don't want. return algebra::math::atanh(v[2] / Fastor::norm(v)); } /// Get a normalized version of the input vector /// /// @param v the input vector -template -ALGEBRA_HOST constexpr Fastor::Tensor normalize( - const Fastor::Tensor &v) { +template +ALGEBRA_HOST constexpr auto normalize( + const Fastor::AbstractTensor &v) { - return (static_cast(1.0) / Fastor::norm(v)) * v; + return v / Fastor::norm(v); } /// Dot product between two input vectors @@ -85,6 +115,18 @@ ALGEBRA_HOST constexpr Fastor::Tensor normalize( /// @param b the second input vector /// /// @return the scalar dot product value +template +ALGEBRA_HOST constexpr auto dot(const Fastor::AbstractTensor &a, + const Fastor::AbstractTensor &b) { + return Fastor::inner(a, b); +} + +/// Dot product between two pure vectors (of type `Tensor`) +/// +/// @param a the first input vector +/// @param b the second input vector +/// +/// @return the scalar dot product value template ALGEBRA_HOST_DEVICE constexpr scalar_t dot( const Fastor::Tensor &a, @@ -92,10 +134,10 @@ ALGEBRA_HOST_DEVICE constexpr scalar_t dot( return Fastor::inner(a, b); } -/// Dot product between Tensor and Tensor +/// Dot product between a vector and a matrix slice /// -/// @param a the first input vector -/// @param b the second input Tensor +/// @param a the first input: a vector (`Tensor`) +/// @param b the second input: a matrix (`Tensor`) /// /// @return the scalar dot product value template @@ -109,10 +151,10 @@ ALGEBRA_HOST constexpr scalar_t dot(const Fastor::Tensor &a, Fastor::Tensor(b(Fastor::fseq<0, N>(), 0))); } -/// Dot product between Tensor and Tensor +/// Dot product between a matrix slice and a vector /// -/// @param a the second input Tensor -/// @param b the first input vector +/// @param a the first input: a matrix (`Tensor`) +/// @param b the second input: a vector (`Tensor`) /// /// @return the scalar dot product value template @@ -123,10 +165,10 @@ ALGEBRA_HOST constexpr scalar_t dot(const Fastor::Tensor &a, b); } -/// Dot product between two Tensor +/// Dot product between two matrix slices /// -/// @param a the second input Tensor -/// @param b the first input Tensor +/// @param a the first input: a matrix (`Tensor`) +/// @param b the second input: a matrix (`Tensor`) /// /// @return the scalar dot product value template @@ -143,23 +185,34 @@ ALGEBRA_HOST constexpr scalar_t dot(const Fastor::Tensor &a, /// @param b the second input vector /// /// @return a vector (expression) representing the cross product -template -ALGEBRA_HOST_DEVICE constexpr Fastor::Tensor cross( - const Fastor::Tensor &a, - const Fastor::Tensor &b) { +template +ALGEBRA_HOST constexpr auto cross( + const Fastor::AbstractTensor &a, + const Fastor::AbstractTensor &b) { return Fastor::cross(a, b); } -/// Cross product between Tensor and Tensor +/// Cross product between two pure vectors (of type `Tensor`) /// /// @param a the first input vector -/// @param b the second input Tensor +/// @param b the second input vector /// -/// @return a vector representing the cross product +/// @return a vector (expression) representing the cross product template -ALGEBRA_HOST constexpr Fastor::Tensor cross( - const Fastor::Tensor &a, - const Fastor::Tensor &b) { +ALGEBRA_HOST constexpr auto cross(const Fastor::Tensor &a, + const Fastor::Tensor &b) { + return Fastor::cross(a, b); +} + +/// Cross product between a vector and a matrix slice +/// +/// @param a the first input: a vector (`Tensor`) +/// @param b the second input: a matrix (`Tensor`) +/// +/// @return a vector (expression) representing the cross product +template +ALGEBRA_HOST constexpr auto cross(const Fastor::Tensor &a, + const Fastor::Tensor &b) { // We need to specify the type of the Tensor slice because Fastor by default // is lazy, so it returns an intermediate type which does not play well with @@ -168,31 +221,29 @@ ALGEBRA_HOST constexpr Fastor::Tensor cross( Fastor::Tensor(b(Fastor::fseq<0, 3>(), 0))); } -/// Cross product between Tensor and Tensor +/// Cross product between a matrix slice and a vector /// -/// @param a the second input Tensor -/// @param b the first input vector +/// @param a the first input: a matrix (`Tensor`) +/// @param b the second input: a vector (`Tensor`) /// -/// @return a vector representing the cross product +/// @return a vector (expression) representing the cross product template -ALGEBRA_HOST constexpr Fastor::Tensor cross( - const Fastor::Tensor &a, - const Fastor::Tensor &b) { +ALGEBRA_HOST constexpr auto cross(const Fastor::Tensor &a, + const Fastor::Tensor &b) { return Fastor::cross(Fastor::Tensor(a(Fastor::fseq<0, 3>(), 0)), b); } -/// Cross product between two Tensor +/// Cross product between two matrix slices /// -/// @param a the second input Tensor -/// @param b the first input Tensor +/// @param a the second input matrix (`Tensor`) +/// @param b the first input matrix (`Tensor`) /// -/// @return a vector representing the cross product +/// @return a vector (expression) representing the cross product template -ALGEBRA_HOST constexpr Fastor::Tensor cross( - const Fastor::Tensor &a, - const Fastor::Tensor &b) { +ALGEBRA_HOST constexpr auto cross(const Fastor::Tensor &a, + const Fastor::Tensor &b) { return Fastor::cross(Fastor::Tensor(a(Fastor::fseq<0, 3>(), 0)), Fastor::Tensor(b(Fastor::fseq<0, 3>(), 0))); diff --git a/tests/fastor/fastor_fastor.cpp b/tests/fastor/fastor_fastor.cpp index cb4a1147..f2cfb61f 100644 --- a/tests/fastor/fastor_fastor.cpp +++ b/tests/fastor/fastor_fastor.cpp @@ -33,9 +33,46 @@ struct test_specialisation_name { } }; +// This test checks to see if the `dot` function can handle when one of its +// operands is a sum or difference of two vectors. We also test to see if the +// `dot` function plays nicely with scalar multiplication (just to be safe). +TYPED_TEST_P(test_host_basics_vector, dot_product_with_ops) { + typename TypeParam::vector3 v1{1.f, 2.f, 3.f}; + typename TypeParam::vector3 v2{3.f, 4.f, 5.f}; + + ASSERT_NEAR(algebra::vector::dot(v1 + v2, v2), 76.f, this->m_epsilon); + ASSERT_NEAR(algebra::vector::dot(v1, v2 - v1), 12.f, this->m_epsilon); + ASSERT_NEAR(algebra::vector::dot(v1 + v2, v1 - v2), -36.f, this->m_epsilon); + ASSERT_NEAR(algebra::vector::dot(v1 + v2, 2 * v2), 152.f, this->m_epsilon); +} + +// This test checks to see if the `cross` function can handle when one of its +// operands is a sum or difference of two vectors. We also test to see if the +// `cross` function plays nicely with scalar multiplication (just to be safe). +TYPED_TEST_P(test_host_basics_vector, cross_product_add_sub) { + typename TypeParam::vector3 v1{1.f, 2.f, 3.f}; + typename TypeParam::vector3 v2{3.f, 4.f, 5.f}; + typename TypeParam::vector3 v3{-6.f, 7.f, -9.f}; + + typename TypeParam::vector3 v = algebra::vector::cross(v1 + v2, v3); + typename TypeParam::vector3 ans{-110.f, -12.f, 64.f}; + + ASSERT_NEAR(v[0], ans[0], this->m_epsilon); + ASSERT_NEAR(v[1], ans[1], this->m_epsilon); + ASSERT_NEAR(v[2], ans[2], this->m_epsilon); + + v = algebra::vector::cross(v3 - 2 * v1, 3 * (v1 + v2)); + ans = {342.f, 12.f, -180.f}; + + ASSERT_NEAR(v[0], ans[0], this->m_epsilon); + ASSERT_NEAR(v[1], ans[1], this->m_epsilon); + ASSERT_NEAR(v[2], ans[2], this->m_epsilon); +} + // Register the tests REGISTER_TYPED_TEST_SUITE_P(test_host_basics_vector, local_vectors, vector3, - getter); + getter, dot_product_with_ops, + cross_product_add_sub); TEST_HOST_BASICS_MATRIX_TESTS(); REGISTER_TYPED_TEST_SUITE_P(test_host_basics_transform, transform3, global_transformations);