From 22268fd999163be0f444186ae68303084b5a041d Mon Sep 17 00:00:00 2001 From: azv Date: Thu, 17 Oct 2013 18:21:03 +0400 Subject: [PATCH] 0024238: Some changes in B-spline surfaces to work with aligned memory, which potentially accelerate calculations with SSE vectorization --- src/BSplSLib/BSplSLib.cdl | 60 +- src/BSplSLib/BSplSLib.cxx | 1062 +++++-------- src/BSplSLib/BSplSLib.lxx | 61 +- src/BSplSLib/BSplSLib_BzSyntaxes.cxx | 84 +- src/Geom/FILES | 3 +- src/Geom/Geom.cdl | 3 +- src/Geom/Geom_BSplineCache.cdl | 108 ++ src/Geom/Geom_BSplineCache.cxx | 24 + src/Geom/Geom_BSplineCache.lxx | 61 + src/Geom/Geom_BSplineSurface.cdl | 17 +- src/Geom/Geom_BSplineSurface.cxx | 217 ++- src/Geom/Geom_BSplineSurface_1.cxx | 192 +-- src/Geom/Geom_BezierSurface.cdl | 14 +- src/Geom/Geom_BezierSurface.cxx | 907 +++++------ src/Geom/Geom_OsculatingSurface.cxx | 375 ++--- src/Geom/Handle_Geom_BSplineCache.hxx | 21 + src/NCollection/FILES | 1 + .../NCollection_AlignAllocator.cxx | 78 + .../NCollection_AlignAllocator.hxx | 81 + src/PLib/PLib.cxx | 1340 ++++------------- 20 files changed, 2023 insertions(+), 2686 deletions(-) create mode 100644 src/Geom/Geom_BSplineCache.cdl create mode 100644 src/Geom/Geom_BSplineCache.cxx create mode 100644 src/Geom/Geom_BSplineCache.lxx create mode 100644 src/Geom/Handle_Geom_BSplineCache.hxx create mode 100644 src/NCollection/NCollection_AlignAllocator.cxx create mode 100644 src/NCollection/NCollection_AlignAllocator.hxx diff --git a/src/BSplSLib/BSplSLib.cdl b/src/BSplSLib/BSplSLib.cdl index 4c07e52312..d219666b31 100755 --- a/src/BSplSLib/BSplSLib.cdl +++ b/src/BSplSLib/BSplSLib.cdl @@ -463,8 +463,8 @@ is UFlatKnots,VFlatKnots : Array1OfReal from TColStd ; Poles : Array2OfPnt from TColgp; Weights : Array2OfReal from TColStd ; - CachePoles : in out Array2OfPnt from TColgp; - CacheWeights : in out Array2OfReal from TColStd); + CacheRational : in out Boolean; + CachePolesWeights : in out Array2OfReal from TColStd); ---Purpose: Perform the evaluation of the Taylor expansion -- of the Bspline normalized between 0 and 1. @@ -477,8 +477,8 @@ is UDegree,VDegree : Integer; UCacheParameter,VCacheParameter : Real; USpanLenght,VSpanLength : Real; - Poles : Array2OfPnt from TColgp ; - Weights : Array2OfReal from TColStd ; + PolesWeightsArray : Array2OfReal from TColStd; + CacheRational : Boolean; Point : out Pnt from gp) ; ---Purpose: Perform the evaluation of the of the cache @@ -494,13 +494,15 @@ is -- effects -- - CoefsD0(U,V : Real; - Poles : Array2OfPnt from TColgp ; - Weights : Array2OfReal from TColStd ; - Point : out Pnt from gp) ; + CoefsD0(U,V : Real; + PolesWeightsArray : Array2OfReal from TColStd; + CacheRational : Boolean; + Point : out Pnt from gp) ; ---Purpose: Calls CacheD0 for Bezier Surfaces Arrays computed with -- the method PolesCoefficients. -- Warning: To be used for BezierSurfaces ONLY!!! + -- Warning: Number of columns of the PolesWeightsArray should be + -- divisible to 4 (three columns for Pole and fourth for Weight) ---C++: inline @@ -508,8 +510,8 @@ is UDegree,VDegree : Integer; UCacheParameter,VCacheParameter : Real; USpanLenght,VSpanLength : Real; - Poles : Array2OfPnt from TColgp ; - Weights : Array2OfReal from TColStd ; + PolesWeightsArray : Array2OfReal from TColStd; + CacheRational : Boolean; Point : out Pnt from gp; VecU, VecV : out Vec from gp) ; @@ -526,14 +528,16 @@ is -- effects -- - CoefsD1(U,V : Real; - Poles : Array2OfPnt from TColgp; - Weights : Array2OfReal from TColStd; - Point : out Pnt from gp; - VecU, VecV : out Vec from gp) ; + CoefsD1(U,V : Real; + PolesWeightsArray : Array2OfReal from TColStd; + CacheRational : Boolean; + Point : out Pnt from gp; + VecU, VecV : out Vec from gp) ; ---Purpose: Calls CacheD0 for Bezier Surfaces Arrays computed with -- the method PolesCoefficients. - -- Warning: To be used for BezierSurfaces ONLY!!! + -- Warning: To be used for BezierSurfaces ONLY!!! + -- Warning: Number of columns of the PolesWeightsArray should be + -- divisible to 4 (three columns for Pole and fourth for Weight) ---C++: inline @@ -541,8 +545,8 @@ is UDegree,VDegree : Integer; UCacheParameter,VCacheParameter : Real; USpanLenght,VSpanLength : Real; - Poles : Array2OfPnt from TColgp ; - Weights : Array2OfReal from TColStd ; + PolesWeightsArray : Array2OfReal from TColStd; + CacheRational : Boolean; Point : out Pnt from gp; VecU, VecV, VecUU, VecUV, VecVV : out Vec from gp) ; @@ -560,25 +564,27 @@ is -- CoefsD2(U,V : Real; - Poles : Array2OfPnt from TColgp ; - Weights : Array2OfReal from TColStd ; + PolesWeightsArray : Array2OfReal from TColStd; + CacheRational : Boolean; Point : out Pnt from gp; VecU, VecV, VecUU, VecUV, VecVV : out Vec from gp) ; ---Purpose: Calls CacheD0 for Bezier Surfaces Arrays computed with -- the method PolesCoefficients. - -- Warning: To be used for BezierSurfaces ONLY!!! + -- Warning: To be used for BezierSurfaces ONLY!!! + -- Warning: Number of columns of the PolesWeightsArray should be + -- divisible to 4 (three columns for Pole and fourth for Weight) ---C++: inline - PolesCoefficients(Poles : Array2OfPnt from TColgp; - CachePoles : in out Array2OfPnt from TColgp); + PolesCoefficients(Poles : Array2OfPnt from TColgp; + CachePolesWeights : in out Array2OfReal from TColStd); ---Purpose: Warning! To be used for BezierSurfaces ONLY!!! ---C++: inline - PolesCoefficients(Poles : Array2OfPnt from TColgp; - Weights : Array2OfReal from TColStd ; - CachePoles : in out Array2OfPnt from TColgp; - CacheWeights : in out Array2OfReal from TColStd) ; + PolesCoefficients(Poles : Array2OfPnt from TColgp; + Weights : Array2OfReal from TColStd ; + CacheRational : in out Boolean; + CachePolesWeights : in out Array2OfReal from TColStd) ; ---Purpose: Encapsulation of BuildCache to perform the -- evaluation of the Taylor expansion for beziersurfaces diff --git a/src/BSplSLib/BSplSLib.cxx b/src/BSplSLib/BSplSLib.cxx index e3d4127378..1f13f3e1da 100755 --- a/src/BSplSLib/BSplSLib.cxx +++ b/src/BSplSLib/BSplSLib.cxx @@ -36,6 +36,7 @@ #include #include #include +#include // for null derivatives static Standard_Real BSplSLib_zero[3] = {0.0, 0.0, 0.0}; @@ -1916,25 +1917,25 @@ void BSplSLib::Unperiodize // stored in homogeneous form //======================================================================= -void BSplSLib::BuildCache -(const Standard_Real U, - const Standard_Real V, - const Standard_Real USpanDomain, - const Standard_Real VSpanDomain, - const Standard_Boolean UPeriodic, - const Standard_Boolean VPeriodic, - const Standard_Integer UDegree, - const Standard_Integer VDegree, - const Standard_Integer UIndex, - const Standard_Integer VIndex, - const TColStd_Array1OfReal& UFlatKnots, - const TColStd_Array1OfReal& VFlatKnots, - const TColgp_Array2OfPnt& Poles, - const TColStd_Array2OfReal& Weights, - TColgp_Array2OfPnt& CachePoles, - TColStd_Array2OfReal& CacheWeights) -{ - Standard_Boolean rational,rational_u,rational_v,flag_u_or_v; +void BSplSLib::BuildCache( + const Standard_Real U, + const Standard_Real V, + const Standard_Real USpanDomain, + const Standard_Real VSpanDomain, + const Standard_Boolean UPeriodic, + const Standard_Boolean VPeriodic, + const Standard_Integer UDegree, + const Standard_Integer VDegree, + const Standard_Integer UIndex, + const Standard_Integer VIndex, + const TColStd_Array1OfReal& UFlatKnots, + const TColStd_Array1OfReal& VFlatKnots, + const TColgp_Array2OfPnt& Poles, + const TColStd_Array2OfReal& Weights, + Standard_Boolean& CacheRational, + TColStd_Array2OfReal& CachePolesWeights) +{ + Standard_Boolean rational,rational_u,rational_v,flag_u_or_v; Standard_Integer kk,d1,d1p1,d2,d2p1,ii,jj,iii,jjj,Index; Standard_Real u1,min_degree_domain,max_degree_domain,f,factor[2],u2; if (&Weights != NULL) @@ -1943,741 +1944,488 @@ void BSplSLib::BuildCache rational_u = rational_v = Standard_False; BSplSLib_DataContainer dc (UDegree, VDegree); flag_u_or_v = - PrepareEval (U, - V, - UIndex, - VIndex, - UDegree, - VDegree, - rational_u, - rational_v, - UPeriodic, - VPeriodic, - Poles, - Weights, - UFlatKnots, - VFlatKnots, - (BSplCLib::NoMults()), - (BSplCLib::NoMults()), - u1, - u2, - d1, - d2, - rational, - dc); + PrepareEval(U, + V, + UIndex, + VIndex, + UDegree, + VDegree, + rational_u, + rational_v, + UPeriodic, + VPeriodic, + Poles, + Weights, + UFlatKnots, + VFlatKnots, + (BSplCLib::NoMults()), + (BSplCLib::NoMults()), + u1, + u2, + d1, + d2, + rational, + dc); d1p1 = d1 + 1; d2p1 = d2 + 1; - if (rational) { - BSplCLib::Bohm(u1,d1,d1,*dc.knots1,4 * d2p1,*dc.poles); - - for (kk = 0; kk <= d1 ; kk++) - BSplCLib::Bohm(u2,d2,d2,*dc.knots2,4,*(dc.poles + kk * 4 * d2p1)); - if (flag_u_or_v) { - min_degree_domain = USpanDomain ; - max_degree_domain = VSpanDomain ; + CacheRational = rational; + + if (rational) + { + BSplCLib::Bohm(u1, d1, d1, *dc.knots1, 4 * d2p1, *dc.poles); + + for (kk = 0; kk <= d1; kk++) + BSplCLib::Bohm(u2, d2, d2, *dc.knots2, 4, *(dc.poles + kk * 4 * d2p1)); + + if (flag_u_or_v) + { + min_degree_domain = USpanDomain; + max_degree_domain = VSpanDomain; } - else { - min_degree_domain = VSpanDomain ; - max_degree_domain = USpanDomain ; + else + { + min_degree_domain = VSpanDomain; + max_degree_domain = USpanDomain; } + factor[0] = 1.0e0 ; - - for (ii = 0 ; ii <= d2 ; ii++) { + for (ii = 0; ii <= d2; ii++) + { iii = ii + 1; - factor[1] = 1.0e0 ; - - for (jj = 0 ; jj <= d1 ; jj++) { - jjj = jj + 1; - Index = jj * d2p1 + ii ; - Index = Index << 2; - gp_Pnt& P = CachePoles(iii,jjj); - f = factor[0] * factor[1]; - P.SetX( f * dc.poles[Index]); Index++; - P.SetY( f * dc.poles[Index]); Index++; - P.SetZ( f * dc.poles[Index]); Index++; - CacheWeights(iii ,jjj) = f * dc.poles[Index] ; - factor[1] *= min_degree_domain / (Standard_Real) (jjj) ; + factor[1] = 1.0e0; + + for (jj = 0; jj <= d1; jj++) + { + jjj = jj * 4; + Index = jj * d2p1 + ii; + Index = Index << 2; + f = factor[0] * factor[1]; + for (kk = 1; kk <= 4; kk++) + { + CachePolesWeights(iii, jjj + kk) = f * dc.poles[Index]; + Index++; + } + factor[1] *= min_degree_domain / (Standard_Real)(jj + 1); } - factor[0] *= max_degree_domain / (Standard_Real) (iii) ; + factor[0] *= max_degree_domain / (Standard_Real)(iii); } } - else { - BSplCLib::Bohm(u1,d1,d1,*dc.knots1,3 * d2p1,*dc.poles); - - for (kk = 0; kk <= d1 ; kk++) - BSplCLib::Bohm(u2,d2,d2,*dc.knots2,3,*(dc.poles + kk * 3 * d2p1)); - if (flag_u_or_v) { - min_degree_domain = USpanDomain ; - max_degree_domain = VSpanDomain ; + else + { + BSplCLib::Bohm(u1, d1, d1, *dc.knots1, 3 * d2p1, *dc.poles); + + for (kk = 0; kk <= d1; kk++) + BSplCLib::Bohm(u2, d2, d2, *dc.knots2, 3, *(dc.poles + kk * 3 * d2p1)); + + if (flag_u_or_v) + { + min_degree_domain = USpanDomain; + max_degree_domain = VSpanDomain; } - else { - min_degree_domain = VSpanDomain ; - max_degree_domain = USpanDomain ; + else + { + min_degree_domain = VSpanDomain; + max_degree_domain = USpanDomain; } - factor[0] = 1.0e0 ; - - for (ii = 0 ; ii <= d2 ; ii++) { + + factor[0] = 1.0e0; + for (ii = 0; ii <= d2; ii++) + { iii = ii + 1; - factor[1] = 1.0e0 ; - - for (jj = 0 ; jj <= d1 ; jj++) { - jjj = jj + 1; - Index = jj * d2p1 + ii ; - Index = (Index << 1) + Index; - gp_Pnt& P = CachePoles(iii,jjj); - f = factor[0] * factor[1]; - P.SetX( f * dc.poles[Index]); Index++; - P.SetY( f * dc.poles[Index]); Index++; - P.SetZ( f * dc.poles[Index]); - factor[1] *= min_degree_domain / (Standard_Real) (jjj) ; - } - factor[0] *= max_degree_domain / (Standard_Real) (iii) ; - } - if (&Weights != NULL) { - // - // means that PrepareEval did found out that the surface was - // locally polynomial but since the surface is constructed - // with some weights we need to set the weight polynomial to constant - // - - for (ii = 1 ; ii <= d2p1 ; ii++) { - - for (jj = 1 ; jj <= d1p1 ; jj++) { - CacheWeights(ii,jj) = 0.0e0 ; - } + factor[1] = 1.0e0; + + for (jj = 0; jj <= d1; jj++) + { + jjj = jj * 4; + Index = jj * d2p1 + ii; + Index = (Index << 1) + Index; + f = factor[0] * factor[1]; + for (kk = 1; kk <= 3; kk++) + { + CachePolesWeights(iii, jjj + kk) = f * dc.poles[Index]; + Index++; + } + CachePolesWeights(iii, jjj + 4) = 0.0e0; + factor[1] *= min_degree_domain / (Standard_Real)(jj + 1); } - CacheWeights(1,1) = 1.0e0 ; + factor[0] *= max_degree_domain / (Standard_Real)(iii); } + CachePolesWeights(1, 4) = 1.0e0; } } + //======================================================================= //function : CacheD0 //purpose : Evaluates the polynomial cache of the Bspline Curve // //======================================================================= -void BSplSLib::CacheD0(const Standard_Real UParameter, - const Standard_Real VParameter, - const Standard_Integer UDegree, - const Standard_Integer VDegree, - const Standard_Real UCacheParameter, - const Standard_Real VCacheParameter, - const Standard_Real USpanLenght, - const Standard_Real VSpanLenght, - const TColgp_Array2OfPnt& PolesArray, - const TColStd_Array2OfReal& WeightsArray, - gp_Pnt& aPoint) +void BSplSLib::CacheD0( + const Standard_Real UParameter, + const Standard_Real VParameter, + const Standard_Integer UDegree, + const Standard_Integer VDegree, + const Standard_Real UCacheParameter, + const Standard_Real VCacheParameter, + const Standard_Real USpanLenght, + const Standard_Real VSpanLenght, + const TColStd_Array2OfReal& PolesWeightsArray, + const Standard_Boolean CacheRational, + gp_Pnt& aPoint) { // // the CacheParameter is where the cache polynomial was evaluated in homogeneous // form // the SpanLenght is the normalizing factor so that the polynomial is between // 0 and 1 - Standard_Integer -// ii, - dimension, - min_degree, - max_degree ; - - Standard_Real - new_parameter[2] , - inverse ; - - Standard_Real * - PArray = (Standard_Real *) - &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ; - Standard_Real * - myPoint = (Standard_Real *) &aPoint ; - if (UDegree <= VDegree) { - min_degree = UDegree ; - max_degree = VDegree ; - new_parameter[1] = (UParameter - UCacheParameter) / USpanLenght ; - new_parameter[0] = (VParameter - VCacheParameter) / VSpanLenght ; - dimension = 3 * (UDegree + 1) ; + Standard_Integer dimension, min_degree, max_degree; + Standard_Real new_parameter[2], inverse; + MEMALIGN(Standard_Real *PArray) = (Standard_Real *) &(PolesWeightsArray(PolesWeightsArray.LowerCol(), PolesWeightsArray.LowerRow())); + Standard_Real *myPoint = (Standard_Real *) &aPoint; + MEMALIGN(Standard_Real local_pole_and_weight[4]); + + if (UDegree <= VDegree) + { + min_degree = UDegree; + max_degree = VDegree; + new_parameter[1] = (UParameter - UCacheParameter) / USpanLenght; + new_parameter[0] = (VParameter - VCacheParameter) / VSpanLenght; + dimension = 4 * (UDegree + 1); } - else { - min_degree = VDegree ; - max_degree = UDegree ; - new_parameter[0] = (UParameter - UCacheParameter) / USpanLenght ; - new_parameter[1] = (VParameter - VCacheParameter) / VSpanLenght ; - dimension = 3 * (VDegree + 1) ; + else + { + min_degree = VDegree; + max_degree = UDegree; + new_parameter[0] = (UParameter - UCacheParameter) / USpanLenght; + new_parameter[1] = (VParameter - VCacheParameter) / VSpanLenght; + dimension = 4 * (VDegree + 1); } - NCollection_LocalArray locpoles(dimension); - + + Standard_Real* locpoles = (Standard_Real*) NCollection_AlignAllocator::Allocate(dimension*sizeof(Standard_Real), DATA_ALIGNMENT); + PLib::NoDerivativeEvalPolynomial(new_parameter[0], - max_degree, - dimension, - max_degree*dimension, - PArray[0], - locpoles[0]) ; - + max_degree, + dimension, + max_degree*dimension, + PArray[0], + locpoles[0]); + PLib::NoDerivativeEvalPolynomial(new_parameter[1], - min_degree, - 3, - (min_degree << 1) + min_degree, - locpoles[0], - myPoint[0]) ; - if (&WeightsArray != NULL) { - dimension = min_degree + 1 ; - Standard_Real * - WArray = (Standard_Real *) - &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ; - PLib::NoDerivativeEvalPolynomial(new_parameter[0], - max_degree, - dimension, - max_degree*dimension, - WArray[0], - locpoles[0]) ; - - PLib::NoDerivativeEvalPolynomial(new_parameter[1], - min_degree, - 1, - min_degree, - locpoles[0], - inverse) ; - inverse = 1.0e0/ inverse ; - - myPoint[0] *= inverse ; - myPoint[1] *= inverse ; - myPoint[2] *= inverse ; + min_degree, + 4, + min_degree * 4, + locpoles[0], + local_pole_and_weight[0]); + + if (CacheRational) + { + inverse = 1.0e0 / local_pole_and_weight[3]; + + myPoint[0] = local_pole_and_weight[0] * inverse; + myPoint[1] = local_pole_and_weight[1] * inverse; + myPoint[2] = local_pole_and_weight[2] * inverse; + } + else + { + myPoint[0] = local_pole_and_weight[0]; + myPoint[1] = local_pole_and_weight[1]; + myPoint[2] = local_pole_and_weight[2]; } + + NCollection_AlignAllocator::Free(locpoles); } + //======================================================================= //function : CacheD1 //purpose : Evaluates the polynomial cache of the Bspline Curve // //======================================================================= -void BSplSLib::CacheD1(const Standard_Real UParameter, - const Standard_Real VParameter, - const Standard_Integer UDegree, - const Standard_Integer VDegree, - const Standard_Real UCacheParameter, - const Standard_Real VCacheParameter, - const Standard_Real USpanLenght, - const Standard_Real VSpanLenght, - const TColgp_Array2OfPnt& PolesArray, - const TColStd_Array2OfReal& WeightsArray, - gp_Pnt& aPoint, - gp_Vec& aVecU, - gp_Vec& aVecV) +void BSplSLib::CacheD1( + const Standard_Real UParameter, + const Standard_Real VParameter, + const Standard_Integer UDegree, + const Standard_Integer VDegree, + const Standard_Real UCacheParameter, + const Standard_Real VCacheParameter, + const Standard_Real USpanLenght, + const Standard_Real VSpanLenght, + const TColStd_Array2OfReal& PolesWeightsArray, + const Standard_Boolean CacheRational, + gp_Pnt& aPoint, + gp_Vec& aVecU, + gp_Vec& aVecV) { - // // the CacheParameter is where the cache polynomial was evaluated in homogeneous // form // the SpanLenght is the normalizing factor so that the polynomial is between // 0 and 1 - Standard_Integer -// ii, -// jj, -// kk, - dimension, - min_degree, - max_degree ; - - Standard_Real - inverse_min, - inverse_max, - new_parameter[2] ; - - Standard_Real * - PArray = (Standard_Real *) - &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ; - Standard_Real local_poles_array[2][2][3], - local_poles_and_weights_array[2][2][4], - local_weights_array[2][2] ; - Standard_Real * my_vec_min, - * my_vec_max, - * my_point ; - my_point = (Standard_Real *) &aPoint ; - // - // initialize in case of rational evaluation - // because RationalDerivative will use all - // the coefficients - // - // - if (&WeightsArray != NULL) { - - local_poles_array [0][0][0] = 0.0e0 ; - local_poles_array [0][0][1] = 0.0e0 ; - local_poles_array [0][0][2] = 0.0e0 ; - local_weights_array [0][0] = 0.0e0 ; - local_poles_and_weights_array[0][0][0] = 0.0e0 ; - local_poles_and_weights_array[0][0][1] = 0.0e0 ; - local_poles_and_weights_array[0][0][2] = 0.0e0 ; - local_poles_and_weights_array[0][0][3] = 0.0e0 ; - - local_poles_array [0][1][0] = 0.0e0 ; - local_poles_array [0][1][1] = 0.0e0 ; - local_poles_array [0][1][2] = 0.0e0 ; - local_weights_array [0][1] = 0.0e0 ; - local_poles_and_weights_array[0][1][0] = 0.0e0 ; - local_poles_and_weights_array[0][1][1] = 0.0e0 ; - local_poles_and_weights_array[0][1][2] = 0.0e0 ; - local_poles_and_weights_array[0][1][3] = 0.0e0 ; - - local_poles_array [1][0][0] = 0.0e0 ; - local_poles_array [1][0][1] = 0.0e0 ; - local_poles_array [1][0][2] = 0.0e0 ; - local_weights_array [1][0] = 0.0e0 ; - local_poles_and_weights_array[1][0][0] = 0.0e0 ; - local_poles_and_weights_array[1][0][1] = 0.0e0 ; - local_poles_and_weights_array[1][0][2] = 0.0e0 ; - local_poles_and_weights_array[1][0][3] = 0.0e0 ; - - local_poles_array [1][1][0] = 0.0e0 ; - local_poles_array [1][1][1] = 0.0e0 ; - local_poles_array [1][1][2] = 0.0e0 ; - local_weights_array [1][1] = 0.0e0 ; - local_poles_and_weights_array[1][1][0] = 0.0e0 ; - local_poles_and_weights_array[1][1][1] = 0.0e0 ; - local_poles_and_weights_array[1][1][2] = 0.0e0 ; - local_poles_and_weights_array[1][1][3] = 0.0e0 ; - } - - if (UDegree <= VDegree) { - min_degree = UDegree ; - max_degree = VDegree ; - inverse_min = 1.0e0/USpanLenght ; - inverse_max = 1.0e0/VSpanLenght ; - new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ; - new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ; - - dimension = 3 * (UDegree + 1) ; - my_vec_min = (Standard_Real *) &aVecU ; - my_vec_max = (Standard_Real *) &aVecV ; + Standard_Integer dimension, min_degree, max_degree, ii; + Standard_Real inverse_min, inverse_max, new_parameter[2]; + MEMALIGN(Standard_Real *PArray) = (Standard_Real *) &(PolesWeightsArray(PolesWeightsArray.LowerCol(), PolesWeightsArray.LowerRow())); + MEMALIGN(Standard_Real local_poles_array[4][3]); + MEMALIGN(Standard_Real local_poles_and_weights_array[4][4]); + Standard_Real *my_vec_min, *my_vec_max, *my_point; + + my_point = (Standard_Real *) &aPoint; + + for (ii = 0; ii < 3; ii++) + { + local_poles_array [3][ii] = 0.0e0; + local_poles_and_weights_array[3][ii] = 0.0e0; } - else { - min_degree = VDegree ; - max_degree = UDegree ; - inverse_min = 1.0e0/VSpanLenght ; - inverse_max = 1.0e0/USpanLenght ; - new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ; - new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ; - dimension = 3 * (VDegree + 1) ; - my_vec_min = (Standard_Real *) &aVecV ; - my_vec_max = (Standard_Real *) &aVecU ; + local_poles_and_weights_array[3][3] = 0.0e0; + + if (UDegree <= VDegree) + { + min_degree = UDegree; + max_degree = VDegree; + inverse_min = 1.0e0 / USpanLenght; + inverse_max = 1.0e0 / VSpanLenght; + new_parameter[0] = (VParameter - VCacheParameter) * inverse_max; + new_parameter[1] = (UParameter - UCacheParameter) * inverse_min; + dimension = 4 * (UDegree + 1); + my_vec_min = (Standard_Real *) &aVecU; + my_vec_max = (Standard_Real *) &aVecV; + } + else + { + min_degree = VDegree; + max_degree = UDegree; + inverse_min = 1.0e0 / VSpanLenght; + inverse_max = 1.0e0 / USpanLenght; + new_parameter[0] = (UParameter - UCacheParameter) * inverse_max; + new_parameter[1] = (VParameter - VCacheParameter) * inverse_min; + dimension = 4 * (VDegree + 1); + my_vec_min = (Standard_Real *) &aVecV; + my_vec_max = (Standard_Real *) &aVecU; } - NCollection_LocalArray locpoles (2 * dimension); - + Standard_Real* locpoles = (Standard_Real*) NCollection_AlignAllocator::Allocate(2*dimension*sizeof(Standard_Real), DATA_ALIGNMENT); + PLib::EvalPolynomial(new_parameter[0], - 1, - max_degree, - dimension, - PArray[0], - locpoles[0]) ; - + 1, + max_degree, + dimension, + PArray[0], + locpoles[0]); + PLib::EvalPolynomial(new_parameter[1], - 1, - min_degree, - 3, - locpoles[0], - local_poles_array[0][0][0]) ; + 1, + min_degree, + 4, + locpoles[0], + local_poles_and_weights_array[0][0]); + PLib::NoDerivativeEvalPolynomial(new_parameter[1], - min_degree, - 3, - (min_degree << 1) + min_degree, - locpoles[dimension], - local_poles_array[1][0][0]) ; - - if (&WeightsArray != NULL) { - dimension = min_degree + 1 ; - Standard_Real * - WArray = (Standard_Real *) - &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ; - PLib::EvalPolynomial(new_parameter[0], - 1, - max_degree, - dimension, - WArray[0], - locpoles[0]) ; - - PLib::EvalPolynomial(new_parameter[1], - 1, - min_degree, - 1, - locpoles[0], - local_weights_array[0][0]) ; - PLib::NoDerivativeEvalPolynomial(new_parameter[1], - min_degree, - 1, - min_degree, - locpoles[dimension], - local_weights_array[1][0]) ; - - local_poles_and_weights_array[0][0][0] = local_poles_array [0][0][0] ; - local_poles_and_weights_array[0][0][1] = local_poles_array [0][0][1] ; - local_poles_and_weights_array[0][0][2] = local_poles_array [0][0][2] ; - local_poles_and_weights_array[0][0][3] = local_weights_array[0][0] ; - - local_poles_and_weights_array[0][1][0] = local_poles_array [0][1][0] ; - local_poles_and_weights_array[0][1][1] = local_poles_array [0][1][1] ; - local_poles_and_weights_array[0][1][2] = local_poles_array [0][1][2] ; - local_poles_and_weights_array[0][1][3] = local_weights_array[0][1] ; - - local_poles_and_weights_array[1][0][0] = local_poles_array [1][0][0] ; - local_poles_and_weights_array[1][0][1] = local_poles_array [1][0][1] ; - local_poles_and_weights_array[1][0][2] = local_poles_array [1][0][2] ; - local_poles_and_weights_array[1][0][3] = local_weights_array[1][0] ; - - local_poles_and_weights_array[1][1][0] = local_poles_array [1][1][0] ; - local_poles_and_weights_array[1][1][1] = local_poles_array [1][1][1] ; - local_poles_and_weights_array[1][1][2] = local_poles_array [1][1][2] ; - local_poles_and_weights_array[1][1][3] = local_weights_array[1][1] ; + min_degree, + 4, + min_degree * 4, + locpoles[dimension], + local_poles_and_weights_array[2][0]); + if (CacheRational) + { BSplSLib::RationalDerivative(1, - 1, - 1, - 1, - local_poles_and_weights_array[0][0][0], - local_poles_array[0][0][0]) ; + 1, + 1, + 1, + local_poles_and_weights_array[0][0], + local_poles_array[0][0]); + + for (ii = 0; ii < 3; ii++) + { + my_point [ii] = local_poles_array [0][ii]; + my_vec_min[ii] = inverse_min * local_poles_array[1][ii]; + my_vec_max[ii] = inverse_max * local_poles_array[2][ii]; + } } - - my_point [0] = local_poles_array [0][0][0] ; - my_vec_min[0] = inverse_min * local_poles_array[0][1][0] ; - my_vec_max[0] = inverse_max * local_poles_array[1][0][0] ; - - my_point [1] = local_poles_array [0][0][1] ; - my_vec_min[1] = inverse_min * local_poles_array[0][1][1] ; - my_vec_max[1] = inverse_max * local_poles_array[1][0][1] ; - - my_point [2] = local_poles_array [0][0][2] ; - my_vec_min[2] = inverse_min * local_poles_array[0][1][2] ; - my_vec_max[2] = inverse_max * local_poles_array[1][0][2] ; + else + { + for (ii = 0; ii < 3; ii++) + { + my_point [ii] = local_poles_and_weights_array [0][ii]; + my_vec_min[ii] = inverse_min * local_poles_and_weights_array[1][ii]; + my_vec_max[ii] = inverse_max * local_poles_and_weights_array[2][ii]; + } + } + + NCollection_AlignAllocator::Free(locpoles); } + //======================================================================= //function : CacheD2 //purpose : Evaluates the polynomial cache of the Bspline Curve // //======================================================================= -void BSplSLib::CacheD2(const Standard_Real UParameter, - const Standard_Real VParameter, - const Standard_Integer UDegree, - const Standard_Integer VDegree, - const Standard_Real UCacheParameter, - const Standard_Real VCacheParameter, - const Standard_Real USpanLenght, - const Standard_Real VSpanLenght, - const TColgp_Array2OfPnt& PolesArray, - const TColStd_Array2OfReal& WeightsArray, - gp_Pnt& aPoint, - gp_Vec& aVecU, - gp_Vec& aVecV, - gp_Vec& aVecUU, - gp_Vec& aVecUV, - gp_Vec& aVecVV) +void BSplSLib::CacheD2( + const Standard_Real UParameter, + const Standard_Real VParameter, + const Standard_Integer UDegree, + const Standard_Integer VDegree, + const Standard_Real UCacheParameter, + const Standard_Real VCacheParameter, + const Standard_Real USpanLength, + const Standard_Real VSpanLength, + const TColStd_Array2OfReal& PolesWeightsArray, + const Standard_Boolean CacheRational, + gp_Pnt& aPoint, + gp_Vec& aVecU, + gp_Vec& aVecV, + gp_Vec& aVecUU, + gp_Vec& aVecUV, + gp_Vec& aVecVV) { // // the CacheParameter is where the cache polynomial was evaluated in homogeneous // form // the SpanLenght is the normalizing factor so that the polynomial is between // 0 and 1 - Standard_Integer - ii, -// jj, - kk, - index, - dimension, - min_degree, - max_degree ; - - Standard_Real - inverse_min, - inverse_max, - new_parameter[2] ; - - Standard_Real * - PArray = (Standard_Real *) - &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ; - Standard_Real local_poles_array[3][3][3], - local_poles_and_weights_array[3][3][4], - local_weights_array[3][3] ; - Standard_Real * my_vec_min, - * my_vec_max, - * my_vec_min_min, - * my_vec_max_max, - * my_vec_min_max, - * my_point ; + Standard_Integer ii, kk, index, dimension, min_degree, max_degree; + Standard_Real inverse_min, inverse_max, new_parameter[2]; + MEMALIGN(Standard_Real *PArray) = (Standard_Real *) &(PolesWeightsArray(PolesWeightsArray.LowerCol(), PolesWeightsArray.LowerRow())); + MEMALIGN(Standard_Real local_poles_array[9][3]); + MEMALIGN(Standard_Real local_poles_and_weights_array[9][4]); + Standard_Real *my_vec_min, *my_vec_max, *my_vec_min_min, *my_vec_max_max, *my_vec_min_max, *my_point; my_point = (Standard_Real *) &aPoint ; - + // // initialize in case the min and max degree are less than 2 // - local_poles_array[0][0][0] = 0.0e0 ; - local_poles_array[0][0][1] = 0.0e0 ; - local_poles_array[0][0][2] = 0.0e0 ; - local_poles_array[0][1][0] = 0.0e0 ; - local_poles_array[0][1][1] = 0.0e0 ; - local_poles_array[0][1][2] = 0.0e0 ; - local_poles_array[0][2][0] = 0.0e0 ; - local_poles_array[0][2][1] = 0.0e0 ; - local_poles_array[0][2][2] = 0.0e0 ; - - local_poles_array[1][0][0] = 0.0e0 ; - local_poles_array[1][0][1] = 0.0e0 ; - local_poles_array[1][0][2] = 0.0e0 ; - local_poles_array[1][1][0] = 0.0e0 ; - local_poles_array[1][1][1] = 0.0e0 ; - local_poles_array[1][1][2] = 0.0e0 ; - local_poles_array[1][2][0] = 0.0e0 ; - local_poles_array[1][2][1] = 0.0e0 ; - local_poles_array[1][2][2] = 0.0e0 ; - - local_poles_array[2][0][0] = 0.0e0 ; - local_poles_array[2][0][1] = 0.0e0 ; - local_poles_array[2][0][2] = 0.0e0 ; - local_poles_array[2][1][0] = 0.0e0 ; - local_poles_array[2][1][1] = 0.0e0 ; - local_poles_array[2][1][2] = 0.0e0 ; - local_poles_array[2][2][0] = 0.0e0 ; - local_poles_array[2][2][1] = 0.0e0 ; - local_poles_array[2][2][2] = 0.0e0 ; - // - // initialize in case of rational evaluation + // and in case of rational evaluation // because RationalDerivative will use all // the coefficients // - // - if (&WeightsArray != NULL) { - - local_poles_and_weights_array[0][0][0] = 0.0e0 ; - local_poles_and_weights_array[0][0][1] = 0.0e0 ; - local_poles_and_weights_array[0][0][2] = 0.0e0 ; - local_poles_and_weights_array[0][1][0] = 0.0e0 ; - local_poles_and_weights_array[0][1][1] = 0.0e0 ; - local_poles_and_weights_array[0][1][2] = 0.0e0 ; - local_poles_and_weights_array[0][2][0] = 0.0e0 ; - local_poles_and_weights_array[0][2][1] = 0.0e0 ; - local_poles_and_weights_array[0][2][2] = 0.0e0 ; - - local_poles_and_weights_array[1][0][0] = 0.0e0 ; - local_poles_and_weights_array[1][0][1] = 0.0e0 ; - local_poles_and_weights_array[1][0][2] = 0.0e0 ; - local_poles_and_weights_array[1][1][0] = 0.0e0 ; - local_poles_and_weights_array[1][1][1] = 0.0e0 ; - local_poles_and_weights_array[1][1][2] = 0.0e0 ; - local_poles_and_weights_array[1][2][0] = 0.0e0 ; - local_poles_and_weights_array[1][2][1] = 0.0e0 ; - local_poles_and_weights_array[1][2][2] = 0.0e0 ; - - local_poles_and_weights_array[2][0][0] = 0.0e0 ; - local_poles_and_weights_array[2][0][1] = 0.0e0 ; - local_poles_and_weights_array[2][0][2] = 0.0e0 ; - local_poles_and_weights_array[2][1][0] = 0.0e0 ; - local_poles_and_weights_array[2][1][1] = 0.0e0 ; - local_poles_and_weights_array[2][1][2] = 0.0e0 ; - local_poles_and_weights_array[2][2][0] = 0.0e0 ; - local_poles_and_weights_array[2][2][1] = 0.0e0 ; - local_poles_and_weights_array[2][2][2] = 0.0e0 ; - - local_poles_and_weights_array[0][0][3] = - local_weights_array[0][0] = 0.0e0 ; - local_poles_and_weights_array[0][1][3] = - local_weights_array[0][1] = 0.0e0 ; - local_poles_and_weights_array[0][2][3] = - local_weights_array[0][2] = 0.0e0 ; - local_poles_and_weights_array[1][0][3] = - local_weights_array[1][0] = 0.0e0 ; - local_poles_and_weights_array[1][1][3] = - local_weights_array[1][1] = 0.0e0 ; - local_poles_and_weights_array[1][2][3] = - local_weights_array[1][2] = 0.0e0 ; - local_poles_and_weights_array[2][0][3] = - local_weights_array[2][0] = 0.0e0 ; - local_poles_and_weights_array[2][1][3] = - local_weights_array[2][1] = 0.0e0 ; - local_poles_and_weights_array[2][2][3] = - local_weights_array[2][2] = 0.0e0 ; - } - - if (UDegree <= VDegree) { - min_degree = UDegree ; - max_degree = VDegree ; - inverse_min = 1.0e0/USpanLenght ; - inverse_max = 1.0e0/VSpanLenght ; - new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ; - new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ; - - dimension = 3 * (UDegree + 1) ; - my_vec_min = (Standard_Real *) &aVecU ; - my_vec_max = (Standard_Real *) &aVecV ; - my_vec_min_min = (Standard_Real *) &aVecUU ; - my_vec_min_max = (Standard_Real *) &aVecUV ; - my_vec_max_max = (Standard_Real *) &aVecVV ; + for (ii = 0; ii < 9; ii++) + { + for (kk = 0; kk < 3; kk++) + local_poles_array[ii][kk] = 0.0e0; + for (kk = 0; kk < 4; kk++) + local_poles_and_weights_array[ii][kk] = 0.0e0; } - else { - min_degree = VDegree ; - max_degree = UDegree ; - inverse_min = 1.0e0/VSpanLenght ; - inverse_max = 1.0e0/USpanLenght ; - new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ; - new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ; - dimension = 3 * (VDegree + 1) ; - my_vec_min = (Standard_Real *) &aVecV ; - my_vec_max = (Standard_Real *) &aVecU ; - my_vec_min_min = (Standard_Real *) &aVecVV ; - my_vec_min_max = (Standard_Real *) &aVecUV ; - my_vec_max_max = (Standard_Real *) &aVecUU ; - } - - NCollection_LocalArray locpoles (3 * dimension); - + + if (UDegree <= VDegree) + { + min_degree = UDegree; + max_degree = VDegree; + inverse_min = 1.0e0 / USpanLength; + inverse_max = 1.0e0 / VSpanLength; + new_parameter[0] = (VParameter - VCacheParameter) * inverse_max; + new_parameter[1] = (UParameter - UCacheParameter) * inverse_min; + dimension = 4 * (UDegree + 1); + my_vec_min = (Standard_Real *) &aVecU; + my_vec_max = (Standard_Real *) &aVecV; + my_vec_min_min = (Standard_Real *) &aVecUU; + my_vec_min_max = (Standard_Real *) &aVecUV; + my_vec_max_max = (Standard_Real *) &aVecVV; + } + else + { + min_degree = VDegree; + max_degree = UDegree; + inverse_min = 1.0e0 / VSpanLength; + inverse_max = 1.0e0 / USpanLength; + new_parameter[0] = (UParameter - UCacheParameter) * inverse_max; + new_parameter[1] = (VParameter - VCacheParameter) * inverse_min; + dimension = 4 * (VDegree + 1); + my_vec_min = (Standard_Real *) &aVecV; + my_vec_max = (Standard_Real *) &aVecU; + my_vec_min_min = (Standard_Real *) &aVecVV; + my_vec_min_max = (Standard_Real *) &aVecUV; + my_vec_max_max = (Standard_Real *) &aVecUU; + } + + Standard_Real* locpoles = (Standard_Real*) NCollection_AlignAllocator::Allocate(3*dimension*sizeof(Standard_Real), DATA_ALIGNMENT); + // // initialize in case min or max degree are less than 2 // Standard_Integer MinIndMax = 2; - if ( max_degree < 2) MinIndMax = max_degree; + if (max_degree < 2) MinIndMax = max_degree; Standard_Integer MinIndMin = 2; - if ( min_degree < 2) MinIndMin = min_degree; - + if (min_degree < 2) MinIndMin = min_degree; + index = MinIndMax * dimension ; - for (ii = MinIndMax ; ii < 3 ; ii++) { - - for (kk = 0 ; kk < dimension ; kk++) { + for (ii = MinIndMax; ii < 3 ; ii++) + { + for (kk = 0 ; kk < dimension ; kk++) + { locpoles[index] = 0.0e0 ; index += 1 ; } } - + PLib::EvalPolynomial(new_parameter[0], - MinIndMax, - max_degree, - dimension, - PArray[0], - locpoles[0]) ; - + MinIndMax, + max_degree, + dimension, + PArray[0], + locpoles[0]); + PLib::EvalPolynomial(new_parameter[1], - MinIndMin, - min_degree, - 3, - locpoles[0], - local_poles_array[0][0][0]) ; + MinIndMin, + min_degree, + 4, + locpoles[0], + local_poles_and_weights_array[0][0]); + PLib::EvalPolynomial(new_parameter[1], - 1, - min_degree, - 3, - locpoles[dimension], - local_poles_array[1][0][0]) ; - + 1, + min_degree, + 4, + locpoles[dimension], + local_poles_and_weights_array[3][0]); + PLib::NoDerivativeEvalPolynomial(new_parameter[1], - min_degree, - 3, - (min_degree << 1) + min_degree, - locpoles[dimension + dimension], - local_poles_array[2][0][0]) ; - - if (&WeightsArray != NULL) { - dimension = min_degree + 1 ; - Standard_Real * - WArray = (Standard_Real *) - &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ; - PLib::EvalPolynomial(new_parameter[0], - MinIndMax, - max_degree, - dimension, - WArray[0], - locpoles[0]) ; - - PLib::EvalPolynomial(new_parameter[1], - MinIndMin, - min_degree, - 1, - locpoles[0], - local_weights_array[0][0]) ; - PLib::EvalPolynomial(new_parameter[1], - 1, - min_degree, - 1, - locpoles[dimension], - local_weights_array[1][0]) ; - PLib::NoDerivativeEvalPolynomial(new_parameter[1], - min_degree, - 1, - min_degree, - locpoles[dimension + dimension], - local_weights_array[2][0]) ; - - - local_poles_and_weights_array[0][0][0] = local_poles_array[0][0][0]; - local_poles_and_weights_array[0][0][1] = local_poles_array[0][0][1]; - local_poles_and_weights_array[0][0][2] = local_poles_array[0][0][2]; - local_poles_and_weights_array[0][1][0] = local_poles_array[0][1][0]; - local_poles_and_weights_array[0][1][1] = local_poles_array[0][1][1]; - local_poles_and_weights_array[0][1][2] = local_poles_array[0][1][2]; - local_poles_and_weights_array[0][2][0] = local_poles_array[0][2][0]; - local_poles_and_weights_array[0][2][1] = local_poles_array[0][2][1]; - local_poles_and_weights_array[0][2][2] = local_poles_array[0][2][2]; - - local_poles_and_weights_array[1][0][0] = local_poles_array[1][0][0]; - local_poles_and_weights_array[1][0][1] = local_poles_array[1][0][1]; - local_poles_and_weights_array[1][0][2] = local_poles_array[1][0][2]; - local_poles_and_weights_array[1][1][0] = local_poles_array[1][1][0]; - local_poles_and_weights_array[1][1][1] = local_poles_array[1][1][1]; - local_poles_and_weights_array[1][1][2] = local_poles_array[1][1][2]; - local_poles_and_weights_array[1][2][0] = local_poles_array[1][2][0]; - local_poles_and_weights_array[1][2][1] = local_poles_array[1][2][1]; - local_poles_and_weights_array[1][2][2] = local_poles_array[1][2][2]; - - local_poles_and_weights_array[2][0][0] = local_poles_array[2][0][0]; - local_poles_and_weights_array[2][0][1] = local_poles_array[2][0][1]; - local_poles_and_weights_array[2][0][2] = local_poles_array[2][0][2]; - local_poles_and_weights_array[2][1][0] = local_poles_array[2][1][0]; - local_poles_and_weights_array[2][1][1] = local_poles_array[2][1][1]; - local_poles_and_weights_array[2][1][2] = local_poles_array[2][1][2]; - local_poles_and_weights_array[2][2][0] = local_poles_array[2][2][0]; - local_poles_and_weights_array[2][2][1] = local_poles_array[2][2][1]; - local_poles_and_weights_array[2][2][2] = local_poles_array[2][2][2]; - - - local_poles_and_weights_array[0][0][3] = local_weights_array[0][0]; - local_poles_and_weights_array[0][1][3] = local_weights_array[0][1]; - local_poles_and_weights_array[0][2][3] = local_weights_array[0][2]; - local_poles_and_weights_array[1][0][3] = local_weights_array[1][0]; - local_poles_and_weights_array[1][1][3] = local_weights_array[1][1]; - local_poles_and_weights_array[1][2][3] = local_weights_array[1][2]; - local_poles_and_weights_array[2][0][3] = local_weights_array[2][0]; - local_poles_and_weights_array[2][1][3] = local_weights_array[2][1]; - local_poles_and_weights_array[2][2][3] = local_weights_array[2][2]; - - BSplSLib::RationalDerivative(2, - 2, - 2, - 2, - local_poles_and_weights_array[0][0][0], - local_poles_array[0][0][0]) ; - } - + min_degree, + 4, + (min_degree << 2), + locpoles[ (dimension<<1) ], + local_poles_and_weights_array[6][0]); + Standard_Real minmin = inverse_min * inverse_min; Standard_Real minmax = inverse_min * inverse_max; Standard_Real maxmax = inverse_max * inverse_max; - - my_point [0] = local_poles_array [0][0][0] ; - my_vec_min [0] = inverse_min * local_poles_array[0][1][0] ; - my_vec_max [0] = inverse_max * local_poles_array[1][0][0] ; - my_vec_min_min[0] = minmin * local_poles_array [0][2][0] ; - my_vec_min_max[0] = minmax * local_poles_array [1][1][0] ; - my_vec_max_max[0] = maxmax * local_poles_array [2][0][0] ; - - my_point [1] = local_poles_array [0][0][1] ; - my_vec_min [1] = inverse_min * local_poles_array[0][1][1] ; - my_vec_max [1] = inverse_max * local_poles_array[1][0][1] ; - my_vec_min_min[1] = minmin * local_poles_array [0][2][1] ; - my_vec_min_max[1] = minmax * local_poles_array [1][1][1] ; - my_vec_max_max[1] = maxmax * local_poles_array [2][0][1] ; - - my_point [2] = local_poles_array [0][0][2] ; - my_vec_min [2] = inverse_min * local_poles_array[0][1][2] ; - my_vec_max [2] = inverse_max * local_poles_array[1][0][2] ; - my_vec_min_min[2] = minmin * local_poles_array [0][2][2] ; - my_vec_min_max[2] = minmax * local_poles_array [1][1][2] ; - my_vec_max_max[2] = maxmax * local_poles_array [2][0][2] ; + + if (CacheRational) + { + BSplSLib::RationalDerivative(2, + 2, + 2, + 2, + local_poles_and_weights_array[0][0], + local_poles_array[0][0]); + + for (ii = 0; ii < 3; ii++) + { + my_point [ii] = local_poles_array [0][ii] ; + my_vec_min [ii] = inverse_min * local_poles_array[1][ii] ; + my_vec_max [ii] = inverse_max * local_poles_array[3][ii] ; + my_vec_min_min[ii] = minmin * local_poles_array [2][ii] ; + my_vec_min_max[ii] = minmax * local_poles_array [4][ii] ; + my_vec_max_max[ii] = maxmax * local_poles_array [6][ii] ; + } + } + else + { + for (ii = 0; ii < 3; ii++) + { + my_point [ii] = local_poles_and_weights_array [0][ii] ; + my_vec_min [ii] = inverse_min * local_poles_and_weights_array[1][ii] ; + my_vec_max [ii] = inverse_max * local_poles_and_weights_array[3][ii] ; + my_vec_min_min[ii] = minmin * local_poles_and_weights_array [2][ii] ; + my_vec_min_max[ii] = minmax * local_poles_and_weights_array [4][ii] ; + my_vec_max_max[ii] = maxmax * local_poles_and_weights_array [6][ii] ; + } + } + + NCollection_AlignAllocator::Free(locpoles); } //======================================================================= diff --git a/src/BSplSLib/BSplSLib.lxx b/src/BSplSLib/BSplSLib.lxx index 5bcba9af8a..d80d9f7718 100755 --- a/src/BSplSLib/BSplSLib.lxx +++ b/src/BSplSLib/BSplSLib.lxx @@ -28,16 +28,17 @@ //purpose : //======================================================================= -inline void BSplSLib::CoefsD0(const Standard_Real U, - const Standard_Real V, - const TColgp_Array2OfPnt& Poles, - const TColStd_Array2OfReal& Weights, - gp_Pnt& Point) +inline void BSplSLib::CoefsD0( + const Standard_Real U, + const Standard_Real V, + const TColStd_Array2OfReal& PolesWeightsArray, + const Standard_Boolean CacheRational, + gp_Pnt& Point) { BSplSLib::CacheD0(U, V, - Poles.RowLength() - 1, Poles.ColLength() - 1, + PolesWeightsArray.RowLength() - 1, (PolesWeightsArray.ColLength()>>2) - 1, 0., 0., 1., 1., - Poles, Weights, Point); + PolesWeightsArray, CacheRational, Point); } //======================================================================= @@ -45,18 +46,19 @@ inline void BSplSLib::CoefsD0(const Standard_Real U, //purpose : //======================================================================= -inline void BSplSLib::CoefsD1(const Standard_Real U, - const Standard_Real V, - const TColgp_Array2OfPnt& Poles, - const TColStd_Array2OfReal& Weights, - gp_Pnt& Point, - gp_Vec& VecU, - gp_Vec& VecV) +inline void BSplSLib::CoefsD1( + const Standard_Real U, + const Standard_Real V, + const TColStd_Array2OfReal& PolesWeightsArray, + const Standard_Boolean CacheRational, + gp_Pnt& Point, + gp_Vec& VecU, + gp_Vec& VecV) { BSplSLib::CacheD1(U, V, - Poles.RowLength() - 1, Poles.ColLength() - 1, + PolesWeightsArray.RowLength() - 1, (PolesWeightsArray.ColLength()>>2) - 1, 0., 0., 1., 1., - Poles, Weights, Point, VecU, VecV); + PolesWeightsArray, CacheRational, Point, VecU, VecV); } //======================================================================= @@ -65,20 +67,20 @@ inline void BSplSLib::CoefsD1(const Standard_Real U, //======================================================================= inline void BSplSLib::CoefsD2(const Standard_Real U, - const Standard_Real V, - const TColgp_Array2OfPnt& Poles, - const TColStd_Array2OfReal& Weights, - gp_Pnt& Point, - gp_Vec& VecU, - gp_Vec& VecV, - gp_Vec& VecUU, - gp_Vec& VecUV, - gp_Vec& VecVV) + const Standard_Real V, + const TColStd_Array2OfReal& PolesWeightsArray, + const Standard_Boolean CacheRational, + gp_Pnt& Point, + gp_Vec& VecU, + gp_Vec& VecV, + gp_Vec& VecUU, + gp_Vec& VecUV, + gp_Vec& VecVV) { BSplSLib::CacheD2(U, V, - Poles.RowLength() - 1, Poles.ColLength() - 1, + PolesWeightsArray.RowLength() - 1, (PolesWeightsArray.ColLength()>>2) - 1, 0., 0., 1., 1., - Poles, Weights, Point, VecU, VecV, + PolesWeightsArray, CacheRational, Point, VecU, VecV, VecUU, VecUV, VecVV); } @@ -89,10 +91,11 @@ inline void BSplSLib::CoefsD2(const Standard_Real U, //======================================================================= inline void BSplSLib::PolesCoefficients(const TColgp_Array2OfPnt& Poles, - TColgp_Array2OfPnt& CachePoles) + TColStd_Array2OfReal& CachePolesWeights) { + Standard_Boolean aCacheRational; BSplSLib::PolesCoefficients(Poles, BSplSLib::NoWeights(), - CachePoles, BSplSLib::NoWeights()); + aCacheRational, CachePolesWeights); } //======================================================================= diff --git a/src/BSplSLib/BSplSLib_BzSyntaxes.cxx b/src/BSplSLib/BSplSLib_BzSyntaxes.cxx index a896b28c34..182062ff1c 100755 --- a/src/BSplSLib/BSplSLib_BzSyntaxes.cxx +++ b/src/BSplSLib/BSplSLib_BzSyntaxes.cxx @@ -34,10 +34,11 @@ //purpose : //======================================================================= -void BSplSLib::PolesCoefficients (const TColgp_Array2OfPnt& Poles, - const TColStd_Array2OfReal& Weights, - TColgp_Array2OfPnt& CachePoles, - TColStd_Array2OfReal& CacheWeights) +void BSplSLib::PolesCoefficients( + const TColgp_Array2OfPnt& Poles, + const TColStd_Array2OfReal& Weights, + Standard_Boolean& CacheRational, + TColStd_Array2OfReal& CachePolesWeights) { Standard_Integer i; Standard_Integer uclas = Poles.ColLength(); @@ -45,53 +46,56 @@ void BSplSLib::PolesCoefficients (const TColgp_Array2OfPnt& Poles, TColStd_Array1OfReal biduflatknots(1,uclas << 1); TColStd_Array1OfReal bidvflatknots(1,vclas << 1); - for(i = 1; i <= uclas; i++) { + for(i = 1; i <= uclas; i++) + { biduflatknots(i ) = 0.; biduflatknots(i + uclas) = 1.; } - for(i = 1; i <= vclas; i++) { + for(i = 1; i <= vclas; i++) + { bidvflatknots(i ) = 0.; bidvflatknots(i + vclas) = 1.; } - if ( uclas > vclas) { + if ( uclas > vclas) + { BSplSLib::BuildCache(0.,0., - 1.,1.,0,0, - uclas - 1,vclas - 1,0,0, - biduflatknots,bidvflatknots, - Poles,Weights, - CachePoles,CacheWeights); + 1.,1.,0,0, + uclas - 1,vclas - 1,0,0, + biduflatknots,bidvflatknots, + Poles,Weights, + CacheRational,CachePolesWeights); } - else { + else + { // BuilCache exige que les resultats soient formates en [MaxCoeff,MinCoeff] - TColgp_Array2OfPnt CPoles (1,vclas, 1, uclas); - TColStd_Array2OfReal CWeights(1,vclas, 1, uclas); - Standard_Integer ii, jj; +//// TColgp_Array2OfPnt CPoles (1,vclas, 1, uclas); +//// TColStd_Array2OfReal CWeights(1,vclas, 1, uclas); + TColStd_Array2OfReal aCachePolesWeights(1, vclas, 1, uclas*4); BSplSLib::BuildCache(0.,0., - 1.,1.,0,0, - uclas - 1,vclas - 1,0,0, - biduflatknots,bidvflatknots, - Poles,Weights, - CPoles,CWeights); - if (&Weights == NULL) { - - for (ii = 1; ii <= uclas; ii++) { - - for (jj = 1; jj <= vclas; jj++) { - CachePoles(ii, jj) = CPoles(jj, ii); - } - } - } - else { - - for (ii = 1; ii <= uclas; ii++) { - - for (jj = 1; jj <= vclas; jj++) { - CachePoles (ii, jj) = CPoles (jj, ii); - CacheWeights(ii, jj) = CWeights(jj, ii); - } - } - } + 1.,1.,0,0, + uclas - 1,vclas - 1,0,0, + biduflatknots,bidvflatknots, + Poles,Weights, + CacheRational,aCachePolesWeights); +//// if (&Weights == NULL) +//// { + Standard_Integer ii, jj, kk; + for (ii = 1; ii <= uclas; ii++) + for (jj = 1; jj <= vclas; jj++) + // The loop boundaries were taken to avoid unnecessary subtractions in array indexes + for (kk = -3; kk <= 0; kk++) + CachePolesWeights(ii, (jj<<2) + kk) = aCachePolesWeights(jj, (ii<<2) + kk); +//// } +//// else +//// { +//// for (ii = 1; ii <= uclas; ii++) +//// for (jj = 1; jj <= vclas; jj++) +//// { +//// CachePoles (ii, jj) = CPoles (jj, ii); +//// CacheWeights(ii, jj) = CWeights(jj, ii); +//// } +//// } } } diff --git a/src/Geom/FILES b/src/Geom/FILES index 644fc70c7a..cbef93b001 100755 --- a/src/Geom/FILES +++ b/src/Geom/FILES @@ -1,6 +1,7 @@ Geom_BSplineCurve_1.cxx Geom_BSplineSurface_1.cxx - +Geom_BSplineCache.lxx +Handle_Geom_BSplineCache.hxx diff --git a/src/Geom/Geom.cdl b/src/Geom/Geom.cdl index aa0217adf9..47c5b55df1 100755 --- a/src/Geom/Geom.cdl +++ b/src/Geom/Geom.cdl @@ -145,11 +145,12 @@ exception UndefinedValue inherits DomainError from Standard; class BezierSurface; class BSplineSurface; class RectangularTrimmedSurface; - + class OffsetSurface; private class OsculatingSurface; + class BSplineCache; -- class PlateSurface ; end Geom; diff --git a/src/Geom/Geom_BSplineCache.cdl b/src/Geom/Geom_BSplineCache.cdl new file mode 100644 index 0000000000..50240d4de7 --- /dev/null +++ b/src/Geom/Geom_BSplineCache.cdl @@ -0,0 +1,108 @@ +-- Created on: 2013-10-17 +-- Created by: AZV +-- Copyright (c) 2013 OPEN CASCADE SAS +-- +-- The content of this file is subject to the Open CASCADE Technology Public +-- License Version 6.5 (the "License"). You may not use the content of this file +-- except in compliance with the License. Please obtain a copy of the License +-- at http://www.opencascade.org and read it completely before using this file. +-- +-- The Initial Developer of the Original Code is Open CASCADE S.A.S., having its +-- main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France. +-- +-- The Original Code and all software distributed under the License is +-- distributed on an "AS IS" basis, without warranty of any kind, and the +-- Initial Developer hereby disclaims all such warranties, including without +-- limitation, any warranties of merchantability, fitness for a particular +-- purpose or non-infringement. Please see the License for the specific terms +-- and conditions governing the rights and limitations under the License. + +class BSplineCache from Geom inherits TShared from MMgt + + ---Purpose: Describes a memory aligned cache array of B-spline and + -- Bezier surfaces poles and weights + +uses Array2OfReal from TColStd + +raises RangeError from Standard, + OutOfRange from Standard, + OutOfMemory from Standard, + DimensionMismatch from Standard + +is + + Create (theRowLower, theRowUpper : Integer from Standard; + theColumnLower, theColumnUpper : Integer from Standard) + returns BSplineCache from Geom + ---Purpose: Creates 2D array of real, + -- address of the first element is aligned on 16 or 32 bytes. + raises RangeError from Standard, + OutOfMemory from Standard; + + Destroy (me : mutable); + ---Level: Advanced + ---Purpose: Frees the allocated memory + ---C++: alias ~ + + ColLength (me) + ---Level: Public + ---Purpose: Return the number of rows in the array. + ---C++: inline + returns Integer from Standard; + + RowLength (me) + ---Level: Public + ---Purpose: Returns the number of columns in the array. + ---C++: inline + returns Integer from Standard; + + Array2 (me) + ---Purpose: Returns the 2D array; the returned array is not modifiable. + ---C++: return const & + ---C++: inline + returns Array2OfReal from TColStd + is static; + + ChangeArray2 (me : mutable) + ---Purpose: Returns a modifiable reference to 2D array. + ---C++: return & + ---C++: inline + returns Array2OfReal from TColStd + is static; + + SetValue (me : mutable; + theRow, theColumn : Integer from Standard; + theValue : Real from Standard) + ---Purpose: Assigns the value to the (, ) item of array. + -- Exceptions + -- Standard_OutOfRange if or lies outside the bounds of this array. + ---C++: inline + raises OutOfRange from Standard; + + + Value (me; + theRow, theCol: Integer from Standard) + ---Level: Public + ---Purpose: Returns the value of the element of index (, ) + ---C++: alias operator() + ---C++: return const & + ---C++: inline + returns Real from Standard + raises OutOfRange from Standard; + + + ChangeValue (me: mutable; + theRow, theColumn: Integer from Standard) + ---Level: Public + ---Purpose: Returns the value of the element of index (, ). + -- Allows to change this value. + ---C++: alias operator() + ---C++: return & + ---C++: inline + returns Real from Standard + raises OutOfRange from Standard; + +fields + myCacheArray : Array2OfReal from TColStd; + +end BSplineCache; diff --git a/src/Geom/Geom_BSplineCache.cxx b/src/Geom/Geom_BSplineCache.cxx new file mode 100644 index 0000000000..4aa36b2b1e --- /dev/null +++ b/src/Geom/Geom_BSplineCache.cxx @@ -0,0 +1,24 @@ +#include + +#ifndef _NCollection_AlignAllocator_HeaderFile +#include +#endif + +IMPLEMENT_STANDARD_RTTIEXT(Geom_BSplineCache, MMgt_TShared) + +Geom_BSplineCache::Geom_BSplineCache(const Standard_Integer theRowLower, const Standard_Integer theRowUpper, + const Standard_Integer theColumnLower, const Standard_Integer theColumnUpper) + : myCacheArray( + *(Standard_Real*)NCollection_AlignAllocator::Allocate( + (theRowUpper-theRowLower+1)*(theColumnUpper-theColumnLower+1)*sizeof(Standard_Real), + DATA_ALIGNMENT), + theRowLower, theRowUpper, theColumnLower, theColumnUpper) + {} + +void Geom_BSplineCache::Destroy() +{ + Standard_Real* aPtr = (Standard_Real*) &(myCacheArray(myCacheArray.LowerCol(), myCacheArray.LowerRow())); + if (aPtr) + NCollection_AlignAllocator::Free(aPtr); +} + diff --git a/src/Geom/Geom_BSplineCache.lxx b/src/Geom/Geom_BSplineCache.lxx new file mode 100644 index 0000000000..f93f2a89a2 --- /dev/null +++ b/src/Geom/Geom_BSplineCache.lxx @@ -0,0 +1,61 @@ +#ifndef _NCollection_AlignAllocator_HeaderFile +#include +#endif + + +Geom_BSplineCache::Geom_BSplineCache(const Standard_Integer theRowLower, const Standard_Integer theRowUpper, + const Standard_Integer theColumnLower, const Standard_Integer theColumnUpper) + : myCacheArray( + *(Standard_Real*)NCollection_AlignAllocator::Allocate( + (theRowUpper-theRowLower+1)*(theColumnUpper-theColumnLower+1)*sizeof(Standard_Real), + DATA_ALIGNMENT), + theRowLower, theRowUpper, theColumnLower, theColumnUpper) + {} + +void Geom_BSplineCache::Destroy() +{ + Standard_Real* aPtr = (Standard_Real*) &(myCacheArray(myCacheArray.LowerCol(), myCacheArray.LowerRow())); + if (aPtr) + NCollection_AlignAllocator::Free(aPtr); +} + +inline const TColStd_Array2OfReal& Geom_BSplineCache::Array2() const +{ + return myCacheArray; +} + +inline TColStd_Array2OfReal& Geom_BSplineCache::ChangeArray2() +{ + return myCacheArray; +} + +inline Standard_Integer Geom_BSplineCache::ColLength() const +{ + return myCacheArray.ColLength(); +} + +inline Standard_Integer Geom_BSplineCache::RowLength() const +{ + return myCacheArray.RowLength(); +} + +inline void Geom_BSplineCache::SetValue(const Standard_Integer theRow, + const Standard_Integer theColumn, + const Standard_Real theValue) +{ + myCacheArray.SetValue(theRow, theColumn, theValue); +} + +inline const Standard_Real& Geom_BSplineCache::Value(const Standard_Integer theRow, + const Standard_Integer theColumn) const +{ + return myCacheArray(theRow, theColumn); +} + +inline Standard_Real& Geom_BSplineCache::ChangeValue(const Standard_Integer theRow, + const Standard_Integer theColumn) +{ + return myCacheArray.ChangeValue(theRow, theColumn); +} + +#include diff --git a/src/Geom/Geom_BSplineSurface.cdl b/src/Geom/Geom_BSplineSurface.cdl index 102802f366..fbec75cda1 100755 --- a/src/Geom/Geom_BSplineSurface.cdl +++ b/src/Geom/Geom_BSplineSurface.cdl @@ -154,7 +154,8 @@ uses Array1OfInteger from TColStd, Curve from Geom, Geometry from Geom, Shape from GeomAbs, - Mutex from Standard + Mutex from Standard, + BSplineCache from Geom raises ConstructionError from Standard, DimensionError from Standard, @@ -1427,7 +1428,7 @@ fields umults : HArray1OfInteger from TColStd; vmults : HArray1OfInteger from TColStd; -- Inplementation of the cache on surfaces - cachepoles : HArray2OfPnt from TColgp; + cachePolesWeights : BSplineCache from Geom; -- Taylor expansion of the poles function, in homogeneous -- form if the curve is rational. The taylor expansion -- is normalized so that the span corresponds to @@ -1458,14 +1459,12 @@ fields -- f (u0,v0) f (u0,v0) ------------- ----------- -- 2 3! -- - -- The size of the array is (1,Max degree) (1, Min degree) + -- The size of the array is (1,Max degree) (1, 4*Min degree) + -- The array consists of numbers which can be composed to 4D row vectors + -- with three components for storing the point and the fourth for the weight -- - cacheweights : HArray2OfReal from TColStd; - -- Taylor expansion of the poles function, in homogeneous - -- form if the curve is rational. The taylor expansion - -- is normalized so that the span corresponds to - -- [0 1]x[0 1]. The Taylor expension of lower degree - -- is stored as consecutive Real in the array as explained above + cacheRational : Boolean; + ucacheparameter : Real ; vcacheparameter : Real ; -- Parameters at which the Taylor expension is stored in diff --git a/src/Geom/Geom_BSplineSurface.cxx b/src/Geom/Geom_BSplineSurface.cxx index e390b7f0d0..2d18e48ffd 100755 --- a/src/Geom/Geom_BSplineSurface.cxx +++ b/src/Geom/Geom_BSplineSurface.cxx @@ -38,6 +38,7 @@ #include #include + //======================================================================= //function : CheckSurfaceData //purpose : Internal use only. @@ -203,10 +204,9 @@ Geom_BSplineSurface::Geom_BSplineSurface vmults->ChangeArray1() = VMults; MinDegree = Min(udeg,vdeg) ; MaxDegree = Max(udeg,vdeg) ; - cachepoles = new TColgp_HArray2OfPnt(1,MaxDegree + 1, - 1,MinDegree + 1) ; - cacheweights.Nullify() ; + cachePolesWeights = new Geom_BSplineCache(1, MaxDegree + 1, 1, (MinDegree + 1) * 4); + ucacheparameter = 0.0e0 ; vcacheparameter = 0.0e0 ; ucachespanlenght = 1.0e0 ; @@ -257,7 +257,7 @@ Geom_BSplineSurface::Geom_BSplineSurface for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) { for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) { if (Weights(i,j) <= gp::Resolution()) - Standard_ConstructionError::Raise("Geom_BSplineSurface"); + Standard_ConstructionError::Raise("Geom_BSplineSurface"); } } @@ -268,19 +268,19 @@ Geom_BSplineSurface::Geom_BSplineSurface // check CheckSurfaceData(Poles, - UKnots , VKnots, - UMults , VMults, - UDegree , VDegree, - UPeriodic, VPeriodic); + UKnots , VKnots, + UMults , VMults, + UDegree , VDegree, + UPeriodic, VPeriodic); // copy arrays poles = new TColgp_HArray2OfPnt(1,Poles.ColLength(), - 1,Poles.RowLength()); + 1,Poles.RowLength()); poles->ChangeArray2() = Poles; weights = new TColStd_HArray2OfReal (1,Poles.ColLength(), - 1,Poles.RowLength()); + 1,Poles.RowLength()); weights->ChangeArray2() = Weights; uknots = new TColStd_HArray1OfReal (1,UKnots.Length()); @@ -296,12 +296,9 @@ Geom_BSplineSurface::Geom_BSplineSurface vmults->ChangeArray1() = VMults; MinDegree = Min(udeg,vdeg) ; MaxDegree = Max(udeg,vdeg) ; - cachepoles = new TColgp_HArray2OfPnt(1,MaxDegree + 1, - 1,MinDegree + 1) ; - if (urational || vrational) { - cacheweights = new TColStd_HArray2OfReal (1,MaxDegree + 1, - 1,MinDegree + 1); - } + + cachePolesWeights = new Geom_BSplineCache(1, MaxDegree + 1, 1, (MinDegree + 1) * 4); + ucacheparameter = 0.0e0 ; vcacheparameter = 0.0e0 ; ucachespanlenght = 1.0e0 ; @@ -1380,49 +1377,31 @@ void Geom_BSplineSurface::ValidateCache(const Standard_Real Uparameter, { Standard_Real NewParameter ; Standard_Integer LocalIndex = 0 ; - Standard_Integer MinDegree, - MaxDegree ; - // - // check if the degree did not change - // + Standard_Integer MinDegree, MaxDegree ; + // check if the degree is changed MinDegree = Min(udeg,vdeg) ; MaxDegree = Max(udeg,vdeg) ; - if (cachepoles->ColLength() < MaxDegree + 1 || - cachepoles->RowLength() < MinDegree + 1) { - cachepoles = new TColgp_HArray2OfPnt(1,MaxDegree + 1, - 1,MinDegree + 1); - } - // - // Verif + poussee pour les poids - // - if (urational || vrational) { - if (cacheweights.IsNull()) { - cacheweights = new TColStd_HArray2OfReal(1,MaxDegree + 1, - 1,MinDegree + 1); - } - else { - if (cacheweights->ColLength() < MaxDegree + 1 || - cacheweights->RowLength() < MinDegree + 1) { - cacheweights = new TColStd_HArray2OfReal(1,MaxDegree + 1, - 1,MinDegree + 1); - } - } + if (cachePolesWeights->ColLength() < (MaxDegree + 1) || + cachePolesWeights->RowLength() < (MinDegree + 1) * 4) + { + cachePolesWeights = new Geom_BSplineCache(1, MaxDegree + 1, 1, (MinDegree + 1) * 4); } BSplCLib::LocateParameter(udeg, - (ufknots->Array1()), - (BSplCLib::NoMults()), - Uparameter, - uperiodic, - LocalIndex, - NewParameter); - ucachespanindex = LocalIndex ; - if (Uparameter == ufknots->Value(LocalIndex + 1)) { - + (ufknots->Array1()), + (BSplCLib::NoMults()), + Uparameter, + uperiodic, + LocalIndex, + NewParameter); + ucachespanindex = LocalIndex; + if (Uparameter == ufknots->Value(LocalIndex + 1)) + { LocalIndex += 1 ; - ucacheparameter = ufknots->Value(LocalIndex) ; - if (LocalIndex == ufknots->Upper() - udeg) { + ucacheparameter = ufknots->Value(LocalIndex); + if (LocalIndex == ufknots->Upper() - udeg) + { // // for the last span if the parameter is outside of // the domain of the curve than use the last knot @@ -1431,31 +1410,34 @@ void Geom_BSplineSurface::ValidateCache(const Standard_Real Uparameter, // the IsCacheValid will know for sure we are extending // the Bspline // - - ucachespanlenght = ufknots->Value(LocalIndex - 1) - ucacheparameter ; + ucachespanlenght = ufknots->Value(LocalIndex - 1) - ucacheparameter; } - else { - ucachespanlenght = ufknots->Value(LocalIndex + 1) - ucacheparameter ; - } - } - else { - ucacheparameter = ufknots->Value(LocalIndex) ; - ucachespanlenght = ufknots->Value(LocalIndex + 1) - ucacheparameter ; + else + { + ucachespanlenght = ufknots->Value(LocalIndex + 1) - ucacheparameter; } + } + else + { + ucacheparameter = ufknots->Value(LocalIndex); + ucachespanlenght = ufknots->Value(LocalIndex + 1) - ucacheparameter; + } LocalIndex = 0 ; BSplCLib::LocateParameter(vdeg, - (vfknots->Array1()), - (BSplCLib::NoMults()), - Vparameter, - vperiodic, - LocalIndex, - NewParameter); - vcachespanindex = LocalIndex ; - if (Vparameter == vfknots->Value(LocalIndex + 1)) { - LocalIndex += 1 ; - vcacheparameter = vfknots->Value(LocalIndex) ; - if (LocalIndex == vfknots->Upper() - vdeg) { + (vfknots->Array1()), + (BSplCLib::NoMults()), + Vparameter, + vperiodic, + LocalIndex, + NewParameter); + vcachespanindex = LocalIndex; + if (Vparameter == vfknots->Value(LocalIndex + 1)) + { + LocalIndex += 1; + vcacheparameter = vfknots->Value(LocalIndex); + if (LocalIndex == vfknots->Upper() - vdeg) + { // // for the last span if the parameter is outside of // the domain of the curve than use the last knot @@ -1464,57 +1446,60 @@ void Geom_BSplineSurface::ValidateCache(const Standard_Real Uparameter, // the IsCacheValid will know for sure we are extending // the Bspline // - - vcachespanlenght = vfknots->Value(LocalIndex - 1) - vcacheparameter ; + vcachespanlenght = vfknots->Value(LocalIndex - 1) - vcacheparameter; } - else { - vcachespanlenght = vfknots->Value(LocalIndex + 1) - vcacheparameter ; + else + { + vcachespanlenght = vfknots->Value(LocalIndex + 1) - vcacheparameter; } } - else { - vcacheparameter = vfknots->Value(LocalIndex) ; - vcachespanlenght = vfknots->Value(LocalIndex + 1) - vcacheparameter ; + else + { + vcacheparameter = vfknots->Value(LocalIndex); + vcachespanlenght = vfknots->Value(LocalIndex + 1) - vcacheparameter; } - - Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2, - uspanlenght_11 = ucachespanlenght/2, - vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2, - vspanlenght_11 = vcachespanlenght/2 ; - if (urational || vrational) { + + Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2; + Standard_Real uspanlenght_11 = ucachespanlenght/2; + Standard_Real vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2; + Standard_Real vspanlenght_11 = vcachespanlenght/2; + if (urational || vrational) + { BSplSLib::BuildCache(uparameter_11, - vparameter_11, - uspanlenght_11, - vspanlenght_11, - uperiodic, - vperiodic, - udeg, - vdeg, - ucachespanindex, - vcachespanindex, - (ufknots->Array1()), - (vfknots->Array1()), - poles->Array2(), - weights->Array2(), - cachepoles->ChangeArray2(), - cacheweights->ChangeArray2()) ; + vparameter_11, + uspanlenght_11, + vspanlenght_11, + uperiodic, + vperiodic, + udeg, + vdeg, + ucachespanindex, + vcachespanindex, + (ufknots->Array1()), + (vfknots->Array1()), + poles->Array2(), + weights->Array2(), + cacheRational, + cachePolesWeights->ChangeArray2()); } - else { + else + { BSplSLib::BuildCache(uparameter_11, - vparameter_11, - uspanlenght_11, - vspanlenght_11, - uperiodic, - vperiodic, - udeg, - vdeg, - ucachespanindex, - vcachespanindex, - (ufknots->Array1()), - (vfknots->Array1()), - poles->Array2(), - *((TColStd_Array2OfReal*) NULL), - cachepoles->ChangeArray2(), - *((TColStd_Array2OfReal*) NULL)) ; + vparameter_11, + uspanlenght_11, + vspanlenght_11, + uperiodic, + vperiodic, + udeg, + vdeg, + ucachespanindex, + vcachespanindex, + (ufknots->Array1()), + (vfknots->Array1()), + poles->Array2(), + *((TColStd_Array2OfReal*) NULL), + cacheRational, + cachePolesWeights->ChangeArray2()); } validcache = 1 ; } diff --git a/src/Geom/Geom_BSplineSurface_1.cxx b/src/Geom/Geom_BSplineSurface_1.cxx index 864c997f00..75024efa36 100755 --- a/src/Geom/Geom_BSplineSurface_1.cxx +++ b/src/Geom/Geom_BSplineSurface_1.cxx @@ -123,39 +123,23 @@ void Geom_BSplineSurface::D0(const Standard_Real U, Standard_Mutex::Sentry aSentry(MySurface->myMutex); if(!IsCacheValid(new_u, new_v)) - MySurface->ValidateCache(new_u, new_v); - - Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2, - uspanlenght_11 = ucachespanlenght/2, - vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2, - vspanlenght_11 = vcachespanlenght/2 ; - if (cacheweights.IsNull()) { - - BSplSLib::CacheD0(new_u, - new_v, - udeg, - vdeg, - uparameter_11, - vparameter_11, - uspanlenght_11, - vspanlenght_11, - cachepoles->Array2(), - *((TColStd_Array2OfReal*) NULL), - P) ; - } - else { - BSplSLib::CacheD0(new_u, - new_v, - udeg, - vdeg, - uparameter_11, - vparameter_11, - uspanlenght_11, - vspanlenght_11, - cachepoles->Array2(), - cacheweights->Array2(), - P) ; - } + MySurface->ValidateCache(new_u, new_v); + + Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2; + Standard_Real uspanlenght_11 = ucachespanlenght/2; + Standard_Real vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2; + Standard_Real vspanlenght_11 = vcachespanlenght/2; + BSplSLib::CacheD0(new_u, + new_v, + udeg, + vdeg, + uparameter_11, + vparameter_11, + uspanlenght_11, + vspanlenght_11, + cachePolesWeights->Array2(), + cacheRational, + P); } //======================================================================= @@ -176,45 +160,26 @@ void Geom_BSplineSurface::D1(const Standard_Real U, Standard_Mutex::Sentry aSentry(MySurface->myMutex); if(!IsCacheValid(new_u, new_v)) - MySurface->ValidateCache(new_u, new_v); + MySurface->ValidateCache(new_u, new_v); - Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2, - uspanlenght_11 = ucachespanlenght/2, - vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2, - vspanlenght_11 = vcachespanlenght/2 ; + Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2; + Standard_Real uspanlenght_11 = ucachespanlenght/2; + Standard_Real vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2; + Standard_Real vspanlenght_11 = vcachespanlenght/2; - if (cacheweights.IsNull()) { - - BSplSLib::CacheD1(new_u, - new_v, - udeg, - vdeg, - uparameter_11, - vparameter_11, - uspanlenght_11, - vspanlenght_11, - cachepoles->Array2(), - *((TColStd_Array2OfReal*) NULL), - P, - D1U, - D1V) ; - } - else { - - BSplSLib::CacheD1(new_u, - new_v, - udeg, - vdeg, - uparameter_11, - vparameter_11, - uspanlenght_11, - vspanlenght_11, - cachepoles->Array2(), - cacheweights->Array2(), - P, - D1U, - D1V) ; - } + BSplSLib::CacheD1(new_u, + new_v, + udeg, + vdeg, + uparameter_11, + vparameter_11, + uspanlenght_11, + vspanlenght_11, + cachePolesWeights->Array2(), + cacheRational, + P, + D1U, + D1V); } //======================================================================= @@ -223,13 +188,13 @@ void Geom_BSplineSurface::D1(const Standard_Real U, //======================================================================= void Geom_BSplineSurface::D2 (const Standard_Real U, - const Standard_Real V, - gp_Pnt& P, - gp_Vec& D1U, - gp_Vec& D1V, - gp_Vec& D2U, - gp_Vec& D2V, - gp_Vec& D2UV) const + const Standard_Real V, + gp_Pnt& P, + gp_Vec& D1U, + gp_Vec& D1V, + gp_Vec& D2U, + gp_Vec& D2V, + gp_Vec& D2UV) const { Standard_Real new_u(U), new_v(V); PeriodicNormalization(new_u, new_v); @@ -238,49 +203,30 @@ void Geom_BSplineSurface::D2 (const Standard_Real U, Standard_Mutex::Sentry aSentry(MySurface->myMutex); if(!IsCacheValid(new_u, new_v)) - MySurface->ValidateCache(new_u, new_v); - - Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2, - uspanlenght_11 = ucachespanlenght/2, - vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2, - vspanlenght_11 = vcachespanlenght/2 ; - if (cacheweights.IsNull()) { - BSplSLib::CacheD2(new_u, - new_v, - udeg, - vdeg, - uparameter_11, - vparameter_11, - uspanlenght_11, - vspanlenght_11, - cachepoles->Array2(), - *((TColStd_Array2OfReal*) NULL), - P, - D1U, - D1V, - D2U, - D2UV, - D2V); - } - else { - BSplSLib::CacheD2(new_u, - new_v, - udeg, - vdeg, - uparameter_11, - vparameter_11, - uspanlenght_11, - vspanlenght_11, - cachepoles->Array2(), - cacheweights->Array2(), - P, - D1U, - D1V, - D2U, - D2UV, - D2V); - } - } + MySurface->ValidateCache(new_u, new_v); + + Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2; + Standard_Real uspanlenght_11 = ucachespanlenght/2; + Standard_Real vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2; + Standard_Real vspanlenght_11 = vcachespanlenght/2; + + BSplSLib::CacheD2(new_u, + new_v, + udeg, + vdeg, + uparameter_11, + vparameter_11, + uspanlenght_11, + vspanlenght_11, + cachePolesWeights->Array2(), + cacheRational, + P, + D1U, + D1V, + D2U, + D2UV, + D2V); +} //======================================================================= //function : D3 @@ -301,9 +247,9 @@ void Geom_BSplineSurface::D3 (const Standard_Real U, gp_Vec& D3UVV) const { BSplSLib::D3(U,V,0,0,POLES,WEIGHTS,UFKNOTS,VFKNOTS,FMULTS,FMULTS, - udeg,vdeg,urational,vrational,uperiodic,vperiodic, - P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV); - } + udeg,vdeg,urational,vrational,uperiodic,vperiodic, + P,D1U,D1V,D2U,D2V,D2UV,D3U,D3V,D3UUV,D3UVV); +} //======================================================================= //function : DN diff --git a/src/Geom/Geom_BezierSurface.cdl b/src/Geom/Geom_BezierSurface.cdl index 0704db91da..6ea9042002 100755 --- a/src/Geom/Geom_BezierSurface.cdl +++ b/src/Geom/Geom_BezierSurface.cdl @@ -105,7 +105,8 @@ uses Array1OfReal from TColStd, Vec from gp, Curve from Geom, Geometry from Geom, - Shape from GeomAbs + Shape from GeomAbs, + BSplineCache from Geom raises ConstructionError from Standard, @@ -668,9 +669,10 @@ is -- Bezier surface. This value is 25. - Create (SurfacePoles, SurfaceCoefficients : HArray2OfPnt from TColgp; - PoleWeights, CoefficientWeights : HArray2OfReal from TColStd; - IsURational, IsVRational : Boolean) + Create (SurfacePoles : HArray2OfPnt from TColgp; + PoleWeights : HArray2OfReal from TColStd; + SurfaceCoeffsAndWeights : BSplineCache from Geom; + IsURational, IsVRational : Boolean) returns mutable BezierSurface is private; @@ -723,8 +725,8 @@ fields poles : HArray2OfPnt from TColgp; weights : HArray2OfReal from TColStd; - coeffs : HArray2OfPnt from TColgp; - wcoeffs : HArray2OfReal from TColStd; + coeffsweights : BSplineCache from Geom; + cacherational : Boolean; ucacheparameter : Real ; vcacheparameter : Real ; diff --git a/src/Geom/Geom_BezierSurface.cxx b/src/Geom/Geom_BezierSurface.cxx index 03c75959fc..854ee84aaf 100755 --- a/src/Geom/Geom_BezierSurface.cxx +++ b/src/Geom/Geom_BezierSurface.cxx @@ -471,14 +471,13 @@ Geom_BezierSurface::Geom_BezierSurface //purpose : //======================================================================= -Geom_BezierSurface::Geom_BezierSurface - (const Handle(TColgp_HArray2OfPnt)& SurfacePoles, - const Handle(TColgp_HArray2OfPnt)& SurfaceCoefs, - const Handle(TColStd_HArray2OfReal)& PoleWeights, - const Handle(TColStd_HArray2OfReal)& CoefWeights, - const Standard_Boolean IsURational, - const Standard_Boolean IsVRational) -:maxderivinvok(Standard_False) +Geom_BezierSurface::Geom_BezierSurface( + const Handle(TColgp_HArray2OfPnt)& SurfacePoles, + const Handle(TColStd_HArray2OfReal)& PoleWeights, + const Handle(Geom_BSplineCache)& SurfaceCoeffsAndWeights, + const Standard_Boolean IsURational, + const Standard_Boolean IsVRational) + :maxderivinvok(Standard_False) { urational = IsURational; vrational = IsVRational; @@ -490,21 +489,17 @@ Geom_BezierSurface::Geom_BezierSurface Standard_Integer NbUPoles = SurfacePoles->ColLength(); Standard_Integer NbVPoles = SurfacePoles->RowLength(); - poles = new TColgp_HArray2OfPnt (1,NbUPoles, - 1,NbVPoles) ; + poles = new TColgp_HArray2OfPnt(1, NbUPoles, 1, NbVPoles); poles->ChangeArray2() = SurfacePoles->Array2(); - coeffs = new TColgp_HArray2OfPnt (1,SurfaceCoefs->ColLength(), - 1,SurfaceCoefs->RowLength()) ; - coeffs->ChangeArray2() = SurfaceCoefs->Array2(); + coeffsweights = new Geom_BSplineCache(1, SurfaceCoeffsAndWeights->ColLength(), + 1, SurfaceCoeffsAndWeights->RowLength()); + coeffsweights->ChangeArray2() = SurfaceCoeffsAndWeights->Array2(); - if ( urational || vrational) { - weights = new TColStd_HArray2OfReal (1,NbUPoles,1,NbVPoles); + if ( urational || vrational) + { + weights = new TColStd_HArray2OfReal(1,NbUPoles,1,NbVPoles); weights->ChangeArray2() = PoleWeights->Array2(); - - wcoeffs = new TColStd_HArray2OfReal (1,SurfaceCoefs->ColLength(), - 1,SurfaceCoefs->RowLength()) ; - wcoeffs->ChangeArray2() = CoefWeights->Array2(); } } @@ -540,7 +535,7 @@ void Geom_BezierSurface::ExchangeUV () TColgp_Array2OfPnt& snpoles = npoles->ChangeArray2(); TColStd_Array2OfReal& snweights = nweights->ChangeArray2(); - + Standard_Integer i, j; for (i = LC; i <= UC; i++) { for (j = LR; j <= UR; j++) { @@ -555,8 +550,10 @@ void Geom_BezierSurface::ExchangeUV () Standard_Boolean temp = urational; urational = vrational; vrational = temp; - coeffs = new TColgp_HArray2OfPnt (LC, UC, LR, UR); - wcoeffs = new TColStd_HArray2OfReal (LC, UC, LR, UR); + + // coeffsweights is an array of components of 4D row vector + // (three for pole point and fourth for weight) + coeffsweights = new Geom_BSplineCache(LC, UC, LR, UR<<2); UpdateCoefficients(); } @@ -648,43 +645,46 @@ void Geom_BezierSurface::Increase (const Standard_Integer UDeg, //purpose : //======================================================================= -void Geom_BezierSurface::InsertPoleColAfter - (const Standard_Integer VIndex, - const TColgp_Array1OfPnt& CPoles) +void Geom_BezierSurface::InsertPoleColAfter( + const Standard_Integer VIndex, + const TColgp_Array1OfPnt& CPoles) { const TColgp_Array2OfPnt & Poles = poles->Array2(); - if (VIndex < 1 || VIndex > Poles.RowLength()) Standard_OutOfRange::Raise(); - if (CPoles.Length() != Poles.ColLength()) { + if (VIndex < 1 || VIndex > Poles.RowLength()) + Standard_OutOfRange::Raise(); + if (CPoles.Length() != Poles.ColLength()) Standard_ConstructionError::Raise(); - } Handle(TColgp_HArray2OfPnt) npoles = - new TColgp_HArray2OfPnt(1,poles->ColLength(),1,poles->RowLength()+1); + new TColgp_HArray2OfPnt(1, poles->ColLength(), 1, poles->RowLength()+1); Handle(TColStd_HArray2OfReal) nweights; - if (urational || vrational) { + if (urational || vrational) + { nweights = - new TColStd_HArray2OfReal(1,poles->ColLength(),1,poles->RowLength()+1); + new TColStd_HArray2OfReal(1, poles->ColLength(), 1, poles->RowLength()+1); - TColStd_Array1OfReal CWeights(nweights->LowerRow(),nweights->UpperRow()); + TColStd_Array1OfReal CWeights(nweights->LowerRow(), nweights->UpperRow()); CWeights.Init(1.); - AddRatPoleCol (poles->Array2(), weights->Array2(), - CPoles, CWeights, VIndex, - npoles->ChangeArray2(), nweights->ChangeArray2()); + AddRatPoleCol(poles->Array2(), weights->Array2(), + CPoles, CWeights, VIndex, + npoles->ChangeArray2(), nweights->ChangeArray2()); } - else { - AddPoleCol (poles->Array2(), - CPoles, VIndex, - npoles->ChangeArray2()); + else + { + AddPoleCol(poles->Array2(), + CPoles, VIndex, + npoles->ChangeArray2()); } poles = npoles; weights = nweights; - coeffs = new TColgp_HArray2OfPnt(1,poles->ColLength(), - 1,poles->RowLength()); - wcoeffs = new TColStd_HArray2OfReal(1,poles->ColLength(), - 1,poles->RowLength()); + + // coeffsweights is an array of components of 4D row vector + // (three for pole point and fourth for weight) + coeffsweights = new Geom_BSplineCache(1, poles->ColLength(), + 1, poles->RowLength()<<2); UpdateCoefficients(); } @@ -693,44 +693,45 @@ void Geom_BezierSurface::InsertPoleColAfter //purpose : //======================================================================= -void Geom_BezierSurface::InsertPoleColAfter - (const Standard_Integer VIndex, - const TColgp_Array1OfPnt& CPoles, - const TColStd_Array1OfReal& CPoleWeights) +void Geom_BezierSurface::InsertPoleColAfter( + const Standard_Integer VIndex, + const TColgp_Array1OfPnt& CPoles, + const TColStd_Array1OfReal& CPoleWeights) { const TColgp_Array2OfPnt & Poles = poles->Array2(); - if (VIndex < 1 || VIndex > Poles.RowLength()) Standard_OutOfRange::Raise(); - if (CPoles.Length() != Poles.ColLength() || - CPoleWeights.Length() != CPoles.Length()) { + if (VIndex < 1 || VIndex > Poles.RowLength()) + Standard_OutOfRange::Raise(); + if (CPoles.Length() != Poles.ColLength() || + CPoleWeights.Length() != CPoles.Length()) Standard_ConstructionError::Raise(); - } + Standard_Integer Index = CPoleWeights.Lower(); - while (Index <= CPoleWeights.Upper()) { - if (CPoleWeights (Index) <= gp::Resolution()) { + while (Index <= CPoleWeights.Upper()) + { + if (CPoleWeights (Index) <= gp::Resolution()) Standard_ConstructionError::Raise(); - } Index++; } Handle(TColgp_HArray2OfPnt) npoles = - new TColgp_HArray2OfPnt(1,poles->ColLength(),1,poles->RowLength()+1); - + new TColgp_HArray2OfPnt(1, poles->ColLength(), 1, poles->RowLength()+1); + Handle(TColStd_HArray2OfReal) nweights = - new TColStd_HArray2OfReal(1,poles->ColLength(),1,poles->RowLength()+1); - - AddRatPoleCol (poles->Array2(), weights->Array2(), - CPoles, CPoleWeights, VIndex, - npoles->ChangeArray2(), nweights->ChangeArray2()); + new TColStd_HArray2OfReal(1, poles->ColLength(), 1, poles->RowLength()+1); + + AddRatPoleCol(poles->Array2(), weights->Array2(), + CPoles, CPoleWeights, VIndex, + npoles->ChangeArray2(), nweights->ChangeArray2()); poles = npoles; weights = nweights; - coeffs = new TColgp_HArray2OfPnt(1,poles->ColLength(), - 1,poles->RowLength()); - wcoeffs = new TColStd_HArray2OfReal(1,poles->ColLength(), - 1,poles->RowLength()); - Rational(weights->Array2(), urational, vrational); + // coeffsweights is an array of components of 4D row vector + // (three for pole point and fourth for weight) + coeffsweights = new Geom_BSplineCache(1, poles->ColLength(), + 1, poles->RowLength()<<2); + Rational(weights->Array2(), urational, vrational); UpdateCoefficients(); } @@ -763,44 +764,48 @@ void Geom_BezierSurface::InsertPoleColBefore //purpose : //======================================================================= -void Geom_BezierSurface::InsertPoleRowAfter (const Standard_Integer UIndex, - const TColgp_Array1OfPnt& CPoles) +void Geom_BezierSurface::InsertPoleRowAfter( + const Standard_Integer UIndex, + const TColgp_Array1OfPnt& CPoles) { const TColgp_Array2OfPnt & Poles = poles->Array2(); - if (UIndex < 1 || UIndex > Poles.ColLength()) Standard_OutOfRange::Raise(); - if (CPoles.Length() != Poles.RowLength()) { + if (UIndex < 1 || UIndex > Poles.ColLength()) + Standard_OutOfRange::Raise(); + if (CPoles.Length() != Poles.RowLength()) Standard_ConstructionError::Raise(); - } Handle(TColgp_HArray2OfPnt) npoles = - new TColgp_HArray2OfPnt(1,poles->ColLength()+1,1,poles->RowLength()); + new TColgp_HArray2OfPnt(1, poles->ColLength()+1, 1, poles->RowLength()); Handle(TColStd_HArray2OfReal) nweights; - if (urational || vrational) { + if (urational || vrational) + { nweights = - new TColStd_HArray2OfReal(1,poles->ColLength()+1,1,poles->RowLength()); + new TColStd_HArray2OfReal(1, poles->ColLength()+1, 1, poles->RowLength()); -// TColStd_Array1OfReal CWeights(nweights->LowerCol(),nweights->UpperCol(), -// 1.0); ??????????? - TColStd_Array1OfReal CWeights(1.0, - nweights->LowerCol(),nweights->UpperCol()); + // TColStd_Array1OfReal CWeights(nweights->LowerCol(),nweights->UpperCol(), + // 1.0); ??????????? + TColStd_Array1OfReal CWeights(1.0, + nweights->LowerCol(), nweights->UpperCol()); - AddRatPoleRow (poles->Array2(), weights->Array2(), - CPoles, CWeights, UIndex, - npoles->ChangeArray2(), nweights->ChangeArray2()); + AddRatPoleRow(poles->Array2(), weights->Array2(), + CPoles, CWeights, UIndex, + npoles->ChangeArray2(), nweights->ChangeArray2()); } - else { - AddPoleRow (poles->Array2(), - CPoles, UIndex, - npoles->ChangeArray2()); + else + { + AddPoleRow(poles->Array2(), + CPoles, UIndex, + npoles->ChangeArray2()); } poles = npoles; weights = nweights; - coeffs = new TColgp_HArray2OfPnt(1,poles->ColLength(), - 1,poles->RowLength()); - wcoeffs = new TColStd_HArray2OfReal(1,poles->ColLength(), - 1,poles->RowLength()); + + // coeffsweights is an array of components of 4D row vector + // (three for pole point and fourth for weight) + coeffsweights = new Geom_BSplineCache(1, poles->ColLength(), + 1, poles->RowLength()<<2); UpdateCoefficients(); } @@ -810,44 +815,45 @@ void Geom_BezierSurface::InsertPoleRowAfter (const Standard_Integer UIndex, //purpose : //======================================================================= -void Geom_BezierSurface::InsertPoleRowAfter - (const Standard_Integer UIndex, - const TColgp_Array1OfPnt& CPoles, - const TColStd_Array1OfReal& CPoleWeights) +void Geom_BezierSurface::InsertPoleRowAfter( + const Standard_Integer UIndex, + const TColgp_Array1OfPnt& CPoles, + const TColStd_Array1OfReal& CPoleWeights) { const TColgp_Array2OfPnt & Poles = poles->Array2(); - if (UIndex < 1 || UIndex > Poles.ColLength()) Standard_OutOfRange::Raise(); + if (UIndex < 1 || UIndex > Poles.ColLength()) + Standard_OutOfRange::Raise(); if (CPoles.Length() != Poles.RowLength() || - CPoleWeights.Length() != CPoles.Length()) { + CPoleWeights.Length() != CPoles.Length()) Standard_ConstructionError::Raise(); - } + Standard_Integer Index = CPoleWeights.Lower(); - while (Index <= CPoleWeights.Upper()) { - if (CPoleWeights(Index) <= gp::Resolution()) { + while (Index <= CPoleWeights.Upper()) + { + if (CPoleWeights(Index) <= gp::Resolution()) Standard_ConstructionError::Raise(); - } Index++; } Handle(TColgp_HArray2OfPnt) npoles = - new TColgp_HArray2OfPnt(1,poles->ColLength()+1,1,poles->RowLength()); - + new TColgp_HArray2OfPnt(1, poles->ColLength()+1, 1, poles->RowLength()); + Handle(TColStd_HArray2OfReal) nweights = - new TColStd_HArray2OfReal(1,poles->ColLength()+1,1,poles->RowLength()); - - AddRatPoleCol (poles->Array2(), weights->Array2(), - CPoles, CPoleWeights, UIndex, - npoles->ChangeArray2(), nweights->ChangeArray2()); + new TColStd_HArray2OfReal(1, poles->ColLength()+1, 1, poles->RowLength()); + + AddRatPoleCol(poles->Array2(), weights->Array2(), + CPoles, CPoleWeights, UIndex, + npoles->ChangeArray2(), nweights->ChangeArray2()); poles = npoles; weights = nweights; - coeffs = new TColgp_HArray2OfPnt(1,poles->ColLength(), - 1,poles->RowLength()); - wcoeffs = new TColStd_HArray2OfReal(1,poles->ColLength(), - 1,poles->RowLength()); - Rational(weights->Array2(), urational, vrational); + // coeffsweights is an array of components of 4D row vector + // (three for pole point and fourth for weight) + coeffsweights = new Geom_BSplineCache(1, poles->ColLength(), + 1, poles->RowLength()<<2); + Rational(weights->Array2(), urational, vrational); UpdateCoefficients(); } @@ -880,38 +886,43 @@ void Geom_BezierSurface::InsertPoleRowBefore //purpose : //======================================================================= -void Geom_BezierSurface::RemovePoleCol (const Standard_Integer VIndex) +void Geom_BezierSurface::RemovePoleCol(const Standard_Integer VIndex) { const TColgp_Array2OfPnt & Poles = poles->Array2(); - if (VIndex < 1 || VIndex > Poles.RowLength()) Standard_OutOfRange::Raise(); - if (Poles.RowLength() <= 2) Standard_ConstructionError::Raise(); + if (VIndex < 1 || VIndex > Poles.RowLength()) + Standard_OutOfRange::Raise(); + if (Poles.RowLength() <= 2) + Standard_ConstructionError::Raise(); Handle(TColgp_HArray2OfPnt) npoles = - new TColgp_HArray2OfPnt(1,poles->ColLength(),1,poles->RowLength()-1); + new TColgp_HArray2OfPnt(1, poles->ColLength(), 1, poles->RowLength()-1); Handle(TColStd_HArray2OfReal) nweights; - if (urational || vrational) { + if (urational || vrational) + { nweights = - new TColStd_HArray2OfReal(1,poles->ColLength(),1,poles->RowLength()-1); + new TColStd_HArray2OfReal(1, poles->ColLength(), 1, poles->RowLength()-1); - DeleteRatPoleCol (poles->Array2(), weights->Array2(), - VIndex, - npoles->ChangeArray2(), nweights->ChangeArray2()); + DeleteRatPoleCol(poles->Array2(), weights->Array2(), + VIndex, + npoles->ChangeArray2(), nweights->ChangeArray2()); // Mise a jour de la rationalite Rational(nweights->Array2(), urational, vrational); } - else { - DeletePoleCol (poles->Array2(), - VIndex, - npoles->ChangeArray2()); + else + { + DeletePoleCol(poles->Array2(), + VIndex, + npoles->ChangeArray2()); } poles = npoles; weights = nweights; - coeffs = new TColgp_HArray2OfPnt(1,poles->ColLength(), - 1,poles->RowLength()); - wcoeffs = new TColStd_HArray2OfReal(1,poles->ColLength(), - 1,poles->RowLength()); + + // coeffsweights is an array of components of 4D row vector + // (three for pole point and fourth for weight) + coeffsweights = new Geom_BSplineCache(1, poles->ColLength(), + 1, poles->RowLength()<<2); UpdateCoefficients(); } @@ -920,39 +931,44 @@ void Geom_BezierSurface::RemovePoleCol (const Standard_Integer VIndex) //purpose : //======================================================================= -void Geom_BezierSurface::RemovePoleRow (const Standard_Integer UIndex) +void Geom_BezierSurface::RemovePoleRow(const Standard_Integer UIndex) { const TColgp_Array2OfPnt & Poles = poles->Array2(); - if (UIndex < 1 || UIndex > Poles.ColLength()) Standard_OutOfRange::Raise(); - if (Poles.ColLength() <= 2) Standard_ConstructionError::Raise(); + if (UIndex < 1 || UIndex > Poles.ColLength()) + Standard_OutOfRange::Raise(); + if (Poles.ColLength() <= 2) + Standard_ConstructionError::Raise(); Handle(TColgp_HArray2OfPnt) npoles = - new TColgp_HArray2OfPnt(1,poles->ColLength()-1,1,poles->RowLength()); + new TColgp_HArray2OfPnt(1, poles->ColLength()-1, 1, poles->RowLength()); Handle(TColStd_HArray2OfReal) nweights; - if (urational || vrational) { + if (urational || vrational) + { nweights = - new TColStd_HArray2OfReal(1,poles->ColLength()-1,1,poles->RowLength()); + new TColStd_HArray2OfReal(1, poles->ColLength()-1, 1, poles->RowLength()); - DeleteRatPoleRow (poles->Array2(), weights->Array2(), - UIndex, - npoles->ChangeArray2(), nweights->ChangeArray2()); + DeleteRatPoleRow(poles->Array2(), weights->Array2(), + UIndex, + npoles->ChangeArray2(), nweights->ChangeArray2()); // Mise a jour de la rationalite Rational(nweights->Array2(), urational, vrational); } - else { - DeletePoleRow (poles->Array2(), - UIndex, - npoles->ChangeArray2()); + else + { + DeletePoleRow(poles->Array2(), + UIndex, + npoles->ChangeArray2()); } poles = npoles; weights = nweights; - coeffs = new TColgp_HArray2OfPnt(1,poles->ColLength(), - 1,poles->RowLength()); - wcoeffs = new TColStd_HArray2OfReal(1,poles->ColLength(), - 1,poles->RowLength()); + + // coeffsweights is an array of components of 4D row vector + // (three for pole point and fourth for weight) + coeffsweights = new Geom_BSplineCache(1, poles->ColLength(), + 1, poles->RowLength()<<2); UpdateCoefficients(); } @@ -961,63 +977,78 @@ void Geom_BezierSurface::RemovePoleRow (const Standard_Integer UIndex) //purpose : //======================================================================= -void Geom_BezierSurface::Segment - (const Standard_Real U1, - const Standard_Real U2, - const Standard_Real V1, - const Standard_Real V2) +void Geom_BezierSurface::Segment( + const Standard_Real U1, + const Standard_Real U2, + const Standard_Real V1, + const Standard_Real V2) { Standard_Boolean rat = (urational || vrational); Handle(TColgp_HArray2OfPnt) Coefs; Handle(TColStd_HArray2OfReal) WCoefs; - - if (validcache == 0) UpdateCoefficients(0., 0.); + if (validcache == 0) + UpdateCoefficients(); + + Standard_Integer ii, jj; + Coefs = new TColgp_HArray2OfPnt(1, UDegree()+1, 1, VDegree()+1); + if (rat) + WCoefs = new TColStd_HArray2OfReal(1, UDegree()+1, 1, VDegree()+1); // Attention si udeg <= vdeg u et v sont intervertis // dans les coeffs, il faut donc tout transposer. - if(UDegree() <= VDegree()) { - Standard_Integer ii, jj; - Coefs = new (TColgp_HArray2OfPnt)(1,UDegree()+1,1,VDegree()+1); - if (rat) { - WCoefs = new (TColStd_HArray2OfReal)(1,UDegree()+1,1,VDegree()+1); - } + if(UDegree() <= VDegree()) + { for (ii=1; ii<=UDegree()+1; ii++) - for (jj=1; jj<=VDegree()+1; jj++) { - Coefs->SetValue(ii, jj, coeffs->Value(jj,ii)); - if (rat) WCoefs->SetValue(ii, jj, wcoeffs->Value(jj,ii)); + for (jj=1; jj<=VDegree()+1; jj++) + { + Coefs->SetValue(ii, jj, gp_Pnt(coeffsweights->Value(jj, (ii<<2)-3), + coeffsweights->Value(jj, (ii<<2)-2), + coeffsweights->Value(jj, (ii<<2)-1))); + if (rat) + WCoefs->SetValue(ii, jj, coeffsweights->Value(jj, ii<<2)); + } + } + else + { + for (ii=1; ii<=UDegree()+1; ii++) + for (jj=1; jj<=VDegree()+1; jj++) + { + Coefs->SetValue(ii, jj, gp_Pnt(coeffsweights->Value(ii, (jj<<2)-3), + coeffsweights->Value(ii, (jj<<2)-2), + coeffsweights->Value(ii, (jj<<2)-1))); + if (rat) + WCoefs->SetValue(ii, jj, coeffsweights->Value(ii, jj<<2)); } } - else { - Coefs = coeffs; - if (rat) {WCoefs = wcoeffs;} - } -// Trim dans la base cannonique et Update des Poles et Coeffs + // Trim dans la base cannonique et Update des Poles et Coeffs -// PMN : tranfo sur les parametres + // PMN : tranfo sur les parametres Standard_Real ufirst = 2*(U1 - 0.5), - ulast = 2*(U2 - 0.5), + ulast = 2*(U2 - 0.5), vfirst = 2*(V1 - 0.5), vlast = 2*(V2 - 0.5); - if (rat) { - PLib::UTrimming (ufirst, ulast, Coefs->ChangeArray2(), - WCoefs->ChangeArray2()); - PLib::VTrimming (vfirst, vlast, Coefs->ChangeArray2(), - WCoefs->ChangeArray2()); + if (rat) + { + PLib::UTrimming(ufirst, ulast, Coefs->ChangeArray2(), + WCoefs->ChangeArray2()); + PLib::VTrimming(vfirst, vlast, Coefs->ChangeArray2(), + WCoefs->ChangeArray2()); PLib::CoefficientsPoles(Coefs->Array2(), - WCoefs->Array2(), - poles->ChangeArray2(), - weights->ChangeArray2()); + WCoefs->Array2(), + poles->ChangeArray2(), + weights->ChangeArray2()); } - else { - PLib::UTrimming (ufirst, ulast, Coefs->ChangeArray2(), - *((TColStd_Array2OfReal*) NULL)); - PLib::VTrimming (vfirst, vlast, Coefs->ChangeArray2(), - *((TColStd_Array2OfReal*) NULL)); - PLib::CoefficientsPoles (Coefs->Array2(), - *((TColStd_Array2OfReal*) NULL), - poles->ChangeArray2(), - *((TColStd_Array2OfReal*) NULL)); + else + { + PLib::UTrimming(ufirst, ulast, Coefs->ChangeArray2(), + *((TColStd_Array2OfReal*) NULL)); + PLib::VTrimming(vfirst, vlast, Coefs->ChangeArray2(), + *((TColStd_Array2OfReal*) NULL)); + PLib::CoefficientsPoles(Coefs->Array2(), + *((TColStd_Array2OfReal*) NULL), + poles->ChangeArray2(), + *((TColStd_Array2OfReal*) NULL)); } UpdateCoefficients(); } @@ -1178,45 +1209,50 @@ void Geom_BezierSurface::SetPoleRow //purpose : //======================================================================= -void Geom_BezierSurface::SetWeight (const Standard_Integer UIndex, - const Standard_Integer VIndex, - const Standard_Real Weight) +void Geom_BezierSurface::SetWeight( + const Standard_Integer UIndex, + const Standard_Integer VIndex, + const Standard_Real Weight) { - // compute new rationality + // compute new rationality Standard_Boolean wasrat = (urational||vrational); - if (!wasrat) { + if (!wasrat) + { // a weight of 1. does not turn to rational - if (Abs(Weight - 1.) <= gp::Resolution()) { + if (Abs(Weight - 1.) <= gp::Resolution()) + { UpdateCoefficients(); //Pour l'appel via SetPole return; } - + // set weights of 1. - weights = new TColStd_HArray2OfReal (1, poles->ColLength(), - 1, poles->RowLength(), 1.); - wcoeffs = new TColStd_HArray2OfReal (1, poles->ColLength(), - 1, poles->RowLength()); + weights = new TColStd_HArray2OfReal(1, poles->ColLength(), + 1, poles->RowLength(), 1.); } TColStd_Array2OfReal & Weights = weights->ChangeArray2(); - if (Weight <= gp::Resolution()) + if (Weight <= gp::Resolution()) Standard_ConstructionError::Raise("Geom_BezierSurface::SetWeight"); if (UIndex < 1 || UIndex > Weights.ColLength() || VIndex < 1 || - VIndex > Weights.RowLength()) Standard_OutOfRange::Raise(); + VIndex > Weights.RowLength()) + Standard_OutOfRange::Raise(); - if (Abs (Weight - Weights (UIndex, VIndex)) > gp::Resolution()) { + if (Abs (Weight - Weights (UIndex, VIndex)) > gp::Resolution()) + { Weights (UIndex, VIndex) = Weight; Rational(Weights, urational, vrational); } - // is it turning into non rational - if (wasrat) { - if (!(urational || vrational)) { + // is it turning into non rational + if (wasrat) + { + if (!(urational || vrational)) + { weights.Nullify(); - wcoeffs.Nullify(); + cacherational = Standard_False; } } @@ -1228,44 +1264,45 @@ void Geom_BezierSurface::SetWeight (const Standard_Integer UIndex, //purpose : //======================================================================= -void Geom_BezierSurface::SetWeightCol - (const Standard_Integer VIndex, - const TColStd_Array1OfReal& CPoleWeights) +void Geom_BezierSurface::SetWeightCol( + const Standard_Integer VIndex, + const TColStd_Array1OfReal& CPoleWeights) { Standard_Integer I; - // compute new rationality + // compute new rationality Standard_Boolean wasrat = (urational||vrational); - if (!wasrat) { + if (!wasrat) + { // set weights of 1. - weights = new TColStd_HArray2OfReal (1, poles->ColLength(), - 1, poles->RowLength(), 1.); - wcoeffs = new TColStd_HArray2OfReal (1, poles->ColLength(), - 1, poles->RowLength()); + weights = new TColStd_HArray2OfReal(1, poles->ColLength(), + 1, poles->RowLength(), 1.); } TColStd_Array2OfReal & Weights = weights->ChangeArray2(); - if (VIndex < 1 || VIndex > Weights.RowLength()) Standard_OutOfRange::Raise(); - - if (CPoleWeights.Length() != Weights.ColLength()) { + if (VIndex < 1 || VIndex > Weights.RowLength()) + Standard_OutOfRange::Raise(); + if (CPoleWeights.Length() != Weights.ColLength()) Standard_ConstructionError::Raise("Geom_BezierSurface::SetWeightCol"); - } I = CPoleWeights.Lower(); - while (I <= CPoleWeights.Upper()) { - if (CPoleWeights(I) <= gp::Resolution()) { + while (I <= CPoleWeights.Upper()) + { + if (CPoleWeights(I) <= gp::Resolution()) Standard_ConstructionError::Raise(); - } - Weights (I, VIndex) = CPoleWeights (I); + + Weights(I, VIndex) = CPoleWeights(I); I++; } - Rational(Weights, urational, vrational); + Rational(Weights, urational, vrational); - // is it turning into non rational - if (wasrat) { - if (!(urational || vrational)) { + // is it turning into non rational + if (wasrat) + { + if (!(urational || vrational)) + { weights.Nullify(); - wcoeffs.Nullify(); + cacherational = Standard_False; } } @@ -1277,19 +1314,18 @@ void Geom_BezierSurface::SetWeightCol //purpose : //======================================================================= -void Geom_BezierSurface::SetWeightRow - (const Standard_Integer UIndex, - const TColStd_Array1OfReal& CPoleWeights) +void Geom_BezierSurface::SetWeightRow( + const Standard_Integer UIndex, + const TColStd_Array1OfReal& CPoleWeights) { Standard_Integer I; - // compute new rationality + // compute new rationality Standard_Boolean wasrat = (urational||vrational); - if (!wasrat) { + if (!wasrat) + { // set weights of 1. - weights = new TColStd_HArray2OfReal (1, poles->ColLength(), - 1, poles->RowLength(), 1.); - wcoeffs = new TColStd_HArray2OfReal (1, poles->ColLength(), - 1, poles->RowLength()); + weights = new TColStd_HArray2OfReal(1, poles->ColLength(), + 1, poles->RowLength(), 1.); } TColStd_Array2OfReal & Weights = weights->ChangeArray2(); @@ -1298,26 +1334,28 @@ void Geom_BezierSurface::SetWeightRow if (CPoleWeights.Lower() < 1 || CPoleWeights.Lower() > Weights.RowLength() || CPoleWeights.Upper() < 1 || - CPoleWeights.Upper() > Weights.RowLength() ) { + CPoleWeights.Upper() > Weights.RowLength() ) Standard_ConstructionError::Raise("Geom_BezierSurface::SetWeightRow"); - } I = CPoleWeights.Lower(); - while (I <= CPoleWeights.Upper()) { - if (CPoleWeights(I) <= gp::Resolution()) { + while (I <= CPoleWeights.Upper()) + { + if (CPoleWeights(I) <= gp::Resolution()) Standard_ConstructionError::Raise(); - } - Weights (UIndex, I) = CPoleWeights (I); + + Weights(UIndex, I) = CPoleWeights(I); I++; } Rational(Weights, urational, vrational); - // is it turning into non rational - if (wasrat) { - if (!(urational || vrational)) { + // is it turning into non rational + if (wasrat) + { + if (!(urational || vrational)) + { weights.Nullify(); - wcoeffs.Nullify(); + cacherational = Standard_False; } } @@ -1449,63 +1487,64 @@ GeomAbs_Shape Geom_BezierSurface::Continuity () const //purpose : //======================================================================= -void Geom_BezierSurface::D0 (const Standard_Real U, - const Standard_Real V, - gp_Pnt& P ) const +void Geom_BezierSurface::D0( + const Standard_Real U, + const Standard_Real V, + gp_Pnt& P ) const { - - if (validcache == 1) { + if (validcache == 1) + { // // XAB : cet algorithme devient instable pour les hauts degres // RBD : Beaucoup moins maintenant avec le calcul d'un nouveau cache // sur [-1,1]. // - Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2, - uspanlenght_11 = ucachespanlenght/2, - vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2, - vspanlenght_11 = vcachespanlenght/2 ; - if (urational || vrational) { - BSplSLib::CacheD0(U, V, UDegree(), VDegree(), - uparameter_11, vparameter_11, - uspanlenght_11, vspanlenght_11, - coeffs->Array2(), - wcoeffs->Array2(), - P); - } - else { - BSplSLib::CacheD0(U, V, UDegree(), VDegree(), - uparameter_11, vparameter_11, - uspanlenght_11, vspanlenght_11, - coeffs->Array2(), - *((TColStd_Array2OfReal*) NULL), - P); - } + Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2; + Standard_Real uspanlenght_11 = ucachespanlenght/2; + Standard_Real vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2; + Standard_Real vspanlenght_11 = vcachespanlenght/2; + + BSplSLib::CacheD0(U, V, UDegree(), VDegree(), + uparameter_11, vparameter_11, + uspanlenght_11, vspanlenght_11, + coeffsweights->Array2(), + cacherational, + P); } - else { - Standard_Real array_u[2] ; - Standard_Real array_v[2] ; - Standard_Integer mult_u[2] ; - Standard_Integer mult_v[2] ; - TColStd_Array1OfReal biduknots(array_u[0],1,2); biduknots(1) = 0.; biduknots(2) = 1.; - TColStd_Array1OfInteger bidumults(mult_u[0],1,2); bidumults.Init(UDegree() + 1); - TColStd_Array1OfReal bidvknots(array_v[0],1,2); bidvknots(1) = 0.; bidvknots(2) = 1.; - TColStd_Array1OfInteger bidvmults(mult_v[0],1,2); bidvmults.Init(VDegree() + 1); - if (urational || vrational) { - BSplSLib::D0(U, V, 1,1,poles->Array2(), - weights->Array2(), - biduknots,bidvknots,bidumults,bidvmults, - UDegree(),VDegree(), - urational,vrational,Standard_False,Standard_False, - P) ; + else + { + Standard_Real array_u[2]; + Standard_Real array_v[2]; + Standard_Integer mult_u[2]; + Standard_Integer mult_v[2]; + TColStd_Array1OfReal biduknots(array_u[0], 1, 2); + biduknots(1) = 0.; + biduknots(2) = 1.; + TColStd_Array1OfInteger bidumults(mult_u[0], 1, 2); + bidumults.Init(UDegree() + 1); + TColStd_Array1OfReal bidvknots(array_v[0], 1, 2); + bidvknots(1) = 0.; + bidvknots(2) = 1.; + TColStd_Array1OfInteger bidvmults(mult_v[0], 1, 2); + bidvmults.Init(VDegree() + 1); + + if (urational || vrational) + { + BSplSLib::D0(U, V, 1, 1, poles->Array2(), + weights->Array2(), + biduknots, bidvknots, bidumults, bidvmults, + UDegree(), VDegree(), + urational, vrational, Standard_False, Standard_False, + P); } - else { - - BSplSLib::D0(U, V, 1,1,poles->Array2(), - *((TColStd_Array2OfReal*) NULL), - biduknots,bidvknots,bidumults,bidvmults, - UDegree(),VDegree(), - urational,vrational,Standard_False,Standard_False, - P) ; + else + { + BSplSLib::D0(U, V, 1, 1, poles->Array2(), + *((TColStd_Array2OfReal*) NULL), + biduknots, bidvknots, bidumults, bidvmults, + UDegree(), VDegree(), + urational, vrational, Standard_False, Standard_False, + P); } } } @@ -1515,65 +1554,65 @@ void Geom_BezierSurface::D0 (const Standard_Real U, //purpose : //======================================================================= -void Geom_BezierSurface::D1 - (const Standard_Real U, - const Standard_Real V, - gp_Pnt& P, - gp_Vec& D1U, - gp_Vec& D1V ) const +void Geom_BezierSurface::D1( + const Standard_Real U, + const Standard_Real V, + gp_Pnt& P, + gp_Vec& D1U, + gp_Vec& D1V ) const { - - if (validcache == 1) { + if (validcache == 1) + { // // XAB : cet algorithme devient instable pour les hauts degres // RBD : Beaucoup moins maintenant avec le calcul d'un nouveau cache // sur [-1,1]. // - Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2, - uspanlenght_11 = ucachespanlenght/2, - vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2, - vspanlenght_11 = vcachespanlenght/2 ; - if (urational || vrational) { - BSplSLib::CacheD1(U, V, UDegree(), VDegree(), - uparameter_11, vparameter_11, - uspanlenght_11, vspanlenght_11, - coeffs->Array2(), - wcoeffs->Array2(), - P, D1U, D1V); - } - else { - BSplSLib::CacheD1(U, V, UDegree(), VDegree(), - uparameter_11, vparameter_11, - uspanlenght_11, vspanlenght_11, - coeffs->Array2(), - *((TColStd_Array2OfReal*) NULL), - P, D1U, D1V); - } + Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2; + Standard_Real uspanlenght_11 = ucachespanlenght/2; + Standard_Real vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2; + Standard_Real vspanlenght_11 = vcachespanlenght/2; + + BSplSLib::CacheD1(U, V, UDegree(), VDegree(), + uparameter_11, vparameter_11, + uspanlenght_11, vspanlenght_11, + coeffsweights->Array2(), + cacherational, + P, D1U, D1V); } - else { - Standard_Real array_u[2] ; - Standard_Real array_v[2] ; - Standard_Integer mult_u[2] ; - Standard_Integer mult_v[2] ; - TColStd_Array1OfReal biduknots(array_u[0],1,2); biduknots(1) = 0.; biduknots(2) = 1.; - TColStd_Array1OfInteger bidumults(mult_u[0],1,2); bidumults.Init(UDegree() + 1); - TColStd_Array1OfReal bidvknots(array_v[0],1,2); bidvknots(1) = 0.; bidvknots(2) = 1.; - TColStd_Array1OfInteger bidvmults(mult_v[0],1,2); bidvmults.Init(VDegree() + 1); - if (urational || vrational) { - BSplSLib::D1(U, V, 1,1,poles->Array2(), - weights->Array2(), - biduknots,bidvknots,bidumults,bidvmults, - UDegree(),VDegree(), - urational,vrational,Standard_False,Standard_False, - P,D1U, D1V) ; + else + { + Standard_Real array_u[2]; + Standard_Real array_v[2]; + Standard_Integer mult_u[2]; + Standard_Integer mult_v[2]; + TColStd_Array1OfReal biduknots(array_u[0], 1, 2); + biduknots(1) = 0.; + biduknots(2) = 1.; + TColStd_Array1OfInteger bidumults(mult_u[0], 1, 2); + bidumults.Init(UDegree() + 1); + TColStd_Array1OfReal bidvknots(array_v[0], 1, 2); + bidvknots(1) = 0.; + bidvknots(2) = 1.; + TColStd_Array1OfInteger bidvmults(mult_v[0], 1, 2); + bidvmults.Init(VDegree() + 1); + + if (urational || vrational) + { + BSplSLib::D1(U, V, 1, 1, poles->Array2(), + weights->Array2(), + biduknots, bidvknots, bidumults, bidvmults, + UDegree(), VDegree(), + urational, vrational, Standard_False, Standard_False, + P, D1U, D1V); } else { - BSplSLib::D1(U, V, 1,1,poles->Array2(), - *((TColStd_Array2OfReal*) NULL), - biduknots,bidvknots,bidumults,bidvmults, - UDegree(),VDegree(), - urational,vrational,Standard_False,Standard_False, - P,D1U, D1V) ; + BSplSLib::D1(U, V, 1, 1, poles->Array2(), + *((TColStd_Array2OfReal*) NULL), + biduknots, bidvknots, bidumults, bidvmults, + UDegree(), VDegree(), + urational, vrational, Standard_False, Standard_False, + P, D1U, D1V); } } } @@ -1583,69 +1622,72 @@ void Geom_BezierSurface::D1 //purpose : //======================================================================= -void Geom_BezierSurface::D2 - (const Standard_Real U, - const Standard_Real V, - gp_Pnt& P, - gp_Vec& D1U, gp_Vec& D1V, - gp_Vec& D2U, gp_Vec& D2V, gp_Vec& D2UV ) const +void Geom_BezierSurface::D2( + const Standard_Real U, + const Standard_Real V, + gp_Pnt& P, + gp_Vec& D1U, + gp_Vec& D1V, + gp_Vec& D2U, + gp_Vec& D2V, + gp_Vec& D2UV ) const { - - if (validcache == 1) { + if (validcache == 1) + { // // XAB : cet algorithme devient instable pour les hauts degres // RBD : Beaucoup moins maintenant avec le calcul d'un nouveau cache // sur [-1,1]. // - Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2, - uspanlenght_11 = ucachespanlenght/2, - vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2, - vspanlenght_11 = vcachespanlenght/2 ; - if (urational || vrational) { - //-- ATTENTION a l'ORDRE d'appel ds BSPLSLIB - BSplSLib::CacheD2(U, V, UDegree(), VDegree(), - uparameter_11, vparameter_11, - uspanlenght_11, vspanlenght_11, - coeffs->Array2(), - wcoeffs->Array2(), - P, D1U, D1V, D2U, D2UV , D2V); - } - else { - //-- ATTENTION a l'ORDRE d'appel ds BSPLSLIB - BSplSLib::CacheD2(U, V, UDegree(), VDegree(), - uparameter_11, vparameter_11, - uspanlenght_11, vspanlenght_11, - coeffs->Array2(), - *((TColStd_Array2OfReal*) NULL), - P, D1U, D1V, D2U, D2UV , D2V); - } + Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2; + Standard_Real uspanlenght_11 = ucachespanlenght/2; + Standard_Real vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2; + Standard_Real vspanlenght_11 = vcachespanlenght/2 ; + + //-- ATTENTION a l'ORDRE d'appel ds BSPLSLIB + BSplSLib::CacheD2(U, V, UDegree(), VDegree(), + uparameter_11, vparameter_11, + uspanlenght_11, vspanlenght_11, + coeffsweights->Array2(), + cacherational, + P, D1U, D1V, D2U, D2UV , D2V); } - else { - Standard_Real array_u[2] ; - Standard_Real array_v[2] ; - Standard_Integer mult_u[2] ; - Standard_Integer mult_v[2] ; - TColStd_Array1OfReal biduknots(array_u[0],1,2); biduknots(1) = 0.; biduknots(2) = 1.; - TColStd_Array1OfInteger bidumults(mult_u[0],1,2); bidumults.Init(UDegree() + 1); - TColStd_Array1OfReal bidvknots(array_v[0],1,2); bidvknots(1) = 0.; bidvknots(2) = 1.; - TColStd_Array1OfInteger bidvmults(mult_v[0],1,2); bidvmults.Init(VDegree() + 1); - if (urational || vrational) { + else + { + Standard_Real array_u[2]; + Standard_Real array_v[2]; + Standard_Integer mult_u[2]; + Standard_Integer mult_v[2]; + TColStd_Array1OfReal biduknots(array_u[0], 1, 2); + biduknots(1) = 0.; + biduknots(2) = 1.; + TColStd_Array1OfInteger bidumults(mult_u[0], 1, 2); + bidumults.Init(UDegree() + 1); + TColStd_Array1OfReal bidvknots(array_v[0], 1, 2); + bidvknots(1) = 0.; + bidvknots(2) = 1.; + TColStd_Array1OfInteger bidvmults(mult_v[0], 1, 2); + bidvmults.Init(VDegree() + 1); + + if (urational || vrational) + { //-- ATTENTION a l'ORDRE d'appel ds BSPLSLIB BSplSLib::D2(U, V, 1,1,poles->Array2(), - weights->Array2(), - biduknots,bidvknots,bidumults,bidvmults, - UDegree(),VDegree(), - urational,vrational,Standard_False,Standard_False, - P,D1U, D1V, D2U, D2V , D2UV) ; + weights->Array2(), + biduknots,bidvknots,bidumults,bidvmults, + UDegree(),VDegree(), + urational,vrational,Standard_False,Standard_False, + P,D1U, D1V, D2U, D2V , D2UV) ; } - else { + else + { //-- ATTENTION a l'ORDRE d'appel ds BSPLSLIB BSplSLib::D2(U, V, 1,1,poles->Array2(), - *((TColStd_Array2OfReal*) NULL), - biduknots,bidvknots,bidumults,bidvmults, - UDegree(),VDegree(), - urational,vrational,Standard_False,Standard_False, - P,D1U, D1V, D2U, D2V, D2UV ) ; + *((TColStd_Array2OfReal*) NULL), + biduknots,bidvknots,bidumults,bidvmults, + UDegree(),VDegree(), + urational,vrational,Standard_False,Standard_False, + P,D1U, D1V, D2U, D2V, D2UV ) ; } } } @@ -2077,7 +2119,7 @@ void Geom_BezierSurface::Resolution(const Standard_Real Tolerance3D, Handle(Geom_Geometry) Geom_BezierSurface::Copy() const { Handle(Geom_BezierSurface) S = new Geom_BezierSurface - (poles, coeffs, weights, wcoeffs, urational, vrational); + (poles, weights, coeffsweights, urational, vrational); return S; } @@ -2086,9 +2128,9 @@ Handle(Geom_Geometry) Geom_BezierSurface::Copy() const //purpose : //======================================================================= -void Geom_BezierSurface::Init - (const Handle(TColgp_HArray2OfPnt)& Poles, - const Handle(TColStd_HArray2OfReal)& Weights) +void Geom_BezierSurface::Init( + const Handle(TColgp_HArray2OfPnt)& Poles, + const Handle(TColStd_HArray2OfReal)& Weights) { Standard_Integer NbUPoles = Poles->ColLength(); Standard_Integer NbVPoles = Poles->RowLength(); @@ -2097,17 +2139,16 @@ void Geom_BezierSurface::Init Standard_Integer mincls = Min(NbUPoles, NbVPoles); // set fields - poles = Poles; - coeffs = new TColgp_HArray2OfPnt (1,maxcls,1,mincls); + poles = Poles; - if (urational || vrational) { + // coeffsweights is an array of components of 4D row vector + // (three for pole point and fourth for weight) + coeffsweights = new Geom_BSplineCache(1, maxcls, 1, mincls<<2); + + if (urational || vrational) weights = Weights; - wcoeffs = new TColStd_HArray2OfReal (1,maxcls,1,mincls); - } - else { + else weights.Nullify(); - wcoeffs.Nullify(); - } UpdateCoefficients(); } @@ -2118,8 +2159,9 @@ void Geom_BezierSurface::Init //purpose : //======================================================================= -void Geom_BezierSurface::UpdateCoefficients(const Standard_Real , - const Standard_Real ) +void Geom_BezierSurface::UpdateCoefficients( + const Standard_Real , + const Standard_Real ) { maxderivinvok = Standard_False; ucacheparameter = 0.; @@ -2129,31 +2171,30 @@ void Geom_BezierSurface::UpdateCoefficients(const Standard_Real , TColStd_Array1OfReal bidvflatknots(BSplCLib::FlatBezierKnots(VDegree()), 1, 2*(VDegree()+1)); - Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2, - uspanlenght_11 = ucachespanlenght/2, - vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2, - vspanlenght_11 = vcachespanlenght/2 ; - - if ( urational || vrational ) { - BSplSLib::BuildCache(uparameter_11,vparameter_11, - uspanlenght_11,vspanlenght_11,0,0, - UDegree(),VDegree(),0,0, - biduflatknots,bidvflatknots, - poles->Array2(), - weights->Array2(), - coeffs->ChangeArray2(), - wcoeffs->ChangeArray2()); - } - else { - BSplSLib::BuildCache(uparameter_11,vparameter_11, - uspanlenght_11,vspanlenght_11,0,0, - UDegree(),VDegree(),0,0, - biduflatknots,bidvflatknots, - poles->Array2(), - *((TColStd_Array2OfReal*) NULL), - coeffs->ChangeArray2(), - *((TColStd_Array2OfReal*) NULL)); - } + Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2; + Standard_Real uspanlenght_11 = ucachespanlenght/2; + Standard_Real vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2; + Standard_Real vspanlenght_11 = vcachespanlenght/2; + + if ( urational || vrational ) + BSplSLib::BuildCache(uparameter_11, vparameter_11, + uspanlenght_11, vspanlenght_11, 0, 0, + UDegree(), VDegree(), 0, 0, + biduflatknots, bidvflatknots, + poles->Array2(), + weights->Array2(), + cacherational, + coeffsweights->ChangeArray2()); + else + BSplSLib::BuildCache(uparameter_11, vparameter_11, + uspanlenght_11, vspanlenght_11, 0, 0, + UDegree(), VDegree(), 0, 0, + biduflatknots, bidvflatknots, + poles->Array2(), + *((TColStd_Array2OfReal*) NULL), + cacherational, + coeffsweights->ChangeArray2()); + validcache = 1; } diff --git a/src/Geom/Geom_OsculatingSurface.cxx b/src/Geom/Geom_OsculatingSurface.cxx index a358db9d3b..34edc33322 100755 --- a/src/Geom/Geom_OsculatingSurface.cxx +++ b/src/Geom/Geom_OsculatingSurface.cxx @@ -439,11 +439,11 @@ Standard_Boolean Geom_OsculatingSurface::VOscSurf //======================================================================= Standard_Boolean Geom_OsculatingSurface::BuildOsculatingSurface - (const Standard_Real Param, - const Standard_Integer SUKnot, - const Standard_Integer SVKnot, - const Handle(Geom_BSplineSurface)& BS, - Handle(Geom_BSplineSurface)& BSpl) const + (const Standard_Real Param, + const Standard_Integer SUKnot, + const Standard_Integer SVKnot, + const Handle(Geom_BSplineSurface)& BS, + Handle(Geom_BSplineSurface)& BSpl) const { Standard_Integer i, j; Standard_Boolean OsculSurf=Standard_True; @@ -452,210 +452,211 @@ Standard_Boolean Geom_OsculatingSurface::BuildOsculatingSurface cout<<"======================================"<UDegree(); vdeg = BS->VDegree(); if( (IsAlongU() && vdeg <=1) || (IsAlongV() && udeg <=1)) { #if defined(DEB) && defined(OCCT_DEVELOPMENT) - cout<<" surface osculatrice nulle "<ChangeValue(i) = i-1; - PolynomialVIntervals->ChangeValue(i) = i-1; - TrueUIntervals->ChangeValue(i) = BS->UKnot(SUKnot+i-1); - TrueVIntervals->ChangeValue(i) = BS->VKnot(SVKnot+i-1); - } - + MinDegree = (Standard_Integer ) Min(udeg,vdeg) ; + MaxDegree = (Standard_Integer ) Max(udeg,vdeg) ; + + // cachepoles is a 2D grid of 4D row vectors (three components for pole and fourth for weight) + TColStd_Array2OfReal cachepoles(1, MaxDegree + 1, 1, (MinDegree + 1) * 4); + Standard_Boolean cacheRational; + // end for cache + + // for polynomial grid + Standard_Integer MaxUDegree, MaxVDegree; + Standard_Integer UContinuity, VContinuity; + + Handle(TColStd_HArray2OfInteger) NumCoeffPerSurface = + new TColStd_HArray2OfInteger(1, 1, 1, 2); + Handle(TColStd_HArray1OfReal) PolynomialUIntervals = + new TColStd_HArray1OfReal(1, 2); + Handle(TColStd_HArray1OfReal) PolynomialVIntervals = + new TColStd_HArray1OfReal(1, 2); + Handle(TColStd_HArray1OfReal) TrueUIntervals = + new TColStd_HArray1OfReal(1, 2); + Handle(TColStd_HArray1OfReal) TrueVIntervals = + new TColStd_HArray1OfReal(1, 2); + MaxUDegree = (Standard_Integer ) udeg; + MaxVDegree = (Standard_Integer ) vdeg; + + for (i=1;i<=2;i++) + { + PolynomialUIntervals->ChangeValue(i) = i-1; + PolynomialVIntervals->ChangeValue(i) = i-1; + TrueUIntervals->ChangeValue(i) = BS->UKnot(SUKnot+i-1); + TrueVIntervals->ChangeValue(i) = BS->VKnot(SVKnot+i-1); + } - Standard_Integer OscUNumCoeff=0, OscVNumCoeff=0; - if (IsAlongU()) - { + Standard_Integer OscUNumCoeff=0, OscVNumCoeff=0; + if (IsAlongU()) + { #if defined(DEB) && defined(OCCT_DEVELOPMENT) - cout<<">>>>>>>>>>> AlongU"<>>>>>>>>>> AlongU"<>>>>>>>>>> AlongV"<>>>>>>>>>> AlongV"<ChangeValue(1,1) = OscUNumCoeff; + NumCoeffPerSurface->ChangeValue(1,2) = OscVNumCoeff; + + Handle(TColStd_HArray1OfReal) Coefficients = new TColStd_HArray1OfReal(1, + NumCoeffPerSurface->Value(1,1)*NumCoeffPerSurface->Value(1,2)*3); + // end for polynomial grid + + // building the cache + Standard_Integer ULocalIndex, VLocalIndex; + Standard_Real ucacheparameter, vcacheparameter,uspanlength, vspanlength; + TColgp_Array2OfPnt NewPoles(1, BS->NbUPoles(), 1, BS->NbVPoles()); + TColStd_Array1OfReal UFlatKnots(1, BS->NbUPoles() + BS->UDegree() + 1); + TColStd_Array1OfReal VFlatKnots(1, BS->NbVPoles() + BS->VDegree() + 1); + BS->Poles(NewPoles); + BS->UKnotSequence(UFlatKnots); + BS->VKnotSequence(VFlatKnots); + + VLocalIndex = 0; + ULocalIndex = 0; + for(j = 1; j <= SVKnot; j++) VLocalIndex += BS->VMultiplicity(j); + for(i = 1; i <= SUKnot; i++) ULocalIndex += BS->UMultiplicity(i); + ucacheparameter = BS->UKnot(SUKnot); + vcacheparameter = BS->VKnot(SVKnot); + vspanlength = BS->VKnot(SVKnot + 1) - BS->VKnot(SVKnot); + uspanlength = BS->UKnot(SUKnot + 1) - BS->UKnot(SUKnot); + + // On se ramene toujours a un parametrage tel que localement ce soit l'iso + // u=0 ou v=0 qui soit degeneree + + Standard_Boolean IsVNegative = Param > vcacheparameter + vspanlength/2; + Standard_Boolean IsUNegative = Param > ucacheparameter + uspanlength/2; + + if (IsAlongU() && (Param > vcacheparameter + vspanlength/2)) + vcacheparameter = vcacheparameter + vspanlength; + if (IsAlongV() && (Param > ucacheparameter + uspanlength/2)) + ucacheparameter = ucacheparameter + uspanlength; + + BSplSLib::BuildCache(ucacheparameter, + vcacheparameter, + uspanlength, + vspanlength, + BS->IsUPeriodic(), + BS->IsVPeriodic(), + BS->UDegree(), + BS->VDegree(), + ULocalIndex, + VLocalIndex, + UFlatKnots, + VFlatKnots, + NewPoles, + BSplSLib::NoWeights(), + cacheRational, + cachepoles); + Standard_Integer m, n, index; + TColgp_Array2OfPnt OscCoeff(1,OscUNumCoeff , 1, OscVNumCoeff); + + if (IsAlongU()) + { + if (udeg > vdeg) + { + for(n = 1; n <= udeg + 1; n++) + for(m = 1; m <= vdeg; m++) + OscCoeff(n,m) = gp_Pnt(cachepoles(n, m*4), cachepoles(n, m*4+1), cachepoles(n, m*4+2)); } - NumCoeffPerSurface->ChangeValue(1,1) = OscUNumCoeff; - NumCoeffPerSurface->ChangeValue(1,2) = OscVNumCoeff; - - Handle(TColStd_HArray1OfReal) Coefficients = new TColStd_HArray1OfReal(1, - NumCoeffPerSurface->Value(1,1)*NumCoeffPerSurface->Value(1,2)*3); -// end for polynomial grid - -// building the cache - Standard_Integer ULocalIndex, VLocalIndex; - Standard_Real ucacheparameter, vcacheparameter,uspanlength, vspanlength; - TColgp_Array2OfPnt NewPoles(1, BS->NbUPoles(), 1, BS->NbVPoles()); - TColStd_Array1OfReal UFlatKnots(1, BS->NbUPoles() + BS->UDegree() + 1); - TColStd_Array1OfReal VFlatKnots(1, BS->NbVPoles() + BS->VDegree() + 1); - BS->Poles(NewPoles); - BS->UKnotSequence(UFlatKnots); - BS->VKnotSequence(VFlatKnots); - - VLocalIndex = 0; - ULocalIndex = 0; - for(j = 1; j <= SVKnot; j++) VLocalIndex += BS->VMultiplicity(j); - for(i = 1; i <= SUKnot; i++) ULocalIndex += BS->UMultiplicity(i); - ucacheparameter = BS->UKnot(SUKnot); - vcacheparameter = BS->VKnot(SVKnot); - vspanlength = BS->VKnot(SVKnot + 1) - BS->VKnot(SVKnot); - uspanlength = BS->UKnot(SUKnot + 1) - BS->UKnot(SUKnot); - - // On se ramene toujours a un parametrage tel que localement ce soit l'iso - // u=0 ou v=0 qui soit degeneree - - Standard_Boolean IsVNegative = Param > vcacheparameter + vspanlength/2; - Standard_Boolean IsUNegative = Param > ucacheparameter + uspanlength/2; - - if (IsAlongU() && (Param > vcacheparameter + vspanlength/2)) - vcacheparameter = vcacheparameter + vspanlength; - if (IsAlongV() && (Param > ucacheparameter + uspanlength/2)) - ucacheparameter = ucacheparameter + uspanlength; - - BSplSLib::BuildCache(ucacheparameter, - vcacheparameter, - uspanlength, - vspanlength, - BS->IsUPeriodic(), - BS->IsVPeriodic(), - BS->UDegree(), - BS->VDegree(), - ULocalIndex, - VLocalIndex, - UFlatKnots, - VFlatKnots, - NewPoles, - BSplSLib::NoWeights(), - cachepoles, - BSplSLib::NoWeights()); - Standard_Integer m, n, index; - TColgp_Array2OfPnt OscCoeff(1,OscUNumCoeff , 1, OscVNumCoeff); - - if (IsAlongU()) - { - if (udeg > vdeg) - { - for(n = 1; n <= udeg + 1; n++) - for(m = 1; m <= vdeg; m++) - OscCoeff(n,m) = cachepoles(n,m+1) ; - } - else - { - for(n = 1; n <= udeg + 1; n++) - for(m = 1; m <= vdeg; m++) - OscCoeff(n,m) = cachepoles(m+1,n) ; - } - if (IsVNegative) PLib::VTrimming(-1,0,OscCoeff,PLib::NoWeights2()); - - index=1; - for(n = 1; n <= udeg + 1; n++) - for(m = 1; m <= vdeg; m++) - { - Coefficients->ChangeValue(index++) = OscCoeff(n,m).X(); - Coefficients->ChangeValue(index++) = OscCoeff(n,m).Y(); - Coefficients->ChangeValue(index++) = OscCoeff(n,m).Z(); - } - } + else + { + for(n = 1; n <= udeg + 1; n++) + for(m = 1; m <= vdeg; m++) + OscCoeff(n,m) = gp_Pnt(cachepoles(m+1, (n-1)*4), cachepoles(m+1, (n-1)*4+1), cachepoles(m+1, (n-1)*4+2)); + } + if (IsVNegative) PLib::VTrimming(-1,0,OscCoeff,PLib::NoWeights2()); + + index=1; + for(n = 1; n <= udeg + 1; n++) + for(m = 1; m <= vdeg; m++) + { + Coefficients->ChangeValue(index++) = OscCoeff(n,m).X(); + Coefficients->ChangeValue(index++) = OscCoeff(n,m).Y(); + Coefficients->ChangeValue(index++) = OscCoeff(n,m).Z(); + } + } - if (IsAlongV()) - { - if (udeg > vdeg) - { - for(n = 1; n <= udeg; n++) - for(m = 1; m <= vdeg + 1; m++) - OscCoeff(n,m) = cachepoles(n+1,m); - } - else - { - for(n = 1; n <= udeg; n++) - for(m = 1; m <= vdeg + 1; m++) - OscCoeff(n,m) = cachepoles(m,n+1); - } - if (IsUNegative) PLib::UTrimming(-1,0,OscCoeff,PLib::NoWeights2()); - index=1; - for(n = 1; n <= udeg; n++) - for(m = 1; m <= vdeg + 1; m++) - { - Coefficients->ChangeValue(index++) = OscCoeff(n,m).X(); - Coefficients->ChangeValue(index++) = OscCoeff(n,m).Y(); - Coefficients->ChangeValue(index++) = OscCoeff(n,m).Z(); - } - } + if (IsAlongV()) + { + if (udeg > vdeg) + { + for(n = 1; n <= udeg; n++) + for(m = 1; m <= vdeg + 1; m++) + OscCoeff(n,m) = gp_Pnt(cachepoles(n+1, (m-1)*4), cachepoles(n+1, (m-1)*4+1), cachepoles(n+1, (m-1)*4+2)); + } + else + { + for(n = 1; n <= udeg; n++) + for(m = 1; m <= vdeg + 1; m++) + OscCoeff(n,m) = gp_Pnt(cachepoles(m, n*4), cachepoles(m, n*4+1), cachepoles(m, n*4+2)); + } + if (IsUNegative) PLib::UTrimming(-1,0,OscCoeff,PLib::NoWeights2()); + index=1; + for(n = 1; n <= udeg; n++) + for(m = 1; m <= vdeg + 1; m++) + { + Coefficients->ChangeValue(index++) = OscCoeff(n,m).X(); + Coefficients->ChangeValue(index++) = OscCoeff(n,m).Y(); + Coefficients->ChangeValue(index++) = OscCoeff(n,m).Z(); + } + } - if (IsAlongU()) MaxVDegree--; - if (IsAlongV()) MaxUDegree--; - UContinuity = - 1; - VContinuity = - 1; - - Convert_GridPolynomialToPoles Data(1,1, - UContinuity, - VContinuity, - MaxUDegree, - MaxVDegree, - NumCoeffPerSurface, - Coefficients, - PolynomialUIntervals, - PolynomialVIntervals, - TrueUIntervals, - TrueVIntervals); - -// Handle(Geom_BSplineSurface) BSpl = - BSpl =new Geom_BSplineSurface(Data.Poles()->Array2(), - Data.UKnots()->Array1(), - Data.VKnots()->Array1(), - Data.UMultiplicities()->Array1(), - Data.VMultiplicities()->Array1(), - Data.UDegree(), - Data.VDegree(), - 0, 0); + if (IsAlongU()) MaxVDegree--; + if (IsAlongV()) MaxUDegree--; + UContinuity = - 1; + VContinuity = - 1; + + Convert_GridPolynomialToPoles Data(1,1, + UContinuity, + VContinuity, + MaxUDegree, + MaxVDegree, + NumCoeffPerSurface, + Coefficients, + PolynomialUIntervals, + PolynomialVIntervals, + TrueUIntervals, + TrueVIntervals); + + // Handle(Geom_BSplineSurface) BSpl = + BSpl = new Geom_BSplineSurface(Data.Poles()->Array2(), + Data.UKnots()->Array1(), + Data.VKnots()->Array1(), + Data.UMultiplicities()->Array1(), + Data.VMultiplicities()->Array1(), + Data.UDegree(), + Data.VDegree(), + 0, 0); #if defined(DEB) && defined(OCCT_DEVELOPMENT) - cout<<"^====================================^"< +#endif +#ifndef _Standard_DefineHandle_HeaderFile +#include +#endif +#ifndef _Handle_MMgt_TShared_HeaderFile +#include +#endif + +class Standard_Transient; +class Handle(Standard_Type); +class Handle(MMgt_TShared); +class Geom_BSplineCache; + +DEFINE_STANDARD_HANDLE(Geom_BSplineCache,MMgt_TShared) + +#endif // _Handle_Geom_BSplineCache_HeaderFile diff --git a/src/NCollection/FILES b/src/NCollection/FILES index 8011e6edea..2ecd6d92a3 100755 --- a/src/NCollection/FILES +++ b/src/NCollection/FILES @@ -7,6 +7,7 @@ NCollection_IncAllocator.cxx NCollection_HeapAllocator.hxx NCollection_HeapAllocator.cxx NCollection_StdAllocator.hxx +NCollection_AlignAllocator.hxx NCollection_ListNode.hxx NCollection_BaseList.hxx diff --git a/src/NCollection/NCollection_AlignAllocator.cxx b/src/NCollection/NCollection_AlignAllocator.cxx new file mode 100644 index 0000000000..5e12fcc455 --- /dev/null +++ b/src/NCollection/NCollection_AlignAllocator.cxx @@ -0,0 +1,78 @@ +// Created on: 2013-10-11 +// Created by: Artem ZHIDKOV +// Copyright (c) 2009-2013 OPEN CASCADE SAS +// +// The content of this file is subject to the Open CASCADE Technology Public +// License Version 6.5 (the "License"). You may not use the content of this file +// except in compliance with the License. Please obtain a copy of the License +// at http://www.opencascade.org and read it completely before using this file. +// +// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its +// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France. +// +// The Original Code and all software distributed under the License is +// distributed on an "AS IS" basis, without warranty of any kind, and the +// Initial Developer hereby disclaims all such warranties, including without +// limitation, any warranties of merchantability, fitness for a particular +// purpose or non-infringement. Please see the License for the specific terms +// and conditions governing the rights and limitations under the License. + +#include +#include +#include + + +#ifdef _MSC_VER + #include +#elif (defined(__GNUC__) && __GNUC__ >= 4 && __GNUC_MINOR__ >= 1) + #include +#else + extern "C" int posix_memalign (void** thePtr, size_t theAlign, size_t theBytesCount); +#endif + + +IMPLEMENT_STANDARD_HANDLE (NCollection_AlignAllocator, NCollection_BaseAllocator) +IMPLEMENT_STANDARD_RTTIEXT(NCollection_AlignAllocator, NCollection_BaseAllocator) + +//======================================================================= +//function : Allocate +//purpose : +//======================================================================= + +void* NCollection_AlignAllocator::Allocate(const Standard_Size& theSize, + const Standard_Size& theAlignment) +{ + void* aResult; +#if defined(_MSC_VER) + aResult = _aligned_malloc(theSize, theAlignment); +#elif (defined(__GNUC__) && __GNUC__ >= 4 && __GNUC_MINOR__ >= 1) + aResult = _mm_malloc(theSize, theAlignment); +#else + if (posix_memalign(&aResult, theAlignment, theSize)) + aResult = NULL; +#endif + + if (aResult == NULL) + { + char aBuffer[96]; + Sprintf (aBuffer, "Failed to allocate %" PRIuPTR " bytes in global dynamic heap", theSize); + Standard_OutOfMemory::Raise (aBuffer); + } + return aResult; +} + +//======================================================================= +//function : Free +//purpose : +//======================================================================= + +void NCollection_AlignAllocator::Free(void* anAddress) +{ +#if defined(_MSC_VER) + _aligned_free(anAddress); +#elif (defined(__GNUC__) && __GNUC__ >= 4 && __GNUC_MINOR__ >= 1) + _mm_free(anAddress); +#else + free(anAddress); +#endif +} diff --git a/src/NCollection/NCollection_AlignAllocator.hxx b/src/NCollection/NCollection_AlignAllocator.hxx new file mode 100644 index 0000000000..2cdff7f4cd --- /dev/null +++ b/src/NCollection/NCollection_AlignAllocator.hxx @@ -0,0 +1,81 @@ +// Created on: 2013-10-11 +// Created by: Artem ZHIDKOV +// Copyright (c) 2002-2013 OPEN CASCADE SAS +// +// The content of this file is subject to the Open CASCADE Technology Public +// License Version 6.5 (the "License"). You may not use the content of this file +// except in compliance with the License. Please obtain a copy of the License +// at http://www.opencascade.org and read it completely before using this file. +// +// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its +// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France. +// +// The Original Code and all software distributed under the License is +// distributed on an "AS IS" basis, without warranty of any kind, and the +// Initial Developer hereby disclaims all such warranties, including without +// limitation, any warranties of merchantability, fitness for a particular +// purpose or non-infringement. Please see the License for the specific terms +// and conditions governing the rights and limitations under the License. + +#ifndef _NCollection_AlignAllocator_HeaderFile +#define _NCollection_AlignAllocator_HeaderFile + +#include + + +//====================================================== +// Definitions for SSE and memory alignment +//====================================================== + +// The macro for alignment of data for SSE and AVX vectorization +#if (defined(__AVX__) || (defined(_M_IX86_FP) && _M_IX86_FP>2)) +# define DATA_ALIGNMENT 32 +# define NEED_DATA_ALIGN +#elif (defined(__SSE__) || defined(__SSE2__) || defined(__SSE3__) || defined(__SSE4_1__) || defined(__SSE4_2__) || defined(_M_X64) || (defined(_M_IX86_FP) && _M_IX86_FP>0)) +# define DATA_ALIGNMENT 16 +# define NEED_DATA_ALIGN +#else +# define DATA_ALIGNMENT 4 +#endif + +// The macro for array (pointer) alignment according DATA_ALIGNMENT +#ifdef NEED_DATA_ALIGN +# ifdef _MSC_VER +# define MEMALIGN(thePtrDef) __declspec(align(DATA_ALIGNMENT)) thePtrDef +# else +# define MEMALIGN(thePtrDef) thePtrDef__attribute__((aligned(DATA_ALIGNMENT))) +# endif +#else +# define MEMALIGN(thePtrDef) thePtrDef +#endif +//====================================================== + + +class Handle_NCollection_AlignAllocator; + +/** \brief Allocates aligned memory + */ +class NCollection_AlignAllocator //: public NCollection_BaseAllocator +{ +public: + Standard_EXPORT static void* Allocate(const Standard_Size& theSize, + const Standard_Size& theAlignment = 4); + Standard_EXPORT static void Free(void * anAddress); + +protected: + //! Constructor - prohibited + NCollection_AlignAllocator(void) {}; + +private: + //! Copy constructor - prohibited + NCollection_AlignAllocator(const NCollection_AlignAllocator&); + +public: +// Declaration of CASCADE RTTI +DEFINE_STANDARD_RTTI (NCollection_AlignAllocator) +}; + +// Definition of HANDLE object using Standard_DefineHandle.hxx +DEFINE_STANDARD_HANDLE (NCollection_AlignAllocator, NCollection_BaseAllocator) + +#endif diff --git a/src/PLib/PLib.cxx b/src/PLib/PLib.cxx index 9bcd5889ac..d5be2c454e 100755 --- a/src/PLib/PLib.cxx +++ b/src/PLib/PLib.cxx @@ -665,1074 +665,300 @@ void PLib::EvalPolynomial(const Standard_Real Par, // where d is the Degree // { - Standard_Integer DegreeDimension = Degree * Dimension; + Standard_Integer kk, jj; + Standard_Real *PA = &PolynomialCoeff + Degree * Dimension; - Standard_Integer jj; - Standard_Real *RA = &Results ; - Standard_Real *PA = &PolynomialCoeff ; - Standard_Real *tmpRA = RA; - Standard_Real *tmpPA = PA + DegreeDimension; - - switch (Dimension) { - - case 1 : { - *tmpRA = *tmpPA; - if (DerivativeRequest > 0 ) { - tmpRA++ ; - Standard_Real *valRA; - Standard_Integer ii, LocalRequest; - Standard_Integer Index1, Index2; - Standard_Integer MaxIndex1, MaxIndex2; - if (DerivativeRequest < Degree) { - LocalRequest = DerivativeRequest; - MaxIndex2 = MaxIndex1 = LocalRequest; - } - else { - LocalRequest = Degree; - MaxIndex2 = MaxIndex1 = Degree; - } - MaxIndex2 --; - - for (ii = 1; ii <= LocalRequest; ii++) { - *tmpRA = 0.0e0; tmpRA ++ ; - } - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA --; - Index1 = MaxIndex1; - Index2 = MaxIndex2; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 --; - Index2 --; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); - } - } - else { - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA --; - *tmpRA = Par * (*tmpRA) + (*tmpPA); - } - } - break; - } - case 2 : { - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; - tmpPA --; - if (DerivativeRequest > 0 ) { - Standard_Real *valRA; - Standard_Integer ii, LocalRequest; - Standard_Integer Index1, Index2; - Standard_Integer MaxIndex1, MaxIndex2; - if (DerivativeRequest < Degree) { - LocalRequest = DerivativeRequest; - MaxIndex2 = MaxIndex1 = LocalRequest << 1; - } - else { - LocalRequest = Degree; - MaxIndex2 = MaxIndex1 = DegreeDimension; - } - MaxIndex2 -= 2; - - for (ii = 1; ii <= LocalRequest; ii++) { - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - } - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= 2; - - Index1 = MaxIndex1; - Index2 = MaxIndex2; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 2; - Index2 -= 2; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 1; - Index2 = MaxIndex2 + 1; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 2; - Index2 -= 2; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); - - tmpPA --; - } - } - else { - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= 2; - tmpRA = RA; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); - tmpPA --; - } - } - break; - } - case 3 : { - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; - tmpPA -= 2; - if (DerivativeRequest > 0 ) { - Standard_Real *valRA; - Standard_Integer ii, LocalRequest; - Standard_Integer Index1, Index2; - Standard_Integer MaxIndex1, MaxIndex2; - if (DerivativeRequest < Degree) { - LocalRequest = DerivativeRequest; - MaxIndex2 = MaxIndex1 = (LocalRequest << 1) + LocalRequest; - } - else { - LocalRequest = Degree; - MaxIndex2 = MaxIndex1 = DegreeDimension; - } - MaxIndex2 -= 3; - - for (ii = 1; ii <= LocalRequest; ii++) { - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - } - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= 3; - - Index1 = MaxIndex1; - Index2 = MaxIndex2; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 3; - Index2 -= 3; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 1; - Index2 = MaxIndex2 + 1; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 3; - Index2 -= 3; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 2; - Index2 = MaxIndex2 + 2; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 3; - Index2 -= 3; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); - - tmpPA -= 2; - } - } - else { - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= 3; - tmpRA = RA; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); - tmpPA -= 2; - } - } - break; - } - case 6 : { - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; - tmpPA -= 5; - if (DerivativeRequest > 0 ) { - Standard_Real *valRA; - Standard_Integer ii, LocalRequest; - Standard_Integer Index1, Index2; - Standard_Integer MaxIndex1, MaxIndex2; - if (DerivativeRequest < Degree) { - LocalRequest = DerivativeRequest; - MaxIndex2 = MaxIndex1 = (LocalRequest << 2) + (LocalRequest << 1); - } - else { - LocalRequest = Degree; - MaxIndex2 = MaxIndex1 = DegreeDimension; - } - MaxIndex2 -= 6; - - for (ii = 1; ii <= LocalRequest; ii++) { - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - } - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= 6; - - Index1 = MaxIndex1; - Index2 = MaxIndex2; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 6; - Index2 -= 6; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 1; - Index2 = MaxIndex2 + 1; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 6; - Index2 -= 6; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 2; - Index2 = MaxIndex2 + 2; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 6; - Index2 -= 6; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 3; - Index2 = MaxIndex2 + 3; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 6; - Index2 -= 6; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 4; - Index2 = MaxIndex2 + 4; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 6; - Index2 -= 6; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 5; - Index2 = MaxIndex2 + 5; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 6; - Index2 -= 6; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); - - tmpPA -= 5; - } - } - else { - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= 6; - tmpRA = RA; - - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); - tmpPA -= 5; - } - } - break; - } - case 9 : { - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; - tmpPA -= 8; - if (DerivativeRequest > 0 ) { - Standard_Real *valRA; - Standard_Integer ii, LocalRequest; - Standard_Integer Index1, Index2; - Standard_Integer MaxIndex1, MaxIndex2; - if (DerivativeRequest < Degree) { - LocalRequest = DerivativeRequest; - MaxIndex2 = MaxIndex1 = (LocalRequest << 3) + LocalRequest; - } - else { - LocalRequest = Degree; - MaxIndex2 = MaxIndex1 = DegreeDimension; - } - MaxIndex2 -= 9; - - for (ii = 1; ii <= LocalRequest; ii++) { - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - } - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= 9; - - Index1 = MaxIndex1; - Index2 = MaxIndex2; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 9; - Index2 -= 9; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 1; - Index2 = MaxIndex2 + 1; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 9; - Index2 -= 9; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 2; - Index2 = MaxIndex2 + 2; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 9; - Index2 -= 9; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 3; - Index2 = MaxIndex2 + 3; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 9; - Index2 -= 9; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 4; - Index2 = MaxIndex2 + 4; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 9; - Index2 -= 9; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 5; - Index2 = MaxIndex2 + 5; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 9; - Index2 -= 9; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 6; - Index2 = MaxIndex2 + 6; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 9; - Index2 -= 9; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 7; - Index2 = MaxIndex2 + 7; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 9; - Index2 -= 9; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 8; - Index2 = MaxIndex2 + 8; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 9; - Index2 -= 9; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); - - tmpPA -= 8; - } - } - else { - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= 9; - tmpRA = RA; - - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; + switch (DerivativeRequest) + { + case 0 : + { + Standard_Real *resD0 = &Results; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; +//#pragma vector aligned + for (kk = 0; kk < Dimension; kk++) + resD0[kk] = PA[kk]; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); - tmpPA -= 8; - } - } - break; - } - case 12 : { - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; - tmpPA -= 11; - if (DerivativeRequest > 0 ) { - Standard_Real *valRA; - Standard_Integer ii, LocalRequest; - Standard_Integer Index1, Index2; - Standard_Integer MaxIndex1, MaxIndex2; - if (DerivativeRequest < Degree) { - LocalRequest = DerivativeRequest; - MaxIndex2 = MaxIndex1 = (LocalRequest << 3) + (LocalRequest << 2); - } - else { - LocalRequest = Degree; - MaxIndex2 = MaxIndex1 = DegreeDimension; - } - MaxIndex2 -= 12; - - for (ii = 1; ii <= LocalRequest; ii++) { - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - } - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= 12; - - Index1 = MaxIndex1; - Index2 = MaxIndex2; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 12; - Index2 -= 12; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 1; - Index2 = MaxIndex2 + 1; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 12; - Index2 -= 12; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 2; - Index2 = MaxIndex2 + 2; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 12; - Index2 -= 12; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 3; - Index2 = MaxIndex2 + 3; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 12; - Index2 -= 12; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 4; - Index2 = MaxIndex2 + 4; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 12; - Index2 -= 12; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 5; - Index2 = MaxIndex2 + 5; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 12; - Index2 -= 12; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 6; - Index2 = MaxIndex2 + 6; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 12; - Index2 -= 12; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 7; - Index2 = MaxIndex2 + 7; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 12; - Index2 -= 12; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 8; - Index2 = MaxIndex2 + 8; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 12; - Index2 -= 12; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 9; - Index2 = MaxIndex2 + 9; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 12; - Index2 -= 12; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 10; - Index2 = MaxIndex2 + 10; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 12; - Index2 -= 12; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 11; - Index2 = MaxIndex2 + 11; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 12; - Index2 -= 12; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); - - tmpPA -= 11; - } + for (jj = 0; jj < Degree; jj++) + { + PA -= Dimension; +//#pragma vector aligned + for (kk = 0; kk < Dimension; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + } + break; } - else { - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= 12; - tmpRA = RA; - - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; + case 1 : + { + Standard_Real *resD0 = &Results; + Standard_Real *resD1 = &Results + Dimension; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); - tmpPA -= 11; - } + switch (Dimension) + { + case 4: + { +//#pragma vector aligned + for (kk = 0; kk < 4; kk++) + resD1[kk] = resD0[kk] = PA[kk]; + PA -= 4; + +//#pragma vector aligned + for (kk = 0; kk < 4; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + + for (jj = 1; jj < Degree; jj++) + { + PA -= 4; + +//#pragma vector aligned + for (kk = 0; kk < 4; kk++) + resD1[kk] = Par * resD1[kk] + resD0[kk]; + +//#pragma vector aligned + for (kk = 0; kk < 4; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + } + break; + } + case 8: + { +//#pragma vector aligned + for (kk = 0; kk < 8; kk++) + resD1[kk] = resD0[kk] = PA[kk]; + PA -= 8; + +//#pragma vector aligned + for (kk = 0; kk < 8; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + + for (jj = 1; jj < Degree; jj++) + { + PA -= 8; + +//#pragma vector aligned + for (kk = 0; kk < 8; kk++) + resD1[kk] = Par * resD1[kk] + resD0[kk]; + +//#pragma vector aligned + for (kk = 0; kk < 8; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + } + break; + } + case 12: + { +//#pragma vector aligned + for (kk = 0; kk < 12; kk++) + resD1[kk] = resD0[kk] = PA[kk]; + PA -= 12; + +//#pragma vector aligned + for (kk = 0; kk < 12; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + + for (jj = 1; jj < Degree; jj++) + { + PA -= 12; + +//#pragma vector aligned + for (kk = 0; kk < 12; kk++) + resD1[kk] = Par * resD1[kk] + resD0[kk]; + +//#pragma vector aligned + for (kk = 0; kk < 12; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + } + break; + } + case 16: + { +//#pragma vector aligned + for (kk = 0; kk < 16; kk++) + resD1[kk] = resD0[kk] = PA[kk]; + PA -= 16; + +//#pragma vector aligned + for (kk = 0; kk < 16; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + + for (jj = 1; jj < Degree; jj++) + { + PA -= 16; + +//#pragma vector aligned + for (kk = 0; kk < 16; kk++) + resD1[kk] = Par * resD1[kk] + resD0[kk]; + +//#pragma vector aligned + for (kk = 0; kk < 16; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + } + break; + } + case 20: + { +//#pragma vector aligned + for (kk = 0; kk < 20; kk++) + resD1[kk] = resD0[kk] = PA[kk]; + PA -= 20; + +//#pragma vector aligned + for (kk = 0; kk < 20; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + + for (jj = 1; jj < Degree; jj++) + { + PA -= 20; + +//#pragma vector aligned + for (kk = 0; kk < 20; kk++) + resD1[kk] = Par * resD1[kk] + resD0[kk]; + +//#pragma vector aligned + for (kk = 0; kk < 20; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + } + break; + } + case 24: + { +//#pragma vector aligned + for (kk = 0; kk < 24; kk++) + resD1[kk] = resD0[kk] = PA[kk]; + PA -= 24; + +//#pragma vector aligned + for (kk = 0; kk < 24; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + + for (jj = 1; jj < Degree; jj++) + { + PA -= 24; + +//#pragma vector aligned + for (kk = 0; kk < 24; kk++) + resD1[kk] = Par * resD1[kk] + resD0[kk]; + +//#pragma vector aligned + for (kk = 0; kk < 24; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + } + break; + } + default: + { +//#pragma vector aligned + for (kk = 0; kk < Dimension; kk++) + resD1[kk] = resD0[kk] = PA[kk]; + PA -= Dimension; + +//#pragma vector aligned + for (kk = 0; kk < Dimension; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + + for (jj = 1; jj < Degree; jj++) + { + PA -= Dimension; + +//#pragma vector aligned + for (kk = 0; kk < Dimension; kk++) + resD1[kk] = Par * resD1[kk] + resD0[kk]; + +//#pragma vector aligned + for (kk = 0; kk < Dimension; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + } + break; + } + } + break; } - break; - break; - } - case 15 : { - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; + case 2 : + { + Standard_Real *resD0 = &Results; + Standard_Real *resD1 = &Results + Dimension; + Standard_Real *resD2 = &Results + 2 * Dimension; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - *tmpRA = *tmpPA; tmpRA++; - tmpPA -= 14; - if (DerivativeRequest > 0 ) { - Standard_Real *valRA; - Standard_Integer ii, LocalRequest; - Standard_Integer Index1, Index2; - Standard_Integer MaxIndex1, MaxIndex2; - if (DerivativeRequest < Degree) { - LocalRequest = DerivativeRequest; - MaxIndex2 = MaxIndex1 = (LocalRequest << 4) - LocalRequest; - } - else { - LocalRequest = Degree; - MaxIndex2 = MaxIndex1 = DegreeDimension; - } - MaxIndex2 -= 15; - - for (ii = 1; ii <= LocalRequest; ii++) { - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - *tmpRA = 0.0e0; tmpRA++; - } - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= 15; - - Index1 = MaxIndex1; - Index2 = MaxIndex2; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 1; - Index2 = MaxIndex2 + 1; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 2; - Index2 = MaxIndex2 + 2; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 3; - Index2 = MaxIndex2 + 3; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 4; - Index2 = MaxIndex2 + 4; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 5; - Index2 = MaxIndex2 + 5; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 6; - Index2 = MaxIndex2 + 6; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 7; - Index2 = MaxIndex2 + 7; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 8; - Index2 = MaxIndex2 + 8; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 9; - Index2 = MaxIndex2 + 9; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 10; - Index2 = MaxIndex2 + 10; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 11; - Index2 = MaxIndex2 + 11; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 12; - Index2 = MaxIndex2 + 12; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 13; - Index2 = MaxIndex2 + 13; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - - Index1 = MaxIndex1 + 14; - Index2 = MaxIndex2 + 14; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= 15; - Index2 -= 15; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); - - tmpPA -= 14; - } + switch (Degree) + { + case 1 : + { + Standard_Real *tmpPA = PA - Dimension; + + for (kk = 0; kk < Dimension; kk++) + { + resD0[kk] = Par * PA[kk] + tmpPA[kk]; + resD1[kk] = PA[kk]; + resD2[kk] = 0.0e0; + } + break; + } + default : + { + for (kk = 0; kk < Dimension; kk++) + { + resD2[kk] = 2. * PA[kk]; + resD1[kk] = Par * PA[kk] + (Par * PA[kk] + PA[kk - Dimension]); + resD0[kk] = Par * (Par * PA[kk] + PA[kk - Dimension]) + PA[kk - 2 * Dimension]; + + } + PA -= 2 * Dimension; + + for (jj = 2; jj < Degree; jj++) + { + PA -= Dimension; + + for (kk = 0; kk < Dimension; kk++) + resD2[kk] = Par * resD2[kk] + 2. * resD1[kk]; + for (kk = 0; kk < Dimension; kk++) + resD1[kk] = Par * resD1[kk] + resD0[kk]; + for (kk = 0; kk < Dimension; kk++) + resD0[kk] = Par * resD0[kk] + PA[kk]; + } + break; + } + } + break; } - else { - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= 15; - tmpRA = RA; - - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; + default : + { + Standard_Integer ii, LocalRequest, LocReqDim; + Standard_Real *RA = &Results; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; + for (kk = 0; kk < Dimension; kk++) + RA[kk] = PA[kk]; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; + if (DerivativeRequest < Degree) + LocalRequest = DerivativeRequest; + else + LocalRequest = Degree; + LocReqDim = LocalRequest * Dimension; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - *tmpRA = Par * (*tmpRA) + (*tmpPA); - tmpPA -= 14; - } - } - break; - } - default : { - Standard_Integer kk ; - for ( kk = 0 ; kk < Dimension ; kk++) { - *tmpRA = *tmpPA; tmpRA++; tmpPA++; - } - tmpPA -= Dimension; - if (DerivativeRequest > 0 ) { - Standard_Real *valRA; - Standard_Integer ii, LocalRequest; - Standard_Integer Index1, Index2; - Standard_Integer MaxIndex1, MaxIndex2; - if (DerivativeRequest < Degree) { - LocalRequest = DerivativeRequest; - MaxIndex2 = MaxIndex1 = LocalRequest * Dimension; - } - else { - LocalRequest = Degree; - MaxIndex2 = MaxIndex1 = DegreeDimension; - } - MaxIndex2 -= Dimension; - - for (ii = 1; ii <= MaxIndex1; ii++) { - *tmpRA = 0.0e0; tmpRA++; - } - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= Dimension; - - for (kk = 0 ; kk < Dimension ; kk++) { - Index1 = MaxIndex1 + kk; - Index2 = MaxIndex2 + kk; - - for (ii = LocalRequest ; ii > 0 ; ii--) { - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ; - Index1 -= Dimension; - Index2 -= Dimension; - } - valRA = &RA[Index1]; - *valRA = Par * (*valRA) + (*tmpPA); tmpPA++; - } - tmpPA -= Dimension; - } - } - else { - - for (jj = Degree ; jj > 0 ; jj--) { - tmpPA -= Dimension; - tmpRA = RA; - - for (kk = 0 ; kk < Dimension ; kk++) { - *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++; - } - tmpPA -= Dimension; - } - } + for (ii = 1; ii <= DerivativeRequest; ii++) + { + RA += Dimension; + for (kk = 0; kk < Dimension; kk++) + RA[kk] = 0.0e0; + } + + for (jj = 0; jj < Degree; jj++) + { + RA = &Results + LocReqDim; + for (ii = LocalRequest; ii > 1; ii--) + { + for (kk = 0; kk < Dimension; kk++) + RA[kk] = Par * RA[kk] + ((double)ii) * RA[kk - Dimension]; + RA -= Dimension; + } + + for (kk = 0; kk < Dimension; kk++) + RA[kk] = Par * RA[kk] + RA[kk - Dimension]; + RA -= Dimension; + PA -= Dimension; + + for (kk = 0; kk < Dimension; kk++) + RA[kk] = Par * RA[kk] + PA[kk]; + } + break; } } } -- 2.39.5