- 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.
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
)
--- /dev/null
+// 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 <PLib_DoubleJacobiPolynomial.hxx>
+#include <PLib_JacobiPolynomial.hxx>
+
+#include <gtest/gtest.h>
+
+#include <Standard_Real.hxx>
+#include <Standard_Integer.hxx>
+#include <TColStd_Array1OfReal.hxx>
+#include <TColStd_HArray1OfReal.hxx>
+#include <GeomAbs_Shape.hxx>
+#include <Precision.hxx>
+
+// 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
--- /dev/null
+// 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 <PLib_HermitJacobi.hxx>
+#include <PLib_JacobiPolynomial.hxx>
+
+#include <gtest/gtest.h>
+
+#include <Standard_Real.hxx>
+#include <Standard_Integer.hxx>
+#include <TColStd_Array1OfReal.hxx>
+#include <GeomAbs_Shape.hxx>
+#include <Precision.hxx>
+
+// 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<Standard_Real> 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<Standard_Real> 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<Standard_Real> 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
--- /dev/null
+// 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 <PLib_JacobiPolynomial.hxx>
+
+#include <gtest/gtest.h>
+
+#include <Standard_Real.hxx>
+#include <Standard_Integer.hxx>
+#include <TColStd_Array1OfReal.hxx>
+#include <TColStd_Array2OfReal.hxx>
+#include <GeomAbs_Shape.hxx>
+#include <Precision.hxx>
+
+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<Standard_Integer> 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<Standard_Real> 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<Standard_Real> 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
--- /dev/null
+// 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 <PLib.hxx>
+
+#include <gtest/gtest.h>
+
+#include <Standard_Real.hxx>
+#include <Standard_Integer.hxx>
+#include <TColStd_Array1OfReal.hxx>
+#include <TColStd_Array2OfReal.hxx>
+#include <TColgp_Array1OfPnt.hxx>
+#include <TColgp_Array1OfPnt2d.hxx>
+#include <TColgp_Array2OfPnt.hxx>
+#include <GeomAbs_Shape.hxx>
+#include <Precision.hxx>
+
+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<TColStd_Array1OfReal*>(nullptr));
+ EXPECT_EQ(PLib::NoWeights2(), static_cast<TColStd_Array2OfReal*>(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
NewV = MaxDegreeV;
math_Vector MaxErr2(1, 2);
+ MaxError = 0.0; // Initialize MaxError
+
//**********************************************************************
//-------------------- Coupure des coefficients ------------------------
//**********************************************************************
projects / occt.git / commitdiff