Add multipole moments, and a unit test for their accuracy

benchmarks/CMakeLists.txt

@@ -17,7 +17,12 @@ if (NOT TARGET symtensor)
     find_package(symtensor REQUIRED)
 endif ()
 
-add_executable(benchmarks symmetricTensor.cpp multipole.cpp)
+add_executable(
+        benchmarks
+        symmetricTensor.cpp
+        multipole.cpp
+        multipoleMoment.cpp
+)
 target_link_libraries(benchmarks PRIVATE
         Catch2::Catch2WithMain
         symtensor::symtensor

benchmarks/multipoleMoment.cpp

@@ -0,0 +1,10 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/benchmark/catch_benchmark_all.hpp>
+
+#include <symtensor/Multipole.h>
+
+using namespace symtensor;
+
+TEST_CASE("benchmark: Multipole moment gravity approximation", "[MultipoleMoment]") {
+    // todo
+}

include/symtensor/Index.h

@@ -67,10 +67,10 @@ namespace symtensor {
         template<>
         struct IndexTypeForDimension<4> {
             enum class type : std::size_t {
-                W = 0,
-                X,
+                X = 0,
                 Y,
-                Z
+                Z,
+                W,
             };
         };
 

include/symtensor/MultipoleBase.h

@@ -116,7 +116,7 @@ namespace symtensor {
         template<std::size_t R>
         using TensorType = std::tuple_element_t<indexForRank(R), TensorTuple>;
 
-        using Ranks = std::integer_sequence<std::size_t, Tensors::Rank...>;
+        using Ranks = std::integer_sequence<std::size_t, Rank<Tensors>()...>;
 
     private:
 
@@ -141,7 +141,8 @@ namespace symtensor {
          * @param types The values for each tensor contained in the multipole.
          *   So long as each tensor type has an implicit constructor, initializer-list syntax works.
          */
-        explicit inline constexpr MultipoleBase(Tensors &&...types) : _tuple(std::forward<Tensors>(types)...) {}
+        explicit inline constexpr MultipoleBase(const Tensors &...types) : _tuple(types...) {}
+        //explicit inline constexpr MultipoleBase(Tensors &&...types) : _tuple(std::forward<Tensors>(types)...) {}
 
         /**
          * @brief Constructor from a single vector

include/symtensor/MultipoleMoment.h

@@ -0,0 +1,111 @@
+/**
+ * @file
+ * @brief Provides the generic multipole type alongside typedefs for common uses.
+ */
+#ifndef SYMTENSOR_MULTIPOLEMOMENT_H
+#define SYMTENSOR_MULTIPOLEMOMENT_H
+
+#include <symtensor/MultipoleBase.h>
+
+namespace symtensor {
+
+    /**
+     * @brief Generic multipole type
+     *
+     * @tparam Order the rank of the highest symmetric tensor
+     * @tparam Tensors prefix of types which make up the lowest ranks of the multipole.
+     *  Must end in a Symmetric tensor -- the final type is "promoted" to produce higher rank tensors
+     *  until the requested Order is reached.
+     */
+    template<std::size_t Order, typename ...Tensors>
+    class MultipoleMoment : public MultipoleBaseFromTuple<
+            MultipoleMoment<Order, Tensors...>,
+            TensorSequence<Order, typename last_type<Tensors...>::Scalar, Tensors...>
+    > {
+
+        using Base = MultipoleBaseFromTuple<
+                MultipoleMoment<Order, Tensors...>,
+                TensorSequence<Order, typename last_type<Tensors...>::Scalar, Tensors...>
+        >;
+
+        using Self = MultipoleMoment;
+
+        using Scalar = typename last_type<Tensors...>::Scalar;
+
+    public:
+
+        using Base::Base;
+
+        template<typename S, indexable Vector>
+        static inline constexpr MultipoleMoment FromPointMass(const S &mass, const Vector &position) {
+            return [&]<std::size_t... i>(std::index_sequence<i...>) constexpr {
+                return Self{
+                        mass,
+                        std::tuple_element_t<i + 1, typename Self::TensorTuple>::CartesianPower(position)...
+                };
+            }(std::make_index_sequence<Self::NumTensors - 1>());
+        }
+
+        inline constexpr const Scalar &scalar() const { return this->template tensor<0>(); }
+
+        inline constexpr Scalar &scalar() { return this->template tensor<0>(); }
+
+        template<typename T>
+        inline constexpr Self &operator/=(const T &scalar) {
+            [&]<std::size_t... i>(std::index_sequence<i...>) constexpr {
+                ((this->template tensor<i + 1>() /= scalar), ...);
+            }(std::make_index_sequence<Self::NumTensors - 1>());
+            return *this;
+        }
+
+        template<typename T>
+        inline constexpr Self &operator*=(const T &scalar) {
+            [&]<std::size_t... i>(std::index_sequence<i...>) constexpr {
+                ((this->template tensor<i + 1>() *= scalar), ...);
+            }(std::make_index_sequence<Self::NumTensors - 1>());
+            return *this;
+        }
+
+
+    };
+
+    template<std::size_t Order>
+    using MultipoleMoment2f = MultipoleMoment<Order, SymmetricTensor2f<1>>;
+
+    template<std::size_t Order>
+    using MultipoleMoment3f = MultipoleMoment<Order, SymmetricTensor3f<1>>;
+
+    template<typename ...Tensors>
+    class Quadrupole : public MultipoleMoment<2, Tensors...> {
+        using Base = MultipoleMoment<2, Tensors...>;
+    public:
+        using Base::Base;
+    };
+
+    /// 2d quadrupole moment with floating point elements
+    using QuadrupoleMoment2f = MultipoleMoment2f<2>;
+
+    /// 3d quadrupole moment with floating point elements
+    using QuadrupoleMoment3f = MultipoleMoment3f<2>;
+
+    /// 2d octupole moment with floating point elements
+    using OctupoleMoment2f = MultipoleMoment2f<3>;
+
+    /// 3d octupole moment with floating point elements
+    using OctupoleMoment3f = MultipoleMoment3f<3>;
+
+    /// 2d hexadecupole moment with floating point elements
+    using HexadecupoleMoment2f = MultipoleMoment2f<4>;
+
+    /// 3d hexadecupole moment with floating point elements
+    using HexadecupoleMoment3f = MultipoleMoment3f<4>;
+
+    /// 2d triacontadyupole moment with floating point elements
+    using TriacontadyupoleMoment2f = MultipoleMoment2f<5>;
+
+    /// 3d triacontadyupole moment with floating point elements
+    using TriacontadyupoleMoment3f = MultipoleMoment3f<5>;
+
+}
+
+#endif //SYMTENSOR_MULTIPOLEMOMENT_H

include/symtensor/SymmetricTensorBase.h

@@ -482,6 +482,10 @@ namespace symtensor {
         /// @copydoc operator*=(const Scalar &)
         inline constexpr Implementation operator*(const Scalar &scalar) const { return Self{*this} *= scalar; }
 
+        friend inline constexpr Implementation operator*(const Scalar &scalar, const Self &tensor) {
+            return tensor * scalar;
+        }
+
         /// @copydoc operator/=(const Scalar &)
         inline constexpr Implementation operator/(const Scalar &scalar) const { return Self{*this} /= scalar; }
 
@@ -524,6 +528,34 @@ namespace symtensor {
             return Self{*this} -= other;
         }
 
+        /// @}
+    public:
+        /// @name Other Tensor-tensor operations
+        /// @{
+
+        template<symmetric_tensor<D, I> OtherTensor>
+        friend auto operator*(
+                const Implementation &lhs,
+                const OtherTensor &rhs
+        ) {
+
+
+            constexpr std::size_t OtherRank = OtherTensor::Rank;
+            constexpr std::size_t ProductRank = Rank - OtherRank;
+            using ProductTensor = ReplaceRank<Implementation, ProductRank>;
+
+            // todo: is this correct?
+            ProductTensor product{};
+            [&]<std::size_t... i>(std::index_sequence<i...>) {
+                ((
+                        product.template at<head<Index, Rank, ProductRank>(lexicographicalIndices<i>())>() +=
+                                lhs.template at<lexicographicalIndices<i>()>() *
+                                rhs.template at<(tail<Index, Rank, OtherRank>(lexicographicalIndices<i>()))>()
+                ), ...);
+            }(std::make_index_sequence<NumValues>());
+            return product;
+        }
+
         /// @}
     public:
         /// @name Comparison operations

include/symtensor/concepts.h

@@ -35,11 +35,10 @@ namespace symtensor {
     //    template<typename T>
     //    concept compile_time_indexable = requires{ compile_time_indexable_helper<T, decltype(std::index_sequence<T::Dimensions>)> };
 
-    template<typename T, std::size_t R, std::size_t D, typename I>
+    template<typename T, std::size_t D, typename I>
     concept symmetric_tensor = requires(T tensor){
 
-//        requires T::Dimensions == D;
-//        requires T::Rank == R;
+        requires T::Dimensions == D;
         requires std::same_as<typename T::Index, I>;
         // todo: check for at<>
     };

include/symtensor/glm.h

@@ -7,6 +7,17 @@
 
 namespace symtensor {
 
+    template<class Implementation, typename Scalar, std::size_t D, typename I>
+    constexpr static inline glm::vec<int(D), Scalar>
+    to_glm(const symtensor::SymmetricTensorBase<Implementation, Scalar, D, 1, I> &vector) {
+
+        glm::vec<int(D), Scalar> result{};
+        for (std::size_t i = 0; i < std::size_t(D); ++i) {
+            result[i] = vector[i];
+        }
+        return result;
+    }
+
     template<class Implementation, typename Scalar, std::size_t D, typename I>
     constexpr static inline glm::mat<int(D), int(D), Scalar>
     to_glm(const symtensor::SymmetricTensorBase<Implementation, Scalar, D, 2, I> &symmetricMatrix) {

include/symtensor/util.h

@@ -20,6 +20,14 @@ namespace symtensor {
             return V * pow(V, P - 1);
     }
 
+    template<std::size_t P, typename T>
+    inline T pow(T value) {
+        if constexpr (P == 0)
+            return T(1);
+        else
+            return value * pow<P - 1>(value);
+    }
+
     constexpr std::size_t pascal(std::size_t x, std::size_t y) {
         if (x == 0 || y == 0) return 1;
         else return pascal(x - 1, y) + pascal(x, y - 1);
@@ -40,26 +48,26 @@ namespace symtensor {
         return array;
     }
 
-    template<typename T, std::size_t N>
+    template<typename T, std::size_t N, std::size_t NewN = N - 1>
     static constexpr auto head(std::array<T, N> array) {
         if constexpr (N == 0)
             // Never try to get the head of an empty array! This makes the compiler very angry.
             return array;
         else {
-            std::array<T, N - 1> copy{};
-            std::copy(array.begin(), array.end() - 1, copy.begin());
+            std::array<T, NewN> copy{};
+            std::copy(array.begin(), array.begin() + NewN, copy.begin());
             return copy;
         }
     }
 
-    template<typename T, std::size_t N>
+    template<typename T, std::size_t N, std::size_t NewN = N - 1>
     static constexpr auto tail(std::array<T, N> array) {
         if constexpr (N == 0)
             // Never try to get the tail of an empty array! This makes the compiler very angry.
             return array;
         else {
-            std::array<T, N - 1> copy{};
-            std::copy(array.begin() + 1, array.end(), copy.begin());
+            std::array<T, NewN> copy{};
+            std::copy(array.end() - NewN, array.end(), copy.begin());
             return copy;
         }
     }

tests/CMakeLists.txt

@@ -22,6 +22,7 @@ add_executable(tests
         util.cpp
         symmetricTensor.cpp
         multipole.cpp
+        multipoleMoment.cpp
         )
 set_property(TARGET tests PROPERTY CXX_STANDARD 20)
 target_link_libraries(