diff --git a/src/AdvApprox/AdvApprox_ApproxAFunction.cxx b/src/AdvApprox/AdvApprox_ApproxAFunction.cxx index 5bbecf7d45..af57b2a019 100755 --- a/src/AdvApprox/AdvApprox_ApproxAFunction.cxx +++ b/src/AdvApprox/AdvApprox_ApproxAFunction.cxx @@ -786,20 +786,20 @@ void AdvApprox_ApproxAFunction::Perform(const Standard_Integer Num1DSS, TColStd_Array1OfReal AverageError(1,myMaxSegments * TotalNumSS) ; - Approximation (TotalDimension, TotalNumSS, LocalDimension, - myFirst, myLast, - *(AdvApprox_EvaluatorFunction*)myEvaluator, - CutTool, - ContinuityOrder, - NumMaxCoeffs, myMaxSegments, - LocalTolerances, code_precis, - NumCurves, // Nombre de courbe en sortie - NumCoeffPerCurvePtr->ChangeArray1(), // Nbre de coeff par courbe - LocalCoefficientsPtr->ChangeArray1(),// Les Coeffs solutions - IntervalsPtr->ChangeArray1(), // La Table des decoupes - ErrorMax, // Majoration de l'erreur - AverageError, // Erreur moyenne constatee - ErrorCode) ; + Approximation ( TotalDimension, + TotalNumSS, LocalDimension, + myFirst, myLast, + *(AdvApprox_EvaluatorFunction*)myEvaluator, + CutTool, ContinuityOrder, NumMaxCoeffs, + myMaxSegments, LocalTolerances, code_precis, + NumCurves, // Nombre de courbe en sortie + NumCoeffPerCurvePtr->ChangeArray1(), // Nbre de coeff par courbe + LocalCoefficientsPtr->ChangeArray1(), // Les Coeffs solutions + IntervalsPtr->ChangeArray1(), // La Table des decoupes + ErrorMax, // Majoration de l'erreur + AverageError, // Erreur moyenne constatee + ErrorCode) ; + if (ErrorCode == 0 || ErrorCode == -1) { // // si tout est OK ou bien on a un resultat dont l une des erreurs max est diff --git a/src/Approx/Approx_SweepApproximation.cxx b/src/Approx/Approx_SweepApproximation.cxx index ff478da58c..2743fd499b 100755 --- a/src/Approx/Approx_SweepApproximation.cxx +++ b/src/Approx/Approx_SweepApproximation.cxx @@ -478,11 +478,11 @@ Standard_Boolean Approx_SweepApproximation::D1(const Standard_Real Param, myPoles->ChangeValue(ii).ChangeCoord() -= Translation.XYZ(); // Homothety on all. - myDPoles->ChangeValue(ii) *= myWeigths->Value(ii); + const Standard_Real aWeight = myWeigths->Value(ii); + myDPoles->ChangeValue(ii) *= aWeight; Vaux.SetXYZ( myPoles->Value(ii).Coord()); myDPoles->ChangeValue(ii) += myDWeigths->Value(ii)*Vaux; - myPoles->ChangeValue(ii).ChangeCoord() - *= myWeigths->Value(ii); // for the cash + myPoles->ChangeValue(ii).ChangeCoord() *= aWeight; // for the cash } diff --git a/src/DrawTrSurf/DrawTrSurf_Drawable.cxx b/src/DrawTrSurf/DrawTrSurf_Drawable.cxx index 672bd4383c..3a966e7ab0 100755 --- a/src/DrawTrSurf/DrawTrSurf_Drawable.cxx +++ b/src/DrawTrSurf/DrawTrSurf_Drawable.cxx @@ -127,21 +127,25 @@ static void PlotCurve (Draw_Display& aDisplay, //======================================================================= void DrawTrSurf_Drawable::DrawCurveOn (Adaptor3d_Curve& C, - Draw_Display& aDisplay) const + Draw_Display& aDisplay) const { gp_Pnt P; - if (myDrawMode == 1) { + if (myDrawMode == 1) + { Standard_Real Fleche = myDeflection/aDisplay.Zoom(); GCPnts_UniformDeflection LineVu(C,Fleche); - if (LineVu.IsDone()) { + if (LineVu.IsDone()) + { aDisplay.MoveTo(LineVu.Value(1)); - for (Standard_Integer i = 2; i <= LineVu.NbPoints(); i++) { - aDisplay.DrawTo(LineVu.Value(i)); + for (Standard_Integer i = 2; i <= LineVu.NbPoints(); i++) + { + aDisplay.DrawTo(LineVu.Value(i)); } - } + } } - else { - Standard_Real j; + else + { + Standard_Integer j; Standard_Integer intrv, nbintv = C.NbIntervals(GeomAbs_CN); TColStd_Array1OfReal TI(1,nbintv+1); C.Intervals(TI,GeomAbs_CN); @@ -150,36 +154,44 @@ void DrawTrSurf_Drawable::DrawCurveOn (Adaptor3d_Curve& C, GeomAbs_CurveType CurvType = C.GetType(); gp_Pnt aPPnt=P, aNPnt; - for (intrv = 1; intrv <= nbintv; intrv++) { + for (intrv = 1; intrv <= nbintv; intrv++) + { Standard_Real t = TI(intrv); Standard_Real step = (TI(intrv+1) - t) / myDiscret; - switch (CurvType) { - case GeomAbs_Line : - break; - case GeomAbs_Circle : - case GeomAbs_Ellipse : - for (j = 1; j < myDiscret; j++) { - t += step; - C.D0(t,P); - aDisplay.DrawTo(P); - } - break; - case GeomAbs_Parabola : - case GeomAbs_Hyperbola : - case GeomAbs_BezierCurve : - case GeomAbs_BSplineCurve : - case GeomAbs_OtherCurve : - for (j = 1; j <= myDiscret/2; j++) { - C.D0 (t+step*2., aNPnt); + switch (CurvType) + { + case GeomAbs_Line: + break; + case GeomAbs_Circle: + case GeomAbs_Ellipse: + for (j = 1; j < myDiscret; j++) + { + t += step; + C.D0(t,P); + aDisplay.DrawTo(P); + } + break; + case GeomAbs_Parabola: + case GeomAbs_Hyperbola: + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + case GeomAbs_OtherCurve: + const Standard_Integer nIter = myDiscret/2; + for (j = 1; j < nIter; j++) + { + const Standard_Real t1 = t+step*2.; + C.D0 (t1, aNPnt); PlotCurve (aDisplay, C, t, step, aPPnt, aNPnt); - aPPnt = aNPnt; - t += step*2.; - } - break; + aPPnt = aNPnt; + t = t1; + } + + break; } C.D0(TI(intrv+1),P); + PlotCurve (aDisplay, C, t, step, aPPnt, P); aDisplay.DrawTo(P); } } diff --git a/src/Extrema/Extrema_CurveTool.cdl b/src/Extrema/Extrema_CurveTool.cdl index 12c5f662b3..13e6726171 100755 --- a/src/Extrema/Extrema_CurveTool.cdl +++ b/src/Extrema/Extrema_CurveTool.cdl @@ -83,6 +83,11 @@ is U : Real from Standard) returns Pnt from gp; ---C++: inline + + D0 (myclass; C : Curve from Adaptor3d; + U : Real from Standard; + P : out Pnt from gp); + ---C++: inline D1 (myclass; C : Curve from Adaptor3d; U : Real from Standard; @@ -90,10 +95,22 @@ is V : out Vec from gp); ---C++: inline - D2 (myclass; C : Curve from Adaptor3d; + D2 (myclass; C : Curve from Adaptor3d; + U : Real from Standard; + P : out Pnt from gp; + V1, V2 : out Vec from gp); + ---C++: inline + + D3 (myclass; C : Curve from Adaptor3d; + U : Real from Standard; + P : out Pnt from gp; + V1, V2, V3 : out Vec from gp); + ---C++: inline + + DN (myclass; C : Curve from Adaptor3d; U : Real from Standard; - P : out Pnt from gp; - V1, V2 : out Vec from gp); + N : Integer from Standard) + returns Vec from gp; ---C++: inline Line(myclass; C : Curve from Adaptor3d) returns Lin from gp; diff --git a/src/Extrema/Extrema_CurveTool.lxx b/src/Extrema/Extrema_CurveTool.lxx index 6f8490430c..85a62994ea 100755 --- a/src/Extrema/Extrema_CurveTool.lxx +++ b/src/Extrema/Extrema_CurveTool.lxx @@ -62,6 +62,17 @@ inline gp_Pnt Extrema_CurveTool::Value(const Adaptor3d_Curve& C, return C.Value(U); } +//======================================================================= +//function : D0 +//purpose : +//======================================================================= + +inline void Extrema_CurveTool::D0(const Adaptor3d_Curve& C, + const Standard_Real U, + gp_Pnt& P) +{ + C.D0(U, P); +} //======================================================================= //function : D1 @@ -90,6 +101,31 @@ inline void Extrema_CurveTool::D2(const Adaptor3d_Curve& C, C.D2(U, P, V1, V2); } +//======================================================================= +//function : D3 +//purpose : +//======================================================================= + + inline void Extrema_CurveTool::D3(const Adaptor3d_Curve& C, + const Standard_Real U, + gp_Pnt& P, + gp_Vec& V1, + gp_Vec& V2, + gp_Vec& V3) +{ + C.D3(U, P, V1, V2, V3); +} + +//======================================================================= +//function : DN +//purpose : +//======================================================================= +inline gp_Vec Extrema_CurveTool::DN(const Adaptor3d_Curve& C, + const Standard_Real U, + const Standard_Integer N) +{ + return C.DN(U, N); +} //======================================================================= diff --git a/src/Extrema/Extrema_FuncExtCC.cdl b/src/Extrema/Extrema_FuncExtCC.cdl index e49d6b7292..45ebdccec7 100755 --- a/src/Extrema/Extrema_FuncExtCC.cdl +++ b/src/Extrema/Extrema_FuncExtCC.cdl @@ -50,7 +50,6 @@ is ---Purpose: SetCurve (me: in out; theRank: Integer; C1: Curve1); - ---C++: inline ---Purpose: SetTolerance (me: in out; theTol: Real); @@ -103,6 +102,14 @@ is ---Purpose: Returns a tolerance specified in the constructor -- or in SetTolerance() method. + SubIntervalInitialize(me: in out; theUfirst, theUlast: Vector); + ---Purpose: Determines of boundaries of subinterval for find of root. + + SearchOfTolerance(me: in out; C: Address from Standard) returns Real from Standard; + ---Purpose: Computes a Tol value. If 1st derivative of curve + -- |D1| + + +static const Standard_Real MinTol = 1.e-20; +static const Standard_Real TolFactor = 1.e-12; +static const Standard_Real MinStep = 1e-7; + +static const Standard_Integer MaxOrder = 3; + + + //============================================================================= +Standard_Real Extrema_FuncExtCC::SearchOfTolerance(const Standard_Address C) + { + const Standard_Integer NPoint = 10; + Standard_Real aStartParam, anEndParam; + + if(C==myC1) + { + aStartParam = myUinfium; + anEndParam = myUsupremum; + } + else if(C==myC2) + { + aStartParam = myVinfium; + anEndParam = myVsupremum; + } + else + { +#ifdef DEB + cout << "+++ Function Extrema_FuncExtCC::SearchOfTolerance(...)" << endl; + cout << "Warning: No curve for tolerance computing!---"< anEndParam) + u = anEndParam; + + Pnt Ptemp; //empty point (is not used below) + Vec VDer; // 1st derivative vector + Tool1::D1(*((Curve1*)C), u, Ptemp, VDer); + Standard_Real vm = VDer.Magnitude(); + if(vm > aMax) + aMax = vm; + } + while(++aNum < NPoint+1); + + return Max(aMax*TolFactor,MinTol); + } -#define Tol 1.e-20 -#define delta 1.e-9 - +//============================================================================= Extrema_FuncExtCC::Extrema_FuncExtCC(const Standard_Real thetol) : myC1 (0), myC2 (0), myTol (thetol) -{ -} + { + math_Vector V1(1,2), V2(1,2); + V1(1) = 0.0; + V2(1) = 0.0; + V1(2) = 0.0; + V2(2) = 0.0; + SubIntervalInitialize(V1, V2); + myMaxDerivOrderC1 = 0; + myTolC1=MinTol; + myMaxDerivOrderC2 = 0; + myTolC2=MinTol; + } +//============================================================================= Extrema_FuncExtCC::Extrema_FuncExtCC (const Curve1& C1, const Curve2& C2, const Standard_Real thetol) : myC1 ((Standard_Address)&C1), myC2 ((Standard_Address)&C2), myTol (thetol) -{ -} + { + math_Vector V1(1,2), V2(1,2); + V1(1) = Tool1::FirstParameter(*((Curve1*)myC1)); + V2(1) = Tool1::LastParameter(*((Curve1*)myC1)); + V1(2) = Tool2::FirstParameter(*((Curve2*)myC2)); + V2(2) = Tool2::LastParameter(*((Curve2*)myC2)); + SubIntervalInitialize(V1, V2); + + switch(Tool1::GetType(*((Curve1*)myC1))) + { + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + case GeomAbs_OtherCurve: + myMaxDerivOrderC1 = MaxOrder; + myTolC1 = SearchOfTolerance((Standard_Address)&C1); + break; + default: + myMaxDerivOrderC1 = 0; + myTolC1=MinTol; + break; + } + + switch(Tool2::GetType(*((Curve2*)myC2))) + { + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + case GeomAbs_OtherCurve: + myMaxDerivOrderC2 = MaxOrder; + myTolC2 = SearchOfTolerance((Standard_Address)&C2); + break; + default: + myMaxDerivOrderC2 = 0; + myTolC2=MinTol; + break; + } + } -Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV, - math_Vector& F) -{ +//============================================================================= +void Extrema_FuncExtCC::SetCurve (const Standard_Integer theRank, const Curve1& C) + { + Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()") + if (theRank == 1) + { + myC1 = (Standard_Address)&C; + switch(/*Tool1::GetType(*((Curve1*)myC1))*/ C.GetType()) + { + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + case GeomAbs_OtherCurve: + myMaxDerivOrderC1 = MaxOrder; + myTolC1 = SearchOfTolerance((Standard_Address)&C); + break; + default: + myMaxDerivOrderC1 = 0; + myTolC1=MinTol; + break; + } + } + else if (theRank == 2) + { + myC2 = (Standard_Address)&C; + switch(/*Tool2::GetType(*((Curve2*)myC2))*/C.GetType()) + { + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + case GeomAbs_OtherCurve: + myMaxDerivOrderC2 = MaxOrder; + myTolC2 = SearchOfTolerance((Standard_Address)&C); + break; + default: + myMaxDerivOrderC2 = 0; + myTolC2=MinTol; + break; + } + } + } + +//============================================================================= +Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV, math_Vector& F) + { myU = UV(1); myV = UV(2); Tool1::D1(*((Curve1*)myC1), myU,myP1,myDu); Tool2::D1(*((Curve2*)myC2), myV,myP2,myDv); + Vec P1P2 (myP1,myP2); Standard_Real Ndu = myDu.Magnitude(); - if (Ndu <= Tol) { - Pnt P1, P2; - P1 = Tool1::Value(*((Curve1*)myC1), myU-delta); - P2 = Tool1::Value(*((Curve1*)myC1), myU+delta); - Vec V(P1,P2); - myDu = V; - Ndu = myDu.Magnitude(); - if (Ndu <= Tol) { - return Standard_False; + + + if(myMaxDerivOrderC1 != 0) + { + if (Ndu <= myTolC1) + { + //Derivative is approximated by Taylor-series + const Standard_Real DivisionFactor = 1.e-3; + Standard_Real du; + if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst())) + du = 0.0; + else + du = myUsupremum-myUinfium; + + const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); + + Standard_Integer n = 1; //Derivative order + Vec V; + Standard_Boolean IsDeriveFound; + + do + { + V = Tool1::DN(*((Curve1*)myC1),myU,++n); + Ndu = V.Magnitude(); + IsDeriveFound = (Ndu > myTolC1); + } + while(!IsDeriveFound && n < myMaxDerivOrderC1); + + if(IsDeriveFound) + { + Standard_Real u; + + if(myU-myUinfium < aDelta) + u = myU+aDelta; + else + u = myU-aDelta; + + Pnt P1, P2; + Tool1::D0(*((Curve1*)myC1),Min(myU, u),P1); + Tool1::D0(*((Curve1*)myC1),Max(myU, u),P2); + + Vec V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + myDu = -V; + else + myDu = V; + }//if(IsDeriveFound) + else + { + //Derivative is approximated by three points + + Pnt Ptemp; //(0,0,0)-coordinate + Pnt P1, P2, P3; + Standard_Boolean IsParameterGrown; + + if(myU-myUinfium < 2*aDelta) + { + Tool1::D0(*((Curve1*)myC1),myU,P1); + Tool1::D0(*((Curve1*)myC1),myU+aDelta,P2); + Tool1::D0(*((Curve1*)myC1),myU+2*aDelta,P3); + IsParameterGrown = Standard_True; + } + else + { + Tool1::D0(*((Curve1*)myC1),myU-2*aDelta,P1); + Tool1::D0(*((Curve1*)myC1),myU-aDelta,P2); + Tool1::D0(*((Curve1*)myC1),myU,P3); + IsParameterGrown = Standard_False; + } + + Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3); + + if(IsParameterGrown) + myDu=-3*V1+4*V2-V3; + else + myDu=V1-4*V2+3*V3; + }//else of if(IsDeriveFound) + Ndu = myDu.Magnitude(); + }//if (Ndu <= myTolC1) condition + }//if(myMaxDerivOrder != 0) + + if (Ndu <= MinTol) + { +#ifdef DEB + cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl; + cout << "Warning: 1st derivative of C1 is equal to zero!---"<= RealLast()) || (myVinfium <= RealFirst())) + dv = 0.0; + else + dv = myVsupremum-myVinfium; + + const Standard_Real aDelta = Max(dv*DivisionFactor,MinStep); + +//Derivative is approximated by Taylor-series + + Standard_Integer n = 1; //Derivative order + Vec V; + Standard_Boolean IsDeriveFound; + + do + { + V = Tool2::DN(*((Curve2*)myC2),myV,++n); + Ndv = V.Magnitude(); + IsDeriveFound = (Ndv > myTolC2); + } + while(!IsDeriveFound && n < myMaxDerivOrderC2); + + if(IsDeriveFound) + { + Standard_Real v; + + if(myV-myVinfium < aDelta) + v = myV+aDelta; + else + v = myV-aDelta; + + Pnt P1, P2; + Tool2::D0(*((Curve2*)myC2),Min(myV, v),P1); + Tool2::D0(*((Curve2*)myC2),Max(myV, v),P2); + + Vec V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + myDv = -V; + else + myDv = V; + }//if(IsDeriveFound) + else + { + //Derivative is approximated by three points + + Pnt Ptemp; //(0,0,0)-coordinate + Pnt P1, P2, P3; + Standard_Boolean IsParameterGrown; + + if(myV-myVinfium < 2*aDelta) + { + Tool2::D0(*((Curve2*)myC2),myV,P1); + Tool2::D0(*((Curve2*)myC2),myV+aDelta,P2); + Tool2::D0(*((Curve2*)myC2),myV+2*aDelta,P3); + IsParameterGrown = Standard_True; + } + else + { + Tool2::D0(*((Curve2*)myC2),myV-2*aDelta,P1); + Tool2::D0(*((Curve2*)myC2),myV-aDelta,P2); + Tool2::D0(*((Curve2*)myC2),myV,P3); + IsParameterGrown = Standard_False; + } + + Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3); + + if(IsParameterGrown) + myDv=-3*V1+4*V2-V3; + else + myDv=V1-4*V2+3*V3; + }//else of if(IsDeriveFound) + + Ndv = myDv.Magnitude(); + }//if (Ndv <= myTolC2) + }//if(myMaxDerivOrder != 0) + + if (Ndv <= MinTol) + { +#ifdef DEB + cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl; + cout << "1st derivative of C2 is equal to zero!---"<= RealLast()) || (myUinfium <= RealFirst())) + du = 0.0; + else + du = myUsupremum-myUinfium; + + const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep); + + Standard_Real dv; + if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst())) + dv = 0.0; + else + dv = myVsupremum-myVinfium; + + const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep); Vec P1P2 (myP1,myP2); - - Standard_Real Ndu = myDu.Magnitude(); - if (Ndu <= Tol) { - Pnt P1, P2; - Vec V1; - Tool1::D1(*((Curve1*)myC1),myU+delta, P2, Duu); - Tool1::D1(*((Curve1*)myC1),myU-delta, P1, V1); - Vec V(P1,P2); - myDu = V; - Duu -= V1; - Ndu = myDu.Magnitude(); - if (Ndu <= Tol) { - return Standard_False; - } - } - - Standard_Real Ndv = myDv.Magnitude(); - if (Ndv <= Tol) { - Pnt P1, P2; - Vec V1; - Tool2::D1(*((Curve2*)myC2),myV+delta, P2, Dvv); - Tool2::D1(*((Curve2*)myC2),myV-delta, P1, V1); - Vec V(P1,P2); - myDv = V; - Dvv -= V1; - Ndv = myDv.Magnitude(); - if (Ndv <= Tol) { - return Standard_False; - } - } - F(1) = P1P2.Dot(myDu)/Ndu; - F(2) = P1P2.Dot(myDv)/Ndv; + if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1)) + { +//Derivative is approximated by three points + + math_Vector FF1(1,2), FF2(1,2), FF3(1,2); + Standard_Real F1, F2, F3; + +/////////////////////////// Search of DF1_u derivative (begin) /////////////////// + if(myU-myUinfium < 2*aDeltaU) + { + F1=F(1); + math_Vector UV2(1,2), UV3(1,2); + UV2(1)=myU+aDeltaU; + UV2(2)=myV; + UV3(1)=myU+2*aDeltaU; + UV3(2)=myV; + if(!((Value(UV2,FF2)) && (Value(UV3,FF3)))) + { +#ifdef DEB + cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; + cout << "There are many points close to singularity points " + "and which have zero-derivative." << endl; + cout << "Try to decrease aDelta variable's value. ---" << endl; +#endif + return Standard_False; + } + + F2 = FF2(1); + F3 = FF3(1); + + Df(1,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU); + }//if(myU-myUinfium < 2*aDeltaU) + else + { + F3 = F(1); + math_Vector UV2(1,2), UV1(1,2); + UV2(1)=myU-aDeltaU; + UV2(2)=myV; + UV1(1)=myU-2*aDeltaU; + UV1(2)=myV; + + if(!((Value(UV2,FF2)) && (Value(UV1,FF1)))) + { +#ifdef DEB + cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; + cout << "There are many points close to singularity points " + "and which have zero-derivative." << endl; + cout << "Try to decrease aDelta variable's value. ---" << endl; +#endif + return Standard_False; + } + + F1 = FF1(1); + F2 = FF2(1); + + Df(1,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU); + }//else of if(myU-myUinfium < 2*aDeltaU) condition +/////////////////////////// Search of DF1_u derivative (end) /////////////////// + + //Return to previous values + myU = myU_old; + myV = myV_old; + +/////////////////////////// Search of DF1_v derivative (begin) /////////////////// + if(myV-myVinfium < 2*aDeltaV) + { + F1=F(1); + math_Vector UV2(1,2), UV3(1,2); + UV2(1)=myU; + UV2(2)=myV+aDeltaV; + UV3(1)=myU; + UV3(2)=myV+2*aDeltaV; + + if(!((Value(UV2,FF2)) && (Value(UV3,FF3)))) + { +#ifdef DEB + cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; + cout << "There are many points close to singularity points " + "and which have zero-derivative." << endl; + cout << "Try to decrease aDelta variable's value. ---" << endl; +#endif + return Standard_False; + } + F2 = FF2(1); + F3 = FF3(1); + + Df(1,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV); + }//if(myV-myVinfium < 2*aDeltaV) + else + { + F3 = F(1); + math_Vector UV2(1,2), UV1(1,2); + UV2(1)=myU; + UV2(2)=myV-aDeltaV; + UV1(1)=myU; + UV1(2)=myV-2*aDeltaV; + if(!((Value(UV2,FF2)) && (Value(UV1,FF1)))) + { +#ifdef DEB + cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; + cout << "There are many points close to singularity points " + "and which have zero-derivative." << endl; + cout << "Try to decrease aDelta variable's value. ---" << endl; +#endif + return Standard_False; + } + + F1 = FF1(1); + F2 = FF2(1); + + Df(1,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV); + }//else of if(myV-myVinfium < 2*aDeltaV) +/////////////////////////// Search of DF1_v derivative (end) /////////////////// + + //Return to previous values + myU = myU_old; + myV = myV_old; + myP1 = myP1_old, myP2 = myP2_old; + myDu = myDu_old, myDv = myDv_old; + }//if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1)) + else + { + const Standard_Real Ndu = myDu.Magnitude(); + Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu)); + Df(1,2) = myDv.Dot(myDu)/Ndu; + }//else of if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1)) + + if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2)) + { +//Derivative is approximated by three points + + math_Vector FF1(1,2), FF2(1,2), FF3(1,2); + Standard_Real F1, F2, F3; + +/////////////////////////// Search of DF2_v derivative (begin) /////////////////// + if(myV-myVinfium < 2*aDeltaV) + { + F1=F(2); + math_Vector UV2(1,2), UV3(1,2); + UV2(1)=myU; + UV2(2)=myV+aDeltaV; + UV3(1)=myU; + UV3(2)=myV+2*aDeltaV; + + if(!((Value(UV2,FF2)) && (Value(UV3,FF3)))) + { +#ifdef DEB + cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; + cout << "There are many points close to singularity points " + "and which have zero-derivative." << endl; + cout << "Try to decrease aDelta variable's value. ---" << endl; +#endif + return Standard_False; + } + + F2 = FF2(2); + F3 = FF3(2); + + Df(2,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV); + + }//if(myV-myVinfium < 2*aDeltaV) + else + { + F3 = F(2); + math_Vector UV2(1,2), UV1(1,2); + UV2(1)=myU; + UV2(2)=myV-aDeltaV; + UV1(1)=myU; + UV1(2)=myV-2*aDeltaV; + + if(!((Value(UV2,FF2)) && (Value(UV1,FF1)))) + { +#ifdef DEB + cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; + cout << "There are many points close to singularity points " + "and which have zero-derivative." << endl; + cout << "Try to decrease aDelta variable's value. ---" << endl; +#endif + return Standard_False; + } + + F1 = FF1(2); + F2 = FF2(2); + + Df(2,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV); + }//else of if(myV-myVinfium < 2*aDeltaV) +/////////////////////////// Search of DF2_v derivative (end) /////////////////// + + //Return to previous values + myU = myU_old; + myV = myV_old; + +/////////////////////////// Search of DF2_u derivative (begin) /////////////////// + if(myU-myUinfium < 2*aDeltaU) + { + F1=F(2); + math_Vector UV2(1,2), UV3(1,2); + UV2(1)=myU+aDeltaU; + UV2(2)=myV; + UV3(1)=myU+2*aDeltaU; + UV3(2)=myV; + if(!((Value(UV2,FF2)) && (Value(UV3,FF3)))) + { +#ifdef DEB + cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; + cout << "There are many points close to singularity points " + "and which have zero-derivative." << endl; + cout << "Try to decrease aDelta variable's value. ---" << endl; +#endif + return Standard_False; + } + + F2 = FF2(2); + F3 = FF3(2); + + Df(2,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU); + }//if(myU-myUinfium < 2*aDelta) + else + { + F3 = F(2); + math_Vector UV2(1,2), UV1(1,2); + UV2(1)=myU-aDeltaU; + UV2(2)=myV; + UV1(1)=myU-2*aDeltaU; + UV1(2)=myV; + + if(!((Value(UV2,FF2)) && (Value(UV1,FF1)))) + { +#ifdef DEB + cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; + cout << "There are many points close to singularity points " + "and which have zero-derivative." << endl; + cout << "Try to decrease aDelta variable's value. ---" << endl; +#endif + return Standard_False; + } + + F1 = FF1(2); + F2 = FF2(2); + + Df(2,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU); + }//else of if(myU-myUinfium < 2*aDeltaU) +/////////////////////////// Search of DF2_u derivative (end) /////////////////// + + //Return to previous values + myU = myU_old; + myV = myV_old; + myP1 = myP1_old; + myP2 = myP2_old; + myDu = myDu_old; + myDv = myDv_old; + }//if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2)) + else + { + Standard_Real Ndv = myDv.Magnitude(); + Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv)); + Df(2,1) = -myDu.Dot(myDv)/Ndv; + }//else of if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2)) - Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu)); - Df(1,2) = myDv.Dot(myDu)/Ndu; - Df(2,1) = -myDu.Dot(myDv)/Ndv; - Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv)); return Standard_True; -} + }//end of function //============================================================================= Standard_Integer Extrema_FuncExtCC::GetStateNumber () @@ -166,11 +730,11 @@ Standard_Integer Extrema_FuncExtCC::GetStateNumber () Vec P1P2 (myP1, myP2); Standard_Real mod = Du.Magnitude(); - if(mod > Tol) { + if(mod > myTolC1) { Du /= mod; } mod = Dv.Magnitude(); - if(mod > Tol) { + if(mod > myTolC2) { Dv /= mod; } @@ -192,9 +756,12 @@ void Extrema_FuncExtCC::Points (const Standard_Integer N, } //============================================================================= - - - - - - +void Extrema_FuncExtCC::SubIntervalInitialize(const math_Vector& theInfBound, + const math_Vector& theSupBound) + { + myUinfium = theInfBound(1); + myUsupremum = theSupBound(1); + myVinfium = theInfBound(2); + myVsupremum = theSupBound(2); + } +//============================================================================= \ No newline at end of file diff --git a/src/Extrema/Extrema_FuncExtCC.lxx b/src/Extrema/Extrema_FuncExtCC.lxx index f66daccdfb..1ff093e5fa 100755 --- a/src/Extrema/Extrema_FuncExtCC.lxx +++ b/src/Extrema/Extrema_FuncExtCC.lxx @@ -16,16 +16,6 @@ // and conditions governing the rights and limitations under the License. //============================================================================= - -inline void Extrema_FuncExtCC::SetCurve (const Standard_Integer theRank, const Curve1& C) -{ - Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()") - if (theRank == 1) {myC1 = (Standard_Address)&C;} - else {myC2 = (Standard_Address)&C;} -} - -//============================================================================= - inline void Extrema_FuncExtCC::SetTolerance (const Standard_Real theTol) { myTol = theTol; diff --git a/src/Extrema/Extrema_FuncExtPC.cdl b/src/Extrema/Extrema_FuncExtPC.cdl index dee78f8310..6a7f9be3d2 100755 --- a/src/Extrema/Extrema_FuncExtPC.cdl +++ b/src/Extrema/Extrema_FuncExtPC.cdl @@ -94,6 +94,14 @@ is TypeMismatch from Standard; -- if N < 1 or N > NbExt(me). + SubIntervalInitialize(me: in out; theUfirst, theUlast: Real from Standard); + ---Purpose: Determines boundaries of subinterval for find of root. + + SearchOfTolerance(me: in out) returns Real from Standard; + ---Purpose: Computes a Tol value. If 1st derivative of curve + -- |D1| #include -#define delta 1.e-9 -#define Tol 1.e-20 +static const Standard_Real TolFactor = 1.e-12; +static const Standard_Real MinTol = 1.e-20; +static const Standard_Real MinStep = 1e-7; + +static const Standard_Integer MaxOrder = 3; + + /*----------------------------------------------------------------------------- Fonction permettant de rechercher une distance extremale entre un point P et @@ -36,6 +41,35 @@ les algorithmes math_FunctionRoot et math_FunctionRoots. = ||D1c(u)|| ** 2 + (C(u)-P).D2c(u) } ----------------------------------------------------------------------------*/ +Standard_Real Extrema_FuncExtPC::SearchOfTolerance() + { + const Standard_Integer NPoint = 10; + const Standard_Real aStep = (myUsupremum - myUinfium)/(Standard_Real)NPoint; + + Standard_Integer aNum = 0; + Standard_Real aMax = -Precision::Infinite(); //Maximum value of 1st derivative + //(it is computed with using NPoint point) + + do + { + Standard_Real u = myUinfium + aNum*aStep; //parameter for every point + if(u > myUsupremum) + u = myUsupremum; + + Pnt Ptemp; //empty point (is not used below) + Vec VDer; // 1st derivative vector + Tool::D1(*((Curve*)myC), u, Ptemp, VDer); + Standard_Real vm = VDer.Magnitude(); + if(vm > aMax) + aMax = vm; + } + while(++aNum < NPoint+1); + + return Max(aMax*TolFactor,MinTol); + + } + +//============================================================================= Extrema_FuncExtPC::Extrema_FuncExtPC(): myU(0.), @@ -44,31 +78,67 @@ myD1f(0.) myPinit = Standard_False; myCinit = Standard_False; myD1Init = Standard_False; + + SubIntervalInitialize(0.0,0.0); + myMaxDerivOrder = 0; + myTol=MinTol; } //============================================================================= Extrema_FuncExtPC::Extrema_FuncExtPC (const Pnt& P, - const Curve& C): -myU(0.), -myD1f(0.) -{ + const Curve& C): myU(0.), myD1f(0.) + { myP = P; myC = (Standard_Address)&C; myPinit = Standard_True; myCinit = Standard_True; myD1Init = Standard_False; -} + + SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)), + Tool::LastParameter(*((Curve*)myC))); + + switch(Tool::GetType(*((Curve*)myC))) + { + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + case GeomAbs_OtherCurve: + myMaxDerivOrder = MaxOrder; + myTol = SearchOfTolerance(); + break; + default: + myMaxDerivOrder = 0; + myTol=MinTol; + break; + } + } //============================================================================= void Extrema_FuncExtPC::Initialize(const Curve& C) -{ + { myC = (Standard_Address)&C; myCinit = Standard_True; myPoint.Clear(); mySqDist.Clear(); myIsMin.Clear(); -} + + SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)), + Tool::LastParameter(*((Curve*)myC))); + + switch(Tool::GetType(*((Curve*)myC))) + { + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + case GeomAbs_OtherCurve: + myMaxDerivOrder = MaxOrder; + myTol = SearchOfTolerance(); + break; + default: + myMaxDerivOrder = 0; + myTol=MinTol; + break; + } + } //============================================================================= @@ -83,30 +153,106 @@ void Extrema_FuncExtPC::SetPoint(const Pnt& P) //============================================================================= + Standard_Boolean Extrema_FuncExtPC::Value (const Standard_Real U, Standard_Real& F) { - if (!myPinit || !myCinit) Standard_TypeMismatch::Raise(); + if (!myPinit || !myCinit) + Standard_TypeMismatch::Raise("No init"); + myU = U; Vec D1c; Tool::D1(*((Curve*)myC),myU,myPc,D1c); Standard_Real Ndu = D1c.Magnitude(); - if (Ndu <= Tol) { // Cas Singulier (PMN 22/04/1998) - Pnt P1, P2; - P2 = Tool::Value(*((Curve*)myC),myU + delta); - P1 = Tool::Value(*((Curve*)myC),myU - delta); - Vec V(P1,P2); - D1c = V; - Ndu = D1c.Magnitude(); - if (Ndu <= Tol) { - Vec aD2; - Tool::D2(*((Curve*)myC),myU,myPc,D1c,aD2); - Ndu = aD2.Magnitude(); - - if(Ndu <= Tol) - return Standard_False; - D1c = aD2; + + if(myMaxDerivOrder != 0) + { + if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998) + { + const Standard_Real DivisionFactor = 1.e-3; + Standard_Real du; + if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst())) + du = 0.0; + else + du = myUsupremum-myUinfium; + + const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); +//Derivative is approximated by Taylor-series + + Standard_Integer n = 1; //Derivative order + Vec V; + Standard_Boolean IsDeriveFound; + + do + { + V = Tool::DN(*((Curve*)myC),myU,++n); + Ndu = V.Magnitude(); + IsDeriveFound = (Ndu > myTol); + } + while(!IsDeriveFound && n < myMaxDerivOrder); + + if(IsDeriveFound) + { + Standard_Real u; + + if(myU-myUinfium < aDelta) + u = myU+aDelta; + else + u = myU-aDelta; + + Pnt P1, P2; + Tool::D0(*((Curve*)myC),Min(myU, u),P1); + Tool::D0(*((Curve*)myC),Max(myU, u),P2); + + Vec V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + D1c = -V; + else + D1c = V; + }//if(IsDeriveFound) + else + { +//Derivative is approximated by three points + + Pnt Ptemp; //(0,0,0)-coordinate + Pnt P1, P2, P3; + Standard_Boolean IsParameterGrown; + + if(myU-myUinfium < 2*aDelta) + { + Tool::D0(*((Curve*)myC),myU,P1); + Tool::D0(*((Curve*)myC),myU+aDelta,P2); + Tool::D0(*((Curve*)myC),myU+2*aDelta,P3); + IsParameterGrown = Standard_True; + } + else + { + Tool::D0(*((Curve*)myC),myU-2*aDelta,P1); + Tool::D0(*((Curve*)myC),myU-aDelta,P2); + Tool::D0(*((Curve*)myC),myU,P3); + IsParameterGrown = Standard_False; + } + + Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3); + + if(IsParameterGrown) + D1c=-3*V1+4*V2-V3; + else + D1c=V1-4*V2+3*V3; + } + Ndu = D1c.Magnitude(); + }//(if (Ndu <= myTol)) condition + }//if(myMaxDerivOrder != 0) + + if (Ndu <= MinTol) + { +#ifdef DEB + cout << "+++Function Extrema_FuncExtPC::Value(...)." << endl; + cout << "Warning: 1st derivative is equal to zero!---"<= RealLast()) || (myUinfium <= RealFirst())) + du = 0.0; + else + du = myUsupremum-myUinfium; + + const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); + + Standard_Real F1, F2, F3; + + if(myU-myUinfium < 2*aDelta) + { + F1=F; + const Standard_Real U1 = myU, U2 = myU+aDelta, U3 = myU+2*aDelta; + + if(!((Value(U2,F2)) && (Value(U3,F3)))) + { +#ifdef DEB + cout << "+++ Function Extrema_FuncExtPC::Values(...)" << endl; + cout << "There are many points close to singularity points " + "and which have zero-derivative." << endl; + cout << "Try to decrease aDelta variable's value. ---" << endl; +#endif + myD1Init = Standard_False; + return Standard_False; + } + +//After calling of Value(...) function variable myU will be redeterminated. +//So we must return it previous value. + D1f=(-3*F1+4*F2-F3)/(2.0*aDelta); + } + else + { + F3 = F; + const Standard_Real U1 = myU-2*aDelta, U2 = myU-aDelta, U3 = myU; + + if(!((Value(U2,F2)) && (Value(U1,F1)))) + { +#ifdef DEB + cout << "+++ Function Extrema_FuncExtPC::Values(...)" << endl; + cout << "There are many points close to singularity points " + "and which have zero-derivative." << endl; + cout << "Try to decrease aDelta variable's value. ---" << endl; +#endif + myD1Init = Standard_False; + return Standard_False; + } +//After calling of Value(...) function variable myU will be redeterminated. +//So we must return it previous value. + D1f=(F1-4*F2+3*F3)/(2.0*aDelta); + } + myU = U; + myPc = myPc_old; + myP = myP_old; + } + else + { + Vec PPc (myP,myPc); + D1f = Ndu + (PPc.Dot(D2c)/Ndu) - F*(D1c.Dot(D2c))/(Ndu*Ndu); } - } - - Vec PPc (myP,myPc); - F = PPc.Dot(D1c)/Ndu; - D1f = Ndu + (PPc.Dot(D2c)/Ndu) - F*(D1c.Dot(D2c))/(Ndu*Ndu); myD1f = D1f; + myD1Init = Standard_True; return Standard_True; -} + } //============================================================================= Standard_Integer Extrema_FuncExtPC::GetStateNumber () @@ -200,3 +414,8 @@ POnC Extrema_FuncExtPC::Point (const Standard_Integer N) const } //============================================================================= +void Extrema_FuncExtPC::SubIntervalInitialize(const Standard_Real theUfirst, const Standard_Real theUlast) + { + myUinfium = theUfirst; + myUsupremum = theUlast; + } \ No newline at end of file diff --git a/src/Extrema/Extrema_GExtCC2d.gxx b/src/Extrema/Extrema_GExtCC2d.gxx index d7e8b84fd4..cf210292cf 100755 --- a/src/Extrema/Extrema_GExtCC2d.gxx +++ b/src/Extrema/Extrema_GExtCC2d.gxx @@ -40,6 +40,20 @@ #include +//======================================================================= +//function : IsParallelDot +//purpose : This function returns True if angle between and +// vectors or between and <-theV2> vectors is less than +//AngTol. +//======================================================================= +static Standard_Boolean IsParallelDot( gp_Vec2d theV1, + gp_Vec2d theV2, + Standard_Real AngTol) + { + return Abs(theV1.Dot(theV2)) >= + theV1.Magnitude()*theV2.Magnitude()*cos(AngTol); + } + Extrema_GExtCC2d::Extrema_GExtCC2d() {} @@ -535,7 +549,8 @@ void Extrema_GExtCC2d::Results(const Extrema_ECC2d& AlgExt, gp_Vec2d v1, v2; Tool1::D1(C1,U,p, v1); Tool2::D1(*((Curve2*)myC),U2,p, v2); - if (v1.IsParallel(v2, Precision::Angular())) { + if (IsParallelDot(v1, v2, Precision::Angular())) + { mynbext++; Val = AlgExt.SquareDistance(i); P1.SetValues(U, P1.Value()); @@ -543,7 +558,7 @@ void Extrema_GExtCC2d::Results(const Extrema_ECC2d& AlgExt, mySqDist.Append(Val); mypoints.Append(P1); mypoints.Append(P2); - } + } // modified by NIZHNY-EAP Thu Jan 27 16:41:00 2000 ___END___ } } diff --git a/src/Extrema/Extrema_GExtPC.gxx b/src/Extrema/Extrema_GExtPC.gxx index ddaf85e83a..54b87d46f2 100755 --- a/src/Extrema/Extrema_GExtPC.gxx +++ b/src/Extrema/Extrema_GExtPC.gxx @@ -131,8 +131,9 @@ void Extrema_GExtPC::Perform(const ThePoint& P) if (myuinf >= anInfToCheck) anInfToCheck = myuinf; if (myusup <= aSupToCheck) aSupToCheck = myusup; if((aSupToCheck - anInfToCheck) <= mytolu) continue; - - if (i != 1) { + + if (i != 1) + { TheCurveTool::D1(*((TheCurve*)myC), myintuinf, PP, V1); s1 = (TheVector(P, PP))*V1; if (s1*s2 < 0.0) { @@ -145,6 +146,7 @@ void Extrema_GExtPC::Perform(const ThePoint& P) TheCurveTool::D1(*((TheCurve*)myC), myintusup, PP, V1); s2 = (TheVector(P, PP))*V1; } + IntervalPerform(P); IntExtIsDone = IntExtIsDone || mydone; } diff --git a/src/Extrema/Extrema_GenExtCC.gxx b/src/Extrema/Extrema_GenExtCC.gxx index cff9b1e548..b24371ce87 100755 --- a/src/Extrema/Extrema_GenExtCC.gxx +++ b/src/Extrema/Extrema_GenExtCC.gxx @@ -184,6 +184,8 @@ b- Calcul des minima: UVinf(2) = aCache2->TrimFirstParameter(); UVsup(1) = aCache1->TrimLastParameter(); UVsup(2) = aCache2->TrimLastParameter(); + + myF.SubIntervalInitialize(UVinf,UVsup); // - des 'bords' du tableau TbDist2 for (NoV = 0; NoV <= aNbV+1; NoV++) { diff --git a/src/Extrema/Extrema_GenExtPC.gxx b/src/Extrema/Extrema_GenExtPC.gxx index 2f354a33a5..2d5a3b929b 100755 --- a/src/Extrema/Extrema_GenExtPC.gxx +++ b/src/Extrema/Extrema_GenExtPC.gxx @@ -171,6 +171,7 @@ Methode: -----------------------------------------------------------------------------*/ { myF.SetPoint(P); + myF.SubIntervalInitialize(myumin,myusup); myDone = Standard_False; math_FunctionRoots S (myF, myumin, myusup, mynbsample, mytolu, mytolF, mytolF); diff --git a/src/Geom/Geom_OffsetCurve.cxx b/src/Geom/Geom_OffsetCurve.cxx index cec8bb11a3..77e3277560 100755 --- a/src/Geom/Geom_OffsetCurve.cxx +++ b/src/Geom/Geom_OffsetCurve.cxx @@ -60,10 +60,11 @@ typedef gp_XYZ XYZ; +//ordre de derivation maximum pour la recherche de la premiere +//derivee non nulle +static const int maxDerivOrder = 3; +static const Standard_Real MinStep = 1e-7; -static const int MaxDegree = 9; - //ordre de derivation maximum pour la recherche de la premiere - //derivee non nulle @@ -297,7 +298,7 @@ void Geom_OffsetCurve::D2 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2) cons //purpose : //======================================================================= -void Geom_OffsetCurve::D3 (const Standard_Real U, Pnt& P, Vec& V1, Vec& V2, Vec& V3) +void Geom_OffsetCurve::D3 (const Standard_Real theU, Pnt& P, Vec& theV1, Vec& V2, Vec& V3) const { @@ -314,22 +315,71 @@ const { // (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir + // (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir + const Standard_Real aTol = gp::Resolution(); + + Standard_Boolean IsDirectionChange = Standard_False; + + basisCurve->D3 (theU, P, theV1, V2, V3); + Vec V4 = basisCurve->DN (theU, 4); + if(theV1.Magnitude() <= aTol) + { + const Standard_Real anUinfium = basisCurve->FirstParameter(); + const Standard_Real anUsupremum = basisCurve->LastParameter(); + + const Standard_Real DivisionFactor = 1.e-3; + Standard_Real du; + if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) + du = 0.0; + else + du = anUsupremum-anUinfium; + + const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); +//Derivative is approximated by Taylor-series + + Standard_Integer anIndex = 1; //Derivative order + Vec V; + + do + { + V = basisCurve->DN(theU,++anIndex); + } + while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder); + + Standard_Real u; + + if(theU-anUinfium < aDelta) + u = theU+aDelta; + else + u = theU-aDelta; + + Pnt P1, P2; + basisCurve->D0(Min(theU, u),P1); + basisCurve->D0(Max(theU, u),P2); + + Vec V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + { + theV1 = -V; + V2 = -basisCurve->DN (theU, anIndex + 1); + V3 = -basisCurve->DN (theU, anIndex + 2); + V4 = -basisCurve->DN (theU, anIndex + 3); + + IsDirectionChange = Standard_True; + } + else + { + theV1 = V; + V2 = basisCurve->DN (theU, anIndex + 1); + V3 = basisCurve->DN (theU, anIndex + 2); + V4 = basisCurve->DN (theU, anIndex + 3); + } + }//if(V1.Magnitude() <= aTol) - basisCurve->D3 (U, P, V1, V2, V3); - Vec V4 = basisCurve->DN (U, 4); - Standard_Integer Index = 2; - while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { - V1 = basisCurve->DN (U, Index); - Index++; - } - if (Index != 2) { - V2 = basisCurve->DN (U, Index); - V3 = basisCurve->DN (U, Index + 1); - V4 = basisCurve->DN (U, Index + 2); - } XYZ OffsetDir = direction.XYZ(); - XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir); + XYZ Ndir = (theV1.XYZ()).Crossed (OffsetDir); XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir); XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir); XYZ D3Ndir = (V4.XYZ()).Crossed (OffsetDir); @@ -351,7 +401,12 @@ const { D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r)); D3Ndir.Multiply (offsetValue/R); + + if(IsDirectionChange) + V3=-V3; + V3.Add (Vec(D3Ndir)); + // V2 = P" (U) : Standard_Real R4 = R2 * R2; D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2)); @@ -362,7 +417,7 @@ const { DNdir.Multiply(R); DNdir.Subtract (Ndir.Multiplied (Dr/R)); DNdir.Multiply (offsetValue/R2); - V1.Add (Vec(DNdir)); + theV1.Add (Vec(DNdir)); } else { // V3 = P"' (U) : @@ -372,7 +427,12 @@ const { D3Ndir.Add (Ndir.Multiplied (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r)); D3Ndir.Multiply (offsetValue); + + if(IsDirectionChange) + V3=-V3; + V3.Add (Vec(D3Ndir)); + // V2 = P" (U) : D2Ndir.Divide (R); D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R3)); @@ -382,12 +442,10 @@ const { // V1 = P' (U) : DNdir.Multiply (offsetValue/R); DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3)); - V1.Add (Vec(DNdir)); + theV1.Add (Vec(DNdir)); } //P (U) : - Ndir.Multiply (offsetValue/R); - Ndir.Add (P.XYZ()); - P.SetXYZ (Ndir); + D0(theU,P); } @@ -396,134 +454,30 @@ const { //purpose : //======================================================================= -Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const { +Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const + { + Standard_RangeError_Raise_if (N < 1, "Exception: " + "Geom_OffsetCurve::DN(...). N<1."); - - if (N < 1) Standard_RangeError::Raise(); - XYZ OffsetDir = direction.XYZ(); - Pnt P; - Vec V1, V2, dummy; - if (N == 1) { - basisCurve->D2 (U, P, V1, V2); - Standard_Integer Index = 2; - while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { - V1 = basisCurve->DN (U, Index); - Index++; - } - if (Index != 2) V2 = basisCurve->DN (U, Index); - XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir); - XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir); - Standard_Real R2 = Ndir.SquareModulus(); - Standard_Real R = Sqrt (R2); - Standard_Real R3 = R * R2; - Standard_Real Dr = Ndir.Dot (DNdir); - if (R3 <= gp::Resolution()) { - if (R2 <= gp::Resolution()) Geom_UndefinedDerivative::Raise(); - Ndir.Multiply (Dr/R); - DNdir.Multiply(R); - DNdir.Subtract (Ndir); - DNdir.Multiply (offsetValue/R2); - V1.Add (Vec(DNdir)); - } - else { - Ndir.Multiply (offsetValue * Dr / R3); - DNdir.Multiply (offsetValue/R); - DNdir.Subtract (Ndir); - V1.Add (Vec(DNdir)); - } - dummy = V1; + gp_Vec VN, Vtemp; + gp_Pnt Ptemp; + switch (N) + { + case 1: + D1( U, Ptemp, VN); + break; + case 2: + D2( U, Ptemp, Vtemp, VN); + break; + case 3: + D3( U, Ptemp, Vtemp, Vtemp, VN); + break; + default: + Standard_NotImplemented::Raise("Exception: " + "Derivative order is greater than 3. Cannot compute of derivative."); } - else if (N == 2) { - Vec V3; - basisCurve->D3 (U, P, V1, V2, V3); - Standard_Integer Index = 2; - while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { - V1 = basisCurve->DN (U, Index); - Index++; - } - if (Index != 2) { - V2 = basisCurve->DN (U, Index); - V3 = basisCurve->DN (U, Index + 1); - } - XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir); - XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir); - XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir); - Standard_Real R2 = Ndir.SquareModulus(); - Standard_Real R = Sqrt (R2); - Standard_Real R3 = R2 * R; Standard_Real R4 = R2 * R2; Standard_Real R5 = R3 * R2; - Standard_Real Dr = Ndir.Dot (DNdir); - Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir); - if (R5 <= gp::Resolution()) { - if (R4 <= gp::Resolution()) Geom_UndefinedDerivative::Raise(); - Ndir.Multiply ((3.0 * Dr * Dr / R4) - (D2r/R2)); - DNdir.Multiply (2.0 * Dr / R2); - D2Ndir.Subtract (DNdir); - D2Ndir.Subtract (Ndir); - D2Ndir.Multiply (offsetValue / R); - V2.Add (Vec(D2Ndir)); - } - else { - Ndir.Multiply ((3.0 * Dr * Dr / R4) - (D2r / R2)); - DNdir.Multiply (2.0 * Dr / R2); - D2Ndir.Divide (R); - D2Ndir.Subtract (DNdir); - D2Ndir.Subtract (Ndir); - D2Ndir.Multiply (offsetValue); - V2.Add (Vec(D2Ndir)); - } - dummy = V2; - } - else if (N == 3) { - Vec V3; - basisCurve->D3 (U, P, V1, V2, V3); - Vec V4 = basisCurve->DN (U, 4); - Standard_Integer Index = 2; - while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { - V1 = basisCurve->DN (U, Index); - Index++; - } - if (Index != 2) { - V2 = basisCurve->DN (U, Index); - V3 = basisCurve->DN (U, Index + 1); - V4 = basisCurve->DN (U, Index + 2); - } - XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir); - XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir); - XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir); - XYZ D3Ndir = (V4.XYZ()).Crossed (OffsetDir); - Standard_Real R2 = Ndir.SquareModulus(); - Standard_Real R = Sqrt (R2); Standard_Real R3 = R2 * R; Standard_Real R4 = R2 * R2; - Standard_Real R5 = R3 * R2; Standard_Real R6 = R3 * R3; Standard_Real R7 = R5 * R2; - Standard_Real Dr = Ndir.Dot (DNdir); - Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir); - Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir); - if (R7 <= gp::Resolution()) { - if (R6 <= gp::Resolution()) Geom_UndefinedDerivative::Raise(); - D2Ndir.Multiply (3.0 * Dr / R2); - DNdir.Multiply (3.0 * ((D2r/R2) + (Dr*Dr)/R4)); - Ndir.Multiply (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r); - D3Ndir.Subtract (D2Ndir); - D3Ndir.Subtract (DNdir); - D3Ndir.Add (Ndir); - D3Ndir.Multiply (offsetValue/R); - V3.Add (Vec(D3Ndir)); - } - else { - D2Ndir.Multiply (3.0 * Dr / R3); - DNdir.Multiplied (3.0 * ((D2r/R3) + (Dr*Dr/R5))); - Ndir.Multiply (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r); - D3Ndir.Divide (R); - D3Ndir.Subtract (D2Ndir); - D3Ndir.Subtract (DNdir); - D3Ndir.Add (Ndir); - D3Ndir.Multiply (offsetValue); - V3.Add (Vec(D3Ndir)); - } - dummy = V3; - } - else { Standard_NotImplemented::Raise(); } - - return dummy; + + return VN; } //======================================================================= @@ -531,24 +485,70 @@ Vec Geom_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const //purpose : //======================================================================= -void Geom_OffsetCurve::D0(const Standard_Real U, - gp_Pnt& P, - gp_Pnt& Pbasis, - gp_Vec& V1basis)const -{ +void Geom_OffsetCurve::D0(const Standard_Real theU, gp_Pnt& theP, + gp_Pnt& thePbasis, gp_Vec& theV1basis)const + { + const Standard_Real aTol = gp::Resolution(); - basisCurve->D1 (U, Pbasis, V1basis); - Standard_Integer Index = 2; - while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { - V1basis = basisCurve->DN (U, Index); - Index++; - } - XYZ Ndir = (V1basis.XYZ()).Crossed (direction.XYZ()); + basisCurve->D1 (theU, thePbasis, theV1basis); + Standard_Real Ndu = theV1basis.Magnitude(); + + if(Ndu <= aTol) + { + const Standard_Real anUinfium = basisCurve->FirstParameter(); + const Standard_Real anUsupremum = basisCurve->LastParameter(); + + const Standard_Real DivisionFactor = 1.e-3; + Standard_Real du; + if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) + du = 0.0; + else + du = anUsupremum-anUinfium; + + const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); +//Derivative is approximated by Taylor-series + + Standard_Integer anIndex = 1; //Derivative order + gp_Vec V; + + do + { + V = basisCurve->DN(theU,++anIndex); + Ndu = V.Magnitude(); + } + while((Ndu <= aTol) && anIndex < maxDerivOrder); + + Standard_Real u; + + if(theU-anUinfium < aDelta) + u = theU+aDelta; + else + u = theU-aDelta; + + gp_Pnt P1, P2; + basisCurve->D0(Min(theU, u),P1); + basisCurve->D0(Max(theU, u),P2); + + gp_Vec V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + theV1basis = -V; + else + theV1basis = V; + + Ndu = theV1basis.Magnitude(); + }//if(Ndu <= aTol) + + XYZ Ndir = (theV1basis.XYZ()).Crossed (direction.XYZ()); Standard_Real R = Ndir.Modulus(); - if (R <= gp::Resolution()) Geom_UndefinedValue::Raise(); + if (R <= gp::Resolution()) + Geom_UndefinedValue::Raise("Exception: Undefined normal vector " + "because tangent vector has zero-magnitude!"); + Ndir.Multiply (offsetValue/R); - Ndir.Add (Pbasis.XYZ()); - P.SetXYZ(Ndir); + Ndir.Add (thePbasis.XYZ()); + theP.SetXYZ(Ndir); } //======================================================================= @@ -556,26 +556,76 @@ void Geom_OffsetCurve::D0(const Standard_Real U, //purpose : //======================================================================= -void Geom_OffsetCurve::D1 ( const Standard_Real U, +void Geom_OffsetCurve::D1 ( const Standard_Real theU, Pnt& P , Pnt& PBasis , - Vec& V1, Vec& V1basis, Vec& V2basis) const { + Vec& theV1, Vec& V1basis, Vec& V2basis) const { // P(u) = p(u) + Offset * Ndir / R // with R = || p' ^ V|| and Ndir = P' ^ direction (local normal direction) // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R)) - basisCurve->D2 (U, PBasis, V1basis, V2basis); - V1 = V1basis; + const Standard_Real aTol = gp::Resolution(); + + basisCurve->D2 (theU, PBasis, V1basis, V2basis); + theV1 = V1basis; Vec V2 = V2basis; - Standard_Integer Index = 2; - while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { - V1 = basisCurve->DN (U, Index); - Index++; - } - if (Index != 2) V2 = basisCurve->DN (U, Index); + + if(theV1.Magnitude() <= aTol) + { + const Standard_Real anUinfium = basisCurve->FirstParameter(); + const Standard_Real anUsupremum = basisCurve->LastParameter(); + + const Standard_Real DivisionFactor = 1.e-3; + Standard_Real du; + if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) + du = 0.0; + else + du = anUsupremum-anUinfium; + + const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); +//Derivative is approximated by Taylor-series + + Standard_Integer anIndex = 1; //Derivative order + Vec V; + + do + { + V = basisCurve->DN(theU,++anIndex); + } + while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder); + + Standard_Real u; + + if(theU-anUinfium < aDelta) + u = theU+aDelta; + else + u = theU-aDelta; + + Pnt P1, P2; + basisCurve->D0(Min(theU, u),P1); + basisCurve->D0(Max(theU, u),P2); + + Vec V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + { + theV1 = -V; + V2 = - basisCurve->DN (theU, anIndex+1); + } + else + { + theV1 = V; + V2 = basisCurve->DN (theU, anIndex+1); + } + + V2basis = V2; + V1basis = theV1; + }//if(theV1.Magnitude() <= aTol) + XYZ OffsetDir = direction.XYZ(); - XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir); + XYZ Ndir = (theV1.XYZ()).Crossed (OffsetDir); XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir); Standard_Real R2 = Ndir.SquareModulus(); Standard_Real R = Sqrt (R2); @@ -587,18 +637,16 @@ void Geom_OffsetCurve::D1 ( const Standard_Real U, DNdir.Multiply(R); DNdir.Subtract (Ndir.Multiplied (Dr/R)); DNdir.Multiply (offsetValue/R2); - V1.Add (Vec(DNdir)); + theV1.Add (Vec(DNdir)); } else { // Same computation as IICURV in EUCLID-IS because the stability is // better DNdir.Multiply (offsetValue/R); DNdir.Subtract (Ndir.Multiplied (offsetValue * Dr/R3)); - V1.Add (Vec(DNdir)); + theV1.Add (Vec(DNdir)); } - Ndir.Multiply (offsetValue/R); - Ndir.Add (PBasis.XYZ()); - P.SetXYZ (Ndir); + D0(theU,P); } @@ -607,9 +655,9 @@ void Geom_OffsetCurve::D1 ( const Standard_Real U, //purpose : //======================================================================= -void Geom_OffsetCurve::D2 (const Standard_Real U, +void Geom_OffsetCurve::D2 (const Standard_Real theU, Pnt& P , Pnt& PBasis , - Vec& V1 , Vec& V2 , + Vec& theV1 , Vec& V2 , Vec& V1basis, Vec& V2basis, Vec& V3basis) const { // P(u) = p(u) + Offset * Ndir / R @@ -620,21 +668,75 @@ void Geom_OffsetCurve::D2 (const Standard_Real U, // P"(u) = p"(u) + (Offset / R) * (D2Ndir/DU - DNdir * (2.0 * Dr/ R**2) + // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2))) - basisCurve->D3 (U, PBasis, V1basis, V2basis, V3basis); - Standard_Integer Index = 2; - V1 = V1basis; + const Standard_Real aTol = gp::Resolution(); + + Standard_Boolean IsDirectionChange = Standard_False; + + basisCurve->D3 (theU, PBasis, V1basis, V2basis, V3basis); + + theV1 = V1basis; V2 = V2basis; Vec V3 = V3basis; - while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { - V1 = basisCurve->DN (U, Index); - Index++; - } - if (Index != 2) { - V2 = basisCurve->DN (U, Index); - V3 = basisCurve->DN (U, Index + 1); - } + + if(theV1.Magnitude() <= aTol) + { + const Standard_Real anUinfium = basisCurve->FirstParameter(); + const Standard_Real anUsupremum = basisCurve->LastParameter(); + + const Standard_Real DivisionFactor = 1.e-3; + Standard_Real du; + if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) + du = 0.0; + else + du = anUsupremum-anUinfium; + + const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); +//Derivative is approximated by Taylor-series + + Standard_Integer anIndex = 1; //Derivative order + Vec V; + + do + { + V = basisCurve->DN(theU,++anIndex); + } + while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder); + + Standard_Real u; + + if(theU-anUinfium < aDelta) + u = theU+aDelta; + else + u = theU-aDelta; + + Pnt P1, P2; + basisCurve->D0(Min(theU, u),P1); + basisCurve->D0(Max(theU, u),P2); + + Vec V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + { + theV1 = -V; + V2 = -basisCurve->DN (theU, anIndex+1); + V3 = -basisCurve->DN (theU, anIndex + 2); + + IsDirectionChange = Standard_True; + } + else + { + theV1 = V; + V2 = basisCurve->DN (theU, anIndex+1); + V3 = basisCurve->DN (theU, anIndex + 2); + } + + V2basis = V2; + V1basis = theV1; + }//if(V1.Magnitude() <= aTol) + XYZ OffsetDir = direction.XYZ(); - XYZ Ndir = (V1.XYZ()).Crossed (OffsetDir); + XYZ Ndir = (theV1.XYZ()).Crossed (OffsetDir); XYZ DNdir = (V2.XYZ()).Crossed (OffsetDir); XYZ D2Ndir = (V3.XYZ()).Crossed (OffsetDir); Standard_Real R2 = Ndir.SquareModulus(); @@ -644,6 +746,7 @@ void Geom_OffsetCurve::D2 (const Standard_Real U, Standard_Real R5 = R3 * R2; Standard_Real Dr = Ndir.Dot (DNdir); Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir); + if (R5 <= gp::Resolution()) { //We try another computation but the stability is not very good //dixit ISG. @@ -653,12 +756,17 @@ void Geom_OffsetCurve::D2 (const Standard_Real U, D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2)); D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2))); D2Ndir.Multiply (offsetValue / R); + + if(IsDirectionChange) + V2=-V2; + V2.Add (Vec(D2Ndir)); + // V1 = P' (U) : DNdir.Multiply(R); DNdir.Subtract (Ndir.Multiplied (Dr/R)); DNdir.Multiply (offsetValue/R2); - V1.Add (Vec(DNdir)); + theV1.Add (Vec(DNdir)); } else { // Same computation as IICURV in EUCLID-IS because the stability is @@ -670,16 +778,19 @@ void Geom_OffsetCurve::D2 (const Standard_Real U, offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3)) ) ); + + if(IsDirectionChange) + V2=-V2; + V2.Add (Vec(D2Ndir)); + // V1 = P' (U) : DNdir.Multiply (offsetValue/R); DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3)); - V1.Add (Vec(DNdir)); + theV1.Add (Vec(DNdir)); } //P (U) : - Ndir.Multiply (offsetValue/R); - Ndir.Add (PBasis.XYZ()); - P.SetXYZ (Ndir); + D0(theU,P); } @@ -717,23 +828,15 @@ Standard_Real Geom_OffsetCurve::Offset () const { return offsetValue; } //purpose : //======================================================================= -void Geom_OffsetCurve::Value ( -const Standard_Real U, Pnt& P, Pnt& PBasis, Vec& V1basis) const { +void Geom_OffsetCurve::Value (const Standard_Real theU, Pnt& theP, + Pnt& thePbasis, Vec& theV1basis) const + { + if (basisCurve->Continuity() == GeomAbs_C0) + Geom_UndefinedValue::Raise("Exception: Basis curve is C0 continuity!"); - if (basisCurve->Continuity() == GeomAbs_C0) Geom_UndefinedValue::Raise(); - basisCurve->D1 (U, PBasis, V1basis); - Standard_Integer Index = 2; - while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { - V1basis = basisCurve->DN (U, Index); - Index++; + basisCurve->D1(theU, thePbasis, theV1basis); + D0(theU,theP); } - XYZ Ndir = (V1basis.XYZ()).Crossed (direction.XYZ()); - Standard_Real R = Ndir.Modulus(); - if (R <= gp::Resolution()) Geom_UndefinedValue::Raise(); - Ndir.Multiply (offsetValue/R); - Ndir.Add (PBasis.XYZ()); - P.SetXYZ (Ndir); -} //======================================================================= diff --git a/src/Geom2d/Geom2d_OffsetCurve.cxx b/src/Geom2d/Geom2d_OffsetCurve.cxx index 583e88cd37..dfd5fc05f0 100755 --- a/src/Geom2d/Geom2d_OffsetCurve.cxx +++ b/src/Geom2d/Geom2d_OffsetCurve.cxx @@ -53,11 +53,10 @@ typedef gp_Trsf2d Trsf2d; typedef gp_XY XY; - -static const int MaxDegree = 9; //ordre de derivation maximum pour la recherche de la premiere //derivee non nulle - +static const int maxDerivOrder = 3; +static const Standard_Real MinStep = 1e-7; //======================================================================= @@ -192,47 +191,141 @@ GeomAbs_Shape Geom2d_OffsetCurve::Continuity () const //purpose : //======================================================================= -void Geom2d_OffsetCurve::D0 (const Standard_Real U, - Pnt2d& P ) const -{ - Vec2d V1; +void Geom2d_OffsetCurve::D0 (const Standard_Real theU, + Pnt2d& theP ) const + { + const Standard_Real aTol = gp::Resolution(); - basisCurve->D1 (U, P, V1); - Standard_Integer Index = 2; - while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { - V1 = basisCurve->DN (U, Index); - Index++; + Vec2d vD1; + + basisCurve->D1 (theU, theP, vD1); + Standard_Real Ndu = vD1.Magnitude(); + + if(Ndu <= aTol) + { + const Standard_Real anUinfium = basisCurve->FirstParameter(); + const Standard_Real anUsupremum = basisCurve->LastParameter(); + + const Standard_Real DivisionFactor = 1.e-3; + Standard_Real du; + if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) + du = 0.0; + else + du = anUsupremum-anUinfium; + + const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); +//Derivative is approximated by Taylor-series + + Standard_Integer anIndex = 1; //Derivative order + Vec2d V; + + do + { + V = basisCurve->DN(theU,++anIndex); + Ndu = V.Magnitude(); + } + while((Ndu <= aTol) && anIndex < maxDerivOrder); + + Standard_Real u; + + if(theU-anUinfium < aDelta) + u = theU+aDelta; + else + u = theU-aDelta; + + Pnt2d P1, P2; + basisCurve->D0(Min(theU, u),P1); + basisCurve->D0(Max(theU, u),P2); + + Vec2d V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + vD1 = -V; + else + vD1 = V; + + Ndu = vD1.Magnitude(); + }//if(Ndu <= aTol) + + if (Ndu <= aTol) + Geom2d_UndefinedValue::Raise("Exception: Undefined normal vector " + "because tangent vector has zero-magnitude!"); + + Standard_Real A = vD1.Y(); + Standard_Real B = - vD1.X(); + A = A * offsetValue/Ndu; + B = B * offsetValue/Ndu; + theP.SetCoord(theP.X() + A, theP.Y() + B); } - Standard_Real A = V1.Y(); - Standard_Real B = - V1.X(); - Standard_Real R = Sqrt(A*A + B * B); - if (R <= gp::Resolution()) Geom2d_UndefinedValue::Raise(); - A = A * offsetValue/R; - B = B * offsetValue/R; - P.SetCoord(P.X() + A, P.Y() + B); -} //======================================================================= //function : D1 //purpose : //======================================================================= - -void Geom2d_OffsetCurve::D1 (const Standard_Real U, Pnt2d& P, Vec2d& V1) const { - +void Geom2d_OffsetCurve::D1 (const Standard_Real theU, Pnt2d& P, Vec2d& theV1) const + { // P(u) = p(u) + Offset * Ndir / R // with R = || p' ^ Z|| and Ndir = P' ^ Z // P'(u) = p'(u) + (Offset / R**2) * (DNdir/DU * R - Ndir * (DR/R)) + const Standard_Real aTol = gp::Resolution(); + Vec2d V2; - basisCurve->D2 (U, P, V1, V2); - Standard_Integer Index = 2; - while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { - V1 = basisCurve->DN (U, Index); - Index++; - } - if (Index != 2) { V2 = basisCurve->DN (U, Index); } - XY Ndir (V1.Y(), -V1.X()); + basisCurve->D2 (theU, P, theV1, V2); + + if(theV1.Magnitude() <= aTol) + { + const Standard_Real anUinfium = basisCurve->FirstParameter(); + const Standard_Real anUsupremum = basisCurve->LastParameter(); + + const Standard_Real DivisionFactor = 1.e-3; + Standard_Real du; + if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) + du = 0.0; + else + du = anUsupremum-anUinfium; + + const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); +//Derivative is approximated by Taylor-series + + Standard_Integer anIndex = 1; //Derivative order + Vec2d V; + + do + { + V = basisCurve->DN(theU,++anIndex); + } + while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder); + + Standard_Real u; + + if(theU-anUinfium < aDelta) + u = theU+aDelta; + else + u = theU-aDelta; + + Pnt2d P1, P2; + basisCurve->D0(Min(theU, u),P1); + basisCurve->D0(Max(theU, u),P2); + + Vec2d V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + { + theV1 = -V; + V2 = - basisCurve->DN (theU, anIndex+1); + } + else + { + theV1 = V; + V2 = basisCurve->DN (theU, anIndex+1); + } + }//if(theV1.Magnitude() <= aTol) + + XY Ndir (theV1.Y(), -theV1.X()); XY DNdir (V2.Y(), -V2.X()); Standard_Real R2 = Ndir.SquareModulus(); Standard_Real R = Sqrt (R2); @@ -244,18 +337,17 @@ void Geom2d_OffsetCurve::D1 (const Standard_Real U, Pnt2d& P, Vec2d& V1) const { DNdir.Multiply(R); DNdir.Subtract (Ndir.Multiplied (Dr/R)); DNdir.Multiply (offsetValue/R2); - V1.Add (Vec2d(DNdir)); + theV1.Add (Vec2d(DNdir)); } else { // Same computation as IICURV in EUCLID-IS because the stability is // better DNdir.Multiply (offsetValue/R); DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3)); - V1.Add (Vec2d(DNdir)); + theV1.Add (Vec2d(DNdir)); } - Ndir.Multiply (offsetValue/R); - Ndir.Add (P.XY()); - P.SetXY (Ndir); + + D0(theU, P); } //======================================================================= @@ -263,9 +355,9 @@ void Geom2d_OffsetCurve::D1 (const Standard_Real U, Pnt2d& P, Vec2d& V1) const { //purpose : //======================================================================= -void Geom2d_OffsetCurve::D2 (const Standard_Real U, +void Geom2d_OffsetCurve::D2 (const Standard_Real theU, Pnt2d& P, - Vec2d& V1, Vec2d& V2) const + Vec2d& theV1, Vec2d& V2) const { // P(u) = p(u) + Offset * Ndir / R // with R = || p' ^ Z|| and Ndir = P' ^ Z @@ -276,17 +368,67 @@ void Geom2d_OffsetCurve::D2 (const Standard_Real U, // Ndir * ( (3.0 * Dr**2 / R**4) - (D2r / R**2))) Vec2d V3; - basisCurve->D3 (U, P, V1, V2, V3); - Standard_Integer Index = 2; - while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { - V1 = basisCurve->DN (U, Index); - Index++; - } - if (Index != 2) { - V2 = basisCurve->DN (U, Index); - V3 = basisCurve->DN (U, Index + 1); - } - XY Ndir (V1.Y(), -V1.X()); + basisCurve->D3 (theU, P, theV1, V2, V3); + + const Standard_Real aTol = gp::Resolution(); + + Standard_Boolean IsDirectionChange = Standard_False; + + if(theV1.Magnitude() <= aTol) + { + const Standard_Real anUinfium = basisCurve->FirstParameter(); + const Standard_Real anUsupremum = basisCurve->LastParameter(); + + const Standard_Real DivisionFactor = 1.e-3; + Standard_Real du; + if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) + du = 0.0; + else + du = anUsupremum-anUinfium; + + const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); +//Derivative is approximated by Taylor-series + + Standard_Integer anIndex = 1; //Derivative order + Vec2d V; + + do + { + V = basisCurve->DN(theU,++anIndex); + } + while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder); + + Standard_Real u; + + if(theU-anUinfium < aDelta) + u = theU+aDelta; + else + u = theU-aDelta; + + Pnt2d P1, P2; + basisCurve->D0(Min(theU, u),P1); + basisCurve->D0(Max(theU, u),P2); + + Vec2d V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + { + theV1 = -V; + V2 = -basisCurve->DN (theU, anIndex+1); + V3 = -basisCurve->DN (theU, anIndex + 2); + + IsDirectionChange = Standard_True; + } + else + { + theV1 = V; + V2 = basisCurve->DN (theU, anIndex+1); + V3 = basisCurve->DN (theU, anIndex + 2); + } + }//if(V1.Magnitude() <= aTol) + + XY Ndir (theV1.Y(), -theV1.X()); XY DNdir (V2.Y(), -V2.X()); XY D2Ndir (V3.Y(), -V3.X()); Standard_Real R2 = Ndir.SquareModulus(); @@ -296,40 +438,54 @@ void Geom2d_OffsetCurve::D2 (const Standard_Real U, Standard_Real R5 = R3 * R2; Standard_Real Dr = Ndir.Dot (DNdir); Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir); - if (R5 <= gp::Resolution()) { + if (R5 <= gp::Resolution()) + { //We try another computation but the stability is not very good //dixit ISG. - if (R4 <= gp::Resolution()) { Geom2d_UndefinedDerivative::Raise(); } + if (R4 <= gp::Resolution()) + { + Geom2d_UndefinedDerivative::Raise(); + } + // V2 = P" (U) : + Standard_Real R4 = R2 * R2; + D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2)); + D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2))); + D2Ndir.Multiply (offsetValue / R); + + if(IsDirectionChange) + V2=-V2; + + V2.Add (Vec2d(D2Ndir)); + + // V1 = P' (U) : + DNdir.Multiply(R); + DNdir.Subtract (Ndir.Multiplied (Dr/R)); + DNdir.Multiply (offsetValue/R2); + theV1.Add (Vec2d(DNdir)); + } + else + { + // Same computation as IICURV in EUCLID-IS because the stability is + // better. // V2 = P" (U) : - Standard_Real R4 = R2 * R2; - D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2)); - D2Ndir.Add (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2))); - D2Ndir.Multiply (offsetValue / R); - V2.Add (Vec2d(D2Ndir)); - // V1 = P' (U) : - DNdir.Multiply(R); - DNdir.Subtract (Ndir.Multiplied (Dr/R)); - DNdir.Multiply (offsetValue/R2); - V1.Add (Vec2d(DNdir)); - } - else { - // Same computation as IICURV in EUCLID-IS because the stability is - // better. - // V2 = P" (U) : D2Ndir.Multiply (offsetValue/R); D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3)); D2Ndir.Add (Ndir.Multiplied (offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3)))); + + if(IsDirectionChange) + V2=-V2; + V2.Add (Vec2d(D2Ndir)); + // V1 = P' (U) - DNdir.Multiply (offsetValue/R); - DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3)); - V1.Add (Vec2d(DNdir)); - } - //P (U) : - Ndir.Multiply (offsetValue/R); - Ndir.Add (P.XY()); - P.SetXY (Ndir); + DNdir.Multiply (offsetValue/R); + DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3)); + theV1.Add (Vec2d(DNdir)); + } + + //P (U) : + D0(theU, P); } @@ -338,9 +494,9 @@ void Geom2d_OffsetCurve::D2 (const Standard_Real U, //purpose : //======================================================================= -void Geom2d_OffsetCurve::D3 (const Standard_Real U, +void Geom2d_OffsetCurve::D3 (const Standard_Real theU, Pnt2d& P, - Vec2d& V1, Vec2d& V2, Vec2d& V3) const { + Vec2d& theV1, Vec2d& V2, Vec2d& V3) const { // P(u) = p(u) + Offset * Ndir / R @@ -356,89 +512,149 @@ void Geom2d_OffsetCurve::D3 (const Standard_Real U, // (D3r/R2) * Ndir + (6.0 * Dr * Dr / R4) * Ndir + // (6.0 * Dr * D2r / R4) * Ndir - (15.0 * Dr* Dr* Dr /R6) * Ndir + const Standard_Real aTol = gp::Resolution(); + + Standard_Boolean IsDirectionChange = Standard_False; + + basisCurve->D3 (theU, P, theV1, V2, V3); + Vec2d V4 = basisCurve->DN (theU, 4); + + if(theV1.Magnitude() <= aTol) + { + const Standard_Real anUinfium = basisCurve->FirstParameter(); + const Standard_Real anUsupremum = basisCurve->LastParameter(); + + const Standard_Real DivisionFactor = 1.e-3; + Standard_Real du; + if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) + du = 0.0; + else + du = anUsupremum-anUinfium; + + const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); +//Derivative is approximated by Taylor-series + + Standard_Integer anIndex = 1; //Derivative order + Vec2d V; + + do + { + V = basisCurve->DN(theU,++anIndex); + } + while((V.Magnitude() <= aTol) && anIndex < maxDerivOrder); + + Standard_Real u; + + if(theU-anUinfium < aDelta) + u = theU+aDelta; + else + u = theU-aDelta; + + Pnt2d P1, P2; + basisCurve->D0(Min(theU, u),P1); + basisCurve->D0(Max(theU, u),P2); + + Vec2d V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + { + theV1 = -V; + V2 = -basisCurve->DN (theU, anIndex + 1); + V3 = -basisCurve->DN (theU, anIndex + 2); + V4 = -basisCurve->DN (theU, anIndex + 3); + + IsDirectionChange = Standard_True; + } + else + { + theV1 = V; + V2 = basisCurve->DN (theU, anIndex + 1); + V3 = basisCurve->DN (theU, anIndex + 2); + V4 = basisCurve->DN (theU, anIndex + 3); + } + }//if(V1.Magnitude() <= aTol) + + XY Ndir (theV1.Y(), -theV1.X()); + XY DNdir (V2.Y(), -V2.X()); + XY D2Ndir (V3.Y(), -V3.X()); + XY D3Ndir (V4.Y(), -V4.X()); + Standard_Real R2 = Ndir.SquareModulus(); + Standard_Real R = Sqrt (R2); + Standard_Real R3 = R2 * R; + Standard_Real R4 = R2 * R2; + Standard_Real R5 = R3 * R2; + Standard_Real R6 = R3 * R3; + Standard_Real R7 = R5 * R2; + Standard_Real Dr = Ndir.Dot (DNdir); + Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir); + Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir); + + if (R7 <= gp::Resolution()) + { + //We try another computation but the stability is not very good + //dixit ISG. + + if (R6 <= gp::Resolution()) + Geom2d_UndefinedDerivative::Raise(); + + // V3 = P"' (U) : + D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R2)); + D3Ndir.Subtract ( + (DNdir.Multiplied ((3.0 * offsetValue) * ((D2r/R2) + (Dr*Dr)/R4)))); + D3Ndir.Add (Ndir.Multiplied ( + (offsetValue * (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r)))); + D3Ndir.Multiply (offsetValue/R); + + if(IsDirectionChange) + V3=-V3; + + V3.Add (Vec2d(D3Ndir)); - basisCurve->D3 (U, P, V1, V2, V3); - Vec2d V4 = basisCurve->DN (U, 4); - Standard_Integer Index = 2; - while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { - V1 = basisCurve->DN (U, Index); - Index++; - } - if (Index != 2) { - V2 = basisCurve->DN (U, Index); - V3 = basisCurve->DN (U, Index + 1); - V4 = basisCurve->DN (U, Index + 2); - } - XY Ndir (V1.Y(), -V1.X()); - XY DNdir (V2.Y(), -V2.X()); - XY D2Ndir (V3.Y(), -V3.X()); - XY D3Ndir (V4.Y(), -V4.X()); - Standard_Real R2 = Ndir.SquareModulus(); - Standard_Real R = Sqrt (R2); - Standard_Real R3 = R2 * R; - Standard_Real R4 = R2 * R2; - Standard_Real R5 = R3 * R2; - Standard_Real R6 = R3 * R3; - Standard_Real R7 = R5 * R2; - Standard_Real Dr = Ndir.Dot (DNdir); - Standard_Real D2r = Ndir.Dot (D2Ndir) + DNdir.Dot (DNdir); - Standard_Real D3r = Ndir.Dot (D3Ndir) + 3.0 * DNdir.Dot (D2Ndir); - if (R7 <= gp::Resolution()) { - //We try another computation but the stability is not very good - //dixit ISG. - if (R6 <= gp::Resolution()) Geom2d_UndefinedDerivative::Raise(); - // V3 = P"' (U) : - D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R2)); - D3Ndir.Subtract ( - (DNdir.Multiplied ((3.0 * offsetValue) * ((D2r/R2) + (Dr*Dr)/R4)))); - D3Ndir.Add (Ndir.Multiplied ( - (offsetValue * (6.0*Dr*Dr/R4 + 6.0*Dr*D2r/R4 - 15.0*Dr*Dr*Dr/R6 - D3r)) - )); - D3Ndir.Multiply (offsetValue/R); - V3.Add (Vec2d(D3Ndir)); - // V2 = P" (U) : - Standard_Real R4 = R2 * R2; - D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2)); - D2Ndir.Subtract (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2))); - D2Ndir.Multiply (offsetValue / R); - V2.Add (Vec2d(D2Ndir)); - // V1 = P' (U) : - DNdir.Multiply(R); - DNdir.Subtract (Ndir.Multiplied (Dr/R)); - DNdir.Multiply (offsetValue/R2); - V1.Add (Vec2d(DNdir)); - } - else { - // Same computation as IICURV in EUCLID-IS because the stability is - // better. - // V3 = P"' (U) : - D3Ndir.Multiply (offsetValue/R); - D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R3)); - D3Ndir.Subtract (DNdir.Multiplied ( - ((3.0 * offsetValue) * ((D2r/R3) + (Dr*Dr)/R5))) ); - D3Ndir.Add (Ndir.Multiplied ( - (offsetValue * (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r)) - )); - V3.Add (Vec2d(D3Ndir)); - // V2 = P" (U) : - D2Ndir.Multiply (offsetValue/R); - D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3)); - D2Ndir.Subtract (Ndir.Multiplied ( - offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3)) - ) - ); - V2.Add (Vec2d(D2Ndir)); - // V1 = P' (U) : - DNdir.Multiply (offsetValue/R); - DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3)); - V1.Add (Vec2d(DNdir)); - } - //P (U) : - Ndir.Multiply (offsetValue/R); - Ndir.Add (P.XY()); - P.SetXY (Ndir); -} + // V2 = P" (U) : + Standard_Real R4 = R2 * R2; + D2Ndir.Subtract (DNdir.Multiplied (2.0 * Dr / R2)); + D2Ndir.Subtract (Ndir.Multiplied (((3.0 * Dr * Dr)/R4) - (D2r/R2))); + D2Ndir.Multiply (offsetValue / R); + V2.Add (Vec2d(D2Ndir)); + // V1 = P' (U) : + DNdir.Multiply(R); + DNdir.Subtract (Ndir.Multiplied (Dr/R)); + DNdir.Multiply (offsetValue/R2); + theV1.Add (Vec2d(DNdir)); + } + else + { + // Same computation as IICURV in EUCLID-IS because the stability is + // better. + // V3 = P"' (U) : + D3Ndir.Multiply (offsetValue/R); + D3Ndir.Subtract (D2Ndir.Multiplied (3.0 * offsetValue * Dr / R3)); + D3Ndir.Subtract (DNdir.Multiplied ( + ((3.0 * offsetValue) * ((D2r/R3) + (Dr*Dr)/R5))) ); + D3Ndir.Add (Ndir.Multiplied ( + (offsetValue * (6.0*Dr*Dr/R5 + 6.0*Dr*D2r/R5 - 15.0*Dr*Dr*Dr/R7 - D3r)))); + + if(IsDirectionChange) + V3=-V3; + + V3.Add (Vec2d(D3Ndir)); + + // V2 = P" (U) : + D2Ndir.Multiply (offsetValue/R); + D2Ndir.Subtract (DNdir.Multiplied (2.0 * offsetValue * Dr / R3)); + D2Ndir.Subtract (Ndir.Multiplied ( + offsetValue * (((3.0 * Dr * Dr) / R5) - (D2r / R3)))); + V2.Add (Vec2d(D2Ndir)); + // V1 = P' (U) : + DNdir.Multiply (offsetValue/R); + DNdir.Subtract (Ndir.Multiplied (offsetValue*Dr/R3)); + theV1.Add (Vec2d(DNdir)); + } + //P (U) : + D0(theU, P); + } //======================================================================= //function : DN @@ -448,7 +664,7 @@ void Geom2d_OffsetCurve::D3 (const Standard_Real U, Vec2d Geom2d_OffsetCurve::DN (const Standard_Real U, const Standard_Integer N) const { - Standard_RangeError_Raise_if (N < 1, "Geom2d_OffsetCurve::DN()"); +Standard_RangeError_Raise_if (N < 1, "Exception: Geom2d_OffsetCurve::DN(). N<1."); gp_Vec2d VN, VBidon; gp_Pnt2d PBidon; @@ -457,7 +673,8 @@ Vec2d Geom2d_OffsetCurve::DN (const Standard_Real U, case 2: D2( U, PBidon, VBidon, VN); break; case 3: D3( U, PBidon, VBidon, VBidon, VN); break; default: - Standard_NotImplemented::Raise(); + Standard_NotImplemented::Raise("Exception: Derivative order is greater than 3. " + "Cannot compute of derivative."); } return VN; @@ -469,25 +686,13 @@ Vec2d Geom2d_OffsetCurve::DN (const Standard_Real U, //purpose : //======================================================================= -void Geom2d_OffsetCurve::Value (const Standard_Real U, - Pnt2d& P, Pnt2d& Pbasis, - Vec2d& V1basis ) const -{ - - basisCurve->D1 (U, Pbasis, V1basis); - Standard_Integer Index = 2; - while (V1basis.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { - V1basis = basisCurve->DN (U, Index); - Index++; +void Geom2d_OffsetCurve::Value (const Standard_Real theU, + Pnt2d& theP, Pnt2d& thePbasis, + Vec2d& theV1basis ) const + { + basisCurve->D1(theU, thePbasis, theV1basis); + D0(theU,theP); } - Standard_Real A = V1basis.Y(); - Standard_Real B = - V1basis.X(); - Standard_Real R = Sqrt(A*A + B * B); - if (R <= gp::Resolution()) Geom2d_UndefinedValue::Raise(); - A = A * offsetValue/R; - B = B * offsetValue/R; - P.SetCoord (A + Pbasis.X(), B + Pbasis.Y()); -} //======================================================================= @@ -509,7 +714,7 @@ void Geom2d_OffsetCurve::D1 (const Standard_Real U, V1 = V1basis; Vec2d V2 = V2basis; Standard_Integer Index = 2; - while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { + while (V1.Magnitude() <= gp::Resolution() && Index <= maxDerivOrder) { V1 = basisCurve->DN (U, Index); Index++; } @@ -567,7 +772,7 @@ void Geom2d_OffsetCurve::D2 (const Standard_Real U, V1 = V1basis; V2 = V2basis; Vec2d V3 = V3basis; - while (V1.Magnitude() <= gp::Resolution() && Index <= MaxDegree) { + while (V1.Magnitude() <= gp::Resolution() && Index <= maxDerivOrder) { V1 = basisCurve->DN (U, Index); Index++; } diff --git a/src/Geom2dInt/Geom2dInt_CurveTool.cdl b/src/Geom2dInt/Geom2dInt_CurveTool.cdl index 9c78ab248c..5b4e8d4ea2 100755 --- a/src/Geom2dInt/Geom2dInt_CurveTool.cdl +++ b/src/Geom2dInt/Geom2dInt_CurveTool.cdl @@ -44,7 +44,7 @@ uses is - TheType(myclass; C: IntCurveCurve) + GetType(myclass; C: IntCurveCurve) ---C++: inline returns CurveType from GeomAbs; @@ -145,6 +145,15 @@ is D2 (myclass; C: IntCurveCurve; U: Real ; P: out Pnt2d; T,N: out Vec2d); ---C++: inline + + D3 (myclass; C: IntCurveCurve; U: Real ; + P: out Pnt2d; T,N,V: out Vec2d); + ---C++: inline + + DN(myclass; C: IntCurveCurve; U: Real ; + N: Integer from Standard) + returns Vec2d; + ---C++: inline NbIntervals(myclass ; C: IntCurveCurve) ---Purpose : output the number of interval of continuity C2 of @@ -165,6 +174,10 @@ is -- used if Type == Composite. ---C++: inline + Degree(myclass; C : IntCurveCurve) returns Integer from Standard; + ---C++: inline + + end CurveTool; diff --git a/src/Geom2dInt/Geom2dInt_CurveTool.lxx b/src/Geom2dInt/Geom2dInt_CurveTool.lxx index 1e16a4046a..79cd30df46 100755 --- a/src/Geom2dInt/Geom2dInt_CurveTool.lxx +++ b/src/Geom2dInt/Geom2dInt_CurveTool.lxx @@ -35,7 +35,7 @@ #define IS_C2_COMPOSITE 0 //============================================================ -inline GeomAbs_CurveType Geom2dInt_CurveTool::TheType(const IntCurveCurve& C) { +inline GeomAbs_CurveType Geom2dInt_CurveTool::GetType(const IntCurveCurve& C) { return(C.GetType()); } //============================================================ @@ -85,6 +85,25 @@ inline void Geom2dInt_CurveTool::D2 (const IntCurveCurve& C, C.D2(U,P,T,N); } + +//============================================================ +inline void Geom2dInt_CurveTool::D3 (const IntCurveCurve& C, + const Standard_Real U, + gp_Pnt2d& P, + gp_Vec2d& T, + gp_Vec2d& N, + gp_Vec2d& V) { + + C.D3(U,P,T,N,V); +} +//============================================================ +inline gp_Vec2d Geom2dInt_CurveTool::DN(const Adaptor2d_Curve2d& C, + const Standard_Real U, + const Standard_Integer N) + { + return C.DN(U,N); + } + //============================================================ inline Standard_Real Geom2dInt_CurveTool::FirstParameter (const IntCurveCurve& C) { return(C.FirstParameter()); @@ -134,4 +153,7 @@ inline Standard_Integer Geom2dInt_CurveTool::NbIntervals(const IntCurveCurve& C) } //============================================================ - + inline Standard_Integer Geom2dInt_CurveTool::Degree(const IntCurveCurve& C) +{ + return C.Degree(); +} diff --git a/src/GeomFill/GeomFill_Frenet.cdl b/src/GeomFill/GeomFill_Frenet.cdl index 2036593a55..31869ee3c1 100755 --- a/src/GeomFill/GeomFill_Frenet.cdl +++ b/src/GeomFill/GeomFill_Frenet.cdl @@ -188,7 +188,19 @@ is -- Warning : It used only for C2 aproximation returns Boolean is private; - + + RotateTrihedron(me; + Tangent : out Vec from gp; + Normal : out Vec from gp; + BiNormal : out Vec from gp; + NewTangent : in Vec from gp) + ---Purpose: revolves the trihedron (which is determined + -- of given "Tangent", "Normal" and "BiNormal" vectors) + -- to coincide "Tangent" and "NewTangent" axes. + returns Boolean from Standard + is private; + + fields P : Pnt from gp; mySngl : HArray1OfReal from TColStd; diff --git a/src/GeomFill/GeomFill_Frenet.cxx b/src/GeomFill/GeomFill_Frenet.cxx index 964f0040d3..8212e121c7 100755 --- a/src/GeomFill/GeomFill_Frenet.cxx +++ b/src/GeomFill/GeomFill_Frenet.cxx @@ -31,8 +31,11 @@ #include #include -#define NullTol 1.e-10 -#define MaxSingular 1.e-5 +static const Standard_Real NullTol = 1.e-10; +static const Standard_Real MaxSingular = 1.e-5; + +static const Standard_Integer maxDerivOrder = 3; + //======================================================================= //function : FDeriv @@ -59,6 +62,30 @@ static gp_Vec DDeriv(const gp_Vec& F, const gp_Vec& DF, const gp_Vec& D2F) return Result; } +//======================================================================= +//function : CosAngle +//purpose : Return a cosine between vectors theV1 and theV2. +//======================================================================= +static Standard_Real CosAngle(const gp_Vec& theV1, const gp_Vec& theV2) + { + const Standard_Real aTol = gp::Resolution(); + + const Standard_Real m1 = theV1.Magnitude(), m2 = theV2.Magnitude(); + if((m1 <= aTol) || (m2 <= aTol)) //Vectors are codirectional + return 1.0; + + Standard_Real aCAng = theV1.Dot(theV2)/(m1*m2); + + if(aCAng > 1.0) + aCAng = 1.0; + + if(aCAng < -1.0) + aCAng = -1.0; + + + return aCAng; + } + //======================================================================= //function : GeomFill_Frenet //purpose : @@ -314,36 +341,205 @@ Handle(GeomFill_TrihedronLaw) GeomFill_Frenet::Copy() const else isSngl = Standard_False; } +//======================================================================= +//function : RotateTrihedron +//purpose : This function revolves the trihedron (which is determined of +// given "Tangent", "Normal" and "BiNormal" vectors) +// to coincide "Tangent" and "NewTangent" axes. +//======================================================================= +Standard_Boolean + GeomFill_Frenet::RotateTrihedron( gp_Vec& Tangent, + gp_Vec& Normal, + gp_Vec& BiNormal, + const gp_Vec& NewTangent) const + { + const Standard_Real anInfCOS = cos(Precision::Angular()); //0.99999995 + const Standard_Real aTol = gp::Resolution(); + + gp_Vec anAxis = Tangent.Crossed(NewTangent); + const Standard_Real NT = anAxis.Magnitude(); + if(NT <= aTol) + //No rotation required + return Standard_True; + else + anAxis /= NT; //Normalization + + const Standard_Real aPx = anAxis.X(), aPy = anAxis.Y(), aPz = anAxis.Z(); + const Standard_Real aCAng = CosAngle(Tangent,NewTangent); //cosine + + const Standard_Real anAddCAng = 1.0 - aCAng; + const Standard_Real aSAng = sqrt(1.0 - aCAng*aCAng); //sine + +//According to rotate direction, sine of rotation angle might be +//positive or negative. +//We can research to choose necessary sign. But I think, it is more +//effectively, to rotate "Tangent" vector in both direction. After that +//we can choose necessary rotation direction in depend of results. + + const gp_Vec aV11(anAddCAng*aPx*aPx+aCAng, anAddCAng*aPx*aPy-aPz*aSAng, anAddCAng*aPx*aPz+aPy*aSAng); + const gp_Vec aV12(anAddCAng*aPx*aPx+aCAng, anAddCAng*aPx*aPy+aPz*aSAng, anAddCAng*aPx*aPz-aPy*aSAng); + const gp_Vec aV21(anAddCAng*aPx*aPy+aPz*aSAng, anAddCAng*aPy*aPy+aCAng, anAddCAng*aPy*aPz-aPx*aSAng); + const gp_Vec aV22(anAddCAng*aPx*aPy-aPz*aSAng, anAddCAng*aPy*aPy+aCAng, anAddCAng*aPy*aPz+aPx*aSAng); + const gp_Vec aV31(anAddCAng*aPx*aPz-aPy*aSAng, anAddCAng*aPy*aPz+aPx*aSAng, anAddCAng*aPz*aPz+aCAng); + const gp_Vec aV32(anAddCAng*aPx*aPz+aPy*aSAng, anAddCAng*aPy*aPz-aPx*aSAng, anAddCAng*aPz*aPz+aCAng); + + gp_Vec aT1(Tangent.Dot(aV11), Tangent.Dot(aV21), Tangent.Dot(aV31)); + gp_Vec aT2(Tangent.Dot(aV12), Tangent.Dot(aV22), Tangent.Dot(aV32)); + + if(CosAngle(aT1,NewTangent) >= CosAngle(aT2,NewTangent)) + { + Tangent = aT1; + Normal = gp_Vec(Normal.Dot(aV11), Normal.Dot(aV21), Normal.Dot(aV31)); + BiNormal = gp_Vec(BiNormal.Dot(aV11), BiNormal.Dot(aV21), BiNormal.Dot(aV31)); + } + else + { + Tangent = aT2; + Normal = gp_Vec(Normal.Dot(aV12), Normal.Dot(aV22), Normal.Dot(aV32)); + BiNormal = gp_Vec(BiNormal.Dot(aV12), BiNormal.Dot(aV22), BiNormal.Dot(aV32)); + } + + return CosAngle(Tangent,NewTangent) >= anInfCOS; + } + + //======================================================================= //function : D0 //purpose : //======================================================================= - Standard_Boolean GeomFill_Frenet::D0(const Standard_Real Param, + Standard_Boolean GeomFill_Frenet::D0(const Standard_Real theParam, gp_Vec& Tangent, gp_Vec& Normal, gp_Vec& BiNormal) { + const Standard_Real aTol = gp::Resolution(); + Standard_Real norm; Standard_Integer Index; Standard_Real Delta = 0.; - if(IsSingular(Param, Index)) - if (SingularD0(Param, Index, Tangent, Normal, BiNormal, Delta)) + if(IsSingular(theParam, Index)) + if (SingularD0(theParam, Index, Tangent, Normal, BiNormal, Delta)) return Standard_True; - Standard_Real theParam = Param + Delta; - myTrimmed->D2(theParam, P, Tangent, BiNormal); - Tangent.Normalize(); - BiNormal = Tangent.Crossed(BiNormal); - norm = BiNormal.Magnitude(); - if (norm <= gp::Resolution()) { - gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent); - BiNormal.SetXYZ(Axe.YDirection().XYZ()); - } - else BiNormal.Normalize(); + Standard_Real aParam = theParam + Delta; + myTrimmed->D2(aParam, P, Tangent, BiNormal); - Normal = BiNormal; - Normal.Cross(Tangent); + const Standard_Real DivisionFactor = 1.e-3; + const Standard_Real anUinfium = myTrimmed->FirstParameter(); + const Standard_Real anUsupremum = myTrimmed->LastParameter(); + const Standard_Real aDelta = (anUsupremum - anUinfium)*DivisionFactor; + Standard_Real Ndu = Tangent.Magnitude(); + +////////////////////////////////////////////////////////////////////////////////////////////////////// + if(Ndu <= aTol) + { + gp_Vec aTn; +//Derivative is approximated by Taylor-series + + Standard_Integer anIndex = 1; //Derivative order + Standard_Boolean isDeriveFound = Standard_False; + + do + { + aTn = myTrimmed->DN(theParam,++anIndex); + Ndu = aTn.Magnitude(); + isDeriveFound = Ndu > aTol; + } + while(!isDeriveFound && anIndex < maxDerivOrder); + + if(isDeriveFound) + { + Standard_Real u; + + if(theParam-anUinfium < aDelta) + u = theParam+aDelta; + else + u = theParam-aDelta; + + gp_Pnt P1, P2; + myTrimmed->D0(Min(theParam, u),P1); + myTrimmed->D0(Max(theParam, u),P2); + + gp_Vec V1(P1,P2); + Standard_Real aDirFactor = aTn.Dot(V1); + + if(aDirFactor < 0.0) + aTn = -aTn; + }//if(IsDeriveFound) + else + { +//Derivative is approximated by three points + + gp_Pnt Ptemp(0.0,0.0,0.0); //(0,0,0)-coordinate + gp_Pnt P1, P2, P3; + Standard_Boolean IsParameterGrown; + + if(theParam-anUinfium < 2*aDelta) + { + myTrimmed->D0(theParam,P1); + myTrimmed->D0(theParam+aDelta,P2); + myTrimmed->D0(theParam+2*aDelta,P3); + IsParameterGrown = Standard_True; + } + else + { + myTrimmed->D0(theParam-2*aDelta,P1); + myTrimmed->D0(theParam-aDelta,P2); + myTrimmed->D0(theParam,P3); + IsParameterGrown = Standard_False; + } + + gp_Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3); + + if(IsParameterGrown) + aTn=-3*V1+4*V2-V3; + else + aTn=V1-4*V2+3*V3; + }//else of "if(IsDeriveFound)" condition + Ndu = aTn.Magnitude(); + gp_Pnt Pt = P; + Standard_Real dPar = 10.0*aDelta; + +//Recursive calling is used for determine of trihedron for +//point, which is near to given. + if(theParam-anUinfium < dPar) + { + if(D0(aParam+dPar,Tangent,Normal,BiNormal) == Standard_False) + return Standard_False; + } + else + { + if(D0(aParam-dPar,Tangent,Normal,BiNormal) == Standard_False) + return Standard_False; + } + + P = Pt; + + if(RotateTrihedron(Tangent,Normal,BiNormal,aTn) == Standard_False) + { +#ifdef DEB + cout << "Cannot coincide two tangents." << endl; +#endif + return Standard_False; + } + }//if(Ndu <= aTol) + else + { + Tangent = Tangent/Ndu; + BiNormal = Tangent.Crossed(BiNormal); + norm = BiNormal.Magnitude(); + if (norm <= gp::Resolution()) + { + gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent); + BiNormal.SetXYZ(Axe.YDirection().XYZ()); + } + else + BiNormal.Normalize(); + + Normal = BiNormal; + Normal.Cross(Tangent); + } return Standard_True; } diff --git a/src/GeomFill/GeomFill_SnglrFunc.cdl b/src/GeomFill/GeomFill_SnglrFunc.cdl index ceb14ac280..67ec8467a5 100755 --- a/src/GeomFill/GeomFill_SnglrFunc.cdl +++ b/src/GeomFill/GeomFill_SnglrFunc.cdl @@ -99,6 +99,27 @@ is --- Purpose : Raised if the continuity of the current interval -- is not C2. is redefined static; + + D3 (me; U : Real; P : out Pnt from gp; V1, V2, V3 : out Vec from gp) + --- Purpose : + -- Returns the point P of parameter U, the first, the second + -- and the third derivative. + raises + DomainError from Standard + --- Purpose : Raised if the continuity of the current interval + -- is not C1. + + is redefined static; + + DN (me; U : Real; N : Integer) returns Vec from gp + --- Purpose : + -- The returned vector gives the value of the derivative for the + -- order of derivation N. + raises + OutOfRange from Standard + --- Purpose : Raised if N < 1. + is redefined static; + Resolution(me; R3d : Real) returns Real ---Purpose : Returns the parametric resolution corresponding diff --git a/src/GeomFill/GeomFill_SnglrFunc.cxx b/src/GeomFill/GeomFill_SnglrFunc.cxx index 859ea9c6bf..de1bedfb18 100755 --- a/src/GeomFill/GeomFill_SnglrFunc.cxx +++ b/src/GeomFill/GeomFill_SnglrFunc.cxx @@ -21,6 +21,7 @@ #include +#include #include GeomFill_SnglrFunc::GeomFill_SnglrFunc(const Handle(Adaptor3d_HCurve)& HC) : @@ -121,6 +122,45 @@ void GeomFill_SnglrFunc::SetRatio(const Standard_Real Ratio) V2 *= ratio; } +void GeomFill_SnglrFunc::D3(const Standard_Real U,gp_Pnt& P,gp_Vec& V1,gp_Vec& V2,gp_Vec& V3) const + { + gp_Vec DC, D2C, D3C, D4C, D5C; + myHCurve->D3(U, P, DC, D2C, D3C); + D4C = myHCurve->DN(U, 4); + D5C = myHCurve->DN(U, 5); + P = gp_Pnt(DC.Crossed(D2C).XYZ()).ChangeCoord()*ratio; + V1 = DC.Crossed(D3C)*ratio; + V2 = (D2C.Crossed(D3C) + DC.Crossed(D4C))*ratio; + V3 = (DC.Crossed(D5C) + D2C.Crossed(D4C)*2)*ratio; + } + +gp_Vec GeomFill_SnglrFunc::DN(const Standard_Real U,const Standard_Integer N) const + { + Standard_RangeError_Raise_if (N < 1, "Exception: Geom2d_OffsetCurve::DN(). N<1."); + + gp_Vec D1C, D2C, D3C; + gp_Pnt C; + + switch(N) + { + case 1: + D1(U,C,D1C); + return D1C; + case 2: + D2(U,C,D1C,D2C); + return D2C; + case 3: + D3(U,C,D1C,D2C,D3C); + return D3C; + default: + Standard_NotImplemented::Raise("Exception: Derivative order is greater than 3. " + "Cannot compute of derivative."); + } + + return gp_Vec(); + + } + Standard_Real GeomFill_SnglrFunc::Resolution(const Standard_Real R3D) const { return Precision::Parametric(R3D); diff --git a/src/GeometryTest/GeometryTest_APICommands.cxx b/src/GeometryTest/GeometryTest_APICommands.cxx index d5fe1ecdfd..6cb85a2fa9 100755 --- a/src/GeometryTest/GeometryTest_APICommands.cxx +++ b/src/GeometryTest/GeometryTest_APICommands.cxx @@ -55,12 +55,12 @@ Standard_IMPORT Draw_Viewer dout; //======================================================================= static Standard_Integer proj (Draw_Interpretor& di, Standard_Integer n, const char** a) -{ - if ( n < 5) { + if ( n < 5) + { cout << " Use proj curve/surf x y z [extrema algo: g(grad)/t(tree)]" << endl; return 1; - } + } gp_Pnt P(Draw::Atof(a[2]),Draw::Atof(a[3]),Draw::Atof(a[4])); @@ -69,70 +69,88 @@ static Standard_Integer proj (Draw_Interpretor& di, Standard_Integer n, const ch Handle(Geom_Curve) GC = DrawTrSurf::GetCurve(a[1]); Handle(Geom_Surface) GS; Extrema_ExtAlgo aProjAlgo = Extrema_ExtAlgo_Grad; + if (n == 6 && a[5][0] == 't') aProjAlgo = Extrema_ExtAlgo_Tree; - if (GC.IsNull()) { + if (GC.IsNull()) + { GS = DrawTrSurf::GetSurface(a[1]); + if (GS.IsNull()) return 1; + Standard_Real U1, U2, V1, V2; GS->Bounds(U1,U2,V1,V2); GeomAPI_ProjectPointOnSurf proj(P,GS,U1,U2,V1,V2,aProjAlgo); + Standard_Real UU,VV; - for ( Standard_Integer i = 1; i <= proj.NbPoints(); i++) { + for ( Standard_Integer i = 1; i <= proj.NbPoints(); i++) + { gp_Pnt P1 = proj.Point(i); - if ( P.Distance(P1) > Precision::Confusion()) { - Handle(Geom_Line) L = new Geom_Line(P,gp_Vec(P,P1)); - Handle(Geom_TrimmedCurve) CT = - new Geom_TrimmedCurve(L, 0., P.Distance(P1)); + if ( P.Distance(P1) > Precision::Confusion()) + { + Handle(Geom_Line) L = new Geom_Line(P,gp_Vec(P,P1)); + Handle(Geom_TrimmedCurve) CT = + new Geom_TrimmedCurve(L, 0., P.Distance(P1)); Sprintf(name,"%s%d","ext_",i); - char* temp = name; // portage WNT - DrawTrSurf::Set(temp, CT); - di << name << " "; - } - else { + char* temp = name; // portage WNT + DrawTrSurf::Set(temp, CT); + di << name << " "; + } + else + { Sprintf(name,"%s%d","ext_",i); - di << name << " "; - char* temp = name; // portage WNT - DrawTrSurf::Set(temp, P1); - proj.Parameters(i,UU,VV); - di << " Le point est sur la surface." << "\n"; - di << " Ses parametres sont: UU = " << UU << "\n"; - di << " VV = " << VV << "\n"; + di << name << " "; + char* temp = name; // portage WNT + DrawTrSurf::Set(temp, P1); + proj.Parameters(i,UU,VV); + di << " Le point est sur la surface." << "\n"; + di << " Ses parametres sont: UU = " << UU << "\n"; + di << " VV = " << VV << "\n"; + } } } - - } - else { + else + { GeomAPI_ProjectPointOnCurve proj(P,GC,GC->FirstParameter(), - GC->LastParameter()); -// Standard_Real UU; - for ( Standard_Integer i = 1; i <= proj.NbPoints(); i++) { + GC->LastParameter()); + + if(proj.NbPoints() == 0) + { + cout << "No project point was found." << endl; + return 0; + } + + for ( Standard_Integer i = 1; i <= proj.NbPoints(); i++) + { gp_Pnt P1 = proj.Point(i); Standard_Real UU = proj.Parameter(i); di << " parameter " << i << " = " << UU << "\n"; - if ( P.Distance(P1) > Precision::Confusion()) { - Handle(Geom_Line) L = new Geom_Line(P,gp_Vec(P,P1)); - Handle(Geom_TrimmedCurve) CT = - new Geom_TrimmedCurve(L, 0., P.Distance(P1)); + if ( P.Distance(P1) > Precision::Confusion()) + { + Handle(Geom_Line) L = new Geom_Line(P,gp_Vec(P,P1)); + Handle(Geom_TrimmedCurve) CT = + new Geom_TrimmedCurve(L, 0., P.Distance(P1)); Sprintf(name,"%s%d","ext_",i); - char* temp = name; // portage WNT - DrawTrSurf::Set(temp, CT); - di << name << " "; - } - else { + char* temp = name; // portage WNT + DrawTrSurf::Set(temp, CT); + di << name << " "; + } + else + { Sprintf(name,"%s%d","ext_",i); - char* temp = name; // portage WNT - DrawTrSurf::Set(temp, P1); - di << name << " "; - UU = proj.Parameter(i); - di << " Le point est sur la courbe." << "\n"; - di << " Son parametre est U = " << UU << "\n"; + char* temp = name; // portage WNT + DrawTrSurf::Set(temp, P1); + di << name << " "; + UU = proj.Parameter(i); + di << " Le point est sur la courbe." << "\n"; + di << " Son parametre est U = " << UU << "\n"; + } } } - } + return 0; } @@ -349,78 +367,133 @@ static Standard_Integer extrema(Draw_Interpretor& di, Standard_Integer n, const } char name[100]; - if ( C1 && C2) { + if ( C1 && C2) + { GeomAPI_ExtremaCurveCurve Ex(GC1,GC2,U1f,U1l,U2f,U2l); - if(!Ex.Extrema().IsParallel()) { - for ( Standard_Integer i = 1; i <= Ex.NbExtrema(); i++) { - gp_Pnt P1,P2; - Ex.Points(i,P1,P2); - if (P1.Distance(P2) < 1.e-16) { - di << "Extrema " << i << " is point : " << P1.X() << " " << P1.Y() << " " << P1.Z() << "\n"; - continue; - } - Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2)); - Handle(Geom_TrimmedCurve) CT = - new Geom_TrimmedCurve(L, 0., P1.Distance(P2)); - Sprintf(name,"%s%d","ext_",i); - char* temp = name; // portage WNT - DrawTrSurf::Set(temp, CT); - di << name << " "; + + if(!Ex.Extrema().IsParallel()) + { + const Standard_Integer aNExtr = Ex.NbExtrema(); + if(aNExtr == 0) + { + di << "No solutions!\n"; + } + else + { + for ( Standard_Integer i = 1; i <= aNExtr; i++) + { + gp_Pnt P1,P2; + Ex.Points(i,P1,P2); + Standard_Real U1,V1; + Ex.Parameters(i,U1,V1); + if (P1.Distance(P2) < 1.e-16) + { + di << "Extrema " << i << " is point : " << P1.X() << " " << P1.Y() << " " << P1.Z() << "\n"; + continue; + } + + Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2)); + Handle(Geom_TrimmedCurve) CT = + new Geom_TrimmedCurve(L, 0., P1.Distance(P2)); +#ifdef DEB + Sprintf(name,"%s%d (U=%f; V=%f)","ext_",i,U1,V1); +#else + Sprintf(name,"%s%d","ext_",i); +#endif + char* temp = name; // portage WNT + DrawTrSurf::Set(temp, CT); + di << name << " "; + } + } + } + else + { + di << "Infinite number of extremas, distance = " << Ex.LowerDistance() << "\n"; } } - else { - di << "Infinite number of extremas, distance = " << Ex.LowerDistance() << "\n"; - } - } - else if ( C1 & S2) { + else if ( C1 & S2) + { GeomAPI_ExtremaCurveSurface Ex(GC1,GS2,U1f,U1l,U2f,U2l,V2f,V2l); - for ( Standard_Integer i = 1; i <= Ex.NbExtrema(); i++) { - gp_Pnt P1,P2; - Ex.Points(i,P1,P2); - if (P1.Distance(P2) < 1.e-16) continue; - Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2)); - Handle(Geom_TrimmedCurve) CT = - new Geom_TrimmedCurve(L, 0., P1.Distance(P2)); - Sprintf(name,"%s%d","ext_",i); - char* temp = name; // portage WNT - DrawTrSurf::Set(temp, CT); - di << name << " "; + + const Standard_Integer aNExtr = Ex.NbExtrema(); + if(aNExtr == 0) + { + di << "No solutions!\n"; + } + else + { + for ( Standard_Integer i = 1; i <= aNExtr; i++) + { + gp_Pnt P1,P2; + Ex.Points(i,P1,P2); + if (P1.Distance(P2) < 1.e-16) continue; + Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2)); + Handle(Geom_TrimmedCurve) CT = + new Geom_TrimmedCurve(L, 0., P1.Distance(P2)); + Sprintf(name,"%s%d","ext_",i); + char* temp = name; // portage WNT + DrawTrSurf::Set(temp, CT); + di << name << " "; + } + } } - } - else if ( S1 & C2) { + else if ( S1 & C2) + { GeomAPI_ExtremaCurveSurface Ex(GC2,GS1,U2f,U2l,U1f,U1l,V1f,V1l); - for ( Standard_Integer i = 1; i <= Ex.NbExtrema(); i++) { - gp_Pnt P1,P2; - Ex.Points(i,P1,P2); - if (P1.Distance(P2) < 1.e-16) continue; - Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2)); - Handle(Geom_TrimmedCurve) CT = - new Geom_TrimmedCurve(L, 0., P1.Distance(P2)); - Sprintf(name,"%s%d","ext_",i); - char* temp = name; // portage WNT - DrawTrSurf::Set(temp, CT); - di << name << " "; + + const Standard_Integer aNExtr = Ex.NbExtrema(); + if(aNExtr == 0) + { + di << "No solutions!\n"; + } + else + { + for ( Standard_Integer i = 1; i <= aNExtr; i++) + { + gp_Pnt P1,P2; + Ex.Points(i,P1,P2); + if (P1.Distance(P2) < 1.e-16) continue; + Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2)); + Handle(Geom_TrimmedCurve) CT = + new Geom_TrimmedCurve(L, 0., P1.Distance(P2)); + Sprintf(name,"%s%d","ext_",i); + char* temp = name; // portage WNT + DrawTrSurf::Set(temp, CT); + di << name << " "; + } + } } - } - else if ( S1 & S2) { + else if ( S1 & S2) + { GeomAPI_ExtremaSurfaceSurface Ex(GS1,GS2,U1f,U1l,V1f,V1l,U2f,U2l,V2f,V2l); - for ( Standard_Integer i = 1; i <= Ex.NbExtrema(); i++) { - gp_Pnt P1,P2; - Ex.Points(i,P1,P2); - if (P1.Distance(P2) < 1.e-16) continue; - Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2)); - Handle(Geom_TrimmedCurve) CT = - new Geom_TrimmedCurve(L, 0., P1.Distance(P2)); - Sprintf(name,"%s%d","ext_",i); - char* temp = name; // portage WNT - DrawTrSurf::Set(temp, CT); - di << name << " "; + + const Standard_Integer aNExtr = Ex.NbExtrema(); + if(aNExtr == 0) + { + di << "No solutions!\n"; + } + else + { + for ( Standard_Integer i = 1; i <= aNExtr; i++) + { + gp_Pnt P1,P2; + Ex.Points(i,P1,P2); + if (P1.Distance(P2) < 1.e-16) + continue; + + Handle(Geom_Line) L = new Geom_Line(P1,gp_Vec(P1,P2)); + Handle(Geom_TrimmedCurve) CT = + new Geom_TrimmedCurve(L, 0., P1.Distance(P2)); + Sprintf(name,"%s%d","ext_",i); + char* temp = name; // portage WNT + DrawTrSurf::Set(temp, CT); + di << name << " "; + } + } } - } return 0; - -} + } //======================================================================= //function : totalextcc diff --git a/src/GeometryTest/GeometryTest_CurveCommands.cxx b/src/GeometryTest/GeometryTest_CurveCommands.cxx index fe4878bdbd..c130e0ff23 100755 --- a/src/GeometryTest/GeometryTest_CurveCommands.cxx +++ b/src/GeometryTest/GeometryTest_CurveCommands.cxx @@ -1329,150 +1329,219 @@ static Standard_Integer surfpoints (Draw_Interpretor& /*di*/, Standard_Integer / //purpose : //======================================================================= static Standard_Integer intersection (Draw_Interpretor& di, Standard_Integer n, const char** a) -{ - if (n < 4) { + { + if (n < 4) return 1; - } + // Handle(Geom_Curve) GC1; Handle(Geom_Surface) GS1 = DrawTrSurf::GetSurface(a[2]); - if (GS1.IsNull()) { + if (GS1.IsNull()) + { GC1 = DrawTrSurf::GetCurve(a[2]); if (GC1.IsNull()) return 1; - } + } + // Handle(Geom_Surface) GS2 = DrawTrSurf::GetSurface(a[3]); - if (GS2.IsNull()) { + if (GS2.IsNull()) return 1; - } + // Standard_Real tol = Precision::Confusion(); - if (n == 5 || n == 9 || n == 13 || n == 17) tol = Draw::Atof(a[n-1]); + if (n == 5 || n == 9 || n == 13 || n == 17) + tol = Draw::Atof(a[n-1]); + // Handle(Geom_Curve) Result; gp_Pnt Point; + // - if (GC1.IsNull()) { + if (GC1.IsNull()) + { GeomInt_IntSS Inters; // // Surface Surface - if (n <= 5) { + if (n <= 5) + { // General case Inters.Perform(GS1,GS2,tol,Standard_True); - } - else if (n == 8 || n == 9 || n == 12 || n == 13 || n == 16 || n == 17) { + } + else if (n == 8 || n == 9 || n == 12 || n == 13 || n == 16 || n == 17) + { Standard_Boolean useStart = Standard_True, useBnd = Standard_True; Standard_Integer ista1=0,ista2=0,ibnd1=0,ibnd2=0; Standard_Real UVsta[4]; Handle(GeomAdaptor_HSurface) AS1,AS2; + // - if (n <= 9) { // user starting point + if (n <= 9) // user starting point + { useBnd = Standard_False; - ista1 = 4; ista2 = 7; - } - else if (n <= 13) { // user bounding + ista1 = 4; + ista2 = 7; + } + else if (n <= 13) // user bounding + { useStart = Standard_False; ibnd1 = 4; ibnd2 = 11; - } - else { // both user starting point and bounding + } + else // both user starting point and bounding + { ista1 = 4; ista2 = 7; ibnd1 = 8; ibnd2 = 15; - } + } + if (useStart) + { for (Standard_Integer i=ista1; i <= ista2; i++) + { UVsta[i-ista1] = Draw::Atof(a[i]); - if (useBnd) { + } + } + + if (useBnd) + { Standard_Real UVbnd[8]; for (Standard_Integer i=ibnd1; i <= ibnd2; i++) UVbnd[i-ibnd1] = Draw::Atof(a[i]); + AS1 = new GeomAdaptor_HSurface(GS1,UVbnd[0],UVbnd[1],UVbnd[2],UVbnd[3]); AS2 = new GeomAdaptor_HSurface(GS2,UVbnd[4],UVbnd[5],UVbnd[6],UVbnd[7]); - } + } + // - if (useStart && !useBnd) { + if (useStart && !useBnd) + { Inters.Perform(GS1,GS2,tol,UVsta[0],UVsta[1],UVsta[2],UVsta[3]); - } - else if (!useStart && useBnd) { + } + else if (!useStart && useBnd) + { Inters.Perform(AS1,AS2,tol); - } - else { + } + else + { Inters.Perform(AS1,AS2,tol,UVsta[0],UVsta[1],UVsta[2],UVsta[3]); - } - }//else if (n == 8 || n == 9 || n == 12 || n == 13 || n == 16 || n == 17) { - else { + } + }//else if (n == 8 || n == 9 || n == 12 || n == 13 || n == 16 || n == 17) + else + { di<<"incorrect number of arguments"<<"\n"; return 1; - } + } + // - if (!Inters.IsDone()) { + if (!Inters.IsDone()) + { + di<<"No intersections found!"<<"\n"; + return 1; - } + } + // char buf[1024]; Standard_Integer i, aNbLines, aNbPoints; - // + + // aNbLines = Inters.NbLines(); - if (aNbLines >= 2) { - for (i=1; i<=aNbLines; ++i) { - Sprintf(buf, "%s_%d",a[1],i); - Result = Inters.Line(i); - const char* temp = buf; - DrawTrSurf::Set(temp,Result); + if (aNbLines >= 2) + { + for (i=1; i<=aNbLines; ++i) + { + Sprintf(buf, "%s_%d",a[1],i); + di << buf << " "; + Result = Inters.Line(i); + const char* temp = buf; + DrawTrSurf::Set(temp,Result); + } } - } - else if (aNbLines == 1) { + else if (aNbLines == 1) + { Result = Inters.Line(1); + Sprintf(buf,"%s",a[1]); + di << buf << " "; DrawTrSurf::Set(a[1],Result); } + // aNbPoints=Inters.NbPoints(); - for (i=1; i<=aNbPoints; ++i) { + for (i=1; i<=aNbPoints; ++i) + { Point=Inters.Point(i); Sprintf(buf,"%s_p_%d",a[1],i); - const char* temp =buf; + di << buf << " "; + const char* temp = buf; DrawTrSurf::Set(temp, Point); - } - }// if (GC1.IsNull()) { - // - else { + } + }// if (GC1.IsNull()) + else + { // Curve Surface GeomAPI_IntCS Inters(GC1,GS2); - // - if (!Inters.IsDone()) return 1; + // + if (!Inters.IsDone()) + { + di<<"No intersections found!"<<"\n"; + return 1; + } + Standard_Integer nblines = Inters.NbSegments(); Standard_Integer nbpoints = Inters.NbPoints(); - if ( (nblines+nbpoints) >= 2) { - char newname[1024]; + + char newname[1024]; + + if ( (nblines+nbpoints) >= 2) + { Standard_Integer i; Standard_Integer Compt = 1; - for (i = 1; i <= nblines; i++, Compt++) { - Sprintf(newname,"%s_%d",a[1],Compt); - Result = Inters.Segment(i); - const char* temp = newname; // pour portage WNT - DrawTrSurf::Set(temp,Result); + + if(nblines >= 1) + cout << " Lines: " << endl; + + for (i = 1; i <= nblines; i++, Compt++) + { + Sprintf(newname,"%s_%d",a[1],Compt); + di << newname << " "; + Result = Inters.Segment(i); + const char* temp = newname; // pour portage WNT + DrawTrSurf::Set(temp,Result); + } + + if(nbpoints >= 1) + cout << " Points: " << endl; + + const Standard_Integer imax = nblines+nbpoints; + + for (/*i = 1*/; i <= imax; i++, Compt++) + { + Sprintf(newname,"%s_%d",a[1],i); + di << newname << " "; + Point = Inters.Point(i); + const char* temp = newname; // pour portage WNT + DrawTrSurf::Set(temp,Point); + } } - for (i = 1; i <= nbpoints; i++, Compt++) { - Sprintf(newname,"%s_%d",a[1],i); - Point = Inters.Point(i); - const char* temp = newname; // pour portage WNT - DrawTrSurf::Set(temp,Point); - } - } - else if (nblines == 1) { + else if (nblines == 1) + { Result = Inters.Segment(1); + Sprintf(newname,"%s",a[1]); + di << newname << " "; DrawTrSurf::Set(a[1],Result); - } - else if (nbpoints == 1) { + } + else if (nbpoints == 1) + { Point = Inters.Point(1); + Sprintf(newname,"%s",a[1]); + di << newname << " "; DrawTrSurf::Set(a[1],Point); + } } - } dout.Flush(); return 0; -} + } //======================================================================= //function : CurveCommands diff --git a/src/GeomliteTest/GeomliteTest_API2dCommands.cxx b/src/GeomliteTest/GeomliteTest_API2dCommands.cxx index f8eeeb69c0..83b40f1d41 100755 --- a/src/GeomliteTest/GeomliteTest_API2dCommands.cxx +++ b/src/GeomliteTest/GeomliteTest_API2dCommands.cxx @@ -245,7 +245,8 @@ static Standard_Integer extrema(Draw_Interpretor& di, Standard_Integer n, const // modified by APV (compilation error - LINUX) // for ( Standard_Integer i = 1; i <= Ex.NbExtrema(); i++) { Standard_Integer i; - for ( i = 1; i <= Ex.NbExtrema(); i++) { + const Standard_Integer aNExtr = Ex.NbExtrema(); + for ( i = 1; i <= aNExtr; i++) { // modified by APV (compilation error - LINUX) gp_Pnt2d P1,P2; @@ -267,7 +268,7 @@ static Standard_Integer extrema(Draw_Interpretor& di, Standard_Integer n, const } } if (i==1) - di << "No decisions "; + di << "No solutions!\n"; return 0; } diff --git a/src/HLRBRep/HLRBRep_CurveTool.cdl b/src/HLRBRep/HLRBRep_CurveTool.cdl index 918d4fbac7..003c3964c3 100755 --- a/src/HLRBRep/HLRBRep_CurveTool.cdl +++ b/src/HLRBRep/HLRBRep_CurveTool.cdl @@ -245,4 +245,8 @@ is NbSamples(myclass; C: Address from Standard) returns Integer from Standard; + Degree(myclass; C: Address from Standard) + returns Integer from Standard; + ---C++: inline + end CurveTool; diff --git a/src/HLRBRep/HLRBRep_CurveTool.lxx b/src/HLRBRep/HLRBRep_CurveTool.lxx index b5ed07f2fe..2065e82f4f 100755 --- a/src/HLRBRep/HLRBRep_CurveTool.lxx +++ b/src/HLRBRep/HLRBRep_CurveTool.lxx @@ -313,3 +313,15 @@ inline Handle(Geom2d_BSplineCurve) inline Standard_Real HLRBRep_CurveTool::EpsX(const Standard_Address C) { return(1e-10); } + + +//======================================================================= +//function : Degree +//purpose : +//======================================================================= + +inline Standard_Integer + HLRBRep_CurveTool::Degree (const Standard_Address C) +{ + return(((HLRBRep_Curve *)C)->Degree()); +} diff --git a/src/IntCurve/IntCurve_IntCurveCurveGen.gxx b/src/IntCurve/IntCurve_IntCurveCurveGen.gxx index 08003db052..9273f42f78 100755 --- a/src/IntCurve/IntCurve_IntCurveCurveGen.gxx +++ b/src/IntCurve/IntCurve_IntCurveCurveGen.gxx @@ -39,7 +39,7 @@ void IntCurve_IntCurveCurveGen::Perform(const TheCurve& C, IntRes2d_Domain D1; Standard_Real TolDomain = Tol; if(Tol #include +static const Standard_Real MinStep = 1.0e-7; -LProp_CLProps::LProp_CLProps (const Curve& C, + +LProp_CLProps::LProp_CLProps (const Curve& C, const Standard_Real U, - const Standard_Integer N, + const Standard_Integer N, const Standard_Real Resolution) - : myCurve(C), - level(N), - cn(4), // cn(Tool::Continuity(C)), RLE - linTol(Resolution), - tangentStatus (LProp_Undecided) + : myCurve(C), myDerOrder(N), myCN(4), + myLinTol(Resolution), myTangentStatus (LProp_Undecided) { - Standard_OutOfRange_Raise_if (N < 0 || N > 3, - "LProp_CLProps::LProp_CLProps()"); - + "LProp_CLProps::LProp_CLProps()"); + SetParameter(U); } -LProp_CLProps::LProp_CLProps (const Curve& C, - const Standard_Integer N, - const Standard_Real Resolution) - : myCurve(C), - u(RealLast()), - level(N), - cn(4), // (Tool::Continuity(C)), RLE - linTol(Resolution), - tangentStatus (LProp_Undecided) - +LProp_CLProps::LProp_CLProps (const Curve& C, const Standard_Integer N, + const Standard_Real Resolution) + : myCurve(C), myU(RealLast()), myDerOrder(N), myCN(4), + myLinTol(Resolution), myTangentStatus (LProp_Undecided) { - - Standard_OutOfRange_Raise_if (N < 0 || N > 3, ""); + Standard_OutOfRange_Raise_if (N < 0 || N > 3, + "LProp_CLProps::LProp_CLProps()"); } LProp_CLProps::LProp_CLProps (const Standard_Integer N, const Standard_Real Resolution) - : u(RealLast()), - level(N), - cn(0), - linTol(Resolution), - tangentStatus (LProp_Undecided) + : myU(RealLast()), myDerOrder(N), myCN(0), myLinTol(Resolution), + myTangentStatus (LProp_Undecided) { - Standard_OutOfRange_Raise_if (N < 0 || N > 3, ""); } void LProp_CLProps::SetParameter(const Standard_Real U) -{ - u = U; - switch (level) { +{ + myU = U; + switch (myDerOrder) + { case 0: - Tool::Value(myCurve, u, pnt); + Tool::Value(myCurve, myU, myPnt); break; case 1: - Tool::D1(myCurve, u, pnt, d[0]); + Tool::D1(myCurve, myU, myPnt, myDerivArr[0]); break; case 2: - Tool::D2(myCurve, u, pnt, d[0], d[1]); + Tool::D2(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1]); break; case 3: - Tool::D3(myCurve, u, pnt, d[0], d[1], d[2]); + Tool::D3(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1], myDerivArr[2]); break; } - tangentStatus = LProp_Undecided; + + myTangentStatus = LProp_Undecided; } -void LProp_CLProps::SetCurve(const Curve& C) { - myCurve = C ; - cn = 4; // Tool::Continuity(C); RLE +void LProp_CLProps::SetCurve(const Curve& C) +{ + myCurve = C ; + myCN = 4; // Tool::Continuity(C); RLE } const Pnt& LProp_CLProps::Value () const -{ - return pnt; +{ + return myPnt; } const Vec& LProp_CLProps::D1 () { - if (level < 1) { - level = 1; - Tool::D1(myCurve, u, pnt, d[0]); + if (myDerOrder < 1) + { + myDerOrder = 1; + Tool::D1(myCurve, myU, myPnt, myDerivArr[0]); } - return d[0]; + + return myDerivArr[0]; } const Vec& LProp_CLProps::D2 () { - if (level < 2) { - level = 2; - Tool::D2(myCurve, u, pnt, d[0], d[1]); + if (myDerOrder < 2) + { + myDerOrder = 2; + Tool::D2(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1]); } - return d[1]; + + return myDerivArr[1]; } const Vec& LProp_CLProps::D3 () { - if (level < 3) { - level = 3; - Tool::D3(myCurve, u, pnt, d[0], d[1], d[2]); + if (myDerOrder < 3) + { + myDerOrder = 3; + Tool::D3(myCurve, myU, myPnt, myDerivArr[0], myDerivArr[1], myDerivArr[2]); } - return d[2]; + + return myDerivArr[2]; } Standard_Boolean LProp_CLProps::IsTangentDefined () { - - if (tangentStatus == LProp_Undefined) { + if (myTangentStatus == LProp_Undefined) return Standard_False; - } - else if (tangentStatus >= LProp_Defined) { + else if (myTangentStatus >= LProp_Defined) return Standard_True; - } // tangentStatus == Lprop_Undecided // we have to calculate the first non null derivative - Standard_Real Tol = linTol * linTol; + const Standard_Real Tol = myLinTol * myLinTol; + Vec V; + Standard_Integer Order = 0; - while (Order < 4) { - Order++; - if(cn >= Order) { - switch(Order) { - case 1 : - V = D1(); - break; - case 2 : - V = D2(); - break; - case 3 : - V = D3(); - break; - }; - if(V.SquareMagnitude() > Tol) { - significantFirstDerivativeOrder = Order; - tangentStatus = LProp_Defined; - return Standard_True; - } - } - else { - tangentStatus = LProp_Undefined; + while (Order++ < 4) + { + if(myCN >= Order) + { + switch(Order) + { + case 1: + V = D1(); + break; + case 2: + V = D2(); + break; + case 3: + V = D3(); + break; + }//switch(Order) + + if(V.SquareMagnitude() > Tol) + { + mySignificantFirstDerivativeOrder = Order; + myTangentStatus = LProp_Defined; + return Standard_True; + }//if(V.SquareMagnitude() > Tol) + }//if(cn >= Order) + else + { + myTangentStatus = LProp_Undefined; return Standard_False; - } - } + }// else of "if(cn >= Order)" condition + }//while (Order < 4) + return Standard_False; } - void LProp_CLProps::Tangent (Dir& D) { + if(!IsTangentDefined()) + LProp_NotDefined::Raise(); + + if(mySignificantFirstDerivativeOrder == 1) + D = Dir(myDerivArr[0]); + else if (mySignificantFirstDerivativeOrder > 1) + { + const Standard_Real DivisionFactor = 1.e-3; + const Standard_Real anUsupremum = Tool::LastParameter(myCurve), + anUinfium = Tool::FirstParameter(myCurve); + + Standard_Real du; + if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) + du = 0.0; + else + du = anUsupremum-anUinfium; + + const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); - if(!IsTangentDefined()) { LProp_NotDefined::Raise(); } - D = Dir(d[significantFirstDerivativeOrder - 1]); + Vec V = myDerivArr[mySignificantFirstDerivativeOrder - 1]; + + Standard_Real u; + + if(myU-anUinfium < aDelta) + u = myU+aDelta; + else + u = myU-aDelta; + + Pnt P1, P2; + Tool::Value(myCurve, Min(myU, u),P1); + Tool::Value(myCurve, Max(myU, u),P2); + + Vec V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + V = -V; + + D = Dir(V); + }//else if (mySignificantFirstDerivativeOrder > 1) } Standard_Real LProp_CLProps::Curvature () { Standard_Boolean isDefined = IsTangentDefined(); LProp_NotDefined_Raise_if(!isDefined, - "LProp_CLProps::CurvatureNotDefined()"); + "LProp_CLProps::CurvatureNotDefined()"); // if the first derivative is null the curvature is infinite. - if(significantFirstDerivativeOrder > 1) return RealLast(); + if(mySignificantFirstDerivativeOrder > 1) + return RealLast(); - Standard_Real Tol = linTol * linTol; - Standard_Real DD1 = d[0].SquareMagnitude(); - Standard_Real DD2 = d[1].SquareMagnitude(); + Standard_Real Tol = myLinTol * myLinTol; + Standard_Real DD1 = myDerivArr[0].SquareMagnitude(); + Standard_Real DD2 = myDerivArr[1].SquareMagnitude(); + // if the second derivative is null the curvature is null. - if (DD2 <= Tol) { - curvature = 0.0; + if (DD2 <= Tol) + { + myCurvature = 0.0; } - else { - Standard_Real N = d[0].CrossSquareMagnitude(d[1]); + else + { + Standard_Real N = myDerivArr[0].CrossSquareMagnitude(myDerivArr[1]); // if d[0] and d[1] are colinear the curvature is null. Standard_Real t = N/(DD1*DD2); - if (t<=Tol) { - curvature = 0.0; + if (t<=Tol) + { + myCurvature = 0.0; } - else { - curvature = sqrt(N) / (DD1*sqrt(DD1)); + else + { + myCurvature = sqrt(N) / (DD1*sqrt(DD1)); } } - return curvature; -} + return myCurvature; +} void LProp_CLProps::Normal (Dir& D) { Standard_Real c = Curvature(); - if(c==RealLast() || Abs(c) <= linTol) { LProp_NotDefined::Raise(); } + if(c==RealLast() || Abs(c) <= myLinTol) + { + LProp_NotDefined::Raise("LProp_CLProps::Normal(...):" + "Curvature is null or infinity"); + } // we used here the following vector relation // a ^ (b ^ c) = b(ac) - c(ab) // Norm = d[0] ^ (d[1] ^ d[0]) - - Vec Norm = d[1] * (d[0] * d[0]) - d[0] * (d[0] * d[1]); + + Vec Norm = myDerivArr[1] * (myDerivArr[0] * myDerivArr[0]) - myDerivArr[0] * (myDerivArr[0] * myDerivArr[1]); D = Dir(Norm); } - void LProp_CLProps::CentreOfCurvature (Pnt& P) { - if(Abs(Curvature()) <= linTol) { LProp_NotDefined::Raise(); } + if(Abs(Curvature()) <= myLinTol) + { + LProp_NotDefined::Raise(); + } // we used here the following vector relation // a ^ (b ^ c) = b(ac) - c(ab) // Norm = d[0] ^ (d[1] ^ d[0]) - Vec Norm = d[1] * (d[0] * d[0]) - d[0] * (d[0] * d[1]); + Vec Norm = myDerivArr[1] * (myDerivArr[0] * myDerivArr[0]) - myDerivArr[0] * (myDerivArr[0] * myDerivArr[1]); Norm.Normalize(); - Norm.Divide(curvature); - P= pnt.Translated(Norm); + Norm.Divide(myCurvature); + P= myPnt.Translated(Norm); } diff --git a/src/LProp/LProp_SLProps.cdl b/src/LProp/LProp_SLProps.cdl index 3e30651074..8ae91900b5 100755 --- a/src/LProp/LProp_SLProps.cdl +++ b/src/LProp/LProp_SLProps.cdl @@ -197,35 +197,35 @@ is fields - surf : Surface; - u : Real; - v : Real; - level : Integer; - cn : Integer; - linTol : Real; + mySurf : Surface; + myU : Real; + myV : Real; + myDerOrder : Integer; + myCN : Integer; + myLinTol : Real; - pnt : Pnt from gp; - d1U : Vec from gp; - d1V : Vec from gp; - d2U : Vec from gp; - d2V : Vec from gp; - dUV : Vec from gp; + myPnt : Pnt from gp; + myD1u : Vec from gp; + myD1v : Vec from gp; + myD2u : Vec from gp; + myD2v : Vec from gp; + myDuv : Vec from gp; - normal : Dir from gp; - minCurv : Real; - maxCurv : Real; - dirMinCurv : Dir from gp; - dirMaxCurv : Dir from gp; - meanCurv : Real; - gausCurv : Real; + myNormal : Dir from gp; + myMinCurv : Real; + myMaxCurv : Real; + myDirMinCurv : Dir from gp; + myDirMaxCurv : Dir from gp; + myMeanCurv : Real; + myGausCurv : Real; - significantFirstUDerivativeOrder : Integer; - significantFirstVDerivativeOrder : Integer; + mySignificantFirstDerivativeOrderU : Integer; + mySignificantFirstDerivativeOrderV : Integer; - uTangentStatus : Status from LProp; - vTangentStatus : Status from LProp; - normalStatus : Status from LProp; - curvatureStatus : Status from LProp; + myUTangentStatus : Status from LProp; + myVTangentStatus : Status from LProp; + myNormalStatus : Status from LProp; + myCurvatureStatus : Status from LProp; end SLProps; diff --git a/src/LProp/LProp_SLProps.gxx b/src/LProp/LProp_SLProps.gxx index 66604b6135..608522e40b 100755 --- a/src/LProp/LProp_SLProps.gxx +++ b/src/LProp/LProp_SLProps.gxx @@ -26,56 +26,62 @@ #include #include +static const Standard_Real MinStep = 1.0e-7; static Standard_Boolean IsTangentDefined (LProp_SLProps& SProp, - const Standard_Integer cn, - const Standard_Real linTol, - const Standard_Integer Derivative, - Standard_Integer& Order, - LProp_Status& Status) + const Standard_Integer cn, + const Standard_Real linTol, + const Standard_Integer Derivative, + Standard_Integer& Order, + LProp_Status& Status) { Standard_Real Tol = linTol * linTol; gp_Vec V[2]; Order = 0; - while (Order < 3) { + + while (Order < 3) + { Order++; - if(cn >= Order) { - switch(Order) { - case 1 : - V[0] = SProp.D1U(); - V[1] = SProp.D1V(); - break; - case 2 : - V[0] = SProp.D2U(); - V[1] = SProp.D2V(); - break; - }; - if(V[Derivative].SquareMagnitude() > Tol) { - Status = LProp_Defined; - return Standard_True; + if(cn >= Order) + { + switch(Order) + { + case 1: + V[0] = SProp.D1U(); + V[1] = SProp.D1V(); + break; + case 2: + V[0] = SProp.D2U(); + V[1] = SProp.D2V(); + break; + }//switch(Order) + + if(V[Derivative].SquareMagnitude() > Tol) + { + Status = LProp_Defined; + return Standard_True; } - } - else { + }//if(cn >= Order) + else + { Status = LProp_Undefined; return Standard_False; } } + return Standard_False; } LProp_SLProps::LProp_SLProps (const Surface& S, - const Standard_Real U, - const Standard_Real V, - const Standard_Integer N, - const Standard_Real Resolution) - : surf(S), - level(N), - cn(4), // (Tool::Continuity(S)), - linTol(Resolution) + const Standard_Real U, + const Standard_Real V, + const Standard_Integer N, + const Standard_Real Resolution) + : mySurf(S),myDerOrder(N), myCN(4), // (Tool::Continuity(S)), + myLinTol(Resolution) { - Standard_OutOfRange_Raise_if(N < 0 || N > 2, - "LProp_SLProps::LProp_SLProps()"); + "LProp_SLProps::LProp_SLProps()"); SetParameters(U, V); } @@ -83,292 +89,328 @@ LProp_SLProps::LProp_SLProps (const Surface& S, LProp_SLProps::LProp_SLProps (const Surface& S, const Standard_Integer N, const Standard_Real Resolution) - : surf(S), - u(RealLast()), v(RealLast()), - level(N), - cn(4), // (Tool::Continuity(S)), - linTol(Resolution), - uTangentStatus (LProp_Undecided), - vTangentStatus (LProp_Undecided), - normalStatus (LProp_Undecided), - curvatureStatus(LProp_Undecided) + : mySurf(S), myU(RealLast()), myV(RealLast()), myDerOrder(N), + myCN(4), // (Tool::Continuity(S)) + myLinTol(Resolution), + myUTangentStatus (LProp_Undecided), + myVTangentStatus (LProp_Undecided), + myNormalStatus (LProp_Undecided), + myCurvatureStatus(LProp_Undecided) { Standard_OutOfRange_Raise_if(N < 0 || N > 2, - "LProp_SLProps::LProp_SLProps()"); + "LProp_SLProps::LProp_SLProps()"); } LProp_SLProps::LProp_SLProps (const Standard_Integer N, - const Standard_Real Resolution) - :u(RealLast()), v(RealLast()), - level(N), - cn(0), - linTol(Resolution), - uTangentStatus (LProp_Undecided), - vTangentStatus (LProp_Undecided), - normalStatus (LProp_Undecided), - curvatureStatus(LProp_Undecided) + const Standard_Real Resolution) + : myU(RealLast()), myV(RealLast()), myDerOrder(N), myCN(0), + myLinTol(Resolution), + myUTangentStatus (LProp_Undecided), + myVTangentStatus (LProp_Undecided), + myNormalStatus (LProp_Undecided), + myCurvatureStatus(LProp_Undecided) { Standard_OutOfRange_Raise_if(N < 0 || N > 2, - "LProp_SLProps::LProp_SLProps() bad level"); + "LProp_SLProps::LProp_SLProps() bad level"); } -void LProp_SLProps::SetSurface (const Surface& S ) { - - surf = S; - cn = 4; // =Tool::Continuity(S); -} - -void LProp_SLProps::SetParameters (const Standard_Real U, const Standard_Real V) +void LProp_SLProps::SetSurface (const Surface& S ) { - u = U; - v = V; - switch (level) { + mySurf = S; + myCN = 4; // =Tool::Continuity(S); +} + +void LProp_SLProps::SetParameters (const Standard_Real U, + const Standard_Real V) +{ + myU = U; + myV = V; + switch (myDerOrder) + { case 0: - Tool::Value(surf, u, v, pnt); + Tool::Value(mySurf, myU, myV, myPnt); break; case 1: - Tool::D1(surf, u, v, pnt, d1U, d1V); + Tool::D1(mySurf, myU, myV, myPnt, myD1u, myD1v); break; case 2: - Tool::D2(surf, u, v, pnt, d1U, d1V, d2U, d2V, dUV); + Tool::D2(mySurf, myU, myV, myPnt, myD1u, myD1v, myD2u, myD2v, myDuv); break; - }; - uTangentStatus = LProp_Undecided; - vTangentStatus = LProp_Undecided; - normalStatus = LProp_Undecided; - curvatureStatus = LProp_Undecided; + } + + myUTangentStatus = LProp_Undecided; + myVTangentStatus = LProp_Undecided; + myNormalStatus = LProp_Undecided; + myCurvatureStatus = LProp_Undecided; } const gp_Pnt& LProp_SLProps::Value() const { - return pnt; + return myPnt; } const gp_Vec& LProp_SLProps::D1U() { - if (level < 1) { - level =1; - Tool::D1(surf,u,v,pnt,d1U,d1V); + if (myDerOrder < 1) + { + myDerOrder =1; + Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v); } - return d1U; + + return myD1u; } const gp_Vec& LProp_SLProps::D1V() { - if (level < 1) { - level =1; - Tool::D1(surf,u,v,pnt,d1U,d1V); + if (myDerOrder < 1) + { + myDerOrder =1; + Tool::D1(mySurf,myU,myV,myPnt,myD1u,myD1v); } - return d1V; + + return myD1v; } const gp_Vec& LProp_SLProps::D2U() { - if (level < 2) { - level =2; - Tool::D2(surf,u,v,pnt,d1U,d1V,d2U,d2V,dUV); + if (myDerOrder < 2) + { + myDerOrder =2; + Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv); } - return d2U; + + return myD2u; } const gp_Vec& LProp_SLProps::D2V() { - if (level < 2) { - level =2; - Tool::D2(surf,u,v,pnt,d1U,d1V,d2U,d2V,dUV); + if (myDerOrder < 2) + { + myDerOrder =2; + Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv); } - return d2V; + + return myD2v; } + const gp_Vec& LProp_SLProps::DUV() { - if (level < 2) { - level =2; - Tool::D2(surf,u,v,pnt,d1U,d1V,d2U,d2V,dUV); + if (myDerOrder < 2) + { + myDerOrder =2; + Tool::D2(mySurf,myU,myV,myPnt,myD1u,myD1v,myD2u,myD2v,myDuv); } - return dUV; + + return myDuv; } Standard_Boolean LProp_SLProps::IsTangentUDefined () { - if (uTangentStatus == LProp_Undefined) { + if (myUTangentStatus == LProp_Undefined) return Standard_False; - } - else if (uTangentStatus >= LProp_Defined) { + else if (myUTangentStatus >= LProp_Defined) return Standard_True; - } + // uTangentStatus == Lprop_Undecided // we have to calculate the first non null U derivative - return IsTangentDefined(*this, cn, linTol, 0, - significantFirstUDerivativeOrder, uTangentStatus); + return IsTangentDefined(*this, myCN, myLinTol, 0, + mySignificantFirstDerivativeOrderU, myUTangentStatus); } -void LProp_SLProps::TangentU (gp_Dir& D) { +void LProp_SLProps::TangentU (gp_Dir& D) +{ + if(!IsTangentUDefined()) + LProp_NotDefined::Raise(); - if(!IsTangentUDefined()) { LProp_NotDefined::Raise(); } - if(significantFirstUDerivativeOrder == 1) { - D = gp_Dir(d1U); - } - else { - D = gp_Dir(d2U); + if(mySignificantFirstDerivativeOrderU == 1) + D = gp_Dir(myD1u); + else + { + const Standard_Real DivisionFactor = 1.e-3; + Standard_Real anUsupremum, anUinfium; + Standard_Real anVsupremum, anVinfium; + Tool::Bounds(mySurf,anUinfium,anVinfium,anUsupremum,anVsupremum); + Standard_Real du; + if((anUsupremum >= RealLast()) || (anUinfium <= RealFirst())) + du = 0.0; + else + du = anUsupremum-anUinfium; + + const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep); + + gp_Vec V = myD2u; + + Standard_Real u; + + if(myU-anUinfium < aDeltaU) + u = myU+aDeltaU; + else + u = myU-aDeltaU; + + gp_Pnt P1, P2; + Tool::Value(mySurf, Min(myU, u),myV,P1); + Tool::Value(mySurf, Max(myU, u),myV,P2); + + gp_Vec V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + V = -V; + + D = gp_Dir(V); } } Standard_Boolean LProp_SLProps::IsTangentVDefined () { - - if (vTangentStatus == LProp_Undefined) { + if (myVTangentStatus == LProp_Undefined) return Standard_False; - } - else if (vTangentStatus >= LProp_Defined) { + else if (myVTangentStatus >= LProp_Defined) return Standard_True; - } + // vTangentStatus == Lprop_Undecided // we have to calculate the first non null V derivative - return IsTangentDefined(*this, cn, linTol, 1, - significantFirstVDerivativeOrder, vTangentStatus); + return IsTangentDefined(*this, myCN, myLinTol, 1, + mySignificantFirstDerivativeOrderV, myVTangentStatus); } -void LProp_SLProps::TangentV (gp_Dir& D) { +void LProp_SLProps::TangentV (gp_Dir& D) +{ + if(!IsTangentVDefined()) + LProp_NotDefined::Raise(); + + if(mySignificantFirstDerivativeOrderV == 1) + D = gp_Dir(myD1v); + else + { + const Standard_Real DivisionFactor = 1.e-3; + Standard_Real anUsupremum, anUinfium; + Standard_Real anVsupremum, anVinfium; + Tool::Bounds(mySurf,anUinfium,anVinfium,anUsupremum,anVsupremum); + + Standard_Real dv; + if((anVsupremum >= RealLast()) || (anVinfium <= RealFirst())) + dv = 0.0; + else + dv = anVsupremum-anVinfium; + + const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep); - if(!IsTangentVDefined()) { LProp_NotDefined::Raise(); } - if(significantFirstVDerivativeOrder == 1) { - D = gp_Dir(d1V); - } - else { - D = gp_Dir(d2V); + gp_Vec V = myD2v; + + Standard_Real v; + + if(myV-anVinfium < aDeltaV) + v = myV+aDeltaV; + else + v = myV-aDeltaV; + + gp_Pnt P1, P2; + Tool::Value(mySurf, myU, Min(myV, v),P1); + Tool::Value(mySurf, myU, Max(myV, v),P2); + + gp_Vec V1(P1,P2); + Standard_Real aDirFactor = V.Dot(V1); + + if(aDirFactor < 0.0) + V = -V; + + D = gp_Dir(V); } } Standard_Boolean LProp_SLProps::IsNormalDefined() { - - if (normalStatus == LProp_Undefined) { + if (myNormalStatus == LProp_Undefined) return Standard_False; - } - else if (normalStatus >= LProp_Defined) { + else if (myNormalStatus >= LProp_Defined) return Standard_True; - } + // status = UnDecided - // first try the standard computation of the normal. CSLib_DerivativeStatus Status; - CSLib::Normal(d1U, d1V, linTol, Status, normal); - if (Status == CSLib_Done ) { - normalStatus = LProp_Computed; + CSLib::Normal(myD1u, myD1v, myLinTol, Status, myNormal); + if (Status == CSLib_Done ) + { + myNormalStatus = LProp_Computed; return Standard_True; } + // else solve the degenerated case only if continuity >= 2 -/* if (cn >= 2) { - if(level < 2) this->D2U(); - Standard_Boolean Done; - CSLib_NormalStatus Stat; - CSLib::Normal(d1U, d1V, d2U, d2V, dUV, linTol, Done, Stat, normal); - if (Done) { - normalStatus = LProp_Computed; - return Standard_True; - } - }*/ - /* - else { - Standard_Integer MaxOrder=3; - CSLib_NormalStatus Stat; - gp_Dir thenormal; - TColgp_Array2OfVec DerNUV(0,MaxOrder,0,MaxOrder); - TColgp_Array2OfVec DerSurf(0,MaxOrder+1,0,MaxOrder+1); - Standard_Integer i,j,OrderU,OrderV; - Standard_Real Umin,Umax,Vmin,Vmax; - - Tool::Bounds(surf, Umin, Vmin, Umax, Vmax); - // Calcul des derivees - for(i=1;i<=MaxOrder+1;i++){ - DerSurf.SetValue(i,0, Tool::DN(surf,u,v,i,0)); - } - - for(i=0;i<=MaxOrder+1;i++) - for(j=1;j<=MaxOrder+1;j++){ - DerSurf.SetValue(i,j, Tool::DN(surf,u,v,i,j)); - } - - for(i=0;i<=MaxOrder;i++) - for(j=0;j<=MaxOrder;j++){ - DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurf)); - } - - CSLib::Normal(MaxOrder,DerNUV, 1.e-9, u,v, - Umin,Umax,Vmin,Vmax, - Stat, thenormal, OrderU, OrderV); - normal.SetXYZ(thenormal.XYZ()); - if (Stat == CSLib_Defined) { - normalStatus = LProp_Computed; - return Standard_True; - } - } - */ - - normalStatus = LProp_Undefined; + myNormalStatus = LProp_Undefined; return Standard_False; } -const gp_Dir& LProp_SLProps::Normal () { +const gp_Dir& LProp_SLProps::Normal () +{ + if(!IsNormalDefined()) + { + LProp_NotDefined::Raise(); + } - if(!IsNormalDefined()) { LProp_NotDefined::Raise(); } - return normal; + return myNormal; } + Standard_Boolean LProp_SLProps::IsCurvatureDefined () { - - if (curvatureStatus == LProp_Undefined) { + if (myCurvatureStatus == LProp_Undefined) return Standard_False; - } - else if (curvatureStatus >= LProp_Defined) { + else if (myCurvatureStatus >= LProp_Defined) return Standard_True; - } - if(cn < 2) { - curvatureStatus = LProp_Undefined; + + if(myCN < 2) + { + myCurvatureStatus = LProp_Undefined; return Standard_False; } + // status = UnDecided - if (!IsNormalDefined()) { - curvatureStatus = LProp_Undefined; + if (!IsNormalDefined()) + { + myCurvatureStatus = LProp_Undefined; return Standard_False; } + // pour eviter un plantage dans le cas du caro pointu // en fait on doit pouvoir calculer les courbure // avoir - if(!IsTangentUDefined() || !IsTangentVDefined()) { - curvatureStatus = LProp_Undefined; + if(!IsTangentUDefined() || !IsTangentVDefined()) + { + myCurvatureStatus = LProp_Undefined; return Standard_False; } // here we compute the curvature features of the surface - gp_Vec Norm (normal); + gp_Vec Norm (myNormal); - Standard_Real E = d1U.SquareMagnitude(); - Standard_Real F = d1U.Dot(d1V); - Standard_Real G = d1V.SquareMagnitude(); + Standard_Real E = myD1u.SquareMagnitude(); + Standard_Real F = myD1u.Dot(myD1v); + Standard_Real G = myD1v.SquareMagnitude(); - if(level < 2) this->D2U(); + if(myDerOrder < 2) + this->D2U(); - Standard_Real L = Norm.Dot(d2U); - Standard_Real M = Norm.Dot(dUV); - Standard_Real N = Norm.Dot(d2V); + Standard_Real L = Norm.Dot(myD2u); + Standard_Real M = Norm.Dot(myDuv); + Standard_Real N = Norm.Dot(myD2v); Standard_Real A = E * M - F * L; Standard_Real B = E * N - G * L; Standard_Real C = F * N - G * M; Standard_Real MaxABC = Max(Max(Abs(A),Abs(B)),Abs(C)); - if (MaxABC < RealEpsilon()) { // ombilic - minCurv = N / G; - maxCurv = minCurv; - dirMinCurv = gp_Dir (d1U); - dirMaxCurv = gp_Dir (d1U.Crossed(Norm)); - meanCurv = minCurv; // (Cmin + Cmax) / 2. - gausCurv = minCurv * minCurv; // (Cmin * Cmax) - curvatureStatus = LProp_Computed; + if (MaxABC < RealEpsilon()) // ombilic + { + myMinCurv = N / G; + myMaxCurv = myMinCurv; + myDirMinCurv = gp_Dir (myD1u); + myDirMaxCurv = gp_Dir (myD1u.Crossed(Norm)); + myMeanCurv = myMinCurv; // (Cmin + Cmax) / 2. + myGausCurv = myMinCurv * myMinCurv; // (Cmin * Cmax) + myCurvatureStatus = LProp_Computed; return Standard_True; } @@ -377,105 +419,124 @@ Standard_Boolean LProp_SLProps::IsCurvatureDefined () C = C / MaxABC; Standard_Real Curv1, Curv2, Root1, Root2; gp_Vec VectCurv1, VectCurv2; - if (Abs(A) > RealEpsilon()) { + + if (Abs(A) > RealEpsilon()) + { math_DirectPolynomialRoots Root (A, B, C); - if(Root.NbSolutions() != 2) { - curvatureStatus = LProp_Undefined; + if(Root.NbSolutions() != 2) + { + myCurvatureStatus = LProp_Undefined; return Standard_False; } - else { - Root1 = Root.Value(1); + else + { + Root1 = Root.Value(1); Root2 = Root.Value(2); Curv1 = ((L * Root1 + 2. * M) * Root1 + N) / ((E * Root1 + 2. * F) * Root1 + G); Curv2 = ((L * Root2 + 2. * M) * Root2 + N) / ((E * Root2 + 2. * F) * Root2 + G); - VectCurv1 = Root1 * d1U + d1V; - VectCurv2 = Root2 * d1U + d1V; + VectCurv1 = Root1 * myD1u + myD1v; + VectCurv2 = Root2 * myD1u + myD1v; } } - else if (Abs(C) > RealEpsilon()) { + else if (Abs(C) > RealEpsilon()) + { math_DirectPolynomialRoots Root(C, B, A); - if((Root.NbSolutions() != 2)) { - curvatureStatus = LProp_Undefined; + + if((Root.NbSolutions() != 2)) + { + myCurvatureStatus = LProp_Undefined; return Standard_False; } - else { + else + { Root1 = Root.Value(1); Root2 = Root.Value(2); Curv1 = ((N * Root1 + 2. * M) * Root1 + L) / ((G * Root1 + 2. * F) * Root1 + E); Curv2 = ((N * Root2 + 2. * M) * Root2 + L) / ((G * Root2 + 2. * F) * Root2 + E); - VectCurv1 = d1U + Root1 * d1V; - VectCurv2 = d1U + Root2 * d1V; + VectCurv1 = myD1u + Root1 * myD1v; + VectCurv2 = myD1u + Root2 * myD1v; } } - else { + else + { Curv1 = L / E; Curv2 = N / G; - VectCurv1 = d1U; - VectCurv2 = d1V; + VectCurv1 = myD1u; + VectCurv2 = myD1v; } - if (Curv1 < Curv2) { - minCurv = Curv1; - maxCurv = Curv2; - dirMinCurv = gp_Dir (VectCurv1); - dirMaxCurv = gp_Dir (VectCurv2); - } - else { - minCurv = Curv2; - maxCurv = Curv1; - dirMinCurv = gp_Dir (VectCurv2); - dirMaxCurv = gp_Dir (VectCurv1); - } - meanCurv = ((N * E) - (2. * M * F) + (L * G)) // voir Farin p.282 - / (2. * ((E * G) - (F * F))); - gausCurv = ((L * N) - (M * M)) - / ((E * G) - (F * F)); - curvatureStatus = LProp_Computed; - return Standard_True; -} + if (Curv1 < Curv2) + { + myMinCurv = Curv1; + myMaxCurv = Curv2; + myDirMinCurv = gp_Dir (VectCurv1); + myDirMaxCurv = gp_Dir (VectCurv2); + } + else + { + myMinCurv = Curv2; + myMaxCurv = Curv1; + myDirMinCurv = gp_Dir (VectCurv2); + myDirMaxCurv = gp_Dir (VectCurv1); + } + + myMeanCurv = ((N * E) - (2. * M * F) + (L * G)) // voir Farin p.282 + / (2. * ((E * G) - (F * F))); + myGausCurv = ((L * N) - (M * M)) + / ((E * G) - (F * F)); + myCurvatureStatus = LProp_Computed; + return Standard_True; +} Standard_Boolean LProp_SLProps::IsUmbilic () { - if(!IsCurvatureDefined()) { LProp_NotDefined::Raise(); } - return Abs(maxCurv - minCurv) < Abs(Epsilon(maxCurv)); + if(!IsCurvatureDefined()) + LProp_NotDefined::Raise(); + + return Abs(myMaxCurv - myMinCurv) < Abs(Epsilon(myMaxCurv)); } Standard_Real LProp_SLProps::MaxCurvature () { - if(!IsCurvatureDefined()) { LProp_NotDefined::Raise(); } - return maxCurv; + if(!IsCurvatureDefined()) + LProp_NotDefined::Raise(); + + return myMaxCurv; } Standard_Real LProp_SLProps::MinCurvature () { - if(!IsCurvatureDefined()) { LProp_NotDefined::Raise(); } - return minCurv; + if(!IsCurvatureDefined()) + LProp_NotDefined::Raise(); + + return myMinCurv; } void LProp_SLProps::CurvatureDirections(gp_Dir& Max, gp_Dir& Min) { + if(!IsCurvatureDefined()) + LProp_NotDefined::Raise(); - if(!IsCurvatureDefined()) { LProp_NotDefined::Raise(); } - Max = dirMaxCurv; - Min = dirMinCurv; + Max = myDirMaxCurv; + Min = myDirMinCurv; } -Standard_Real LProp_SLProps::MeanCurvature () { - - if(!IsCurvatureDefined()) { LProp_NotDefined::Raise(); } - return meanCurv; +Standard_Real LProp_SLProps::MeanCurvature () +{ + if(!IsCurvatureDefined()) + LProp_NotDefined::Raise(); + return myMeanCurv; } -Standard_Real LProp_SLProps::GaussianCurvature () { +Standard_Real LProp_SLProps::GaussianCurvature () +{ + if(!IsCurvatureDefined()) + LProp_NotDefined::Raise(); - if(!IsCurvatureDefined()) { LProp_NotDefined::Raise(); } - return gausCurv; + return myGausCurv; } - - - diff --git a/src/gp/gp_Vec.cdl b/src/gp/gp_Vec.cdl index 6bebd43b53..291d592278 100755 --- a/src/gp/gp_Vec.cdl +++ b/src/gp/gp_Vec.cdl @@ -151,7 +151,7 @@ is raises VectorWithNullMagnitude is static; - + IsParallel (me; Other : Vec; AngularTolerance : Real) returns Boolean --- Purpose : -- Returns True if Angle(, Other) <= AngularTolerance or @@ -164,7 +164,7 @@ is raises VectorWithNullMagnitude is static; - + Angle (me; Other : Vec) returns Real --- Purpose : -- Computes the angular value between and diff --git a/src/gp/gp_Vec2d.cdl b/src/gp/gp_Vec2d.cdl index d77d101d3a..6d2dd32d5a 100755 --- a/src/gp/gp_Vec2d.cdl +++ b/src/gp/gp_Vec2d.cdl @@ -136,7 +136,7 @@ is raises VectorWithNullMagnitude is static; - + IsParallel (me; Other : Vec2d; AngularTolerance : Real) returns Boolean ---C++: inline @@ -212,6 +212,11 @@ is --- Purpose : Computes the scalar product ---C++: alias operator * + GetNormal(me) returns Vec2d is static; + ---C++: inline + -- Purpose : Returns a vector {Y(), -X()} which + -- is normal to given. + Multiply (me : in out; Scalar : Real) is static; ---C++: inline ---C++: alias operator *= diff --git a/src/gp/gp_Vec2d.lxx b/src/gp/gp_Vec2d.lxx index fa6e3e4ca4..8585c6fe65 100755 --- a/src/gp/gp_Vec2d.lxx +++ b/src/gp/gp_Vec2d.lxx @@ -244,3 +244,8 @@ inline gp_Vec2d operator* (const Standard_Real Scalar, const gp_Vec2d& V) { return V.Multiplied(Scalar); } +inline gp_Vec2d gp_Vec2d::GetNormal() const + { + return gp_Vec2d(this->Y(), (-1)*this->X()); + } + diff --git a/tests/bugs/modalg_5/bug23706_1 b/tests/bugs/modalg_5/bug23706_1 new file mode 100755 index 0000000000..cc98cbe935 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_1 @@ -0,0 +1,32 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +bsplinecurve r2 4 3 1 5 2 1 3 5 0 8 0 1 2 8 2 1 4 8 3 1 4 8 3 1 6 8 4 1 10 8 10 1 +mkedge spine r2 +wire spine spine +circle profile 0 8 0 1 0 1 0.2 +mkedge profile profile +wire profile profile +mkplane profile profile +pipe p spine profile +explode p f +mksurface ss1 p_1 +mksurface ss2 p_2 +mksurface ss3 p_3 +offset o1 ss1 0.1 +offset o2 ss2 0.1 +offset o3 ss3 0.1 + +mkface res o2 +set info [sprops res] +regexp {Mass +: +([-0-9.+eE]+)} $info full sq +set sq_check 9.00819 + +if { [expr 1.*abs($sq_check - $sq)/$sq_check] > 0.01 } { + puts "Error : The square of result shape is $sq" +} diff --git a/tests/bugs/modalg_5/bug23706_10 b/tests/bugs/modalg_5/bug23706_10 new file mode 100755 index 0000000000..491b9f083e --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_10 @@ -0,0 +1,17 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +bsplinecurve r3 2 6 1 3 2 1 3 1 4 1 5 1 6 3 2 5 3 1 3 7 3 1 4 8 3 1 4 8 3 1 4 8 3 1 5 9 3 1 9 7 3 1 +bsplinecurve r4 2 6 2 3 2.5 1 3 1 3.5 1 4 1 4.5 3 -1 2 3 1 1 11 3 1 3 9 3 1 3 9 3 1 3 9 3 1 5 7 3 1 7 4 3 1 +set info [extrema r3 r4] + +if { [regexp "Extrema 3 is point : 4 8 3" $info] != 1 || [regexp "Extrema 8 is point : 4 8 3" $info] != 1 } { + puts "Error : Point of extrema is wrong" +} else { + puts "OK: Point of extrema is valid" +} diff --git a/tests/bugs/modalg_5/bug23706_11 b/tests/bugs/modalg_5/bug23706_11 new file mode 100755 index 0000000000..207023b2a8 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_11 @@ -0,0 +1,17 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +bsplinecurve r1 2 5 1 3 2 1 3 1 4 1 5 3 2 5 3 1 3 7 3 1 4 8 3 1 4 8 3 1 5 9 3 1 9 7 3 1 +bsplinecurve r2 2 5 2 3 2.5 1 3 1 3.5 1 4 3 -1 2 3 1 1 11 3 1 3 9 3 1 3 9 3 1 3 9 3 1 5 7 3 1 7 4 3 1 +set info [extrema r1 r2] + +if { [llength $info] != 3 } { + puts "Error : Extrema is wrong" +} else { + puts "OK: Extrema is valid" +} diff --git a/tests/bugs/modalg_5/bug23706_12 b/tests/bugs/modalg_5/bug23706_12 new file mode 100755 index 0000000000..743e4e01ac --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_12 @@ -0,0 +1,18 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +bsplinecurve r9 2 6 1 3 2 1 3 1 4 1 5 1 6 3 4 -3 3 1 6 8 3 1 10 11 3 1 10 11 3 1 10 11 3 1 14 14 3 1 5 8 3 1 +bsplinecurve r10 2 6 2 3 2.5 1 3 1 3.5 1 4 1 4.5 3 5 20 3 1 8 15 3 1 12 18 3 1 12 18 3 1 12 18 3 1 16 21 3 1 7 12 3 1 + +set info [extrema r9 r10] + +if { [llength $info] != 6 } { + puts "Error : Extrema is wrong" +} else { + puts "OK: Extrema is valid" +} diff --git a/tests/bugs/modalg_5/bug23706_13 b/tests/bugs/modalg_5/bug23706_13 new file mode 100755 index 0000000000..6a1f221f34 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_13 @@ -0,0 +1,26 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +2dbsplinecurve b3 2 6 1 3 2 1 3 1 4 1 5 1 6 3 2 5 1 3 7 1 4 8 1 4 8 1 4 8 1 5 9 1 9 7 1 +2dbsplinecurve b4 2 6 2 3 2.5 1 3 1 3.5 1 4 1 4.5 3 -1 2 1 1 11 1 3 9 1 3 9 1 3 9 1 5 7 1 7 4 1 + +set info [2dextrema b3 b4] +set status 0 +for { set i 1 } { $i <= 15 } { incr i 1 } { + regexp "dist $i: +(\[-0-9.+eE\]+)" $info full pp + if { $pp != 1.4142135623730951 } { + puts "Error : Extrema is wrong on dist $i" + set status 1 + } +} + +if { $status != 0 } { + puts "Error : Extrema is wrong" +} else { + puts "OK: Extrema is valid" +} diff --git a/tests/bugs/modalg_5/bug23706_14 b/tests/bugs/modalg_5/bug23706_14 new file mode 100755 index 0000000000..f15381d2b1 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_14 @@ -0,0 +1,49 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +2dbsplinecurve b9 2 8 1 2 2 1 3 1 4 1 5 1 6 1 7 1 8 2 4 -3 1 6 8 1 10 11 1 10 11 1 10 11 1 14 14 1 5 8 1 +2dbsplinecurve b10 2 8 2 2 2.5 1 3 1 3.5 1 4 1 4.5 1 5 1 5.5 2 5 20 1 8 15 1 12 18 1 12 18 1 12 18 1 16 21 1 7 12 1 +set info [2dextrema b9 b10] + +set status 0 +for { set i 1 } { $i <= 9 } { incr i } { + regexp "dist $i: +(\[-0-9.+eE\]+)" $info full pp1 + if { $pp1 != 7.2801098892805181 } { + puts "Error : Extrema is wrong on dist $i" + set status 1 + } +} + +for { set j 11 } { $j <= 19 } { incr j 1 } { + regexp "dist $j: +(\[-0-9.+eE\]+)" $info full pp2 + if { $pp2 != 7.2801098892805181 } { + puts "Error : Extrema is wrong on dist $j" + set status 1 + } +} + +regexp {dist 10: +([-0-9.+eE]+)} $info full pp3 +regexp {dist 20: +([-0-9.+eE]+)} $info full pp4 +regexp {dist 21: +([-0-9.+eE]+)} $info full pp5 +set pp_c 4.316921907096102 +set pp_ch 4.3169219070961038 + +if { $pp3 != $pp_c } { + puts "Error : Extrema is wrong on dist 10" + set status 1 +} +if { $pp4 != $pp_ch || $pp5 != $pp_ch} { + puts "Error : Extrema is wrong on dist 20 or 21" + set status 1 +} + +if { $status != 0 } { + puts "Error : Extrema is wrong" +} else { + puts "OK: Extrema is valid" +} diff --git a/tests/bugs/modalg_5/bug23706_15 b/tests/bugs/modalg_5/bug23706_15 new file mode 100755 index 0000000000..c5a7f0672e --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_15 @@ -0,0 +1,34 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +2dbsplinecurve b7 2 5 1 3 2 1 3 1 4 1 5 3 4 -3 1 6 8 1 10 11 1 10 11 1 14 14 1 5 8 1 +2dbsplinecurve b8 2 5 2 3 2.5 1 3 1 3.5 1 4 3 5 20 1 8 15 1 12 18 1 12 18 1 16 21 1 7 12 1 +set info [2dextrema b7 b8] + +set status 0 +for { set i 2 } { $i <= 5 } { incr i } { + regexp "dist $i: +(\[-0-9.+eE\]+)" $info full pp1 + if { $pp1 !=4.3624023150195192 } { + puts "Error : Extrema is wrong on dist $i" + set status 1 + } +} + +regexp {dist 1: +([-0-9.+eE]+)} $info full pp2 +set pp_ch 4.3624023150195184 + +if { $pp2 != $pp_ch } { + puts "Error : Extrema is wrong on dist 1" + set status 1 +} + +if { $status != 0 } { + puts "Error : Extrema is wrong" +} else { + puts "OK: Extrema is valid" +} diff --git a/tests/bugs/modalg_5/bug23706_16 b/tests/bugs/modalg_5/bug23706_16 new file mode 100644 index 0000000000..9fd030f097 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_16 @@ -0,0 +1,22 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.999999652077201 +set y 5.0000000062915735 +set z 5.00002142991819367 +set pp_ch 0.9991079538920743 + +restore [locate_data_file bug23706_c.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp +if { $pp != $pp_ch } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_17 b/tests/bugs/modalg_5/bug23706_17 new file mode 100644 index 0000000000..8928dff363 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_17 @@ -0,0 +1,22 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.99999991301930024 +set y 5.00000000157289337 +set z 5.00000535747954842 +set pp_ch 0.99955486819730277 + +restore [locate_data_file bug23706_c.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp +if { $pp != $pp_ch } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_19 b/tests/bugs/modalg_5/bug23706_19 new file mode 100644 index 0000000000..90a6521433 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_19 @@ -0,0 +1,22 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.99999999837571056 +set y 5.0000000000293724 +set z 5.0000001000463034 +set pp_ch 0.99993927567416474 + +restore [locate_data_file bug23706_c.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp +if { $pp != $pp_ch } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_2 b/tests/bugs/modalg_5/bug23706_2 new file mode 100755 index 0000000000..10cc9a1f43 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_2 @@ -0,0 +1,26 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +bsplinecurve r2 4 3 1 5 2 1 3 5 0 8 0 1 2 8 2 1 4 8 3 1 4 8 3 1 6 8 4 1 10 8 10 1 +mkedge spine r2 +wire spine spine +mksweep spine +addsweep spine -R +buildsweep spine -R +explode spine f +mksurface ss spine_1 +offset o1 ss 2 + +mkface res o1 +set info [sprops res] +regexp {Mass +: +([-0-9.+eE]+)} $info full sq +set sq_check 254.476 + +if { [expr 1.*abs($sq_check - $sq)/$sq_check] > 0.01 } { + puts "Error : The square of result shape is $sq" +} diff --git a/tests/bugs/modalg_5/bug23706_20 b/tests/bugs/modalg_5/bug23706_20 new file mode 100644 index 0000000000..a40965a0e7 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_20 @@ -0,0 +1,22 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.9999965207720098 +set y 5.0000000629157348 +set z 5.0002142991819367 +set pp_ch 0.99715423329884789 + +restore [locate_data_file bug23706_c2.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp +if { $pp != $pp_ch } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_21 b/tests/bugs/modalg_5/bug23706_21 new file mode 100644 index 0000000000..7f2cf8bf13 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_21 @@ -0,0 +1,22 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.999999652077201 +set y 5.0000000062915735 +set z 5.00002142991819367 +set pp_ch 0.99910795390933105 + +restore [locate_data_file bug23706_c2.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp +if { $pp != $pp_ch } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_22 b/tests/bugs/modalg_5/bug23706_22 new file mode 100644 index 0000000000..b95aa011f4 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_22 @@ -0,0 +1,22 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.99999991301930024 +set y 5.00000000157289337 +set z 5.00000535747954842 +set pp_ch 0.99955486819834238 + +restore [locate_data_file bug23706_c2.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp +if { $pp != $pp_ch } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_24 b/tests/bugs/modalg_5/bug23706_24 new file mode 100644 index 0000000000..8a2672b716 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_24 @@ -0,0 +1,22 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.99999999837571056 +set y 5.0000000000293724 +set z 5.0000001000463034 +set pp_ch 0.9999392756740122 + +restore [locate_data_file bug23706_c2.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp +if { $pp != $pp_ch } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_26 b/tests/bugs/modalg_5/bug23706_26 new file mode 100644 index 0000000000..bed1d5bfce --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_26 @@ -0,0 +1,25 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 11.0 +set y -6.0 +set z 5.0 +set pp_ch1 0.22894170490369878 +set pp_ch2 1.7710582950963012 + +restore [locate_data_file bug23706_c03.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} + diff --git a/tests/bugs/modalg_5/bug23706_27 b/tests/bugs/modalg_5/bug23706_27 new file mode 100644 index 0000000000..506d935440 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_27 @@ -0,0 +1,24 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.0 +set y 6.0 +set z -3.0 +set pp_ch1 1 +set pp_ch2 1 + +restore [locate_data_file bug23706_c03.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_28 b/tests/bugs/modalg_5/bug23706_28 new file mode 100644 index 0000000000..05e74545cc --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_28 @@ -0,0 +1,26 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 5.0 +set y 8.0 +set z -2.0 +set pp_ch1 1 +set pp_ch2 1 +set pp_ch3 1.1865241781930462 + +restore [locate_data_file bug23706_c03.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_29 b/tests/bugs/modalg_5/bug23706_29 new file mode 100644 index 0000000000..2d1bc735d0 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_29 @@ -0,0 +1,26 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x -4.0 +set y 4.0 +set z 1.0 +set pp_ch1 0 +set pp_ch2 1 +set pp_ch3 1 + +restore [locate_data_file bug23706_c03.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_3 b/tests/bugs/modalg_5/bug23706_3 new file mode 100755 index 0000000000..71c5592d2c --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_3 @@ -0,0 +1,23 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +restore [locate_data_file bug23706_c2.draw] c +mkedge e1 c +wire w e1 +polyline pl 3.9 4.9 5 4 5.1 5 4.1 3.9 5 3.9 4.9 5 +mkplane pl pl +pipe result w pl +don result +fit +checkshape result + +set square 1.86489 + + + + diff --git a/tests/bugs/modalg_5/bug23706_31 b/tests/bugs/modalg_5/bug23706_31 new file mode 100644 index 0000000000..6f53a40d4c --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_31 @@ -0,0 +1,25 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 11.0 +set y -6.0 +set z 5.0 +set pp_ch1 0.22894170490369881 +set pp_ch2 1.7710582950963012 + +restore [locate_data_file bug23706_c04.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} + diff --git a/tests/bugs/modalg_5/bug23706_32 b/tests/bugs/modalg_5/bug23706_32 new file mode 100644 index 0000000000..6f7bda2fd4 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_32 @@ -0,0 +1,24 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.0 +set y 6.0 +set z -3.0 +set pp_ch1 1 +set pp_ch2 1 + +restore [locate_data_file bug23706_c04.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_33 b/tests/bugs/modalg_5/bug23706_33 new file mode 100644 index 0000000000..2ca8763abe --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_33 @@ -0,0 +1,26 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 5.0 +set y 8.0 +set z -2.0 +set pp_ch1 0.81347582180695399 +set pp_ch2 1 +set pp_ch3 1 + +restore [locate_data_file bug23706_c04.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_34 b/tests/bugs/modalg_5/bug23706_34 new file mode 100644 index 0000000000..b14a05f809 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_34 @@ -0,0 +1,26 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x -4.0 +set y 4.0 +set z 1.0 +set pp_ch1 1 +set pp_ch2 1 +set pp_ch3 2 + +restore [locate_data_file bug23706_c04.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_36 b/tests/bugs/modalg_5/bug23706_36 new file mode 100644 index 0000000000..22bc5dd56e --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_36 @@ -0,0 +1,25 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 11.0 +set y -6.0 +set z 5.0 +set pp_ch1 0.22894170490369881 +set pp_ch2 1.7205732840814361 + +restore [locate_data_file bug23706_c05.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} + diff --git a/tests/bugs/modalg_5/bug23706_37 b/tests/bugs/modalg_5/bug23706_37 new file mode 100644 index 0000000000..4ffd361388 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_37 @@ -0,0 +1,24 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.0 +set y 6.0 +set z -3.0 +set pp_ch1 1 +set pp_ch2 1 + +restore [locate_data_file bug23706_c05.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_38 b/tests/bugs/modalg_5/bug23706_38 new file mode 100644 index 0000000000..d55ae36f0b --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_38 @@ -0,0 +1,26 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 5.0 +set y 8.0 +set z -2.0 +set pp_ch1 1 +set pp_ch2 1 +set pp_ch3 1.0371228345434986 + +restore [locate_data_file bug23706_c05.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_39 b/tests/bugs/modalg_5/bug23706_39 new file mode 100644 index 0000000000..bd679f6ca2 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_39 @@ -0,0 +1,26 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x -4.0 +set y 4.0 +set z 1.0 +set pp_ch1 0 +set pp_ch2 1 +set pp_ch3 1 + +restore [locate_data_file bug23706_c05.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_4 b/tests/bugs/modalg_5/bug23706_4 new file mode 100755 index 0000000000..e84fa2b97c --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_4 @@ -0,0 +1,37 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +2dbsplinecurve c1 2 5 0 3 0.2 1 0.3 1 0.4 1 0.5 3 2 0 1 3 -1 1 5 5 1 5 5 1 6 8 1 4 7 1 +2dbsplinecurve c3 2 4 1 3 2 1 3 1 5 3 6 3 1 5.001 5.01 1 5.001 5.01 1 3 9 1 2 11 1 +set info [2dintersect c1 c3] + +if { [regexp "Intersection point 1" $info] != 1 } { + puts "Error : Intersection should have two points" +} else { + regexp {Intersection point 1 +: +([-0-9.+eE]+)} $info full p11t + regexp {Intersection point 1 +: +[-0-9.+eE]+ +([-0-9.+eE]+)} $info full p12t +} + +if { [regexp "Intersection point 2" $info] != 1 } { + puts "Error : Intersection should have two points" +} else { + regexp {Intersection point 2 +: +([-0-9.+eE]+)} $info full p21t + regexp {Intersection point 2 +: +[-0-9.+eE]+ +([-0-9.+eE]+)} $info full p22t +} + +set p11 [expr round($p11t*10000)] +set p12 [expr round($p12t*10000)] +set p21 [expr round($p21t*10000)] +set p22 [expr round($p22t*10000)] + +if { ${p11} != 50024 || ${p12} != 50072 || ${p21} != 40024 || ${p22} != 70012 } { + puts "Error : Points of intersection have wrong coordinates" +} else { + puts "OK: Points of intersection are right" +} + diff --git a/tests/bugs/modalg_5/bug23706_40 b/tests/bugs/modalg_5/bug23706_40 new file mode 100644 index 0000000000..00fe14a61d --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_40 @@ -0,0 +1,22 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x -3.0 +set y 15.0 +set z -9.0 +set pp_ch 0.31967360381308058 + +restore [locate_data_file bug23706_c07.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp +if { $pp != $pp_ch } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_41 b/tests/bugs/modalg_5/bug23706_41 new file mode 100644 index 0000000000..73d21aef3e --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_41 @@ -0,0 +1,23 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 11.0 +set y -6.0 +set z 5.0 +set pp_ch 1.7205732840814361 + +restore [locate_data_file bug23706_c07.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp +if { $pp != $pp_ch } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} + diff --git a/tests/bugs/modalg_5/bug23706_42 b/tests/bugs/modalg_5/bug23706_42 new file mode 100644 index 0000000000..7d133ce83b --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_42 @@ -0,0 +1,24 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.0 +set y 6.0 +set z -3.0 +set pp_ch1 1 +set pp_ch2 1 + +restore [locate_data_file bug23706_c07.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_43 b/tests/bugs/modalg_5/bug23706_43 new file mode 100644 index 0000000000..bc482f13f8 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_43 @@ -0,0 +1,26 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 5.0 +set y 8.0 +set z -2.0 +set pp_ch1 1 +set pp_ch2 1 +set pp_ch3 1.0371228345434986 + +restore [locate_data_file bug23706_c07.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_44 b/tests/bugs/modalg_5/bug23706_44 new file mode 100644 index 0000000000..749b807d91 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_44 @@ -0,0 +1,26 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x -4.0 +set y 4.0 +set z 1.0 +set pp_ch1 0.034819847916144751 +set pp_ch2 1 +set pp_ch3 1 + +restore [locate_data_file bug23706_c07.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_45 b/tests/bugs/modalg_5/bug23706_45 new file mode 100644 index 0000000000..902082288d --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_45 @@ -0,0 +1,22 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x -3.0 +set y 15.0 +set z -9.0 +set pp_ch 0.51963826852813999 + +restore [locate_data_file bug23706_c08.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp +if { $pp != $pp_ch } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_46 b/tests/bugs/modalg_5/bug23706_46 new file mode 100644 index 0000000000..d5166afa33 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_46 @@ -0,0 +1,24 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 11.0 +set y -6.0 +set z 5.0 +set pp_ch1 0.48292970726418566 +set pp_ch2 1.7205732840814361 + +restore [locate_data_file bug23706_c08.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_47 b/tests/bugs/modalg_5/bug23706_47 new file mode 100644 index 0000000000..abaa114e9c --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_47 @@ -0,0 +1,24 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.0 +set y 6.0 +set z -3.0 +set pp_ch1 1 +set pp_ch2 1 + +restore [locate_data_file bug23706_c08.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_48 b/tests/bugs/modalg_5/bug23706_48 new file mode 100644 index 0000000000..8bf2f89506 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_48 @@ -0,0 +1,26 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 5.0 +set y 8.0 +set z -2.0 +set pp_ch1 1 +set pp_ch2 1 +set pp_ch3 1.0371228345434986 + +restore [locate_data_file bug23706_c08.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_49 b/tests/bugs/modalg_5/bug23706_49 new file mode 100644 index 0000000000..cd75d4d348 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_49 @@ -0,0 +1,26 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x -4.0 +set y 4.0 +set z 1.0 +set pp_ch1 0.087689905182099182 +set pp_ch2 1 +set pp_ch3 1 + +restore [locate_data_file bug23706_c08.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 } { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_5 b/tests/bugs/modalg_5/bug23706_5 new file mode 100755 index 0000000000..19d169fb6b --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_5 @@ -0,0 +1,34 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +2dbsplinecurve c1 2 5 0 3 0.2 1 0.3 1 0.4 1 0.5 3 2 0 1 3 -1 1 5 5 1 5 5 1 6 8 1 4 7 1 +2dbsplinecurve c2 2 4 1 3 2 1 3 1 5 3 6 3 1 5 5 1 5 5 1 3 9 1 2 11 1 +set info [2dintersect c1 c2] + +if { [regexp "Intersection point 1" $info] != 1 } { + puts "Error : Intersection should have two points" +} else { + regexp {Intersection point 1 +: +([-0-9.+eE]+)} $info full p11 + regexp {Intersection point 1 +: +[-0-9.+eE]+ +([-0-9.+eE]+)} $info full p12 +} + +if { [regexp "Intersection point 2" $info] != 1 } { + puts "Error : Intersection should have two points" +} else { + regexp {Intersection point 2 +: +([-0-9.+eE]+)} $info full p21t + regexp {Intersection point 2 +: +[-0-9.+eE]+ +([-0-9.+eE]+)} $info full p22t +} + +set p21 [expr int($p21t)] +set p22 [expr int($p22t)] + +if { ${p11} != 5 || ${p12} != 5 || ${p21} != 4 || ${p22} != 7 } { + puts "Error : Points of intersection have wrong coordinates" +} else { + puts "OK: Points of intersection are right" +} diff --git a/tests/bugs/modalg_5/bug23706_50 b/tests/bugs/modalg_5/bug23706_50 new file mode 100644 index 0000000000..18eedbbb98 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_50 @@ -0,0 +1,40 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.0 +set y 3.0 +set z 2.0 +set pp_ch1 2.261838779028444 +set pp_ch2 2.7514388736312116 +set pp_ch3 3.5195936992321921 +set pp_ch4 3.9600115496393977 +set pp_ch5 5.4999999987220543 +set pp_ch6 6.8388132447593541 +set pp_ch7 7.8261046366621292 + +restore [locate_data_file bug23706_c11.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4 +regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5 +regexp {parameter 6 += +([-0-9.+eE]+)} $info full pp6 +regexp {parameter 7 += +([-0-9.+eE]+)} $info full pp7 +if { $pp1 != $pp_ch1 || + $pp2 != $pp_ch2 || + $pp3 != $pp_ch3 || + $pp4 != $pp_ch4 || + $pp5 != $pp_ch5 || + $pp6 != $pp_ch6 || + $pp7 != $pp_ch7} { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_51 b/tests/bugs/modalg_5/bug23706_51 new file mode 100644 index 0000000000..1c2938bc87 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_51 @@ -0,0 +1,34 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 5.0 +set y 7.0 +set z 8.0 +set pp_ch1 2.8126840147763663 +set pp_ch2 3.5195936992321926 +set pp_ch3 3.9600115496393977 +set pp_ch4 5.4999999987220543 +set pp_ch5 7.2883607799598096 +set pp_ch6 1 + +restore [locate_data_file bug23706_c11.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4 +regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5 +regexp {parameter 6 += +([-0-9.+eE]+)} $info full pp6 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || + $pp3 != $pp_ch3 || $pp4 != $pp_ch4 || + $pp5 != $pp_ch5 || $pp6 != $pp_ch6} { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_52 b/tests/bugs/modalg_5/bug23706_52 new file mode 100644 index 0000000000..fd19bca607 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_52 @@ -0,0 +1,30 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 11.0 +set y -2.0 +set z -2.0 +set pp_ch1 2.9473269594602054 +set pp_ch2 4.4416670680933228 +set pp_ch3 5.4999999987220543 +set pp_ch4 6.6582576262308306 +set pp_ch5 7.7414573419084736 + +restore [locate_data_file bug23706_c11.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4 +regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 || $pp4 != $pp_ch4 || $pp5 != $pp_ch5} { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_53 b/tests/bugs/modalg_5/bug23706_53 new file mode 100644 index 0000000000..4fe42247be --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_53 @@ -0,0 +1,40 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.0 +set y 3.0 +set z 2.0 +set pp_ch1 1.1738953633378706 +set pp_ch2 2.1611867552406454 +set pp_ch3 3.5000000012779413 +set pp_ch4 5.0399884503606023 +set pp_ch5 5.4804063007678074 +set pp_ch6 6.2485611263687888 +set pp_ch7 6.7381612209715556 + +restore [locate_data_file bug23706_c12.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4 +regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5 +regexp {parameter 6 += +([-0-9.+eE]+)} $info full pp6 +regexp {parameter 7 += +([-0-9.+eE]+)} $info full pp7 +if { $pp1 != $pp_ch1 || + $pp2 != $pp_ch2 || + $pp3 != $pp_ch3 || + $pp4 != $pp_ch4 || + $pp5 != $pp_ch5 || + $pp6 != $pp_ch6 || + $pp7 != $pp_ch7} { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_54 b/tests/bugs/modalg_5/bug23706_54 new file mode 100644 index 0000000000..e3cbb582ba --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_54 @@ -0,0 +1,34 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 5.0 +set y 7.0 +set z 8.0 +set pp_ch1 1.7116392200401909 +set pp_ch2 3.5000000012779413 +set pp_ch3 5.0399884503606023 +set pp_ch4 5.4804063007678074 +set pp_ch5 6.1873159852236332 +set pp_ch6 8 + +restore [locate_data_file bug23706_c12.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4 +regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5 +regexp {parameter 6 += +([-0-9.+eE]+)} $info full pp6 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || + $pp3 != $pp_ch3 || $pp4 != $pp_ch4 || + $pp5 != $pp_ch5 || $pp6 != $pp_ch6} { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_55 b/tests/bugs/modalg_5/bug23706_55 new file mode 100644 index 0000000000..833dd50804 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_55 @@ -0,0 +1,30 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 11.0 +set y -2.0 +set z -2.0 +set pp_ch1 1.2585426580915264 +set pp_ch2 2.3417423737691694 +set pp_ch3 3.499999996505935 +set pp_ch4 4.5583329319066772 +set pp_ch5 6.052673040539795 + +restore [locate_data_file bug23706_c12.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4 +regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 || $pp4 != $pp_ch4 || $pp5 != $pp_ch5} { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_56 b/tests/bugs/modalg_5/bug23706_56 new file mode 100644 index 0000000000..585a783c77 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_56 @@ -0,0 +1,34 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.0 +set y 3.0 +set z 2.0 +set pp_ch1 1.8318851868378956 +set pp_ch2 3.0397214383562297 +set pp_ch3 5.5 +set pp_ch4 6.8388132447593541 +set pp_ch5 7.8261046366621292 + +restore [locate_data_file bug23706_c13.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4 +regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5 +if { $pp1 != $pp_ch1 || + $pp2 != $pp_ch2 || + $pp3 != $pp_ch3 || + $pp4 != $pp_ch4 || + $pp5 != $pp_ch5} { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_57 b/tests/bugs/modalg_5/bug23706_57 new file mode 100644 index 0000000000..d075ef7dd6 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_57 @@ -0,0 +1,29 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 5.0 +set y 7.0 +set z 8.0 +set pp_ch1 3.0397214383562297 +set pp_ch2 5.5 +set pp_ch3 7.2883607799598096 +set pp_ch4 1 + +restore [locate_data_file bug23706_c13.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || + $pp3 != $pp_ch3 || $pp4 != $pp_ch4} { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_58 b/tests/bugs/modalg_5/bug23706_58 new file mode 100644 index 0000000000..2be2338ea2 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_58 @@ -0,0 +1,30 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 11.0 +set y -2.0 +set z -2.0 +set pp_ch1 2.2389225099869194 +set pp_ch2 3.219764556283669 +set pp_ch3 5.5 +set pp_ch4 6.6582576262308306 +set pp_ch5 7.7414573419084736 + +restore [locate_data_file bug23706_c13.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4 +regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 || $pp4 != $pp_ch4 || $pp5 != $pp_ch5} { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_59 b/tests/bugs/modalg_5/bug23706_59 new file mode 100644 index 0000000000..58ded812e2 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_59 @@ -0,0 +1,34 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 3.0 +set y 3.0 +set z 2.0 +set pp_ch1 1.1738953633378706 +set pp_ch2 2.1611867552406454 +set pp_ch3 3.5000000000000004 +set pp_ch4 5.9602785616437703 +set pp_ch5 7.1681148131621049 + +restore [locate_data_file bug23706_c14.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4 +regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5 +if { $pp1 != $pp_ch1 || + $pp2 != $pp_ch2 || + $pp3 != $pp_ch3 || + $pp4 != $pp_ch4 || + $pp5 != $pp_ch5} { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_6 b/tests/bugs/modalg_5/bug23706_6 new file mode 100755 index 0000000000..8314b56cc6 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_6 @@ -0,0 +1,18 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +pload XSDRAW +2dbsplinecurve cc 3 2 0 4 1 4 -1 -1 1 0 -1 1 0 0 1 0 0 1 +offset2dcurve o cc .5 +set info [length o] +regexp {The length o is+ +([-0-9.+eE]+)} $info full ll +set ll_check 2.3717833300483151 + +if { [expr 1.*abs($ll_check - $ll)/$ll_check] > 0.01 } { + puts "Error : The lenght of result shape is $ll" +} diff --git a/tests/bugs/modalg_5/bug23706_60 b/tests/bugs/modalg_5/bug23706_60 new file mode 100644 index 0000000000..5ec0c89a2a --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_60 @@ -0,0 +1,29 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 5.0 +set y 7.0 +set z 8.0 +set pp_ch1 1.7116392200401909 +set pp_ch2 3.5000000000000004 +set pp_ch3 5.9602785616437703 +set pp_ch4 8 + +restore [locate_data_file bug23706_c14.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || + $pp3 != $pp_ch3 || $pp4 != $pp_ch4} { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_61 b/tests/bugs/modalg_5/bug23706_61 new file mode 100644 index 0000000000..8057244123 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_61 @@ -0,0 +1,30 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +#################################### +## Cannot project point on curve +#################################### + +set x 11.0 +set y -2.0 +set z -2.0 +set pp_ch1 1.2585426580915264 +set pp_ch2 2.3417423737691694 +set pp_ch3 3.5000000000000004 +set pp_ch4 5.7802354437163306 +set pp_ch5 6.761077490013081 + +restore [locate_data_file bug23706_c14.draw] c +set info [proj c $x $y $z] + +regexp {parameter 1 += +([-0-9.+eE]+)} $info full pp1 +regexp {parameter 2 += +([-0-9.+eE]+)} $info full pp2 +regexp {parameter 3 += +([-0-9.+eE]+)} $info full pp3 +regexp {parameter 4 += +([-0-9.+eE]+)} $info full pp4 +regexp {parameter 5 += +([-0-9.+eE]+)} $info full pp5 +if { $pp1 != $pp_ch1 || $pp2 != $pp_ch2 || $pp3 != $pp_ch3 || $pp4 != $pp_ch4 || $pp5 != $pp_ch5} { + puts "Error : Projection is not correct" +} else { + puts "OK: Projection is correct" +} diff --git a/tests/bugs/modalg_5/bug23706_7 b/tests/bugs/modalg_5/bug23706_7 new file mode 100755 index 0000000000..ae4941f0a6 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_7 @@ -0,0 +1,19 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +pload XSDRAW +bsplinecurve cc 3 2 0 4 1 4 -1 -1 2 1 0 -1 2 1 0 0 2 1 0 0 2 1 +point pp 0 0 1 +offsetcurve o cc .5 pp +set info [length o] +regexp {The length o is+ +([-0-9.+eE]+)} $info full ll +set ll_check 2.3717833300483151 + +if { [expr 1.*abs($ll_check - $ll)/$ll_check] > 0.01 } { + puts "Error : The lenght of result shape is $ll" +} diff --git a/tests/bugs/modalg_5/bug23706_8 b/tests/bugs/modalg_5/bug23706_8 new file mode 100755 index 0000000000..2f7a3c6a86 --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_8 @@ -0,0 +1,25 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +bsplinecurve r2 2 4 1 3 2 2 3 1 4 3 0 8 0 1 2 8 2 1 4 8 3 1 4 8 3 1 6 8 4 1 10 8 10 1 +mkedge spine r2 +wire spine spine +circle profile 0 8 0 1 0 1 0.51 +mkedge profile profile +wire profile profile +mkplane profile profile +pipe p spine profile +explode p f +mksurface ss p_2 +bsplinecurve r3 2 5 1 3 2 1 3 1 4 1 5 3 9 7 3 1 5.447213595499958 9 2.105572809000084 1 4.223606797749979 8 2.552786404500042 1 4.223606797749979 8 2.552786404500042 1 3 7 3 1 2 5 3 1 +set info [intersect res r3 ss] +if { [regexp "res_1" $info] != 1 || [regexp "res_2" $info] != 1 || [regexp "res_3" $info] != 0 } { + puts "Error : Wrong number of points of intersection. Should contain two points" +} else { + puts "OK: Number of points of intersection is valid" +} diff --git a/tests/bugs/modalg_5/bug23706_9 b/tests/bugs/modalg_5/bug23706_9 new file mode 100755 index 0000000000..bb7d45db3f --- /dev/null +++ b/tests/bugs/modalg_5/bug23706_9 @@ -0,0 +1,19 @@ +puts "============" +puts "OCC23706" +puts "============" +puts "" +######################################################################### +# Cannot project point on curve +######################################################################### + +pload XSDRAW +2dbsplinecurve c1 2 5 0 3 0.2 1 0.3 1 0.4 1 0.5 3 2 0 1 3 -1 1 5 5 1 5 5 1 6 8 1 4 7 1 +offset2dcurve o1 c1 2 + +set info [length o1] +regexp {The length o1 is+ +([-0-9.+eE]+)} $info full ll +set ll_check 19.244437838214424 + +if { [expr 1.*abs($ll_check - $ll)/$ll_check] > 0.01 } { + puts "Error : The lenght of result shape is $ll" +} diff --git a/tests/bugs/moddata_3/bug23706 b/tests/bugs/moddata_3/bug23706 new file mode 100644 index 0000000000..42092736f6 --- /dev/null +++ b/tests/bugs/moddata_3/bug23706 @@ -0,0 +1,20 @@ +puts "========" +puts "OCC23706" +puts "========" +puts "" +############################### +## Cannot project point 3.9999965207720098 5.0000000629157348 5.0002142991819367 on curve from attached file. +############################### + +set BugNumber OCC23706 + +restore [locate_data_file bug23706_c.draw] c +proj c 3.9999965207720098 5.0000000629157348 5.0002142991819367 + +if {[isdraw ext_1] == 1} { + dump ext_1 + puts "${BugNumber} OK" +} else { + puts "Faulty ${BugNumber}" +} + \ No newline at end of file