Skip to content
Snippets Groups Projects
Commit ae74f3e3 authored by Emmanuel Labourasse's avatar Emmanuel Labourasse
Browse files

Add tests for TaylorPolynomial

parent 74b23a60
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
tests/test_TaylorPolynomial.cpp 0 → 100644
+ 331
− 0
View file @ ae74f3e3
#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}));
// }
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment