lili/OCCT
1
0
Fork 0
mirror of https://github.com/Open-Cascade-SAS/OCCT.git synced 2026-09-06 20:58:02 +08:00
Code Issues Packages Projects Releases Wiki Activity

Modeling Data - Optimize BSplCLib interpolation and blend evaluation erformance (#1082)

Browse Source 
Profiling identified several performance bottlenecks in the BSpline interpolation and blend surface computation pipeline. This commit addresses them through four categories of optimization:

1. Static initialization for GeomFill convertors: the monomial-to-BSpline conversion matrices in GeomFill_QuasiAngularConvertor::Init() and GeomFill_PolynomialConvertor::Init() are mathematical constants that were recomputed on every call via Convert_CompPolynomialToPoles. Now computed once via static lambda-initialized locals.

2. Stack allocation for small matrices/arrays: InterpolationMatrix in BSplCLib::Interpolate, aBSplineBasis in BuildBSpMatrix, and parameters/contact_array in Convert_CompPolynomialToPoles::Perform now use stack buffers when sizes fit, avoiding heap allocation.

3. Raw pointer access in hot loops: replaced multi-layer accessor chains (math_Matrix::Value -> math_DoubleTab::Value -> NCollection_Array2::Value -> NCollection_Array1::at with bounds checks) with direct pointer arithmetic in EvalBsplineBasis, FactorBandedMatrix, BuildBSpMatrix, SolveBandedSystem, and math_VectorBase operations (Multiply, TMultiply, Multiplied, Norm, Norm2).

4. Eliminated redundant recomputation: cached AdvApprox_ApproxAFunction:: NbPoles() results in Approx_SweepApproximation, Approx_CurveOnSurface, and Approx_Curve2d instead of recomputing BSplCLib::NbPoles in inner loops. Cached math_FunctionSetRoot solver in BRepBlend_AppFuncRoot to avoid repeated construction/destruction per SearchPoint call.

Also fixed undefined behavior in BSplCLib::NbPoles where pointer arithmetic created a pointer before the array start (pmu -= f).
This commit is contained in:
Pasukhin Dmitry
2026-02-13 21:35:03 +00:00
committed by GitHub
parent a9f45dfa15
commit ea9443d154
11 changed files with 305 additions and 259 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/FoundationClasses/TKMath/BSplCLib/BSplCLib.cxx
+37 -22
View File
@@ -361,17 +361,16 @@ int BSplCLib::MinKnotMult(const Array1OfInteger& Mults, const int FromK1, const
int BSplCLib::NbPoles(const int Degree, const bool Periodic, const NCollection_Array1<int>& Mults)
{
int i, sigma = 0;
int f = Mults.Lower();
int l = Mults.Upper();
const int f = Mults.Lower();
const int l = Mults.Upper();
const int* pmu = &Mults(f);
pmu -= f;
int Mf = pmu[f];
int Ml = pmu[l];
const int Mf = pmu[0];
const int Ml = pmu[l - f];
if (Mf <= 0)
return 0;
if (Ml <= 0)
return 0;
int sigma;
if (Periodic)
{
if (Mf > Degree)
@@ -384,7 +383,7 @@ int BSplCLib::NbPoles(const int Degree, const bool Periodic, const NCollection_A
}
else
{
int Deg1 = Degree + 1;
const int Deg1 = Degree + 1;
if (Mf > Deg1)
return 0;
if (Ml > Deg1)
@@ -392,7 +391,7 @@ int BSplCLib::NbPoles(const int Degree, const bool Periodic, const NCollection_A
sigma = Mf + Ml - Deg1;
}
for (i = f + 1; i < l; i++)
for (int i = 1; i < l - f; i++)
{
if (pmu[i] <= 0)
return 0;
@@ -3008,39 +3007,45 @@ int BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
return 1;
}
double* PolesArray = &Array;
double* PolesArray = &Array;
const int aLRow = Matrix.LowerRow();
const int aNCols = Matrix.ColNumber();
const double* aMatData = &Matrix(aLRow, 1);
for (int ii = Matrix.LowerRow() + 1; ii <= Matrix.UpperRow(); ii++)
for (int ii = aLRow + 1; ii <= Matrix.UpperRow(); ii++)
{
const int aMinIndex =
(ii - LowerBandWidth >= Matrix.LowerRow() ? ii - LowerBandWidth : Matrix.LowerRow());
const int aMinIndex = (ii - LowerBandWidth >= aLRow ? ii - LowerBandWidth : aLRow);
const double* aRowII = aMatData + (ii - aLRow) * aNCols;
for (int jj = aMinIndex; jj < ii; jj++)
{
const double aCoeff = aRowII[jj - ii + LowerBandWidth];
for (int kk = 0; kk < ArrayDimension; kk++)
{
PolesArray[(ii - 1) * ArrayDimension + kk] +=
PolesArray[(jj - 1) * ArrayDimension + kk] * Matrix(ii, jj - ii + LowerBandWidth + 1);
PolesArray[(jj - 1) * ArrayDimension + kk] * aCoeff;
}
}
}
for (int ii = Matrix.UpperRow(); ii >= Matrix.LowerRow(); ii--)
for (int ii = Matrix.UpperRow(); ii >= aLRow; ii--)
{
const int aMaxIndex =
(ii + UpperBandWidth <= Matrix.UpperRow() ? ii + UpperBandWidth : Matrix.UpperRow());
const double* aRowII = aMatData + (ii - aLRow) * aNCols;
for (int jj = aMaxIndex; jj > ii; jj--)
{
const double aCoeff = aRowII[jj - ii + LowerBandWidth];
for (int kk = 0; kk < ArrayDimension; kk++)
{
PolesArray[(ii - 1) * ArrayDimension + kk] -=
PolesArray[(jj - 1) * ArrayDimension + kk] * Matrix(ii, jj - ii + LowerBandWidth + 1);
PolesArray[(jj - 1) * ArrayDimension + kk] * aCoeff;
}
}
// Fixing a bug PRO18577 to avoid division by zero
const double aDivisor = Matrix(ii, LowerBandWidth + 1);
const double aDivisor = aRowII[LowerBandWidth];
constexpr double THE_TOLERANCE = 1.0e-16;
if (std::abs(aDivisor) <= THE_TOLERANCE)
{
@@ -3147,9 +3152,14 @@ void BSplCLib::Interpolate(const int Degree,
double& Poles,
int& InversionProblem)
{
int ErrorCode, UpperBandWidth, LowerBandWidth;
// double *PolesArray = &Poles ;
math_Matrix InterpolationMatrix(1, Parameters.Length(), 1, 2 * Degree + 1);
int ErrorCode, UpperBandWidth, LowerBandWidth;
constexpr int THE_STACK_LIMIT = 2048;
const int aMatSize = Parameters.Length() * (2 * Degree + 1);
double aStackBuf[THE_STACK_LIMIT];
math_Matrix InterpolationMatrix =
(aMatSize <= THE_STACK_LIMIT)
? math_Matrix(aStackBuf, 1, Parameters.Length(), 1, 2 * Degree + 1)
: math_Matrix(1, Parameters.Length(), 1, 2 * Degree + 1);
ErrorCode = BSplCLib::BuildBSpMatrix(Parameters,
ContactOrderArray,
FlatKnots,
@@ -3187,9 +3197,14 @@ void BSplCLib::Interpolate(const int Degree,
double& Weights,
int& InversionProblem)
{
int ErrorCode, UpperBandWidth, LowerBandWidth;
math_Matrix InterpolationMatrix(1, Parameters.Length(), 1, 2 * Degree + 1);
int ErrorCode, UpperBandWidth, LowerBandWidth;
constexpr int THE_STACK_LIMIT = 2048;
const int aMatSize = Parameters.Length() * (2 * Degree + 1);
double aStackBuf[THE_STACK_LIMIT];
math_Matrix InterpolationMatrix =
(aMatSize <= THE_STACK_LIMIT)
? math_Matrix(aStackBuf, 1, Parameters.Length(), 1, 2 * Degree + 1)
: math_Matrix(1, Parameters.Length(), 1, 2 * Degree + 1);
ErrorCode = BSplCLib::BuildBSpMatrix(Parameters,
ContactOrderArray,
FlatKnots,
src/FoundationClasses/TKMath/BSplCLib/BSplCLib_2.cxx
+60 -43
View File
@@ -316,12 +316,20 @@ int BSplCLib::BuildBSpMatrix(const NCollection_Array1<double>& Parameters,
return 1;
}
math_Matrix aBSplineBasis(1, aMaxOrder, 1, aMaxOrder);
double aBasisBuf[aMaxOrder * aMaxOrder];
math_Matrix aBSplineBasis(aBasisBuf, 1, aMaxOrder, 1, aMaxOrder);
// Zero the entire matrix once instead of per-row zero-fill loops.
Matrix.Init(0.0);
const int aMatLRow = Matrix.LowerRow();
const int aMatNCols = Matrix.ColNumber();
double* aMatData = &Matrix(aMatLRow, 1);
for (int i = Parameters.Lower(); i <= Parameters.Upper(); i++)
{
int aFirstNonZeroIndex = 0;
const int anErrorCode = BSplCLib::EvalBsplineBasis(ContactOrderArray(i),
const int aContactOrder = ContactOrderArray(i);
const int anErrorCode = BSplCLib::EvalBsplineBasis(aContactOrder,
anOrder,
FlatKnots,
Parameters(i),
@@ -332,19 +340,12 @@ int BSplCLib::BuildBSpMatrix(const NCollection_Array1<double>& Parameters,
return 2;
}
int anIndex = LowerBandWidth + 1 + aFirstNonZeroIndex - i;
for (int j = 1; j < anIndex; j++)
int anIndex = LowerBandWidth + 1 + aFirstNonZeroIndex - i;
double* aRowData = aMatData + (i - aMatLRow) * aMatNCols;
const double* aBasisSrc = aBasisBuf + aContactOrder * aMaxOrder;
for (int j = 0; j < anOrder; j++)
{
Matrix.Value(i, j) = 0.0;
}
for (int j = 1; j <= anOrder; j++)
{
Matrix.Value(i, anIndex) = aBSplineBasis(ContactOrderArray(i) + 1, j);
anIndex += 1;
}
for (int j = anIndex; j <= aBandWidth; j++)
{
Matrix.Value(i, j) = 0.0;
aRowData[anIndex + j - 1] = aBasisSrc[j];
}
}
@@ -361,14 +362,20 @@ int BSplCLib::FactorBandedMatrix(math_Matrix& Matrix,
const int aBandWidth = UpperBandWidth + LowerBandWidth + 1;
PivotIndexProblem = 0;
for (int i = Matrix.LowerRow() + 1; i <= Matrix.UpperRow(); i++)
const int aLRow = Matrix.LowerRow();
const int aNCols = Matrix.ColNumber();
double* aData = &Matrix(aLRow, 1);
for (int i = aLRow + 1; i <= Matrix.UpperRow(); i++)
{
const int aMinIndex = (LowerBandWidth - i + 2 >= 1 ? LowerBandWidth - i + 2 : 1);
double* aRowI = aData + (i - aLRow) * aNCols;
for (int j = aMinIndex; j <= LowerBandWidth; j++)
{
const int anIndex = i - LowerBandWidth + j - 1;
const double aPivot = Matrix(anIndex, LowerBandWidth + 1);
const int anIndex = i - LowerBandWidth + j - 1;
const double* aRowIdx = aData + (anIndex - aLRow) * aNCols;
const double aPivot = aRowIdx[LowerBandWidth];
if (std::abs(aPivot) <= RealSmall())
{
PivotIndexProblem = anIndex;
@@ -376,12 +383,12 @@ int BSplCLib::FactorBandedMatrix(math_Matrix& Matrix,
}
const double anInverse = -1.0 / aPivot;
Matrix(i, j) = Matrix(i, j) * anInverse;
aRowI[j - 1] = aRowI[j - 1] * anInverse;
const int aMaxIndex = aBandWidth + anIndex - i;
for (int k = j + 1; k <= aMaxIndex; k++)
{
Matrix(i, k) += Matrix(i, j) * Matrix(anIndex, k + i - anIndex);
aRowI[k - 1] += aRowI[j - 1] * aRowIdx[k + i - anIndex - 1];
}
}
}
@@ -448,8 +455,17 @@ int BSplCLib::EvalBsplineBasis(const int DerivativeReque
FirstNonZeroBsplineIndex = ii - Order + 1;
BsplineBasis(1, 1) = 1.0;
aLocalRequest = DerivativeRequest;
// Use raw pointers for BsplineBasis and FlatKnots to bypass accessor overhead
// (math_Matrix::Value -> math_DoubleTab::Value -> NCollection_Array2::Value ->
// NCollection_Array1::at with bounds checks and DYLD stubs in tight loops).
const int aBasisNCols = BsplineBasis.ColNumber();
double* aBasisData = &BsplineBasis(1, 1);
const double* aKnotsData = &FlatKnots(FlatKnots.Lower());
constexpr double aResolution = gp::Resolution();
ii -= FlatKnots.Lower(); // rebase to zero-based indexing into aKnotsData
aBasisData[0] = 1.0;
aLocalRequest = DerivativeRequest;
if (DerivativeRequest >= Order)
{
aLocalRequest = Order - 1;
@@ -457,21 +473,21 @@ int BSplCLib::EvalBsplineBasis(const int DerivativeReque
for (int qq = 2; qq <= Order - aLocalRequest; qq++)
{
BsplineBasis(1, qq) = 0.0;
aBasisData[qq - 1] = 0.0;
for (int pp = 1; pp <= qq - 1; pp++)
{
const double aScale = FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1);
if (std::abs(aScale) < gp::Resolution())
const double aScale = aKnotsData[ii + pp] - aKnotsData[ii - qq + pp + 1];
if (std::abs(aScale) < aResolution)
{
return 2;
}
const double aFactor = (Parameter - FlatKnots(ii - qq + pp + 1)) / aScale;
const double aSaved = aFactor * BsplineBasis(1, pp);
BsplineBasis(1, pp) *= (1.0 - aFactor);
BsplineBasis(1, pp) += BsplineBasis(1, qq);
BsplineBasis(1, qq) = aSaved;
const double aFactor = (Parameter - aKnotsData[ii - qq + pp + 1]) / aScale;
const double aSaved = aFactor * aBasisData[pp - 1];
aBasisData[pp - 1] *= (1.0 - aFactor);
aBasisData[pp - 1] += aBasisData[qq - 1];
aBasisData[qq - 1] = aSaved;
}
}
@@ -479,37 +495,38 @@ int BSplCLib::EvalBsplineBasis(const int DerivativeReque
{
for (int pp = 1; pp <= qq - 1; pp++)
{
BsplineBasis(Order - qq + 2, pp) = BsplineBasis(1, pp);
aBasisData[(Order - qq + 1) * aBasisNCols + (pp - 1)] = aBasisData[pp - 1];
}
BsplineBasis(1, qq) = 0.0;
aBasisData[qq - 1] = 0.0;
for (int ss = Order - aLocalRequest + 1; ss <= qq; ss++)
{
BsplineBasis(Order - ss + 2, qq) = 0.0;
aBasisData[(Order - ss + 1) * aBasisNCols + (qq - 1)] = 0.0;
}
for (int pp = 1; pp <= qq - 1; pp++)
{
const double aScale = FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1);
if (std::abs(aScale) < gp::Resolution())
const double aScale = aKnotsData[ii + pp] - aKnotsData[ii - qq + pp + 1];
if (std::abs(aScale) < aResolution)
{
return 2;
}
const double anInverse = 1.0 / aScale;
const double aFactor = (Parameter - FlatKnots(ii - qq + pp + 1)) * anInverse;
double aSaved = aFactor * BsplineBasis(1, pp);
BsplineBasis(1, pp) *= (1.0 - aFactor);
BsplineBasis(1, pp) += BsplineBasis(1, qq);
BsplineBasis(1, qq) = aSaved;
const double aFactor = (Parameter - aKnotsData[ii - qq + pp + 1]) * anInverse;
double aSaved = aFactor * aBasisData[pp - 1];
aBasisData[pp - 1] *= (1.0 - aFactor);
aBasisData[pp - 1] += aBasisData[qq - 1];
aBasisData[qq - 1] = aSaved;
const double aLocalInverse = static_cast<double>(qq - 1) * anInverse;
for (int ss = Order - aLocalRequest + 1; ss <= qq; ss++)
{
aSaved = aLocalInverse * BsplineBasis(Order - ss + 2, pp);
BsplineBasis(Order - ss + 2, pp) *= -aLocalInverse;
BsplineBasis(Order - ss + 2, pp) += BsplineBasis(Order - ss + 2, qq);
BsplineBasis(Order - ss + 2, qq) = aSaved;
double* aRowS = aBasisData + (Order - ss + 1) * aBasisNCols;
aSaved = aLocalInverse * aRowS[pp - 1];
aRowS[pp - 1] *= -aLocalInverse;
aRowS[pp - 1] += aRowS[qq - 1];
aRowS[qq - 1] = aSaved;
}
}
}
src/FoundationClasses/TKMath/Convert/Convert_CompPolynomialToPoles.cxx
+5 -2
View File
@@ -197,7 +197,10 @@ void Convert_CompPolynomialToPoles::Perform(const int Nu
myFlatKnots = NCollection_Array1<double>(1, num_flat_knots);
BSplCLib::KnotSequence(myKnots, myMults, myDegree, false, myFlatKnots);
NCollection_Array1<double> parameters(1, num_poles);
constexpr int THE_MAX_POLES = 128;
double aParamBuf[THE_MAX_POLES];
int aContactBuf[THE_MAX_POLES];
NCollection_Array1<double> parameters(aParamBuf[0], 1, num_poles, num_poles <= THE_MAX_POLES);
BSplCLib::BuildSchoenbergPoints(myDegree, myFlatKnots, parameters);
myPoles = NCollection_Array2<double>(1, num_poles, 1, Dimension);
index = 2;
@@ -205,7 +208,7 @@ void Convert_CompPolynomialToPoles::Perform(const int Nu
Pindex = PolynomialIntervals.LowerRow();
poles_array = (double*)&myPoles.ChangeValue(1, 1);
NCollection_Array1<int> contact_array(1, num_poles);
NCollection_Array1<int> contact_array(aContactBuf[0], 1, num_poles, num_poles <= THE_MAX_POLES);
poles_index = 0;
for (ii = 1; ii <= num_poles; ii++, poles_index += Dimension)
src/FoundationClasses/TKMath/math/math_VectorBase.lxx
+81 -69
View File
@@ -101,66 +101,60 @@ void math_VectorBase<TheItemType>::SetLower(const int theLower)
template <typename TheItemType>
double math_VectorBase<TheItemType>::Norm() const
{
// 4-way unrolled accumulation for better vectorization
double aSum1 = 0.0, aSum2 = 0.0, aSum3 = 0.0, aSum4 = 0.0;
int anIndex = Lower();
int anUpper = Upper();
int anUpper4 = anUpper - 3;
const int aLen = Length();
const TheItemType* aPtr = &Array(Lower());
double aSum1 = 0.0, aSum2 = 0.0, aSum3 = 0.0, aSum4 = 0.0;
int i = 0;
int aLen4 = aLen - 3;
// Process 4 elements at a time
for (; anIndex <= anUpper4; anIndex += 4)
for (; i < aLen4; i += 4)
{
const double aVal0 = static_cast<double>(Array(anIndex));
const double aVal1 = static_cast<double>(Array(anIndex + 1));
const double aVal2 = static_cast<double>(Array(anIndex + 2));
const double aVal3 = static_cast<double>(Array(anIndex + 3));
const double aVal0 = static_cast<double>(aPtr[i]);
const double aVal1 = static_cast<double>(aPtr[i + 1]);
const double aVal2 = static_cast<double>(aPtr[i + 2]);
const double aVal3 = static_cast<double>(aPtr[i + 3]);
aSum1 += aVal0 * aVal0;
aSum2 += aVal1 * aVal1;
aSum3 += aVal2 * aVal2;
aSum4 += aVal3 * aVal3;
}
// Process remaining elements
for (; anIndex <= anUpper; ++anIndex)
for (; i < aLen; ++i)
{
const double aVal = static_cast<double>(Array(anIndex));
const double aVal = static_cast<double>(aPtr[i]);
aSum1 += aVal * aVal;
}
// Combine partial sums (pairwise for better numerical stability)
return std::sqrt((aSum1 + aSum2) + (aSum3 + aSum4));
}
template <typename TheItemType>
double math_VectorBase<TheItemType>::Norm2() const
{
// 4-way unrolled accumulation for better vectorization
double aSum1 = 0.0, aSum2 = 0.0, aSum3 = 0.0, aSum4 = 0.0;
int anIndex = Lower();
int anUpper = Upper();
int anUpper4 = anUpper - 3;
const int aLen = Length();
const TheItemType* aPtr = &Array(Lower());
double aSum1 = 0.0, aSum2 = 0.0, aSum3 = 0.0, aSum4 = 0.0;
int i = 0;
int aLen4 = aLen - 3;
// Process 4 elements at a time
for (; anIndex <= anUpper4; anIndex += 4)
for (; i < aLen4; i += 4)
{
const double aVal0 = static_cast<double>(Array(anIndex));
const double aVal1 = static_cast<double>(Array(anIndex + 1));
const double aVal2 = static_cast<double>(Array(anIndex + 2));
const double aVal3 = static_cast<double>(Array(anIndex + 3));
const double aVal0 = static_cast<double>(aPtr[i]);
const double aVal1 = static_cast<double>(aPtr[i + 1]);
const double aVal2 = static_cast<double>(aPtr[i + 2]);
const double aVal3 = static_cast<double>(aPtr[i + 3]);
aSum1 += aVal0 * aVal0;
aSum2 += aVal1 * aVal1;
aSum3 += aVal2 * aVal2;
aSum4 += aVal3 * aVal3;
}
// Process remaining elements
for (; anIndex <= anUpper; ++anIndex)
for (; i < aLen; ++i)
{
const double aVal = static_cast<double>(Array(anIndex));
const double aVal = static_cast<double>(aPtr[i]);
aSum1 += aVal * aVal;
}
// Combine partial sums (pairwise for better numerical stability)
return (aSum1 + aSum2) + (aSum3 + aSum4);
}
@@ -446,17 +440,22 @@ void math_VectorBase<TheItemType>::Multiply(const math_Matrix&
(Length() != theLeft.RowNumber()) || (theLeft.ColNumber() != theRight.Length()),
"math_VectorBase::Multiply() - input matrix and /or vector have wrong dimensions");
int Index = Lower();
for (int I = theLeft.LowerRow(); I <= theLeft.UpperRow(); I++)
// result[r] = sum_c theLeft(r, c) * theRight[c]
const int aNRows = theLeft.RowNumber();
const int aNCols = theLeft.ColNumber();
const double* aMatData = &theLeft(theLeft.LowerRow(), theLeft.LowerCol());
const TheItemType* aRightPtr = &theRight.Array(theRight.Lower());
TheItemType* aResPtr = &Array(Lower());
for (int r = 0; r < aNRows; r++)
{
Array(Index) = 0.0;
int K = theRight.Lower();
for (int J = theLeft.LowerCol(); J <= theLeft.UpperCol(); J++)
TheItemType aSum = 0.0;
const double* aRowData = aMatData + r * aNCols;
for (int c = 0; c < aNCols; c++)
{
Array(Index) = Array(Index) + theLeft(I, J) * theRight.Array(K);
K++;
aSum += aRowData[c] * aRightPtr[c];
}
Index++;
aResPtr[r] = aSum;
}
}
@@ -468,17 +467,21 @@ void math_VectorBase<TheItemType>::Multiply(const math_VectorBase<TheItemType>&
(Length() != theRight.ColNumber()) || (theLeft.Length() != theRight.RowNumber()),
"math_VectorBase::Multiply() - input matrix and /or vector have wrong dimensions");
int Index = Lower();
for (int J = theRight.LowerCol(); J <= theRight.UpperCol(); J++)
// result[c] = sum_r theLeft[r] * theRight(r, c)
const int aNCols = theRight.ColNumber();
const int aNRows = theRight.RowNumber();
const double* aMatData = &theRight(theRight.LowerRow(), theRight.LowerCol());
const TheItemType* aLeftPtr = &theLeft.Array(theLeft.Lower());
TheItemType* aResPtr = &Array(Lower());
for (int c = 0; c < aNCols; c++)
{
Array(Index) = 0.0;
int K = theLeft.Lower();
for (int I = theRight.LowerRow(); I <= theRight.UpperRow(); I++)
TheItemType aSum = 0.0;
for (int r = 0; r < aNRows; r++)
{
Array(Index) = Array(Index) + theLeft.Array(K) * theRight(I, J);
K++;
aSum += aLeftPtr[r] * aMatData[r * aNCols + c];
}
Index++;
aResPtr[c] = aSum;
}
}
@@ -490,17 +493,21 @@ void math_VectorBase<TheItemType>::TMultiply(const math_Matrix&
(Length() != theTLeft.ColNumber()) || (theTLeft.RowNumber() != theRight.Length()),
"math_VectorBase::TMultiply() - input matrix and /or vector have wrong dimensions");
int Index = Lower();
for (int I = theTLeft.LowerCol(); I <= theTLeft.UpperCol(); I++)
// result[c] = sum_r theTLeft(r, c) * theRight[r] (transpose-multiply)
const int aNCols = theTLeft.ColNumber();
const int aNRows = theTLeft.RowNumber();
const double* aMatData = &theTLeft(theTLeft.LowerRow(), theTLeft.LowerCol());
const TheItemType* aRightPtr = &theRight.Array(theRight.Lower());
TheItemType* aResPtr = &Array(Lower());
for (int c = 0; c < aNCols; c++)
{
Array(Index) = 0.0;
int K = theRight.Lower();
for (int J = theTLeft.LowerRow(); J <= theTLeft.UpperRow(); J++)
TheItemType aSum = 0.0;
for (int r = 0; r < aNRows; r++)
{
Array(Index) = Array(Index) + theTLeft(J, I) * theRight.Array(K);
K++;
aSum += aMatData[r * aNCols + c] * aRightPtr[r];
}
Index++;
aResPtr[c] = aSum;
}
}
@@ -512,17 +519,22 @@ void math_VectorBase<TheItemType>::TMultiply(const math_VectorBase<TheItemType>&
(Length() != theTRight.RowNumber()) || (theLeft.Length() != theTRight.ColNumber()),
"math_VectorBase::TMultiply() - input matrix and /or vector have wrong dimensions");
int Index = Lower();
for (int J = theTRight.LowerRow(); J <= theTRight.UpperRow(); J++)
// result[r] = sum_c theLeft[c] * theTRight(r, c)
const int aNCols = theTRight.ColNumber();
const int aNRows = theTRight.RowNumber();
const double* aMatData = &theTRight(theTRight.LowerRow(), theTRight.LowerCol());
const TheItemType* aLeftPtr = &theLeft.Array(theLeft.Lower());
TheItemType* aResPtr = &Array(Lower());
for (int r = 0; r < aNRows; r++)
{
Array(Index) = 0.0;
int K = theLeft.Lower();
for (int I = theTRight.LowerCol(); I <= theTRight.UpperCol(); I++)
TheItemType aSum = 0.0;
const double* aRowData = aMatData + r * aNCols;
for (int c = 0; c < aNCols; c++)
{
Array(Index) = Array(Index) + theLeft.Array(K) * theTRight(J, I);
K++;
aSum += aLeftPtr[c] * aRowData[c];
}
Index++;
aResPtr[r] = aSum;
}
}
@@ -530,19 +542,19 @@ template <typename TheItemType>
TheItemType math_VectorBase<TheItemType>::Multiplied(
const math_VectorBase<TheItemType>& theRight) const
{
double Result = 0;
Standard_DimensionError_Raise_if(
Length() != theRight.Length(),
"math_VectorBase::Multiplied() - input vector has wrong dimensions");
int I = theRight.Lower();
for (int Index = Lower(); Index <= Upper(); Index++)
const int aLen = Length();
const TheItemType* aLeftPtr = &Array(Lower());
const TheItemType* aRightPtr = &theRight.Array(theRight.Lower());
TheItemType aResult = 0;
for (int i = 0; i < aLen; i++)
{
Result = Result + Array(Index) * theRight.Array(I);
I++;
aResult += aLeftPtr[i] * aRightPtr[i];
}
return Result;
return aResult;
}
template <typename TheItemType>
src/ModelingAlgorithms/TKFillet/BRepBlend/BRepBlend_AppFuncRoot.cxx
+11 -7
View File
@@ -75,8 +75,14 @@ BRepBlend_AppFuncRoot::BRepBlend_AppFuncRoot(occ::handle<BRepBlend_Line>& Line,
{
myBary.SetCoord(0, 0, 0);
}
mySolver = std::make_unique<math_FunctionSetRoot>(Func, myTolerance, 30);
}
//==================================================================================================
BRepBlend_AppFuncRoot::~BRepBlend_AppFuncRoot() = default;
//================================================================================
// Function: D0
// Purpose : Calculation of section for v = Param, if calculation fails
@@ -337,28 +343,26 @@ bool BRepBlend_AppFuncRoot::SearchPoint(Blend_AppFunction& Func,
// (2) Calculation of the solution ------------------------
Func.Set(Param);
Func.GetBounds(X1, X2);
math_FunctionSetRoot rsnld(Func, myTolerance, 30);
mySolver->Perform(Func, XInit, X1, X2);
rsnld.Perform(Func, XInit, X1, X2);
if (!rsnld.IsDone())
if (!mySolver->IsDone())
{
#ifdef BREPBLEND_DEB
std::cout << "AppFunc : RNLD Not done en t = " << Param << std::endl;
#endif
return false;
}
rsnld.Root(Sol);
mySolver->Root(Sol);
// (3) Storage of the point
Point(Func, Param, Sol, Pnt);
// (4) Insertion of the point if the calculation seems long.
if ((!Trouve) && (rsnld.NbIterations() > 3))
if ((!Trouve) && (mySolver->NbIterations() > 3))
{
#ifdef OCCT_DEBUG
std::cout << "Evaluation in t = " << Param << "given" << std::endl;
rsnld.Dump(std::cout);
mySolver->Dump(std::cout);
#endif
myLine->InsertBefore(Index + 1, Pnt);
}
src/ModelingAlgorithms/TKFillet/BRepBlend/BRepBlend_AppFuncRoot.hxx
+16 -9
View File
@@ -30,8 +30,12 @@
#include <gp_Vec2d.hxx>
#include <Standard_Integer.hxx>
#include <GeomAbs_Shape.hxx>
#include <memory>
class BRepBlend_Line;
class Blend_AppFunction;
class math_FunctionSetRoot;
//! Function to approximate by AppSurface
class BRepBlend_AppFuncRoot : public Approx_SweepFunction
@@ -152,6 +156,8 @@ public:
DEFINE_STANDARD_RTTIEXT(BRepBlend_AppFuncRoot, Approx_SweepFunction)
Standard_EXPORT ~BRepBlend_AppFuncRoot();
protected:
Standard_EXPORT BRepBlend_AppFuncRoot(occ::handle<BRepBlend_Line>& Line,
Blend_AppFunction& Func,
@@ -166,15 +172,16 @@ private:
const int LastIndex,
int& ParamIndex) const;
occ::handle<BRepBlend_Line> myLine;
void* myFunc;
math_Vector myTolerance;
Blend_Point myPnt;
gp_Pnt myBary;
math_Vector X1;
math_Vector X2;
math_Vector XInit;
math_Vector Sol;
occ::handle<BRepBlend_Line> myLine;
void* myFunc;
math_Vector myTolerance;
Blend_Point myPnt;
gp_Pnt myBary;
math_Vector X1;
math_Vector X2;
math_Vector XInit;
math_Vector Sol;
std::unique_ptr<math_FunctionSetRoot> mySolver;
};
#endif // _BRepBlend_AppFuncRoot_HeaderFile
src/ModelingAlgorithms/TKGeomAlgo/GeomFill/GeomFill_PolynomialConvertor.cxx
+41 -57
View File
@@ -23,9 +23,7 @@
#include <Standard_Integer.hxx>
#include <StdFail_NotDone.hxx>
#include <NCollection_Array1.hxx>
#include <NCollection_HArray1.hxx>
#include <NCollection_Array2.hxx>
#include <NCollection_HArray2.hxx>
#include <math_Vector.hxx>
GeomFill_PolynomialConvertor::GeomFill_PolynomialConvertor()
@@ -43,63 +41,49 @@ bool GeomFill_PolynomialConvertor::Initialized() const
void GeomFill_PolynomialConvertor::Init()
{
if (myinit)
return; // On n'initialise qu'une fois
int ii, jj;
double terme;
math_Matrix H(1, Ordre, 1, Ordre), B(1, Ordre, 1, Ordre);
occ::handle<NCollection_HArray1<double>> Coeffs =
new (NCollection_HArray1<double>)(1, Ordre * Ordre),
TrueInter = new (NCollection_HArray1<double>)(1, 2);
occ::handle<NCollection_HArray2<double>> Poles1d =
new (NCollection_HArray2<double>)(1, Ordre, 1, Ordre),
Inter = new (NCollection_HArray2<double>)(1, 1, 1, 2);
// Calcul de B
Inter->SetValue(1, 1, -1);
Inter->SetValue(1, 2, 1);
TrueInter->SetValue(1, -1);
TrueInter->SetValue(2, 1);
Coeffs->Init(0);
for (ii = 1; ii <= Ordre; ii++)
{
Coeffs->SetValue(ii + (ii - 1) * Ordre, 1);
}
// Convertion ancienne formules
occ::handle<NCollection_HArray1<int>> Ncf = new (NCollection_HArray1<int>)(1, 1);
Ncf->Init(Ordre);
Convert_CompPolynomialToPoles AConverter(1, 1, 8, 8, Ncf, Coeffs, Inter, TrueInter);
/* Convert_CompPolynomialToPoles
AConverter(8, Ordre-1, Ordre-1,
Coeffs,
Inter,
TrueInter); En attente du bon Geomlite*/
Poles1d = new NCollection_HArray2<double>(AConverter.Poles());
for (jj = 1; jj <= Ordre; jj++)
{
for (ii = 1; ii <= Ordre; ii++)
{
terme = Poles1d->Value(ii, jj);
if (std::abs(terme - 1) < 1.e-9)
terme = 1; // petite retouche
if (std::abs(terme + 1) < 1.e-9)
terme = -1;
B(ii, jj) = terme;
}
}
// Calcul de H
myinit = PLib::HermiteCoefficients(-1, 1, Ordre / 2 - 1, Ordre / 2 - 1, H);
H.Transpose();
if (!myinit)
return;
// reste l'essentiel
BH = B * H;
// BH = B * H where B is the monomial-to-BSpline conversion matrix on [-1,1], degree 7,
// and H is the Hermite coefficients matrix. Both are mathematical constants computed once.
static const math_Matrix THE_BH_MATRIX = []() {
constexpr int anOrdre = 8;
NCollection_Array1<double> aCoeffs(1, anOrdre * anOrdre);
NCollection_Array1<double> anInter(1, 2), aTrueInter(1, 2);
anInter.SetValue(1, -1);
anInter.SetValue(2, 1);
aTrueInter.SetValue(1, -1);
aTrueInter.SetValue(2, 1);
aCoeffs.Init(0);
for (int ii = 1; ii <= anOrdre; ii++)
aCoeffs.SetValue(ii + (ii - 1) * anOrdre, 1);
Convert_CompPolynomialToPoles aConverter(anOrdre,
anOrdre - 1,
anOrdre - 1,
aCoeffs,
anInter,
aTrueInter);
const NCollection_Array2<double>& aPoles = aConverter.Poles();
math_Matrix aB(1, anOrdre, 1, anOrdre);
for (int jj = 1; jj <= anOrdre; jj++)
for (int ii = 1; ii <= anOrdre; ii++)
{
double aTerm = aPoles.Value(ii, jj);
if (std::abs(aTerm - 1) < 1.e-9)
aTerm = 1;
if (std::abs(aTerm + 1) < 1.e-9)
aTerm = -1;
aB(ii, jj) = aTerm;
}
math_Matrix aH(1, anOrdre, 1, anOrdre);
PLib::HermiteCoefficients(-1, 1, anOrdre / 2 - 1, anOrdre / 2 - 1, aH);
aH.Transpose();
return math_Matrix(aB * aH);
}();
BH = THE_BH_MATRIX;
myinit = true;
}
void GeomFill_PolynomialConvertor::Section(const gp_Pnt& FirstPnt,
src/ModelingAlgorithms/TKGeomAlgo/GeomFill/GeomFill_QuasiAngularConvertor.cxx
+35 -34
View File
@@ -23,7 +23,6 @@
#include <StdFail_NotDone.hxx>
#include <NCollection_Array1.hxx>
#include <NCollection_Array2.hxx>
#include <NCollection_HArray2.hxx>
#define NullAngle 1.e-6
@@ -72,43 +71,45 @@ bool GeomFill_QuasiAngularConvertor::Initialized() const
void GeomFill_QuasiAngularConvertor::Init()
{
if (myinit)
return; // On n'initialise qu'une fois
int ii, jj, Ordre = 7;
double terme;
NCollection_Array1<double> Coeffs(1, Ordre * Ordre), TrueInter(1, 2), Inter(1, 2);
occ::handle<NCollection_HArray2<double>> Poles1d =
new (NCollection_HArray2<double>)(1, Ordre, 1, Ordre);
return;
// Calcul de B
Inter.SetValue(1, -1);
Inter.SetValue(2, 1);
TrueInter.SetValue(1, -1);
TrueInter.SetValue(2, 1);
// B is the monomial-to-BSpline conversion matrix on [-1,1], degree 6.
// It is a mathematical constant computed once.
static const math_Matrix THE_BASIS = []() {
constexpr int anOrdre = 7;
NCollection_Array1<double> aCoeffs(1, anOrdre * anOrdre);
NCollection_Array1<double> anInter(1, 2), aTrueInter(1, 2);
anInter.SetValue(1, -1);
anInter.SetValue(2, 1);
aTrueInter.SetValue(1, -1);
aTrueInter.SetValue(2, 1);
aCoeffs.Init(0);
for (int ii = 1; ii <= anOrdre; ii++)
aCoeffs.SetValue(ii + (ii - 1) * anOrdre, 1);
Coeffs.Init(0);
for (ii = 1; ii <= Ordre; ii++)
{
Coeffs.SetValue(ii + (ii - 1) * Ordre, 1);
}
Convert_CompPolynomialToPoles aConverter(anOrdre,
anOrdre - 1,
anOrdre - 1,
aCoeffs,
anInter,
aTrueInter);
const NCollection_Array2<double>& aPoles = aConverter.Poles();
math_Matrix aResult(1, anOrdre, 1, anOrdre);
for (int jj = 1; jj <= anOrdre; jj++)
for (int ii = 1; ii <= anOrdre; ii++)
{
double aTerm = aPoles.Value(ii, jj);
if (std::abs(aTerm - 1) < 1.e-9)
aTerm = 1;
if (std::abs(aTerm + 1) < 1.e-9)
aTerm = -1;
aResult(ii, jj) = aTerm;
}
return aResult;
}();
// Convertion
Convert_CompPolynomialToPoles AConverter(Ordre, Ordre - 1, Ordre - 1, Coeffs, Inter, TrueInter);
Poles1d = new NCollection_HArray2<double>(AConverter.Poles());
B = THE_BASIS;
for (jj = 1; jj <= Ordre; jj++)
{
for (ii = 1; ii <= Ordre; ii++)
{
terme = Poles1d->Value(ii, jj);
if (std::abs(terme - 1) < 1.e-9)
terme = 1; // petite retouche
if (std::abs(terme + 1) < 1.e-9)
terme = -1;
B(ii, jj) = terme;
}
}
// Init des polynomes
Vx.Init(0);
Vx(1) = 1;
Vy.Init(0);
src/ModelingData/TKGeomBase/Approx/Approx_Curve2d.cxx
+5 -4
View File
@@ -149,12 +149,13 @@ Approx_Curve2d::Approx_Curve2d(const occ::handle<Adaptor2d_Curve2d>& C2D,
if (myHasResult)
{
NCollection_Array1<gp_Pnt2d> Poles2d(1, aApprox.NbPoles());
NCollection_Array1<double> Poles1dU(1, aApprox.NbPoles());
const int aNbPoles = aApprox.NbPoles();
NCollection_Array1<gp_Pnt2d> Poles2d(1, aNbPoles);
NCollection_Array1<double> Poles1dU(1, aNbPoles);
aApprox.Poles1d(1, Poles1dU);
NCollection_Array1<double> Poles1dV(1, aApprox.NbPoles());
NCollection_Array1<double> Poles1dV(1, aNbPoles);
aApprox.Poles1d(2, Poles1dV);
for (int i = 1; i <= aApprox.NbPoles(); i++)
for (int i = 1; i <= aNbPoles; i++)
Poles2d.SetValue(i, gp_Pnt2d(Poles1dU.Value(i), Poles1dV.Value(i)));
occ::handle<NCollection_HArray1<double>> Knots = aApprox.Knots();
src/ModelingData/TKGeomBase/Approx/Approx_CurveOnSurface.cxx
+6 -5
View File
@@ -502,21 +502,22 @@ void Approx_CurveOnSurface::Perform(const int theMaxSegments,
occ::handle<NCollection_HArray1<int>> Mults = aApprox.Multiplicities();
int Degree = aApprox.Degree();
const int aNbPoles = aApprox.NbPoles();
if (!theOnly2d)
{
NCollection_Array1<gp_Pnt> Poles(1, aApprox.NbPoles());
NCollection_Array1<gp_Pnt> Poles(1, aNbPoles);
aApprox.Poles(1, Poles);
myCurve3d = new Geom_BSplineCurve(Poles, Knots->Array1(), Mults->Array1(), Degree);
myError3d = aApprox.MaxError(3, 1);
}
if (!theOnly3d)
{
NCollection_Array1<gp_Pnt2d> Poles2d(1, aApprox.NbPoles());
NCollection_Array1<double> Poles1dU(1, aApprox.NbPoles());
NCollection_Array1<gp_Pnt2d> Poles2d(1, aNbPoles);
NCollection_Array1<double> Poles1dU(1, aNbPoles);
aApprox.Poles1d(1, Poles1dU);
NCollection_Array1<double> Poles1dV(1, aApprox.NbPoles());
NCollection_Array1<double> Poles1dV(1, aNbPoles);
aApprox.Poles1d(2, Poles1dV);
for (int i = 1; i <= aApprox.NbPoles(); i++)
for (int i = 1; i <= aNbPoles; i++)
Poles2d.SetValue(i, gp_Pnt2d(Poles1dU.Value(i), Poles1dV.Value(i)));
myCurve2d = new Geom2d_BSplineCurve(Poles2d, Knots->Array1(), Mults->Array1(), Degree);
src/ModelingData/TKGeomBase/Approx/Approx_SweepApproximation.cxx
+8 -7
View File
@@ -303,12 +303,13 @@ void Approx_SweepApproximation::Approximation(
// --> Fill Champs of the surface ----
int ii, jj;
vdeg = Approx.Degree();
vdeg = Approx.Degree();
const int aNbPoles = Approx.NbPoles();
// Unfortunately Adv_Approx stores the transposition of the required
// so, writing tabPoles = Approx.Poles() will give an erroneous result
// It is only possible to allocate and recopy term by term...
tabPoles = new (NCollection_HArray2<gp_Pnt>)(1, Num3DSS, 1, Approx.NbPoles());
tabWeights = new (NCollection_HArray2<double>)(1, Num3DSS, 1, Approx.NbPoles());
tabPoles = new (NCollection_HArray2<gp_Pnt>)(1, Num3DSS, 1, aNbPoles);
tabWeights = new (NCollection_HArray2<double>)(1, Num3DSS, 1, aNbPoles);
if (Num1DSS == Num3DSS)
{
@@ -316,7 +317,7 @@ void Approx_SweepApproximation::Approximation(
gp_Pnt P;
for (ii = 1; ii <= Num3DSS; ii++)
{
for (jj = 1; jj <= Approx.NbPoles(); jj++)
for (jj = 1; jj <= aNbPoles; jj++)
{
P = Approx.Poles()->Value(jj, ii);
wpoid = Approx.Poles1d()->Value(jj, ii);
@@ -332,7 +333,7 @@ void Approx_SweepApproximation::Approximation(
tabWeights->Init(1);
for (ii = 1; ii <= Num3DSS; ii++)
{
for (jj = 1; jj <= Approx.NbPoles(); jj++)
for (jj = 1; jj <= aNbPoles; jj++)
{
tabPoles->SetValue(ii, jj, Approx.Poles()->Value(jj, ii));
}
@@ -355,10 +356,10 @@ void Approx_SweepApproximation::Approximation(
{
TrsfInv = AAffin->Value(ii).Inverted();
occ::handle<NCollection_HArray1<gp_Pnt2d>> P2d =
new (NCollection_HArray1<gp_Pnt2d>)(1, Approx.NbPoles());
new (NCollection_HArray1<gp_Pnt2d>)(1, aNbPoles);
Approx.Poles2d(ii, P2d->ChangeArray1());
// do not forget to apply inverted homothety.
for (jj = 1; jj <= Approx.NbPoles(); jj++)
for (jj = 1; jj <= aNbPoles; jj++)
{
TrsfInv.Transforms(P2d->ChangeValue(jj).ChangeCoord());
}