From c63628e84569cd2d89da418017d662e8cc39ba48 Mon Sep 17 00:00:00 2001 From: nbv Date: Fri, 4 Oct 2013 18:26:23 +0400 Subject: [PATCH] 0024211: Definition of Basic Runtime Check parameter causes regression in debug mode Out of ChoixRef array boundaries. Uninitialized variable in IntCurve_IntPolyPolyGen::findIntersect(...) function. Handling of infinity numbers in sprops command is added. test (CPU-limit) --- src/GProp/GProp_SGProps.gxx | 1060 ++++++++++++++++++---------- src/IntImp/IntImp_Int2S.gxx | 108 +-- src/IntWalk/IntWalk_PWalking_1.gxx | 4 +- tests/bugs/modalg_1/bug19793_2 | 2 +- 4 files changed, 741 insertions(+), 433 deletions(-) diff --git a/src/GProp/GProp_SGProps.gxx b/src/GProp/GProp_SGProps.gxx index abf1b6a23e..099a6c0edc 100755 --- a/src/GProp/GProp_SGProps.gxx +++ b/src/GProp/GProp_SGProps.gxx @@ -102,10 +102,97 @@ static HMath_Vector IxyU = new math_Vector(1,SM,0.0); static HMath_Vector IxzU = new math_Vector(1,SM,0.0); static HMath_Vector IyzU = new math_Vector(1,SM,0.0); +static inline Standard_Real MultiplicationInf(Standard_Real theMA, Standard_Real theMB) +{ + if((theMA == 0.0) || (theMB == 0.0)) //strictly zerro (without any tolerances) + return 0.0; + + if(Precision::IsPositiveInfinite(theMA)) + { + if(theMB < 0.0) + return -Precision::Infinite(); + else + return Precision::Infinite(); + } + + if(Precision::IsPositiveInfinite(theMB)) + { + if(theMA < 0.0) + return -Precision::Infinite(); + else + return Precision::Infinite(); + } + + if(Precision::IsNegativeInfinite(theMA)) + { + if(theMB < 0.0) + return +Precision::Infinite(); + else + return -Precision::Infinite(); + } + + if(Precision::IsNegativeInfinite(theMB)) + { + if(theMA < 0.0) + return +Precision::Infinite(); + else + return -Precision::Infinite(); + } + + return (theMA * theMB); +} + +static inline Standard_Real AdditionInf(Standard_Real theMA, Standard_Real theMB) +{ + if(Precision::IsPositiveInfinite(theMA)) + { + if(Precision::IsNegativeInfinite(theMB)) + return 0.0; + else + return Precision::Infinite(); + } + + if(Precision::IsPositiveInfinite(theMB)) + { + if(Precision::IsNegativeInfinite(theMA)) + return 0.0; + else + return Precision::Infinite(); + } + + if(Precision::IsNegativeInfinite(theMA)) + { + if(Precision::IsPositiveInfinite(theMB)) + return 0.0; + else + return -Precision::Infinite(); + } + + if(Precision::IsNegativeInfinite(theMB)) + { + if(Precision::IsPositiveInfinite(theMA)) + return 0.0; + else + return -Precision::Infinite(); + } + + return (theMA + theMB); +} + +static inline Standard_Real Multiplication(Standard_Real theMA, Standard_Real theMB) +{ + return (theMA * theMB); +} + +static inline Standard_Real Addition(Standard_Real theMA, Standard_Real theMB) +{ + return (theMA + theMB); +} + static Standard_Integer FillIntervalBounds(Standard_Real A, Standard_Real B, const TColStd_Array1OfReal& Knots, - HMath_Vector& VA, + HMath_Vector& VA, HMath_Vector& VB) { Standard_Integer i = 1, iEnd = Knots.Upper(), j = 1, k = 1; @@ -129,14 +216,14 @@ static Standard_Integer FillIntervalBounds(Standard_Real A, } static inline Standard_Integer MaxSubs(Standard_Integer n, Standard_Integer coeff = SUBS_POWER){ -// return n = IntegerLast()/coeff < n? IntegerLast(): n*coeff + 1; + // return n = IntegerLast()/coeff < n? IntegerLast(): n*coeff + 1; return Min((n * coeff + 1),SM); } static Standard_Integer LFillIntervalBounds(Standard_Real A, Standard_Real B, const TColStd_Array1OfReal& Knots, - const Standard_Integer NumSubs) + const Standard_Integer NumSubs) { Standard_Integer iEnd = Knots.Upper(), jEnd = L1->Upper(); if(iEnd - 1 > jEnd){ @@ -162,7 +249,7 @@ static Standard_Integer LFillIntervalBounds(Standard_Real A, static Standard_Integer UFillIntervalBounds(Standard_Real A, Standard_Real B, const TColStd_Array1OfReal& Knots, - const Standard_Integer NumSubs) + const Standard_Integer NumSubs) { Standard_Integer iEnd = Knots.Upper(), jEnd = U1->Upper(); if(iEnd - 1 > jEnd){ @@ -190,10 +277,16 @@ static Standard_Real CCompute(Face& S, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, - const Standard_Real EpsDim, - const Standard_Boolean isErrorCalculation, + const Standard_Real EpsDim, + const Standard_Boolean isErrorCalculation, const Standard_Boolean isVerifyComputation) { + Standard_Real (*FuncAdd)(Standard_Real, Standard_Real); + Standard_Real (*FuncMul)(Standard_Real, Standard_Real); + + FuncAdd = Addition; + FuncMul = Multiplication; + Standard_Boolean isNaturalRestriction = S.NaturalRestriction(); Standard_Integer NumSubs = SUBS_POWER; @@ -206,7 +299,17 @@ static Standard_Real CCompute(Face& S, //Face parametrization in U and V direction Standard_Real BV1, BV2, v; Standard_Real BU1, BU2, u1, u2, um, ur, u; - S.Bounds (BU1, BU2, BV1, BV2); u1 = BU1; + S.Bounds (BU1, BU2, BV1, BV2); + u1 = BU1; + + if(Precision::IsInfinite(BU1) || Precision::IsInfinite(BU2) || + Precision::IsInfinite(BV1) || Precision::IsInfinite(BV2)) + { + FuncAdd = AdditionInf; + FuncMul = MultiplicationInf; + } + + //location point used to compute the inertia Standard_Real xloc, yloc, zloc; loc.Coord (xloc, yloc, zloc); // use member of parent class @@ -221,7 +324,7 @@ static Standard_Real CCompute(Face& S, Standard_Real Dul; // Dul = Du / Dl Standard_Real CDim[2], CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz; Standard_Real LocDim[2], LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy, LocIxz, LocIyz; - + Standard_Real ErrorU, ErrorL, ErrorLMax = 0.0, Eps=0.0, EpsL=0.0, EpsU=0.0; Standard_Integer iD = 0, NbLSubs, iLS, iLSubEnd, iGL, iGLEnd, NbLGaussP[2], LRange[2], iL, kL, kLEnd, IL, JL; @@ -237,191 +340,327 @@ static Standard_Real CCompute(Face& S, NbUGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbUGaussP[0])); math::GaussPoints(NbUGaussP[0],UGaussP0); math::GaussWeights(NbUGaussP[0],UGaussW0); math::GaussPoints(NbUGaussP[1],UGaussP1); math::GaussWeights(NbUGaussP[1],UGaussW1); - + NbUSubs = S.SUIntSubs(); TColStd_Array1OfReal UKnots(1,NbUSubs+1); S.UKnots(UKnots); - - while (isNaturalRestriction || D.More()) { - if(isNaturalRestriction){ + while (isNaturalRestriction || D.More()) + { + if(isNaturalRestriction) + { NbLGaussP[0] = Min(2*NbUGaussP[0],math::GaussPointsMax()); - }else{ + } + else + { S.Load(D.Value()); ++iD; NbLGaussP[0] = S.LIntOrder(EpsDim); } - NbLGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbLGaussP[0])); math::GaussPoints(NbLGaussP[0],LGaussP0); math::GaussWeights(NbLGaussP[0],LGaussW0); math::GaussPoints(NbLGaussP[1],LGaussP1); math::GaussWeights(NbLGaussP[1],LGaussW1); - + NbLSubs = isNaturalRestriction? S.SVIntSubs(): S.LIntSubs(); TColStd_Array1OfReal LKnots(1,NbLSubs+1); - if(isNaturalRestriction){ + if(isNaturalRestriction) + { S.VKnots(LKnots); l1 = BV1; l2 = BV2; - }else{ + } + else + { S.LKnots(LKnots); l1 = S.FirstParameter(); l2 = S.LastParameter(); } + ErrorL = 0.0; kLEnd = 1; JL = 0; //OCC503(apo): if(Abs(l2-l1) < EPS_PARAM) continue; - if(Abs(l2-l1) > EPS_PARAM) { + if(Abs(l2-l1) > EPS_PARAM) + { iLSubEnd = LFillIntervalBounds(l1, l2, LKnots, NumSubs); LMaxSubs = MaxSubs(iLSubEnd); - if(LMaxSubs > DimL.Vector()->Upper()) LMaxSubs = DimL.Vector()->Upper(); - DimL.Init(0.0,1,LMaxSubs); ErrL.Init(0.0,1,LMaxSubs); ErrUL.Init(0.0,1,LMaxSubs); - do{// while: L - if(++JL > iLSubEnd){ - LRange[0] = IL = ErrL->Max(); LRange[1] = JL; - L1(JL) = (L1(IL) + L2(IL))/2.0; L2(JL) = L2(IL); L2(IL) = L1(JL); - }else LRange[0] = IL = JL; - if(JL == LMaxSubs || Abs(L2(JL) - L1(JL)) < EPS_PARAM) - if(kLEnd == 1){ - DimL(JL) = ErrL(JL) = IxL(JL) = IyL(JL) = IzL(JL) = - IxxL(JL) = IyyL(JL) = IzzL(JL) = IxyL(JL) = IxzL(JL) = IyzL(JL) = 0.0; - }else{ - JL--; - EpsL = ErrorL; Eps = EpsL/0.9; - break; - } - else - for(kL=0; kL < kLEnd; kL++){ - iLS = LRange[kL]; - lm = 0.5*(L2(iLS) + L1(iLS)); - lr = 0.5*(L2(iLS) - L1(iLS)); - CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0; - for(iGL=0; iGL < iGLEnd; iGL++){// - CDim[iGL] = 0.0; - for(iL=1; iL<=NbLGaussP[iGL]; iL++){ - l = lm + lr*(*LGaussP[iGL])(iL); - if(isNaturalRestriction){ - v = l; u2 = BU2; Dul = (*LGaussW[iGL])(iL); - }else{ - S.D12d (l, Puv, Vuv); - Dul = Vuv.Y()*(*LGaussW[iGL])(iL); // Dul = Du / Dl - if(Abs(Dul) < EPS_PARAM) continue; - v = Puv.Y(); u2 = Puv.X(); - //Check on cause out off bounds of value current parameter - if(v < BV1) v = BV1; else if(v > BV2) v = BV2; - if(u2 < BU1) u2 = BU1; else if(u2 > BU2) u2 = BU2; - } - ErrUL(iLS) = 0.0; - kUEnd = 1; JU = 0; - if(Abs(u2-u1) < EPS_PARAM) continue; - iUSubEnd = UFillIntervalBounds(u1, u2, UKnots, NumSubs); - UMaxSubs = MaxSubs(iUSubEnd); - if(UMaxSubs > DimU.Vector()->Upper()) UMaxSubs = DimU.Vector()->Upper(); - DimU.Init(0.0,1,UMaxSubs); ErrU.Init(0.0,1,UMaxSubs); ErrorU = 0.0; - do{//while: U - if(++JU > iUSubEnd){ - URange[0] = IU = ErrU->Max(); URange[1] = JU; - U1(JU) = (U1(IU)+U2(IU))/2.0; U2(JU) = U2(IU); U2(IU) = U1(JU); - }else URange[0] = IU = JU; - if(JU == UMaxSubs || Abs(U2(JU) - U1(JU)) < EPS_PARAM) - if(kUEnd == 1){ - DimU(JU) = ErrU(JU) = IxU(JU) = IyU(JU) = IzU(JU) = - IxxU(JU) = IyyU(JU) = IzzU(JU) = IxyU(JU) = IxzU(JU) = IyzU(JU) = 0.0; - }else{ - JU--; - EpsU = ErrorU; Eps = EpsU*Abs((u2-u1)*Dul)/0.1; EpsL = 0.9*Eps; - break; - } - else - for(kU=0; kU < kUEnd; kU++){ - iUS = URange[kU]; - um = 0.5*(U2(iUS) + U1(iUS)); - ur = 0.5*(U2(iUS) - U1(iUS)); - LocIx = LocIy = LocIz = LocIxx = LocIyy = LocIzz = LocIxy = LocIxz = LocIyz = 0.0; - iGUEnd = iGLEnd - iGL; - for(iGU=0; iGU < iGUEnd; iGU++){// - LocDim[iGU] = 0.0; - for(iU=1; iU<=NbUGaussP[iGU]; iU++){ - u = um + ur*(*UGaussP[iGU])(iU); - S.Normal(u, v, Ps, VNor); - ds = VNor.Magnitude(); //Jacobien(x,y,z) -> (u,v)=||n|| - ds *= (*UGaussW[iGU])(iU); - LocDim[iGU] += ds; - if(iGU > 0) continue; - Ps.Coord(x, y, z); - x -= xloc; y -= yloc; z -= zloc; - LocIx += x*ds; LocIy += y*ds; LocIz += z*ds; - LocIxy += x*y*ds; LocIyz += y*z*ds; LocIxz += x*z*ds; - x *= x; y *= y; z *= z; - LocIxx += (y+z)*ds; LocIyy += (x+z)*ds; LocIzz += (x+y)*ds; - }//for: iU - }//for: iGU - DimU(iUS) = LocDim[0]*ur; - if(iGL > 0) continue; - ErrU(iUS) = Abs(LocDim[1]-LocDim[0])*ur; - IxU(iUS) = LocIx*ur; IyU(iUS) = LocIy*ur; IzU(iUS) = LocIz*ur; - IxxU(iUS) = LocIxx*ur; IyyU(iUS) = LocIyy*ur; IzzU(iUS) = LocIzz*ur; - IxyU(iUS) = LocIxy*ur; IxzU(iUS) = LocIxz*ur; IyzU(iUS) = LocIyz*ur; - }//for: kU (iUS) - if(JU == iUSubEnd) kUEnd = 2; - if(kUEnd == 2) ErrorU = ErrU(ErrU->Max()); - }while((ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1); - for(i=1; i<=JU; i++) CDim[iGL] += DimU(i)*Dul; - if(iGL > 0) continue; - ErrUL(iLS) = ErrorU*Abs((u2-u1)*Dul); - for(i=1; i<=JU; i++){ - CIx += IxU(i)*Dul; CIy += IyU(i)*Dul; CIz += IzU(i)*Dul; - CIxx += IxxU(i)*Dul; CIyy += IyyU(i)*Dul; CIzz += IzzU(i)*Dul; - CIxy += IxyU(i)*Dul; CIxz += IxzU(i)*Dul; CIyz += IyzU(i)*Dul; - } - }//for: iL - }//for: iGL - DimL(iLS) = CDim[0]*lr; - if(iGLEnd == 2) ErrL(iLS) = Abs(CDim[1]-CDim[0])*lr + ErrUL(iLS); - IxL(iLS) = CIx*lr; IyL(iLS) = CIy*lr; IzL(iLS) = CIz*lr; - IxxL(iLS) = CIxx*lr; IyyL(iLS) = CIyy*lr; IzzL(iLS) = CIzz*lr; - IxyL(iLS) = CIxy*lr; IxzL(iLS) = CIxz*lr; IyzL(iLS) = CIyz*lr; - }//for: (kL)iLS - // Calculate/correct epsilon of computation by current value of Dim - //That is need for not spend time for - if(JL == iLSubEnd){ - kLEnd = 2; - Standard_Real DDim = 0.0; - for(i=1; i<=JL; i++) DDim += DimL(i); - DDim = Abs(DDim*EpsDim); - if(DDim > Eps) { - Eps = DDim; EpsL = 0.9*Eps; - } - } - if(kLEnd == 2) ErrorL = ErrL(ErrL->Max()); - }while((ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1); - for(i=1; i<=JL; i++){ - Dim += DimL(i); - Ix += IxL(i); Iy += IyL(i); Iz += IzL(i); - Ixx += IxxL(i); Iyy += IyyL(i); Izz += IzzL(i); - Ixy += IxyL(i); Ixz += IxzL(i); Iyz += IyzL(i); - } - ErrorLMax = Max(ErrorLMax, ErrorL); - } - if(isNaturalRestriction) break; - D.Next(); - } - if(Abs(Dim) >= EPS_DIM){ - Ix /= Dim; Iy /= Dim; Iz /= Dim; - g.SetCoord (Ix, Iy, Iz); - }else{ - Dim =0.0; - g.SetCoord (0., 0.,0.); - } - inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz), - gp_XYZ (-Ixy, Iyy, -Iyz), - gp_XYZ (-Ixz, -Iyz, Izz)); + if(LMaxSubs > DimL.Vector()->Upper()) + LMaxSubs = DimL.Vector()->Upper(); - if(iGLEnd == 2) Eps = Dim != 0.0? ErrorLMax/Abs(Dim): 0.0; - else Eps = EpsDim; - return Eps; + DimL.Init(0.0,1,LMaxSubs); ErrL.Init(0.0,1,LMaxSubs); ErrUL.Init(0.0,1,LMaxSubs); + + do{// while: L + if(++JL > iLSubEnd) + { + LRange[0] = IL = ErrL->Max(); + LRange[1] = JL; + L1(JL) = (L1(IL) + L2(IL))/2.0; + L2(JL) = L2(IL); + L2(IL) = L1(JL); + } + else + LRange[0] = IL = JL; + + if(JL == LMaxSubs || Abs(L2(JL) - L1(JL)) < EPS_PARAM) + if(kLEnd == 1) + { + DimL(JL) = ErrL(JL) = IxL(JL) = IyL(JL) = IzL(JL) = + IxxL(JL) = IyyL(JL) = IzzL(JL) = IxyL(JL) = IxzL(JL) = IyzL(JL) = 0.0; + }else{ + JL--; + EpsL = ErrorL; Eps = EpsL/0.9; + break; + } + else + for(kL=0; kL < kLEnd; kL++) + { + iLS = LRange[kL]; + lm = 0.5*(L2(iLS) + L1(iLS)); + lr = 0.5*(L2(iLS) - L1(iLS)); + CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0; + + for(iGL=0; iGL < iGLEnd; iGL++) + { + CDim[iGL] = 0.0; + for(iL=1; iL<=NbLGaussP[iGL]; iL++) + { + l = lm + lr*(*LGaussP[iGL])(iL); + if(isNaturalRestriction) + { + v = l; u2 = BU2; Dul = (*LGaussW[iGL])(iL); + } + else + { + S.D12d (l, Puv, Vuv); + Dul = Vuv.Y()*(*LGaussW[iGL])(iL); // Dul = Du / Dl + if(Abs(Dul) < EPS_PARAM) + continue; + + v = Puv.Y(); + u2 = Puv.X(); + //Check on cause out off bounds of value current parameter + if(v < BV1) + v = BV1; + else if(v > BV2) + v = BV2; + + if(u2 < BU1) + u2 = BU1; + else if(u2 > BU2) + u2 = BU2; + } + + ErrUL(iLS) = 0.0; + kUEnd = 1; + JU = 0; + + if(Abs(u2-u1) < EPS_PARAM) + continue; + + iUSubEnd = UFillIntervalBounds(u1, u2, UKnots, NumSubs); + UMaxSubs = MaxSubs(iUSubEnd); + if(UMaxSubs > DimU.Vector()->Upper()) + UMaxSubs = DimU.Vector()->Upper(); + + DimU.Init(0.0,1,UMaxSubs); ErrU.Init(0.0,1,UMaxSubs); ErrorU = 0.0; + + do{//while: U + if(++JU > iUSubEnd) + { + URange[0] = IU = ErrU->Max(); + URange[1] = JU; + U1(JU) = (U1(IU)+U2(IU))/2.0; + U2(JU) = U2(IU); + U2(IU) = U1(JU); + } + else + URange[0] = IU = JU; + + if(JU == UMaxSubs || Abs(U2(JU) - U1(JU)) < EPS_PARAM) + if(kUEnd == 1) + { + DimU(JU) = ErrU(JU) = IxU(JU) = IyU(JU) = IzU(JU) = + IxxU(JU) = IyyU(JU) = IzzU(JU) = IxyU(JU) = IxzU(JU) = IyzU(JU) = 0.0; + }else + { + JU--; + EpsU = ErrorU; Eps = EpsU*Abs((u2-u1)*Dul)/0.1; EpsL = 0.9*Eps; + break; + } + else + for(kU=0; kU < kUEnd; kU++) + { + iUS = URange[kU]; + um = 0.5*(U2(iUS) + U1(iUS)); + ur = 0.5*(U2(iUS) - U1(iUS)); + LocIx = LocIy = LocIz = LocIxx = LocIyy = LocIzz = LocIxy = LocIxz = LocIyz = 0.0; + iGUEnd = iGLEnd - iGL; + for(iGU=0; iGU < iGUEnd; iGU++) + { + LocDim[iGU] = 0.0; + for(iU=1; iU <= NbUGaussP[iGU]; iU++) + { + u = um + ur*(*UGaussP[iGU])(iU); + S.Normal(u, v, Ps, VNor); + ds = VNor.Magnitude(); //Jacobien(x,y,z) -> (u,v)=||n|| + ds *= (*UGaussW[iGU])(iU); + LocDim[iGU] += ds; + + if(iGU > 0) + continue; + + Ps.Coord(x, y, z); + x = FuncAdd(x, -xloc); + y = FuncAdd(y, -yloc); + z = FuncAdd(z, -zloc); + + const Standard_Real XdS = FuncMul(x, ds); + const Standard_Real YdS = FuncMul(y, ds); + const Standard_Real ZdS = FuncMul(z, ds); + + LocIx = FuncAdd(LocIx, XdS); + LocIy = FuncAdd(LocIy, YdS); + LocIz = FuncAdd(LocIz, ZdS); + LocIxy = FuncAdd(LocIxy, FuncMul(x, YdS)); + LocIyz = FuncAdd(LocIyz, FuncMul(y, ZdS)); + LocIxz = FuncAdd(LocIxz, FuncMul(x, ZdS)); + x = Precision::IsInfinite(x) ? Precision::Infinite() : x*x; + y = Precision::IsInfinite(y) ? Precision::Infinite() : y*y; + z = Precision::IsInfinite(z) ? Precision::Infinite() : z*z; + LocIxx = FuncAdd(LocIxx, FuncAdd(YdS, ZdS)); + LocIyy = FuncAdd(LocIyy, FuncAdd(XdS, ZdS)); + LocIzz = FuncAdd(LocIzz, FuncAdd(XdS, YdS)); + }//for: iU + }//for: iGU + + DimU(iUS) = FuncMul(LocDim[0],ur); + if(iGL > 0) + continue; + + ErrU(iUS) = FuncMul(Abs(LocDim[1]-LocDim[0]), ur); + IxU(iUS) = FuncMul(LocIx, ur); + IyU(iUS) = FuncMul(LocIy, ur); + IzU(iUS) = FuncMul(LocIz, ur); + IxxU(iUS) = FuncMul(LocIxx, ur); + IyyU(iUS) = FuncMul(LocIyy, ur); + IzzU(iUS) = FuncMul(LocIzz, ur); + IxyU(iUS) = FuncMul(LocIxy, ur); + IxzU(iUS) = FuncMul(LocIxz, ur); + IyzU(iUS) = FuncMul(LocIyz, ur); + }//for: kU (iUS) + + if(JU == iUSubEnd) + kUEnd = 2; + + if(kUEnd == 2) + ErrorU = ErrU(ErrU->Max()); + }while((ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1); + + for(i=1; i<=JU; i++) + CDim[iGL] = FuncAdd(CDim[iGL], FuncMul(DimU(i), Dul)); + + if(iGL > 0) + continue; + + ErrUL(iLS) = ErrorU*Abs((u2-u1)*Dul); + for(i=1; i<=JU; i++) + { + CIx = FuncAdd(CIx, FuncMul(IxU(i), Dul)); + CIy = FuncAdd(CIy, FuncMul(IyU(i), Dul)); + CIz = FuncAdd(CIz, FuncMul(IzU(i), Dul)); + CIxx = FuncAdd(CIxx, FuncMul(IxxU(i), Dul)); + CIyy = FuncAdd(CIyy, FuncMul(IyyU(i), Dul)); + CIzz = FuncAdd(CIzz, FuncMul(IzzU(i), Dul)); + CIxy = FuncAdd(CIxy, FuncMul(IxyU(i), Dul)); + CIxz = FuncAdd(CIxz, FuncMul(IxzU(i), Dul)); + CIyz = FuncAdd(CIyz, FuncMul(IyzU(i), Dul)); + } + }//for: iL + }//for: iGL + + DimL(iLS) = FuncMul(CDim[0], lr); + if(iGLEnd == 2) + ErrL(iLS) = FuncAdd(FuncMul(Abs(CDim[1]-CDim[0]),lr), ErrUL(iLS)); + + IxL(iLS) = FuncMul(CIx, lr); + IyL(iLS) = FuncMul(CIy, lr); + IzL(iLS) = FuncMul(CIz, lr); + IxxL(iLS) = FuncMul(CIxx, lr); + IyyL(iLS) = FuncMul(CIyy, lr); + IzzL(iLS) = FuncMul(CIzz, lr); + IxyL(iLS) = FuncMul(CIxy, lr); + IxzL(iLS) = FuncMul(CIxz, lr); + IyzL(iLS) = FuncMul(CIyz, lr); + }//for: (kL)iLS + // Calculate/correct epsilon of computation by current value of Dim + //That is need for not spend time for + if(JL == iLSubEnd) + { + kLEnd = 2; + Standard_Real DDim = 0.0; + for(i=1; i<=JL; i++) + DDim += DimL(i); + + DDim = Abs(DDim*EpsDim); + if(DDim > Eps) + { + Eps = DDim; + EpsL = 0.9*Eps; + } + } + + if(kLEnd == 2) + ErrorL = ErrL(ErrL->Max()); + }while((ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1); + + for(i=1; i<=JL; i++) + { + Dim = FuncAdd(Dim, DimL(i)); + Ix = FuncAdd(Ix, IxL(i)); + Iy = FuncAdd(Iy, IyL(i)); + Iz = FuncAdd(Iz, IzL(i)); + Ixx = FuncAdd(Ixx, IxxL(i)); + Iyy = FuncAdd(Iyy, IyyL(i)); + Izz = FuncAdd(Izz, IzzL(i)); + Ixy = FuncAdd(Ixy, IxyL(i)); + Ixz = FuncAdd(Ixz, IxzL(i)); + Iyz = FuncAdd(Iyz, IyzL(i)); + } + + ErrorLMax = Max(ErrorLMax, ErrorL); + } + + if(isNaturalRestriction) + break; + + D.Next(); + } + + if(Abs(Dim) >= EPS_DIM) + { + Ix /= Dim; + Iy /= Dim; + Iz /= Dim; + g.SetCoord (Ix, Iy, Iz); + } + else + { + Dim =0.0; + g.SetCoord (0., 0.,0.); + } + + inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz), + gp_XYZ (-Ixy, Iyy, -Iyz), + gp_XYZ (-Ixz, -Iyz, Izz)); + + if(iGLEnd == 2) + Eps = Dim != 0.0? ErrorLMax/Abs(Dim): 0.0; + else + Eps = EpsDim; + + return Eps; } static Standard_Real Compute(Face& S, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, - Standard_Real EpsDim) + Standard_Real EpsDim) { Standard_Boolean isErrorCalculation = 0.0 > EpsDim || EpsDim < 0.001? 1: 0; Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0; @@ -431,7 +670,7 @@ static Standard_Real Compute(Face& S, const gp_Pnt& loc, Standard_Real& Dim, gp_ } static Standard_Real Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, - Standard_Real EpsDim) + Standard_Real EpsDim) { Standard_Boolean isErrorCalculation = 0.0 > EpsDim || EpsDim < 0.001? 1: 0; Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0; @@ -439,260 +678,315 @@ static Standard_Real Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Rea return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation); } -static void Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& dim, gp_Pnt& g, gp_Mat& inertia){ - Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz; - dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0; +static void Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& dim, gp_Pnt& g, gp_Mat& inertia) +{ + Standard_Real (*FuncAdd)(Standard_Real, Standard_Real); + Standard_Real (*FuncMul)(Standard_Real, Standard_Real); - Standard_Real x, y, z; - Standard_Integer NbCGaussgp_Pnts = 0; + FuncAdd = Addition; + FuncMul = Multiplication; - Standard_Real l1, l2, lm, lr, l; //boundary curve parametrization - Standard_Real v1, v2, vm, vr, v; //Face parametrization in v direction - Standard_Real u1, u2, um, ur, u; - Standard_Real ds; //Jacobien (x, y, z) -> (u, v) = ||n|| + Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz; + dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0; - gp_Pnt P; //On the Face - gp_Vec VNor; + Standard_Real x, y, z; + Standard_Integer NbCGaussgp_Pnts = 0; - gp_Pnt2d Puv; //On the boundary curve u-v - gp_Vec2d Vuv; - Standard_Real Dul; // Dul = Du / Dl - Standard_Real CArea, CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz; - Standard_Real LocArea, LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy, - LocIxz, LocIyz; + Standard_Real l1, l2, lm, lr, l; //boundary curve parametrization + Standard_Real v1, v2, vm, vr, v; //Face parametrization in v direction + Standard_Real u1, u2, um, ur, u; + Standard_Real ds; //Jacobien (x, y, z) -> (u, v) = ||n|| + + gp_Pnt P; //On the Face + gp_Vec VNor; + + gp_Pnt2d Puv; //On the boundary curve u-v + gp_Vec2d Vuv; + Standard_Real Dul; // Dul = Du / Dl + Standard_Real CArea, CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz; + Standard_Real LocArea, LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy, + LocIxz, LocIyz; - S.Bounds (u1, u2, v1, v2); + S.Bounds (u1, u2, v1, v2); - Standard_Integer NbUGaussgp_Pnts = Min(S.UIntegrationOrder (), - math::GaussPointsMax()); - Standard_Integer NbVGaussgp_Pnts = Min(S.VIntegrationOrder (), - math::GaussPointsMax()); - - Standard_Integer NbGaussgp_Pnts = Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts); - - //Number of Gauss points for the integration - //on the Face - math_Vector GaussSPV (1, NbGaussgp_Pnts); - math_Vector GaussSWV (1, NbGaussgp_Pnts); - math::GaussPoints (NbGaussgp_Pnts,GaussSPV); - math::GaussWeights (NbGaussgp_Pnts,GaussSWV); + if(Precision::IsInfinite(u1) || Precision::IsInfinite(u2) || + Precision::IsInfinite(v1) || Precision::IsInfinite(v2)) + { + FuncAdd = AdditionInf; + FuncMul = MultiplicationInf; + } - //location point used to compute the inertia - Standard_Real xloc, yloc, zloc; - loc.Coord (xloc, yloc, zloc); + Standard_Integer NbUGaussgp_Pnts = Min(S.UIntegrationOrder (), + math::GaussPointsMax()); + Standard_Integer NbVGaussgp_Pnts = Min(S.VIntegrationOrder (), + math::GaussPointsMax()); - while (D.More()) { + Standard_Integer NbGaussgp_Pnts = Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts); - S.Load(D.Value()); - NbCGaussgp_Pnts = Min(S.IntegrationOrder (), math::GaussPointsMax()); - - math_Vector GaussCP (1, NbCGaussgp_Pnts); - math_Vector GaussCW (1, NbCGaussgp_Pnts); - math::GaussPoints (NbCGaussgp_Pnts,GaussCP); - math::GaussWeights (NbCGaussgp_Pnts,GaussCW); + //Number of Gauss points for the integration + //on the Face + math_Vector GaussSPV (1, NbGaussgp_Pnts); + math_Vector GaussSWV (1, NbGaussgp_Pnts); + math::GaussPoints (NbGaussgp_Pnts,GaussSPV); + math::GaussWeights (NbGaussgp_Pnts,GaussSWV); - CArea = 0.0; - CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0; - l1 = S.FirstParameter (); - l2 = S.LastParameter (); - lm = 0.5 * (l2 + l1); - lr = 0.5 * (l2 - l1); - Puv = S.Value2d (lm); - vm = Puv.Y(); - Puv = S.Value2d (lr); - vr = Puv.Y(); + //location point used to compute the inertia + Standard_Real xloc, yloc, zloc; + loc.Coord (xloc, yloc, zloc); - for (Standard_Integer i = 1; i <= NbCGaussgp_Pnts; i++) { - l = lm + lr * GaussCP (i); - S.D12d(l, Puv, Vuv); - v = Puv.Y(); - u2 = Puv.X(); - Dul = Vuv.Y(); - Dul *= GaussCW (i); - um = 0.5 * (u2 + u1); - ur = 0.5 * (u2 - u1); - LocArea = LocIx = LocIy = LocIz = LocIxx = LocIyy = LocIzz = + while (D.More()) { + + S.Load(D.Value()); + NbCGaussgp_Pnts = Min(S.IntegrationOrder (), math::GaussPointsMax()); + + math_Vector GaussCP (1, NbCGaussgp_Pnts); + math_Vector GaussCW (1, NbCGaussgp_Pnts); + math::GaussPoints (NbCGaussgp_Pnts,GaussCP); + math::GaussWeights (NbCGaussgp_Pnts,GaussCW); + + CArea = 0.0; + CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0; + l1 = S.FirstParameter (); + l2 = S.LastParameter (); + lm = 0.5 * (l2 + l1); + lr = 0.5 * (l2 - l1); + + Puv = S.Value2d (lm); + vm = Puv.Y(); + Puv = S.Value2d (lr); + vr = Puv.Y(); + + for (Standard_Integer i = 1; i <= NbCGaussgp_Pnts; i++) { + l = lm + lr * GaussCP (i); + S.D12d(l, Puv, Vuv); + v = Puv.Y(); + u2 = Puv.X(); + Dul = Vuv.Y(); + Dul *= GaussCW (i); + um = 0.5 * (u2 + u1); + ur = 0.5 * (u2 - u1); + LocArea = LocIx = LocIy = LocIz = LocIxx = LocIyy = LocIzz = LocIxy = LocIxz = LocIyz = 0.0; - for (Standard_Integer j = 1; j <= NbGaussgp_Pnts; j++) { - u = um + ur * GaussSPV (j); - S.Normal (u, v, P, VNor); - ds = VNor.Magnitude(); //normal.Magnitude - ds = ds * Dul * GaussSWV (j); - LocArea += ds; - P.Coord (x, y, z); - x -= xloc; - y -= yloc; - z -= zloc; - LocIx += x * ds; - LocIy += y * ds; - LocIz += z * ds; - LocIxy += x * y * ds; - LocIyz += y * z * ds; - LocIxz += x * z * ds; - x *= x; - y *= y; - z *= z; - LocIxx += (y + z) * ds; - LocIyy += (x + z) * ds; - LocIzz += (x + y) * ds; - } - CArea += LocArea * ur; - CIx += LocIx * ur; - CIy += LocIy * ur; - CIz += LocIz * ur; - CIxx += LocIxx * ur; - CIyy += LocIyy * ur; - CIzz += LocIzz * ur; - CIxy += LocIxy * ur; - CIxz += LocIxz * ur; - CIyz += LocIyz * ur; + for (Standard_Integer j = 1; j <= NbGaussgp_Pnts; j++) { + u = FuncAdd(um, FuncMul(ur, GaussSPV (j))); + S.Normal (u, v, P, VNor); + ds = VNor.Magnitude(); //normal.Magnitude + ds = FuncMul(ds, Dul) * GaussSWV (j); + LocArea = FuncAdd(LocArea, ds); + P.Coord (x, y, z); + + x = FuncAdd(x, -xloc); + y = FuncAdd(y, -yloc); + z = FuncAdd(z, -zloc); + + const Standard_Real XdS = FuncMul(x, ds); + const Standard_Real YdS = FuncMul(y, ds); + const Standard_Real ZdS = FuncMul(z, ds); + + LocIx = FuncAdd(LocIx, XdS); + LocIy = FuncAdd(LocIy, YdS); + LocIz = FuncAdd(LocIz, ZdS); + LocIxy = FuncAdd(LocIxy, FuncMul(x, YdS)); + LocIyz = FuncAdd(LocIyz, FuncMul(y, ZdS)); + LocIxz = FuncAdd(LocIxz, FuncMul(x, ZdS)); + x = Precision::IsInfinite(x) ? Precision::Infinite() : x*x; + y = Precision::IsInfinite(y) ? Precision::Infinite() : y*y; + z = Precision::IsInfinite(z) ? Precision::Infinite() : z*z; + LocIxx = FuncAdd(LocIxx, FuncAdd(YdS, ZdS)); + LocIyy = FuncAdd(LocIyy, FuncAdd(XdS, ZdS)); + LocIzz = FuncAdd(LocIzz, FuncAdd(XdS, YdS)); } - dim += CArea * lr; - Ix += CIx * lr; - Iy += CIy * lr; - Iz += CIz * lr; - Ixx += CIxx * lr; - Iyy += CIyy * lr; - Izz += CIzz * lr; - Ixy += CIxy * lr; - Ixz += CIxz * lr; - Iyz += CIyz * lr; - D.Next(); - } - if (Abs(dim) >= EPS_DIM) { - Ix /= dim; - Iy /= dim; - Iz /= dim; - g.SetCoord (Ix, Iy, Iz); - } - else { - dim =0.; - g.SetCoord (0., 0.,0.); - } - inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz), - gp_XYZ (-Ixy, Iyy, -Iyz), - gp_XYZ (-Ixz, -Iyz, Izz)); + + CArea = FuncAdd(CArea, FuncMul(LocArea, ur)); + CIx = FuncAdd(CIx, FuncMul(LocIx, ur)); + CIy = FuncAdd(CIy, FuncMul(LocIy, ur)); + CIz = FuncAdd(CIz, FuncMul(LocIz, ur)); + CIxx = FuncAdd(CIxx, FuncMul(LocIxx, ur)); + CIyy = FuncAdd(CIyy, FuncMul(LocIyy, ur)); + CIzz = FuncAdd(CIzz, FuncMul(LocIzz, ur)); + CIxy = FuncAdd(CIxy, FuncMul(LocIxy, ur)); + CIxz = FuncAdd(CIxz, FuncMul(LocIxz, ur)); + CIyz = FuncAdd(CIyz, FuncMul(LocIyz, ur)); + } + + dim = FuncAdd(dim, FuncMul(CArea, lr)); + Ix = FuncAdd(Ix, FuncMul(CIx, lr)); + Iy = FuncAdd(Iy, FuncMul(CIy, lr)); + Iz = FuncAdd(Iz, FuncMul(CIz, lr)); + Ixx = FuncAdd(Ixx, FuncMul(CIxx, lr)); + Iyy = FuncAdd(Iyy, FuncMul(CIyy, lr)); + Izz = FuncAdd(Izz, FuncMul(CIzz, lr)); + Ixy = FuncAdd(Ixy, FuncMul(CIxy, lr)); + Ixz = FuncAdd(Iyz, FuncMul(CIxz, lr)); + Iyz = FuncAdd(Ixz, FuncMul(CIyz, lr)); + D.Next(); + } + + if (Abs(dim) >= EPS_DIM) { + Ix /= dim; + Iy /= dim; + Iz /= dim; + g.SetCoord (Ix, Iy, Iz); + } + else { + dim =0.; + g.SetCoord (0., 0.,0.); + } + + inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz), + gp_XYZ (-Ixy, Iyy, -Iyz), + gp_XYZ (-Ixz, -Iyz, Izz)); } +static void Compute(const Face& S, + const gp_Pnt& loc, + Standard_Real& dim, + gp_Pnt& g, + gp_Mat& inertia) +{ + Standard_Real (*FuncAdd)(Standard_Real, Standard_Real); + Standard_Real (*FuncMul)(Standard_Real, Standard_Real); - -static void Compute(const Face& S, const gp_Pnt& loc, Standard_Real& dim, gp_Pnt& g, gp_Mat& inertia){ - Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz; - dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0; + FuncAdd = Addition; + FuncMul = Multiplication; - Standard_Real LowerU, UpperU, LowerV, UpperV; - S.Bounds (LowerU, UpperU, LowerV, UpperV); - Standard_Integer UOrder = Min(S.UIntegrationOrder (), - math::GaussPointsMax()); - Standard_Integer VOrder = Min(S.VIntegrationOrder (), - math::GaussPointsMax()); - gp_Pnt P; - gp_Vec VNor; - Standard_Real dsi, ds; - Standard_Real ur, um, u, vr, vm, v; - Standard_Real x, y, z; - Standard_Real Ixi, Iyi, Izi, Ixxi, Iyyi, Izzi, Ixyi, Ixzi, Iyzi; - Standard_Real xloc, yloc, zloc; - loc.Coord (xloc, yloc, zloc); + Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz; + dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0; - Standard_Integer i, j; - math_Vector GaussPU (1, UOrder); //gauss points and weights - math_Vector GaussWU (1, UOrder); - math_Vector GaussPV (1, VOrder); - math_Vector GaussWV (1, VOrder); + Standard_Real LowerU, UpperU, LowerV, UpperV; + S.Bounds (LowerU, UpperU, LowerV, UpperV); - //Recuperation des points de Gauss dans le fichier GaussPoints. - math::GaussPoints (UOrder,GaussPU); - math::GaussWeights (UOrder,GaussWU); - math::GaussPoints (VOrder,GaussPV); - math::GaussWeights (VOrder,GaussWV); + if(Precision::IsInfinite(LowerU) || Precision::IsInfinite(UpperU) || + Precision::IsInfinite(LowerV) || Precision::IsInfinite(UpperV)) + { + FuncAdd = AdditionInf; + FuncMul = MultiplicationInf; + } - // Calcul des integrales aux points de gauss : - um = 0.5 * (UpperU + LowerU); - vm = 0.5 * (UpperV + LowerV); - ur = 0.5 * (UpperU - LowerU); - vr = 0.5 * (UpperV - LowerV); + Standard_Integer UOrder = Min(S.UIntegrationOrder (), + math::GaussPointsMax()); + Standard_Integer VOrder = Min(S.VIntegrationOrder (), + math::GaussPointsMax()); + gp_Pnt P; + gp_Vec VNor; + Standard_Real dsi, ds; + Standard_Real ur, um, u, vr, vm, v; + Standard_Real x, y, z; + Standard_Real Ixi, Iyi, Izi, Ixxi, Iyyi, Izzi, Ixyi, Ixzi, Iyzi; + Standard_Real xloc, yloc, zloc; + loc.Coord (xloc, yloc, zloc); - for (j = 1; j <= VOrder; j++) { - v = vm + vr * GaussPV (j); - dsi = Ixi = Iyi = Izi = Ixxi = Iyyi = Izzi = Ixyi = Ixzi = Iyzi = 0.0; + Standard_Integer i, j; + math_Vector GaussPU (1, UOrder); //gauss points and weights + math_Vector GaussWU (1, UOrder); + math_Vector GaussPV (1, VOrder); + math_Vector GaussWV (1, VOrder); - for (i = 1; i <= UOrder; i++) { - u = um + ur * GaussPU (i); - S.Normal (u, v, P, VNor); - ds = VNor.Magnitude() * GaussWU (i); - P.Coord (x, y, z); - x -= xloc; - y -= yloc; - z -= zloc; - dsi += ds; - Ixi += x * ds; - Iyi += y * ds; - Izi += z * ds; - Ixyi += x * y * ds; - Iyzi += y * z * ds; - Ixzi += x * z * ds; - x *= x; - y *= y; - z *= z; - Ixxi += (y + z) * ds; - Iyyi += (x + z) * ds; - Izzi += (x + y) * ds; - } - dim += dsi * GaussWV (j); - Ix += Ixi * GaussWV (j); - Iy += Iyi * GaussWV (j); - Iz += Izi * GaussWV (j); - Ixx += Ixxi * GaussWV (j); - Iyy += Iyyi * GaussWV (j); - Izz += Izzi * GaussWV (j); - Ixy += Ixyi * GaussWV (j); - Iyz += Iyzi * GaussWV (j); - Ixz += Ixzi * GaussWV (j); - } - vr *= ur; - Ixx *= vr; - Iyy *= vr; - Izz *= vr; - Ixy *= vr; - Ixz *= vr; - Iyz *= vr; - if (Abs(dim) >= EPS_DIM) { - Ix /= dim; - Iy /= dim; - Iz /= dim; - dim *= vr; - g.SetCoord (Ix, Iy, Iz); - } - else { - dim =0.; - g.SetCoord (0.,0.,0.); - } - inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz), - gp_XYZ (-Ixy, Iyy, -Iyz), - gp_XYZ (-Ixz, -Iyz, Izz)); + //Recuperation des points de Gauss dans le fichier GaussPoints. + math::GaussPoints (UOrder,GaussPU); + math::GaussWeights (UOrder,GaussWU); + math::GaussPoints (VOrder,GaussPV); + math::GaussWeights (VOrder,GaussWV); + + // Calcul des integrales aux points de gauss : + um = 0.5 * FuncAdd(UpperU, LowerU); + vm = 0.5 * FuncAdd(UpperV, LowerV); + ur = 0.5 * FuncAdd(UpperU, -LowerU); + vr = 0.5 * FuncAdd(UpperV, -LowerV); + + for (j = 1; j <= VOrder; j++) { + v = FuncAdd(vm, FuncMul(vr, GaussPV(j))); + dsi = Ixi = Iyi = Izi = Ixxi = Iyyi = Izzi = Ixyi = Ixzi = Iyzi = 0.0; + + for (i = 1; i <= UOrder; i++) { + u = FuncAdd(um, FuncMul(ur, GaussPU (i))); + S.Normal (u, v, P, VNor); + ds = FuncMul(VNor.Magnitude(), GaussWU (i)); + P.Coord (x, y, z); + + x = FuncAdd(x, -xloc); + y = FuncAdd(y, -yloc); + z = FuncAdd(z, -zloc); + + dsi = FuncAdd(dsi, ds); + + const Standard_Real XdS = FuncMul(x, ds); + const Standard_Real YdS = FuncMul(y, ds); + const Standard_Real ZdS = FuncMul(z, ds); + + Ixi = FuncAdd(Ixi, XdS); + Iyi = FuncAdd(Iyi, YdS); + Izi = FuncAdd(Izi, ZdS); + Ixyi = FuncAdd(Ixyi, FuncMul(x, YdS)); + Iyzi = FuncAdd(Iyzi, FuncMul(y, ZdS)); + Ixzi = FuncAdd(Ixzi, FuncMul(x, ZdS)); + x = Precision::IsInfinite(x) ? Precision::Infinite() : x*x; + y = Precision::IsInfinite(y) ? Precision::Infinite() : y*y; + z = Precision::IsInfinite(z) ? Precision::Infinite() : z*z; + Ixxi = FuncAdd(Ixxi, FuncAdd(YdS, ZdS)); + Iyyi = FuncAdd(Iyyi, FuncAdd(XdS, ZdS)); + Izzi = FuncAdd(Izzi, FuncAdd(XdS, YdS)); + } + + dim = FuncAdd(dim, FuncMul(dsi, GaussWV (j))); + Ix = FuncAdd(Ix, FuncMul(Ixi, GaussWV (j))); + Iy = FuncAdd(Iy, FuncMul(Iyi, GaussWV (j))); + Iz = FuncAdd(Iz, FuncMul(Izi, GaussWV (j))); + Ixx = FuncAdd(Ixx, FuncMul(Ixxi, GaussWV (j))); + Iyy = FuncAdd(Iyy, FuncMul(Iyyi, GaussWV (j))); + Izz = FuncAdd(Izz, FuncMul(Izzi, GaussWV (j))); + Ixy = FuncAdd(Ixy, FuncMul(Ixyi, GaussWV (j))); + Iyz = FuncAdd(Iyz, FuncMul(Iyzi, GaussWV (j))); + Ixz = FuncAdd(Ixz, FuncMul(Ixzi, GaussWV (j))); + } + + vr = FuncMul(vr, ur); + Ixx = FuncMul(vr, Ixx); + Iyy = FuncMul(vr, Iyy); + Izz = FuncMul(vr, Izz); + Ixy = FuncMul(vr, Ixy); + Ixz = FuncMul(vr, Ixz); + Iyz = FuncMul(vr, Iyz); + + if (Abs(dim) >= EPS_DIM) + { + Ix /= dim; + Iy /= dim; + Iz /= dim; + dim *= vr; + g.SetCoord (Ix, Iy, Iz); + } + else + { + dim =0.; + g.SetCoord (0.,0.,0.); + } + + inertia = gp_Mat (gp_XYZ ( Ixx, -Ixy, -Ixz), + gp_XYZ (-Ixy, Iyy, -Iyz), + gp_XYZ (-Ixz, -Iyz, Izz)); } GProp_SGProps::GProp_SGProps(){} GProp_SGProps::GProp_SGProps (const Face& S, - const gp_Pnt& SLocation - ) + const gp_Pnt& SLocation + ) { - SetLocation(SLocation); - Perform(S); + SetLocation(SLocation); + Perform(S); } GProp_SGProps::GProp_SGProps (Face& S, Domain& D, - const gp_Pnt& SLocation - ) + const gp_Pnt& SLocation + ) { - SetLocation(SLocation); - Perform(S,D); + SetLocation(SLocation); + Perform(S,D); } GProp_SGProps::GProp_SGProps(Face& S, const gp_Pnt& SLocation, const Standard_Real Eps){ diff --git a/src/IntImp/IntImp_Int2S.gxx b/src/IntImp/IntImp_Int2S.gxx index 73a4dc5ea7..a3d45db05f 100755 --- a/src/IntImp/IntImp_Int2S.gxx +++ b/src/IntImp/IntImp_Int2S.gxx @@ -30,7 +30,7 @@ #include #include #include - + //Standard_IMPORT extern IntImp_ConstIsoparametric *ChoixRef; @@ -41,72 +41,72 @@ IntImp_Int2S::IntImp_Int2S() { IntImp_Int2S::IntImp_Int2S(const ThePSurface& surf1, - const ThePSurface& surf2, - const Standard_Real TolTangency ) : - done(Standard_True), - empty(Standard_True), - myZerParFunc(surf1,surf2), - tol(TolTangency*TolTangency) + const ThePSurface& surf2, + const Standard_Real TolTangency ) : +done(Standard_True), +empty(Standard_True), +myZerParFunc(surf1,surf2), +tol(TolTangency*TolTangency) { ua0 = ThePSurfaceTool::FirstUParameter(surf1); //-- ThePSurfaceTool::UIntervalFirst(surf1); va0 = ThePSurfaceTool::FirstVParameter(surf1); //-- ThePSurfaceTool::VIntervalFirst(surf1); ua1 = ThePSurfaceTool::LastUParameter(surf1); //-- ThePSurfaceTool::UIntervalLast(surf1); va1 = ThePSurfaceTool::LastVParameter(surf1); //-- ThePSurfaceTool::VIntervalLast(surf1); - + ub0 = ThePSurfaceTool::FirstUParameter(surf2); //-- ThePSurfaceTool::UIntervalFirst(surf2); vb0 = ThePSurfaceTool::FirstVParameter(surf2); //-- ThePSurfaceTool::VIntervalFirst(surf2); ub1 = ThePSurfaceTool::LastUParameter(surf2); //-- ThePSurfaceTool::UIntervalLast(surf2); vb1 = ThePSurfaceTool::LastVParameter(surf2); //-- ThePSurfaceTool::VIntervalLast(surf2); - + ures1 = ThePSurfaceTool::UResolution(surf1,Precision::Confusion()); vres1 = ThePSurfaceTool::VResolution(surf1,Precision::Confusion()); - + ures2 = ThePSurfaceTool::UResolution(surf2,Precision::Confusion()); vres2 = ThePSurfaceTool::VResolution(surf2,Precision::Confusion()); } IntImp_Int2S::IntImp_Int2S(const TColStd_Array1OfReal& Param, - const ThePSurface& surf1, - const ThePSurface& surf2, - const Standard_Real TolTangency ) : - done(Standard_True), - empty(Standard_True), - myZerParFunc(surf1,surf2), - tol(TolTangency*TolTangency) + const ThePSurface& surf1, + const ThePSurface& surf2, + const Standard_Real TolTangency ) : +done(Standard_True), +empty(Standard_True), +myZerParFunc(surf1,surf2), +tol(TolTangency*TolTangency) { math_FunctionSetRoot Rsnld(myZerParFunc,15); //-- Modif lbr 18 MAI ????????????? ua0 = ThePSurfaceTool::FirstUParameter(surf1); //-- ThePSurfaceTool::UIntervalFirst(surf1); va0 = ThePSurfaceTool::FirstVParameter(surf1); //-- ThePSurfaceTool::VIntervalFirst(surf1); ua1 = ThePSurfaceTool::LastUParameter(surf1); //-- ThePSurfaceTool::UIntervalLast(surf1); va1 = ThePSurfaceTool::LastVParameter(surf1); //-- ThePSurfaceTool::VIntervalLast(surf1); - + ub0 = ThePSurfaceTool::FirstUParameter(surf2); //-- ThePSurfaceTool::UIntervalFirst(surf2); vb0 = ThePSurfaceTool::FirstVParameter(surf2); //-- ThePSurfaceTool::VIntervalFirst(surf2); ub1 = ThePSurfaceTool::LastUParameter(surf2); //-- ThePSurfaceTool::UIntervalLast(surf2); vb1 = ThePSurfaceTool::LastVParameter(surf2); //-- ThePSurfaceTool::VIntervalLast(surf2); - + ures1 = ThePSurfaceTool::UResolution(surf1,Precision::Confusion()); vres1 = ThePSurfaceTool::VResolution(surf1,Precision::Confusion()); - + ures2 = ThePSurfaceTool::UResolution(surf2,Precision::Confusion()); vres2 = ThePSurfaceTool::VResolution(surf2,Precision::Confusion()); Perform(Param,Rsnld); } -IntImp_ConstIsoparametric IntImp_Int2S:: - Perform(const TColStd_Array1OfReal& Param, - math_FunctionSetRoot& Rsnld, - const IntImp_ConstIsoparametric ChoixIso ) +IntImp_ConstIsoparametric IntImp_Int2S:: Perform(const TColStd_Array1OfReal& Param, + math_FunctionSetRoot& Rsnld, + const IntImp_ConstIsoparametric ChoixIso) { Standard_Real BornInfBuf[3], BornSupBuf[3], ToleranceBuf[3], UVapBuf[3]; Standard_Real UvresBuf[4]; - math_Vector BornInf (BornInfBuf, 1, 3), BornSup (BornSupBuf, 1, 3), Tolerance (ToleranceBuf, 1, 3), UVap (UVapBuf, 1, 3); + math_Vector BornInf (BornInfBuf, 1, 3), BornSup (BornSupBuf, 1, 3), + Tolerance (ToleranceBuf, 1, 3), UVap (UVapBuf, 1, 3); TColStd_Array1OfReal Uvres (UvresBuf[0], 1, 4); IntImp_ConstIsoparametric BestChoix; - myZerParFunc.ComputeParameters(ChoixIso,Param,UVap, - BornInf,BornSup,Tolerance); + + myZerParFunc.ComputeParameters(ChoixIso,Param,UVap,BornInf,BornSup,Tolerance); Rsnld.SetTolerance(Tolerance); Rsnld.Perform(myZerParFunc,UVap,BornInf,BornSup); BestChoix = ChoixIso; @@ -118,9 +118,9 @@ IntImp_ConstIsoparametric IntImp_Int2S:: tangent = myZerParFunc.IsTangent(UVap,Uvres,BestChoix); pint.SetValue(myZerParFunc.Point(),Uvres(1),Uvres(2),Uvres(3),Uvres(4)); if (!tangent) { - d3d = myZerParFunc.Direction(); - d2d1 = myZerParFunc.DirectionOnS1(); - d2d2 = myZerParFunc.DirectionOnS2(); + d3d = myZerParFunc.Direction(); + d2d1 = myZerParFunc.DirectionOnS1(); + d2d2 = myZerParFunc.DirectionOnS2(); } } else { @@ -133,9 +133,8 @@ IntImp_ConstIsoparametric IntImp_Int2S:: return BestChoix; } -IntImp_ConstIsoparametric IntImp_Int2S:: Perform( - const TColStd_Array1OfReal& Param, - math_FunctionSetRoot& Rsnld) +IntImp_ConstIsoparametric IntImp_Int2S:: Perform(const TColStd_Array1OfReal& Param, + math_FunctionSetRoot& Rsnld) { gp_Vec DPUV[4]; gp_Pnt P1,P2; @@ -157,15 +156,20 @@ IntImp_ConstIsoparametric IntImp_Int2S:: Perform( Epsuv[2] = ThePSurfaceTool::UResolution(Caro2,Precision::Confusion()); Epsuv[3] = ThePSurfaceTool::VResolution(Caro2,Precision::Confusion()); - for (Standard_Integer j=0;j<=3;j++) UVd[j] = Param(j+1); + for (Standard_Integer j=0;j<=3;j++) + UVd[j] = Param(j+1); empty = Standard_True; Standard_Boolean Tangent = IntImp_ComputeTangence(DPUV,Epsuv,UVd,ChoixIso); - if (Tangent) return BestChoix; + if (Tangent) + return BestChoix; + Standard_Integer i=0; IntImp_ConstIsoparametric CurrentChoix = BestChoix; //-- Modif 17 Mai 93 - while (empty && i<= 3) { + + while (empty && i<= 3) + { CurrentChoix = Perform(Param,Rsnld,ChoixIso[i]); if(!empty) { BestChoix = CurrentChoix; @@ -174,7 +178,7 @@ IntImp_ConstIsoparametric IntImp_Int2S:: Perform( } if (!empty) { // verifier que l on ne deborde pas les frontieres pint.Parameters(Duv(1),Duv(2),Duv(3),Duv(4)); - + UVd[0] = ua0; //-- ThePSurfaceTool::UIntervalFirst(Caro1); UVd[1] = va0; //-- ThePSurfaceTool::VIntervalFirst(Caro1); UVf[0] = ua1; //-- ThePSurfaceTool::UIntervalLast(Caro1); @@ -184,7 +188,7 @@ IntImp_ConstIsoparametric IntImp_Int2S:: Perform( UVd[3] = vb0; //-- ubThePSurfaceTool::VIntervalFirst(Caro2); UVf[2] = ub1; //-- ThePSurfaceTool::UIntervalLast(Caro2); UVf[3] = vb1; //-- ThePSurfaceTool::VIntervalLast(Caro2); - + Standard_Integer Nc,Iiso; if (Duv(1) <= UVd[0]-Epsuv[0]) { Duv(1) = UVd[0]; @@ -233,18 +237,28 @@ IntImp_ConstIsoparametric IntImp_Int2S:: Perform( if (!empty) { // verification si l on ne deborde pas sur le carreau // reciproque Nc =3-Nc; - if (Duv(Nc) <= UVd[Nc-1]-Epsuv[Nc-1]) Duv(Nc)=UVd[Nc-1]; - else if (Duv(Nc) >=UVf[Nc-1]+Epsuv[Nc-1]) Duv(Nc)=UVf[Nc-1]; - else if (Duv(Nc+1) <= UVd[Nc]) { - Nc = Nc +1; - Duv(Nc)=UVd[Nc-1]; + if (Duv(Nc) <= UVd[Nc-1]-Epsuv[Nc-1]) + Duv(Nc)=UVd[Nc-1]; + else if (Duv(Nc) >=UVf[Nc-1]+Epsuv[Nc-1]) + Duv(Nc)=UVf[Nc-1]; + else if (Duv(Nc+1) <= UVd[Nc]) + { + Nc = Nc + 1; + Duv(Nc)=UVd[Nc-1]; } - else if (Duv(Nc+1) >=UVf[Nc]) { - Nc = Nc +1; - Duv(Nc)=UVf[Nc-1]; + else if (Duv(Nc+1) >=UVf[Nc]) + { + Nc = Nc + 1; + Duv(Nc)=UVf[Nc-1]; } - else return BestChoix; + else + return BestChoix; + empty = Standard_True; + + if(Nc == 4) + Nc = 0; + BestChoix = ChoixRef[Nc]; //en attendant BestChoix = Perform(Duv,Rsnld,BestChoix); } diff --git a/src/IntWalk/IntWalk_PWalking_1.gxx b/src/IntWalk/IntWalk_PWalking_1.gxx index 84a8023af5..f4c6c9d95e 100755 --- a/src/IntWalk/IntWalk_PWalking_1.gxx +++ b/src/IntWalk/IntWalk_PWalking_1.gxx @@ -379,11 +379,9 @@ Standard_Boolean IntWalk_PWalking::PerformFirstPoint (const TColStd_Array1OfRea close = Standard_False; // Standard_Integer i; - Standard_Real aTmp; TColStd_Array1OfReal Param(1,4); // for (i=1; i<=4; ++i) { - aTmp = ParDep(i); Param(i) = ParDep(i); } //-- calculate the first solution point @@ -393,9 +391,11 @@ Standard_Boolean IntWalk_PWalking::PerformFirstPoint (const TColStd_Array1OfRea if (!myIntersectionOn2S.IsDone()) { return Standard_False; } + if (myIntersectionOn2S.IsEmpty()) { return Standard_False; } + FirstPoint = myIntersectionOn2S.Point(); return Standard_True; } diff --git a/tests/bugs/modalg_1/bug19793_2 b/tests/bugs/modalg_1/bug19793_2 index 8f85c95d5d..55dcd43366 100755 --- a/tests/bugs/modalg_1/bug19793_2 +++ b/tests/bugs/modalg_1/bug19793_2 @@ -9,7 +9,7 @@ puts "" # Fuse problem of symetrical shapes. Appendix for NPAL19789 ####################################################################### -cpulimit 500 +cpulimit 1000 #cpulimit 4500 set BugNumber OCC19793