Add function definitions for function calls involving intermediate Fastor types (#162)

Use abstract tensor types in the Fastor plugin

Author: wermos

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,39 +22,65 @@
 
 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 <concepts::scalar scalar_t, auto N>
-requires(N >= 2) ALGEBRA_HOST_DEVICE
-    constexpr auto phi(const Fastor::Tensor<scalar_t, N> &v) {
+template <typename Derived, auto N>
+ALGEBRA_HOST_DEVICE constexpr auto phi(
+    const Fastor::AbstractTensor<Derived, N> &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<scalar_t, SIZE>`.
+
+  // 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 <concepts::scalar scalar_t, auto N>
-requires(N >= 3) ALGEBRA_HOST constexpr scalar_t
-    theta(const Fastor::Tensor<scalar_t, N> &v) noexcept {
-
+template <typename Derived, auto N>
+ALGEBRA_HOST constexpr auto theta(
+    const Fastor::AbstractTensor<Derived, N> &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<scalar_t, SIZE>`.
+
+  // 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 <concepts::scalar scalar_t, auto N>
-requires(N >= 2) ALGEBRA_HOST constexpr scalar_t
-    perp(const Fastor::Tensor<scalar_t, N> &v) noexcept {
-
+template <typename Derived, auto N>
+ALGEBRA_HOST constexpr auto perp(
+    const Fastor::AbstractTensor<Derived, N> &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<scalar_t, SIZE>`.
+  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>())));
 }
 
 /// This method retrieves the norm of a vector, no dimension restriction
 ///
 /// @param v the input vector
-template <concepts::scalar scalar_t, auto N>
-ALGEBRA_HOST constexpr scalar_t norm(const Fastor::Tensor<scalar_t, N> &v) {
+template <typename Derived, auto N>
+ALGEBRA_HOST constexpr auto norm(const Fastor::AbstractTensor<Derived, N> &v) {
 
   return Fastor::norm(v);
 }
@@ -62,21 +89,24 @@ ALGEBRA_HOST constexpr scalar_t norm(const Fastor::Tensor<scalar_t, N> &v) {
 /// rows >= 3
 ///
 /// @param v the input vector
-template <concepts::scalar scalar_t, auto N>
-requires(N >= 3) ALGEBRA_HOST constexpr scalar_t
-    eta(const Fastor::Tensor<scalar_t, N> &v) noexcept {
+template <typename Derived, auto N>
+ALGEBRA_HOST constexpr auto eta(
+    const Fastor::AbstractTensor<Derived, N> &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 <concepts::scalar scalar_t, auto N>
-ALGEBRA_HOST constexpr Fastor::Tensor<scalar_t, N> normalize(
-    const Fastor::Tensor<scalar_t, N> &v) {
+template <typename Derived, auto N>
+ALGEBRA_HOST constexpr auto normalize(
+    const Fastor::AbstractTensor<Derived, N> &v) {
 
-  return (static_cast<scalar_t>(1.0) / Fastor::norm(v)) * v;
+  return v / Fastor::norm(v);
 }
 
 /// Dot product between two input vectors
@@ -85,17 +115,29 @@ ALGEBRA_HOST constexpr Fastor::Tensor<scalar_t, N> normalize(
 /// @param b the second input vector
 ///
 /// @return the scalar dot product value
+template <typename Derived0, auto N0, typename Derived1, auto N1>
+ALGEBRA_HOST constexpr auto dot(const Fastor::AbstractTensor<Derived0, N0> &a,
+                                const Fastor::AbstractTensor<Derived1, N1> &b) {
+  return Fastor::inner(a, b);
+}
+
+/// Dot product between two pure vectors (of type `Tensor<scalar_t, N>`)
+///
+/// @param a the first input vector
+/// @param b the second input vector
+///
+/// @return the scalar dot product value
 template <concepts::scalar scalar_t, auto N>
 ALGEBRA_HOST_DEVICE constexpr scalar_t dot(
     const Fastor::Tensor<scalar_t, N> &a,
     const Fastor::Tensor<scalar_t, N> &b) {
   return Fastor::inner(a, b);
 }
 
-/// Dot product between Tensor<scalar_t, N> and Tensor<scalar_t, N, 1>
+/// Dot product between a vector and a matrix slice
 ///
-/// @param a the first input vector
-/// @param b the second input Tensor<scalar_t, N, 1>
+/// @param a the first input: a vector (`Tensor<scalar_t, N>`)
+/// @param b the second input: a matrix (`Tensor<scalar_t, N, 1>`)
 ///
 /// @return the scalar dot product value
 template <concepts::scalar scalar_t, auto N>
@@ -109,10 +151,10 @@ ALGEBRA_HOST constexpr scalar_t dot(const Fastor::Tensor<scalar_t, N> &a,
                        Fastor::Tensor<scalar_t, N>(b(Fastor::fseq<0, N>(), 0)));
 }
 
-/// Dot product between Tensor<scalar_t, N> and Tensor<scalar_t, N, 1>
+/// Dot product between a matrix slice and a vector
 ///
-/// @param a the second input Tensor<scalar_t, N, 1>
-/// @param b the first input vector
+/// @param a the first input: a matrix (`Tensor<scalar_t, N, 1>`)
+/// @param b the second input: a vector (`Tensor<scalar_t, N>`)
 ///
 /// @return the scalar dot product value
 template <concepts::scalar scalar_t, auto N>
@@ -123,10 +165,10 @@ ALGEBRA_HOST constexpr scalar_t dot(const Fastor::Tensor<scalar_t, N, 1> &a,
                        b);
 }
 
-/// Dot product between two Tensor<scalar_t, 3, 1>
+/// Dot product between two matrix slices
 ///
-/// @param a the second input Tensor<scalar_t, 3, 1>
-/// @param b the first input Tensor<scalar_t, 3, 1>
+/// @param a the first input: a matrix (`Tensor<scalar_t, N, 1>`)
+/// @param b the second input: a matrix (`Tensor<scalar_t, N, 1>`)
 ///
 /// @return the scalar dot product value
 template <concepts::scalar scalar_t, auto N>
@@ -143,23 +185,34 @@ ALGEBRA_HOST constexpr scalar_t dot(const Fastor::Tensor<scalar_t, N, 1> &a,
 /// @param b the second input vector
 ///
 /// @return a vector (expression) representing the cross product
-template <concepts::scalar scalar_t>
-ALGEBRA_HOST_DEVICE constexpr Fastor::Tensor<scalar_t, 3> cross(
-    const Fastor::Tensor<scalar_t, 3> &a,
-    const Fastor::Tensor<scalar_t, 3> &b) {
+template <typename Derived0, auto N0, typename Derived1, auto N1>
+ALGEBRA_HOST constexpr auto cross(
+    const Fastor::AbstractTensor<Derived0, N0> &a,
+    const Fastor::AbstractTensor<Derived1, N1> &b) {
   return Fastor::cross(a, b);
 }
 
-/// Cross product between Tensor<scalar_t, 3> and Tensor<scalar_t, 3, 1>
+/// Cross product between two pure vectors (of type `Tensor<scalar_t, 3>`)
 ///
 /// @param a the first input vector
-/// @param b the second input Tensor<scalar_t, 3, 1>
+/// @param b the second input vector
 ///
-/// @return a vector representing the cross product
+/// @return a vector (expression) representing the cross product
 template <concepts::scalar scalar_t>
-ALGEBRA_HOST constexpr Fastor::Tensor<scalar_t, 3> cross(
-    const Fastor::Tensor<scalar_t, 3> &a,
-    const Fastor::Tensor<scalar_t, 3, 1> &b) {
+ALGEBRA_HOST constexpr auto cross(const Fastor::Tensor<scalar_t, 3> &a,
+                                  const Fastor::Tensor<scalar_t, 3> &b) {
+  return Fastor::cross(a, b);
+}
+
+/// Cross product between a vector and a matrix slice
+///
+/// @param a the first input: a vector (`Tensor<scalar_t, 3>`)
+/// @param b the second input: a matrix (`Tensor<scalar_t, 3, 1>`)
+///
+/// @return a vector (expression) representing the cross product
+template <concepts::scalar scalar_t>
+ALGEBRA_HOST constexpr auto cross(const Fastor::Tensor<scalar_t, 3> &a,
+                                  const Fastor::Tensor<scalar_t, 3, 1> &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<scalar_t, 3> cross(
                        Fastor::Tensor<scalar_t, 3>(b(Fastor::fseq<0, 3>(), 0)));
 }
 
-/// Cross product between Tensor<scalar_t, 3> and Tensor<scalar_t, 3, 1>
+/// Cross product between a matrix slice and a vector
 ///
-/// @param a the second input Tensor<scalar_t, 3, 1>
-/// @param b the first input vector
+/// @param a the first input: a matrix (`Tensor<scalar_t, 3, 1>`)
+/// @param b the second input: a vector (`Tensor<scalar_t, 3>`)
 ///
-/// @return a vector representing the cross product
+/// @return a vector (expression) representing the cross product
 template <concepts::scalar scalar_t>
-ALGEBRA_HOST constexpr Fastor::Tensor<scalar_t, 3> cross(
-    const Fastor::Tensor<scalar_t, 3, 1> &a,
-    const Fastor::Tensor<scalar_t, 3> &b) {
+ALGEBRA_HOST constexpr auto cross(const Fastor::Tensor<scalar_t, 3, 1> &a,
+                                  const Fastor::Tensor<scalar_t, 3> &b) {
 
   return Fastor::cross(Fastor::Tensor<scalar_t, 3>(a(Fastor::fseq<0, 3>(), 0)),
                        b);
 }
 
-/// Cross product between two Tensor<scalar_t, 3, 1>
+/// Cross product between two matrix slices
 ///
-/// @param a the second input Tensor<scalar_t, 3, 1>
-/// @param b the first input Tensor<scalar_t, 3, 1>
+/// @param a the second input matrix (`Tensor<scalar_t, 3, 1>`)
+/// @param b the first input matrix (`Tensor<scalar_t, 3, 1>`)
 ///
-/// @return a vector representing the cross product
+/// @return a vector (expression) representing the cross product
 template <concepts::scalar scalar_t>
-ALGEBRA_HOST constexpr Fastor::Tensor<scalar_t, 3> cross(
-    const Fastor::Tensor<scalar_t, 3, 1> &a,
-    const Fastor::Tensor<scalar_t, 3, 1> &b) {
+ALGEBRA_HOST constexpr auto cross(const Fastor::Tensor<scalar_t, 3, 1> &a,
+                                  const Fastor::Tensor<scalar_t, 3, 1> &b) {
 
   return Fastor::cross(Fastor::Tensor<scalar_t, 3>(a(Fastor::fseq<0, 3>(), 0)),
                        Fastor::Tensor<scalar_t, 3>(b(Fastor::fseq<0, 3>(), 0)));

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);