From 7fdd2b62f9566a1e484b2fd59ba8cccc43bc8647 Mon Sep 17 00:00:00 2001 From: Pasukhin Dmitry Date: Sun, 7 Sep 2025 19:17:13 +0100 Subject: [PATCH] Testing - Add unit tests for PLib functionality (#705) - Introduced comprehensive unit tests for the Jacobi polynomial implementation in PLib, covering constructors, edge cases, Gauss integration points, weights, and basis function evaluations. - Added tests for basic utility functions in PLib, including pole conversion and binomial coefficient calculations. - Implemented checks for Hermite interpolation and polynomial evaluation with derivatives. - Enhanced error handling and edge case testing for small and large coefficients. - Initialised MaxError in PLib_DoubleJacobiPolynomial to ensure consistent behaviour during degree reduction. --- .../TKMath/GTests/FILES.cmake | 4 + .../PLib_DoubleJacobiPolynomial_Test.cxx | 402 +++++++++++++++++ .../TKMath/GTests/PLib_HermitJacobi_Test.cxx | 310 +++++++++++++ .../GTests/PLib_JacobiPolynomial_Test.cxx | 366 +++++++++++++++ .../TKMath/GTests/PLib_Test.cxx | 418 ++++++++++++++++++ .../PLib/PLib_DoubleJacobiPolynomial.cxx | 2 + 6 files changed, 1502 insertions(+) create mode 100644 src/FoundationClasses/TKMath/GTests/PLib_DoubleJacobiPolynomial_Test.cxx create mode 100644 src/FoundationClasses/TKMath/GTests/PLib_HermitJacobi_Test.cxx create mode 100644 src/FoundationClasses/TKMath/GTests/PLib_JacobiPolynomial_Test.cxx create mode 100644 src/FoundationClasses/TKMath/GTests/PLib_Test.cxx diff --git a/src/FoundationClasses/TKMath/GTests/FILES.cmake b/src/FoundationClasses/TKMath/GTests/FILES.cmake index 4620e54c43..ea90fb0e64 100644 --- a/src/FoundationClasses/TKMath/GTests/FILES.cmake +++ b/src/FoundationClasses/TKMath/GTests/FILES.cmake @@ -34,4 +34,8 @@ set(OCCT_TKMath_GTests_FILES math_TrigonometricFunctionRoots_Test.cxx math_Uzawa_Test.cxx math_Vector_Test.cxx + PLib_Test.cxx + PLib_JacobiPolynomial_Test.cxx + PLib_HermitJacobi_Test.cxx + PLib_DoubleJacobiPolynomial_Test.cxx ) diff --git a/src/FoundationClasses/TKMath/GTests/PLib_DoubleJacobiPolynomial_Test.cxx b/src/FoundationClasses/TKMath/GTests/PLib_DoubleJacobiPolynomial_Test.cxx new file mode 100644 index 0000000000..64ba4e0833 --- /dev/null +++ b/src/FoundationClasses/TKMath/GTests/PLib_DoubleJacobiPolynomial_Test.cxx @@ -0,0 +1,402 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +// Test fixture for PLib_DoubleJacobiPolynomial tests +class PLibDoubleJacobiPolynomialTest : public ::testing::Test +{ +protected: + void SetUp() override + { + // Create Jacobi polynomials for testing + myJacobiU = new PLib_JacobiPolynomial(10, GeomAbs_C0); + myJacobiV = new PLib_JacobiPolynomial(8, GeomAbs_C1); + } + + void TearDown() override + { + myJacobiU.Nullify(); + myJacobiV.Nullify(); + } + + Handle(PLib_JacobiPolynomial) myJacobiU; + Handle(PLib_JacobiPolynomial) myJacobiV; + + // Helper to create test coefficients + void createTestCoefficients(TColStd_Array1OfReal& theCoeffs, + Standard_Integer theDimension, + Standard_Integer theDegreeU, + Standard_Integer theDegreeV) + { + Standard_Integer anIndex = theCoeffs.Lower(); // Start from array lower bound + for (Standard_Integer j = 0; j <= theDegreeV; j++) + { + for (Standard_Integer i = 0; i <= theDegreeU; i++) + { + for (Standard_Integer d = 1; d <= theDimension; d++) + { + theCoeffs(anIndex) = Sin(anIndex * 0.1 + i * 0.2 + j * 0.3); + anIndex++; + } + } + } + } +}; + +// Test default constructor +TEST_F(PLibDoubleJacobiPolynomialTest, DefaultConstructor) +{ + PLib_DoubleJacobiPolynomial aDoubleJac; + + // After default construction, accessors should return null handles + EXPECT_TRUE(aDoubleJac.U().IsNull()) << "U polynomial should be null after default construction"; + EXPECT_TRUE(aDoubleJac.V().IsNull()) << "V polynomial should be null after default construction"; + EXPECT_TRUE(aDoubleJac.TabMaxU().IsNull()) << "TabMaxU should be null after default construction"; + EXPECT_TRUE(aDoubleJac.TabMaxV().IsNull()) << "TabMaxV should be null after default construction"; +} + +// Test parameterized constructor +TEST_F(PLibDoubleJacobiPolynomialTest, ParameterizedConstructor) +{ + ASSERT_FALSE(myJacobiU.IsNull()) << "JacobiU should not be null"; + ASSERT_FALSE(myJacobiV.IsNull()) << "JacobiV should not be null"; + + PLib_DoubleJacobiPolynomial aDoubleJac(myJacobiU, myJacobiV); + + // After construction with parameters, accessors should return valid handles + EXPECT_FALSE(aDoubleJac.U().IsNull()) << "U polynomial should not be null"; + EXPECT_FALSE(aDoubleJac.V().IsNull()) << "V polynomial should not be null"; + EXPECT_FALSE(aDoubleJac.TabMaxU().IsNull()) << "TabMaxU should not be null"; + EXPECT_FALSE(aDoubleJac.TabMaxV().IsNull()) << "TabMaxV should not be null"; + + // Test that the returned handles are correct + EXPECT_EQ(aDoubleJac.U().get(), myJacobiU.get()) << "U polynomial handle should match"; + EXPECT_EQ(aDoubleJac.V().get(), myJacobiV.get()) << "V polynomial handle should match"; +} + +// Test MaxErrorU calculation +TEST_F(PLibDoubleJacobiPolynomialTest, MaxErrorU) +{ + PLib_DoubleJacobiPolynomial aDoubleJac(myJacobiU, myJacobiV); + + const Standard_Integer aDimension = 2; + const Standard_Integer aDegreeU = 6; + const Standard_Integer aDegreeV = 5; + const Standard_Integer aDJacCoeff = (aDegreeU + 1) * aDimension; + + // JacCoeff array must be sized based on WorkDegrees + const Standard_Integer aWorkDegreeU = myJacobiU->WorkDegree(); + const Standard_Integer aWorkDegreeV = myJacobiV->WorkDegree(); + const Standard_Integer aCoeffCount = (aWorkDegreeU + 1) * (aWorkDegreeV + 1) * aDimension; + + TColStd_Array1OfReal aJacCoeff(0, aCoeffCount - 1); + for (Standard_Integer i = aJacCoeff.Lower(); i <= aJacCoeff.Upper(); i++) + { + aJacCoeff(i) = Sin(i * 0.1); + } + + EXPECT_NO_THROW({ + Standard_Real aMaxErrU = + aDoubleJac.MaxErrorU(aDimension, aDegreeU, aDegreeV, aDJacCoeff, aJacCoeff); + EXPECT_GE(aMaxErrU, 0.0) << "MaxErrorU should be non-negative"; + EXPECT_FALSE(Precision::IsInfinite(aMaxErrU)) << "MaxErrorU should be finite"; + }) << "MaxErrorU calculation failed"; +} + +// Test MaxErrorV calculation +TEST_F(PLibDoubleJacobiPolynomialTest, MaxErrorV) +{ + PLib_DoubleJacobiPolynomial aDoubleJac(myJacobiU, myJacobiV); + + const Standard_Integer aDimension = 2; + const Standard_Integer aDegreeU = 6; + const Standard_Integer aDegreeV = 5; + const Standard_Integer aDJacCoeff = (aDegreeU + 1) * aDimension; + + // JacCoeff array must be sized based on WorkDegrees + const Standard_Integer aWorkDegreeU = myJacobiU->WorkDegree(); + const Standard_Integer aWorkDegreeV = myJacobiV->WorkDegree(); + const Standard_Integer aCoeffCount = (aWorkDegreeU + 1) * (aWorkDegreeV + 1) * aDimension; + + TColStd_Array1OfReal aJacCoeff(0, aCoeffCount - 1); + for (Standard_Integer i = aJacCoeff.Lower(); i <= aJacCoeff.Upper(); i++) + { + aJacCoeff(i) = Sin(i * 0.1); + } + + EXPECT_NO_THROW({ + Standard_Real aMaxErrV = + aDoubleJac.MaxErrorV(aDimension, aDegreeU, aDegreeV, aDJacCoeff, aJacCoeff); + EXPECT_GE(aMaxErrV, 0.0) << "MaxErrorV should be non-negative"; + EXPECT_FALSE(Precision::IsInfinite(aMaxErrV)) << "MaxErrorV should be finite"; + }) << "MaxErrorV calculation failed"; +} + +// Test general MaxError calculation +TEST_F(PLibDoubleJacobiPolynomialTest, MaxError) +{ + PLib_DoubleJacobiPolynomial aDoubleJac(myJacobiU, myJacobiV); + + const Standard_Integer aDimension = 1; + const Standard_Integer aMinDegreeU = 2, aMaxDegreeU = 5; + const Standard_Integer aMinDegreeV = 4, aMaxDegreeV = 7; // MinDegreeV must be >= MinV=4 + const Standard_Integer aDJacCoeff = (aMaxDegreeU + 1) * aDimension; + const Standard_Real anError = 1e-6; + + // JacCoeff array must be sized based on WorkDegrees + const Standard_Integer aWorkDegreeU = myJacobiU->WorkDegree(); + const Standard_Integer aWorkDegreeV = myJacobiV->WorkDegree(); + const Standard_Integer aCoeffCount = (aWorkDegreeU + 1) * (aWorkDegreeV + 1) * aDimension; + + TColStd_Array1OfReal aJacCoeff(0, aCoeffCount - 1); + for (Standard_Integer i = aJacCoeff.Lower(); i <= aJacCoeff.Upper(); i++) + { + aJacCoeff(i) = Sin(i * 0.1); + } + + EXPECT_NO_THROW({ + Standard_Real aMaxErr = aDoubleJac.MaxError(aDimension, + aMinDegreeU, + aMaxDegreeU, + aMinDegreeV, + aMaxDegreeV, + aDJacCoeff, + aJacCoeff, + anError); + EXPECT_GE(aMaxErr, 0.0) << "MaxError should be non-negative"; + EXPECT_FALSE(Precision::IsInfinite(aMaxErr)) << "MaxError should be finite"; + }) << "MaxError calculation failed"; +} + +// Test degree reduction +TEST_F(PLibDoubleJacobiPolynomialTest, ReduceDegree) +{ + PLib_DoubleJacobiPolynomial aDoubleJac(myJacobiU, myJacobiV); + + const Standard_Integer aDimension = 1; + const Standard_Integer aMinDegreeU = 2, aMaxDegreeU = 6; // Must be >= MinU=2 + const Standard_Integer aMinDegreeV = 4, aMaxDegreeV = 7; // Must be >= MinV=4 + const Standard_Integer aDJacCoeff = (aMaxDegreeU + 1) * aDimension; + const Standard_Real anEpmsCut = 1e-8; + + // JacCoeff array must be sized based on WorkDegrees and use 0-based indexing + const Standard_Integer aWorkDegreeU = myJacobiU->WorkDegree(); + const Standard_Integer aWorkDegreeV = myJacobiV->WorkDegree(); + const Standard_Integer aCoeffCount = (aWorkDegreeU + 1) * (aWorkDegreeV + 1) * aDimension; + + TColStd_Array1OfReal aJacCoeff(0, aCoeffCount - 1); + for (Standard_Integer i = aJacCoeff.Lower(); i <= aJacCoeff.Upper(); i++) + { + aJacCoeff(i) = 1.0 / (i + 10); // Small decreasing values + } + + Standard_Real aMaxError = -1.0; + Standard_Integer aNewDegreeU = -1, aNewDegreeV = -1; + + EXPECT_NO_THROW({ + aDoubleJac.ReduceDegree(aDimension, + aMinDegreeU, + aMaxDegreeU, + aMinDegreeV, + aMaxDegreeV, + aDJacCoeff, + aJacCoeff, + anEpmsCut, + aMaxError, + aNewDegreeU, + aNewDegreeV); + }) << "ReduceDegree failed"; + + // Verify results are reasonable + EXPECT_LE(aNewDegreeU, aMaxDegreeU) << "New U degree should not exceed max"; + EXPECT_GE(aNewDegreeU, aMinDegreeU) << "New U degree should be at least min"; + EXPECT_LE(aNewDegreeV, aMaxDegreeV) << "New V degree should not exceed max"; + EXPECT_GE(aNewDegreeV, aMinDegreeV) << "New V degree should be at least min"; + EXPECT_GE(aMaxError, 0.0) << "Max error should be non-negative"; + EXPECT_FALSE(Precision::IsInfinite(aMaxError)) << "Max error should be finite"; +} + +// Test average error calculation +TEST_F(PLibDoubleJacobiPolynomialTest, AverageError) +{ + PLib_DoubleJacobiPolynomial aDoubleJac(myJacobiU, myJacobiV); + + const Standard_Integer aDimension = 2; + const Standard_Integer aDegreeU = 4; + const Standard_Integer aDegreeV = 3; + const Standard_Integer aDJacCoeff = 0; // For 0-based arrays, start offset should be 0 + + // JacCoeff array must be sized based on WorkDegrees and use 0-based indexing + const Standard_Integer aWorkDegreeU = myJacobiU->WorkDegree(); + const Standard_Integer aWorkDegreeV = myJacobiV->WorkDegree(); + const Standard_Integer aCoeffCount = (aWorkDegreeU + 1) * (aWorkDegreeV + 1) * aDimension; + + TColStd_Array1OfReal aJacCoeff(0, aCoeffCount - 1); + for (Standard_Integer i = aJacCoeff.Lower(); i <= aJacCoeff.Upper(); i++) + { + aJacCoeff(i) = Sin(i * 0.1); + } + + EXPECT_NO_THROW({ + Standard_Real aAvgErr = + aDoubleJac.AverageError(aDimension, aDegreeU, aDegreeV, aDJacCoeff, aJacCoeff); + EXPECT_GE(aAvgErr, 0.0) << "Average error should be non-negative"; + EXPECT_FALSE(Precision::IsInfinite(aAvgErr)) << "Average error should be finite"; + }) << "AverageError calculation failed"; +} + +// Test coefficient conversion to canonical base +TEST_F(PLibDoubleJacobiPolynomialTest, CoefficientConversion) +{ + PLib_DoubleJacobiPolynomial aDoubleJac(myJacobiU, myJacobiV); + + const Standard_Integer aDimension = 3; + const Standard_Integer aDegreeU = 3; + const Standard_Integer aDegreeV = 2; + + // JacCoeff array must be sized based on WorkDegrees, not input degrees + const Standard_Integer aWorkDegreeU = myJacobiU->WorkDegree(); + const Standard_Integer aWorkDegreeV = myJacobiV->WorkDegree(); + const Standard_Integer aJacCoeffSize = (aWorkDegreeU + 1) * (aWorkDegreeV + 1) * aDimension; + + // JacCoeff uses 0-based indexing as per the implementation + TColStd_Array1OfReal aJacCoeff(0, aJacCoeffSize - 1); + // Initialize with test data - use simple pattern for valid coefficients + for (Standard_Integer i = aJacCoeff.Lower(); i <= aJacCoeff.Upper(); i++) + { + aJacCoeff(i) = Sin(i * 0.1); + } + + // Output coefficients array is sized based on actual degrees + const Standard_Integer aCoeffCount = (aDegreeU + 1) * (aDegreeV + 1) * aDimension; + TColStd_Array1OfReal aCoefficients(0, aCoeffCount - 1); + + EXPECT_NO_THROW({ + aDoubleJac.WDoubleJacobiToCoefficients(aDimension, + aDegreeU, + aDegreeV, + aJacCoeff, + aCoefficients); + }) << "Coefficient conversion failed"; + + // Verify output coefficients are finite + for (Standard_Integer i = aCoefficients.Lower(); i <= aCoefficients.Upper(); i++) + { + EXPECT_FALSE(Precision::IsInfinite(aCoefficients(i))) + << "Converted coefficient should be finite at index " << i; + } +} + +// Test with mismatched degrees +TEST_F(PLibDoubleJacobiPolynomialTest, MismatchedDegrees) +{ + Handle(PLib_JacobiPolynomial) aJacU_Low = new PLib_JacobiPolynomial(5, GeomAbs_C0); + Handle(PLib_JacobiPolynomial) aJacV_High = new PLib_JacobiPolynomial(20, GeomAbs_C1); + + EXPECT_NO_THROW({ + PLib_DoubleJacobiPolynomial aDoubleJac(aJacU_Low, aJacV_High); + + // Should handle mismatched degrees gracefully + EXPECT_FALSE(aDoubleJac.U().IsNull()); + EXPECT_FALSE(aDoubleJac.V().IsNull()); + }) << "Constructor should handle mismatched degrees"; +} + +// Test accessor consistency +TEST_F(PLibDoubleJacobiPolynomialTest, AccessorConsistency) +{ + PLib_DoubleJacobiPolynomial aDoubleJac(myJacobiU, myJacobiV); + + // Test that accessors return consistent results + Handle(PLib_JacobiPolynomial) aReturnedU = aDoubleJac.U(); + Handle(PLib_JacobiPolynomial) aReturnedV = aDoubleJac.V(); + Handle(TColStd_HArray1OfReal) aTabMaxU = aDoubleJac.TabMaxU(); + Handle(TColStd_HArray1OfReal) aTabMaxV = aDoubleJac.TabMaxV(); + + // Multiple calls should return same handles + EXPECT_EQ(aReturnedU.get(), aDoubleJac.U().get()) << "U accessor should be consistent"; + EXPECT_EQ(aReturnedV.get(), aDoubleJac.V().get()) << "V accessor should be consistent"; + EXPECT_EQ(aTabMaxU.get(), aDoubleJac.TabMaxU().get()) << "TabMaxU accessor should be consistent"; + EXPECT_EQ(aTabMaxV.get(), aDoubleJac.TabMaxV().get()) << "TabMaxV accessor should be consistent"; + + // Verify TabMax arrays are properly sized and contain valid data + if (!aTabMaxU.IsNull()) + { + const TColStd_Array1OfReal& anArrU = aTabMaxU->Array1(); + for (Standard_Integer i = anArrU.Lower(); i <= anArrU.Upper(); i++) + { + EXPECT_GT(anArrU(i), 0.0) << "TabMaxU values should be positive"; + EXPECT_FALSE(Precision::IsInfinite(anArrU(i))) << "TabMaxU values should be finite"; + } + } + + if (!aTabMaxV.IsNull()) + { + const TColStd_Array1OfReal& anArrV = aTabMaxV->Array1(); + for (Standard_Integer i = anArrV.Lower(); i <= anArrV.Upper(); i++) + { + EXPECT_GT(anArrV(i), 0.0) << "TabMaxV values should be positive"; + EXPECT_FALSE(Precision::IsInfinite(anArrV(i))) << "TabMaxV values should be finite"; + } + } +} + +// Stress test with large degrees +TEST_F(PLibDoubleJacobiPolynomialTest, StressTest) +{ + Handle(PLib_JacobiPolynomial) aJacU_Large = new PLib_JacobiPolynomial(25, GeomAbs_C1); + Handle(PLib_JacobiPolynomial) aJacV_Large = new PLib_JacobiPolynomial(20, GeomAbs_C2); + + EXPECT_NO_THROW({ + PLib_DoubleJacobiPolynomial aDoubleJac(aJacU_Large, aJacV_Large); + + // Test basic operations don't crash with large degrees + const Standard_Integer aDimension = 1; + const Standard_Integer aDegreeU = 15; + const Standard_Integer aDegreeV = 12; + + // Array must be sized based on WorkDegrees for proper method operation + const Standard_Integer aWorkDegreeU = aJacU_Large->WorkDegree(); + const Standard_Integer aWorkDegreeV = aJacV_Large->WorkDegree(); + const Standard_Integer aCoeffCount = (aWorkDegreeU + 1) * (aWorkDegreeV + 1) * aDimension; + const Standard_Integer aDJacCoeff = 0; // For 0-based arrays + + TColStd_Array1OfReal aJacCoeff(0, aCoeffCount - 1); + for (Standard_Integer i = aJacCoeff.Lower(); i <= aJacCoeff.Upper(); i++) + { + aJacCoeff(i) = Sin(i * 0.1); + } + + Standard_Real aMaxErrU = + aDoubleJac.MaxErrorU(aDimension, aDegreeU, aDegreeV, aDJacCoeff, aJacCoeff); + Standard_Real aMaxErrV = + aDoubleJac.MaxErrorV(aDimension, aDegreeU, aDegreeV, aDJacCoeff, aJacCoeff); + Standard_Real aAvgErr = + aDoubleJac.AverageError(aDimension, aDegreeU, aDegreeV, aDJacCoeff, aJacCoeff); + + EXPECT_GE(aMaxErrU, 0.0) << "MaxErrorU should be non-negative"; + EXPECT_GE(aMaxErrV, 0.0) << "MaxErrorV should be non-negative"; + EXPECT_GE(aAvgErr, 0.0) << "AverageError should be non-negative"; + }) << "Large degree operations should complete without crashing"; +} \ No newline at end of file diff --git a/src/FoundationClasses/TKMath/GTests/PLib_HermitJacobi_Test.cxx b/src/FoundationClasses/TKMath/GTests/PLib_HermitJacobi_Test.cxx new file mode 100644 index 0000000000..a2b4f68f9b --- /dev/null +++ b/src/FoundationClasses/TKMath/GTests/PLib_HermitJacobi_Test.cxx @@ -0,0 +1,310 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include +#include + +#include + +#include +#include +#include +#include +#include + +// Test fixture for PLib_HermitJacobi tests +class PLibHermitJacobiTest : public ::testing::Test +{ +protected: + void SetUp() override + { + // Common test setup + } + + void TearDown() override + { + // Cleanup + } + + // Helper to create valid HermitJacobi polynomial instances + Handle(PLib_HermitJacobi) createHermitJacobi(Standard_Integer theDegree, + GeomAbs_Shape theConstraint) + { + // Ensure degree is within valid range + EXPECT_LE(theDegree, 30) << "Degree too high for HermitJacobi polynomial"; + EXPECT_GE(theDegree, 0) << "Degree must be non-negative"; + + return new PLib_HermitJacobi(theDegree, theConstraint); + } +}; + +// Test constructor and basic properties +TEST_F(PLibHermitJacobiTest, ConstructorAndBasicProperties) +{ + // Test with different constraint orders + Handle(PLib_HermitJacobi) aHermC0 = createHermitJacobi(10, GeomAbs_C0); + Handle(PLib_HermitJacobi) aHermC1 = createHermitJacobi(15, GeomAbs_C1); + Handle(PLib_HermitJacobi) aHermC2 = createHermitJacobi(20, GeomAbs_C2); + + ASSERT_FALSE(aHermC0.IsNull()) << "Failed to create C0 HermitJacobi polynomial"; + ASSERT_FALSE(aHermC1.IsNull()) << "Failed to create C1 HermitJacobi polynomial"; + ASSERT_FALSE(aHermC2.IsNull()) << "Failed to create C2 HermitJacobi polynomial"; + + // Test WorkDegree property + EXPECT_EQ(aHermC0->WorkDegree(), 10); + EXPECT_EQ(aHermC1->WorkDegree(), 15); + EXPECT_EQ(aHermC2->WorkDegree(), 20); + + // Test NivConstr property + EXPECT_EQ(aHermC0->NivConstr(), 0); + EXPECT_EQ(aHermC1->NivConstr(), 1); + EXPECT_EQ(aHermC2->NivConstr(), 2); +} + +// Test basis function evaluation D0 +TEST_F(PLibHermitJacobiTest, BasisFunctionD0) +{ + Handle(PLib_HermitJacobi) aHerm = createHermitJacobi(6, GeomAbs_C0); + + TColStd_Array1OfReal aBasisValue(0, aHerm->WorkDegree()); + + // Test at various parameter values + std::vector aTestParams = {-1.0, -0.5, 0.0, 0.5, 1.0}; + + for (Standard_Real aU : aTestParams) + { + EXPECT_NO_THROW({ aHerm->D0(aU, aBasisValue); }) << "D0 evaluation failed at U=" << aU; + + // Basic sanity checks + for (Standard_Integer i = aBasisValue.Lower(); i <= aBasisValue.Upper(); i++) + { + EXPECT_FALSE(Precision::IsInfinite(aBasisValue(i))) + << "Basis value should be finite at index " << i << ", U=" << aU; + } + } +} + +// Test basis function evaluation with derivatives +TEST_F(PLibHermitJacobiTest, BasisFunctionDerivatives) +{ + Handle(PLib_HermitJacobi) aHerm = createHermitJacobi(8, GeomAbs_C1); + + TColStd_Array1OfReal aBasisValue(0, aHerm->WorkDegree()); + TColStd_Array1OfReal aBasisD1(0, aHerm->WorkDegree()); + TColStd_Array1OfReal aBasisD2(0, aHerm->WorkDegree()); + TColStd_Array1OfReal aBasisD3(0, aHerm->WorkDegree()); + + Standard_Real aU = 0.5; // Test at middle point + + // Test D1 + EXPECT_NO_THROW({ aHerm->D1(aU, aBasisValue, aBasisD1); }) << "D1 evaluation failed"; + + // Test D2 + EXPECT_NO_THROW({ aHerm->D2(aU, aBasisValue, aBasisD1, aBasisD2); }) << "D2 evaluation failed"; + + // Test D3 + EXPECT_NO_THROW({ aHerm->D3(aU, aBasisValue, aBasisD1, aBasisD2, aBasisD3); }) + << "D3 evaluation failed"; + + // Verify all values are finite + for (Standard_Integer i = aBasisValue.Lower(); i <= aBasisValue.Upper(); i++) + { + EXPECT_FALSE(Precision::IsInfinite(aBasisValue(i))) << "Basis value should be finite at " << i; + EXPECT_FALSE(Precision::IsInfinite(aBasisD1(i))) + << "First derivative should be finite at " << i; + EXPECT_FALSE(Precision::IsInfinite(aBasisD2(i))) + << "Second derivative should be finite at " << i; + EXPECT_FALSE(Precision::IsInfinite(aBasisD3(i))) + << "Third derivative should be finite at " << i; + } +} + +// Test coefficient conversion +TEST_F(PLibHermitJacobiTest, CoefficientConversion) +{ + const Standard_Integer aWorkDegree = 6; // Use smaller degree that works well with ToCoefficients + Handle(PLib_HermitJacobi) aHerm = createHermitJacobi(aWorkDegree, GeomAbs_C0); + + const Standard_Integer aDimension = 1; + const Standard_Integer aDegree = + aHerm->WorkDegree() - 2 * (aHerm->NivConstr() + 1); // Use computational degree + + // Create test HermitJacobi coefficients with proper size + // ToCoefficients expects arrays sized based on the degree and dimension + Standard_Integer aHermJacSize = (aDegree + 1) * aDimension; + Standard_Integer aCoeffSize = (aDegree + 1) * aDimension; + + // Use 0-based arrays to match the ToCoefficients indexing expectations + TColStd_Array1OfReal aHermJacCoeff(0, aHermJacSize - 1); + for (Standard_Integer i = aHermJacCoeff.Lower(); i <= aHermJacCoeff.Upper(); i++) + { + aHermJacCoeff(i) = Sin(i * 0.3); // Some test values + } + + TColStd_Array1OfReal aCoefficients(0, aCoeffSize - 1); + + EXPECT_NO_THROW({ aHerm->ToCoefficients(aDimension, aDegree, aHermJacCoeff, aCoefficients); }) + << "Coefficient conversion failed"; + + // Verify output is finite + for (Standard_Integer i = aCoefficients.Lower(); i <= aCoefficients.Upper(); i++) + { + EXPECT_FALSE(Precision::IsInfinite(aCoefficients(i))) + << "Converted coefficient should be finite at index " << i; + } +} + +// Test degree reduction +TEST_F(PLibHermitJacobiTest, DegreeReduction) +{ + Handle(PLib_HermitJacobi) aHerm = createHermitJacobi(10, GeomAbs_C0); + + const Standard_Integer aDimension = 1; + const Standard_Integer aMaxDegree = 8; + const Standard_Real aTol = 1e-6; + + // Create test coefficients - must be sized for full WorkDegree + const Standard_Integer aWorkDegree = aHerm->WorkDegree(); + TColStd_Array1OfReal aCoeff(1, (aWorkDegree + 1) * aDimension); + for (Standard_Integer i = aCoeff.Lower(); i <= aCoeff.Upper(); i++) + { + aCoeff(i) = 1.0 / (i + 1); // Decreasing coefficients to allow reduction + } + + Standard_Integer aNewDegree = -1; + Standard_Real aMaxError = -1.0; + + EXPECT_NO_THROW({ + aHerm->ReduceDegree(aDimension, aMaxDegree, aTol, aCoeff.ChangeValue(1), aNewDegree, aMaxError); + }) << "Degree reduction failed"; + + // Verify results are reasonable + EXPECT_LE(aNewDegree, aMaxDegree) << "New degree should not exceed max degree"; + EXPECT_GE(aNewDegree, 0) << "New degree should be non-negative"; + EXPECT_GE(aMaxError, 0.0) << "Max error should be non-negative"; + EXPECT_FALSE(Precision::IsInfinite(aMaxError)) << "Max error should be finite"; +} + +// Test error estimation +TEST_F(PLibHermitJacobiTest, ErrorEstimation) +{ + Handle(PLib_HermitJacobi) aHerm = createHermitJacobi(8, GeomAbs_C1); + + const Standard_Integer aDimension = 1; + + // Create test coefficients + TColStd_Array1OfReal aCoeff(1, 10 * aDimension); + for (Standard_Integer i = aCoeff.Lower(); i <= aCoeff.Upper(); i++) + { + aCoeff(i) = 1.0 / (i + 1); + } + + Standard_Integer aNewDegree = 6; // Reduced from original + + // Test MaxError + Standard_Real aMaxErr = -1.0; + EXPECT_NO_THROW({ aMaxErr = aHerm->MaxError(aDimension, aCoeff.ChangeValue(1), aNewDegree); }) + << "MaxError calculation failed"; + + EXPECT_GE(aMaxErr, 0.0) << "Max error should be non-negative"; + EXPECT_FALSE(Precision::IsInfinite(aMaxErr)) << "Max error should be finite"; + + // Test AverageError + Standard_Real aAvgErr = -1.0; + EXPECT_NO_THROW({ aAvgErr = aHerm->AverageError(aDimension, aCoeff.ChangeValue(1), aNewDegree); }) + << "AverageError calculation failed"; + + EXPECT_GE(aAvgErr, 0.0) << "Average error should be non-negative"; + EXPECT_FALSE(Precision::IsInfinite(aAvgErr)) << "Average error should be finite"; + + // Average error should typically be <= max error + EXPECT_LE(aAvgErr, aMaxErr + Precision::Confusion()) + << "Average error should not exceed max error significantly"; +} + +// Test extreme parameter values +TEST_F(PLibHermitJacobiTest, ExtremeParameterValues) +{ + Handle(PLib_HermitJacobi) aHerm = createHermitJacobi(10, GeomAbs_C2); + + TColStd_Array1OfReal aBasisValue(0, aHerm->WorkDegree()); + + // Test with boundary values + std::vector aExtremeParams = {-0.99999, -1e-12, 1e-12, 0.99999}; + + for (Standard_Real aU : aExtremeParams) + { + EXPECT_NO_THROW({ aHerm->D0(aU, aBasisValue); }) + << "Extreme parameter U=" << aU << " should not crash"; + + // Check that results are finite + for (Standard_Integer i = aBasisValue.Lower(); i <= aBasisValue.Upper(); i++) + { + EXPECT_FALSE(Precision::IsInfinite(aBasisValue(i))) + << "Basis value should be finite at extreme parameter U=" << aU; + } + } +} + +// Test consistency between different derivative orders +TEST_F(PLibHermitJacobiTest, DerivativeConsistency) +{ + Handle(PLib_HermitJacobi) aHerm = createHermitJacobi(6, GeomAbs_C2); + + TColStd_Array1OfReal aBasisValue1(0, aHerm->WorkDegree()); + TColStd_Array1OfReal aBasisD1_1(0, aHerm->WorkDegree()); + + TColStd_Array1OfReal aBasisValue2(0, aHerm->WorkDegree()); + TColStd_Array1OfReal aBasisD1_2(0, aHerm->WorkDegree()); + TColStd_Array1OfReal aBasisD2_2(0, aHerm->WorkDegree()); + + Standard_Real aU = 0.3; + + // Get values from D1 and D2 calls + aHerm->D1(aU, aBasisValue1, aBasisD1_1); + aHerm->D2(aU, aBasisValue2, aBasisD1_2, aBasisD2_2); + + // Values and first derivatives should be consistent + for (Standard_Integer i = aBasisValue1.Lower(); i <= aBasisValue1.Upper(); i++) + { + EXPECT_NEAR(aBasisValue1(i), aBasisValue2(i), Precision::Confusion()) + << "Function values should be consistent between D1 and D2 calls at index " << i; + EXPECT_NEAR(aBasisD1_1(i), aBasisD1_2(i), Precision::Confusion()) + << "First derivatives should be consistent between D1 and D2 calls at index " << i; + } +} + +// Performance test with maximum degree +TEST_F(PLibHermitJacobiTest, PerformanceTest) +{ + Handle(PLib_HermitJacobi) aHermMax = createHermitJacobi(30, GeomAbs_C2); + + TColStd_Array1OfReal aBasisValue(0, aHermMax->WorkDegree()); + TColStd_Array1OfReal aBasisD1(0, aHermMax->WorkDegree()); + TColStd_Array1OfReal aBasisD2(0, aHermMax->WorkDegree()); + TColStd_Array1OfReal aBasisD3(0, aHermMax->WorkDegree()); + + // Test that operations complete in reasonable time + std::vector aTestParams = {-0.8, -0.5, 0.0, 0.5, 0.8}; + + for (Standard_Real aU : aTestParams) + { + EXPECT_NO_THROW({ + aHermMax->D0(aU, aBasisValue); + aHermMax->D1(aU, aBasisValue, aBasisD1); + aHermMax->D2(aU, aBasisValue, aBasisD1, aBasisD2); + aHermMax->D3(aU, aBasisValue, aBasisD1, aBasisD2, aBasisD3); + }) << "High degree operations should complete without crashing at U=" + << aU; + } +} \ No newline at end of file diff --git a/src/FoundationClasses/TKMath/GTests/PLib_JacobiPolynomial_Test.cxx b/src/FoundationClasses/TKMath/GTests/PLib_JacobiPolynomial_Test.cxx new file mode 100644 index 0000000000..4e757c9003 --- /dev/null +++ b/src/FoundationClasses/TKMath/GTests/PLib_JacobiPolynomial_Test.cxx @@ -0,0 +1,366 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace +{ + +} // namespace + +// Test fixture for PLib_JacobiPolynomial tests +class PLibJacobiPolynomialTest : public ::testing::Test +{ +protected: + void SetUp() override + { + // Common test setup + } + + void TearDown() override + { + // Cleanup + } + + // Helper to create valid Jacobi polynomial instances + Handle(PLib_JacobiPolynomial) createJacobiPolynomial(Standard_Integer theDegree, + GeomAbs_Shape theConstraint) + { + // Ensure degree is within valid range (typically <= 30) + EXPECT_LE(theDegree, 30) << "Degree too high for Jacobi polynomial"; + EXPECT_GE(theDegree, 0) << "Degree must be non-negative"; + + return new PLib_JacobiPolynomial(theDegree, theConstraint); + } +}; + +// Test constructor and basic properties +TEST_F(PLibJacobiPolynomialTest, ConstructorAndBasicProperties) +{ + // Test with different constraint orders + Handle(PLib_JacobiPolynomial) aJacC0 = createJacobiPolynomial(10, GeomAbs_C0); + Handle(PLib_JacobiPolynomial) aJacC1 = createJacobiPolynomial(15, GeomAbs_C1); + Handle(PLib_JacobiPolynomial) aJacC2 = createJacobiPolynomial(20, GeomAbs_C2); + + ASSERT_FALSE(aJacC0.IsNull()) << "Failed to create C0 Jacobi polynomial"; + ASSERT_FALSE(aJacC1.IsNull()) << "Failed to create C1 Jacobi polynomial"; + ASSERT_FALSE(aJacC2.IsNull()) << "Failed to create C2 Jacobi polynomial"; + + // Test WorkDegree property + EXPECT_EQ(aJacC0->WorkDegree(), 10); + EXPECT_EQ(aJacC1->WorkDegree(), 15); + EXPECT_EQ(aJacC2->WorkDegree(), 20); + + // Test NivConstr property + EXPECT_EQ(aJacC0->NivConstr(), 0); + EXPECT_EQ(aJacC1->NivConstr(), 1); + EXPECT_EQ(aJacC2->NivConstr(), 2); +} + +// Test constructor with edge cases +TEST_F(PLibJacobiPolynomialTest, ConstructorEdgeCases) +{ + // Test minimum degree + Handle(PLib_JacobiPolynomial) aJacMin = createJacobiPolynomial(0, GeomAbs_C0); + EXPECT_FALSE(aJacMin.IsNull()); + EXPECT_EQ(aJacMin->WorkDegree(), 0); + + // Test maximum recommended degree + Handle(PLib_JacobiPolynomial) aJacMax = createJacobiPolynomial(30, GeomAbs_C2); + EXPECT_FALSE(aJacMax.IsNull()); + EXPECT_EQ(aJacMax->WorkDegree(), 30); + + // Test reasonable high degrees + Handle(PLib_JacobiPolynomial) aJacHigh = createJacobiPolynomial(25, GeomAbs_C0); + EXPECT_GT(aJacHigh->WorkDegree(), 20) << "High degree should be supported"; +} + +// Test Gauss integration points +TEST_F(PLibJacobiPolynomialTest, GaussIntegrationPoints) +{ + Handle(PLib_JacobiPolynomial) aJac = createJacobiPolynomial(10, GeomAbs_C0); + + // Test various numbers of Gauss points (only valid values supported by OCCT) + std::vector aGaussNumbers = {8, 10, 15, 20, 25, 30, 40, 50, 61}; + + for (Standard_Integer aNbGauss : aGaussNumbers) + { + if (aNbGauss > aJac->WorkDegree()) + { + TColStd_Array1OfReal aPoints(0, aNbGauss / 2); + + EXPECT_NO_THROW({ aJac->Points(aNbGauss, aPoints); }) + << "Points calculation failed for " << aNbGauss << " Gauss points"; + + // Verify points are in valid range and ordered + for (Standard_Integer i = aPoints.Lower(); i <= aPoints.Upper(); i++) + { + if (i == 0 && aNbGauss % 2 == 0) + { + // For even number of Gauss points, TabPoints(0) is set to UNDEFINED (-999) + EXPECT_EQ(aPoints(i), -999) << "TabPoints(0) should be UNDEFINED for even NbGaussPoints"; + } + else if (i == 0 && aNbGauss % 2 == 1) + { + // For odd number of Gauss points, TabPoints(0) should be 0 + EXPECT_EQ(aPoints(i), 0.0) << "TabPoints(0) should be 0 for odd NbGaussPoints"; + } + else + { + // Other points should be positive and <= 1 + EXPECT_GT(aPoints(i), 0.0) << "Gauss point should be positive"; + EXPECT_LE(aPoints(i), 1.0) << "Gauss point should be <= 1"; + + if (i > aPoints.Lower() && i > 0) + { + EXPECT_GE(aPoints(i), aPoints(i - 1)) << "Gauss points should be ordered"; + } + } + } + } + } +} + +// Test Gauss integration weights +TEST_F(PLibJacobiPolynomialTest, GaussIntegrationWeights) +{ + Handle(PLib_JacobiPolynomial) aJac = createJacobiPolynomial(8, GeomAbs_C1); + + Standard_Integer aNbGauss = 15; // Must be > degree for valid computation + TColStd_Array2OfReal aWeights(0, aNbGauss / 2, 0, aJac->WorkDegree()); + + aJac->Weights(aNbGauss, aWeights); + + // Basic sanity checks on weights - the array is 2D with specific bounds + EXPECT_EQ(aWeights.LowerRow(), 0) << "Lower row should be 0"; + EXPECT_EQ(aWeights.UpperRow(), aNbGauss / 2) << "Upper row mismatch"; + EXPECT_EQ(aWeights.LowerCol(), 0) << "Lower col should be 0"; + EXPECT_EQ(aWeights.UpperCol(), aJac->WorkDegree()) << "Upper col should match work degree"; + + for (Standard_Integer i = aWeights.LowerRow(); i <= aWeights.UpperRow(); i++) + { + for (Standard_Integer j = aWeights.LowerCol(); j <= aWeights.UpperCol(); j++) + { + EXPECT_FALSE(Precision::IsInfinite(aWeights(i, j))) + << "Weight should be finite at (" << i << "," << j << ")"; + } + } +} + +// Test MaxValue computation +TEST_F(PLibJacobiPolynomialTest, MaxValue) +{ + Handle(PLib_JacobiPolynomial) aJac = createJacobiPolynomial(10, GeomAbs_C0); + + Standard_Integer aTabSize = aJac->WorkDegree() - 2 * (aJac->NivConstr() + 1); + if (aTabSize > 0) + { + TColStd_Array1OfReal aTabMax(0, aTabSize); + + aJac->MaxValue(aTabMax); + + // Verify all max values are positive (they represent maximum absolute values) + for (Standard_Integer i = aTabMax.Lower(); i <= aTabMax.Upper(); i++) + { + EXPECT_GT(aTabMax(i), 0.0) << "Max value should be positive at index " << i; + EXPECT_FALSE(Precision::IsInfinite(aTabMax(i))) + << "Max value should be finite at index " << i; + } + } +} + +// Test basis function evaluation D0 +TEST_F(PLibJacobiPolynomialTest, BasisFunctionD0) +{ + Handle(PLib_JacobiPolynomial) aJac = createJacobiPolynomial(6, GeomAbs_C0); + + TColStd_Array1OfReal aBasisValue(0, aJac->WorkDegree()); + + // Test at various parameter values + std::vector aTestParams = {-1.0, -0.5, 0.0, 0.5, 1.0}; + + for (Standard_Real aU : aTestParams) + { + aJac->D0(aU, aBasisValue); + + // Basic sanity checks + for (Standard_Integer i = aBasisValue.Lower(); i <= aBasisValue.Upper(); i++) + { + EXPECT_FALSE(Precision::IsInfinite(aBasisValue(i))) + << "Basis value should be finite at index " << i << ", U=" << aU; + } + } +} + +// Test basis function evaluation with derivatives +TEST_F(PLibJacobiPolynomialTest, BasisFunctionDerivatives) +{ + Handle(PLib_JacobiPolynomial) aJac = createJacobiPolynomial(8, GeomAbs_C1); + + TColStd_Array1OfReal aBasisValue(0, aJac->WorkDegree()); + TColStd_Array1OfReal aBasisD1(0, aJac->WorkDegree()); + TColStd_Array1OfReal aBasisD2(0, aJac->WorkDegree()); + TColStd_Array1OfReal aBasisD3(0, aJac->WorkDegree()); + + Standard_Real aU = 0.5; // Test at middle point + + // Test D1, D2, D3 evaluations + aJac->D1(aU, aBasisValue, aBasisD1); + aJac->D2(aU, aBasisValue, aBasisD1, aBasisD2); + aJac->D3(aU, aBasisValue, aBasisD1, aBasisD2, aBasisD3); + + // Verify all values are finite + for (Standard_Integer i = aBasisValue.Lower(); i <= aBasisValue.Upper(); i++) + { + EXPECT_FALSE(Precision::IsInfinite(aBasisValue(i))) << "Basis value should be finite at " << i; + EXPECT_FALSE(Precision::IsInfinite(aBasisD1(i))) + << "First derivative should be finite at " << i; + EXPECT_FALSE(Precision::IsInfinite(aBasisD2(i))) + << "Second derivative should be finite at " << i; + EXPECT_FALSE(Precision::IsInfinite(aBasisD3(i))) + << "Third derivative should be finite at " << i; + } +} + +// Test coefficient conversion +TEST_F(PLibJacobiPolynomialTest, CoefficientConversion) +{ + const Standard_Integer aWorkDegree = 6; // Use smaller degree that works well with ToCoefficients + Handle(PLib_JacobiPolynomial) aJac = createJacobiPolynomial(aWorkDegree, GeomAbs_C0); + + const Standard_Integer aDimension = 1; + const Standard_Integer aDegree = + aJac->WorkDegree() - 2 * (aJac->NivConstr() + 1); // Use computational degree + + // Create test Jacobi coefficients with proper size + // ToCoefficients expects arrays sized based on the degree and dimension + // Analysis shows we need arrays that can handle the indexing patterns in the method + Standard_Integer aJacSize = (aDegree + 1) * aDimension; + Standard_Integer aCoeffSize = (aDegree + 1) * aDimension; + + // Use 0-based arrays to match the ToCoefficients indexing expectations + TColStd_Array1OfReal aJacCoeff(0, aJacSize - 1); + for (Standard_Integer i = aJacCoeff.Lower(); i <= aJacCoeff.Upper(); i++) + { + aJacCoeff(i) = Sin(i * 0.1); // Some test values + } + + TColStd_Array1OfReal aCoefficients(0, aCoeffSize - 1); + + aJac->ToCoefficients(aDimension, aDegree, aJacCoeff, aCoefficients); + + // Verify output is finite + for (Standard_Integer i = aCoefficients.Lower(); i <= aCoefficients.Upper(); i++) + { + EXPECT_FALSE(Precision::IsInfinite(aCoefficients(i))) + << "Converted coefficient should be finite at index " << i; + } +} + +// Test degree reduction +TEST_F(PLibJacobiPolynomialTest, DegreeReduction) +{ + Handle(PLib_JacobiPolynomial) aJac = createJacobiPolynomial(10, GeomAbs_C0); + + const Standard_Integer aDimension = 1; + const Standard_Integer aMaxDegree = 8; + const Standard_Real aTol = 1e-6; + + // Create test coefficients - must be sized for full WorkDegree + const Standard_Integer aWorkDegree = aJac->WorkDegree(); + TColStd_Array1OfReal aCoeff(1, (aWorkDegree + 1) * aDimension); + for (Standard_Integer i = aCoeff.Lower(); i <= aCoeff.Upper(); i++) + { + aCoeff(i) = 1.0 / (i + 1); // Decreasing coefficients to allow reduction + } + + Standard_Integer aNewDegree = -1; + Standard_Real aMaxError = -1.0; + + aJac->ReduceDegree(aDimension, aMaxDegree, aTol, aCoeff.ChangeValue(1), aNewDegree, aMaxError); + + // Verify results are reasonable + EXPECT_LE(aNewDegree, aMaxDegree) << "New degree should not exceed max degree"; + EXPECT_GE(aNewDegree, 0) << "New degree should be non-negative"; + EXPECT_GE(aMaxError, 0.0) << "Max error should be non-negative"; + EXPECT_FALSE(Precision::IsInfinite(aMaxError)) << "Max error should be finite"; +} + +// Test error estimation +TEST_F(PLibJacobiPolynomialTest, ErrorEstimation) +{ + Handle(PLib_JacobiPolynomial) aJac = createJacobiPolynomial(8, GeomAbs_C1); + + const Standard_Integer aDimension = 1; + + // Create test coefficients + TColStd_Array1OfReal aCoeff(1, 10 * aDimension); + for (Standard_Integer i = aCoeff.Lower(); i <= aCoeff.Upper(); i++) + { + aCoeff(i) = 1.0 / (i + 1); + } + + Standard_Integer aNewDegree = 6; // Reduced from original + + // Test MaxError + Standard_Real aMaxErr = aJac->MaxError(aDimension, aCoeff.ChangeValue(1), aNewDegree); + + EXPECT_GE(aMaxErr, 0.0) << "Max error should be non-negative"; + EXPECT_FALSE(Precision::IsInfinite(aMaxErr)) << "Max error should be finite"; + + // Test AverageError + Standard_Real aAvgErr = aJac->AverageError(aDimension, aCoeff.ChangeValue(1), aNewDegree); + + EXPECT_GE(aAvgErr, 0.0) << "Average error should be non-negative"; + EXPECT_FALSE(Precision::IsInfinite(aAvgErr)) << "Average error should be finite"; + + // Average error should typically be <= max error + EXPECT_LE(aAvgErr, aMaxErr + Precision::Confusion()) + << "Average error should not exceed max error significantly"; +} + +// Performance and stress tests +TEST_F(PLibJacobiPolynomialTest, StressTests) +{ + // Test with maximum degree + Handle(PLib_JacobiPolynomial) aJacMax = createJacobiPolynomial(30, GeomAbs_C2); + + // Test that basic operations work with high degrees + TColStd_Array1OfReal aBasisValue(0, aJacMax->WorkDegree()); + + aJacMax->D0(0.0, aBasisValue); + aJacMax->D0(0.5, aBasisValue); + aJacMax->D0(1.0, aBasisValue); + + // Test with extreme parameter values + std::vector aExtremeParams = {-0.99999, -1e-10, 1e-10, 0.99999}; + + for (Standard_Real aU : aExtremeParams) + { + aJacMax->D0(aU, aBasisValue); + // Verify basis values are finite + for (Standard_Integer i = aBasisValue.Lower(); i <= aBasisValue.Upper(); i++) + { + EXPECT_FALSE(Precision::IsInfinite(aBasisValue(i))) << "Basis value should be finite"; + } + } +} \ No newline at end of file diff --git a/src/FoundationClasses/TKMath/GTests/PLib_Test.cxx b/src/FoundationClasses/TKMath/GTests/PLib_Test.cxx new file mode 100644 index 0000000000..4809c85aef --- /dev/null +++ b/src/FoundationClasses/TKMath/GTests/PLib_Test.cxx @@ -0,0 +1,418 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace +{ +// Helper function for comparing points with tolerance +void checkPointsEqual(const gp_Pnt& thePnt1, + const gp_Pnt& thePnt2, + const Standard_Real theTolerance = Precision::Confusion()) +{ + EXPECT_NEAR(thePnt1.X(), thePnt2.X(), theTolerance) << "X coordinates differ"; + EXPECT_NEAR(thePnt1.Y(), thePnt2.Y(), theTolerance) << "Y coordinates differ"; + EXPECT_NEAR(thePnt1.Z(), thePnt2.Z(), theTolerance) << "Z coordinates differ"; +} + +// Helper function for comparing 2D points with tolerance +void checkPoint2dEqual(const gp_Pnt2d& thePnt1, + const gp_Pnt2d& thePnt2, + const Standard_Real theTolerance = Precision::Confusion()) +{ + EXPECT_NEAR(thePnt1.X(), thePnt2.X(), theTolerance) << "X coordinates differ"; + EXPECT_NEAR(thePnt1.Y(), thePnt2.Y(), theTolerance) << "Y coordinates differ"; +} +} // namespace + +// Test class for PLib tests +class PLibTest : public ::testing::Test +{ +protected: + void SetUp() override + { + // Setup common test data + } + + void TearDown() override + { + // Cleanup + } +}; + +// Tests for basic utility functions +TEST_F(PLibTest, NoWeightsPointers) +{ + // Test that NoWeights() returns NULL as expected + EXPECT_EQ(PLib::NoWeights(), static_cast(nullptr)); + EXPECT_EQ(PLib::NoWeights2(), static_cast(nullptr)); +} + +// Tests for 3D pole conversion functions +TEST_F(PLibTest, SetGetPoles3D) +{ + // Create test data + TColgp_Array1OfPnt aPoles(1, 3); + aPoles(1) = gp_Pnt(1.0, 2.0, 3.0); + aPoles(2) = gp_Pnt(4.0, 5.0, 6.0); + aPoles(3) = gp_Pnt(7.0, 8.0, 9.0); + + // Test SetPoles and GetPoles without weights + TColStd_Array1OfReal aFP(1, 9); // 3 poles * 3 coordinates + PLib::SetPoles(aPoles, aFP); + + TColgp_Array1OfPnt aResultPoles(1, 3); + PLib::GetPoles(aFP, aResultPoles); + + // Verify results + for (Standard_Integer i = 1; i <= 3; i++) + { + checkPointsEqual(aPoles(i), aResultPoles(i)); + } +} + +TEST_F(PLibTest, SetGetPoles3DWithWeights) +{ + // Create test data + TColgp_Array1OfPnt aPoles(1, 2); + aPoles(1) = gp_Pnt(1.0, 2.0, 3.0); + aPoles(2) = gp_Pnt(4.0, 5.0, 6.0); + + TColStd_Array1OfReal aWeights(1, 2); + aWeights(1) = 0.5; + aWeights(2) = 2.0; + + // Test SetPoles and GetPoles with weights + TColStd_Array1OfReal aFP(1, 8); // 2 poles * (3 coords + 1 weight) + PLib::SetPoles(aPoles, aWeights, aFP); + + TColgp_Array1OfPnt aResultPoles(1, 2); + TColStd_Array1OfReal aResultWeights(1, 2); + PLib::GetPoles(aFP, aResultPoles, aResultWeights); + + // Verify results + for (Standard_Integer i = 1; i <= 2; i++) + { + checkPointsEqual(aPoles(i), aResultPoles(i)); + EXPECT_NEAR(aWeights(i), aResultWeights(i), Precision::Confusion()); + } +} + +// Tests for 2D pole conversion functions +TEST_F(PLibTest, SetGetPoles2D) +{ + // Create test data + TColgp_Array1OfPnt2d aPoles(1, 3); + aPoles(1) = gp_Pnt2d(1.0, 2.0); + aPoles(2) = gp_Pnt2d(3.0, 4.0); + aPoles(3) = gp_Pnt2d(5.0, 6.0); + + // Test SetPoles and GetPoles without weights + TColStd_Array1OfReal aFP(1, 6); // 3 poles * 2 coordinates + PLib::SetPoles(aPoles, aFP); + + TColgp_Array1OfPnt2d aResultPoles(1, 3); + PLib::GetPoles(aFP, aResultPoles); + + // Verify results + for (Standard_Integer i = 1; i <= 3; i++) + { + checkPoint2dEqual(aPoles(i), aResultPoles(i)); + } +} + +TEST_F(PLibTest, SetGetPoles2DWithWeights) +{ + // Create test data + TColgp_Array1OfPnt2d aPoles(1, 2); + aPoles(1) = gp_Pnt2d(1.0, 2.0); + aPoles(2) = gp_Pnt2d(3.0, 4.0); + + TColStd_Array1OfReal aWeights(1, 2); + aWeights(1) = 0.8; + aWeights(2) = 1.5; + + // Test SetPoles and GetPoles with weights + TColStd_Array1OfReal aFP(1, 6); // 2 poles * (2 coords + 1 weight) + PLib::SetPoles(aPoles, aWeights, aFP); + + TColgp_Array1OfPnt2d aResultPoles(1, 2); + TColStd_Array1OfReal aResultWeights(1, 2); + PLib::GetPoles(aFP, aResultPoles, aResultWeights); + + // Verify results + for (Standard_Integer i = 1; i <= 2; i++) + { + checkPoint2dEqual(aPoles(i), aResultPoles(i)); + EXPECT_NEAR(aWeights(i), aResultWeights(i), Precision::Confusion()); + } +} + +// Test for zero weights handling (safety test) +TEST_F(PLibTest, ZeroWeightsHandling) +{ + // Create test data with zero weight - this should not crash + TColgp_Array1OfPnt2d aPoles(1, 2); + aPoles(1) = gp_Pnt2d(1.0, 2.0); + aPoles(2) = gp_Pnt2d(3.0, 4.0); + + TColStd_Array1OfReal aWeights(1, 2); + aWeights(1) = 1.0; + aWeights(2) = 0.0; // Zero weight - potential division by zero + + TColStd_Array1OfReal aFP(1, 6); + PLib::SetPoles(aPoles, aWeights, aFP); + + // This test should complete without crashing + // The behavior with zero weights may vary, but it shouldn't crash + EXPECT_TRUE(true); // We mainly test that no crash occurs +} + +// Tests for Binomial coefficient function +TEST_F(PLibTest, BinomialCoefficient) +{ + // Test known binomial coefficients + EXPECT_NEAR(PLib::Bin(0, 0), 1.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(5, 0), 1.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(5, 5), 1.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(5, 1), 5.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(5, 2), 10.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(5, 3), 10.0, Precision::Confusion()); + EXPECT_NEAR(PLib::Bin(10, 3), 120.0, Precision::Confusion()); + + // Test symmetry property: C(n,k) = C(n,n-k) + for (Standard_Integer n = 1; n <= 10; n++) + { + for (Standard_Integer k = 0; k <= n; k++) + { + EXPECT_NEAR(PLib::Bin(n, k), PLib::Bin(n, n - k), Precision::Confusion()) + << "Binomial coefficient symmetry failed for C(" << n << "," << k << ")"; + } + } +} + +// Tests for polynomial evaluation +TEST_F(PLibTest, EvalPolynomial) +{ + // Test simple polynomial evaluation: f(x) = 1 + 2x + 3x^2 + const Standard_Integer aDegree = 2; + const Standard_Integer aDimension = 1; + const Standard_Integer aDerivOrder = 0; + + TColStd_Array1OfReal aCoeffs(1, 3); + aCoeffs(1) = 1.0; // constant term + aCoeffs(2) = 2.0; // linear term + aCoeffs(3) = 3.0; // quadratic term + + TColStd_Array1OfReal aResults(1, 1); + + // Test at x = 0: f(0) = 1 + PLib::EvalPolynomial(0.0, + aDerivOrder, + aDegree, + aDimension, + aCoeffs.ChangeValue(1), + aResults.ChangeValue(1)); + EXPECT_NEAR(aResults(1), 1.0, Precision::Confusion()); + + // Test at x = 1: f(1) = 1 + 2 + 3 = 6 + PLib::EvalPolynomial(1.0, + aDerivOrder, + aDegree, + aDimension, + aCoeffs.ChangeValue(1), + aResults.ChangeValue(1)); + EXPECT_NEAR(aResults(1), 6.0, Precision::Confusion()); + + // Test at x = 2: f(2) = 1 + 4 + 12 = 17 + PLib::EvalPolynomial(2.0, + aDerivOrder, + aDegree, + aDimension, + aCoeffs.ChangeValue(1), + aResults.ChangeValue(1)); + EXPECT_NEAR(aResults(1), 17.0, Precision::Confusion()); +} + +// Tests for polynomial evaluation with derivatives +TEST_F(PLibTest, EvalPolynomialWithDerivatives) +{ + // Test polynomial f(x) = 1 + 2x + 3x^2 + // f'(x) = 2 + 6x + // f''(x) = 6 + const Standard_Integer aDegree = 2; + const Standard_Integer aDimension = 1; + const Standard_Integer aDerivOrder = 2; + + TColStd_Array1OfReal aCoeffs(1, 3); + aCoeffs(1) = 1.0; + aCoeffs(2) = 2.0; + aCoeffs(3) = 3.0; + + TColStd_Array1OfReal aResults(1, 3); // value + 1st + 2nd derivative + + // Test at x = 1 + PLib::EvalPolynomial(1.0, + aDerivOrder, + aDegree, + aDimension, + aCoeffs.ChangeValue(1), + aResults.ChangeValue(1)); + + EXPECT_NEAR(aResults(1), 6.0, Precision::Confusion()); // f(1) = 6 + EXPECT_NEAR(aResults(2), 8.0, Precision::Confusion()); // f'(1) = 2 + 6 = 8 + EXPECT_NEAR(aResults(3), 6.0, Precision::Confusion()); // f''(1) = 6 +} + +// Tests for constraint order conversion functions +TEST_F(PLibTest, ConstraintOrderConversion) +{ + // Test NivConstr function + EXPECT_EQ(PLib::NivConstr(GeomAbs_C0), 0); + EXPECT_EQ(PLib::NivConstr(GeomAbs_C1), 1); + EXPECT_EQ(PLib::NivConstr(GeomAbs_C2), 2); + + // Test ConstraintOrder function + EXPECT_EQ(PLib::ConstraintOrder(0), GeomAbs_C0); + EXPECT_EQ(PLib::ConstraintOrder(1), GeomAbs_C1); + EXPECT_EQ(PLib::ConstraintOrder(2), GeomAbs_C2); + + // Test round-trip consistency + for (Standard_Integer i = 0; i <= 2; i++) + { + GeomAbs_Shape aShape = PLib::ConstraintOrder(i); + Standard_Integer aLevel = PLib::NivConstr(aShape); + EXPECT_EQ(aLevel, i) << "Round-trip conversion failed for level " << i; + } +} + +// Test for Hermite interpolation +TEST_F(PLibTest, HermiteInterpolate) +{ + + const Standard_Integer aDimension = 1; + const Standard_Real aFirstParam = 0.0; + const Standard_Real aLastParam = 1.0; + const Standard_Integer aFirstOrder = 1; // value + 1st derivative + const Standard_Integer aLastOrder = 1; // value + 1st derivative + + // Define constraints: f(0) = 0, f'(0) = 1, f(1) = 1, f'(1) = 0 + TColStd_Array2OfReal aFirstConstr(1, aDimension, 0, aFirstOrder); + aFirstConstr(1, 0) = 0.0; // f(0) = 0 + aFirstConstr(1, 1) = 1.0; // f'(0) = 1 + + TColStd_Array2OfReal aLastConstr(1, aDimension, 0, aLastOrder); + aLastConstr(1, 0) = 1.0; // f(1) = 1 + aLastConstr(1, 1) = 0.0; // f'(1) = 0 + + const Standard_Integer aCoeffCount = aFirstOrder + aLastOrder + 2; + + TColStd_Array1OfReal aCoeffs(0, + aCoeffCount + - 1); // 0-based indexing as expected by HermiteInterpolate + + // Perform Hermite interpolation + Standard_Boolean aResult = PLib::HermiteInterpolate(aDimension, + aFirstParam, + aLastParam, + aFirstOrder, + aLastOrder, + aFirstConstr, + aLastConstr, + aCoeffs); + + EXPECT_TRUE(aResult) << "Hermite interpolation failed"; + + if (aResult) + { + + // Verify that the resulting polynomial satisfies the constraints + TColStd_Array1OfReal aResults(1, 2); // value + 1st derivative + + // Check at first parameter + PLib::EvalPolynomial(aFirstParam, + 1, + aCoeffCount - 1, + aDimension, + aCoeffs.ChangeValue(0), + aResults.ChangeValue(1)); + EXPECT_NEAR(aResults(1), 0.0, Precision::Confusion()) << "f(0) constraint not satisfied"; + EXPECT_NEAR(aResults(2), 1.0, Precision::Confusion()) << "f'(0) constraint not satisfied"; + + // Check at last parameter + PLib::EvalPolynomial(aLastParam, + 1, + aCoeffCount - 1, + aDimension, + aCoeffs.ChangeValue(0), + aResults.ChangeValue(1)); + EXPECT_NEAR(aResults(1), 1.0, Precision::Confusion()) << "f(1) constraint not satisfied"; + EXPECT_NEAR(aResults(2), 0.0, Precision::Confusion()) << "f'(1) constraint not satisfied"; + } +} + +// Edge case tests +TEST_F(PLibTest, EdgeCases) +{ + // Test with very small coefficients + TColStd_Array1OfReal aSmallCoeffs(1, 3); + aSmallCoeffs(1) = 1.0e-12; + aSmallCoeffs(2) = 1.0e-12; + aSmallCoeffs(3) = 1.0e-12; + + TColStd_Array1OfReal aResults(1, 1); + + // This should not crash even with very small coefficients + EXPECT_NO_THROW( + { PLib::EvalPolynomial(1.0, 0, 2, 1, aSmallCoeffs.ChangeValue(1), aResults.ChangeValue(1)); }); + + // Test with large coefficients + TColStd_Array1OfReal aLargeCoeffs(1, 3); + aLargeCoeffs(1) = 1.0e10; + aLargeCoeffs(2) = 1.0e10; + aLargeCoeffs(3) = 1.0e10; + + EXPECT_NO_THROW( + { PLib::EvalPolynomial(1.0, 0, 2, 1, aLargeCoeffs.ChangeValue(1), aResults.ChangeValue(1)); }); +} + +// Test for Jacobi parameter calculation +TEST_F(PLibTest, JacobiParameters) +{ + Standard_Integer aNbGaussPoints, aWorkDegree; + + // Test various constraint orders and codes + PLib::JacobiParameters(GeomAbs_C0, 10, 1, aNbGaussPoints, aWorkDegree); + EXPECT_GT(aNbGaussPoints, 0) << "Number of Gauss points should be positive"; + EXPECT_GT(aWorkDegree, 0) << "Work degree should be positive"; + + PLib::JacobiParameters(GeomAbs_C1, 15, 2, aNbGaussPoints, aWorkDegree); + EXPECT_GT(aNbGaussPoints, 0); + EXPECT_GT(aWorkDegree, 0); + + PLib::JacobiParameters(GeomAbs_C2, 20, 3, aNbGaussPoints, aWorkDegree); + EXPECT_GT(aNbGaussPoints, 0); + EXPECT_GT(aWorkDegree, 0); +} \ No newline at end of file diff --git a/src/FoundationClasses/TKMath/PLib/PLib_DoubleJacobiPolynomial.cxx b/src/FoundationClasses/TKMath/PLib/PLib_DoubleJacobiPolynomial.cxx index cffc9e07c9..ee08f5caf7 100644 --- a/src/FoundationClasses/TKMath/PLib/PLib_DoubleJacobiPolynomial.cxx +++ b/src/FoundationClasses/TKMath/PLib/PLib_DoubleJacobiPolynomial.cxx @@ -173,6 +173,8 @@ void PLib_DoubleJacobiPolynomial::ReduceDegree(const Standard_Integer Dimen NewV = MaxDegreeV; math_Vector MaxErr2(1, 2); + MaxError = 0.0; // Initialize MaxError + //********************************************************************** //-------------------- Coupure des coefficients ------------------------ //**********************************************************************