1 // Created on: 1992-05-07
2 // Created by: Jacques GOUSSARD
3 // Copyright (c) 1992-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
17 #include <IntAna_ListOfCurve.hxx>
18 #include <IntAna_ListIteratorOfListOfCurve.hxx>
19 #include <IntPatch_WLine.hxx>
21 #include <math_Matrix.hxx>
23 //If Abs(a) <= aNulValue then it is considered that a = 0.
24 static const Standard_Real aNulValue = 1.0e-11;
29 Standard_Boolean ExploreCurve(const gp_Cylinder& aCy,
32 const Standard_Real aTol,
33 IntAna_ListOfCurve& aLC);
35 Standard_Boolean IsToReverse(const gp_Cylinder& Cy1,
36 const gp_Cylinder& Cy2,
37 const Standard_Real Tol);
39 static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
40 const Standard_Real theUlTarget,
41 Standard_Real& theUGiven,
42 const Standard_Real theTol2D,
43 const Standard_Real thePeriod,
44 const Standard_Boolean theFlForce);
46 static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
47 const IntSurf_Quadric& theQuad2,
48 const Handle(IntSurf_LineOn2S)& theLine,
49 const stCoeffsValue& theCoeffs,
50 const Standard_Integer theWLIndex,
51 const Standard_Integer theMinNbPoints,
52 const Standard_Integer theStartPointOnLine,
53 const Standard_Integer theEndPointOnLine,
54 const Standard_Real theU2f,
55 const Standard_Real theU2l,
56 const Standard_Real theTol2D,
57 const Standard_Real thePeriodOfSurf2,
58 const Standard_Boolean isTheReverse);
60 //=======================================================================
62 //purpose : Replaces theParMIN = MIN(theParMIN, theParMAX),
63 // theParMAX = MAX(theParMIN, theParMAX).
64 //=======================================================================
65 static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX)
67 if(theParMIN > theParMAX)
69 const Standard_Real aux = theParMAX;
70 theParMAX = theParMIN;
75 //=======================================================================
76 //function : VBoundaryPrecise
77 //purpose : By default, we shall consider, that V1 and V2 will increase
78 // if U1 increases. But if it is not, new V1set and/or V2set
79 // must be computed as [V1current - DeltaV1] (analogically
80 // for V2). This function processes this case.
81 //=======================================================================
82 static void VBoundaryPrecise( const math_Matrix& theMatr,
83 const Standard_Real theV1AfterDecrByDelta,
84 const Standard_Real theV2AfterDecrByDelta,
85 const Standard_Real theV1f,
86 const Standard_Real theV2f,
87 Standard_Real& theV1Set,
88 Standard_Real& theV2Set)
90 //Now we are going to define if V1 (and V2) increases
91 //(or decreases) when U1 will increase.
92 const Standard_Integer aNbDim = 3;
93 math_Matrix aSyst(1, aNbDim, 1, aNbDim);
95 aSyst.SetCol(1, theMatr.Col(1));
96 aSyst.SetCol(2, theMatr.Col(2));
97 aSyst.SetCol(3, theMatr.Col(4));
99 const Standard_Real aDet = aSyst.Determinant();
101 aSyst.SetCol(1, theMatr.Col(3));
102 const Standard_Real aDet1 = aSyst.Determinant();
104 aSyst.SetCol(1, theMatr.Col(1));
105 aSyst.SetCol(2, theMatr.Col(3));
107 const Standard_Real aDet2 = aSyst.Determinant();
111 theV1Set = Max(theV1AfterDecrByDelta, theV1f);
116 theV2Set = Max(theV2AfterDecrByDelta, theV2f);
120 //=======================================================================
121 //function : DeltaU1Computing
122 //purpose : Computes new step for U1 parameter.
123 //=======================================================================
125 Standard_Boolean DeltaU1Computing(const math_Matrix& theSyst,
126 const math_Vector& theFree,
127 Standard_Real& theDeltaU1Found)
129 Standard_Real aDet = theSyst.Determinant();
131 if(Abs(aDet) > aNulValue)
133 math_Matrix aSyst1(theSyst);
134 aSyst1.SetCol(2, theFree);
136 theDeltaU1Found = Abs(aSyst1.Determinant()/aDet);
137 return Standard_True;
140 return Standard_False;
143 //=======================================================================
144 //function : StepComputing
148 // theMatr must have 3*5-dimension strictly.
150 // {a11*V1+a12*V2+a13*dU1+a14*dU2=b1;
151 // {a21*V1+a22*V2+a23*dU1+a24*dU2=b2;
152 // {a31*V1+a32*V2+a33*dU1+a34*dU2=b3;
153 // theMatr must be following:
154 // (a11 a12 a13 a14 b1)
155 // (a21 a22 a23 a24 b2)
156 // (a31 a32 a33 a34 b3)
157 //=======================================================================
158 static Standard_Boolean StepComputing(const math_Matrix& theMatr,
159 const Standard_Real theV1Cur,
160 const Standard_Real theV2Cur,
161 const Standard_Real theDeltaV1,
162 const Standard_Real theDeltaV2,
163 const Standard_Real theV1f,
164 const Standard_Real theV1l,
165 const Standard_Real theV2f,
166 const Standard_Real theV2l,
167 Standard_Real& theDeltaU1Found/*,
168 Standard_Real& theDeltaU2Found,
169 Standard_Real& theV1Found,
170 Standard_Real& theV2Found*/)
172 Standard_Boolean isSuccess = Standard_False;
173 theDeltaU1Found/* = theDeltaU2Found*/ = RealLast();
174 //theV1Found = theV1set;
175 //theV2Found = theV2Set;
176 const Standard_Integer aNbDim = 3;
178 math_Matrix aSyst(1, aNbDim, 1, aNbDim);
179 math_Vector aFree(1, aNbDim);
181 //By default, increasing V1(U1) and V2(U1) functions is
183 Standard_Real aV1Set = Min(theV1Cur + theDeltaV1, theV1l),
184 aV2Set = Min(theV2Cur + theDeltaV2, theV2l);
186 //However, what is indead?
187 VBoundaryPrecise( theMatr, theV1Cur - theDeltaV1, theV2Cur - theDeltaV2,
188 theV1f, theV2f, aV1Set, aV2Set);
190 aSyst.SetCol(2, theMatr.Col(3));
191 aSyst.SetCol(3, theMatr.Col(4));
193 for(Standard_Integer i = 0; i < 2; i++)
197 aSyst.SetCol(1, theMatr.Col(2));
198 aFree.Set(1, aNbDim, theMatr.Col(5)-aV1Set*theMatr.Col(1));
201 {//i==1 => V2 is known
202 aSyst.SetCol(1, theMatr.Col(1));
203 aFree.Set(1, aNbDim, theMatr.Col(5)-aV2Set*theMatr.Col(2));
206 Standard_Real aNewDU = theDeltaU1Found;
207 if(DeltaU1Computing(aSyst, aFree, aNewDU))
209 isSuccess = Standard_True;
210 if(aNewDU < theDeltaU1Found)
212 theDeltaU1Found = aNewDU;
219 aFree = theMatr.Col(5) - aV1Set*theMatr.Col(1) - aV2Set*theMatr.Col(2);
220 math_Matrix aSyst1(1, aNbDim, 1, 2);
221 aSyst1.SetCol(1, aSyst.Col(2));
222 aSyst1.SetCol(2, aSyst.Col(3));
224 //Now we have overdetermined system.
226 const Standard_Real aDet1 = theMatr(1,3)*theMatr(2,4) - theMatr(2,3)*theMatr(1,4);
227 const Standard_Real aDet2 = theMatr(1,3)*theMatr(3,4) - theMatr(3,3)*theMatr(1,4);
228 const Standard_Real aDet3 = theMatr(2,3)*theMatr(3,4) - theMatr(3,3)*theMatr(2,4);
229 const Standard_Real anAbsD1 = Abs(aDet1);
230 const Standard_Real anAbsD2 = Abs(aDet2);
231 const Standard_Real anAbsD3 = Abs(aDet3);
233 if(anAbsD1 >= anAbsD2)
235 if(anAbsD1 >= anAbsD3)
238 if(anAbsD1 <= aNulValue)
241 theDeltaU1Found = Abs(aFree(1)*theMatr(2,4) - aFree(2)*theMatr(1,4))/anAbsD1;
242 isSuccess = Standard_True;
247 if(anAbsD3 <= aNulValue)
250 theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
251 isSuccess = Standard_True;
256 if(anAbsD2 >= anAbsD3)
259 if(anAbsD2 <= aNulValue)
262 theDeltaU1Found = Abs(aFree(1)*theMatr(3,4) - aFree(3)*theMatr(1,4))/anAbsD2;
263 isSuccess = Standard_True;
268 if(anAbsD3 <= aNulValue)
271 theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
272 isSuccess = Standard_True;
280 //=======================================================================
281 //function : ProcessBounds
283 //=======================================================================
284 void ProcessBounds(const Handle(IntPatch_ALine)& alig, //-- ligne courante
285 const IntPatch_SequenceOfLine& slin,
286 const IntSurf_Quadric& Quad1,
287 const IntSurf_Quadric& Quad2,
288 Standard_Boolean& procf,
289 const gp_Pnt& ptf, //-- Debut Ligne Courante
290 const Standard_Real first, //-- Paramf
291 Standard_Boolean& procl,
292 const gp_Pnt& ptl, //-- Fin Ligne courante
293 const Standard_Real last, //-- Paraml
294 Standard_Boolean& Multpoint,
295 const Standard_Real Tol)
297 Standard_Integer j,k;
298 Standard_Real U1,V1,U2,V2;
299 IntPatch_Point ptsol;
302 if (procf && procl) {
303 j = slin.Length() + 1;
310 //-- On prend les lignes deja enregistrees
312 while (j <= slin.Length()) {
313 if(slin.Value(j)->ArcType() == IntPatch_Analytic) {
314 const Handle(IntPatch_ALine)& aligold = *((Handle(IntPatch_ALine)*)&slin.Value(j));
317 //-- On prend les vertex des lignes deja enregistrees
319 while (k <= aligold->NbVertex()) {
320 ptsol = aligold->Vertex(k);
322 d=ptf.Distance(ptsol.Value());
324 if (!ptsol.IsMultiple()) {
325 //-- le point ptsol (de aligold) est declare multiple sur aligold
326 Multpoint = Standard_True;
327 ptsol.SetMultiple(Standard_True);
328 aligold->Replace(k,ptsol);
330 ptsol.SetParameter(first);
331 alig->AddVertex(ptsol);
332 alig->SetFirstPoint(alig->NbVertex());
333 procf = Standard_True;
335 //-- On restore le point avec son parametre sur aligold
336 ptsol = aligold->Vertex(k);
341 if (ptl.Distance(ptsol.Value()) <= Tol) {
342 if (!ptsol.IsMultiple()) {
343 Multpoint = Standard_True;
344 ptsol.SetMultiple(Standard_True);
345 aligold->Replace(k,ptsol);
347 ptsol.SetParameter(last);
348 alig->AddVertex(ptsol);
349 alig->SetLastPoint(alig->NbVertex());
350 procl = Standard_True;
352 //-- On restore le point avec son parametre sur aligold
353 ptsol = aligold->Vertex(k);
357 if (procf && procl) {
358 k = aligold->NbVertex()+1;
364 if (procf && procl) {
372 if (!procf && !procl) {
373 Quad1.Parameters(ptf,U1,V1);
374 Quad2.Parameters(ptf,U2,V2);
375 ptsol.SetValue(ptf,Tol,Standard_False);
376 ptsol.SetParameters(U1,V1,U2,V2);
377 ptsol.SetParameter(first);
378 if (ptf.Distance(ptl) <= Tol) {
379 ptsol.SetMultiple(Standard_True); // a voir
380 Multpoint = Standard_True; // a voir de meme
381 alig->AddVertex(ptsol);
382 alig->SetFirstPoint(alig->NbVertex());
384 ptsol.SetParameter(last);
385 alig->AddVertex(ptsol);
386 alig->SetLastPoint(alig->NbVertex());
389 alig->AddVertex(ptsol);
390 alig->SetFirstPoint(alig->NbVertex());
391 Quad1.Parameters(ptl,U1,V1);
392 Quad2.Parameters(ptl,U2,V2);
393 ptsol.SetValue(ptl,Tol,Standard_False);
394 ptsol.SetParameters(U1,V1,U2,V2);
395 ptsol.SetParameter(last);
396 alig->AddVertex(ptsol);
397 alig->SetLastPoint(alig->NbVertex());
401 Quad1.Parameters(ptf,U1,V1);
402 Quad2.Parameters(ptf,U2,V2);
403 ptsol.SetValue(ptf,Tol,Standard_False);
404 ptsol.SetParameters(U1,V1,U2,V2);
405 ptsol.SetParameter(first);
406 alig->AddVertex(ptsol);
407 alig->SetFirstPoint(alig->NbVertex());
410 Quad1.Parameters(ptl,U1,V1);
411 Quad2.Parameters(ptl,U2,V2);
412 ptsol.SetValue(ptl,Tol,Standard_False);
413 ptsol.SetParameters(U1,V1,U2,V2);
414 ptsol.SetParameter(last);
415 alig->AddVertex(ptsol);
416 alig->SetLastPoint(alig->NbVertex());
419 //=======================================================================
422 //=======================================================================
423 Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
424 const IntSurf_Quadric& Quad2,
425 const Standard_Real Tol,
426 Standard_Boolean& Empty,
427 Standard_Boolean& Same,
428 Standard_Boolean& Multpoint,
429 IntPatch_SequenceOfLine& slin,
430 IntPatch_SequenceOfPoint& spnt)
433 IntPatch_Point ptsol;
437 IntSurf_TypeTrans trans1,trans2;
438 IntAna_ResultType typint;
443 gp_Cylinder Cy1(Quad1.Cylinder());
444 gp_Cylinder Cy2(Quad2.Cylinder());
446 IntAna_QuadQuadGeo inter(Cy1,Cy2,Tol);
450 return Standard_False;
453 typint = inter.TypeInter();
454 Standard_Integer NbSol = inter.NbSolutions();
455 Empty = Standard_False;
456 Same = Standard_False;
462 Empty = Standard_True;
468 Same = Standard_True;
474 gp_Pnt psol(inter.Point(1));
475 Standard_Real U1,V1,U2,V2;
476 Quad1.Parameters(psol,U1,V1);
477 Quad2.Parameters(psol,U2,V2);
478 ptsol.SetValue(psol,Tol,Standard_True);
479 ptsol.SetParameters(U1,V1,U2,V2);
488 { // Cylinders are tangent to each other by line
489 linsol = inter.Line(1);
490 ptref = linsol.Location();
491 gp_Dir crb1(gp_Vec(ptref,Cy1.Location()));
492 gp_Dir crb2(gp_Vec(ptref,Cy2.Location()));
493 gp_Vec norm1(Quad1.Normale(ptref));
494 gp_Vec norm2(Quad2.Normale(ptref));
495 IntSurf_Situation situcyl1;
496 IntSurf_Situation situcyl2;
498 if (crb1.Dot(crb2) < 0.)
499 { // centre de courbures "opposes"
500 if (norm1.Dot(crb1) > 0.)
502 situcyl2 = IntSurf_Inside;
506 situcyl2 = IntSurf_Outside;
509 if (norm2.Dot(crb2) > 0.)
511 situcyl1 = IntSurf_Inside;
515 situcyl1 = IntSurf_Outside;
520 if (Cy1.Radius() < Cy2.Radius())
522 if (norm1.Dot(crb1) > 0.)
524 situcyl2 = IntSurf_Inside;
528 situcyl2 = IntSurf_Outside;
531 if (norm2.Dot(crb2) > 0.)
533 situcyl1 = IntSurf_Outside;
537 situcyl1 = IntSurf_Inside;
542 if (norm1.Dot(crb1) > 0.)
544 situcyl2 = IntSurf_Outside;
548 situcyl2 = IntSurf_Inside;
551 if (norm2.Dot(crb2) > 0.)
553 situcyl1 = IntSurf_Inside;
557 situcyl1 = IntSurf_Outside;
562 Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_True, situcyl1, situcyl2);
567 for (i=1; i <= NbSol; i++)
569 linsol = inter.Line(i);
570 ptref = linsol.Location();
571 gp_Vec lsd = linsol.Direction();
572 Standard_Real qwe = lsd.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
575 trans1 = IntSurf_Out;
578 else if (qwe <-0.00000001)
581 trans2 = IntSurf_Out;
585 trans1=trans2=IntSurf_Undecided;
588 Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
599 IntPatch_Point pmult1, pmult2;
601 elipsol = inter.Ellipse(1);
603 gp_Pnt pttang1(ElCLib::Value(0.5*M_PI, elipsol));
604 gp_Pnt pttang2(ElCLib::Value(1.5*M_PI, elipsol));
606 Multpoint = Standard_True;
607 pmult1.SetValue(pttang1,Tol,Standard_True);
608 pmult2.SetValue(pttang2,Tol,Standard_True);
609 pmult1.SetMultiple(Standard_True);
610 pmult2.SetMultiple(Standard_True);
612 Standard_Real oU1,oV1,oU2,oV2;
613 Quad1.Parameters(pttang1,oU1,oV1);
614 Quad2.Parameters(pttang1,oU2,oV2);
615 pmult1.SetParameters(oU1,oV1,oU2,oV2);
617 Quad1.Parameters(pttang2,oU1,oV1);
618 Quad2.Parameters(pttang2,oU2,oV2);
619 pmult2.SetParameters(oU1,oV1,oU2,oV2);
621 // on traite la premiere ellipse
623 //-- Calcul de la Transition de la ligne
624 ElCLib::D1(0.,elipsol,ptref,Tgt);
625 Standard_Real qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
628 trans1 = IntSurf_Out;
631 else if (qwe<-0.00000001)
634 trans2 = IntSurf_Out;
638 trans1=trans2=IntSurf_Undecided;
641 //-- Transition calculee au point 0 -> Trans2 , Trans1
642 //-- car ici, on devarit calculer en PI
643 Handle(IntPatch_GLine) glig = new IntPatch_GLine(elipsol,Standard_False,trans2,trans1);
646 Standard_Real aU1, aV1, aU2, aV2;
648 gp_Pnt aP (ElCLib::Value(0., elipsol));
650 aIP.SetValue(aP,Tol,Standard_False);
651 aIP.SetMultiple(Standard_False);
653 Quad1.Parameters(aP, aU1, aV1);
654 Quad2.Parameters(aP, aU2, aV2);
655 aIP.SetParameters(aU1, aV1, aU2, aV2);
657 aIP.SetParameter(0.);
658 glig->AddVertex(aIP);
659 glig->SetFirstPoint(1);
661 aIP.SetParameter(2.*M_PI);
662 glig->AddVertex(aIP);
663 glig->SetLastPoint(2);
666 pmult1.SetParameter(0.5*M_PI);
667 glig->AddVertex(pmult1);
669 pmult2.SetParameter(1.5*M_PI);
670 glig->AddVertex(pmult2);
675 //-- Transitions calculee au point 0 OK
677 // on traite la deuxieme ellipse
678 elipsol = inter.Ellipse(2);
680 Standard_Real param1 = ElCLib::Parameter(elipsol,pttang1);
681 Standard_Real param2 = ElCLib::Parameter(elipsol,pttang2);
682 Standard_Real parampourtransition = 0.0;
685 pmult1.SetParameter(0.5*M_PI);
686 pmult2.SetParameter(1.5*M_PI);
687 parampourtransition = M_PI;
690 pmult1.SetParameter(1.5*M_PI);
691 pmult2.SetParameter(0.5*M_PI);
692 parampourtransition = 0.0;
695 //-- Calcul des transitions de ligne pour la premiere ligne
696 ElCLib::D1(parampourtransition,elipsol,ptref,Tgt);
697 qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
700 trans1 = IntSurf_Out;
703 else if(qwe< -0.00000001)
706 trans2 = IntSurf_Out;
710 trans1=trans2=IntSurf_Undecided;
713 //-- La transition a ete calculee sur un point de cette ligne
714 glig = new IntPatch_GLine(elipsol,Standard_False,trans1,trans2);
717 Standard_Real aU1, aV1, aU2, aV2;
719 gp_Pnt aP (ElCLib::Value(0., elipsol));
721 aIP.SetValue(aP,Tol,Standard_False);
722 aIP.SetMultiple(Standard_False);
724 Quad1.Parameters(aP, aU1, aV1);
725 Quad2.Parameters(aP, aU2, aV2);
726 aIP.SetParameters(aU1, aV1, aU2, aV2);
728 aIP.SetParameter(0.);
729 glig->AddVertex(aIP);
730 glig->SetFirstPoint(1);
732 aIP.SetParameter(2.*M_PI);
733 glig->AddVertex(aIP);
734 glig->SetLastPoint(2);
737 glig->AddVertex(pmult1);
738 glig->AddVertex(pmult2);
744 case IntAna_NoGeometricSolution:
746 Standard_Boolean bReverse;
747 Standard_Real U1,V1,U2,V2;
750 bReverse=IsToReverse(Cy1, Cy2, Tol);
753 Cy2=Quad1.Cylinder();
754 Cy1=Quad2.Cylinder();
757 IntAna_IntQuadQuad anaint(Cy1,Cy2,Tol);
758 if (!anaint.IsDone())
760 return Standard_False;
763 if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0)
765 Empty = Standard_True;
769 NbSol = anaint.NbPnt();
770 for (i = 1; i <= NbSol; i++)
772 psol = anaint.Point(i);
773 Quad1.Parameters(psol,U1,V1);
774 Quad2.Parameters(psol,U2,V2);
775 ptsol.SetValue(psol,Tol,Standard_True);
776 ptsol.SetParameters(U1,V1,U2,V2);
780 gp_Pnt ptvalid, ptf, ptl;
783 Standard_Real first,last,para;
784 IntAna_Curve curvsol;
785 Standard_Boolean tgfound;
786 Standard_Integer kount;
788 NbSol = anaint.NbCurve();
789 for (i = 1; i <= NbSol; i++)
791 curvsol = anaint.Curve(i);
792 curvsol.Domain(first,last);
793 ptf = curvsol.Value(first);
794 ptl = curvsol.Value(last);
798 tgfound = Standard_False;
802 para = (1.123*first + para)/2.123;
803 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
811 Handle(IntPatch_ALine) alig;
814 Standard_Real qwe = tgvalid.DotCross( Quad2.Normale(ptvalid),
815 Quad1.Normale(ptvalid));
818 trans1 = IntSurf_Out;
821 else if(qwe<-0.00000001)
824 trans2 = IntSurf_Out;
828 trans1=trans2=IntSurf_Undecided;
830 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
834 alig = new IntPatch_ALine(curvsol,Standard_False);
835 //-- cout << "Transition indeterminee" << endl;
838 Standard_Boolean TempFalse1 = Standard_False;
839 Standard_Boolean TempFalse2 = Standard_False;
841 ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1,ptf,first,
842 TempFalse2,ptl,last,Multpoint,Tol);
850 return Standard_False;
853 return Standard_True;
856 //=======================================================================
857 //function : ShortCosForm
858 //purpose : Represents theCosFactor*cosA+theSinFactor*sinA as
859 // theCoeff*cos(A-theAngle) if it is possibly (all angles are
861 //=======================================================================
862 static void ShortCosForm( const Standard_Real theCosFactor,
863 const Standard_Real theSinFactor,
864 Standard_Real& theCoeff,
865 Standard_Real& theAngle)
867 theCoeff = sqrt(theCosFactor*theCosFactor+theSinFactor*theSinFactor);
869 if(IsEqual(theCoeff, 0.0))
875 theAngle = acos(Abs(theCosFactor/theCoeff));
877 if(theSinFactor > 0.0)
879 if(IsEqual(theCosFactor, 0.0))
883 else if(theCosFactor < 0.0)
885 theAngle = M_PI-theAngle;
888 else if(IsEqual(theSinFactor, 0.0))
890 if(theCosFactor < 0.0)
895 if(theSinFactor < 0.0)
897 if(theCosFactor > 0.0)
899 theAngle = 2.0*M_PI-theAngle;
901 else if(IsEqual(theCosFactor, 0.0))
903 theAngle = 3.0*M_PI/2.0;
905 else if(theCosFactor < 0.0)
907 theAngle = M_PI+theAngle;
919 //Stores equations coefficients
922 stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&);
932 Standard_Real mK21; //sinU2
933 Standard_Real mK11; //sinU1
934 Standard_Real mL21; //cosU2
935 Standard_Real mL11; //cosU1
936 Standard_Real mM1; //Free member
938 Standard_Real mK22; //sinU2
939 Standard_Real mK12; //sinU1
940 Standard_Real mL22; //cosU2
941 Standard_Real mL12; //cosU1
942 Standard_Real mM2; //Free member
950 Standard_Real mPSIV1;
952 Standard_Real mPSIV2;
960 stCoeffsValue::stCoeffsValue( const gp_Cylinder& theCyl1,
961 const gp_Cylinder& theCyl2):
962 mVecA1(-theCyl1.Radius()*theCyl1.XAxis().Direction().XYZ()),
963 mVecA2(theCyl2.Radius()*theCyl2.XAxis().Direction().XYZ()),
964 mVecB1(-theCyl1.Radius()*theCyl1.YAxis().Direction().XYZ()),
965 mVecB2(theCyl2.Radius()*theCyl2.YAxis().Direction().XYZ()),
966 mVecC1(theCyl1.Axis().Direction().XYZ()),
967 mVecC2(theCyl2.Axis().Direction().XYZ().Reversed()),
968 mVecD(theCyl2.Location().XYZ() - theCyl1.Location().XYZ())
970 enum CoupleOfEquation
976 }aFoundCouple = COENONE;
979 Standard_Real aDetV1V2 = 0.0;
981 const Standard_Real aDelta1 = mVecC1(1)*mVecC2(2)-mVecC1(2)*mVecC2(1); //1-2
982 const Standard_Real aDelta2 = mVecC1(2)*mVecC2(3)-mVecC1(3)*mVecC2(2); //2-3
983 const Standard_Real aDelta3 = mVecC1(1)*mVecC2(3)-mVecC1(3)*mVecC2(1); //1-3
984 const Standard_Real anAbsD1 = Abs(aDelta1); //1-2
985 const Standard_Real anAbsD2 = Abs(aDelta2); //2-3
986 const Standard_Real anAbsD3 = Abs(aDelta3); //1-3
988 if(anAbsD1 >= anAbsD2)
990 if(anAbsD3 > anAbsD1)
992 aFoundCouple = COE13;
997 aFoundCouple = COE12;
1003 if(anAbsD3 > anAbsD2)
1005 aFoundCouple = COE13;
1010 aFoundCouple = COE23;
1015 if(Abs(aDetV1V2) < aNulValue)
1017 Standard_Failure::Raise("Error. Exception in divide by zerro (IntCyCyTrim)!!!!");
1020 switch(aFoundCouple)
1026 math_Vector aVTemp(mVecA1);
1027 mVecA1(1) = aVTemp(2);
1028 mVecA1(2) = aVTemp(3);
1029 mVecA1(3) = aVTemp(1);
1032 mVecA2(1) = aVTemp(2);
1033 mVecA2(2) = aVTemp(3);
1034 mVecA2(3) = aVTemp(1);
1037 mVecB1(1) = aVTemp(2);
1038 mVecB1(2) = aVTemp(3);
1039 mVecB1(3) = aVTemp(1);
1042 mVecB2(1) = aVTemp(2);
1043 mVecB2(2) = aVTemp(3);
1044 mVecB2(3) = aVTemp(1);
1047 mVecC1(1) = aVTemp(2);
1048 mVecC1(2) = aVTemp(3);
1049 mVecC1(3) = aVTemp(1);
1052 mVecC2(1) = aVTemp(2);
1053 mVecC2(2) = aVTemp(3);
1054 mVecC2(3) = aVTemp(1);
1057 mVecD(1) = aVTemp(2);
1058 mVecD(2) = aVTemp(3);
1059 mVecD(3) = aVTemp(1);
1065 math_Vector aVTemp = mVecA1;
1066 mVecA1(2) = aVTemp(3);
1067 mVecA1(3) = aVTemp(2);
1070 mVecA2(2) = aVTemp(3);
1071 mVecA2(3) = aVTemp(2);
1074 mVecB1(2) = aVTemp(3);
1075 mVecB1(3) = aVTemp(2);
1078 mVecB2(2) = aVTemp(3);
1079 mVecB2(3) = aVTemp(2);
1082 mVecC1(2) = aVTemp(3);
1083 mVecC1(3) = aVTemp(2);
1086 mVecC2(2) = aVTemp(3);
1087 mVecC2(3) = aVTemp(2);
1090 mVecD(2) = aVTemp(3);
1091 mVecD(3) = aVTemp(2);
1099 //------- For V1 (begin)
1101 mK21 = (mVecC2(2)*mVecB2(1)-mVecC2(1)*mVecB2(2))/aDetV1V2;
1103 mK11 = (mVecC2(2)*mVecB1(1)-mVecC2(1)*mVecB1(2))/aDetV1V2;
1105 mL21 = (mVecC2(2)*mVecA2(1)-mVecC2(1)*mVecA2(2))/aDetV1V2;
1107 mL11 = (mVecC2(2)*mVecA1(1)-mVecC2(1)*mVecA1(2))/aDetV1V2;
1109 mM1 = (mVecC2(2)*mVecD(1)-mVecC2(1)*mVecD(2))/aDetV1V2;
1110 //------- For V1 (end)
1112 //------- For V2 (begin)
1114 mK22 = (mVecC1(1)*mVecB2(2)-mVecC1(2)*mVecB2(1))/aDetV1V2;
1116 mK12 = (mVecC1(1)*mVecB1(2)-mVecC1(2)*mVecB1(1))/aDetV1V2;
1118 mL22 = (mVecC1(1)*mVecA2(2)-mVecC1(2)*mVecA2(1))/aDetV1V2;
1120 mL12 = (mVecC1(1)*mVecA1(2)-mVecC1(2)*mVecA1(1))/aDetV1V2;
1122 mM2 = (mVecC1(1)*mVecD(2)-mVecC1(2)*mVecD(1))/aDetV1V2;
1123 //------- For V1 (end)
1125 ShortCosForm(mL11, mK11, mK1, mFIV1);
1126 ShortCosForm(mL21, mK21, mL1, mPSIV1);
1127 ShortCosForm(mL12, mK12, mK2, mFIV2);
1128 ShortCosForm(mL22, mK22, mL2, mPSIV2);
1130 const Standard_Real aA1=mVecC1(3)*mK21+mVecC2(3)*mK22-mVecB2(3), //sinU2
1131 aA2=mVecC1(3)*mL21+mVecC2(3)*mL22-mVecA2(3), //cosU2
1132 aB1=mVecB1(3)-mVecC1(3)*mK11-mVecC2(3)*mK12, //sinU1
1133 aB2=mVecA1(3)-mVecC1(3)*mL11-mVecC2(3)*mL12; //cosU1
1135 mC =mVecD(3) - mVecC1(3)*mM1 -mVecC2(3)*mM2; //Free
1137 Standard_Real aA = 0.0;
1139 ShortCosForm(aB2,aB1,mB,mFI1);
1140 ShortCosForm(aA2,aA1,aA,mFI2);
1146 //=======================================================================
1147 //function : CylCylMonotonicity
1148 //purpose : Determines, if U2(U1) function is increasing.
1149 //=======================================================================
1150 static Standard_Boolean CylCylMonotonicity( const Standard_Real theU1par,
1151 const Standard_Integer theWLIndex,
1152 const stCoeffsValue& theCoeffs,
1153 const Standard_Real thePeriod,
1154 Standard_Boolean& theIsIncreasing)
1156 // U2(U1) = FI2 + (+/-)acos(B*cos(U1 - FI1) + C);
1158 //Func. U2(X1) = FI2 + X1;
1159 //Func. X1(X2) = anArccosFactor*X2;
1160 //Func. X2(X3) = acos(X3);
1161 //Func. X3(X4) = B*X4 + C;
1162 //Func. X4(U1) = cos(U1 - FI1).
1165 //U2(X1) is always increasing.
1166 //X1(X2) is increasing if anArccosFactor > 0.0 and
1167 //is decreasing otherwise.
1168 //X2(X3) is always decreasing.
1169 //Therefore, U2(X3) is decreasing if anArccosFactor > 0.0 and
1170 //is increasing otherwise.
1171 //X3(X4) is increasing if B > 0 and is decreasing otherwise.
1172 //X4(U1) is increasing if U1 - FI1 in [PI, 2*PI) and
1173 //is decreasing U1 - FI1 in [0, PI) or if (U1 - FI1 == 2PI).
1174 //After that, we can predict behaviour of U2(U1) function:
1175 //if it is increasing or decreasing.
1177 //For "+/-" sign. If isPlus == TRUE, "+" is chosen, otherwise, we choose "-".
1178 Standard_Boolean isPlus = Standard_False;
1183 isPlus = Standard_True;
1186 isPlus = Standard_False;
1189 //Standard_Failure::Raise("Error. Range Error!!!!");
1190 return Standard_False;
1193 Standard_Real aU1Temp = theU1par - theCoeffs.mFI1;
1194 InscribePoint(0, thePeriod, aU1Temp, 0.0, thePeriod, Standard_False);
1196 theIsIncreasing = Standard_True;
1198 if(((M_PI - aU1Temp) < RealSmall()) && (aU1Temp < thePeriod))
1200 theIsIncreasing = Standard_False;
1203 if(theCoeffs.mB < 0.0)
1205 theIsIncreasing = !theIsIncreasing;
1210 theIsIncreasing = !theIsIncreasing;
1213 return Standard_True;
1216 static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
1217 const Standard_Integer theWLIndex,
1218 const stCoeffsValue& theCoeffs,
1219 Standard_Real& theU2)
1221 Standard_Real aSign = 0.0;
1232 //Standard_Failure::Raise("Error. Range Error!!!!");
1233 return Standard_False;
1236 Standard_Real anArg = theCoeffs.mB *
1237 cos(theU1par - theCoeffs.mFI1) +
1240 if(aNulValue > 1.0 - anArg)
1243 if(anArg + 1.0 < aNulValue)
1246 theU2 = theCoeffs.mFI2 + aSign*acos(anArg);
1248 return Standard_True;
1251 static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1,
1252 const Standard_Real theU2,
1253 const stCoeffsValue& theCoeffs,
1254 Standard_Real& theV1,
1255 Standard_Real& theV2)
1257 theV1 = theCoeffs.mK21 * sin(theU2) +
1258 theCoeffs.mK11 * sin(theU1) +
1259 theCoeffs.mL21 * cos(theU2) +
1260 theCoeffs.mL11 * cos(theU1) + theCoeffs.mM1;
1262 theV2 = theCoeffs.mK22 * sin(theU2) +
1263 theCoeffs.mK12 * sin(theU1) +
1264 theCoeffs.mL22 * cos(theU2) +
1265 theCoeffs.mL12 * cos(theU1) + theCoeffs.mM2;
1267 return Standard_True;
1271 static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
1272 const Standard_Integer theWLIndex,
1273 const stCoeffsValue& theCoeffs,
1274 Standard_Real& theU2,
1275 Standard_Real& theV1,
1276 Standard_Real& theV2)
1278 if(!CylCylComputeParameters(theU1par, theWLIndex, theCoeffs, theU2))
1279 return Standard_False;
1281 if(!CylCylComputeParameters(theU1par, theU2, theCoeffs, theV1, theV2))
1282 return Standard_False;
1284 return Standard_True;
1287 //=======================================================================
1288 //function : SearchOnVBounds
1290 //=======================================================================
1291 static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType,
1292 const stCoeffsValue& theCoeffs,
1293 const Standard_Real theVzad,
1294 const Standard_Real theVInit,
1295 const Standard_Real theInitU2,
1296 const Standard_Real theInitMainVar,
1297 Standard_Real& theMainVariableValue)
1299 const Standard_Integer aNbDim = 3;
1300 const Standard_Real aMaxError = 4.0*M_PI; // two periods
1302 theMainVariableValue = theInitMainVar;
1303 const Standard_Real aTol2 = 1.0e-18;
1304 Standard_Real aMainVarPrev = theInitMainVar, aU2Prev = theInitU2, anOtherVar = theVInit;
1306 //Structure of aMatr:
1307 // C_{1}*U_{1} & C_{2}*U_{2} & C_{3}*V_{*},
1308 //where C_{1}, C_{2} and C_{3} are math_Vector.
1309 math_Matrix aMatr(1, aNbDim, 1, aNbDim);
1311 Standard_Real anError = RealLast();
1312 Standard_Real anErrorPrev = anError;
1313 Standard_Integer aNbIter = 0;
1316 if(++aNbIter > 1000)
1317 return Standard_False;
1319 const Standard_Real aSinU1 = sin(aMainVarPrev),
1320 aCosU1 = cos(aMainVarPrev),
1321 aSinU2 = sin(aU2Prev),
1322 aCosU2 = cos(aU2Prev);
1324 math_Vector aVecFreeMem = (theCoeffs.mVecA2 * aU2Prev +
1325 theCoeffs.mVecB2) * aSinU2 -
1326 (theCoeffs.mVecB2 * aU2Prev -
1327 theCoeffs.mVecA2) * aCosU2 +
1328 (theCoeffs.mVecA1 * aMainVarPrev +
1329 theCoeffs.mVecB1) * aSinU1 -
1330 (theCoeffs.mVecB1 * aMainVarPrev -
1331 theCoeffs.mVecA1) * aCosU1 +
1334 math_Vector aMSum(1, 3);
1339 aMatr.SetCol(3, theCoeffs.mVecC2);
1340 aMSum = theCoeffs.mVecC1 * theVzad;
1341 aVecFreeMem -= aMSum;
1342 aMSum += theCoeffs.mVecC2*anOtherVar;
1346 aMatr.SetCol(3, theCoeffs.mVecC1);
1347 aMSum = theCoeffs.mVecC2 * theVzad;
1348 aVecFreeMem -= aMSum;
1349 aMSum += theCoeffs.mVecC1*anOtherVar;
1353 return Standard_False;
1356 aMatr.SetCol(1, theCoeffs.mVecA1 * aSinU1 - theCoeffs.mVecB1 * aCosU1);
1357 aMatr.SetCol(2, theCoeffs.mVecA2 * aSinU2 - theCoeffs.mVecB2 * aCosU2);
1359 Standard_Real aDetMainSyst = aMatr.Determinant();
1361 if(Abs(aDetMainSyst) < aNulValue)
1363 return Standard_False;
1366 math_Matrix aM1(aMatr), aM2(aMatr), aM3(aMatr);
1367 aM1.SetCol(1, aVecFreeMem);
1368 aM2.SetCol(2, aVecFreeMem);
1369 aM3.SetCol(3, aVecFreeMem);
1371 const Standard_Real aDetMainVar = aM1.Determinant();
1372 const Standard_Real aDetVar1 = aM2.Determinant();
1373 const Standard_Real aDetVar2 = aM3.Determinant();
1375 Standard_Real aDelta = aDetMainVar/aDetMainSyst-aMainVarPrev;
1377 if(Abs(aDelta) > aMaxError)
1378 return Standard_False;
1380 anError = aDelta*aDelta;
1381 aMainVarPrev += aDelta;
1384 aDelta = aDetVar1/aDetMainSyst-aU2Prev;
1386 if(Abs(aDelta) > aMaxError)
1387 return Standard_False;
1389 anError += aDelta*aDelta;
1393 aDelta = aDetVar2/aDetMainSyst-anOtherVar;
1394 anError += aDelta*aDelta;
1395 anOtherVar += aDelta;
1397 if(anError > anErrorPrev)
1398 {//Method diverges. Keep the best result
1399 const Standard_Real aSinU1 = sin(aMainVarPrev),
1400 aCosU1 = cos(aMainVarPrev),
1401 aSinU2 = sin(aU2Prev),
1402 aCosU2 = cos(aU2Prev);
1403 aMSum -= (theCoeffs.mVecA1*aCosU1 +
1404 theCoeffs.mVecB1*aSinU1 +
1405 theCoeffs.mVecA2*aCosU2 +
1406 theCoeffs.mVecB2*aSinU2 +
1408 const Standard_Real aSQNorm = aMSum.Norm2();
1409 return (aSQNorm < aTol2);
1413 theMainVariableValue = aMainVarPrev;
1416 anErrorPrev = anError;
1418 while(anError > aTol2);
1420 theMainVariableValue = aMainVarPrev;
1422 return Standard_True;
1425 //=======================================================================
1426 //function : InscribePoint
1427 //purpose : If theFlForce==TRUE theUGiven will be changed forcefully
1428 // even if theUGiven is already inscibed in the boundary
1429 // (if it is possible; i.e. if new theUGiven is inscribed
1430 // in the boundary, too).
1431 //=======================================================================
1432 Standard_Boolean InscribePoint( const Standard_Real theUfTarget,
1433 const Standard_Real theUlTarget,
1434 Standard_Real& theUGiven,
1435 const Standard_Real theTol2D,
1436 const Standard_Real thePeriod,
1437 const Standard_Boolean theFlForce)
1439 if((theUfTarget - theUGiven <= theTol2D) &&
1440 (theUGiven - theUlTarget <= theTol2D))
1441 {//it has already inscribed
1450 Standard_Real anUtemp = theUGiven + thePeriod;
1451 if((theUfTarget - anUtemp <= theTol2D) &&
1452 (anUtemp - theUlTarget <= theTol2D))
1454 theUGiven = anUtemp;
1455 return Standard_True;
1458 anUtemp = theUGiven - thePeriod;
1459 if((theUfTarget - anUtemp <= theTol2D) &&
1460 (anUtemp - theUlTarget <= theTol2D))
1462 theUGiven = anUtemp;
1466 return Standard_True;
1469 if(IsEqual(thePeriod, 0.0))
1470 {//it cannot be inscribed
1471 return Standard_False;
1474 Standard_Real aDU = theUGiven - theUfTarget;
1481 while(((theUGiven - theUfTarget)*aDU < 0.0) && !((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D)))
1486 return ((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D));
1489 //=======================================================================
1490 //function : InscribeInterval
1491 //purpose : Adjusts theUfGiven and after that fits theUlGiven to result
1492 //=======================================================================
1493 static Standard_Boolean InscribeInterval(const Standard_Real theUfTarget,
1494 const Standard_Real theUlTarget,
1495 Standard_Real& theUfGiven,
1496 Standard_Real& theUlGiven,
1497 const Standard_Real theTol2D,
1498 const Standard_Real thePeriod)
1500 Standard_Real anUpar = theUfGiven;
1501 const Standard_Real aDelta = theUlGiven - theUfGiven;
1502 if(InscribePoint(theUfTarget, theUlTarget, anUpar,
1503 theTol2D, thePeriod, Standard_False))
1505 theUfGiven = anUpar;
1506 theUlGiven = theUfGiven + aDelta;
1510 anUpar = theUlGiven;
1511 if(InscribePoint(theUfTarget, theUlTarget, anUpar,
1512 theTol2D, thePeriod, Standard_False))
1514 theUlGiven = anUpar;
1515 theUfGiven = theUlGiven - aDelta;
1519 return Standard_False;
1523 return Standard_True;
1526 //=======================================================================
1527 //function : InscribeAndSortArray
1528 //purpose : Sort from Min to Max value
1529 //=======================================================================
1530 static void InscribeAndSortArray( Standard_Real theArr[],
1531 const Standard_Integer theNOfMembers,
1532 const Standard_Real theUf,
1533 const Standard_Real theUl,
1534 const Standard_Real theTol2D,
1535 const Standard_Real thePeriod)
1537 for(Standard_Integer i = 0; i < theNOfMembers; i++)
1539 if(Precision::IsInfinite(theArr[i]))
1542 theArr[i] = -theArr[i];
1547 InscribePoint(theUf, theUl, theArr[i], theTol2D, thePeriod, Standard_False);
1549 for(Standard_Integer j = i - 1; j >= 0; j--)
1552 if(theArr[j+1] < theArr[j])
1554 Standard_Real anUtemp = theArr[j+1];
1555 theArr[j+1] = theArr[j];
1556 theArr[j] = anUtemp;
1563 //=======================================================================
1564 //function : AddPointIntoWL
1565 //purpose : Surf1 is a surface, whose U-par is variable.
1566 //=======================================================================
1567 static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
1568 const IntSurf_Quadric& theQuad2,
1569 const stCoeffsValue& theCoeffs,
1570 const Standard_Boolean isTheReverse,
1571 const Standard_Boolean isThePrecise,
1572 const gp_Pnt2d& thePntOnSurf1,
1573 const gp_Pnt2d& thePntOnSurf2,
1574 const Standard_Real theUfSurf1,
1575 const Standard_Real theUlSurf1,
1576 const Standard_Real theUfSurf2,
1577 const Standard_Real theUlSurf2,
1578 const Standard_Real theVfSurf1,
1579 const Standard_Real theVlSurf1,
1580 const Standard_Real theVfSurf2,
1581 const Standard_Real theVlSurf2,
1582 const Standard_Real thePeriodOfSurf1,
1583 const Handle(IntSurf_LineOn2S)& theLine,
1584 const Standard_Integer theWLIndex,
1585 const Standard_Real theTol3D,
1586 const Standard_Real theTol2D,
1587 const Standard_Boolean theFlForce,
1588 const Standard_Boolean theFlBefore = Standard_False)
1590 gp_Pnt aPt1(theQuad1.Value(thePntOnSurf1.X(), thePntOnSurf1.Y())),
1591 aPt2(theQuad2.Value(thePntOnSurf2.X(), thePntOnSurf2.Y()));
1593 Standard_Real aU1par = thePntOnSurf1.X();
1594 if(!InscribePoint(theUfSurf1, theUlSurf1, aU1par, theTol2D,
1595 thePeriodOfSurf1, theFlForce))
1596 return Standard_False;
1598 Standard_Real aU2par = thePntOnSurf2.X();
1599 if(!InscribePoint(theUfSurf2, theUlSurf2, aU2par, theTol2D,
1600 thePeriodOfSurf1, Standard_False))
1601 return Standard_False;
1603 Standard_Real aV1par = thePntOnSurf1.Y();
1604 if((aV1par - theVlSurf1 > theTol2D) || (theVfSurf1 - aV1par > theTol2D))
1605 return Standard_False;
1607 Standard_Real aV2par = thePntOnSurf2.Y();
1608 if((aV2par - theVlSurf2 > theTol2D) || (theVfSurf2 - aV2par > theTol2D))
1609 return Standard_False;
1611 IntSurf_PntOn2S aPnt;
1615 aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
1621 aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
1626 Standard_Integer aNbPnts = theLine->NbPoints();
1629 Standard_Real aUl = 0.0, aVl = 0.0;
1630 const IntSurf_PntOn2S aPlast = theLine->Value(aNbPnts);
1632 aPlast.ParametersOnS2(aUl, aVl);
1634 aPlast.ParametersOnS1(aUl, aVl);
1636 if(!theFlBefore && (aU1par <= aUl))
1637 {//Parameter value must be increased if theFlBefore == FALSE.
1638 return Standard_False;
1641 //theTol2D is minimal step along parameter changed.
1642 //Therefore, if we apply this minimal step two
1643 //neighbour points will be always "same". Consequently,
1644 //we should reduce tolerance for IsSame checking.
1645 const Standard_Real aDTol = 1.0-Epsilon(1.0);
1646 if(aPnt.IsSame(aPlast, theTol3D*aDTol, theTol2D*aDTol))
1648 theLine->RemovePoint(aNbPnts);
1655 return Standard_True;
1657 //Try to precise existing WLine
1658 aNbPnts = theLine->NbPoints();
1661 Standard_Real aU1 = 0.0, aU2 = 0.0, aU3 = 0.0, aV = 0.0;
1664 theLine->Value(aNbPnts).ParametersOnS2(aU3, aV);
1665 theLine->Value(aNbPnts-1).ParametersOnS2(aU2, aV);
1666 theLine->Value(aNbPnts-2).ParametersOnS2(aU1, aV);
1670 theLine->Value(aNbPnts).ParametersOnS1(aU3, aV);
1671 theLine->Value(aNbPnts-1).ParametersOnS1(aU2, aV);
1672 theLine->Value(aNbPnts-2).ParametersOnS1(aU1, aV);
1675 const Standard_Real aStepPrev = aU2-aU1;
1676 const Standard_Real aStep = aU3-aU2;
1678 const Standard_Integer aDeltaStep = RealToInt(aStepPrev/aStep);
1680 if((1 < aDeltaStep) && (aDeltaStep < 2000))
1682 SeekAdditionalPoints( theQuad1, theQuad2, theLine,
1683 theCoeffs, theWLIndex, aDeltaStep, aNbPnts-2,
1684 aNbPnts-1, theUfSurf2, theUlSurf2,
1685 theTol2D, thePeriodOfSurf1, isTheReverse);
1689 return Standard_True;
1692 //=======================================================================
1693 //function : AddBoundaryPoint
1695 //=======================================================================
1696 static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1,
1697 const IntSurf_Quadric& theQuad2,
1698 const Handle(IntPatch_WLine)& theWL,
1699 const stCoeffsValue& theCoeffs,
1700 const Bnd_Box2d& theUVSurf1,
1701 const Bnd_Box2d& theUVSurf2,
1702 const Standard_Real theTol3D,
1703 const Standard_Real theTol2D,
1704 const Standard_Real thePeriod,
1705 const Standard_Real theU1,
1706 const Standard_Real theU2,
1707 const Standard_Real theV1,
1708 const Standard_Real theV1Prev,
1709 const Standard_Real theV2,
1710 const Standard_Real theV2Prev,
1711 const Standard_Integer theWLIndex,
1712 const Standard_Boolean isTheReverse,
1713 const Standard_Boolean theFlForce,
1714 Standard_Boolean& isTheFound1,
1715 Standard_Boolean& isTheFound2)
1717 Standard_Real aUSurf1f = 0.0, //const
1721 Standard_Real aUSurf2f = 0.0, //const
1726 theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
1727 theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
1729 SearchBoundType aTS1 = SearchNONE, aTS2 = SearchNONE;
1730 Standard_Real aV1zad = aVSurf1f, aV2zad = aVSurf2f;
1732 Standard_Real anUpar1 = theU1, anUpar2 = theU1;
1733 Standard_Real aVf = theV1, aVl = theV1Prev;
1735 if( (Abs(aVf-aVSurf1f) <= theTol2D) ||
1736 ((aVf-aVSurf1f)*(aVl-aVSurf1f) <= 0.0))
1740 isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1f, theV2, theU2, theU1, anUpar1);
1742 else if((Abs(aVf-aVSurf1l) <= theTol2D) ||
1743 ((aVf-aVSurf1l)*(aVl-aVSurf1l) <= 0.0))
1747 isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1l, theV2, theU2, theU1, anUpar1);
1753 if((Abs(aVf-aVSurf2f) <= theTol2D) ||
1754 ((aVf-aVSurf2f)*(aVl-aVSurf2f) <= 0.0))
1758 isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2f, theV1, theU2, theU1, anUpar2);
1760 else if((Abs(aVf-aVSurf2l) <= theTol2D) ||
1761 ((aVf-aVSurf2l)*(aVl-aVSurf2l) <= 0.0))
1765 isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2l, theV1, theU2, theU1, anUpar2);
1768 if(!isTheFound1 && !isTheFound2)
1769 return Standard_True;
1771 if(anUpar1 < anUpar2)
1775 Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
1776 if(!CylCylComputeParameters(anUpar1, theWLIndex, theCoeffs, aU2, aV1, aV2))
1778 isTheFound1 = Standard_False;
1781 if(aTS1 == SearchV1)
1784 if(aTS1 == SearchV2)
1787 if(isTheFound1 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False,
1788 gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2),
1789 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
1790 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod,
1791 theWL->Curve(), theWLIndex, theTol3D,
1792 theTol2D, theFlForce))
1794 isTheFound1 = Standard_False;
1800 Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
1801 if(!CylCylComputeParameters(anUpar2, theWLIndex, theCoeffs, aU2, aV1, aV2))
1803 isTheFound2 = Standard_False;
1806 if(aTS2 == SearchV1)
1809 if(aTS2 == SearchV2)
1812 if(isTheFound2 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False,
1813 gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2),
1814 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
1815 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod,
1816 theWL->Curve(), theWLIndex, theTol3D,
1817 theTol2D, theFlForce))
1819 isTheFound2 = Standard_False;
1827 Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
1828 if(!CylCylComputeParameters(anUpar2, theWLIndex, theCoeffs, aU2, aV1, aV2))
1830 isTheFound2 = Standard_False;
1833 if(aTS2 == SearchV1)
1836 if(aTS2 == SearchV2)
1839 if(isTheFound2 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False,
1840 gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2),
1841 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
1842 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod,
1843 theWL->Curve(), theWLIndex, theTol3D,
1844 theTol2D, theFlForce))
1846 isTheFound2 = Standard_False;
1852 Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
1853 if(!CylCylComputeParameters(anUpar1, theWLIndex, theCoeffs, aU2, aV1, aV2))
1855 isTheFound1 = Standard_False;
1858 if(aTS1 == SearchV1)
1861 if(aTS1 == SearchV2)
1864 if(isTheFound1 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False,
1865 gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2),
1866 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
1867 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod,
1868 theWL->Curve(), theWLIndex, theTol3D,
1869 theTol2D, theFlForce))
1871 isTheFound1 = Standard_False;
1876 return Standard_True;
1879 //=======================================================================
1880 //function : SeekAdditionalPoints
1882 //=======================================================================
1883 static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
1884 const IntSurf_Quadric& theQuad2,
1885 const Handle(IntSurf_LineOn2S)& theLine,
1886 const stCoeffsValue& theCoeffs,
1887 const Standard_Integer theWLIndex,
1888 const Standard_Integer theMinNbPoints,
1889 const Standard_Integer theStartPointOnLine,
1890 const Standard_Integer theEndPointOnLine,
1891 const Standard_Real theU2f,
1892 const Standard_Real theU2l,
1893 const Standard_Real theTol2D,
1894 const Standard_Real thePeriodOfSurf2,
1895 const Standard_Boolean isTheReverse)
1897 if(theLine.IsNull())
1900 Standard_Integer aNbPoints = theEndPointOnLine - theStartPointOnLine + 1;
1901 if(aNbPoints >= theMinNbPoints)
1906 Standard_Real aMinDeltaParam = theTol2D;
1909 Standard_Real u1 = 0.0, v1 = 0.0, u2 = 0.0, v2 = 0.0;
1913 theLine->Value(theStartPointOnLine).ParametersOnS2(u1, v1);
1914 theLine->Value(theEndPointOnLine).ParametersOnS2(u2, v2);
1918 theLine->Value(theStartPointOnLine).ParametersOnS1(u1, v1);
1919 theLine->Value(theEndPointOnLine).ParametersOnS1(u2, v2);
1922 aMinDeltaParam = Max(Abs(u2 - u1)/IntToReal(theMinNbPoints), aMinDeltaParam);
1925 Standard_Integer aLastPointIndex = theEndPointOnLine;
1926 Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0;
1928 Standard_Integer aNbPointsPrev = 0;
1929 while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
1931 aNbPointsPrev = aNbPoints;
1932 for(Standard_Integer fp = theStartPointOnLine, lp = 0; fp < aLastPointIndex; fp = lp + 1)
1934 Standard_Real U1f = 0.0, V1f = 0.0; //first point in 1st suraface
1935 Standard_Real U1l = 0.0, V1l = 0.0; //last point in 1st suraface
1937 Standard_Real U2f = 0.0, V2f = 0.0; //first point in 2nd suraface
1938 Standard_Real U2l = 0.0, V2l = 0.0; //last point in 2nd suraface
1944 theLine->Value(fp).ParametersOnS2(U1f, V1f);
1945 theLine->Value(lp).ParametersOnS2(U1l, V1l);
1947 theLine->Value(fp).ParametersOnS1(U2f, V2f);
1948 theLine->Value(lp).ParametersOnS1(U2l, V2l);
1952 theLine->Value(fp).ParametersOnS1(U1f, V1f);
1953 theLine->Value(lp).ParametersOnS1(U1l, V1l);
1955 theLine->Value(fp).ParametersOnS2(U2f, V2f);
1956 theLine->Value(lp).ParametersOnS2(U2l, V2l);
1959 if(Abs(U1l - U1f) <= aMinDeltaParam)
1961 //Step is minimal. It is not necessary to divide it.
1965 U1prec = 0.5*(U1f+U1l);
1967 if(!CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec))
1970 InscribePoint(theU2f, theU2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False);
1972 const gp_Pnt aP1(theQuad1.Value(U1prec, V1prec)), aP2(theQuad2.Value(U2prec, V2prec));
1973 const gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
1976 //cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << endl;
1979 IntSurf_PntOn2S anIP;
1982 anIP.SetValue(aPInt, U2prec, V2prec, U1prec, V1prec);
1986 anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec);
1989 theLine->InsertBefore(lp, anIP);
1995 if(aNbPoints >= theMinNbPoints)
2002 //=======================================================================
2003 //function : CriticalPointsComputing
2005 //=======================================================================
2006 static void CriticalPointsComputing(const stCoeffsValue& theCoeffs,
2007 const Standard_Real theUSurf1f,
2008 const Standard_Real theUSurf1l,
2009 const Standard_Real theUSurf2f,
2010 const Standard_Real theUSurf2l,
2011 const Standard_Real thePeriod,
2012 const Standard_Real theTol2D,
2013 const Standard_Integer theNbCritPointsMax,
2014 Standard_Real theU1crit[])
2017 theU1crit[1] = thePeriod;
2018 theU1crit[2] = theUSurf1f;
2019 theU1crit[3] = theUSurf1l;
2021 const Standard_Real aCOS = cos(theCoeffs.mFI2);
2022 const Standard_Real aBSB = Abs(theCoeffs.mB);
2023 if((theCoeffs.mC - aBSB <= aCOS) && (aCOS <= theCoeffs.mC + aBSB))
2025 Standard_Real anArg = (aCOS - theCoeffs.mC) / theCoeffs.mB;
2031 theU1crit[4] = -acos(anArg) + theCoeffs.mFI1;
2032 theU1crit[5] = acos(anArg) + theCoeffs.mFI1;
2035 Standard_Real aSf = cos(theUSurf2f - theCoeffs.mFI2);
2036 Standard_Real aSl = cos(theUSurf2l - theCoeffs.mFI2);
2039 theU1crit[6] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? -acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : -Precision::Infinite();
2040 theU1crit[7] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? -acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : Precision::Infinite();
2041 theU1crit[8] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : -Precision::Infinite();
2042 theU1crit[9] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : Precision::Infinite();
2044 //preparative treatment of array
2045 InscribeAndSortArray(theU1crit, theNbCritPointsMax, 0.0, thePeriod, theTol2D, thePeriod);
2046 for(Standard_Integer i = 1; i < theNbCritPointsMax; i++)
2048 Standard_Real &a = theU1crit[i],
2049 &b = theU1crit[i-1];
2050 const Standard_Real aRemain = fmod(Abs(a - b), thePeriod); // >= 0, because Abs(a - b) >= 0
2051 if((Abs(a - b) < theTol2D) || (aRemain < theTol2D) || (Abs(aRemain - thePeriod) < theTol2D))
2054 b = Precision::Infinite();
2059 //=======================================================================
2060 //function : IntCyCyTrim
2062 //=======================================================================
2063 Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
2064 const IntSurf_Quadric& theQuad2,
2065 const Standard_Real theTol3D,
2066 const Standard_Real theTol2D,
2067 const Bnd_Box2d& theUVSurf1,
2068 const Bnd_Box2d& theUVSurf2,
2069 const Standard_Boolean isTheReverse,
2070 Standard_Boolean& isTheEmpty,
2071 IntPatch_SequenceOfLine& theSlin,
2072 IntPatch_SequenceOfPoint& theSPnt)
2074 Standard_Real aUSurf1f = 0.0, //const
2078 Standard_Real aUSurf2f = 0.0, //const
2083 theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
2084 theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
2086 const gp_Cylinder& aCyl1 = theQuad1.Cylinder(),
2087 aCyl2 = theQuad2.Cylinder();
2089 IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
2091 if (!anInter.IsDone())
2093 return Standard_False;
2096 IntAna_ResultType aTypInt = anInter.TypeInter();
2098 if(aTypInt != IntAna_NoGeometricSolution)
2099 { //It is not necessary (because result is an analytic curve) or
2100 //it is impossible to make Walking-line.
2102 return Standard_False;
2107 const Standard_Integer aNbMaxPoints = 2000;
2108 const Standard_Integer aNbMinPoints = 200;
2109 const Standard_Integer aNbPoints = Min(Max(aNbMinPoints,
2110 RealToInt(20.0*aCyl1.Radius())), aNbMaxPoints);
2111 const Standard_Real aPeriod = 2.0*M_PI;
2112 const Standard_Real aStepMin = theTol2D,
2113 aStepMax = (aUSurf1l-aUSurf1f)/IntToReal(aNbPoints);
2114 const Standard_Integer aNbWLines = 2;
2116 const stCoeffsValue anEquationCoeffs(aCyl1, aCyl2);
2119 const Standard_Integer aNbOfBoundaries = 2;
2120 Standard_Real aU1f[aNbOfBoundaries] = {-Precision::Infinite(), -Precision::Infinite()};
2121 Standard_Real aU1l[aNbOfBoundaries] = {Precision::Infinite(), Precision::Infinite()};
2123 if(anEquationCoeffs.mB > 0.0)
2125 if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) < -1.0)
2126 {//There is NOT intersection
2127 return Standard_True;
2129 else if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
2131 aU1f[0] = anEquationCoeffs.mFI1;
2132 aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
2134 else if((1 + anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
2135 (anEquationCoeffs.mB <= 1 - anEquationCoeffs.mC))
2137 Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
2143 const Standard_Real aDAngle = acos(anArg);
2144 //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
2145 aU1f[0] = anEquationCoeffs.mFI1;
2146 aU1l[0] = aDAngle + anEquationCoeffs.mFI1;
2147 aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
2148 aU1l[1] = aPeriod + anEquationCoeffs.mFI1;
2150 else if((1 - anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
2151 (anEquationCoeffs.mB <= 1 + anEquationCoeffs.mC))
2153 Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
2159 const Standard_Real aDAngle = acos(anArg);
2160 //U=[aDAngle;2*PI-aDAngle]+aFI1
2162 aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
2163 aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
2165 else if(anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
2167 Standard_Real anArg1 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB,
2168 anArg2 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
2179 const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
2180 //(U=[aDAngle1;aDAngle2]+aFI1) ||
2181 //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
2183 aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1;
2184 aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1;
2185 aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1;
2186 aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1;
2190 Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
2193 else if(anEquationCoeffs.mB < 0.0)
2195 if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) > 1.0)
2196 {//There is NOT intersection
2197 return Standard_True;
2199 else if(-anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
2201 aU1f[0] = anEquationCoeffs.mFI1;
2202 aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
2204 else if((-anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
2205 ( anEquationCoeffs.mB <= anEquationCoeffs.mC - 1))
2207 Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
2213 const Standard_Real aDAngle = acos(anArg);
2214 //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
2216 aU1f[0] = anEquationCoeffs.mFI1;
2217 aU1l[0] = aDAngle + anEquationCoeffs.mFI1;
2218 aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
2219 aU1l[1] = aPeriod + anEquationCoeffs.mFI1;
2221 else if((anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
2222 (anEquationCoeffs.mB <= -anEquationCoeffs.mB - 1))
2224 Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
2230 const Standard_Real aDAngle = acos(anArg);
2231 //U=[aDAngle;2*PI-aDAngle]+aFI1
2233 aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
2234 aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
2236 else if(-anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
2238 Standard_Real anArg1 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB,
2239 anArg2 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
2250 const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
2251 //(U=[aDAngle1;aDAngle2]+aFI1) ||
2252 //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
2254 aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1;
2255 aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1;
2256 aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1;
2257 aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1;
2261 Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
2266 Standard_Failure::Raise("Error. Exception. Unhandled case (B-parameter computation)!");
2269 for(Standard_Integer i = 0; i < aNbOfBoundaries; i++)
2271 if(Precision::IsInfinite(aU1f[i]) && Precision::IsInfinite(aU1l[i]))
2274 InscribeInterval(aUSurf1f, aUSurf1l, aU1f[i], aU1l[i], theTol2D, aPeriod);
2277 if( !Precision::IsInfinite(aU1f[0]) && !Precision::IsInfinite(aU1f[1]) &&
2278 !Precision::IsInfinite(aU1l[0]) && !Precision::IsInfinite(aU1l[1]))
2280 if( ((aU1f[1] <= aU1l[0]) || (aU1l[1] <= aU1l[0])) &&
2281 ((aU1f[0] <= aU1l[1]) || (aU1l[0] <= aU1l[1])))
2282 {//Join all intervals to one
2283 aU1f[0] = Min(aU1f[0], aU1f[1]);
2284 aU1l[0] = Max(aU1l[0], aU1l[1]);
2286 aU1f[1] = -Precision::Infinite();
2287 aU1l[1] = Precision::Infinite();
2292 //[0...1] - in these points parameter U1 goes through
2293 // the seam-edge of the first cylinder.
2294 //[2...3] - First and last U1 parameter.
2295 //[4...5] - in these points parameter U2 goes through
2296 // the seam-edge of the second cylinder.
2297 //[6...9] - in these points an intersection line goes through
2298 // U-boundaries of the second surface.
2299 const Standard_Integer aNbCritPointsMax = 10;
2300 Standard_Real anU1crit[aNbCritPointsMax] = {Precision::Infinite(),
2301 Precision::Infinite(),
2302 Precision::Infinite(),
2303 Precision::Infinite(),
2304 Precision::Infinite(),
2305 Precision::Infinite(),
2306 Precision::Infinite(),
2307 Precision::Infinite(),
2308 Precision::Infinite(),
2309 Precision::Infinite()};
2311 CriticalPointsComputing(anEquationCoeffs, aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
2312 aPeriod, theTol2D, aNbCritPointsMax, anU1crit);
2315 //Getting Walking-line
2319 WLFStatus_Absent = 0,
2320 WLFStatus_Exist = 1,
2321 WLFStatus_Broken = 2
2324 for(Standard_Integer aCurInterval = 0; aCurInterval < aNbOfBoundaries; aCurInterval++)
2326 if(Precision::IsInfinite(aU1f[aCurInterval]) && Precision::IsInfinite(aU1l[aCurInterval]))
2329 Standard_Boolean isAddedIntoWL[aNbWLines];
2330 for(Standard_Integer i = 0; i < aNbWLines; i++)
2331 isAddedIntoWL[i] = Standard_False;
2333 Standard_Real anUf = aU1f[aCurInterval], anUl = aU1l[aCurInterval];
2334 const Standard_Boolean isDeltaPeriod = IsEqual(anUl-anUf, aPeriod);
2336 //Inscribe and sort critical points
2337 InscribeAndSortArray(anU1crit, aNbCritPointsMax, anUf, anUl, theTol2D, aPeriod);
2341 Standard_Real aU2[aNbWLines], aV1[aNbWLines], aV2[aNbWLines];
2342 WLFStatus aWLFindStatus[aNbWLines];
2343 Standard_Real aV1Prev[aNbWLines], aV2Prev[aNbWLines];
2344 Standard_Real anUexpect[aNbWLines];
2345 Standard_Boolean isAddingWLEnabled[aNbWLines];
2347 Handle(IntSurf_LineOn2S) aL2S[aNbWLines];
2348 Handle(IntPatch_WLine) aWLine[aNbWLines];
2349 for(Standard_Integer i = 0; i < aNbWLines; i++)
2351 aL2S[i] = new IntSurf_LineOn2S();
2352 aWLine[i] = new IntPatch_WLine(aL2S[i], Standard_False);
2353 aWLFindStatus[i] = WLFStatus_Absent;
2354 isAddingWLEnabled[i] = Standard_True;
2355 aU2[i] = aV1[i] = aV2[i] = 0.0;
2356 aV1Prev[i] = aV2Prev[i] = 0.0;
2357 anUexpect[i] = anUf;
2360 Standard_Real anU1 = anUf;
2362 Standard_Real aCriticalDelta[aNbCritPointsMax];
2363 for(Standard_Integer i = 0; i < aNbCritPointsMax; i++)
2364 aCriticalDelta[i] = anU1 - anU1crit[i];
2366 Standard_Boolean isFirst = Standard_True;
2370 for(Standard_Integer i = 0; i < aNbCritPointsMax; i++)
2372 if((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0)
2376 for(Standard_Integer i = 0; i < aNbWLines; i++)
2378 aWLFindStatus[i] = WLFStatus_Broken;
2379 anUexpect[i] = anU1;
2386 if(IsEqual(anU1, anUl))
2388 for(Standard_Integer i = 0; i < aNbWLines; i++)
2390 aWLFindStatus[i] = WLFStatus_Broken;
2391 anUexpect[i] = anU1;
2395 //if isAddedIntoWL* == TRUE WLine contains only one point
2396 //(which was end point of previous WLine). If we will
2397 //add point found on the current step WLine will contain only
2398 //two points. At that both these points will be equal to the
2399 //points found earlier. Therefore, new WLine will repeat
2400 //already existing WLine. Consequently, it is necessary
2401 //to forbid building new line in this case.
2403 isAddingWLEnabled[i] = (!isAddedIntoWL[i]);
2407 isAddingWLEnabled[i] = ((theTol2D >= (anUexpect[i] - anU1)) ||
2408 (aWLFindStatus[i] == WLFStatus_Absent));
2410 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2414 for(Standard_Integer i = 0; i < aNbWLines; i++)
2416 isAddingWLEnabled[i] = ((theTol2D >= (anUexpect[i] - anU1)) ||
2417 (aWLFindStatus[i] == WLFStatus_Absent));
2418 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2421 for(Standard_Integer i = 0; i < aNbWLines; i++)
2423 const Standard_Integer aNbPntsWL = aWLine[i].IsNull() ? 0 :
2424 aWLine[i]->Curve()->NbPoints();
2426 CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]);
2428 InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False);
2431 {//the line has not contained any points yet
2432 if(((aUSurf2f + aPeriod - aUSurf2l) <= 2.0*theTol2D) &&
2433 ((Abs(aU2[i] - aUSurf2f) < theTol2D) ||
2434 (Abs(aU2[i]-aUSurf2l) < theTol2D)))
2436 //In this case aU2[i] can have two values: current aU2[i] or
2437 //aU2[i]+aPeriod (aU2[i]-aPeriod). It is necessary to choose
2440 Standard_Boolean isIncreasing = Standard_True;
2441 CylCylMonotonicity(anU1, i, anEquationCoeffs, aPeriod, isIncreasing);
2443 //If U2(U1) is increasing and U2 is considered to be equal aUSurf2l
2444 //then after the next step (when U1 will be increased) U2 will be
2445 //increased too. And we will go out of surface boundary.
2446 //Therefore, If U2(U1) is increasing then U2 must be equal aUSurf2f.
2447 //Analogically, if U2(U1) is decreasing.
2460 if(((aUSurf2l - aUSurf2f) >= aPeriod) &&
2461 ((Abs(aU2[i] - aUSurf2f) < theTol2D) ||
2462 (Abs(aU2[i]-aUSurf2l) < theTol2D)))
2464 Standard_Real aU2prev = 0.0, aV2prev = 0.0;
2466 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS1(aU2prev, aV2prev);
2468 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS2(aU2prev, aV2prev);
2470 if(2.0*Abs(aU2prev - aU2[i]) > aPeriod)
2472 if(aU2prev > aU2[i])
2480 if( (aWLFindStatus[i] == WLFStatus_Broken) ||
2481 (aWLFindStatus[i] == WLFStatus_Absent))
2482 {//Begin and end of WLine must be on boundary point
2483 //or on seam-edge strictly (if it is possible).
2484 if(Abs(aU2[i]) <= theTol2D)
2486 else if(Abs(aU2[i] - aPeriod) <= theTol2D)
2488 else if(Abs(aU2[i] - aUSurf2f) <= theTol2D)
2490 else if(Abs(aU2[i] - aUSurf2l) <= theTol2D)
2494 CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs, aV1[i], aV2[i]);
2498 aV1Prev[i] = aV1[i];
2499 aV2Prev[i] = aV2[i];
2501 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2503 isFirst = Standard_False;
2505 //Looking for points into WLine
2506 Standard_Boolean isBroken = Standard_False;
2507 for(Standard_Integer i = 0; i < aNbWLines; i++)
2509 if(!isAddingWLEnabled[i])
2511 Standard_Boolean isBoundIntersect = Standard_False;
2512 if( (Abs(aV1[i] - aVSurf1f) <= theTol2D) ||
2513 ((aV1[i]-aVSurf1f)*(aV1Prev[i]-aVSurf1f) < 0.0))
2515 isBoundIntersect = Standard_True;
2517 else if( (Abs(aV1[i] - aVSurf1l) <= theTol2D) ||
2518 ( (aV1[i]-aVSurf1l)*(aV1Prev[i]-aVSurf1l) < 0.0))
2520 isBoundIntersect = Standard_True;
2522 else if( (Abs(aV2[i] - aVSurf2f) <= theTol2D) ||
2523 ( (aV2[i]-aVSurf2f)*(aV2Prev[i]-aVSurf2f) < 0.0))
2525 isBoundIntersect = Standard_True;
2527 else if( (Abs(aV2[i] - aVSurf2l) <= theTol2D) ||
2528 ( (aV2[i]-aVSurf2l)*(aV2Prev[i]-aVSurf2l) < 0.0))
2530 isBoundIntersect = Standard_True;
2533 aV1Prev[i] = aV1[i];
2534 aV2Prev[i] = aV2[i];
2536 if(aWLFindStatus[i] == WLFStatus_Broken)
2537 isBroken = Standard_True;
2539 if(!isBoundIntersect)
2545 anUexpect[i] = anU1;
2549 const Standard_Boolean isInscribe =
2550 ((aUSurf2f-aU2[i]) <= theTol2D) && ((aU2[i]-aUSurf2l) <= theTol2D) &&
2551 ((aVSurf1f - aV1[i]) <= theTol2D) && ((aV1[i] - aVSurf1l) <= theTol2D) &&
2552 ((aVSurf2f - aV2[i]) <= theTol2D) && ((aV2[i] - aVSurf2l) <= theTol2D);
2554 //isVIntersect == TRUE if intersection line intersects two (!)
2555 //V-bounds of cylinder (1st or 2nd - no matter)
2556 const Standard_Boolean isVIntersect =
2557 ( ((aVSurf1f-aV1[i])*(aVSurf1f-aV1Prev[i]) < RealSmall()) &&
2558 ((aVSurf1l-aV1[i])*(aVSurf1l-aV1Prev[i]) < RealSmall())) ||
2559 ( ((aVSurf2f-aV2[i])*(aVSurf2f-aV2Prev[i]) < RealSmall()) &&
2560 ((aVSurf2l-aV2[i])*(aVSurf2l-aV2Prev[i]) < RealSmall()));
2563 //isFound1 == TRUE if intersection line intersects V-bounds
2564 // (First or Last - no matter) of the 1st cylynder
2565 //isFound2 == TRUE if intersection line intersects V-bounds
2566 // (First or Last - no matter) of the 2nd cylynder
2567 Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
2568 Standard_Boolean isForce = Standard_False;
2570 if((aWLFindStatus[i] == WLFStatus_Absent))
2572 if(((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1-aUSurf1l) < theTol2D))
2574 isForce = Standard_True;
2578 AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs,
2579 theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod,
2580 anU1, aU2[i], aV1[i], aV1Prev[i],
2581 aV2[i], aV2Prev[i], i, isTheReverse,
2582 isForce, isFound1, isFound2);
2584 const Standard_Boolean isPrevVBound = !isVIntersect &&
2585 ((Abs(aV1Prev[i] - aVSurf1f) <= theTol2D) ||
2586 (Abs(aV1Prev[i] - aVSurf1l) <= theTol2D) ||
2587 (Abs(aV2Prev[i] - aVSurf2f) <= theTol2D) ||
2588 (Abs(aV2Prev[i] - aVSurf2l) <= theTol2D));
2590 if((aWLFindStatus[i] == WLFStatus_Exist) && (isFound1 || isFound2) && !isPrevVBound)
2592 aWLFindStatus[i] = WLFStatus_Broken; //start a new line
2596 if((aWLFindStatus[i] == WLFStatus_Absent) && (isFound1 || isFound2))
2598 aWLFindStatus[i] = WLFStatus_Exist;
2601 if(( aWLFindStatus[i] != WLFStatus_Broken) || (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl))
2603 if(aWLine[i]->NbPnts() > 0)
2605 Standard_Real aU2p = 0.0, aV2p = 0.0;
2607 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
2609 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
2611 const Standard_Real aDelta = aU2[i] - aU2p;
2613 if(2*Abs(aDelta) > aPeriod)
2626 if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse, Standard_True,
2627 gp_Pnt2d(anU1, aV1[i]), gp_Pnt2d(aU2[i], aV2[i]),
2628 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
2629 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod,
2630 aWLine[i]->Curve(), i, theTol3D, theTol2D, isForce))
2632 if(aWLFindStatus[i] == WLFStatus_Absent)
2634 aWLFindStatus[i] = WLFStatus_Exist;
2637 else if(!isFound1 && !isFound2)
2638 {//We do not add any point while doing this iteration
2639 if(aWLFindStatus[i] == WLFStatus_Exist)
2641 aWLFindStatus[i] = WLFStatus_Broken;
2647 {//We do not add any point while doing this iteration
2648 if(aWLFindStatus[i] == WLFStatus_Exist)
2650 aWLFindStatus[i] = WLFStatus_Broken;
2654 aV1Prev[i] = aV1[i];
2655 aV2Prev[i] = aV2[i];
2657 if(aWLFindStatus[i] == WLFStatus_Broken)
2658 isBroken = Standard_True;
2659 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2662 {//current lines are filled. Go to the next lines
2665 Standard_Boolean isAdded = Standard_True;
2667 for(Standard_Integer i = 0; i < aNbWLines; i++)
2669 if(isAddingWLEnabled[i])
2674 isAdded = Standard_False;
2676 Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
2678 AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs,
2679 theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod,
2680 anU1, aU2[i], aV1[i], aV1Prev[i],
2681 aV2[i], aV2Prev[i], i, isTheReverse,
2682 Standard_False, isFound1, isFound2);
2684 if(isFound1 || isFound2)
2686 isAdded = Standard_True;
2689 if(aWLine[i]->NbPnts() > 0)
2691 Standard_Real aU2p = 0.0, aV2p = 0.0;
2693 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
2695 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
2697 const Standard_Real aDelta = aU2[i] - aU2p;
2699 if(2*Abs(aDelta) > aPeriod)
2712 if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse,
2713 Standard_True, gp_Pnt2d(anU1, aV1[i]),
2714 gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
2715 aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
2716 aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
2717 i, theTol3D, theTol2D, Standard_False))
2719 isAdded = Standard_True;
2725 Standard_Real anUmaxAdded = RealFirst();
2726 for(Standard_Integer i = 0; i < aNbWLines; i++)
2728 Standard_Real aU1c = 0.0, aV1c = 0.0;
2730 aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS2(aU1c, aV1c);
2732 aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS1(aU1c, aV1c);
2734 anUmaxAdded = Max(anUmaxAdded, aU1c);
2737 for(Standard_Integer i = 0; i < aNbWLines; i++)
2739 if(isAddingWLEnabled[i])
2744 CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs, aU2[i], aV1[i], aV2[i]);
2746 AddPointIntoWL( theQuad1, theQuad2, anEquationCoeffs, isTheReverse,
2747 Standard_True, gp_Pnt2d(anUmaxAdded, aV1[i]),
2748 gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
2749 aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
2750 aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
2751 i, theTol3D, theTol2D, Standard_False);
2761 const Standard_Real aDeltaV1 = (aVSurf1l - aVSurf1f)/IntToReal(aNbPoints),
2762 aDeltaV2 = (aVSurf2l - aVSurf2f)/IntToReal(aNbPoints);
2764 math_Matrix aMatr(1, 3, 1, 5);
2766 Standard_Real aMinUexp = RealLast();
2768 for(Standard_Integer i = 0; i < aNbWLines; i++)
2770 if(theTol2D < (anUexpect[i] - anU1))
2775 if(aWLFindStatus[i] == WLFStatus_Absent)
2777 anUexpect[i] += aStepMax;
2778 aMinUexp = Min(aMinUexp, anUexpect[i]);
2782 Standard_Real aStepTmp = aStepMax;
2784 const Standard_Real aSinU1 = sin(anU1),
2786 aSinU2 = sin(aU2[i]),
2787 aCosU2 = cos(aU2[i]);
2789 aMatr.SetCol(1, anEquationCoeffs.mVecC1);
2790 aMatr.SetCol(2, anEquationCoeffs.mVecC2);
2791 aMatr.SetCol(3, anEquationCoeffs.mVecA1*aSinU1 - anEquationCoeffs.mVecB1*aCosU1);
2792 aMatr.SetCol(4, anEquationCoeffs.mVecA2*aSinU2 - anEquationCoeffs.mVecB2*aCosU2);
2793 aMatr.SetCol(5, anEquationCoeffs.mVecA1*aCosU1 + anEquationCoeffs.mVecB1*aSinU1 +
2794 anEquationCoeffs.mVecA2*aCosU2 + anEquationCoeffs.mVecB2*aSinU2 +
2795 anEquationCoeffs.mVecD);
2797 if(!StepComputing(aMatr, aV1[i], aV2[i], aDeltaV1, aDeltaV2, aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aStepTmp))
2799 //To avoid cycling-up
2800 anUexpect[i] += aStepMax;
2801 aMinUexp = Min(aMinUexp, anUexpect[i]);
2806 if(aStepTmp < aStepMin)
2807 aStepTmp = aStepMin;
2809 if(aStepTmp > aStepMax)
2810 aStepTmp = aStepMax;
2812 anUexpect[i] = anU1 + aStepTmp;
2813 aMinUexp = Min(aMinUexp, anUexpect[i]);
2819 if(Precision::PConfusion() >= (anUl - anU1))
2824 for(Standard_Integer i = 0; i < aNbWLines; i++)
2826 if(aWLine[i]->NbPnts() != 1)
2827 isAddedIntoWL[i] = Standard_False;
2830 {//strictly equal. Tolerance is considered above.
2831 anUexpect[i] = anUl;
2836 for(Standard_Integer i = 0; i < aNbWLines; i++)
2838 if((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
2840 isTheEmpty = Standard_False;
2841 Standard_Real u1, v1, u2, v2;
2842 aWLine[i]->Point(1).Parameters(u1, v1, u2, v2);
2844 aP.SetParameter(u1);
2845 aP.SetParameters(u1, v1, u2, v2);
2846 aP.SetTolerance(theTol3D);
2847 aP.SetValue(aWLine[i]->Point(1).Value());
2851 else if(aWLine[i]->NbPnts() > 1)
2853 Standard_Boolean isGood = Standard_True;
2855 if(aWLine[i]->NbPnts() == 2)
2857 const IntSurf_PntOn2S& aPf = aWLine[i]->Point(1);
2858 const IntSurf_PntOn2S& aPl = aWLine[i]->Point(2);
2860 if(aPf.IsSame(aPl, Precision::Confusion()))
2861 isGood = Standard_False;
2866 isTheEmpty = Standard_False;
2867 isAddedIntoWL[i] = Standard_True;
2868 SeekAdditionalPoints( theQuad1, theQuad2, aWLine[i]->Curve(),
2869 anEquationCoeffs, i, aNbPoints, 1,
2870 aWLine[i]->NbPnts(), aUSurf2f, aUSurf2l,
2871 theTol2D, aPeriod, isTheReverse);
2873 aWLine[i]->ComputeVertexParameters(theTol3D);
2874 theSlin.Append(aWLine[i]);
2879 isAddedIntoWL[i] = Standard_False;
2883 //aWLine[i]->Dump();
2890 //Delete the points in theSPnt, which
2891 //lie at least in one of the line in theSlin.
2892 for(Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
2894 for(Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); aNbLin++)
2896 const Handle(IntPatch_WLine)& aWLine1 =
2897 Handle(IntPatch_WLine)::DownCast(theSlin.Value(aNbLin));
2899 const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1);
2900 const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aWLine1->NbPnts());
2902 const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNbPnt).PntOn2S();
2903 if( aPntCur.IsSame(aPntFWL1, Precision::Confusion()) ||
2904 aPntCur.IsSame(aPntLWL1, Precision::Confusion()))
2906 theSPnt.Remove(aNbPnt);
2913 const Standard_Real aDU = aStepMin + Epsilon(aStepMin);
2914 //Try to add new points in the neighbourhood of existing point
2915 for(Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
2917 const IntPatch_Point& aPnt2S = theSPnt.Value(aNbPnt);
2919 Standard_Real u1, v1, u2, v2;
2920 aPnt2S.Parameters(u1, v1, u2, v2);
2922 Handle(IntSurf_LineOn2S) aL2S = new IntSurf_LineOn2S();
2923 Handle(IntPatch_WLine) aWLine = new IntPatch_WLine(aL2S, Standard_False);
2924 aWLine->Curve()->Add(aPnt2S.PntOn2S());
2926 //Define the index of WLine, which lies the point aPnt2S in.
2927 Standard_Real anUf = 0.0, anUl = 0.0, aCurU2 = 0.0;
2928 Standard_Integer anIndex = 0;
2931 anUf = Max(u2 - aStepMax, aUSurf1f);
2937 anUf = Max(u1 - aStepMax, aUSurf1f);
2941 Standard_Real aDelta = RealLast();
2942 for (Standard_Integer i = 0; i < aNbWLines; i++)
2944 Standard_Real anU2t = 0.0;
2945 if(!CylCylComputeParameters(anUl, i, anEquationCoeffs, anU2t))
2948 const Standard_Real aDU = Abs(anU2t - aCurU2);
2956 //Try to fill aWLine by additional points
2957 while(anUl - anUf > RealSmall())
2959 Standard_Real anU2 = 0.0, anV1 = 0.0, anV2 = 0.0;
2960 Standard_Boolean isDone =
2961 CylCylComputeParameters(anUf, anIndex, anEquationCoeffs,
2970 if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse, Standard_True,
2971 gp_Pnt2d(anUf, anV1), gp_Pnt2d(anU2, anV2),
2972 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
2973 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l,
2974 aPeriod, aWLine->Curve(), anIndex, theTol3D,
2975 theTol2D, Standard_False, Standard_True))
2981 if(aWLine->NbPnts() > 1)
2983 SeekAdditionalPoints( theQuad1, theQuad2, aWLine->Curve(),
2984 anEquationCoeffs, anIndex, aNbMinPoints,
2985 1, aWLine->NbPnts(), aUSurf2f, aUSurf2l,
2986 theTol2D, aPeriod, isTheReverse);
2988 aWLine->ComputeVertexParameters(theTol3D);
2989 theSlin.Append(aWLine);
2991 theSPnt.Remove(aNbPnt);
2996 return Standard_True;
2999 //=======================================================================
3000 //function : IntCySp
3002 //=======================================================================
3003 Standard_Boolean IntCySp(const IntSurf_Quadric& Quad1,
3004 const IntSurf_Quadric& Quad2,
3005 const Standard_Real Tol,
3006 const Standard_Boolean Reversed,
3007 Standard_Boolean& Empty,
3008 Standard_Boolean& Multpoint,
3009 IntPatch_SequenceOfLine& slin,
3010 IntPatch_SequenceOfPoint& spnt)
3015 IntSurf_TypeTrans trans1,trans2;
3016 IntAna_ResultType typint;
3017 IntPatch_Point ptsol;
3024 Cy = Quad1.Cylinder();
3025 Sp = Quad2.Sphere();
3028 Cy = Quad2.Cylinder();
3029 Sp = Quad1.Sphere();
3031 IntAna_QuadQuadGeo inter(Cy,Sp,Tol);
3033 if (!inter.IsDone()) {return Standard_False;}
3035 typint = inter.TypeInter();
3036 Standard_Integer NbSol = inter.NbSolutions();
3037 Empty = Standard_False;
3043 Empty = Standard_True;
3049 gp_Pnt psol(inter.Point(1));
3050 Standard_Real U1,V1,U2,V2;
3051 Quad1.Parameters(psol,U1,V1);
3052 Quad2.Parameters(psol,U2,V2);
3053 ptsol.SetValue(psol,Tol,Standard_True);
3054 ptsol.SetParameters(U1,V1,U2,V2);
3061 cirsol = inter.Circle(1);
3064 ElCLib::D1(0.,cirsol,ptref,Tgt);
3067 gp_Vec TestCurvature(ptref,Sp.Location());
3068 gp_Vec Normsp,Normcyl;
3070 Normcyl = Quad1.Normale(ptref);
3071 Normsp = Quad2.Normale(ptref);
3074 Normcyl = Quad2.Normale(ptref);
3075 Normsp = Quad1.Normale(ptref);
3078 IntSurf_Situation situcyl;
3079 IntSurf_Situation situsp;
3081 if (Normcyl.Dot(TestCurvature) > 0.) {
3082 situsp = IntSurf_Outside;
3083 if (Normsp.Dot(Normcyl) > 0.) {
3084 situcyl = IntSurf_Inside;
3087 situcyl = IntSurf_Outside;
3091 situsp = IntSurf_Inside;
3092 if (Normsp.Dot(Normcyl) > 0.) {
3093 situcyl = IntSurf_Outside;
3096 situcyl = IntSurf_Inside;
3099 Handle(IntPatch_GLine) glig;
3101 glig = new IntPatch_GLine(cirsol, Standard_True, situcyl, situsp);
3104 glig = new IntPatch_GLine(cirsol, Standard_True, situsp, situcyl);
3109 if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
3110 trans1 = IntSurf_Out;
3111 trans2 = IntSurf_In;
3114 trans1 = IntSurf_In;
3115 trans2 = IntSurf_Out;
3117 Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3120 cirsol = inter.Circle(2);
3121 ElCLib::D1(0.,cirsol,ptref,Tgt);
3122 Standard_Real qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
3123 if(qwe> 0.0000001) {
3124 trans1 = IntSurf_Out;
3125 trans2 = IntSurf_In;
3127 else if(qwe<-0.0000001) {
3128 trans1 = IntSurf_In;
3129 trans2 = IntSurf_Out;
3132 trans1=trans2=IntSurf_Undecided;
3134 glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3140 case IntAna_NoGeometricSolution:
3143 Standard_Real U1,V1,U2,V2;
3144 IntAna_IntQuadQuad anaint(Cy,Sp,Tol);
3145 if (!anaint.IsDone()) {
3146 return Standard_False;
3149 if (anaint.NbPnt()==0 && anaint.NbCurve()==0) {
3150 Empty = Standard_True;
3154 NbSol = anaint.NbPnt();
3155 for (i = 1; i <= NbSol; i++) {
3156 psol = anaint.Point(i);
3157 Quad1.Parameters(psol,U1,V1);
3158 Quad2.Parameters(psol,U2,V2);
3159 ptsol.SetValue(psol,Tol,Standard_True);
3160 ptsol.SetParameters(U1,V1,U2,V2);
3164 gp_Pnt ptvalid,ptf,ptl;
3166 Standard_Real first,last,para;
3167 IntAna_Curve curvsol;
3168 Standard_Boolean tgfound;
3169 Standard_Integer kount;
3171 NbSol = anaint.NbCurve();
3172 for (i = 1; i <= NbSol; i++) {
3173 curvsol = anaint.Curve(i);
3174 curvsol.Domain(first,last);
3175 ptf = curvsol.Value(first);
3176 ptl = curvsol.Value(last);
3180 tgfound = Standard_False;
3183 para = (1.123*first + para)/2.123;
3184 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
3187 tgfound = kount > 5;
3190 Handle(IntPatch_ALine) alig;
3192 Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
3193 Quad1.Normale(ptvalid));
3194 if(qwe> 0.00000001) {
3195 trans1 = IntSurf_Out;
3196 trans2 = IntSurf_In;
3198 else if(qwe<-0.00000001) {
3199 trans1 = IntSurf_In;
3200 trans2 = IntSurf_Out;
3203 trans1=trans2=IntSurf_Undecided;
3205 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
3208 alig = new IntPatch_ALine(curvsol,Standard_False);
3210 Standard_Boolean TempFalse1a = Standard_False;
3211 Standard_Boolean TempFalse2a = Standard_False;
3213 //-- ptf et ptl : points debut et fin de alig
3215 ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1a,ptf,first,
3216 TempFalse2a,ptl,last,Multpoint,Tol);
3218 } //-- boucle sur les lignes
3219 } //-- solution avec au moins une lihne
3225 return Standard_False;
3228 return Standard_True;
3230 //=======================================================================
3231 //function : IntCyCo
3233 //=======================================================================
3234 Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
3235 const IntSurf_Quadric& Quad2,
3236 const Standard_Real Tol,
3237 const Standard_Boolean Reversed,
3238 Standard_Boolean& Empty,
3239 Standard_Boolean& Multpoint,
3240 IntPatch_SequenceOfLine& slin,
3241 IntPatch_SequenceOfPoint& spnt)
3244 IntPatch_Point ptsol;
3248 IntSurf_TypeTrans trans1,trans2;
3249 IntAna_ResultType typint;
3256 Cy = Quad1.Cylinder();
3260 Cy = Quad2.Cylinder();
3263 IntAna_QuadQuadGeo inter(Cy,Co,Tol);
3265 if (!inter.IsDone()) {return Standard_False;}
3267 typint = inter.TypeInter();
3268 Standard_Integer NbSol = inter.NbSolutions();
3269 Empty = Standard_False;
3273 case IntAna_Empty : {
3274 Empty = Standard_True;
3278 case IntAna_Point :{
3279 gp_Pnt psol(inter.Point(1));
3280 Standard_Real U1,V1,U2,V2;
3281 Quad1.Parameters(psol,U1,V1);
3282 Quad1.Parameters(psol,U2,V2);
3283 ptsol.SetValue(psol,Tol,Standard_True);
3284 ptsol.SetParameters(U1,V1,U2,V2);
3289 case IntAna_Circle: {
3295 for(j=1; j<=2; ++j) {
3296 cirsol = inter.Circle(j);
3297 ElCLib::D1(0.,cirsol,ptref,Tgt);
3298 qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
3299 if(qwe> 0.00000001) {
3300 trans1 = IntSurf_Out;
3301 trans2 = IntSurf_In;
3303 else if(qwe<-0.00000001) {
3304 trans1 = IntSurf_In;
3305 trans2 = IntSurf_Out;
3308 trans1=trans2=IntSurf_Undecided;
3310 Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3316 case IntAna_NoGeometricSolution: {
3318 Standard_Real U1,V1,U2,V2;
3319 IntAna_IntQuadQuad anaint(Cy,Co,Tol);
3320 if (!anaint.IsDone()) {
3321 return Standard_False;
3324 if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) {
3325 Empty = Standard_True;
3328 NbSol = anaint.NbPnt();
3329 for (i = 1; i <= NbSol; i++) {
3330 psol = anaint.Point(i);
3331 Quad1.Parameters(psol,U1,V1);
3332 Quad2.Parameters(psol,U2,V2);
3333 ptsol.SetValue(psol,Tol,Standard_True);
3334 ptsol.SetParameters(U1,V1,U2,V2);
3338 gp_Pnt ptvalid, ptf, ptl;
3341 Standard_Real first,last,para;
3342 Standard_Boolean tgfound,firstp,lastp,kept;
3343 Standard_Integer kount;
3346 //IntAna_Curve curvsol;
3348 IntAna_ListOfCurve aLC;
3349 IntAna_ListIteratorOfListOfCurve aIt;
3352 NbSol = anaint.NbCurve();
3353 for (i = 1; i <= NbSol; ++i) {
3354 kept = Standard_False;
3355 //curvsol = anaint.Curve(i);
3358 ExploreCurve(Cy, Co, aC, 10.*Tol, aLC);
3360 aIt.Initialize(aLC);
3361 for (; aIt.More(); aIt.Next()) {
3362 IntAna_Curve& curvsol=aIt.Value();
3364 curvsol.Domain(first, last);
3365 firstp = !curvsol.IsFirstOpen();
3366 lastp = !curvsol.IsLastOpen();
3368 ptf = curvsol.Value(first);
3371 ptl = curvsol.Value(last);
3375 tgfound = Standard_False;
3378 para = (1.123*first + para)/2.123;
3379 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
3382 tgfound = kount > 5;
3385 Handle(IntPatch_ALine) alig;
3387 Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
3388 Quad1.Normale(ptvalid));
3389 if(qwe> 0.00000001) {
3390 trans1 = IntSurf_Out;
3391 trans2 = IntSurf_In;
3393 else if(qwe<-0.00000001) {
3394 trans1 = IntSurf_In;
3395 trans2 = IntSurf_Out;
3398 trans1=trans2=IntSurf_Undecided;
3400 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
3401 kept = Standard_True;
3404 ptvalid = curvsol.Value(para);
3405 alig = new IntPatch_ALine(curvsol,Standard_False);
3406 kept = Standard_True;
3407 //-- cout << "Transition indeterminee" << endl;
3410 Standard_Boolean Nfirstp = !firstp;
3411 Standard_Boolean Nlastp = !lastp;
3412 ProcessBounds(alig,slin,Quad1,Quad2,Nfirstp,ptf,first,
3413 Nlastp,ptl,last,Multpoint,Tol);
3416 } // for (; aIt.More(); aIt.Next())
3417 } // for (i = 1; i <= NbSol; ++i)
3423 return Standard_False;
3425 } // switch (typint)
3427 return Standard_True;
3429 //=======================================================================
3430 //function : ExploreCurve
3432 //=======================================================================
3433 Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy,
3436 const Standard_Real aTol,
3437 IntAna_ListOfCurve& aLC)
3440 Standard_Boolean bFind=Standard_False;
3441 Standard_Real aTheta, aT1, aT2, aDst;
3451 aC.Domain(aT1, aT2);
3454 aDst=aPx.Distance(aPapx);
3459 aDst=aPx.Distance(aPapx);
3464 bFind=aC.FindParameter(aPapx, aTheta);
3469 aPx=aC.Value(aTheta);
3470 aDst=aPx.Distance(aPapx);
3475 // need to be splitted at aTheta
3476 IntAna_Curve aC1, aC2;
3479 aC1.SetDomain(aT1, aTheta);
3481 aC2.SetDomain(aTheta, aT2);
3489 //=======================================================================
3490 //function : IsToReverse
3492 //=======================================================================
3493 Standard_Boolean IsToReverse(const gp_Cylinder& Cy1,
3494 const gp_Cylinder& Cy2,
3495 const Standard_Real Tol)
3497 Standard_Boolean bRet;
3498 Standard_Real aR1,aR2, dR, aSc1, aSc2;
3500 bRet=Standard_False;
3512 gp_Dir aDZ(0.,0.,1.);
3514 const gp_Dir& aDir1=Cy1.Axis().Direction();
3520 const gp_Dir& aDir2=Cy2.Axis().Direction();