dd evaluation feature with multi-D Horner method

src/analysis/PolynomialP.hpp

@@ -2,6 +2,7 @@
 #define POLYNOMIALP_HPP
 
 #include <algebra/TinyVector.hpp>
+#include <utils/Exceptions.hpp>
 
 template <size_t N, size_t Dimension>
 class PolynomialP
@@ -141,7 +142,7 @@ class PolynomialP
 
   PUGS_INLINE
   constexpr double
-  coef(const TinyVector<Dimension, size_t> relative_pos) const
+  operator[](const TinyVector<Dimension, size_t> relative_pos) const
   {
     size_t total_degree = 0;
     for (size_t i = 0; i < Dimension; ++i) {
@@ -153,9 +154,63 @@ class PolynomialP
                                 "degree of the polynomial");
     TinyVector<size_coef> absolute_coefs = this->coefficients();
     size_t absolute_position             = 0;
+    if constexpr (Dimension == 1) {
+      absolute_position = relative_pos[0];
+    } else if constexpr (Dimension == 2) {
+      size_t total_degree = relative_pos[0] + relative_pos[1];
+      absolute_position   = total_degree * (total_degree + 1) / 2 + relative_pos[1];
+    } else {
+      throw UnexpectedError("Not yet Available in 3D");
+      // static_assert(Dimension == 3);
+      // size_t total_degree = relative_pos[0] + relative_pos[1] + relative_pos[2];
+      // return (total_degree + 1) * (total_degree + 2) * (total_degree + 3) / 6 + relative_pos[1];
+      // return (N + 1) * (N + 2) * (N + 3) / 6;
+    }
+
     return absolute_coefs[absolute_position];
   }
 
+  PUGS_INLINE
+  constexpr double
+  operator()(const TinyVector<Dimension> x) const
+  {
+    const auto& P = *this;
+    double value  = 0.;
+    if constexpr (Dimension == 1) {
+      value = m_coefficients[N];
+      for (size_t i = N; i > 0; --i) {
+        value *= x;
+        value += m_coefficients[i - 1];
+      }
+    } else if constexpr (Dimension == 2) {
+      TinyVector<Dimension, size_t> relative_pos(0, N);
+      value = P[relative_pos];
+      for (size_t i = N; i > 0; --i) {
+        value *= x[1];
+        relative_pos  = TinyVector<2, size_t>(N - i + 1, i - 1);
+        double valuex = P[relative_pos];
+        for (size_t j = N - i + 1; j > 0; --j) {
+          valuex *= x[0];
+          relative_pos = TinyVector<2, size_t>(j - 1, i - 1);
+          valuex += P[relative_pos];
+        }
+        value += valuex;
+      }
+      // relative_pos  = TinyVector<2, size_t>(N, 0);
+      // double valuex = P[relative_pos];
+      // for (size_t j = N; j > 0; --j) {
+      //   valuex *= x[0];
+      //   relative_pos = TinyVector<2, size_t>(j - 1, 0);
+      //   valuex += P[relative_pos];
+      // }
+      // value += valuex;
+    } else {
+      throw UnexpectedError("Not yet Available in 3D");
+    }
+
+    return value;
+  }
+
   PUGS_INLINE constexpr PolynomialP() noexcept = default;
   ~PolynomialP()                               = default;
 };

tests/test_PolynomialP.cpp

@@ -87,7 +87,39 @@ TEST_CASE("PolynomialP", "[analysis]")
     REQUIRE(Q == (2 * P));
   }
 
-  // SECTION("product")
+  SECTION("access_coef")
+  {
+    TinyVector<6> coef(1, -2, 10, 7, 2, 9);
+    TinyVector<10> coef2(2, -4, -1, -3, 3, -5, -6, 0, 1, 7);
+    PolynomialP<2, 2> P(coef);
+    PolynomialP<3, 2> Q(coef2);
+    TinyVector<2, size_t> relative_coef(1, 1);
+    TinyVector<2, size_t> relative_coef2(1, 2);
+    REQUIRE(P[relative_coef] == 2);
+    REQUIRE(Q[relative_coef2] == 1);
+    REQUIRE(Q[relative_coef] == 3);
+  }
+
+  SECTION("evaluation")
+  {
+    TinyVector<6> coef(1, -2, 10, 7, 2, 9);
+    TinyVector<10> coef2(2, -4, -1, -3, 3, -5, -6, 0, 1, 7);
+    PolynomialP<2, 2> P(coef);
+    PolynomialP<3, 2> Q(coef2);
+    TinyVector<6> coefx(1, -2, 0, 7, 0, 0);
+    TinyVector<6> coefy(2, 0, -2, 0, 0, 7);
+    TinyVector<2> pos(1, -1);
+    TinyVector<2> pos2(-1, 2);
+    PolynomialP<2, 2> Px(coefx);
+    PolynomialP<2, 2> Py(coefy);
+    REQUIRE(Px(pos) == 6);
+    REQUIRE(Py(pos) == 11);
+    REQUIRE(Px(pos2) == 10);
+    REQUIRE(P(pos2) == 62);
+    REQUIRE(Q(pos) == -20);
+  }
+
+  // // SECTION("product")
   // {
   //   Polynomial<2> P(2, 3, 4);
   //   Polynomial<3> Q(1, 2, -1, 1);
@@ -118,12 +150,6 @@ TEST_CASE("PolynomialP", "[analysis]")
   //   REQUIRE(Polynomial<2>{0, 1, 0} == R);
   // }
 
-  // SECTION("evaluation")
-  // {
-  //   Polynomial<2> P(2, -3, 4);
-  //   REQUIRE(P(3) == 29);
-  // }
-
   // SECTION("primitive")
   // {
   //   Polynomial<2> P(2, -3, 4);