1
0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-09-13 14:27:08 +03:00

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.
This commit is contained in:
Pasukhin Dmitry
2025-09-07 19:17:13 +01:00
committed by GitHub
parent 80ce783ca2
commit 7fdd2b62f9
6 changed files with 1502 additions and 0 deletions

View File

@@ -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
)

View File

@@ -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 <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";
}

View File

@@ -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 <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;
}
}

View File

@@ -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 <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";
}
}
}

View File

@@ -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 <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);
}

View File

@@ -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 ------------------------
//**********************************************************************