Normal error fixed

src/analysis/PolynomialP.hpp

@@ -109,6 +109,7 @@ class PolynomialP
     }
     return *this;
   }
+
   PUGS_INLINE constexpr PolynomialP&
   operator+=(const PolynomialP& Q)
   {
@@ -170,6 +171,65 @@ class PolynomialP
     return absolute_coefs[absolute_position];
   }
 
+  PUGS_INLINE
+  constexpr double
+  operator[](const TinyVector<Dimension, size_t> relative_pos)
+  {
+    size_t total_degree = 0;
+    for (size_t i = 0; i < Dimension; ++i) {
+      Assert((relative_pos[i] <= N),
+             "You are looking for a coefficient of order greater than the degree of the polynomial");
+      total_degree += relative_pos[i];
+    }
+    Assert((total_degree <= N), "The sum of the degrees of the coefficient you are looking for is greater than the "
+                                "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
+  absolute_position(const TinyVector<Dimension, size_t> relative_pos)
+  {
+    size_t total_degree = 0;
+    for (size_t i = 0; i < Dimension; ++i) {
+      Assert((relative_pos[i] <= N),
+             "You are looking for a coefficient of order greater than the degree of the polynomial");
+      total_degree += relative_pos[i];
+    }
+    Assert((total_degree <= N), "The sum of the degrees of the coefficient you are looking for is greater than the "
+                                "degree of the polynomial");
+    size_t abs_pos = 0;
+    if constexpr (Dimension == 1) {
+      abs_pos = relative_pos[0];
+    } else if constexpr (Dimension == 2) {
+      size_t total_degree = relative_pos[0] + relative_pos[1];
+      abs_pos             = 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 abs_pos;
+  }
+
   PUGS_INLINE
   constexpr double
   operator()(const TinyVector<Dimension> x) const
@@ -187,23 +247,15 @@ class PolynomialP
       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);
+        relative_pos  = TinyVector<Dimension, 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);
+          relative_pos = TinyVector<Dimension, 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");
     }
@@ -211,23 +263,62 @@ class PolynomialP
     return value;
   }
 
-  PUGS_INLINE constexpr PolynomialP() noexcept = default;
-  ~PolynomialP()                               = default;
-};
+  PUGS_INLINE size_t
+  find_size_coef(const size_t degree)
+  {
+    if constexpr (Dimension == 1) {
+      return degree + 1;
+    } else if constexpr (Dimension == 2) {
+      return (degree + 1) * (degree + 2) / 2;
+    } else {
+      static_assert(Dimension == 3);
+      return (degree + 1) * (degree + 2) * (degree + 3) / 6;
+    }
+  }
 
-//   // evaluation using Horner's method https://en.wikipedia.org/wiki/Horner's_method
-//   PUGS_INLINE
-//   constexpr double
-//   evaluate(const double& x) const
+  // PUGS_INLINE constexpr auto
+  // derivative(const PolynomialP& P, const size_t var)
   // {
-//     TinyVector<N + 1> coefs = this->coefficients();
-//     double bcoef            = coefs[N];
-//     for (size_t i = N; i > 0; --i) {
-//       bcoef *= x;
-//       bcoef += coefs[i - 1];
+  //   TinyVector<size_coef> coefs(zero);
+  //   PolynomialP<N, Dimension> Q(coefs);
+  //   if constexpr (N == 0) {
+  //     return Q;
+  //   } else {
+  //     Assert(var < Dimension,
+  //            "You can not derive a polynomial with respect to a variable of rank greater than the dimension");
+  //     if constexpr (Dimension == 1) {
+  //       for (size_t i = 0; i < size_coef; ++i) {
+  //         coefs[i] = double(i + 1) * P.coefficients()[i + 1];
+  //       }
+  //       return Q;
+  //     } else if constexpr (Dimension == 2) {
+  //       if (var == 0) {
+  //         for (size_t i = 0; i < N; ++i) {
+  //           for (size_t j = 0; j < N; ++i) {
+  //             TinyVector<Dimension, size_t> relative_pos(i, j);
+  //             TinyVector<Dimension, size_t> relative_posp(i + 1, j);
+  //             Q[relative_pos] = double(i + 1) * P[relative_posp];
   //           }
-//     return bcoef;
   //         }
+  //       } else {
+  //         for (size_t i = 0; i < N; ++i) {
+  //           for (size_t j = 0; j < N; ++i) {
+  //             TinyVector<Dimension, size_t> relative_pos(i, j);
+  //             TinyVector<Dimension, size_t> relative_posp(i, j + 1);
+  //             Q[relative_pos] = double(j + 1) * P[relative_posp];
+  //           }
+  //         }
+  //       }
+  //       return Q;
+  //     } else {
+  //       throw UnexpectedError("Not yet Available in 3D");
+  //     }
+  //   }
+  // }
+
+  PUGS_INLINE constexpr PolynomialP() noexcept = default;
+  ~PolynomialP()                               = default;
+};
 
 //   template <size_t M>
 //   PUGS_INLINE constexpr friend void
@@ -378,21 +469,6 @@ class PolynomialP
 //   return integral;
 // }
 
-// template <size_t N>
-// PUGS_INLINE constexpr auto
-// derivative(const PolynomialP<N>& P)
-// {
-//   if constexpr (N == 0) {
-//     return PolynomialP<0>(0);
-//   } else {
-//     TinyVector<N> coefs;
-//     for (size_t i = 0; i < N; ++i) {
-//       coefs[i] = double(i + 1) * P.coefficients()[i + 1];
-//     }
-//     return PolynomialP<N - 1>(coefs);
-//   }
-// }
-
 // template <size_t N>
 // PUGS_INLINE constexpr PolynomialP<N>
 // lagrangePolynomialP(const TinyVector<N + 1> zeros, const size_t k)

tests/test_PolynomialP.cpp

@@ -116,7 +116,7 @@ TEST_CASE("PolynomialP", "[analysis]")
     REQUIRE(Py(pos) == 11);
     REQUIRE(Px(pos2) == 10);
     REQUIRE(P(pos2) == 62);
-    REQUIRE(Q(pos) == -20);
+    REQUIRE(Q(pos) == -24);
   }
 
   // // SECTION("product")