From d30895f5da57fcf4b409163725c35ef2c823fd43 Mon Sep 17 00:00:00 2001 From: nbv Date: Fri, 23 Sep 2016 17:24:52 +0300 Subject: [PATCH] 0023178: Intersection of cylinders fails to produce results 1. Unification of trimmed and not-trimmed cylinders processing (IntPatch_Intersection::GeomGeomPerfomTrimSurf() method has been removed). 2. Interface of IntPatch_ImpImpIntersection::Perform(...) method has been changed. 3. Now, WLine purging is forbidden for Geom-Geom-Intersection. 4. Bnd_Range class has been created. See Bnd_Range.hxx for detail information. 5. Algorithm of AddBoundaryPoint function has been improved in order to obtain intersection points in both boundaries (VFirst and VLast of every surface). 6. Earlier, method Geom2dConvert::ConcatG1(...) increased resulted B-spline degree (in case of not succession of previous iteration). Now increased value has been limited by Geom2d_BSplineCurve::MaxDegree() value (max degree = 25). 7. Algorithm of B-spline closure definition has been changed in the methods Geom2dConvert::C0BSplineToC1BSplineCurve(...) and Geom2dConvert::C0BSplineToArrayOfC1BSplineCurve(...). Creation of test case for this issue. Adjusting test cases according to their new behavior. Small correction in the code according to KGV's remark. --- src/BSplCLib/BSplCLib_2.cxx | 7 +- src/Bnd/Bnd_Range.cxx | 37 + src/Bnd/Bnd_Range.hxx | 124 ++ src/Bnd/FILES | 2 + src/Geom2dConvert/Geom2dConvert.cxx | 49 +- src/IntPatch/IntPatch_ImpImpIntersection.hxx | 13 +- .../IntPatch_ImpImpIntersection_1.gxx | 13 +- .../IntPatch_ImpImpIntersection_2.gxx | 116 +- .../IntPatch_ImpImpIntersection_4.gxx | 1685 +++++++++-------- src/IntPatch/IntPatch_Intersection.cxx | 139 +- src/IntPatch/IntPatch_Intersection.hxx | 8 - src/IntPatch/IntPatch_WLineTool.cxx | 159 +- src/IntPatch/IntPatch_WLineTool.hxx | 3 + tests/blend/simple/Q6 | 2 +- tests/bugs/modalg_5/bug25742_2 | 19 +- tests/bugs/modalg_6/bug23178 | 43 + tests/bugs/modalg_6/bug27310_2 | 2 +- tests/bugs/modalg_6/bug27766 | 2 +- 18 files changed, 1351 insertions(+), 1072 deletions(-) create mode 100644 src/Bnd/Bnd_Range.cxx create mode 100644 src/Bnd/Bnd_Range.hxx create mode 100644 tests/bugs/modalg_6/bug23178 diff --git a/src/BSplCLib/BSplCLib_2.cxx b/src/BSplCLib/BSplCLib_2.cxx index 41b0835a30..6d7cbad466 100644 --- a/src/BSplCLib/BSplCLib_2.cxx +++ b/src/BSplCLib/BSplCLib_2.cxx @@ -31,8 +31,6 @@ #include -static const Standard_Integer aGlobalMaxDegree = 25; - //======================================================================= //struct : BSplCLib_DataContainer //purpose: Auxiliary structure providing buffers for poles and knots used in @@ -44,8 +42,7 @@ struct BSplCLib_DataContainer BSplCLib_DataContainer(Standard_Integer Degree) { (void)Degree; // avoid compiler warning - Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() || - BSplCLib::MaxDegree() > aGlobalMaxDegree, + Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree(), "BSplCLib: bspline degree is greater than maximum supported"); } @@ -326,7 +323,7 @@ BSplCLib::BuildBSpMatrix(const TColStd_Array1OfReal& Parameters, ReturnCode = 0, FirstNonZeroBsplineIndex, BandWidth, - MaxOrder = aGlobalMaxDegree+1, + MaxOrder = BSplCLib::MaxDegree() + 1, Order ; math_Matrix BSplineBasis(1, MaxOrder, diff --git a/src/Bnd/Bnd_Range.cxx b/src/Bnd/Bnd_Range.cxx new file mode 100644 index 0000000000..d79b211dea --- /dev/null +++ b/src/Bnd/Bnd_Range.cxx @@ -0,0 +1,37 @@ +// Created on: 2016-06-07 +// Created by: Nikolai BUKHALOV +// Copyright (c) 2016 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + + +#include + +//======================================================================= +//function : Common +//purpose : +//======================================================================= +void Bnd_Range::Common(const Bnd_Range& theOther) +{ + if(theOther.IsVoid()) + { + SetVoid(); + } + + if(IsVoid()) + { + return; + } + + myFirst = Max(myFirst, theOther.myFirst); + myLast = Min(myLast, theOther.myLast); +} diff --git a/src/Bnd/Bnd_Range.hxx b/src/Bnd/Bnd_Range.hxx new file mode 100644 index 0000000000..439008f1f8 --- /dev/null +++ b/src/Bnd/Bnd_Range.hxx @@ -0,0 +1,124 @@ +// Created on: 2016-06-07 +// Created by: Nikolai BUKHALOV +// Copyright (c) 2016 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#ifndef _Bnd_Range_HeaderFile +#define _Bnd_Range_HeaderFile + +#include +#include + +//! This class describes a range in 1D space restricted +//! by two real values. +//! A range can be void indicating there is no point included in the range. + +class Bnd_Range +{ +public: + + //! Default constructor. Creates VOID range. + Bnd_Range() : myFirst(0.0), myLast(-1.0) + { + }; + + //! Constructor. Never creates VOID range. + Bnd_Range(const Standard_Real theMin, const Standard_Real theMax) : + myFirst(theMin), myLast(theMax) + { + if(myLast < myFirst) + Standard_ConstructionError::Raise("Last < First"); + }; + + //! Replaces with common-part of and theOther + Standard_EXPORT void Common(const Bnd_Range& theOther); + + //! Extends to include theParameter + void Add(const Standard_Real theParameter) + { + if(IsVoid()) + { + myFirst = myLast = theParameter; + return; + } + + myFirst = Min(myFirst, theParameter); + myLast = Max(myLast, theParameter); + } + + //! Obtain MIN boundary of . + //! If is VOID the method returns false. + Standard_Boolean GetMin(Standard_Real& thePar) const + { + if(IsVoid()) + { + return Standard_False; + } + + thePar = myFirst; + return Standard_True; + } + + //! Obtain MAX boundary of . + //! If is VOID the method returns false. + Standard_Boolean GetMAX(Standard_Real& thePar) const + { + if(IsVoid()) + { + return Standard_False; + } + + thePar = myLast; + return Standard_True; + } + + //! Returns range value (MAX-MIN). Returns negative value for VOID range. + Standard_Real Delta() const + { + return (myLast - myFirst); + } + + //! Is initialized. + Standard_Boolean IsVoid() const + { + return (myLast < myFirst); + } + + //! Initializes by default parameters. Makes VOID. + void SetVoid() + { + myLast = -1.0; + myFirst = 0.0; + } + + //! Extends this to the given value (in both side) + void Enlarge(const Standard_Real theDelta) + { + if (IsVoid()) + { + return; + } + + myFirst -= theDelta; + myLast += theDelta; + } + +private: + //! Start of range + Standard_Real myFirst; + + //! End of range + Standard_Real myLast; +}; + +#endif \ No newline at end of file diff --git a/src/Bnd/FILES b/src/Bnd/FILES index 6c20debbee..845b5257e8 100644 --- a/src/Bnd/FILES +++ b/src/Bnd/FILES @@ -24,6 +24,8 @@ Bnd_Box2d.hxx Bnd_HArray1OfBox.hxx Bnd_HArray1OfBox2d.hxx Bnd_HArray1OfSphere.hxx +Bnd_Range.cxx +Bnd_Range.hxx Bnd_SeqOfBox.hxx Bnd_Sphere.cxx Bnd_Sphere.hxx diff --git a/src/Geom2dConvert/Geom2dConvert.cxx b/src/Geom2dConvert/Geom2dConvert.cxx index 8b952ac081..9190434bc3 100644 --- a/src/Geom2dConvert/Geom2dConvert.cxx +++ b/src/Geom2dConvert/Geom2dConvert.cxx @@ -1025,6 +1025,7 @@ void Geom2dConvert::ConcatG1(TColGeom2d_Array1OfBSplineCurve& ArrayOf for (j=1;j<=nb_curve-1;j++){ //boucle secondaire a l'interieur de chaque groupe Curve1=ArrayOfCurves(j); if ( (j==(nb_curve-1)) &&(Need2DegRepara(ArrayOfCurves))){ + const Standard_Integer aNewCurveDegree = Min(2 * Curve1->Degree(), Geom2d_BSplineCurve::MaxDegree()); Curve2->D1(Curve2->LastParameter(),Pint,Vec1); Curve1->D1(Curve1->FirstParameter(),Pint,Vec2); lambda=Vec2.Magnitude()/Vec1.Magnitude(); @@ -1074,7 +1075,7 @@ void Geom2dConvert::ConcatG1(TColGeom2d_Array1OfBSplineCurve& ArrayOf Curve1FlatKnots, Curve1Poles, FlatKnots, - 2*Curve1->Degree(), + aNewCurveDegree, NewPoles, Status ); @@ -1084,7 +1085,7 @@ void Geom2dConvert::ConcatG1(TColGeom2d_Array1OfBSplineCurve& ArrayOf Curve1FlatKnots, Curve1Weights, FlatKnots, - 2*Curve1->Degree(), + aNewCurveDegree, NewWeights, Status ); @@ -1110,7 +1111,7 @@ void Geom2dConvert::ConcatG1(TColGeom2d_Array1OfBSplineCurve& ArrayOf for (ii=1;ii<=NewPoles.Length();ii++) for (jj=1;jj<=2;jj++) NewPoles(ii).SetCoord(jj,NewPoles(ii).Coord(jj)/NewWeights(ii)); - Curve1= new Geom2d_BSplineCurve(NewPoles,NewWeights,KnotC1,KnotC1Mults,2*Curve1->Degree()); + Curve1 = new Geom2d_BSplineCurve(NewPoles, NewWeights, KnotC1, KnotC1Mults, aNewCurveDegree); } Geom2dConvert_CompCurveToBSplineCurve C(Curve2); fusion=C.Add(Curve1, @@ -1276,6 +1277,8 @@ void Geom2dConvert::ConcatC1(TColGeom2d_Array1OfBSplineCurve& ArrayOf else Curve1=ArrayOfCurves(j); + const Standard_Integer aNewCurveDegree = Min(2 * Curve1->Degree(), Geom2d_BSplineCurve::MaxDegree()); + if (j==0) //initialisation en debut de groupe Curve2=Curve1; else{ @@ -1315,7 +1318,7 @@ void Geom2dConvert::ConcatC1(TColGeom2d_Array1OfBSplineCurve& ArrayOf TColStd_Array1OfReal FlatKnots(1,Curve1FlatKnots.Length()+(Curve1->Degree()*Curve1->NbKnots())); BSplCLib::KnotSequence(KnotC1,KnotC1Mults,FlatKnots); - TColgp_Array1OfPnt2d NewPoles(1,FlatKnots.Length()-(2*Curve1->Degree()+1)); + TColgp_Array1OfPnt2d NewPoles(1, FlatKnots.Length() - (aNewCurveDegree + 1)); Standard_Integer Status; TColStd_Array1OfReal Curve1Weights(1,Curve1->NbPoles()); Curve1->Weights(Curve1Weights); @@ -1330,18 +1333,18 @@ void Geom2dConvert::ConcatC1(TColGeom2d_Array1OfBSplineCurve& ArrayOf Curve1FlatKnots, Curve1Poles, FlatKnots, - 2*Curve1->Degree(), + aNewCurveDegree, NewPoles, Status ); - TColStd_Array1OfReal NewWeights(1,FlatKnots.Length()-(2*Curve1->Degree()+1)); + TColStd_Array1OfReal NewWeights(1, FlatKnots.Length() - (aNewCurveDegree + 1)); // BSplCLib::FunctionReparameterise(reparameterise_evaluator, BSplCLib::FunctionReparameterise(ev, Curve1->Degree(), Curve1FlatKnots, Curve1Weights, FlatKnots, - 2*Curve1->Degree(), + aNewCurveDegree, NewWeights, Status ); @@ -1349,7 +1352,7 @@ void Geom2dConvert::ConcatC1(TColGeom2d_Array1OfBSplineCurve& ArrayOf for (jj=1;jj<=2;jj++) NewPoles(ii).SetCoord(jj,NewPoles(ii).Coord(jj)/NewWeights(ii)); } - Curve1= new Geom2d_BSplineCurve(NewPoles,NewWeights,KnotC1,KnotC1Mults,2*Curve1->Degree()); + Curve1 = new Geom2d_BSplineCurve(NewPoles, NewWeights, KnotC1, KnotC1Mults, aNewCurveDegree); } Geom2dConvert_CompCurveToBSplineCurve C(Curve2); fusion=C.Add(Curve1, @@ -1420,7 +1423,7 @@ void Geom2dConvert::C0BSplineToC1BSplineCurve(Handle(Geom2d_BSplineCurve)& BS, Standard_Integer i,j,nbcurveC1=1; Standard_Real U1,U2; Standard_Boolean closed_flag = Standard_False ; - gp_Pnt2d point; + gp_Pnt2d point1, point2; gp_Vec2d V1,V2; Standard_Boolean fusion; @@ -1453,15 +1456,20 @@ void Geom2dConvert::C0BSplineToC1BSplineCurve(Handle(Geom2d_BSplineCurve)& BS, BSbis->Segment(U1,U2); ArrayOfCurves(i)=BSbis; } + + const Standard_Real anAngularToler = 1.0e-7; Handle(TColStd_HArray1OfInteger) ArrayOfIndices; Handle(TColGeom2d_HArray1OfBSplineCurve) ArrayOfConcatenated; - BS->D1(BS->FirstParameter(),point,V1); //a verifier - BS->D1(BS->LastParameter(),point,V2); + BS->D1(BS->FirstParameter(),point1,V1); //a verifier + BS->D1(BS->LastParameter(),point2,V2); - if ((BS->IsClosed())&&(V1.IsParallel(V2,Precision::Confusion()))) - closed_flag = Standard_True ; - + if ((point1.SquareDistance(point2) < gp::Resolution()) && + (V1.IsParallel(V2, anAngularToler))) + { + closed_flag = Standard_True; + } + Geom2dConvert::ConcatC1(ArrayOfCurves, ArrayOfToler, ArrayOfIndices, @@ -1510,7 +1518,7 @@ void Geom2dConvert::C0BSplineToArrayOfC1BSplineCurve(const Handle(Geom2d_BSpline Standard_Integer i,j,nbcurveC1=1; Standard_Real U1,U2; Standard_Boolean closed_flag = Standard_False ; - gp_Pnt2d point; + gp_Pnt2d point1, point2; gp_Vec2d V1,V2; // Standard_Boolean fusion; @@ -1544,11 +1552,14 @@ void Geom2dConvert::C0BSplineToArrayOfC1BSplineCurve(const Handle(Geom2d_BSpline Handle(TColStd_HArray1OfInteger) ArrayOfIndices; - BS->D1(BS->FirstParameter(),point,V1); - BS->D1(BS->LastParameter(),point,V2); + BS->D1(BS->FirstParameter(),point1,V1); + BS->D1(BS->LastParameter(),point2,V2); - if ((BS->IsClosed())&&(V1.IsParallel(V2,AngularTolerance))) - closed_flag = Standard_True ; + if (((point1.SquareDistance(point2) < gp::Resolution())) && + (V1.IsParallel(V2, AngularTolerance))) + { + closed_flag = Standard_True; + } Geom2dConvert::ConcatC1(ArrayOfCurves, ArrayOfToler, diff --git a/src/IntPatch/IntPatch_ImpImpIntersection.hxx b/src/IntPatch/IntPatch_ImpImpIntersection.hxx index 23153ebe3e..ece87d6db0 100644 --- a/src/IntPatch/IntPatch_ImpImpIntersection.hxx +++ b/src/IntPatch/IntPatch_ImpImpIntersection.hxx @@ -62,7 +62,14 @@ public: //! When intersection result returns IntPatch_RLine and another //! IntPatch_Line (not restriction) we (in case of theIsReqToKeepRLine==TRUE) //! will always keep both lines even if they are coincided. - Standard_EXPORT void Perform (const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_TopolTool)& D1, const Handle(Adaptor3d_HSurface)& S2, const Handle(Adaptor3d_TopolTool)& D2, const Standard_Real TolArc, const Standard_Real TolTang, const Standard_Boolean isTheTrimmed = Standard_False, const Standard_Boolean theIsReqToKeepRLine = Standard_False); + Standard_EXPORT void Perform (const Handle(Adaptor3d_HSurface)& S1, + const Handle(Adaptor3d_TopolTool)& D1, + const Handle(Adaptor3d_HSurface)& S2, + const Handle(Adaptor3d_TopolTool)& D2, + const Standard_Real TolArc, + const Standard_Real TolTang, + const Standard_Boolean theIsReqToKeepRLine = + Standard_False); //! Returns True if the calculus was succesfull. Standard_Boolean IsDone() const; @@ -71,8 +78,8 @@ public: Standard_Boolean IsEmpty() const; //! Returns True if the two patches are considered as - //! entierly tangent, i-e every restriction arc of one - //! patch is inside the geometric base of the otehr patch. + //! entirely tangent, i.e every restriction arc of one + //! patch is inside the geometric base of the other patch. Standard_Boolean TangentFaces() const; //! Returns True when the TangentFaces returns True and the diff --git a/src/IntPatch/IntPatch_ImpImpIntersection_1.gxx b/src/IntPatch/IntPatch_ImpImpIntersection_1.gxx index 18146a1684..e2f049e236 100644 --- a/src/IntPatch/IntPatch_ImpImpIntersection_1.gxx +++ b/src/IntPatch/IntPatch_ImpImpIntersection_1.gxx @@ -70,16 +70,7 @@ static void ProcessBounds(const Handle(IntPatch_ALine)&, const Standard_Real); -static Standard_Boolean IntCyCy(const IntSurf_Quadric&, - const IntSurf_Quadric&, - const Standard_Real, - Standard_Boolean&, - Standard_Boolean&, - Standard_Boolean&, - IntPatch_SequenceOfLine&, - IntPatch_SequenceOfPoint&); - -static Standard_Boolean IntCyCyTrim(const IntSurf_Quadric& theQuad1, +static Standard_Boolean IntCyCy(const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, const Standard_Real theTol3D, const Standard_Real theTol2D, @@ -87,6 +78,8 @@ static Standard_Boolean IntCyCyTrim(const IntSurf_Quadric& theQuad1, const Bnd_Box2d& theUVSurf2, const Standard_Boolean isTheReverse, Standard_Boolean& isTheEmpty, + Standard_Boolean& isTheSameSurface, + Standard_Boolean& isTheMultiplePoint, IntPatch_SequenceOfLine& theSlin, IntPatch_SequenceOfPoint& theSPnt); diff --git a/src/IntPatch/IntPatch_ImpImpIntersection_2.gxx b/src/IntPatch/IntPatch_ImpImpIntersection_2.gxx index b1d3dfcbfd..5972adfd17 100644 --- a/src/IntPatch/IntPatch_ImpImpIntersection_2.gxx +++ b/src/IntPatch/IntPatch_ImpImpIntersection_2.gxx @@ -14,6 +14,8 @@ // Alternatively, this file may be used under the terms of Open CASCADE // commercial license or contractual agreement. +#include + static Standard_Integer SetQuad(const Handle(Adaptor3d_HSurface)& theS, GeomAbs_SurfaceType& theTS, @@ -40,7 +42,7 @@ IntPatch_ImpImpIntersection::IntPatch_ImpImpIntersection const Standard_Real TolTang, const Standard_Boolean theIsReqToKeepRLine) { - Perform(S1,D1,S2,D2,TolArc,TolTang, Standard_False, theIsReqToKeepRLine); + Perform(S1,D1,S2,D2,TolArc,TolTang, theIsReqToKeepRLine); } //======================================================================= //function : Perform @@ -52,13 +54,14 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_TopolTool)& D2, const Standard_Real TolArc, const Standard_Real TolTang, - const Standard_Boolean isTheTrimmed, - const Standard_Boolean theIsReqToKeepRLine) { + const Standard_Boolean theIsReqToKeepRLine) +{ done = Standard_False; - Standard_Boolean isTrimmed = isTheTrimmed; spnt.Clear(); slin.Clear(); + Standard_Boolean isPostProcessingRequired = Standard_True; + empt = Standard_True; tgte = Standard_False; oppo = Standard_False; @@ -149,59 +152,53 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)& S1, case 22: { // Cylinder/Cylinder Standard_Boolean isDONE = Standard_False; - - if(!isTrimmed) + Bnd_Box2d aBox1, aBox2; + + const Standard_Real aU1f = S1->FirstUParameter(); + Standard_Real aU1l = S1->LastUParameter(); + const Standard_Real aU2f = S2->FirstUParameter(); + Standard_Real aU2l = S2->LastUParameter(); + + const Standard_Real anUperiod = 2.0*M_PI; + + if(aU1l - aU1f > anUperiod) + aU1l = aU1f + anUperiod; + + if(aU2l - aU2f > anUperiod) + aU2l = aU2f + anUperiod; + + aBox1.Add(gp_Pnt2d(aU1f, S1->FirstVParameter())); + aBox1.Add(gp_Pnt2d(aU1l, S1->LastVParameter())); + aBox2.Add(gp_Pnt2d(aU2f, S2->FirstVParameter())); + aBox2.Add(gp_Pnt2d(aU2l, S2->LastVParameter())); + + // Resolution is too big if the cylinder radius is + // too small. Therefore, we shall bounde its value above. + // Here, we use simple constant. + const Standard_Real a2DTol = Min(1.0e-4, Min( S1->UResolution(TolTang), + S2->UResolution(TolTang))); + + //The bigger range the bigger nunber of points in Walking-line (WLine) + //we will be able to add and, consequently, we will obtain more + //precise intersection line. + //Every point of WLine is determined as function from U1-parameter, + //where U1 is U-parameter on 1st quadric. + //Therefore, we should use quadric with bigger range as 1st parameter + //in IntCyCy() function. + //On the other hand, there is no point in reversing in case of + //analytical intersection (when result is line, ellipse, point...). + //This result is independent of the arguments order. + Standard_Boolean isReversed = ((aU2l - aU2f) < (aU1l - aU1f)); + + if(isReversed) { - isDONE = IntCyCy(quad1, quad2, TolTang, empt, - SameSurf, multpoint, slin, spnt); + isDONE = IntCyCy(quad2, quad1, TolTang, a2DTol, aBox2, aBox1, + Standard_True, empt, SameSurf, multpoint, slin, spnt); } else { - Bnd_Box2d aBox1, aBox2; - - const Standard_Real aU1f = S1->FirstUParameter(); - Standard_Real aU1l = S1->LastUParameter(); - const Standard_Real aU2f = S2->FirstUParameter(); - Standard_Real aU2l = S2->LastUParameter(); - - const Standard_Real anUperiod = 2.0*M_PI; - - if(aU1l - aU1f > anUperiod) - aU1l = aU1f + anUperiod; - - if(aU2l - aU2f > anUperiod) - aU2l = aU2f + anUperiod; - - aBox1.Add(gp_Pnt2d(aU1f, S1->FirstVParameter())); - aBox1.Add(gp_Pnt2d(aU1l, S1->LastVParameter())); - aBox2.Add(gp_Pnt2d(aU2f, S2->FirstVParameter())); - aBox2.Add(gp_Pnt2d(aU2l, S2->LastVParameter())); - - // Resolution is too big if the cylinder radius is - // too small. Therefore, we shall bounde its value above. - // Here, we use simple constant. - const Standard_Real a2DTol = Min(1.0e-4, Min( S1->UResolution(TolTang), - S2->UResolution(TolTang))); - - Standard_Boolean isReversed = ((aU2l - aU2f) < (aU1l - aU1f)); - - if(isReversed) - { - isDONE = IntCyCyTrim(quad2, quad1, TolTang, a2DTol, aBox2, aBox1, - isReversed, empt, slin, spnt); - } - else - { - isDONE = IntCyCyTrim(quad1, quad2, TolTang, a2DTol, aBox1, aBox2, - isReversed, empt, slin, spnt); - } - - if(!isDONE) - { - isDONE = IntCyCy(quad1, quad2, TolTang, empt, - SameSurf, multpoint, slin, spnt); - isTrimmed = Standard_False; - } + isDONE = IntCyCy(quad1, quad2, TolTang, a2DTol, aBox1, aBox2, + Standard_False, empt, SameSurf, multpoint, slin, spnt); } if (!isDONE) @@ -210,6 +207,17 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)& S1, } bEmpty = empt; + if(!slin.IsEmpty()) + { + const Handle(IntPatch_WLine)& aWLine = + Handle(IntPatch_WLine)::DownCast(slin.Value(1)); + + if(!aWLine.IsNull()) + {//No geometric solution + isPostProcessingRequired = Standard_False; + } + } + break; } // @@ -302,7 +310,7 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)& S1, } // - if(!isTrimmed) + if(isPostProcessingRequired) { if (!SameSurf) { AFunc.SetQuadric(quad2); diff --git a/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx b/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx index cee3358877..f5aca30301 100644 --- a/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx +++ b/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx @@ -15,17 +15,20 @@ // commercial license or contractual agreement. #include -#include -#include +#include #include +#include #include - +#include #include //If Abs(a) <= aNulValue then it is considered that a = 0. static const Standard_Real aNulValue = 1.0e-11; -struct stCoeffsValue; +static void ShortCosForm( const Standard_Real theCosFactor, + const Standard_Real theSinFactor, + Standard_Real& theCoeff, + Standard_Real& theAngle); // static Standard_Boolean ExploreCurve(const gp_Cylinder& aCy, @@ -33,10 +36,6 @@ static IntAna_Curve& aC, const Standard_Real aTol, IntAna_ListOfCurve& aLC); -static - Standard_Boolean IsToReverse(const gp_Cylinder& Cy1, - const gp_Cylinder& Cy2, - const Standard_Real Tol); static Standard_Boolean InscribePoint(const Standard_Real theUfTarget, const Standard_Real theUlTarget, @@ -45,16 +44,399 @@ static Standard_Boolean InscribePoint(const Standard_Real theUfTarget, const Standard_Real thePeriod, const Standard_Boolean theFlForce); + +class ComputationMethods +{ +public: + //Stores equations coefficients + struct stCoeffsValue + { + stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&); + + math_Vector mVecA1; + math_Vector mVecA2; + math_Vector mVecB1; + math_Vector mVecB2; + math_Vector mVecC1; + math_Vector mVecC2; + math_Vector mVecD; + + Standard_Real mK21; //sinU2 + Standard_Real mK11; //sinU1 + Standard_Real mL21; //cosU2 + Standard_Real mL11; //cosU1 + Standard_Real mM1; //Free member + + Standard_Real mK22; //sinU2 + Standard_Real mK12; //sinU1 + Standard_Real mL22; //cosU2 + Standard_Real mL12; //cosU1 + Standard_Real mM2; //Free member + + Standard_Real mK1; + Standard_Real mL1; + Standard_Real mK2; + Standard_Real mL2; + + Standard_Real mFIV1; + Standard_Real mPSIV1; + Standard_Real mFIV2; + Standard_Real mPSIV2; + + Standard_Real mB; + Standard_Real mC; + Standard_Real mFI1; + Standard_Real mFI2; + }; + + + //! Determines, if U2(U1) function is increasing. + static Standard_Boolean CylCylMonotonicity(const Standard_Real theU1par, + const Standard_Integer theWLIndex, + const stCoeffsValue& theCoeffs, + const Standard_Real thePeriod, + Standard_Boolean& theIsIncreasing); + + //! Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0, + //! esimates the tolerance of U2-computing (estimation result is + //! assigned to *theDelta value). + static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par, + const Standard_Integer theWLIndex, + const stCoeffsValue& theCoeffs, + Standard_Real& theU2, + Standard_Real* const theDelta = 0); + + static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1, + const Standard_Real theU2, + const stCoeffsValue& theCoeffs, + Standard_Real& theV1, + Standard_Real& theV2); + + static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par, + const Standard_Integer theWLIndex, + const stCoeffsValue& theCoeffs, + Standard_Real& theU2, + Standard_Real& theV1, + Standard_Real& theV2); + +}; + +ComputationMethods::stCoeffsValue::stCoeffsValue(const gp_Cylinder& theCyl1, + const gp_Cylinder& theCyl2): + mVecA1(-theCyl1.Radius()*theCyl1.XAxis().Direction().XYZ()), + mVecA2(theCyl2.Radius()*theCyl2.XAxis().Direction().XYZ()), + mVecB1(-theCyl1.Radius()*theCyl1.YAxis().Direction().XYZ()), + mVecB2(theCyl2.Radius()*theCyl2.YAxis().Direction().XYZ()), + mVecC1(theCyl1.Axis().Direction().XYZ()), + mVecC2(theCyl2.Axis().Direction().XYZ().Reversed()), + mVecD(theCyl2.Location().XYZ() - theCyl1.Location().XYZ()) +{ + enum CoupleOfEquation + { + COENONE = 0, + COE12 = 1, + COE23 = 2, + COE13 = 3 + }aFoundCouple = COENONE; + + + Standard_Real aDetV1V2 = 0.0; + + const Standard_Real aDelta1 = mVecC1(1)*mVecC2(2)-mVecC1(2)*mVecC2(1); //1-2 + const Standard_Real aDelta2 = mVecC1(2)*mVecC2(3)-mVecC1(3)*mVecC2(2); //2-3 + const Standard_Real aDelta3 = mVecC1(1)*mVecC2(3)-mVecC1(3)*mVecC2(1); //1-3 + const Standard_Real anAbsD1 = Abs(aDelta1); //1-2 + const Standard_Real anAbsD2 = Abs(aDelta2); //2-3 + const Standard_Real anAbsD3 = Abs(aDelta3); //1-3 + + if(anAbsD1 >= anAbsD2) + { + if(anAbsD3 > anAbsD1) + { + aFoundCouple = COE13; + aDetV1V2 = aDelta3; + } + else + { + aFoundCouple = COE12; + aDetV1V2 = aDelta1; + } + } + else + { + if(anAbsD3 > anAbsD2) + { + aFoundCouple = COE13; + aDetV1V2 = aDelta3; + } + else + { + aFoundCouple = COE23; + aDetV1V2 = aDelta2; + } + } + + // In point of fact, every determinant (aDelta1, aDelta2 and aDelta3) is + // cross-product between directions (i.e. sine of angle). + // If sine is too small then sine is (approx.) equal to angle itself. + // Therefore, in this case we should compare sine with angular tolerance. + // This constant is used for check if axes are parallel (see constructor + // AxeOperator::AxeOperator(...) in IntAna_QuadQuadGeo.cxx file). + if(Abs(aDetV1V2) < Precision::Angular()) + { + Standard_Failure::Raise("Error. Exception in divide by zerro (IntCyCyTrim)!!!!"); + } + + switch(aFoundCouple) + { + case COE12: + break; + case COE23: + { + math_Vector aVTemp(mVecA1); + mVecA1(1) = aVTemp(2); + mVecA1(2) = aVTemp(3); + mVecA1(3) = aVTemp(1); + + aVTemp = mVecA2; + mVecA2(1) = aVTemp(2); + mVecA2(2) = aVTemp(3); + mVecA2(3) = aVTemp(1); + + aVTemp = mVecB1; + mVecB1(1) = aVTemp(2); + mVecB1(2) = aVTemp(3); + mVecB1(3) = aVTemp(1); + + aVTemp = mVecB2; + mVecB2(1) = aVTemp(2); + mVecB2(2) = aVTemp(3); + mVecB2(3) = aVTemp(1); + + aVTemp = mVecC1; + mVecC1(1) = aVTemp(2); + mVecC1(2) = aVTemp(3); + mVecC1(3) = aVTemp(1); + + aVTemp = mVecC2; + mVecC2(1) = aVTemp(2); + mVecC2(2) = aVTemp(3); + mVecC2(3) = aVTemp(1); + + aVTemp = mVecD; + mVecD(1) = aVTemp(2); + mVecD(2) = aVTemp(3); + mVecD(3) = aVTemp(1); + + break; + } + case COE13: + { + math_Vector aVTemp = mVecA1; + mVecA1(2) = aVTemp(3); + mVecA1(3) = aVTemp(2); + + aVTemp = mVecA2; + mVecA2(2) = aVTemp(3); + mVecA2(3) = aVTemp(2); + + aVTemp = mVecB1; + mVecB1(2) = aVTemp(3); + mVecB1(3) = aVTemp(2); + + aVTemp = mVecB2; + mVecB2(2) = aVTemp(3); + mVecB2(3) = aVTemp(2); + + aVTemp = mVecC1; + mVecC1(2) = aVTemp(3); + mVecC1(3) = aVTemp(2); + + aVTemp = mVecC2; + mVecC2(2) = aVTemp(3); + mVecC2(3) = aVTemp(2); + + aVTemp = mVecD; + mVecD(2) = aVTemp(3); + mVecD(3) = aVTemp(2); + + break; + } + default: + break; + } + + //------- For V1 (begin) + //sinU2 + mK21 = (mVecC2(2)*mVecB2(1)-mVecC2(1)*mVecB2(2))/aDetV1V2; + //sinU1 + mK11 = (mVecC2(2)*mVecB1(1)-mVecC2(1)*mVecB1(2))/aDetV1V2; + //cosU2 + mL21 = (mVecC2(2)*mVecA2(1)-mVecC2(1)*mVecA2(2))/aDetV1V2; + //cosU1 + mL11 = (mVecC2(2)*mVecA1(1)-mVecC2(1)*mVecA1(2))/aDetV1V2; + //Free member + mM1 = (mVecC2(2)*mVecD(1)-mVecC2(1)*mVecD(2))/aDetV1V2; + //------- For V1 (end) + + //------- For V2 (begin) + //sinU2 + mK22 = (mVecC1(1)*mVecB2(2)-mVecC1(2)*mVecB2(1))/aDetV1V2; + //sinU1 + mK12 = (mVecC1(1)*mVecB1(2)-mVecC1(2)*mVecB1(1))/aDetV1V2; + //cosU2 + mL22 = (mVecC1(1)*mVecA2(2)-mVecC1(2)*mVecA2(1))/aDetV1V2; + //cosU1 + mL12 = (mVecC1(1)*mVecA1(2)-mVecC1(2)*mVecA1(1))/aDetV1V2; + //Free member + mM2 = (mVecC1(1)*mVecD(2)-mVecC1(2)*mVecD(1))/aDetV1V2; + //------- For V1 (end) + + ShortCosForm(mL11, mK11, mK1, mFIV1); + ShortCosForm(mL21, mK21, mL1, mPSIV1); + ShortCosForm(mL12, mK12, mK2, mFIV2); + ShortCosForm(mL22, mK22, mL2, mPSIV2); + + const Standard_Real aA1=mVecC1(3)*mK21+mVecC2(3)*mK22-mVecB2(3), //sinU2 + aA2=mVecC1(3)*mL21+mVecC2(3)*mL22-mVecA2(3), //cosU2 + aB1=mVecB1(3)-mVecC1(3)*mK11-mVecC2(3)*mK12, //sinU1 + aB2=mVecA1(3)-mVecC1(3)*mL11-mVecC2(3)*mL12; //cosU1 + + mC =mVecD(3) - mVecC1(3)*mM1 -mVecC2(3)*mM2; //Free + + Standard_Real aA = 0.0; + + ShortCosForm(aB2,aB1,mB,mFI1); + ShortCosForm(aA2,aA1,aA,mFI2); + + mB /= aA; + mC /= aA; +} + +class WorkWithBoundaries +{ +public: + enum SearchBoundType + { + SearchNONE = 0, + SearchV1 = 1, + SearchV2 = 2 + }; + + struct StPInfo + { + StPInfo() + { + mySurfID = 0; + myU1 = RealLast(); + myV1 = RealLast(); + myU2 = RealLast(); + myV2 = RealLast(); + } + + //Equal to 0 for 1st surface non-zerro for 2nd one. + Standard_Integer mySurfID; + + Standard_Real myU1; + Standard_Real myV1; + Standard_Real myU2; + Standard_Real myV2; + + bool operator>(const StPInfo& theOther) const + { + return myU1 > theOther.myU1; + } + + bool operator<(const StPInfo& theOther) const + { + return myU1 < theOther.myU1; + } + + bool operator==(const StPInfo& theOther) const + { + return myU1 == theOther.myU1; + } + }; + + WorkWithBoundaries(const IntSurf_Quadric& theQuad1, + const IntSurf_Quadric& theQuad2, + const ComputationMethods::stCoeffsValue& theCoeffs, + const Bnd_Box2d& theUVSurf1, + const Bnd_Box2d& theUVSurf2, + const Standard_Integer theNbWLines, + const Standard_Real thePeriod, + const Standard_Real theTol3D, + const Standard_Real theTol2D, + const Standard_Boolean isTheReverse) : + myQuad1(theQuad1), myQuad2(theQuad2), myCoeffs(theCoeffs), + myUVSurf1(theUVSurf1), myUVSurf2(theUVSurf2), myNbWLines(theNbWLines), + myPeriod(thePeriod), myTol3D(theTol3D), myTol2D(theTol2D), + myIsReverse(isTheReverse) + { + }; + + void AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL, + const Standard_Real theU1, + const Standard_Real theU2, + const Standard_Real theV1, + const Standard_Real theV1Prev, + const Standard_Real theV2, + const Standard_Real theV2Prev, + const Standard_Integer theWLIndex, + const Standard_Boolean theFlForce, + Standard_Boolean& isTheFound1, + Standard_Boolean& isTheFound2) const; + + Standard_Boolean BoundariesComputing(Standard_Real theU1f[], + Standard_Real theU1l[]) const; + + void BoundaryEstimation(const gp_Cylinder& theCy1, + const gp_Cylinder& theCy2, + Bnd_Range& theOutBoxS1, + Bnd_Range& theOutBoxS2) const; + +protected: + + Standard_Boolean + SearchOnVBounds(const SearchBoundType theSBType, + const Standard_Real theVzad, + const Standard_Real theVInit, + const Standard_Real theInitU2, + const Standard_Real theInitMainVar, + Standard_Real& theMainVariableValue) const; + + const WorkWithBoundaries& operator=(const WorkWithBoundaries&); + +private: + friend class ComputationMethods; + + const IntSurf_Quadric& myQuad1; + const IntSurf_Quadric& myQuad2; + const ComputationMethods::stCoeffsValue& myCoeffs; + const Bnd_Box2d& myUVSurf1; + const Bnd_Box2d& myUVSurf2; + const Standard_Integer myNbWLines; + const Standard_Real myPeriod; + const Standard_Real myTol3D; + const Standard_Real myTol2D; + const Standard_Boolean myIsReverse; +}; + +static + Standard_Boolean ExploreCurve(const gp_Cylinder& aCy, + const gp_Cone& aCo, + IntAna_Curve& aC, + const Standard_Real aTol, + IntAna_ListOfCurve& aLC); + static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, const Handle(IntSurf_LineOn2S)& theLine, - const stCoeffsValue& theCoeffs, + const ComputationMethods::stCoeffsValue& theCoeffs, const Standard_Integer theWLIndex, const Standard_Integer theMinNbPoints, const Standard_Integer theStartPointOnLine, const Standard_Integer theEndPointOnLine, - const Standard_Real theU2f, - const Standard_Real theU2l, const Standard_Real theTol2D, const Standard_Real thePeriodOfSurf2, const Standard_Boolean isTheReverse); @@ -74,6 +456,29 @@ static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX) } } +//======================================================================= +//function : ExtremaLineLine +//purpose : Computes extrema between the given lines. Returns parameters +// on correspond curve. +//======================================================================= +static inline void ExtremaLineLine(const gp_Ax1& theC1, + const gp_Ax1& theC2, + const Standard_Real theCosA, + const Standard_Real theSqSinA, + Standard_Real& thePar1, + Standard_Real& thePar2) +{ + const gp_Dir &aD1 = theC1.Direction(), + &aD2 = theC2.Direction(); + + const gp_XYZ aL1L2 = theC2.Location().XYZ() - theC1.Location().XYZ(); + const Standard_Real aD1L = aD1.XYZ().Dot(aL1L2), + aD2L = aD2.XYZ().Dot(aL1L2); + + thePar1 = (aD1L - theCosA * aD2L) / theSqSinA; + thePar2 = (theCosA * aD1L - aD2L) / theSqSinA; +} + //======================================================================= //function : VBoundaryPrecise //purpose : By default, we shall consider, that V1 and V2 will increase @@ -165,7 +570,7 @@ static Standard_Boolean StepComputing(const math_Matrix& theMatr, Standard_Real& theV1Found, Standard_Real& theV2Found*/) { -#ifdef OCCT_DEBUG +#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG bool flShow = false; if(flShow) @@ -426,22 +831,24 @@ void ProcessBounds(const Handle(IntPatch_ALine)& alig, //-- ligne coura alig->SetLastPoint(alig->NbVertex()); } } + //======================================================================= -//function : IntCyCy +//function : CyCyAnalyticalIntersect //purpose : //======================================================================= -Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, - const IntSurf_Quadric& Quad2, - const Standard_Real Tol, - Standard_Boolean& Empty, - Standard_Boolean& Same, - Standard_Boolean& Multpoint, - IntPatch_SequenceOfLine& slin, - IntPatch_SequenceOfPoint& spnt) - +Standard_Boolean CyCyAnalyticalIntersect( const IntSurf_Quadric& Quad1, + const IntSurf_Quadric& Quad2, + const IntAna_QuadQuadGeo& theInter, + const Standard_Real Tol, + const Standard_Boolean isTheReverse, + Standard_Boolean& Empty, + Standard_Boolean& Same, + Standard_Boolean& Multpoint, + IntPatch_SequenceOfLine& slin, + IntPatch_SequenceOfPoint& spnt) { IntPatch_Point ptsol; - + Standard_Integer i; IntSurf_TypeTrans trans1,trans2; @@ -453,39 +860,43 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, gp_Cylinder Cy1(Quad1.Cylinder()); gp_Cylinder Cy2(Quad2.Cylinder()); - IntAna_QuadQuadGeo inter(Cy1,Cy2,Tol); - - if (!inter.IsDone()) - { - return Standard_False; - } - - typint = inter.TypeInter(); - Standard_Integer NbSol = inter.NbSolutions(); + typint = theInter.TypeInter(); + Standard_Integer NbSol = theInter.NbSolutions(); Empty = Standard_False; Same = Standard_False; switch (typint) - { + { case IntAna_Empty: { Empty = Standard_True; - } + } break; case IntAna_Same: - { + { Same = Standard_True; - } + } break; - + case IntAna_Point: { - gp_Pnt psol(inter.Point(1)); - Standard_Real U1,V1,U2,V2; - Quad1.Parameters(psol,U1,V1); - Quad2.Parameters(psol,U2,V2); + gp_Pnt psol(theInter.Point(1)); ptsol.SetValue(psol,Tol,Standard_True); + + Standard_Real U1,V1,U2,V2; + + if(isTheReverse) + { + Quad2.Parameters(psol, U1, V1); + Quad1.Parameters(psol, U2, V2); + } + else + { + Quad1.Parameters(psol, U1, V1); + Quad2.Parameters(psol, U2, V2); + } + ptsol.SetParameters(U1,V1,U2,V2); spnt.Append(ptsol); } @@ -496,10 +907,14 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, gp_Pnt ptref; if (NbSol == 1) { // Cylinders are tangent to each other by line - linsol = inter.Line(1); + linsol = theInter.Line(1); ptref = linsol.Location(); + + //Radius-vectors gp_Dir crb1(gp_Vec(ptref,Cy1.Location())); gp_Dir crb2(gp_Vec(ptref,Cy2.Location())); + + //outer normal lines gp_Vec norm1(Quad1.Normale(ptref)); gp_Vec norm2(Quad2.Normale(ptref)); IntSurf_Situation situcyl1; @@ -507,6 +922,11 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, if (crb1.Dot(crb2) < 0.) { // centre de courbures "opposes" + //ATTENTION!!! + // Normal and Radius-vector of the 1st(!) cylinder + // is used for judging what the situation of the 2nd(!) + // cylinder is. + if (norm1.Dot(crb1) > 0.) { situcyl2 = IntSurf_Inside; @@ -526,7 +946,7 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, } } else - { + { if (Cy1.Radius() < Cy2.Radius()) { if (norm1.Dot(crb1) > 0.) @@ -576,9 +996,11 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, { for (i=1; i <= NbSol; i++) { - linsol = inter.Line(i); + linsol = theInter.Line(i); ptref = linsol.Location(); gp_Vec lsd = linsol.Direction(); + + //Theoretically, qwe = +/- 1.0. Standard_Real qwe = lsd.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)); if (qwe >0.00000001) { @@ -607,31 +1029,52 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, gp_Vec Tgt; gp_Pnt ptref; IntPatch_Point pmult1, pmult2; - - elipsol = inter.Ellipse(1); - + + elipsol = theInter.Ellipse(1); + gp_Pnt pttang1(ElCLib::Value(0.5*M_PI, elipsol)); gp_Pnt pttang2(ElCLib::Value(1.5*M_PI, elipsol)); - + Multpoint = Standard_True; pmult1.SetValue(pttang1,Tol,Standard_True); pmult2.SetValue(pttang2,Tol,Standard_True); pmult1.SetMultiple(Standard_True); pmult2.SetMultiple(Standard_True); - + Standard_Real oU1,oV1,oU2,oV2; - Quad1.Parameters(pttang1,oU1,oV1); - Quad2.Parameters(pttang1,oU2,oV2); + + if(isTheReverse) + { + Quad2.Parameters(pttang1,oU1,oV1); + Quad1.Parameters(pttang1,oU2,oV2); + } + else + { + Quad1.Parameters(pttang1,oU1,oV1); + Quad2.Parameters(pttang1,oU2,oV2); + } + pmult1.SetParameters(oU1,oV1,oU2,oV2); - Quad1.Parameters(pttang2,oU1,oV1); - Quad2.Parameters(pttang2,oU2,oV2); - pmult2.SetParameters(oU1,oV1,oU2,oV2); - - // on traite la premiere ellipse + if(isTheReverse) + { + Quad2.Parameters(pttang2,oU1,oV1); + Quad1.Parameters(pttang2,oU2,oV2); + } + else + { + Quad1.Parameters(pttang2,oU1,oV1); + Quad2.Parameters(pttang2,oU2,oV2); + } + pmult2.SetParameters(oU1,oV1,oU2,oV2); + + // on traite la premiere ellipse + //-- Calcul de la Transition de la ligne ElCLib::D1(0.,elipsol,ptref,Tgt); + + //Theoretically, qwe = +/- |Tgt|. Standard_Real qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)); if (qwe>0.00000001) { @@ -644,7 +1087,7 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, trans2 = IntSurf_Out; } else - { + { trans1=trans2=IntSurf_Undecided; } @@ -660,8 +1103,18 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, aIP.SetValue(aP,Tol,Standard_False); aIP.SetMultiple(Standard_False); // - Quad1.Parameters(aP, aU1, aV1); - Quad2.Parameters(aP, aU2, aV2); + + if(isTheReverse) + { + Quad2.Parameters(aP, aU1, aV1); + Quad1.Parameters(aP, aU2, aV2); + } + else + { + Quad1.Parameters(aP, aU1, aV1); + Quad2.Parameters(aP, aU2, aV2); + } + aIP.SetParameters(aU1, aV1, aU2, aV2); // aIP.SetParameter(0.); @@ -685,7 +1138,7 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, //-- Transitions calculee au point 0 OK // // on traite la deuxieme ellipse - elipsol = inter.Ellipse(2); + elipsol = theInter.Ellipse(2); Standard_Real param1 = ElCLib::Parameter(elipsol,pttang1); Standard_Real param2 = ElCLib::Parameter(elipsol,pttang2); @@ -704,6 +1157,8 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, //-- Calcul des transitions de ligne pour la premiere ligne ElCLib::D1(parampourtransition,elipsol,ptref,Tgt); + + //Theoretically, qwe = +/- |Tgt|. qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)); if(qwe> 0.00000001) { @@ -731,8 +1186,18 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, aIP.SetValue(aP,Tol,Standard_False); aIP.SetMultiple(Standard_False); // - Quad1.Parameters(aP, aU1, aV1); - Quad2.Parameters(aP, aU2, aV2); + + if(isTheReverse) + { + Quad2.Parameters(aP, aU1, aV1); + Quad1.Parameters(aP, aU2, aV2); + } + else + { + Quad1.Parameters(aP, aU1, aV1); + Quad2.Parameters(aP, aU2, aV2); + } + aIP.SetParameters(aU1, aV1, aU2, aV2); // aIP.SetParameter(0.); @@ -751,111 +1216,14 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, } break; + case IntAna_Parabola: + case IntAna_Hyperbola: + Standard_Failure::Raise("IntCyCy(): Wrong intersection type!"); + + case IntAna_Circle: + // Circle is useful when we will work with trimmed surfaces + // (two cylinders can be tangent by their basises, e.g. circle) case IntAna_NoGeometricSolution: - { - Standard_Boolean bReverse; - Standard_Real U1,V1,U2,V2; - gp_Pnt psol; - // - bReverse=IsToReverse(Cy1, Cy2, Tol); - if (bReverse) - { - Cy2=Quad1.Cylinder(); - Cy1=Quad2.Cylinder(); - } - // - IntAna_IntQuadQuad anaint(Cy1,Cy2,Tol); - if (!anaint.IsDone()) - { - return Standard_False; - } - - if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) - { - Empty = Standard_True; - } - else - { - NbSol = anaint.NbPnt(); - for (i = 1; i <= NbSol; i++) - { - psol = anaint.Point(i); - Quad1.Parameters(psol,U1,V1); - Quad2.Parameters(psol,U2,V2); - ptsol.SetValue(psol,Tol,Standard_True); - ptsol.SetParameters(U1,V1,U2,V2); - spnt.Append(ptsol); - } - - gp_Pnt ptvalid, ptf, ptl; - gp_Vec tgvalid; - - Standard_Real first,last,para; - IntAna_Curve curvsol; - Standard_Boolean tgfound; - Standard_Integer kount; - - NbSol = anaint.NbCurve(); - for (i = 1; i <= NbSol; i++) - { - curvsol = anaint.Curve(i); - curvsol.Domain(first,last); - ptf = curvsol.Value(first); - ptl = curvsol.Value(last); - - para = last; - kount = 1; - tgfound = Standard_False; - - while (!tgfound) - { - para = (1.123*first + para)/2.123; - tgfound = curvsol.D1u(para,ptvalid,tgvalid); - if (!tgfound) - { - kount ++; - tgfound = kount > 5; - } - } - - Handle(IntPatch_ALine) alig; - if (kount <= 5) - { - Standard_Real qwe = tgvalid.DotCross( Quad2.Normale(ptvalid), - Quad1.Normale(ptvalid)); - if(qwe>0.00000001) - { - trans1 = IntSurf_Out; - trans2 = IntSurf_In; - } - else if(qwe<-0.00000001) - { - trans1 = IntSurf_In; - trans2 = IntSurf_Out; - } - else - { - trans1=trans2=IntSurf_Undecided; - } - alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2); - } - else - { - alig = new IntPatch_ALine(curvsol,Standard_False); - //-- cout << "Transition indeterminee" << endl; - } - - Standard_Boolean TempFalse1 = Standard_False; - Standard_Boolean TempFalse2 = Standard_False; - - ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1,ptf,first, - TempFalse2,ptl,last,Multpoint,Tol); - slin.Append(alig); - } - } - } - break; - default: return Standard_False; } @@ -919,255 +1287,15 @@ static void ShortCosForm( const Standard_Real theCosFactor, } } -enum SearchBoundType -{ - SearchNONE = 0, - SearchV1 = 1, - SearchV2 = 2 -}; - -//Stores equations coefficients -struct stCoeffsValue -{ - stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&); - - math_Vector mVecA1; - math_Vector mVecA2; - math_Vector mVecB1; - math_Vector mVecB2; - math_Vector mVecC1; - math_Vector mVecC2; - math_Vector mVecD; - - Standard_Real mK21; //sinU2 - Standard_Real mK11; //sinU1 - Standard_Real mL21; //cosU2 - Standard_Real mL11; //cosU1 - Standard_Real mM1; //Free member - - Standard_Real mK22; //sinU2 - Standard_Real mK12; //sinU1 - Standard_Real mL22; //cosU2 - Standard_Real mL12; //cosU1 - Standard_Real mM2; //Free member - - Standard_Real mK1; - Standard_Real mL1; - Standard_Real mK2; - Standard_Real mL2; - - Standard_Real mFIV1; - Standard_Real mPSIV1; - Standard_Real mFIV2; - Standard_Real mPSIV2; - - Standard_Real mB; - Standard_Real mC; - Standard_Real mFI1; - Standard_Real mFI2; -}; - -stCoeffsValue::stCoeffsValue( const gp_Cylinder& theCyl1, - const gp_Cylinder& theCyl2): - mVecA1(-theCyl1.Radius()*theCyl1.XAxis().Direction().XYZ()), - mVecA2(theCyl2.Radius()*theCyl2.XAxis().Direction().XYZ()), - mVecB1(-theCyl1.Radius()*theCyl1.YAxis().Direction().XYZ()), - mVecB2(theCyl2.Radius()*theCyl2.YAxis().Direction().XYZ()), - mVecC1(theCyl1.Axis().Direction().XYZ()), - mVecC2(theCyl2.Axis().Direction().XYZ().Reversed()), - mVecD(theCyl2.Location().XYZ() - theCyl1.Location().XYZ()) -{ - enum CoupleOfEquation - { - COENONE = 0, - COE12 = 1, - COE23 = 2, - COE13 = 3 - }aFoundCouple = COENONE; - - - Standard_Real aDetV1V2 = 0.0; - - const Standard_Real aDelta1 = mVecC1(1)*mVecC2(2)-mVecC1(2)*mVecC2(1); //1-2 - const Standard_Real aDelta2 = mVecC1(2)*mVecC2(3)-mVecC1(3)*mVecC2(2); //2-3 - const Standard_Real aDelta3 = mVecC1(1)*mVecC2(3)-mVecC1(3)*mVecC2(1); //1-3 - const Standard_Real anAbsD1 = Abs(aDelta1); //1-2 - const Standard_Real anAbsD2 = Abs(aDelta2); //2-3 - const Standard_Real anAbsD3 = Abs(aDelta3); //1-3 - - if(anAbsD1 >= anAbsD2) - { - if(anAbsD3 > anAbsD1) - { - aFoundCouple = COE13; - aDetV1V2 = aDelta3; - } - else - { - aFoundCouple = COE12; - aDetV1V2 = aDelta1; - } - } - else - { - if(anAbsD3 > anAbsD2) - { - aFoundCouple = COE13; - aDetV1V2 = aDelta3; - } - else - { - aFoundCouple = COE23; - aDetV1V2 = aDelta2; - } - } - - // In point of fact, every determinant (aDelta1, aDelta2 and aDelta3) is - // cross-product between directions (i.e. sine of angle). - // If sine is too small then sine is (approx.) equal to angle itself. - // Therefore, in this case we should compare sine with angular tolerance. - // This constant is used for check if axes are parallel (see constructor - // AxeOperator::AxeOperator(...) in IntAna_QuadQuadGeo.cxx file). - if(Abs(aDetV1V2) < Precision::Angular()) - { - Standard_Failure::Raise("Error. Exception in divide by zerro (IntCyCyTrim)!!!!"); - } - - switch(aFoundCouple) - { - case COE12: - break; - case COE23: - { - math_Vector aVTemp(mVecA1); - mVecA1(1) = aVTemp(2); - mVecA1(2) = aVTemp(3); - mVecA1(3) = aVTemp(1); - - aVTemp = mVecA2; - mVecA2(1) = aVTemp(2); - mVecA2(2) = aVTemp(3); - mVecA2(3) = aVTemp(1); - - aVTemp = mVecB1; - mVecB1(1) = aVTemp(2); - mVecB1(2) = aVTemp(3); - mVecB1(3) = aVTemp(1); - - aVTemp = mVecB2; - mVecB2(1) = aVTemp(2); - mVecB2(2) = aVTemp(3); - mVecB2(3) = aVTemp(1); - - aVTemp = mVecC1; - mVecC1(1) = aVTemp(2); - mVecC1(2) = aVTemp(3); - mVecC1(3) = aVTemp(1); - - aVTemp = mVecC2; - mVecC2(1) = aVTemp(2); - mVecC2(2) = aVTemp(3); - mVecC2(3) = aVTemp(1); - - aVTemp = mVecD; - mVecD(1) = aVTemp(2); - mVecD(2) = aVTemp(3); - mVecD(3) = aVTemp(1); - - break; - } - case COE13: - { - math_Vector aVTemp = mVecA1; - mVecA1(2) = aVTemp(3); - mVecA1(3) = aVTemp(2); - - aVTemp = mVecA2; - mVecA2(2) = aVTemp(3); - mVecA2(3) = aVTemp(2); - - aVTemp = mVecB1; - mVecB1(2) = aVTemp(3); - mVecB1(3) = aVTemp(2); - - aVTemp = mVecB2; - mVecB2(2) = aVTemp(3); - mVecB2(3) = aVTemp(2); - - aVTemp = mVecC1; - mVecC1(2) = aVTemp(3); - mVecC1(3) = aVTemp(2); - - aVTemp = mVecC2; - mVecC2(2) = aVTemp(3); - mVecC2(3) = aVTemp(2); - - aVTemp = mVecD; - mVecD(2) = aVTemp(3); - mVecD(3) = aVTemp(2); - - break; - } - default: - break; - } - - //------- For V1 (begin) - //sinU2 - mK21 = (mVecC2(2)*mVecB2(1)-mVecC2(1)*mVecB2(2))/aDetV1V2; - //sinU1 - mK11 = (mVecC2(2)*mVecB1(1)-mVecC2(1)*mVecB1(2))/aDetV1V2; - //cosU2 - mL21 = (mVecC2(2)*mVecA2(1)-mVecC2(1)*mVecA2(2))/aDetV1V2; - //cosU1 - mL11 = (mVecC2(2)*mVecA1(1)-mVecC2(1)*mVecA1(2))/aDetV1V2; - //Free member - mM1 = (mVecC2(2)*mVecD(1)-mVecC2(1)*mVecD(2))/aDetV1V2; - //------- For V1 (end) - - //------- For V2 (begin) - //sinU2 - mK22 = (mVecC1(1)*mVecB2(2)-mVecC1(2)*mVecB2(1))/aDetV1V2; - //sinU1 - mK12 = (mVecC1(1)*mVecB1(2)-mVecC1(2)*mVecB1(1))/aDetV1V2; - //cosU2 - mL22 = (mVecC1(1)*mVecA2(2)-mVecC1(2)*mVecA2(1))/aDetV1V2; - //cosU1 - mL12 = (mVecC1(1)*mVecA1(2)-mVecC1(2)*mVecA1(1))/aDetV1V2; - //Free member - mM2 = (mVecC1(1)*mVecD(2)-mVecC1(2)*mVecD(1))/aDetV1V2; - //------- For V1 (end) - - ShortCosForm(mL11, mK11, mK1, mFIV1); - ShortCosForm(mL21, mK21, mL1, mPSIV1); - ShortCosForm(mL12, mK12, mK2, mFIV2); - ShortCosForm(mL22, mK22, mL2, mPSIV2); - - const Standard_Real aA1=mVecC1(3)*mK21+mVecC2(3)*mK22-mVecB2(3), //sinU2 - aA2=mVecC1(3)*mL21+mVecC2(3)*mL22-mVecA2(3), //cosU2 - aB1=mVecB1(3)-mVecC1(3)*mK11-mVecC2(3)*mK12, //sinU1 - aB2=mVecA1(3)-mVecC1(3)*mL11-mVecC2(3)*mL12; //cosU1 - - mC =mVecD(3) - mVecC1(3)*mM1 -mVecC2(3)*mM2; //Free - - Standard_Real aA = 0.0; - - ShortCosForm(aB2,aB1,mB,mFI1); - ShortCosForm(aA2,aA1,aA,mFI2); - - mB /= aA; - mC /= aA; -} - //======================================================================= //function : CylCylMonotonicity //purpose : Determines, if U2(U1) function is increasing. //======================================================================= -static Standard_Boolean CylCylMonotonicity( const Standard_Real theU1par, - const Standard_Integer theWLIndex, - const stCoeffsValue& theCoeffs, - const Standard_Real thePeriod, - Standard_Boolean& theIsIncreasing) +Standard_Boolean ComputationMethods::CylCylMonotonicity(const Standard_Real theU1par, + const Standard_Integer theWLIndex, + const stCoeffsValue& theCoeffs, + const Standard_Real thePeriod, + Standard_Boolean& theIsIncreasing) { // U2(U1) = FI2 + (+/-)acos(B*cos(U1 - FI1) + C); //Accordingly, @@ -1235,11 +1363,11 @@ static Standard_Boolean CylCylMonotonicity( const Standard_Real theU1par, // esimates the tolerance of U2-computing (estimation result is // assigned to *theDelta value). //======================================================================= -static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par, +Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par, const Standard_Integer theWLIndex, const stCoeffsValue& theCoeffs, Standard_Real& theU2, - Standard_Real* const theDelta = 0) + Standard_Real* const theDelta) { //This formula is got from some experience and can be changed. const Standard_Real aTol0 = Min(10.0*Epsilon(1.0)*theCoeffs.mB, aNulValue); @@ -1308,7 +1436,7 @@ static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par, return Standard_True; } -static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1, +Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1, const Standard_Real theU2, const stCoeffsValue& theCoeffs, Standard_Real& theV1, @@ -1328,7 +1456,7 @@ static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1, } -static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par, +Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par, const Standard_Integer theWLIndex, const stCoeffsValue& theCoeffs, Standard_Real& theU2, @@ -1348,13 +1476,13 @@ static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par, //function : SearchOnVBounds //purpose : //======================================================================= -static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType, - const stCoeffsValue& theCoeffs, +Standard_Boolean WorkWithBoundaries:: + SearchOnVBounds(const SearchBoundType theSBType, const Standard_Real theVzad, const Standard_Real theVInit, const Standard_Real theInitU2, const Standard_Real theInitMainVar, - Standard_Real& theMainVariableValue) + Standard_Real& theMainVariableValue) const { const Standard_Integer aNbDim = 3; const Standard_Real aMaxError = 4.0*M_PI; // two periods @@ -1381,40 +1509,40 @@ static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType, aSinU2 = sin(aU2Prev), aCosU2 = cos(aU2Prev); - math_Vector aVecFreeMem = (theCoeffs.mVecA2 * aU2Prev + - theCoeffs.mVecB2) * aSinU2 - - (theCoeffs.mVecB2 * aU2Prev - - theCoeffs.mVecA2) * aCosU2 + - (theCoeffs.mVecA1 * aMainVarPrev + - theCoeffs.mVecB1) * aSinU1 - - (theCoeffs.mVecB1 * aMainVarPrev - - theCoeffs.mVecA1) * aCosU1 + - theCoeffs.mVecD; + math_Vector aVecFreeMem = (myCoeffs.mVecA2 * aU2Prev + + myCoeffs.mVecB2) * aSinU2 - + (myCoeffs.mVecB2 * aU2Prev - + myCoeffs.mVecA2) * aCosU2 + + (myCoeffs.mVecA1 * aMainVarPrev + + myCoeffs.mVecB1) * aSinU1 - + (myCoeffs.mVecB1 * aMainVarPrev - + myCoeffs.mVecA1) * aCosU1 + + myCoeffs.mVecD; math_Vector aMSum(1, 3); switch(theSBType) { case SearchV1: - aMatr.SetCol(3, theCoeffs.mVecC2); - aMSum = theCoeffs.mVecC1 * theVzad; + aMatr.SetCol(3, myCoeffs.mVecC2); + aMSum = myCoeffs.mVecC1 * theVzad; aVecFreeMem -= aMSum; - aMSum += theCoeffs.mVecC2*anOtherVar; + aMSum += myCoeffs.mVecC2*anOtherVar; break; case SearchV2: - aMatr.SetCol(3, theCoeffs.mVecC1); - aMSum = theCoeffs.mVecC2 * theVzad; + aMatr.SetCol(3, myCoeffs.mVecC1); + aMSum = myCoeffs.mVecC2 * theVzad; aVecFreeMem -= aMSum; - aMSum += theCoeffs.mVecC1*anOtherVar; + aMSum += myCoeffs.mVecC1*anOtherVar; break; default: return Standard_False; } - aMatr.SetCol(1, theCoeffs.mVecA1 * aSinU1 - theCoeffs.mVecB1 * aCosU1); - aMatr.SetCol(2, theCoeffs.mVecA2 * aSinU2 - theCoeffs.mVecB2 * aCosU2); + aMatr.SetCol(1, myCoeffs.mVecA1 * aSinU1 - myCoeffs.mVecB1 * aCosU1); + aMatr.SetCol(2, myCoeffs.mVecA2 * aSinU2 - myCoeffs.mVecB2 * aCosU2); Standard_Real aDetMainSyst = aMatr.Determinant(); @@ -1460,11 +1588,11 @@ static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType, aCosU1Last = cos(aMainVarPrev), aSinU2Last = sin(aU2Prev), aCosU2Last = cos(aU2Prev); - aMSum -= (theCoeffs.mVecA1*aCosU1Last + - theCoeffs.mVecB1*aSinU1Last + - theCoeffs.mVecA2*aCosU2Last + - theCoeffs.mVecB2*aSinU2Last + - theCoeffs.mVecD); + aMSum -= (myCoeffs.mVecA1*aCosU1Last + + myCoeffs.mVecB1*aSinU1Last + + myCoeffs.mVecA2*aCosU2Last + + myCoeffs.mVecB2*aSinU2Last + + myCoeffs.mVecD); const Standard_Real aSQNorm = aMSum.Norm2(); return (aSQNorm < aTol2); } @@ -1531,27 +1659,10 @@ Standard_Boolean InscribePoint( const Standard_Real theUfTarget, return Standard_True; } - if(IsEqual(thePeriod, 0.0)) - {//it cannot be inscribed - return Standard_False; - } + const Standard_Real aUf = theUfTarget - theTol2D; + const Standard_Real aUl = aUf + thePeriod; - //Make theUGiven nearer to theUfTarget (in order to avoid - //excess iterations) - theUGiven += thePeriod*Floor((theUfTarget-theUGiven)/thePeriod); - Standard_Real aDU = theUGiven - theUfTarget; - - if(aDU > 0.0) - aDU = -thePeriod; - else - aDU = +thePeriod; - - while(((theUGiven - theUfTarget)*aDU < 0.0) && - !((theUfTarget - theUGiven <= theTol2D) && - (theUGiven - theUlTarget <= theTol2D))) - { - theUGiven += aDU; - } + theUGiven = ElCLib::InPeriod(theUGiven, aUf, aUl); return ((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D)); @@ -1625,7 +1736,7 @@ static Standard_Boolean ExcludeNearElements(Standard_Real theArr[], if((anA-anB) < theTol) { if((anB != 0.0) && (anB != theUSurf1f) && (anB != theUSurf1l)) - anA = (anA + anB)/2.0; + anA = (anA + anB)/2.0; else anA = anB; @@ -1645,7 +1756,7 @@ static Standard_Boolean ExcludeNearElements(Standard_Real theArr[], //======================================================================= static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, - const stCoeffsValue& theCoeffs, + const ComputationMethods::stCoeffsValue& theCoeffs, const Standard_Boolean isTheReverse, const Standard_Boolean isThePrecise, const gp_Pnt2d& thePntOnSurf1, @@ -1713,8 +1824,21 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1, aPlast.ParametersOnS1(aUl, aVl); if(!theFlBefore && (aU1par <= aUl)) - {//Parameter value must be increased if theFlBefore == FALSE. - return Standard_False; + { + //Parameter value must be increased if theFlBefore == FALSE. + + aU1par += thePeriodOfSurf1; + + //The condition is as same as in + //InscribePoint(...) function + if((theUfSurf1 - aU1par > theTol2D) || + (aU1par - theUlSurf1 > theTol2D)) + { + //New aU1par is out of target interval. + //Go back to old value. + + return Standard_False; + } } //theTol2D is minimal step along parameter changed. @@ -1760,8 +1884,7 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1, { SeekAdditionalPoints( theQuad1, theQuad2, theLine, theCoeffs, theWLIndex, aDeltaStep, aNbPnts-2, - aNbPnts-1, theUfSurf2, theUlSurf2, - theTol2D, thePeriodOfSurf1, isTheReverse); + aNbPnts-1, theTol2D, thePeriodOfSurf1, isTheReverse); } } @@ -1772,15 +1895,7 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1, //function : AddBoundaryPoint //purpose : //======================================================================= -static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1, - const IntSurf_Quadric& theQuad2, - const Handle(IntPatch_WLine)& theWL, - const stCoeffsValue& theCoeffs, - const Bnd_Box2d& theUVSurf1, - const Bnd_Box2d& theUVSurf2, - const Standard_Real theTol3D, - const Standard_Real theTol2D, - const Standard_Real thePeriod, +void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL, const Standard_Real theU1, const Standard_Real theU2, const Standard_Real theV1, @@ -1788,10 +1903,9 @@ static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1, const Standard_Real theV2, const Standard_Real theV2Prev, const Standard_Integer theWLIndex, - const Standard_Boolean isTheReverse, const Standard_Boolean theFlForce, Standard_Boolean& isTheFound1, - Standard_Boolean& isTheFound2) + Standard_Boolean& isTheFound2) const { Standard_Real aUSurf1f = 0.0, //const aUSurf1l = 0.0, @@ -1802,176 +1916,97 @@ static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1, aVSurf2f = 0.0, aVSurf2l = 0.0; - theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l); - theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l); + myUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l); + myUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l); + + const Standard_Integer aSize = 4; + const Standard_Real anArrVzad[aSize] = {aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l}; - SearchBoundType aTS1 = SearchNONE, aTS2 = SearchNONE; - Standard_Real aV1zad = aVSurf1f, aV2zad = aVSurf2f; + StPInfo aUVPoint[aSize]; - Standard_Real anUpar1 = theU1, anUpar2 = theU1; - Standard_Real aVf = theV1, aVl = theV1Prev; - - if( (Abs(aVf-aVSurf1f) <= theTol2D) || - ((aVf-aVSurf1f)*(aVl-aVSurf1f) <= 0.0)) + for(Standard_Integer anIDSurf = 0; anIDSurf < 4; anIDSurf+=2) { - aTS1 = SearchV1; - aV1zad = aVSurf1f; - isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1f, theV2, theU2, theU1, anUpar1); - } - else if((Abs(aVf-aVSurf1l) <= theTol2D) || - ((aVf-aVSurf1l)*(aVl-aVSurf1l) <= 0.0)) - { - aTS1 = SearchV1; - aV1zad = aVSurf1l; - isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1l, theV2, theU2, theU1, anUpar1); - } + const Standard_Real aVf = (anIDSurf == 0) ? theV1 : theV2, + aVl = (anIDSurf == 0) ? theV1Prev : theV2Prev; - aVf = theV2; - aVl = theV2Prev; + const SearchBoundType aTS = (anIDSurf == 0) ? SearchV1 : SearchV2; - if((Abs(aVf-aVSurf2f) <= theTol2D) || - ((aVf-aVSurf2f)*(aVl-aVSurf2f) <= 0.0)) - { - aTS2 = SearchV2; - aV2zad = aVSurf2f; - isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2f, theV1, theU2, theU1, anUpar2); - } - else if((Abs(aVf-aVSurf2l) <= theTol2D) || - ((aVf-aVSurf2l)*(aVl-aVSurf2l) <= 0.0)) - { - aTS2 = SearchV2; - aV2zad = aVSurf2l; - isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2l, theV1, theU2, theU1, anUpar2); - } - - if(!isTheFound1 && !isTheFound2) - return Standard_True; - - //We increase U1 parameter. Therefore, added point must have U1 parameter less than - //or equal to current (conditions "(anUpar1 < theU1)" and "(anUpar2 < theU1)"). - - if(anUpar1 < anUpar2) - { - if(isTheFound1 && (anUpar1 < theU1)) + for(Standard_Integer anIDBound = 0; anIDBound < 2; anIDBound++) { - Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2; - if(!CylCylComputeParameters(anUpar1, theWLIndex, theCoeffs, aU2, aV1, aV2)) - { - isTheFound1 = Standard_False; - } - - if(aTS1 == SearchV1) - aV1 = aV1zad; - if(aTS1 == SearchV2) - aV2 = aV2zad; + const Standard_Integer anIndex = anIDSurf+anIDBound; - if(isTheFound1 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False, - gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2), - aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l, - aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod, - theWL->Curve(), theWLIndex, theTol3D, - theTol2D, theFlForce)) + aUVPoint[anIndex].mySurfID = anIDSurf; + + if((Abs(aVf-anArrVzad[anIndex]) > myTol2D) && + ((aVf-anArrVzad[anIndex])*(aVl-anArrVzad[anIndex]) > 0.0)) { - isTheFound1 = Standard_False; + continue; } + + const Standard_Boolean aRes = SearchOnVBounds(aTS, anArrVzad[anIndex], + (anIDSurf == 0) ? theV2 : theV1, + theU2, theU1, + aUVPoint[anIndex].myU1); + + if(!aRes || aUVPoint[anIndex].myU1 >= theU1) + { + aUVPoint[anIndex].myU1 = RealLast(); + continue; + } + else + { + Standard_Real &aU1 = aUVPoint[anIndex].myU1, + &aU2 = aUVPoint[anIndex].myU2, + &aV1 = aUVPoint[anIndex].myV1, + &aV2 = aUVPoint[anIndex].myV2; + aU2 = theU2; + aV1 = theV1; + aV2 = theV2; + + if(!ComputationMethods:: + CylCylComputeParameters(aU1, theWLIndex, myCoeffs, aU2, aV1, aV2)) + { + aU1 = RealLast(); + continue; + } + + if(aTS == SearchV1) + aV1 = anArrVzad[anIndex]; + else //if(aTS[anIndex] == SearchV2) + aV2 = anArrVzad[anIndex]; + } + } + } + + std::sort(aUVPoint, aUVPoint+aSize); + + isTheFound1 = isTheFound2 = Standard_False; + + for(Standard_Integer i = 0; i < aSize; i++) + { + if(aUVPoint[i].myU1 == RealLast()) + break; + + if(!AddPointIntoWL(myQuad1, myQuad2, myCoeffs, myIsReverse, Standard_False, + gp_Pnt2d(aUVPoint[i].myU1, aUVPoint[i].myV1), + gp_Pnt2d(aUVPoint[i].myU2, aUVPoint[i].myV2), + aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l, + aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, myPeriod, + theWL->Curve(), theWLIndex, myTol3D, myTol2D, theFlForce)) + { + continue; + } + + if(aUVPoint[i].mySurfID == 0) + { + isTheFound1 = Standard_True; } else { - isTheFound1 = Standard_False; - } - - if(isTheFound2 && (anUpar2 < theU1)) - { - Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2; - if(!CylCylComputeParameters(anUpar2, theWLIndex, theCoeffs, aU2, aV1, aV2)) - { - isTheFound2 = Standard_False; - } - - if(aTS2 == SearchV1) - aV1 = aV1zad; - - if(aTS2 == SearchV2) - aV2 = aV2zad; - - if(isTheFound2 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False, - gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2), - aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l, - aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod, - theWL->Curve(), theWLIndex, theTol3D, - theTol2D, theFlForce)) - { - isTheFound2 = Standard_False; - } - } - else - { - isTheFound2 = Standard_False; + isTheFound2 = Standard_True; } } - else - { - if(isTheFound2 && (anUpar2 < theU1)) - { - Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2; - if(!CylCylComputeParameters(anUpar2, theWLIndex, theCoeffs, aU2, aV1, aV2)) - { - isTheFound2 = Standard_False; - } - - if(aTS2 == SearchV1) - aV1 = aV1zad; - - if(aTS2 == SearchV2) - aV2 = aV2zad; - - if(isTheFound2 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False, - gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2), - aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l, - aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod, - theWL->Curve(), theWLIndex, theTol3D, - theTol2D, theFlForce)) - { - isTheFound2 = Standard_False; - } - } - else - { - isTheFound2 = Standard_False; - } - - if(isTheFound1 && (anUpar1 < theU1)) - { - Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2; - if(!CylCylComputeParameters(anUpar1, theWLIndex, theCoeffs, aU2, aV1, aV2)) - { - isTheFound1 = Standard_False; - } - - if(aTS1 == SearchV1) - aV1 = aV1zad; - - if(aTS1 == SearchV2) - aV2 = aV2zad; - - if(isTheFound1 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False, - gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2), - aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l, - aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod, - theWL->Curve(), theWLIndex, theTol3D, - theTol2D, theFlForce)) - { - isTheFound1 = Standard_False; - } - } - else - { - isTheFound1 = Standard_False; - } - } - - return Standard_True; } //======================================================================= @@ -1981,13 +2016,11 @@ static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1, static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, const Handle(IntSurf_LineOn2S)& theLine, - const stCoeffsValue& theCoeffs, + const ComputationMethods::stCoeffsValue& theCoeffs, const Standard_Integer theWLIndex, const Standard_Integer theMinNbPoints, const Standard_Integer theStartPointOnLine, const Standard_Integer theEndPointOnLine, - const Standard_Real theU2f, - const Standard_Real theU2l, const Standard_Real theTol2D, const Standard_Real thePeriodOfSurf2, const Standard_Boolean isTheReverse) @@ -2062,16 +2095,23 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1, U1prec = 0.5*(U1f+U1l); - if(!CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec)) + if(!ComputationMethods:: + CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec)) + { continue; + } - InscribePoint(theU2f, theU2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False); + MinMax(U2f, U2l); + if(!InscribePoint(U2f, U2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False)) + { + continue; + } const gp_Pnt aP1(theQuad1.Value(U1prec, V1prec)), aP2(theQuad2.Value(U2prec, V2prec)); const gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ())); -#ifdef OCCT_DEBUG - //cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << endl; +#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG + std::cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << std::endl; #endif IntSurf_PntOn2S anIP; @@ -2102,27 +2142,24 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1, //purpose : Computes true domain of future intersection curve (allows // avoiding excess iterations) //======================================================================= -//======================================================================= -static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, - const Standard_Real thePeriod, - Standard_Real theU1f[], - Standard_Real theU1l[]) +Standard_Boolean WorkWithBoundaries::BoundariesComputing(Standard_Real theU1f[], + Standard_Real theU1l[]) const { - if(theCoeffs.mB > 0.0) + if(myCoeffs.mB > 0.0) { - if(theCoeffs.mB + Abs(theCoeffs.mC) < -1.0) + if(myCoeffs.mB + Abs(myCoeffs.mC) < -1.0) {//There is NOT intersection return Standard_False; } - else if(theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0) + else if(myCoeffs.mB + Abs(myCoeffs.mC) <= 1.0) {//U=[0;2*PI]+aFI1 - theU1f[0] = theCoeffs.mFI1; - theU1l[0] = thePeriod + theCoeffs.mFI1; + theU1f[0] = myCoeffs.mFI1; + theU1l[0] = myPeriod + myCoeffs.mFI1; } - else if((1 + theCoeffs.mC <= theCoeffs.mB) && - (theCoeffs.mB <= 1 - theCoeffs.mC)) + else if((1 + myCoeffs.mC <= myCoeffs.mB) && + (myCoeffs.mB <= 1 - myCoeffs.mC)) { - Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB; + Standard_Real anArg = -(myCoeffs.mC + 1) / myCoeffs.mB; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) @@ -2130,15 +2167,15 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, const Standard_Real aDAngle = acos(anArg); //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1) - theU1f[0] = theCoeffs.mFI1; - theU1l[0] = aDAngle + theCoeffs.mFI1; - theU1f[1] = thePeriod - aDAngle + theCoeffs.mFI1; - theU1l[1] = thePeriod + theCoeffs.mFI1; + theU1f[0] = myCoeffs.mFI1; + theU1l[0] = aDAngle + myCoeffs.mFI1; + theU1f[1] = myPeriod - aDAngle + myCoeffs.mFI1; + theU1l[1] = myPeriod + myCoeffs.mFI1; } - else if((1 - theCoeffs.mC <= theCoeffs.mB) && - (theCoeffs.mB <= 1 + theCoeffs.mC)) + else if((1 - myCoeffs.mC <= myCoeffs.mB) && + (myCoeffs.mB <= 1 + myCoeffs.mC)) { - Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB; + Standard_Real anArg = (1 - myCoeffs.mC) / myCoeffs.mB; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) @@ -2147,13 +2184,13 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, const Standard_Real aDAngle = acos(anArg); //U=[aDAngle;2*PI-aDAngle]+aFI1 - theU1f[0] = aDAngle + theCoeffs.mFI1; - theU1l[0] = thePeriod - aDAngle + theCoeffs.mFI1; + theU1f[0] = aDAngle + myCoeffs.mFI1; + theU1l[0] = myPeriod - aDAngle + myCoeffs.mFI1; } - else if(theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0) + else if(myCoeffs.mB - Abs(myCoeffs.mC) >= 1.0) { - Standard_Real anArg1 = (1 - theCoeffs.mC) / theCoeffs.mB, - anArg2 = -(theCoeffs.mC + 1) / theCoeffs.mB; + Standard_Real anArg1 = (1 - myCoeffs.mC) / myCoeffs.mB, + anArg2 = -(myCoeffs.mC + 1) / myCoeffs.mB; if(anArg1 > 1.0) anArg1 = 1.0; if(anArg1 < -1.0) @@ -2168,31 +2205,31 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, //(U=[aDAngle1;aDAngle2]+aFI1) || //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1) - theU1f[0] = aDAngle1 + theCoeffs.mFI1; - theU1l[0] = aDAngle2 + theCoeffs.mFI1; - theU1f[1] = thePeriod - aDAngle2 + theCoeffs.mFI1; - theU1l[1] = thePeriod - aDAngle1 + theCoeffs.mFI1; + theU1f[0] = aDAngle1 + myCoeffs.mFI1; + theU1l[0] = aDAngle2 + myCoeffs.mFI1; + theU1f[1] = myPeriod - aDAngle2 + myCoeffs.mFI1; + theU1l[1] = myPeriod - aDAngle1 + myCoeffs.mFI1; } else { return Standard_False; } } - else if(theCoeffs.mB < 0.0) + else if(myCoeffs.mB < 0.0) { - if(theCoeffs.mB + Abs(theCoeffs.mC) > 1.0) + if(myCoeffs.mB + Abs(myCoeffs.mC) > 1.0) {//There is NOT intersection return Standard_False; } - else if(-theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0) + else if(-myCoeffs.mB + Abs(myCoeffs.mC) <= 1.0) {//U=[0;2*PI]+aFI1 - theU1f[0] = theCoeffs.mFI1; - theU1l[0] = thePeriod + theCoeffs.mFI1; + theU1f[0] = myCoeffs.mFI1; + theU1l[0] = myPeriod + myCoeffs.mFI1; } - else if((-theCoeffs.mC - 1 <= theCoeffs.mB) && - ( theCoeffs.mB <= theCoeffs.mC - 1)) + else if((-myCoeffs.mC - 1 <= myCoeffs.mB) && + ( myCoeffs.mB <= myCoeffs.mC - 1)) { - Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB; + Standard_Real anArg = (1 - myCoeffs.mC) / myCoeffs.mB; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) @@ -2201,15 +2238,15 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, const Standard_Real aDAngle = acos(anArg); //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1) - theU1f[0] = theCoeffs.mFI1; - theU1l[0] = aDAngle + theCoeffs.mFI1; - theU1f[1] = thePeriod - aDAngle + theCoeffs.mFI1; - theU1l[1] = thePeriod + theCoeffs.mFI1; + theU1f[0] = myCoeffs.mFI1; + theU1l[0] = aDAngle + myCoeffs.mFI1; + theU1f[1] = myPeriod - aDAngle + myCoeffs.mFI1; + theU1l[1] = myPeriod + myCoeffs.mFI1; } - else if((theCoeffs.mC - 1 <= theCoeffs.mB) && - (theCoeffs.mB <= -theCoeffs.mB - 1)) + else if((myCoeffs.mC - 1 <= myCoeffs.mB) && + (myCoeffs.mB <= -myCoeffs.mB - 1)) { - Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB; + Standard_Real anArg = -(myCoeffs.mC + 1) / myCoeffs.mB; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) @@ -2218,13 +2255,13 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, const Standard_Real aDAngle = acos(anArg); //U=[aDAngle;2*PI-aDAngle]+aFI1 - theU1f[0] = aDAngle + theCoeffs.mFI1; - theU1l[0] = thePeriod - aDAngle + theCoeffs.mFI1; + theU1f[0] = aDAngle + myCoeffs.mFI1; + theU1l[0] = myPeriod - aDAngle + myCoeffs.mFI1; } - else if(-theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0) + else if(-myCoeffs.mB - Abs(myCoeffs.mC) >= 1.0) { - Standard_Real anArg1 = -(theCoeffs.mC + 1) / theCoeffs.mB, - anArg2 = (1 - theCoeffs.mC) / theCoeffs.mB; + Standard_Real anArg1 = -(myCoeffs.mC + 1) / myCoeffs.mB, + anArg2 = (1 - myCoeffs.mC) / myCoeffs.mB; if(anArg1 > 1.0) anArg1 = 1.0; if(anArg1 < -1.0) @@ -2239,10 +2276,10 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, //(U=[aDAngle1;aDAngle2]+aFI1) || //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1) - theU1f[0] = aDAngle1 + theCoeffs.mFI1; - theU1l[0] = aDAngle2 + theCoeffs.mFI1; - theU1f[1] = thePeriod - aDAngle2 + theCoeffs.mFI1; - theU1l[1] = thePeriod - aDAngle1 + theCoeffs.mFI1; + theU1f[0] = aDAngle1 + myCoeffs.mFI1; + theU1l[0] = aDAngle2 + myCoeffs.mFI1; + theU1f[1] = myPeriod - aDAngle2 + myCoeffs.mFI1; + theU1l[1] = myPeriod - aDAngle1 + myCoeffs.mFI1; } else { @@ -2262,7 +2299,7 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs, //purpose : theNbCritPointsMax contains true number of critical points. // It must be initialized correctly before function calling //======================================================================= -static void CriticalPointsComputing(const stCoeffsValue& theCoeffs, +static void CriticalPointsComputing(const ComputationMethods::stCoeffsValue& theCoeffs, const Standard_Real theUSurf1f, const Standard_Real theUSurf1l, const Standard_Real theUSurf2f, @@ -2272,13 +2309,13 @@ static void CriticalPointsComputing(const stCoeffsValue& theCoeffs, Standard_Integer& theNbCritPointsMax, Standard_Real theU1crit[]) { - //[0...1] - in these points parameter U1 goes through - // the seam-edge of the first cylinder. - //[2...3] - First and last U1 parameter. - //[4...5] - in these points parameter U2 goes through - // the seam-edge of the second cylinder. - //[6...9] - in these points an intersection line goes through - // U-boundaries of the second surface. + //[0...1] - in these points parameter U1 goes through + // the seam-edge of the first cylinder. + //[2...3] - First and last U1 parameter. + //[4...5] - in these points parameter U2 goes through + // the seam-edge of the second cylinder. + //[6...9] - in these points an intersection line goes through + // U-boundaries of the second surface. //[10...11] - Boundary of monotonicity interval of U2(U1) function // (see CylCylMonotonicity() function) @@ -2385,20 +2422,108 @@ static void CriticalPointsComputing(const stCoeffsValue& theCoeffs, } //======================================================================= -//function : IntCyCyTrim +//function : BoundaryEstimation +//purpose : Rough estimation of the parameter range. +//======================================================================= +void WorkWithBoundaries::BoundaryEstimation(const gp_Cylinder& theCy1, + const gp_Cylinder& theCy2, + Bnd_Range& theOutBoxS1, + Bnd_Range& theOutBoxS2) const +{ + const gp_Dir &aD1 = theCy1.Axis().Direction(), + &aD2 = theCy2.Axis().Direction(); + const Standard_Real aR1 = theCy1.Radius(), + aR2 = theCy2.Radius(); + + //Let consider a parallelogram. Its edges are parallel to aD1 and aD2. + //Its altitudes are equal to 2*aR1 and 2*aR2 (diameters of the cylinders). + //In fact, this parallelogram is a projection of the cylinders to the plane + //created by the intersected axes aD1 and aD2 (if the axes are skewed then + //one axis can be translated by parallel shifting till intersection). + + const Standard_Real aCosA = aD1.Dot(aD2); + const Standard_Real aSqSinA = 1.0 - aCosA * aCosA; + + //If sine is small then it can be compared with angle. + if (aSqSinA < Precision::Angular()*Precision::Angular()) + return; + + //Half of delta V. Delta V is a distance between + //projections of two opposite parallelogram vertices + //(joined by the maximal diagonal) to the cylinder axis. + const Standard_Real aSinA = sqrt(aSqSinA); + const Standard_Real aHDV1 = (aR1 * aCosA + aR2)/aSinA, + aHDV2 = (aR2 * aCosA + aR1)/aSinA; + +#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG + //The code in this block is created for test only.It is stupidly to create + //OCCT-test for the method, which will be changed possibly never. + std::cout << "Reference values: aHDV1 = " << aHDV1 << "; aHDV2 = " << aHDV2 << std::endl; +#endif + + //V-parameters of intersection point of the axes (in case of skewed axes, + //see comment above). + Standard_Real aV01 = 0.0, aV02 = 0.0; + ExtremaLineLine(theCy1.Axis(), theCy2.Axis(), aCosA, aSqSinA, aV01, aV02); + + theOutBoxS1.Add(aV01 - aHDV1); + theOutBoxS1.Add(aV01 + aHDV1); + + theOutBoxS2.Add(aV02 - aHDV2); + theOutBoxS2.Add(aV02 + aHDV2); + + theOutBoxS1.Enlarge(Precision::Confusion()); + theOutBoxS2.Enlarge(Precision::Confusion()); + + Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0; + myUVSurf1.Get(aU1, aV1, aU2, aV2); + theOutBoxS1.Common(Bnd_Range(aV1, aV2)); + + myUVSurf2.Get(aU1, aV1, aU2, aV2); + theOutBoxS2.Common(Bnd_Range(aV1, aV2)); +} + +//======================================================================= +//function : IntCyCy //purpose : //======================================================================= -Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, - const IntSurf_Quadric& theQuad2, - const Standard_Real theTol3D, - const Standard_Real theTol2D, - const Bnd_Box2d& theUVSurf1, - const Bnd_Box2d& theUVSurf2, - const Standard_Boolean isTheReverse, - Standard_Boolean& isTheEmpty, - IntPatch_SequenceOfLine& theSlin, - IntPatch_SequenceOfPoint& theSPnt) +Standard_Boolean IntCyCy( const IntSurf_Quadric& theQuad1, + const IntSurf_Quadric& theQuad2, + const Standard_Real theTol3D, + const Standard_Real theTol2D, + const Bnd_Box2d& theUVSurf1, + const Bnd_Box2d& theUVSurf2, + const Standard_Boolean isTheReverse, + Standard_Boolean& isTheEmpty, + Standard_Boolean& isTheSameSurface, + Standard_Boolean& isTheMultiplePoint, + IntPatch_SequenceOfLine& theSlin, + IntPatch_SequenceOfPoint& theSPnt) { + isTheEmpty = Standard_True; + isTheSameSurface = Standard_False; + isTheMultiplePoint = Standard_False; + theSlin.Clear(); + theSPnt.Clear(); + + const gp_Cylinder &aCyl1 = theQuad1.Cylinder(), + &aCyl2 = theQuad2.Cylinder(); + + IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D); + + if (!anInter.IsDone()) + { + return Standard_False; + } + + if(anInter.TypeInter() != IntAna_NoGeometricSolution) + { + return CyCyAnalyticalIntersect( theQuad1, theQuad2, anInter, + theTol3D, isTheReverse, isTheEmpty, + isTheSameSurface, isTheMultiplePoint, + theSlin, theSPnt); + } + Standard_Real aUSurf1f = 0.0, //const aUSurf1l = 0.0, aVSurf1f = 0.0, @@ -2411,27 +2536,6 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l); theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l); - const gp_Cylinder& aCyl1 = theQuad1.Cylinder(), - aCyl2 = theQuad2.Cylinder(); - - IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D); - - if (!anInter.IsDone()) - { - return Standard_False; - } - - IntAna_ResultType aTypInt = anInter.TypeInter(); - - if(aTypInt != IntAna_NoGeometricSolution) - { //It is not necessary (because result is an analytic curve) or - //it is impossible to make Walking-line. - - return Standard_False; - } - - theSlin.Clear(); - theSPnt.Clear(); const Standard_Integer aNbMaxPoints = 2000; const Standard_Integer aNbMinPoints = 200; const Standard_Integer aNbPoints = Min(Max(aNbMinPoints, @@ -2444,14 +2548,23 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, const Standard_Integer aNbWLines = 2; - const stCoeffsValue anEquationCoeffs(aCyl1, aCyl2); + const ComputationMethods::stCoeffsValue anEquationCoeffs(aCyl1, aCyl2); + + const WorkWithBoundaries aBoundWork(theQuad1, theQuad2, anEquationCoeffs, + theUVSurf1, theUVSurf2, aNbWLines, + aPeriod, theTol3D, theTol2D, isTheReverse); + + Bnd_Range aRangeS1, aRangeS2; + aBoundWork.BoundaryEstimation(aCyl1, aCyl2, aRangeS1, aRangeS2); + if (aRangeS1.IsVoid() || aRangeS2.IsVoid()) + return Standard_True; //Boundaries const Standard_Integer aNbOfBoundaries = 2; Standard_Real aU1f[aNbOfBoundaries] = {-Precision::Infinite(), -Precision::Infinite()}; Standard_Real aU1l[aNbOfBoundaries] = {Precision::Infinite(), Precision::Infinite()}; - if(!BoundariesComputing(anEquationCoeffs, aPeriod, aU1f, aU1l)) + if(!aBoundWork.BoundariesComputing(aU1f, aU1l)) return Standard_True; for(Standard_Integer i = 0; i < aNbOfBoundaries; i++) @@ -2620,7 +2733,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, //or on seam-edge strictly (if it is possible). Standard_Real aTol = theTol2D; - CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i], &aTol); + ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs, + aU2[i], &aTol); InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False); aTol = Max(aTol, theTol2D); @@ -2636,7 +2750,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, } else { - CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]); + ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]); InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False); } @@ -2651,7 +2765,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, //correct value. Standard_Boolean isIncreasing = Standard_True; - CylCylMonotonicity(anU1+aStepMin, i, anEquationCoeffs, aPeriod, isIncreasing); + ComputationMethods::CylCylMonotonicity(anU1+aStepMin, i, anEquationCoeffs, + aPeriod, isIncreasing); //If U2(U1) is increasing and U2 is considered to be equal aUSurf2l //then after the next step (when U1 will be increased) U2 will be @@ -2690,7 +2805,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, } } - CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs, aV1[i], aV2[i]); + ComputationMethods::CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs, + aV1[i], aV2[i]); if(isFirst) { @@ -2771,11 +2887,9 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, } } - AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs, - theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod, - anU1, aU2[i], aV1[i], aV1Prev[i], - aV2[i], aV2Prev[i], i, isTheReverse, - isForce, isFound1, isFound2); + aBoundWork.AddBoundaryPoint(aWLine[i], anU1, aU2[i], aV1[i], aV1Prev[i], + aV2[i], aV2Prev[i], i, isForce, + isFound1, isFound2); const Standard_Boolean isPrevVBound = !isVIntersect && ((Abs(aV1Prev[i] - aVSurf1f) <= theTol2D) || @@ -2783,7 +2897,6 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, (Abs(aV2Prev[i] - aVSurf2f) <= theTol2D) || (Abs(aV2Prev[i] - aVSurf2l) <= theTol2D)); - aV1Prev[i] = aV1[i]; aV2Prev[i] = aV2[i]; @@ -2798,7 +2911,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, aWLFindStatus[i] = WLFStatus_Exist; } - if(( aWLFindStatus[i] != WLFStatus_Broken) || (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl)) + if( (aWLFindStatus[i] != WLFStatus_Broken) || + (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl)) { if(aWLine[i]->NbPnts() > 0) { @@ -2872,11 +2986,9 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False; - AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs, - theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod, - anU1, aU2[i], aV1[i], aV1Prev[i], - aV2[i], aV2Prev[i], i, isTheReverse, - Standard_False, isFound1, isFound2); + aBoundWork.AddBoundaryPoint(aWLine[i], anU1, aU2[i], aV1[i], aV1Prev[i], + aV2[i], aV2Prev[i], i, + Standard_False, isFound1, isFound2); if(isFound1 || isFound2) { @@ -2952,7 +3064,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, continue; } - CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs, aU2[i], aV1[i], aV2[i]); + ComputationMethods::CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs, + aU2[i], aV1[i], aV2[i]); AddPointIntoWL( theQuad1, theQuad2, anEquationCoeffs, isTheReverse, Standard_True, gp_Pnt2d(anUmaxAdded, aV1[i]), @@ -2969,8 +3082,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, //Step computing { - const Standard_Real aDeltaV1 = (aVSurf1l - aVSurf1f)/IntToReal(aNbPoints), - aDeltaV2 = (aVSurf2l - aVSurf2f)/IntToReal(aNbPoints); + const Standard_Real aDeltaV1 = aRangeS1.Delta()/IntToReal(aNbPoints), + aDeltaV2 = aRangeS2.Delta()/IntToReal(aNbPoints); math_Matrix aMatr(1, 3, 1, 5); @@ -3078,8 +3191,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, isAddedIntoWL[i] = Standard_True; SeekAdditionalPoints( theQuad1, theQuad2, aWLine[i]->Curve(), anEquationCoeffs, i, aNbPoints, 1, - aWLine[i]->NbPnts(), aUSurf2f, aUSurf2l, - theTol2D, aPeriod, isTheReverse); + aWLine[i]->NbPnts(), theTol2D, aPeriod, + isTheReverse); aWLine[i]->ComputeVertexParameters(theTol3D); theSlin.Append(aWLine[i]); @@ -3090,7 +3203,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, isAddedIntoWL[i] = Standard_False; } -#ifdef OCCT_DEBUG +#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG //aWLine[i]->Dump(); #endif } @@ -3104,7 +3217,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, { for(Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); aNbLin++) { - Handle(IntPatch_WLine) aWLine1 (Handle(IntPatch_WLine)::DownCast(theSlin.Value(aNbLin))); + Handle(IntPatch_WLine) aWLine1 (Handle(IntPatch_WLine):: + DownCast(theSlin.Value(aNbLin))); const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1); const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aWLine1->NbPnts()); @@ -3152,7 +3266,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, for (Standard_Integer i = 0; i < aNbWLines; i++) { Standard_Real anU2t = 0.0; - if(!CylCylComputeParameters(anUl, i, anEquationCoeffs, anU2t)) + if(!ComputationMethods::CylCylComputeParameters(anUl, i, anEquationCoeffs, anU2t)) continue; const Standard_Real aDU2 = Abs(anU2t - aCurU2); @@ -3168,7 +3282,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, { Standard_Real anU2 = 0.0, anV1 = 0.0, anV2 = 0.0; Standard_Boolean isDone = - CylCylComputeParameters(anUf, anIndex, anEquationCoeffs, + ComputationMethods::CylCylComputeParameters(anUf, anIndex, anEquationCoeffs, anU2, anV1, anV2); anUf += aDU; @@ -3177,8 +3291,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, continue; } - if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse, Standard_True, - gp_Pnt2d(anUf, anV1), gp_Pnt2d(anU2, anV2), + if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse, + Standard_True, gp_Pnt2d(anUf, anV1), gp_Pnt2d(anU2, anV2), aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod, aWLine->Curve(), anIndex, theTol3D, @@ -3192,8 +3306,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, { SeekAdditionalPoints( theQuad1, theQuad2, aWLine->Curve(), anEquationCoeffs, anIndex, aNbMinPoints, - 1, aWLine->NbPnts(), aUSurf2f, aUSurf2l, - theTol2D, aPeriod, isTheReverse); + 1, aWLine->NbPnts(), theTol2D, aPeriod, + isTheReverse); aWLine->ComputeVertexParameters(theTol3D); theSlin.Append(aWLine); @@ -3614,7 +3728,7 @@ Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1, ptvalid = curvsol.Value(para); alig = new IntPatch_ALine(curvsol,Standard_False); kept = Standard_True; - //-- cout << "Transition indeterminee" << endl; + //-- std::cout << "Transition indeterminee" << std::endl; } if (kept) { Standard_Boolean Nfirstp = !firstp; @@ -3696,46 +3810,3 @@ Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy, // return bFind; } -//======================================================================= -//function : IsToReverse -//purpose : -//======================================================================= -Standard_Boolean IsToReverse(const gp_Cylinder& Cy1, - const gp_Cylinder& Cy2, - const Standard_Real Tol) -{ - Standard_Boolean bRet; - Standard_Real aR1,aR2, dR, aSc1, aSc2; - // - bRet=Standard_False; - // - aR1=Cy1.Radius(); - aR2=Cy2.Radius(); - dR=aR1-aR2; - if (dR<0.) { - dR=-dR; - } - if (dR>Tol) { - bRet=aR1>aR2; - } - else { - gp_Dir aDZ(0.,0.,1.); - // - const gp_Dir& aDir1=Cy1.Axis().Direction(); - aSc1=aDir1*aDZ; - if (aSc1<0.) { - aSc1=-aSc1; - } - // - const gp_Dir& aDir2=Cy2.Axis().Direction(); - aSc2=aDir2*aDZ; - if (aSc2<0.) { - aSc2=-aSc2; - } - // - if (aSc2DomainIsInfinite() || theD2->DomainIsInfinite()) - { - GeomGeomPerfom( theS1, theD1, theS2, theD2, TolArc, - TolTang, ListOfPnts, RestrictLine, - typs1, typs2, theIsReqToKeepRLine); - } - else - { - GeomGeomPerfomTrimSurf( theS1, theD1, theS2, theD2, - TolArc, TolTang, ListOfPnts, RestrictLine, - typs1, typs2, theIsReqToKeepRLine); - } + GeomGeomPerfom( theS1, theD1, theS2, theD2, TolArc, + TolTang, ListOfPnts, RestrictLine, + typs1, typs2, theIsReqToKeepRLine); } else { @@ -1182,16 +1173,8 @@ void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)& theS1, } else if(ts1 == 1) { - if(theD1->DomainIsInfinite() || theD2->DomainIsInfinite()) - { - GeomGeomPerfom(theS1, theD1, theS2, theD2, TolArc, - TolTang, ListOfPnts, RestrictLine, typs1, typs2, theIsReqToKeepRLine); - } - else - { - GeomGeomPerfomTrimSurf(theS1, theD1, theS2, theD2, - TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2, theIsReqToKeepRLine); - } + GeomGeomPerfom(theS1, theD1, theS2, theD2, TolArc, + TolTang, ListOfPnts, RestrictLine, typs1, typs2, theIsReqToKeepRLine); } if(!theIsReqToPostWLProc) @@ -1438,7 +1421,36 @@ void IntPatch_Intersection::GeomGeomPerfom(const Handle(Adaptor3d_HSurface)& the slin.Append(wlin); } else + { + if(line->ArcType() == IntPatch_Walking) + { + Handle(IntPatch_WLine)::DownCast(line)->EnablePurging(Standard_False); + } + slin.Append(line); + } + } + + for (Standard_Integer i = 1; i <= interii.NbPnts(); i++) + { + spnt.Append(interii.Point(i)); + } + + if((theTyps1 == GeomAbs_Cylinder) && (theTyps2 == GeomAbs_Cylinder)) + { + IntPatch_WLineTool::JoinWLines( slin, spnt, TolTang, + theS1->IsUPeriodic()? theS1->UPeriod() : 0.0, + theS2->IsUPeriodic()? theS2->UPeriod() : 0.0, + theS1->IsVPeriodic()? theS1->VPeriod() : 0.0, + theS2->IsVPeriodic()? theS2->VPeriod() : 0.0, + theS1->FirstUParameter(), + theS1->LastUParameter(), + theS1->FirstVParameter(), + theS1->LastVParameter(), + theS2->FirstUParameter(), + theS2->LastUParameter(), + theS2->FirstVParameter(), + theS2->LastVParameter()); } if(isQuadSet) @@ -1472,11 +1484,6 @@ void IntPatch_Intersection::GeomGeomPerfom(const Handle(Adaptor3d_HSurface)& the theS2->IsVPeriodic()? theS2->VPeriod() : 0.0, aBx1, aBx2); } - - for (Standard_Integer i = 1; i <= interii.NbPnts(); i++) - { - spnt.Append(interii.Point(i)); - } } } else @@ -1568,84 +1575,6 @@ void IntPatch_Intersection:: } } -//======================================================================= -//function : GeomGeomPerfomTrimSurf -//purpose : This function returns ready walking-line (which is not need -// in convertation) as an intersection line between two -// trimmed surfaces. -//======================================================================= -void IntPatch_Intersection:: - GeomGeomPerfomTrimSurf( const Handle(Adaptor3d_HSurface)& theS1, - const Handle(Adaptor3d_TopolTool)& theD1, - const Handle(Adaptor3d_HSurface)& theS2, - const Handle(Adaptor3d_TopolTool)& theD2, - const Standard_Real theTolArc, - const Standard_Real theTolTang, - IntSurf_ListOfPntOn2S& theListOfPnts, - const Standard_Boolean RestrictLine, - const GeomAbs_SurfaceType theTyps1, - const GeomAbs_SurfaceType theTyps2, - const Standard_Boolean theIsReqToKeepRLine) -{ - IntSurf_Quadric Quad1,Quad2; - - if((theTyps1 == GeomAbs_Cylinder) && (theTyps2 == GeomAbs_Cylinder)) - { - IntPatch_ImpImpIntersection anInt; - anInt.Perform(theS1, theD1, theS2, theD2, myTolArc, - myTolTang, Standard_True, theIsReqToKeepRLine); - - done = anInt.IsDone(); - - if(done) - { - empt = anInt.IsEmpty(); - if (!empt) - { - tgte = anInt.TangentFaces(); - if (tgte) - oppo = anInt.OppositeFaces(); - - const Standard_Integer aNbLin = anInt.NbLines(); - const Standard_Integer aNbPts = anInt.NbPnts(); - - for(Standard_Integer aLID = 1; aLID <= aNbLin; aLID++) - { - const Handle(IntPatch_Line)& aLine = anInt.Line(aLID); - slin.Append(aLine); - } - - for(Standard_Integer aPID = 1; aPID <= aNbPts; aPID++) - { - const IntPatch_Point& aPoint = anInt.Point(aPID); - spnt.Append(aPoint); - } - - IntPatch_WLineTool::JoinWLines( slin, spnt, theTolTang, - theS1->IsUPeriodic()? theS1->UPeriod() : 0.0, - theS2->IsUPeriodic()? theS2->UPeriod() : 0.0, - theS1->IsVPeriodic()? theS1->VPeriod() : 0.0, - theS2->IsVPeriodic()? theS2->VPeriod() : 0.0, - theS1->FirstUParameter(), - theS1->LastUParameter(), - theS1->FirstVParameter(), - theS1->LastVParameter(), - theS2->FirstUParameter(), - theS2->LastUParameter(), - theS2->FirstVParameter(), - theS2->LastVParameter()); - } - } - } - else - { - GeomGeomPerfom(theS1, theD1, theS2, theD2, - theTolArc, theTolTang, theListOfPnts, - RestrictLine, theTyps1, theTyps2, theIsReqToKeepRLine); - } -} - - void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_TopolTool)& D1, const Handle(Adaptor3d_HSurface)& S2, diff --git a/src/IntPatch/IntPatch_Intersection.hxx b/src/IntPatch/IntPatch_Intersection.hxx index cfe6b037bb..96a4dac4b4 100644 --- a/src/IntPatch/IntPatch_Intersection.hxx +++ b/src/IntPatch/IntPatch_Intersection.hxx @@ -167,14 +167,6 @@ private: //! will always keep both lines even if they are coincided. Standard_EXPORT void GeomGeomPerfom (const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_TopolTool)& D1, const Handle(Adaptor3d_HSurface)& S2, const Handle(Adaptor3d_TopolTool)& D2, const Standard_Real TolArc, const Standard_Real TolTang, IntSurf_ListOfPntOn2S& LOfPnts, const Standard_Boolean RestrictLine, const GeomAbs_SurfaceType typs1, const GeomAbs_SurfaceType typs2, const Standard_Boolean theIsReqToKeepRLine = Standard_False); - //! Flag theIsReqToKeepRLine has been enterred only for - //! compatibility with TopOpeBRep package. It shall be deleted - //! after deleting TopOpeBRep. - //! When intersection result returns IntPatch_RLine and another - //! IntPatch_Line (not restriction) we (in case of theIsReqToKeepRLine==TRUE) - //! will always keep both lines even if they are coincided. - Standard_EXPORT void GeomGeomPerfomTrimSurf (const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_TopolTool)& D1, const Handle(Adaptor3d_HSurface)& S2, const Handle(Adaptor3d_TopolTool)& D2, const Standard_Real TolArc, const Standard_Real TolTang, IntSurf_ListOfPntOn2S& LOfPnts, const Standard_Boolean RestrictLine, const GeomAbs_SurfaceType typs1, const GeomAbs_SurfaceType typs2, const Standard_Boolean theIsReqToKeepRLine = Standard_False); - Standard_EXPORT void GeomParamPerfom (const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_TopolTool)& D1, const Handle(Adaptor3d_HSurface)& S2, const Handle(Adaptor3d_TopolTool)& D2, const Standard_Boolean isNotAnalitical, const GeomAbs_SurfaceType typs1, const GeomAbs_SurfaceType typs2); diff --git a/src/IntPatch/IntPatch_WLineTool.cxx b/src/IntPatch/IntPatch_WLineTool.cxx index ca38c1c2b1..bb09d22da2 100644 --- a/src/IntPatch/IntPatch_WLineTool.cxx +++ b/src/IntPatch/IntPatch_WLineTool.cxx @@ -17,6 +17,9 @@ #include #include +// It is pure empirical value. +const Standard_Real IntPatch_WLineTool::myMaxConcatAngle = M_PI/6; + //Bit-mask is used for information about //the operation made in //IntPatch_WLineTool::ExtendTwoWlinesToEachOther() method. @@ -830,27 +833,26 @@ static Standard_Boolean IsOutOfDomain(const Bnd_Box2d& theBoxS1, } //======================================================================= -//function : CheckArguments +//function : CheckArgumentsToExtend //purpose : Check if extending is possible // (see IntPatch_WLineTool::ExtendTwoWlinesToEachOther) //======================================================================= -static Standard_Boolean CheckArguments(const IntSurf_Quadric& theS1, - const IntSurf_Quadric& theS2, - const IntSurf_PntOn2S& thePtWL1, - const IntSurf_PntOn2S& thePtWL2, - IntSurf_PntOn2S& theNewPoint, - const gp_Vec& theVec1, - const gp_Vec& theVec2, - const gp_Vec& theVec3, - const Bnd_Box2d& theBoxS1, - const Bnd_Box2d& theBoxS2, - const Standard_Real theToler3D, - const Standard_Real theU1Period, - const Standard_Real theU2Period, - const Standard_Real theV1Period, - const Standard_Real theV2Period) +Standard_Boolean CheckArgumentsToExtend(const IntSurf_Quadric& theS1, + const IntSurf_Quadric& theS2, + const IntSurf_PntOn2S& thePtWL1, + const IntSurf_PntOn2S& thePtWL2, + IntSurf_PntOn2S& theNewPoint, + const gp_Vec& theVec1, + const gp_Vec& theVec2, + const gp_Vec& theVec3, + const Bnd_Box2d& theBoxS1, + const Bnd_Box2d& theBoxS2, + const Standard_Real theToler3D, + const Standard_Real theU1Period, + const Standard_Real theU2Period, + const Standard_Real theV1Period, + const Standard_Real theV2Period) { - const Standard_Real aMaxAngle = M_PI/6; //30 degree const Standard_Real aSqToler = theToler3D*theToler3D; if(theVec3.SquareMagnitude() <= aSqToler) @@ -858,9 +860,9 @@ static Standard_Boolean CheckArguments(const IntSurf_Quadric& theS1, return Standard_False; } - if((theVec1.Angle(theVec2) > aMaxAngle) || - (theVec1.Angle(theVec3) > aMaxAngle) || - (theVec2.Angle(theVec3) > aMaxAngle)) + if((theVec1.Angle(theVec2) > IntPatch_WLineTool::myMaxConcatAngle) || + (theVec1.Angle(theVec3) > IntPatch_WLineTool::myMaxConcatAngle) || + (theVec2.Angle(theVec3) > IntPatch_WLineTool::myMaxConcatAngle)) { return Standard_False; } @@ -907,6 +909,20 @@ static Standard_Boolean CheckArguments(const IntSurf_Quadric& theS1, return Standard_True; } +//======================================================================= +//function : CheckArgumentsToJoin +//purpose : Check if joining is possible +// (see IntPatch_WLineTool::JoinWLines) +//======================================================================= +Standard_Boolean CheckArgumentsToJoin(const gp_Vec& theVec1, + const gp_Vec& theVec2) +{ + // [0, PI] - range + const Standard_Real anAngle = theVec1.Angle(theVec2); + + return (anAngle < IntPatch_WLineTool::myMaxConcatAngle); +} + //======================================================================= //function : ExtendTwoWLFirstFirst //purpose : Performs extending theWLine1 and theWLine2 through their @@ -932,10 +948,10 @@ static void ExtendTwoWLFirstFirst(const IntSurf_Quadric& theS1, Standard_Boolean &theHasBeenJoined) { IntSurf_PntOn2S aPOn2S; - if(!CheckArguments(theS1, theS2, thePtWL1, thePtWL2, aPOn2S, - theVec1, theVec2, theVec3, - theBoxS1, theBoxS2, theToler3D, - theU1Period, theU2Period, theV1Period, theV2Period)) + if(!CheckArgumentsToExtend(theS1, theS2, thePtWL1, thePtWL2, aPOn2S, + theVec1, theVec2, theVec3, + theBoxS1, theBoxS2, theToler3D, + theU1Period, theU2Period, theV1Period, theV2Period)) { return; } @@ -1006,10 +1022,10 @@ static void ExtendTwoWLFirstLast(const IntSurf_Quadric& theS1, Standard_Boolean &theHasBeenJoined) { IntSurf_PntOn2S aPOn2S; - if(!CheckArguments(theS1, theS2, thePtWL1, thePtWL2, aPOn2S, - theVec1, theVec2, theVec3, - theBoxS1, theBoxS2, theToler3D, - theU1Period, theU2Period, theV1Period, theV2Period)) + if(!CheckArgumentsToExtend(theS1, theS2, thePtWL1, thePtWL2, aPOn2S, + theVec1, theVec2, theVec3, + theBoxS1, theBoxS2, theToler3D, + theU1Period, theU2Period, theV1Period, theV2Period)) { return; } @@ -1078,10 +1094,10 @@ static void ExtendTwoWLLastFirst(const IntSurf_Quadric& theS1, Standard_Boolean &theHasBeenJoined) { IntSurf_PntOn2S aPOn2S; - if(!CheckArguments(theS1, theS2, thePtWL1, thePtWL2, aPOn2S, - theVec1, theVec2, theVec3, - theBoxS1, theBoxS2, theToler3D, - theU1Period, theU2Period, theV1Period, theV2Period)) + if(!CheckArgumentsToExtend(theS1, theS2, thePtWL1, thePtWL2, aPOn2S, + theVec1, theVec2, theVec3, + theBoxS1, theBoxS2, theToler3D, + theU1Period, theU2Period, theV1Period, theV2Period)) { return; } @@ -1146,10 +1162,10 @@ static void ExtendTwoWLLastLast(const IntSurf_Quadric& theS1, Standard_Boolean &theHasBeenJoined) { IntSurf_PntOn2S aPOn2S; - if(!CheckArguments(theS1, theS2, thePtWL1, thePtWL2, aPOn2S, - theVec1, theVec2, theVec3, - theBoxS1, theBoxS2, theToler3D, - theU1Period, theU2Period, theV1Period, theV2Period)) + if(!CheckArgumentsToExtend(theS1, theS2, thePtWL1, thePtWL2, aPOn2S, + theVec1, theVec2, theVec3, + theBoxS1, theBoxS2, theToler3D, + theU1Period, theU2Period, theV1Period, theV2Period)) { return; } @@ -1387,10 +1403,19 @@ void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin, { const IntSurf_PntOn2S& aPt1 = aWLine1->Point(2); const IntSurf_PntOn2S& aPt2 = aWLine2->Point(2); - if(!IsSeamOrBound(aPt1, aPt2, aPntFWL1, theU1Period, theU2Period, - theV1Period, theV2Period, theUfSurf1, theUlSurf1, - theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2, - theVfSurf2, theVlSurf2)) + Standard_Boolean aCond = + CheckArgumentsToJoin(gp_Vec(aPntFWL1.Value(), aPt1.Value()), + gp_Vec(aPt2.Value(), aPntFWL2.Value())); + + aCond = aCond && !IsSeamOrBound(aPt1, aPt2, aPntFWL1, + theU1Period, theU2Period, + theV1Period, theV2Period, + theUfSurf1, theUlSurf1, + theVfSurf1, theVlSurf1, + theUfSurf2, theUlSurf2, + theVfSurf2, theVlSurf2); + + if(aCond) { aWLine1->ClearVertexes(); for(Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++) @@ -1413,10 +1438,20 @@ void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin, { const IntSurf_PntOn2S& aPt1 = aWLine1->Point(2); const IntSurf_PntOn2S& aPt2 = aWLine2->Point(aNbPntsWL2-1); - if(!IsSeamOrBound(aPt1, aPt2, aPntFWL1, theU1Period, theU2Period, - theV1Period, theV2Period, theUfSurf1, theUlSurf1, - theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2, - theVfSurf2, theVlSurf2)) + + Standard_Boolean aCond = + CheckArgumentsToJoin(gp_Vec(aPntFWL1.Value(), aPt1.Value()), + gp_Vec(aPt2.Value(), aPntLWL2.Value())); + + aCond = aCond && !IsSeamOrBound(aPt1, aPt2, aPntFWL1, + theU1Period, theU2Period, + theV1Period, theV2Period, + theUfSurf1, theUlSurf1, + theVfSurf1, theVlSurf1, + theUfSurf2, theUlSurf2, + theVfSurf2, theVlSurf2); + + if(aCond) { aWLine1->ClearVertexes(); for(Standard_Integer aNPt = aNbPntsWL2; aNPt >= 1; aNPt--) @@ -1439,10 +1474,20 @@ void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin, { const IntSurf_PntOn2S& aPt1 = aWLine1->Point(aNbPntsWL1-1); const IntSurf_PntOn2S& aPt2 = aWLine2->Point(2); - if(!IsSeamOrBound(aPt1, aPt2, aPntLWL1, theU1Period, theU2Period, - theV1Period, theV2Period, theUfSurf1, theUlSurf1, - theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2, - theVfSurf2, theVlSurf2)) + + Standard_Boolean aCond = + CheckArgumentsToJoin(gp_Vec(aPt1.Value(), aPntLWL1.Value()), + gp_Vec(aPntFWL2.Value(), aPt2.Value())); + + aCond = aCond && !IsSeamOrBound(aPt1, aPt2, aPntLWL1, + theU1Period, theU2Period, + theV1Period, theV2Period, + theUfSurf1, theUlSurf1, + theVfSurf1, theVlSurf1, + theUfSurf2, theUlSurf2, + theVfSurf2, theVlSurf2); + + if(aCond) { aWLine1->ClearVertexes(); for(Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++) @@ -1465,10 +1510,20 @@ void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin, { const IntSurf_PntOn2S& aPt1 = aWLine1->Point(aNbPntsWL1-1); const IntSurf_PntOn2S& aPt2 = aWLine2->Point(aNbPntsWL2-1); - if(!IsSeamOrBound(aPt1, aPt2, aPntLWL1, theU1Period, theU2Period, - theV1Period, theV2Period, theUfSurf1, theUlSurf1, - theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2, - theVfSurf2, theVlSurf2)) + + Standard_Boolean aCond = + CheckArgumentsToJoin(gp_Vec(aPt1.Value(), aPntLWL1.Value()), + gp_Vec(aPntLWL2.Value(), aPt2.Value())); + + aCond = aCond && !IsSeamOrBound(aPt1, aPt2, aPntLWL1, + theU1Period, theU2Period, + theV1Period, theV2Period, + theUfSurf1, theUlSurf1, + theVfSurf1, theVlSurf1, + theUfSurf2, theUlSurf2, + theVfSurf2, theVlSurf2); + + if(aCond) { aWLine1->ClearVertexes(); for(Standard_Integer aNPt = aNbPntsWL2; aNPt >= 1; aNPt--) diff --git a/src/IntPatch/IntPatch_WLineTool.hxx b/src/IntPatch/IntPatch_WLineTool.hxx index 90b0d5d71f..02e4584a10 100644 --- a/src/IntPatch/IntPatch_WLineTool.hxx +++ b/src/IntPatch/IntPatch_WLineTool.hxx @@ -96,6 +96,9 @@ public: const Standard_Real theV2Period, const Bnd_Box2d& theBoxS1, const Bnd_Box2d& theBoxS2); + +//! Max angle to concatenate two WLines to avoid result with C0-continuity + static const Standard_Real myMaxConcatAngle; }; #endif \ No newline at end of file diff --git a/tests/blend/simple/Q6 b/tests/blend/simple/Q6 index 1f3e4adb13..1c2a537fa4 100644 --- a/tests/blend/simple/Q6 +++ b/tests/blend/simple/Q6 @@ -3,4 +3,4 @@ tscale s 0 0 0 SCALE explode s E blend result s SCALE*2 s_5 -checkprops result -s 1.65391e+08 -eps 0.1 +checkprops result -s 1.65391e+008 diff --git a/tests/bugs/modalg_5/bug25742_2 b/tests/bugs/modalg_5/bug25742_2 index 7e4b81ef8b..6dad90a0a2 100755 --- a/tests/bugs/modalg_5/bug25742_2 +++ b/tests/bugs/modalg_5/bug25742_2 @@ -31,16 +31,19 @@ explode b1 f explode b2 f smallview -donly b1_4 b2_1 +donly b1_4 fit +display b2_1 - dchrono h reset dchrono h start bopcurves b1_4 b2_1 -2d dchrono h stop + +checkview -screenshot -2d -path ${imagedir}/${test_image}_1.png + set q [dchrono h show] regexp {CPU user time: ([-0-9.+eE]+) seconds} $q full z @@ -56,10 +59,10 @@ if { $z > ${max_time} } { mksurface s1 b1_4 mksurface s2 b2_1 -dchrono h2 stop -set q2 [dchrono h2 show] +dchrono h2 reset +dchrono h2 start -set CurveNumb [intersect i s1 s2] +set CurveNumb [intersect resi s1 s2] dchrono h2 stop set q2 [dchrono h2 show] @@ -79,4 +82,8 @@ if { [llength ${CurveNumb}] < 1 } { puts "OK : Good intersection" } -checkview -screenshot -2d -path ${imagedir}/${test_image}.png +don resi* +fit +display s1 s2 + +checkview -screenshot -2d -path ${imagedir}/${test_image}_2.png diff --git a/tests/bugs/modalg_6/bug23178 b/tests/bugs/modalg_6/bug23178 new file mode 100644 index 0000000000..a16ee989dc --- /dev/null +++ b/tests/bugs/modalg_6/bug23178 @@ -0,0 +1,43 @@ +puts "================" +puts "OCC23178" +puts "================" +puts "" +####################################################################### +# Intersection of cylinders fails to produce results +####################################################################### + +set GoodNbCurv 6 + +foreach c [directory result*] { + unset $c +} + +cylinder s1 4.999999999813, 1.35836143386791, 0.20975890778078 -1, 0, 0 0, 0, 1 10 +cylinder s2 9.999999999813, 9.1401960568287, 0.11837300583107 0, -0.75870713028692, 0.651431877061437 0, -0.651431877061437, -0.75870713028692 5 + +intersect result s1 s2 + +foreach c [directory result*] { + bounds $c U1 U2 + + if {[dval U2-U1] < 1.0e-9} { + puts "Error: Wrong curve's range!" + } + + xdistcs $c s1 U1 U2 10 2.0e-7 + xdistcs $c s2 U1 U2 10 2.0e-7 +} + +set NbCurv [llength [directory result*]] + +if { $NbCurv == $GoodNbCurv } { + puts "OK: Number of curves is good!" +} else { + puts "Error: $GoodNbCurv is expected but $NbCurv is found!" +} + +smallview +don result* +fit +disp s1 s2 +checkview -screenshot -2d -path ${imagedir}/${test_image}.png diff --git a/tests/bugs/modalg_6/bug27310_2 b/tests/bugs/modalg_6/bug27310_2 index fe5faee133..8c578725be 100644 --- a/tests/bugs/modalg_6/bug27310_2 +++ b/tests/bugs/modalg_6/bug27310_2 @@ -6,7 +6,7 @@ puts "" # Huge tolerance obtained in the result of intersection of two cylindrical faces ################################################# -set ExpTol 0.00027580830315682321 +set ExpTol 6.9230965382090094e-006 set GoodNbCurv 2 restore [locate_data_file OCC496a.brep] a diff --git a/tests/bugs/modalg_6/bug27766 b/tests/bugs/modalg_6/bug27766 index d992b88ed2..c039cd8efc 100644 --- a/tests/bugs/modalg_6/bug27766 +++ b/tests/bugs/modalg_6/bug27766 @@ -7,7 +7,7 @@ puts "" ################################################# set ExpTol 1.0e-7 -set GoodNbCurv 5 +set GoodNbCurv 3 restore [locate_data_file bug27761_c1.brep] c1 restore [locate_data_file bug27761_c2.brep] c2