diff --git a/tests/test_TaylorPolynomial.cpp b/tests/test_TaylorPolynomial.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..eec530eafeef396d12ff3c651881fa53c30249e8
--- /dev/null
+++ b/tests/test_TaylorPolynomial.cpp
@@ -0,0 +1,331 @@
+#include <Kokkos_Core.hpp>
+#include <algebra/TinyMatrix.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <analysis/SquareGaussQuadrature.hpp>
+#include <analysis/TaylorPolynomial.hpp>
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/Types.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class TaylorPolynomial<0, 2>;
+template class TaylorPolynomial<1, 2>;
+template class TaylorPolynomial<2, 2>;
+template class TaylorPolynomial<3, 2>;
+
+// clazy:excludeall=non-pod-global-static
+TinyVector<2> x0{1, -1};
+TinyVector<3> x1{1, -1, 1};
+
+TEST_CASE("TaylorPolynomial", "[analysis]")
+{
+ SECTION("construction")
+ {
+ TinyVector<6> coef(1, 2, 3, 4, 5, 6);
+ REQUIRE_NOTHROW(TaylorPolynomial<2, 2>(coef, x0));
+ }
+
+ SECTION("degree")
+ {
+ TinyVector<3> coef(1, 2, 3);
+ TaylorPolynomial<1, 2> P(coef, x0);
+ REQUIRE(P.degree() == 1);
+ REQUIRE(P.dim() == 2);
+ }
+
+ SECTION("equality")
+ {
+ TinyVector<6> coef(1, 2, 3, 4, 5, 6);
+ TinyVector<3> coef2(1, 2, 3);
+ TinyVector<6> coef3(1, 2, 3, 3, 3, 3);
+ TaylorPolynomial<2, 2> P(coef, x0);
+ TaylorPolynomial<2, 2> Q(coef, x0);
+ TaylorPolynomial<2, 2> R(coef3, x0);
+
+ REQUIRE(P == Q);
+ REQUIRE(P != R);
+ }
+
+ SECTION("addition")
+ {
+ TinyVector<6> coef(1, 2, 3, 4, 5, 6);
+ TinyVector<6> coef2(1, 2, 3, -2, -1, -3);
+ TinyVector<6> coef3(2, 4, 6, 2, 4, 3);
+ TaylorPolynomial<2, 2> P(coef, x0);
+ TaylorPolynomial<2, 2> Q(coef2, x0);
+ TaylorPolynomial<2, 2> R(coef3, x0);
+ REQUIRE(R == (P + Q));
+ REQUIRE((P + Q) == R);
+ }
+
+ SECTION("opposed")
+ {
+ TinyVector<6> coef(1, 2, 3, 4, 5, 6);
+ TinyVector<6> coef2(-1, -2, -3, -4, -5, -6);
+ TaylorPolynomial<2, 2> P(coef, x0);
+ REQUIRE(-P == TaylorPolynomial<2, 2>(coef2, x0));
+ }
+
+ SECTION("difference")
+ {
+ TinyVector<6> coef(1, 2, 3, 4, 5, 6);
+ TinyVector<6> coef2(1, 2, 3, -2, -1, -3);
+ TinyVector<6> coef3(0, 0, 0, 6, 6, 9);
+ TaylorPolynomial<2, 2> P(coef, x0);
+ TaylorPolynomial<2, 2> Q(coef2, x0);
+ TaylorPolynomial<2, 2> R(coef3, x0);
+ R = P - Q;
+ REQUIRE(R == TaylorPolynomial<2, 2>(coef3, x0));
+ }
+
+ SECTION("product_by_scalar")
+ {
+ TinyVector<6> coef(1, 2, 3, 4, 5, 6);
+ TinyVector<6> coef2(2, 4, 6, 8, 10, 12);
+ TaylorPolynomial<2, 2> P(coef, x0);
+ TaylorPolynomial<2, 2> Q(coef2, x0);
+ REQUIRE(Q == (P * 2));
+ REQUIRE(Q == (2 * P));
+ }
+
+ 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);
+ TinyVector<10> coef3(2, -4, -1, -3, 3, -5, -6, 0, 1, 7);
+ TinyVector<20> coef4(2, -4, -1, -3, 3, -5, -6, -2, 1, 7, 3, -2, 1, 2.5, 6, -9, 0.5, 4, -5, -8);
+ TaylorPolynomial<2, 2> P(coef, x0);
+ TaylorPolynomial<3, 2> Q(coef2, x0);
+ TaylorPolynomial<2, 3> R(coef3, x1);
+ TaylorPolynomial<3, 3> S(coef4, x1);
+ TinyVector<2, size_t> relative_coef(1, 1);
+ TinyVector<2, size_t> relative_coef2(1, 2);
+ TinyVector<3, size_t> relative_coef3(1, 0, 1);
+ TinyVector<3, size_t> relative_coef3b(0, 0, 2);
+ TinyVector<3, size_t> relative_coef4(1, 1, 1);
+ TinyVector<3, size_t> relative_coef5(0, 2, 1);
+ REQUIRE(P[relative_coef] == 2);
+ REQUIRE(Q[relative_coef2] == 1);
+ REQUIRE(Q[relative_coef] == 3);
+ REQUIRE(R[relative_coef3] == -6);
+ REQUIRE(R[relative_coef3b] == 7);
+ REQUIRE(S[relative_coef4] == 6);
+ REQUIRE(S[relative_coef5] == 4);
+ }
+
+ SECTION("evaluation")
+ {
+ TinyVector<6> coef(1, -2, 10, 7, 2, 9);
+ TinyVector<10> coef2(2, -4, -1, -3, 3, -5, -6, 0, 1, 7);
+ TaylorPolynomial<2, 2> P(coef, x0);
+ TaylorPolynomial<3, 2> Q(coef2, x0);
+ 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);
+ TaylorPolynomial<2, 2> Px(coefx, x0);
+ TaylorPolynomial<2, 2> Py(coefy, x0);
+ REQUIRE(Px(pos) == 1);
+ REQUIRE(Py(pos) == 2);
+ REQUIRE(Px(pos2) == 33);
+ REQUIRE(P(pos2) == 132);
+ REQUIRE(Q(pos) == 2);
+ }
+
+ SECTION("derivation")
+ {
+ TinyVector<6> coef(1, -2, 10, 7, 2, 9);
+ TinyVector<6> coef2(-2, 14, 2, 0, 0, 0);
+ TinyVector<6> coef3(10, 2, 18, 0, 0, 0);
+ TinyVector<10> coef4(2, -4, -1, -3, 3, -5, -6, 1, 1, 7);
+ TinyVector<10> coef5(-1, -5, 2, 1, 0, 0, 0, 0, 0, 0);
+ TaylorPolynomial<2, 2> P(coef, x0);
+ TaylorPolynomial<2, 2> Q(coef2, x0);
+ TaylorPolynomial<2, 2> R(coef3, x0);
+ TaylorPolynomial<2, 3> S(coef4, x1);
+ TaylorPolynomial<2, 3> T(coef5, x1);
+
+ REQUIRE(Q == P.derivative(0));
+ REQUIRE(R == P.derivative(1));
+ REQUIRE(T == S.derivative(1));
+ }
+
+ SECTION("integrate")
+ {
+ TinyVector<6> coef(1, -2, 10, 7, 2, 9);
+ std::array<TinyVector<2>, 4> positions;
+ std::array<TinyVector<2>, 3> positions2;
+ std::array<TinyVector<2>, 4> positions3;
+ std::array<TinyVector<2>, 2> positions4;
+ positions[0] = TinyVector<2>{0, 0};
+ positions[1] = TinyVector<2>{0, 0.5};
+ positions[2] = TinyVector<2>{0.3, 0.7};
+ positions[3] = TinyVector<2>{0.4, 0.1};
+ positions2[0] = TinyVector<2>{0, 0};
+ positions2[1] = TinyVector<2>{0, 0.5};
+ positions2[2] = TinyVector<2>{0.3, 0.7};
+ positions3[0] = TinyVector<2>{0, 0};
+ positions3[1] = TinyVector<2>{1, 0};
+ positions3[2] = TinyVector<2>{1, 1};
+ positions3[3] = TinyVector<2>{0, 1};
+ positions4[0] = TinyVector<2>{0, 0.5};
+ positions4[1] = TinyVector<2>{0.3, 0.7};
+ TaylorPolynomial<2, 2> P(coef, x0);
+ auto p1 = [](const TinyVector<2>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ return 1 - 2. * (x - 1.) + 10 * (y + 1) + 7 * (x - 1.) * (x - 1) + 2 * (x - 1) * (y + 1) + 9 * (y + 1) * (y + 1);
+ };
+ const QuadratureFormula<2>& l3 = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(3));
+ auto point_list = l3.pointList();
+ auto weight_list = l3.weightList();
+ SquareTransformation<2> s{positions[0], positions[1], positions[2], positions[3]};
+ auto value = weight_list[0] * s.jacobianDeterminant(point_list[0]) * p1(s(point_list[0]));
+ for (size_t i = 1; i < weight_list.size(); ++i) {
+ value += weight_list[i] * s.jacobianDeterminant(point_list[i]) * p1(s(point_list[i]));
+ }
+ SquareTransformation<2> s3{positions3[0], positions3[1], positions3[2], positions3[3]};
+ auto value3 = weight_list[0] * s3.jacobianDeterminant(point_list[0]) * p1(s3(point_list[0]));
+ for (size_t i = 1; i < weight_list.size(); ++i) {
+ value3 += weight_list[i] * s3.jacobianDeterminant(point_list[i]) * p1(s3(point_list[i]));
+ }
+ const QuadratureFormula<2>& t3 = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(3));
+ auto point_list2 = t3.pointList();
+ auto weight_list2 = t3.weightList();
+ TriangleTransformation<2> t{positions2[0], positions2[1], positions2[2]};
+ auto value2 = weight_list2[0] * p1(t(point_list2[0]));
+ for (size_t i = 1; i < weight_list2.size(); ++i) {
+ value2 += weight_list2[i] * p1(t(point_list2[i]));
+ }
+ const QuadratureFormula<1>& l1 = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2));
+ const LineTransformation<2> u{positions4[0], positions4[1]};
+ double value4 = 0.;
+ for (size_t i = 0; i < l1.numberOfPoints(); ++i) {
+ value4 += l1.weight(i) * u.velocityNorm() * p1(u(l1.point(i)));
+ }
+
+ REQUIRE(value == Catch::Approx(integrate(P, positions)));
+ REQUIRE(value2 == Catch::Approx(integrate(P, positions2)));
+ REQUIRE(value3 == Catch::Approx(integrate(P, positions3)));
+ REQUIRE(value4 == Catch::Approx(integrate(P, positions4)));
+ }
+
+ // // SECTION("product")
+ // {
+ // Polynomial<2> P(2, 3, 4);
+ // Polynomial<3> Q(1, 2, -1, 1);
+ // Polynomial<4 > R;
+ // Polynomial<5> S;
+ // R = P;
+ // S = P;
+ // S *= Q;
+ // REQUIRE(Polynomial<5>(2, 7, 8, 7, -1, 4) == (P * Q));
+ // REQUIRE(Polynomial<5>(2, 7, 8, 7, -1, 4) == S);
+ // // REQUIRE_THROWS_AS(R *= Q, AssertError);
+ // }
+
+ // SECTION("divide")
+ // {
+ // Polynomial<2> P(1, 0, 1);
+ // Polynomial<1> Q(0, 1);
+ // Polynomial<1> Q1(0, 1);
+
+ // Polynomial<2> R;
+ // Polynomial<2> S;
+ // REQUIRE(P.realDegree() == 2);
+ // REQUIRE(Q.realDegree() == 1);
+ // REQUIRE(Q1.realDegree() == 1);
+
+ // divide(P, Q1, R, S);
+ // REQUIRE(Polynomial<2>{1, 0, 0} == S);
+ // REQUIRE(Polynomial<2>{0, 1, 0} == R);
+ // }
+
+ // SECTION("primitive")
+ // {
+ // Polynomial<2> P(2, -3, 4);
+ // TinyVector<4> coefs = zero;
+ // Polynomial<3> Q(coefs);
+ // Q = primitive(P);
+ // Polynomial<3> R(0, 2, -3. / 2, 4. / 3);
+ // REQUIRE(Q == R);
+ // }
+
+ // SECTION("integrate")
+ // {
+ // Polynomial<2> P(2, -3, 3);
+ // double xinf = -1;
+ // double xsup = 1;
+ // double result = integrate(P, xinf, xsup);
+ // REQUIRE(result == 6);
+ // result = symmetricIntegrate(P, 2);
+ // REQUIRE(result == 24);
+ // }
+
+ // SECTION("derivative")
+ // {
+ // Polynomial<2> P(2, -3, 3);
+ // Polynomial<1> Q = derivative(P);
+ // REQUIRE(Q == Polynomial<1>(-3, 6));
+
+ // Polynomial<0> P2(3);
+
+ // Polynomial<0> R(0);
+ // REQUIRE(derivative(P2) == R);
+ // }
+
+ // SECTION("affectation")
+ // {
+ // Polynomial<2> Q(2, -3, 3);
+ // Polynomial<4> R(2, -3, 3, 0, 0);
+ // Polynomial<4> P(0, 1, 2, 3, 3);
+ // P = Q;
+ // REQUIRE(P == R);
+ // }
+
+ // SECTION("affectation addition")
+ // {
+ // Polynomial<2> Q(2, -3, 3);
+ // Polynomial<4> R(2, -2, 5, 3, 3);
+ // Polynomial<4> P(0, 1, 2, 3, 3);
+ // P += Q;
+ // REQUIRE(P == R);
+ // }
+
+ // SECTION("power")
+ // {
+ // Polynomial<2> P(2, -3, 3);
+ // Polynomial<4> R(4, -12, 21, -18, 9);
+ // Polynomial<1> Q(0, 2);
+ // Polynomial<2> S = Q.pow<2>(2);
+ // REQUIRE(P.pow<2>(2) == R);
+ // REQUIRE(S == Polynomial<2>(0, 0, 4));
+ // }
+
+ // SECTION("composition")
+ // {
+ // Polynomial<2> P(2, -3, 3);
+ // Polynomial<1> Q(0, 2);
+ // Polynomial<2> R(2, -1, 3);
+ // Polynomial<2> S(1, 2, 2);
+ // REQUIRE(P.compose(Q) == Polynomial<2>(2, -6, 12));
+ // REQUIRE(P.compose2(Q) == Polynomial<2>(2, -6, 12));
+ // REQUIRE(R(S) == Polynomial<4>(4, 10, 22, 24, 12));
+ // }
+
+ // SECTION("Lagrange polynomial")
+ // {
+ // Polynomial<1> S(0.5, -0.5);
+ // Polynomial<1> Q;
+ // Q = lagrangePolynomial<1>(TinyVector<2>{-1, 1}, 0);
+ // REQUIRE(S == Q);
+ // Polynomial<2> P(0, -0.5, 0.5);
+ // Polynomial<2> R;
+ // R = lagrangePolynomial<2>(TinyVector<3>{-1, 0, 1}, 0);
+ // REQUIRE(R == P);
+ // const std::array<Polynomial<2>, 3> basis = lagrangeBasis(TinyVector<3>{-1, 0, 1});
+ // REQUIRE(lagrangeToCanonical(TinyVector<3>{1, 0, 1}, basis) == Polynomial<2>(TinyVector<3>{0, 0, 1}));
+ // }
+}