projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0026193: Incomplete intersection curve
[occt.git] / src / IntPatch / IntPatch_ImpImpIntersection_4.gxx
CommitLineData
b311480e 1// Created on: 1992-05-07
2// Created by: Jacques GOUSSARD
3// Copyright (c) 1992-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 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
973c2be1 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
17#include <IntAna_ListOfCurve.hxx>
18#include <IntAna_ListIteratorOfListOfCurve.hxx>
ecc4f148 19#include <IntPatch_WLine.hxx>
20
b5ef9d91 21#include <math_Matrix.hxx>
22
23//If Abs(a) <= aNulValue then it is considered that a = 0.
24static const Standard_Real aNulValue = 1.0e-11;
25
26struct stCoeffsValue;
fa0291ff 27//
7fd59977 28static
29 Standard_Boolean ExploreCurve(const gp_Cylinder& aCy,
30 const gp_Cone& aCo,
31 IntAna_Curve& aC,
32 const Standard_Real aTol,
33 IntAna_ListOfCurve& aLC);
34static
35 Standard_Boolean IsToReverse(const gp_Cylinder& Cy1,
36 const gp_Cylinder& Cy2,
37 const Standard_Real Tol);
38
b5ef9d91 39static 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);
45
46static 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);
59
60//=======================================================================
61//function : MinMax
62//purpose : Replaces theParMIN = MIN(theParMIN, theParMAX),
9e20ed57 63// theParMAX = MAX(theParMIN, theParMAX).
b5ef9d91 64//=======================================================================
9e20ed57 65static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX)
66{
67 if(theParMIN > theParMAX)
68 {
69 const Standard_Real aux = theParMAX;
70 theParMAX = theParMIN;
71 theParMIN = aux;
72 }
73}
74
b5ef9d91 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//=======================================================================
82static 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)
89{
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);
94
95 aSyst.SetCol(1, theMatr.Col(1));
96 aSyst.SetCol(2, theMatr.Col(2));
97 aSyst.SetCol(3, theMatr.Col(4));
98
99 const Standard_Real aDet = aSyst.Determinant();
100
101 aSyst.SetCol(1, theMatr.Col(3));
102 const Standard_Real aDet1 = aSyst.Determinant();
103
104 aSyst.SetCol(1, theMatr.Col(1));
105 aSyst.SetCol(2, theMatr.Col(3));
106
107 const Standard_Real aDet2 = aSyst.Determinant();
108
109 if(aDet*aDet1 > 0.0)
110 {
111 theV1Set = Max(theV1AfterDecrByDelta, theV1f);
112 }
113
114 if(aDet*aDet2 > 0.0)
115 {
116 theV2Set = Max(theV2AfterDecrByDelta, theV2f);
117 }
118}
119
120//=======================================================================
121//function : DeltaU1Computing
122//purpose : Computes new step for U1 parameter.
123//=======================================================================
124static inline
125 Standard_Boolean DeltaU1Computing(const math_Matrix& theSyst,
126 const math_Vector& theFree,
127 Standard_Real& theDeltaU1Found)
128{
129 Standard_Real aDet = theSyst.Determinant();
130
131 if(Abs(aDet) > aNulValue)
132 {
133 math_Matrix aSyst1(theSyst);
134 aSyst1.SetCol(2, theFree);
135
136 theDeltaU1Found = Abs(aSyst1.Determinant()/aDet);
137 return Standard_True;
138 }
139
140 return Standard_False;
141}
142
143//=======================================================================
144//function : StepComputing
145//purpose :
146//
147//Attention!!!:
148// theMatr must have 3*5-dimension strictly.
149// For system
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//=======================================================================
158static 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*/)
171{
172 Standard_Boolean isSuccess = Standard_False;
173 theDeltaU1Found/* = theDeltaU2Found*/ = RealLast();
174 //theV1Found = theV1set;
175 //theV2Found = theV2Set;
176 const Standard_Integer aNbDim = 3;
177
178 math_Matrix aSyst(1, aNbDim, 1, aNbDim);
179 math_Vector aFree(1, aNbDim);
180
181 //By default, increasing V1(U1) and V2(U1) functions is
182 //considered
183 Standard_Real aV1Set = Min(theV1Cur + theDeltaV1, theV1l),
184 aV2Set = Min(theV2Cur + theDeltaV2, theV2l);
185
186 //However, what is indead?
187 VBoundaryPrecise( theMatr, theV1Cur - theDeltaV1, theV2Cur - theDeltaV2,
188 theV1f, theV2f, aV1Set, aV2Set);
189
190 aSyst.SetCol(2, theMatr.Col(3));
191 aSyst.SetCol(3, theMatr.Col(4));
192
193 for(Standard_Integer i = 0; i < 2; i++)
194 {
195 if(i == 0)
196 {//V1 is known
197 aSyst.SetCol(1, theMatr.Col(2));
198 aFree.Set(1, aNbDim, theMatr.Col(5)-aV1Set*theMatr.Col(1));
199 }
200 else
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));
204 }
205
206 Standard_Real aNewDU = theDeltaU1Found;
207 if(DeltaU1Computing(aSyst, aFree, aNewDU))
208 {
209 isSuccess = Standard_True;
210 if(aNewDU < theDeltaU1Found)
211 {
212 theDeltaU1Found = aNewDU;
213 }
214 }
215 }
216
217 if(!isSuccess)
218 {
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));
223
224 //Now we have overdetermined system.
225
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);
232
233 if(anAbsD1 >= anAbsD2)
234 {
235 if(anAbsD1 >= anAbsD3)
236 {
237 //Det1
238 if(anAbsD1 <= aNulValue)
239 return isSuccess;
240
241 theDeltaU1Found = Abs(aFree(1)*theMatr(2,4) - aFree(2)*theMatr(1,4))/anAbsD1;
242 isSuccess = Standard_True;
243 }
244 else
245 {
246 //Det3
247 if(anAbsD3 <= aNulValue)
248 return isSuccess;
249
250 theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
251 isSuccess = Standard_True;
252 }
253 }
254 else
255 {
256 if(anAbsD2 >= anAbsD3)
257 {
258 //Det2
259 if(anAbsD2 <= aNulValue)
260 return isSuccess;
261
262 theDeltaU1Found = Abs(aFree(1)*theMatr(3,4) - aFree(3)*theMatr(1,4))/anAbsD2;
263 isSuccess = Standard_True;
264 }
265 else
266 {
267 //Det3
268 if(anAbsD3 <= aNulValue)
269 return isSuccess;
270
271 theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
272 isSuccess = Standard_True;
273 }
274 }
275 }
276
277 return isSuccess;
278}
9e20ed57 279
7fd59977 280//=======================================================================
281//function : ProcessBounds
282//purpose :
283//=======================================================================
284void 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)
296{
297 Standard_Integer j,k;
298 Standard_Real U1,V1,U2,V2;
299 IntPatch_Point ptsol;
300 Standard_Real d;
301
302 if (procf && procl) {
303 j = slin.Length() + 1;
304 }
305 else {
306 j = 1;
307 }
308
309
310 //-- On prend les lignes deja enregistrees
311
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));
315 k = 1;
316
317 //-- On prend les vertex des lignes deja enregistrees
318
319 while (k <= aligold->NbVertex()) {
320 ptsol = aligold->Vertex(k);
321 if (!procf) {
322 d=ptf.Distance(ptsol.Value());
323 if (d <= Tol) {
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);
329 }
330 ptsol.SetParameter(first);
331 alig->AddVertex(ptsol);
332 alig->SetFirstPoint(alig->NbVertex());
333 procf = Standard_True;
334
335 //-- On restore le point avec son parametre sur aligold
336 ptsol = aligold->Vertex(k);
337
338 }
339 }
340 if (!procl) {
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);
346 }
347 ptsol.SetParameter(last);
348 alig->AddVertex(ptsol);
349 alig->SetLastPoint(alig->NbVertex());
350 procl = Standard_True;
351
352 //-- On restore le point avec son parametre sur aligold
353 ptsol = aligold->Vertex(k);
354
355 }
356 }
357 if (procf && procl) {
358 k = aligold->NbVertex()+1;
359 }
360 else {
361 k = k+1;
362 }
363 }
364 if (procf && procl) {
365 j = slin.Length()+1;
366 }
367 else {
368 j = j+1;
369 }
370 }
371 }
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());
383
384 ptsol.SetParameter(last);
385 alig->AddVertex(ptsol);
386 alig->SetLastPoint(alig->NbVertex());
387 }
388 else {
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());
398 }
399 }
400 else if (!procf) {
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());
408 }
409 else if (!procl) {
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());
417 }
418}
419//=======================================================================
420//function : IntCyCy
421//purpose :
422//=======================================================================
423Standard_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)
431
432{
433 IntPatch_Point ptsol;
434
435 Standard_Integer i;
436
437 IntSurf_TypeTrans trans1,trans2;
438 IntAna_ResultType typint;
439
440 gp_Elips elipsol;
441 gp_Lin linsol;
442
443 gp_Cylinder Cy1(Quad1.Cylinder());
444 gp_Cylinder Cy2(Quad2.Cylinder());
445
446 IntAna_QuadQuadGeo inter(Cy1,Cy2,Tol);
447
ecc4f148 448 if (!inter.IsDone())
449 {
450 return Standard_False;
451 }
7fd59977 452
453 typint = inter.TypeInter();
454 Standard_Integer NbSol = inter.NbSolutions();
455 Empty = Standard_False;
456 Same = Standard_False;
457
ecc4f148 458 switch (typint)
459 {
460 case IntAna_Empty:
7fd59977 461 {
462 Empty = Standard_True;
463 }
464 break;
465
466 case IntAna_Same:
467 {
468 Same = Standard_True;
469 }
470 break;
ecc4f148 471
472 case IntAna_Point:
7fd59977 473 {
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);
480 spnt.Append(ptsol);
481 }
482 break;
483
484 case IntAna_Line:
485 {
486 gp_Pnt ptref;
ecc4f148 487 if (NbSol == 1)
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;
497
498 if (crb1.Dot(crb2) < 0.)
499 { // centre de courbures "opposes"
500 if (norm1.Dot(crb1) > 0.)
501 {
502 situcyl2 = IntSurf_Inside;
503 }
504 else
505 {
506 situcyl2 = IntSurf_Outside;
507 }
508
509 if (norm2.Dot(crb2) > 0.)
510 {
511 situcyl1 = IntSurf_Inside;
512 }
513 else
514 {
515 situcyl1 = IntSurf_Outside;
516 }
517 }
518 else
519 {
520 if (Cy1.Radius() < Cy2.Radius())
521 {
522 if (norm1.Dot(crb1) > 0.)
523 {
524 situcyl2 = IntSurf_Inside;
525 }
526 else
527 {
528 situcyl2 = IntSurf_Outside;
529 }
530
531 if (norm2.Dot(crb2) > 0.)
532 {
533 situcyl1 = IntSurf_Outside;
534 }
535 else
536 {
537 situcyl1 = IntSurf_Inside;
538 }
539 }
540 else
541 {
542 if (norm1.Dot(crb1) > 0.)
543 {
544 situcyl2 = IntSurf_Outside;
545 }
546 else
547 {
548 situcyl2 = IntSurf_Inside;
549 }
550
551 if (norm2.Dot(crb2) > 0.)
552 {
553 situcyl1 = IntSurf_Inside;
554 }
555 else
556 {
557 situcyl1 = IntSurf_Outside;
558 }
559 }
560 }
561
562 Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_True, situcyl1, situcyl2);
563 slin.Append(glig);
7fd59977 564 }
ecc4f148 565 else
566 {
567 for (i=1; i <= NbSol; i++)
568 {
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));
573 if (qwe >0.00000001)
574 {
575 trans1 = IntSurf_Out;
576 trans2 = IntSurf_In;
577 }
578 else if (qwe <-0.00000001)
579 {
580 trans1 = IntSurf_In;
581 trans2 = IntSurf_Out;
582 }
583 else
584 {
585 trans1=trans2=IntSurf_Undecided;
586 }
587
588 Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
589 slin.Append(glig);
590 }
7fd59977 591 }
592 }
593 break;
594
595 case IntAna_Ellipse:
596 {
597 gp_Vec Tgt;
598 gp_Pnt ptref;
fa0291ff 599 IntPatch_Point pmult1, pmult2;
600
7fd59977 601 elipsol = inter.Ellipse(1);
fa0291ff 602
603 gp_Pnt pttang1(ElCLib::Value(0.5*M_PI, elipsol));
604 gp_Pnt pttang2(ElCLib::Value(1.5*M_PI, elipsol));
7fd59977 605
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);
611
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);
616
617 Quad1.Parameters(pttang2,oU1,oV1);
618 Quad2.Parameters(pttang2,oU2,oV2);
619 pmult2.SetParameters(oU1,oV1,oU2,oV2);
fa0291ff 620
7fd59977 621 // on traite la premiere ellipse
622
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));
ecc4f148 626 if (qwe>0.00000001)
627 {
628 trans1 = IntSurf_Out;
629 trans2 = IntSurf_In;
7fd59977 630 }
ecc4f148 631 else if (qwe<-0.00000001)
632 {
633 trans1 = IntSurf_In;
634 trans2 = IntSurf_Out;
7fd59977 635 }
ecc4f148 636 else
637 {
638 trans1=trans2=IntSurf_Undecided;
7fd59977 639 }
ecc4f148 640
7fd59977 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);
fa0291ff 644 //
645 {
ecc4f148 646 Standard_Real aU1, aV1, aU2, aV2;
647 IntPatch_Point aIP;
648 gp_Pnt aP (ElCLib::Value(0., elipsol));
649 //
650 aIP.SetValue(aP,Tol,Standard_False);
651 aIP.SetMultiple(Standard_False);
652 //
653 Quad1.Parameters(aP, aU1, aV1);
654 Quad2.Parameters(aP, aU2, aV2);
655 aIP.SetParameters(aU1, aV1, aU2, aV2);
656 //
657 aIP.SetParameter(0.);
658 glig->AddVertex(aIP);
659 glig->SetFirstPoint(1);
660 //
661 aIP.SetParameter(2.*M_PI);
662 glig->AddVertex(aIP);
663 glig->SetLastPoint(2);
fa0291ff 664 }
665 //
666 pmult1.SetParameter(0.5*M_PI);
7fd59977 667 glig->AddVertex(pmult1);
fa0291ff 668 //
c6541a0c 669 pmult2.SetParameter(1.5*M_PI);
7fd59977 670 glig->AddVertex(pmult2);
fa0291ff 671
672 //
7fd59977 673 slin.Append(glig);
674
675 //-- Transitions calculee au point 0 OK
fa0291ff 676 //
7fd59977 677 // on traite la deuxieme ellipse
7fd59977 678 elipsol = inter.Ellipse(2);
679
680 Standard_Real param1 = ElCLib::Parameter(elipsol,pttang1);
681 Standard_Real param2 = ElCLib::Parameter(elipsol,pttang2);
ecc4f148 682 Standard_Real parampourtransition = 0.0;
683 if (param1 < param2)
684 {
685 pmult1.SetParameter(0.5*M_PI);
686 pmult2.SetParameter(1.5*M_PI);
687 parampourtransition = M_PI;
7fd59977 688 }
689 else {
ecc4f148 690 pmult1.SetParameter(1.5*M_PI);
691 pmult2.SetParameter(0.5*M_PI);
692 parampourtransition = 0.0;
7fd59977 693 }
694
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));
ecc4f148 698 if(qwe> 0.00000001)
699 {
700 trans1 = IntSurf_Out;
701 trans2 = IntSurf_In;
7fd59977 702 }
ecc4f148 703 else if(qwe< -0.00000001)
704 {
705 trans1 = IntSurf_In;
706 trans2 = IntSurf_Out;
7fd59977 707 }
ecc4f148 708 else
709 {
710 trans1=trans2=IntSurf_Undecided;
7fd59977 711 }
ecc4f148 712
7fd59977 713 //-- La transition a ete calculee sur un point de cette ligne
714 glig = new IntPatch_GLine(elipsol,Standard_False,trans1,trans2);
fa0291ff 715 //
716 {
ecc4f148 717 Standard_Real aU1, aV1, aU2, aV2;
718 IntPatch_Point aIP;
719 gp_Pnt aP (ElCLib::Value(0., elipsol));
720 //
721 aIP.SetValue(aP,Tol,Standard_False);
722 aIP.SetMultiple(Standard_False);
723 //
724 Quad1.Parameters(aP, aU1, aV1);
725 Quad2.Parameters(aP, aU2, aV2);
726 aIP.SetParameters(aU1, aV1, aU2, aV2);
727 //
728 aIP.SetParameter(0.);
729 glig->AddVertex(aIP);
730 glig->SetFirstPoint(1);
731 //
732 aIP.SetParameter(2.*M_PI);
733 glig->AddVertex(aIP);
734 glig->SetLastPoint(2);
7fd59977 735 }
fa0291ff 736 //
7fd59977 737 glig->AddVertex(pmult1);
fa0291ff 738 glig->AddVertex(pmult2);
739 //
7fd59977 740 slin.Append(glig);
741 }
742 break;
7fd59977 743
744 case IntAna_NoGeometricSolution:
745 {
746 Standard_Boolean bReverse;
747 Standard_Real U1,V1,U2,V2;
748 gp_Pnt psol;
749 //
750 bReverse=IsToReverse(Cy1, Cy2, Tol);
ecc4f148 751 if (bReverse)
752 {
753 Cy2=Quad1.Cylinder();
754 Cy1=Quad2.Cylinder();
7fd59977 755 }
756 //
757 IntAna_IntQuadQuad anaint(Cy1,Cy2,Tol);
ecc4f148 758 if (!anaint.IsDone())
759 {
760 return Standard_False;
7fd59977 761 }
762
ecc4f148 763 if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0)
764 {
765 Empty = Standard_True;
7fd59977 766 }
ecc4f148 767 else
768 {
769 NbSol = anaint.NbPnt();
770 for (i = 1; i <= NbSol; i++)
771 {
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);
777 spnt.Append(ptsol);
778 }
779
780 gp_Pnt ptvalid, ptf, ptl;
781 gp_Vec tgvalid;
782
783 Standard_Real first,last,para;
784 IntAna_Curve curvsol;
785 Standard_Boolean tgfound;
786 Standard_Integer kount;
787
788 NbSol = anaint.NbCurve();
789 for (i = 1; i <= NbSol; i++)
790 {
791 curvsol = anaint.Curve(i);
792 curvsol.Domain(first,last);
793 ptf = curvsol.Value(first);
794 ptl = curvsol.Value(last);
795
796 para = last;
797 kount = 1;
798 tgfound = Standard_False;
799
800 while (!tgfound)
801 {
802 para = (1.123*first + para)/2.123;
803 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
804 if (!tgfound)
805 {
806 kount ++;
807 tgfound = kount > 5;
808 }
809 }
810
811 Handle(IntPatch_ALine) alig;
812 if (kount <= 5)
813 {
814 Standard_Real qwe = tgvalid.DotCross( Quad2.Normale(ptvalid),
815 Quad1.Normale(ptvalid));
816 if(qwe>0.00000001)
817 {
818 trans1 = IntSurf_Out;
819 trans2 = IntSurf_In;
820 }
821 else if(qwe<-0.00000001)
822 {
823 trans1 = IntSurf_In;
824 trans2 = IntSurf_Out;
825 }
826 else
827 {
828 trans1=trans2=IntSurf_Undecided;
829 }
830 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
831 }
832 else
833 {
834 alig = new IntPatch_ALine(curvsol,Standard_False);
835 //-- cout << "Transition indeterminee" << endl;
836 }
837
838 Standard_Boolean TempFalse1 = Standard_False;
839 Standard_Boolean TempFalse2 = Standard_False;
840
841 ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1,ptf,first,
842 TempFalse2,ptl,last,Multpoint,Tol);
843 slin.Append(alig);
844 }
7fd59977 845 }
846 }
847 break;
848
849 default:
ecc4f148 850 return Standard_False;
851 }
852
853 return Standard_True;
854}
855
856//=======================================================================
857//function : ShortCosForm
858//purpose : Represents theCosFactor*cosA+theSinFactor*sinA as
859// theCoeff*cos(A-theAngle) if it is possibly (all angles are
860// in radians).
861//=======================================================================
862static void ShortCosForm( const Standard_Real theCosFactor,
863 const Standard_Real theSinFactor,
864 Standard_Real& theCoeff,
865 Standard_Real& theAngle)
866{
867 theCoeff = sqrt(theCosFactor*theCosFactor+theSinFactor*theSinFactor);
868 theAngle = 0.0;
869 if(IsEqual(theCoeff, 0.0))
870 {
871 theAngle = 0.0;
872 return;
873 }
874
875 theAngle = acos(Abs(theCosFactor/theCoeff));
876
877 if(theSinFactor > 0.0)
878 {
879 if(IsEqual(theCosFactor, 0.0))
880 {
881 theAngle = M_PI/2.0;
882 }
883 else if(theCosFactor < 0.0)
884 {
885 theAngle = M_PI-theAngle;
886 }
887 }
888 else if(IsEqual(theSinFactor, 0.0))
889 {
890 if(theCosFactor < 0.0)
891 {
892 theAngle = M_PI;
893 }
894 }
895 if(theSinFactor < 0.0)
896 {
897 if(theCosFactor > 0.0)
898 {
899 theAngle = 2.0*M_PI-theAngle;
900 }
901 else if(IsEqual(theCosFactor, 0.0))
902 {
903 theAngle = 3.0*M_PI/2.0;
904 }
905 else if(theCosFactor < 0.0)
906 {
907 theAngle = M_PI+theAngle;
908 }
909 }
910}
911
912enum SearchBoundType
913{
914 SearchNONE = 0,
915 SearchV1 = 1,
916 SearchV2 = 2
917};
918
919//Stores equations coefficients
920struct stCoeffsValue
921{
922 stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&);
923
b5ef9d91 924 math_Vector mVecA1;
925 math_Vector mVecA2;
926 math_Vector mVecB1;
927 math_Vector mVecB2;
928 math_Vector mVecC1;
929 math_Vector mVecC2;
930 math_Vector mVecD;
ecc4f148 931
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
937
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
943
944 Standard_Real mK1;
945 Standard_Real mL1;
946 Standard_Real mK2;
947 Standard_Real mL2;
948
949 Standard_Real mFIV1;
950 Standard_Real mPSIV1;
951 Standard_Real mFIV2;
952 Standard_Real mPSIV2;
953
954 Standard_Real mB;
955 Standard_Real mC;
956 Standard_Real mFI1;
957 Standard_Real mFI2;
958};
959
960stCoeffsValue::stCoeffsValue( const gp_Cylinder& theCyl1,
b5ef9d91 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())
ecc4f148 969{
ecc4f148 970 enum CoupleOfEquation
971 {
972 COENONE = 0,
973 COE12 = 1,
974 COE23 = 2,
975 COE13 = 3
976 }aFoundCouple = COENONE;
977
978
7c32c7c4 979 Standard_Real aDetV1V2 = 0.0;
980
b5ef9d91 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
7c32c7c4 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
987
988 if(anAbsD1 >= anAbsD2)
ecc4f148 989 {
7c32c7c4 990 if(anAbsD3 > anAbsD1)
ecc4f148 991 {
7c32c7c4 992 aFoundCouple = COE13;
993 aDetV1V2 = aDelta3;
ecc4f148 994 }
995 else
996 {
7c32c7c4 997 aFoundCouple = COE12;
998 aDetV1V2 = aDelta1;
ecc4f148 999 }
1000 }
1001 else
1002 {
7c32c7c4 1003 if(anAbsD3 > anAbsD2)
1004 {
1005 aFoundCouple = COE13;
1006 aDetV1V2 = aDelta3;
1007 }
1008 else
1009 {
1010 aFoundCouple = COE23;
1011 aDetV1V2 = aDelta2;
1012 }
1013 }
1014
1015 if(Abs(aDetV1V2) < aNulValue)
1016 {
1017 Standard_Failure::Raise("Error. Exception in divide by zerro (IntCyCyTrim)!!!!");
ecc4f148 1018 }
1019
1020 switch(aFoundCouple)
1021 {
1022 case COE12:
1023 break;
1024 case COE23:
1025 {
b5ef9d91 1026 math_Vector aVTemp(mVecA1);
1027 mVecA1(1) = aVTemp(2);
1028 mVecA1(2) = aVTemp(3);
1029 mVecA1(3) = aVTemp(1);
ecc4f148 1030
1031 aVTemp = mVecA2;
b5ef9d91 1032 mVecA2(1) = aVTemp(2);
1033 mVecA2(2) = aVTemp(3);
1034 mVecA2(3) = aVTemp(1);
ecc4f148 1035
1036 aVTemp = mVecB1;
b5ef9d91 1037 mVecB1(1) = aVTemp(2);
1038 mVecB1(2) = aVTemp(3);
1039 mVecB1(3) = aVTemp(1);
ecc4f148 1040
1041 aVTemp = mVecB2;
b5ef9d91 1042 mVecB2(1) = aVTemp(2);
1043 mVecB2(2) = aVTemp(3);
1044 mVecB2(3) = aVTemp(1);
ecc4f148 1045
1046 aVTemp = mVecC1;
b5ef9d91 1047 mVecC1(1) = aVTemp(2);
1048 mVecC1(2) = aVTemp(3);
1049 mVecC1(3) = aVTemp(1);
ecc4f148 1050
1051 aVTemp = mVecC2;
b5ef9d91 1052 mVecC2(1) = aVTemp(2);
1053 mVecC2(2) = aVTemp(3);
1054 mVecC2(3) = aVTemp(1);
ecc4f148 1055
1056 aVTemp = mVecD;
b5ef9d91 1057 mVecD(1) = aVTemp(2);
1058 mVecD(2) = aVTemp(3);
1059 mVecD(3) = aVTemp(1);
ecc4f148 1060
1061 break;
1062 }
1063 case COE13:
1064 {
b5ef9d91 1065 math_Vector aVTemp = mVecA1;
1066 mVecA1(2) = aVTemp(3);
1067 mVecA1(3) = aVTemp(2);
ecc4f148 1068
1069 aVTemp = mVecA2;
b5ef9d91 1070 mVecA2(2) = aVTemp(3);
1071 mVecA2(3) = aVTemp(2);
ecc4f148 1072
1073 aVTemp = mVecB1;
b5ef9d91 1074 mVecB1(2) = aVTemp(3);
1075 mVecB1(3) = aVTemp(2);
ecc4f148 1076
1077 aVTemp = mVecB2;
b5ef9d91 1078 mVecB2(2) = aVTemp(3);
1079 mVecB2(3) = aVTemp(2);
ecc4f148 1080
1081 aVTemp = mVecC1;
b5ef9d91 1082 mVecC1(2) = aVTemp(3);
1083 mVecC1(3) = aVTemp(2);
ecc4f148 1084
1085 aVTemp = mVecC2;
b5ef9d91 1086 mVecC2(2) = aVTemp(3);
1087 mVecC2(3) = aVTemp(2);
ecc4f148 1088
1089 aVTemp = mVecD;
b5ef9d91 1090 mVecD(2) = aVTemp(3);
1091 mVecD(3) = aVTemp(2);
ecc4f148 1092
1093 break;
1094 }
1095 default:
1096 break;
1097 }
1098
1099 //------- For V1 (begin)
1100 //sinU2
b5ef9d91 1101 mK21 = (mVecC2(2)*mVecB2(1)-mVecC2(1)*mVecB2(2))/aDetV1V2;
ecc4f148 1102 //sinU1
b5ef9d91 1103 mK11 = (mVecC2(2)*mVecB1(1)-mVecC2(1)*mVecB1(2))/aDetV1V2;
ecc4f148 1104 //cosU2
b5ef9d91 1105 mL21 = (mVecC2(2)*mVecA2(1)-mVecC2(1)*mVecA2(2))/aDetV1V2;
ecc4f148 1106 //cosU1
b5ef9d91 1107 mL11 = (mVecC2(2)*mVecA1(1)-mVecC2(1)*mVecA1(2))/aDetV1V2;
ecc4f148 1108 //Free member
b5ef9d91 1109 mM1 = (mVecC2(2)*mVecD(1)-mVecC2(1)*mVecD(2))/aDetV1V2;
ecc4f148 1110 //------- For V1 (end)
1111
1112 //------- For V2 (begin)
1113 //sinU2
b5ef9d91 1114 mK22 = (mVecC1(1)*mVecB2(2)-mVecC1(2)*mVecB2(1))/aDetV1V2;
ecc4f148 1115 //sinU1
b5ef9d91 1116 mK12 = (mVecC1(1)*mVecB1(2)-mVecC1(2)*mVecB1(1))/aDetV1V2;
ecc4f148 1117 //cosU2
b5ef9d91 1118 mL22 = (mVecC1(1)*mVecA2(2)-mVecC1(2)*mVecA2(1))/aDetV1V2;
ecc4f148 1119 //cosU1
b5ef9d91 1120 mL12 = (mVecC1(1)*mVecA1(2)-mVecC1(2)*mVecA1(1))/aDetV1V2;
ecc4f148 1121 //Free member
b5ef9d91 1122 mM2 = (mVecC1(1)*mVecD(2)-mVecC1(2)*mVecD(1))/aDetV1V2;
ecc4f148 1123 //------- For V1 (end)
1124
1125 ShortCosForm(mL11, mK11, mK1, mFIV1);
1126 ShortCosForm(mL21, mK21, mL1, mPSIV1);
1127 ShortCosForm(mL12, mK12, mK2, mFIV2);
1128 ShortCosForm(mL22, mK22, mL2, mPSIV2);
1129
b5ef9d91 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
ecc4f148 1134
b5ef9d91 1135 mC =mVecD(3) - mVecC1(3)*mM1 -mVecC2(3)*mM2; //Free
ecc4f148 1136
1137 Standard_Real aA = 0.0;
1138
1139 ShortCosForm(aB2,aB1,mB,mFI1);
1140 ShortCosForm(aA2,aA1,aA,mFI2);
1141
1142 mB /= aA;
1143 mC /= aA;
1144}
1145
b5ef9d91 1146//=======================================================================
1147//function : CylCylMonotonicity
1148//purpose : Determines, if U2(U1) function is increasing.
1149//=======================================================================
1150static 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)
1155{
1156 // U2(U1) = FI2 + (+/-)acos(B*cos(U1 - FI1) + C);
1157 //Accordingly,
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).
1163 //
1164 //Consequently,
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.
1176
1177 //For "+/-" sign. If isPlus == TRUE, "+" is chosen, otherwise, we choose "-".
1178 Standard_Boolean isPlus = Standard_False;
1179
1180 switch(theWLIndex)
1181 {
1182 case 0:
1183 isPlus = Standard_True;
1184 break;
1185 case 1:
1186 isPlus = Standard_False;
1187 break;
1188 default:
1189 //Standard_Failure::Raise("Error. Range Error!!!!");
1190 return Standard_False;
1191 }
1192
1193 Standard_Real aU1Temp = theU1par - theCoeffs.mFI1;
1194 InscribePoint(0, thePeriod, aU1Temp, 0.0, thePeriod, Standard_False);
1195
1196 theIsIncreasing = Standard_True;
1197
1198 if(((M_PI - aU1Temp) < RealSmall()) && (aU1Temp < thePeriod))
1199 {
1200 theIsIncreasing = Standard_False;
1201 }
1202
1203 if(theCoeffs.mB < 0.0)
1204 {
1205 theIsIncreasing = !theIsIncreasing;
1206 }
1207
1208 if(!isPlus)
1209 {
1210 theIsIncreasing = !theIsIncreasing;
1211 }
1212
1213 return Standard_True;
1214}
1215
1216static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
1217 const Standard_Integer theWLIndex,
1218 const stCoeffsValue& theCoeffs,
1219 Standard_Real& theU2)
1220{
1221 Standard_Real aSign = 0.0;
1222
1223 switch(theWLIndex)
1224 {
1225 case 0:
1226 aSign = 1.0;
1227 break;
1228 case 1:
1229 aSign = -1.0;
1230 break;
1231 default:
1232 //Standard_Failure::Raise("Error. Range Error!!!!");
1233 return Standard_False;
1234 }
1235
1236 Standard_Real anArg = theCoeffs.mB *
1237 cos(theU1par - theCoeffs.mFI1) +
1238 theCoeffs.mC;
1239
1240 if(aNulValue > 1.0 - anArg)
1241 anArg = 1.0;
1242
1243 if(anArg + 1.0 < aNulValue)
1244 anArg = -1.0;
1245
1246 theU2 = theCoeffs.mFI2 + aSign*acos(anArg);
1247
1248 return Standard_True;
1249}
1250
1251static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1,
1252 const Standard_Real theU2,
1253 const stCoeffsValue& theCoeffs,
1254 Standard_Real& theV1,
1255 Standard_Real& theV2)
1256{
1257 theV1 = theCoeffs.mK21 * sin(theU2) +
1258 theCoeffs.mK11 * sin(theU1) +
1259 theCoeffs.mL21 * cos(theU2) +
1260 theCoeffs.mL11 * cos(theU1) + theCoeffs.mM1;
1261
1262 theV2 = theCoeffs.mK22 * sin(theU2) +
1263 theCoeffs.mK12 * sin(theU1) +
1264 theCoeffs.mL22 * cos(theU2) +
1265 theCoeffs.mL12 * cos(theU1) + theCoeffs.mM2;
1266
1267 return Standard_True;
1268}
1269
1270
1271static 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)
1277{
1278 if(!CylCylComputeParameters(theU1par, theWLIndex, theCoeffs, theU2))
1279 return Standard_False;
1280
1281 if(!CylCylComputeParameters(theU1par, theU2, theCoeffs, theV1, theV2))
1282 return Standard_False;
1283
1284 return Standard_True;
1285}
1286
ecc4f148 1287//=======================================================================
1288//function : SearchOnVBounds
1289//purpose :
1290//=======================================================================
1291static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType,
1292 const stCoeffsValue& theCoeffs,
1293 const Standard_Real theVzad,
b5ef9d91 1294 const Standard_Real theVInit,
ecc4f148 1295 const Standard_Real theInitU2,
1296 const Standard_Real theInitMainVar,
1297 Standard_Real& theMainVariableValue)
1298{
b5ef9d91 1299 const Standard_Integer aNbDim = 3;
ecc4f148 1300 const Standard_Real aMaxError = 4.0*M_PI; // two periods
1301
b5ef9d91 1302 theMainVariableValue = theInitMainVar;
1303 const Standard_Real aTol2 = 1.0e-18;
1304 Standard_Real aMainVarPrev = theInitMainVar, aU2Prev = theInitU2, anOtherVar = theVInit;
ecc4f148 1305
b5ef9d91 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);
1310
ecc4f148 1311 Standard_Real anError = RealLast();
b5ef9d91 1312 Standard_Real anErrorPrev = anError;
ecc4f148 1313 Standard_Integer aNbIter = 0;
1314 do
1315 {
1316 if(++aNbIter > 1000)
1317 return Standard_False;
1318
b5ef9d91 1319 const Standard_Real aSinU1 = sin(aMainVarPrev),
1320 aCosU1 = cos(aMainVarPrev),
1321 aSinU2 = sin(aU2Prev),
1322 aCosU2 = cos(aU2Prev);
ecc4f148 1323
b5ef9d91 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 +
1332 theCoeffs.mVecD;
1333
1334 math_Vector aMSum(1, 3);
ecc4f148 1335
1336 switch(theSBType)
1337 {
1338 case SearchV1:
b5ef9d91 1339 aMatr.SetCol(3, theCoeffs.mVecC2);
1340 aMSum = theCoeffs.mVecC1 * theVzad;
1341 aVecFreeMem -= aMSum;
1342 aMSum += theCoeffs.mVecC2*anOtherVar;
ecc4f148 1343 break;
1344
1345 case SearchV2:
b5ef9d91 1346 aMatr.SetCol(3, theCoeffs.mVecC1);
1347 aMSum = theCoeffs.mVecC2 * theVzad;
1348 aVecFreeMem -= aMSum;
1349 aMSum += theCoeffs.mVecC1*anOtherVar;
ecc4f148 1350 break;
1351
1352 default:
1353 return Standard_False;
1354 }
1355
b5ef9d91 1356 aMatr.SetCol(1, theCoeffs.mVecA1 * aSinU1 - theCoeffs.mVecB1 * aCosU1);
1357 aMatr.SetCol(2, theCoeffs.mVecA2 * aSinU2 - theCoeffs.mVecB2 * aCosU2);
1358
1359 Standard_Real aDetMainSyst = aMatr.Determinant();
ecc4f148 1360
b5ef9d91 1361 if(Abs(aDetMainSyst) < aNulValue)
ecc4f148 1362 {
1363 return Standard_False;
1364 }
1365
b5ef9d91 1366 math_Matrix aM1(aMatr), aM2(aMatr), aM3(aMatr);
1367 aM1.SetCol(1, aVecFreeMem);
1368 aM2.SetCol(2, aVecFreeMem);
1369 aM3.SetCol(3, aVecFreeMem);
ecc4f148 1370
b5ef9d91 1371 const Standard_Real aDetMainVar = aM1.Determinant();
1372 const Standard_Real aDetVar1 = aM2.Determinant();
1373 const Standard_Real aDetVar2 = aM3.Determinant();
ecc4f148 1374
1375 Standard_Real aDelta = aDetMainVar/aDetMainSyst-aMainVarPrev;
1376
1377 if(Abs(aDelta) > aMaxError)
1378 return Standard_False;
1379
1380 anError = aDelta*aDelta;
1381 aMainVarPrev += aDelta;
1382
1383 ///
1384 aDelta = aDetVar1/aDetMainSyst-aU2Prev;
1385
1386 if(Abs(aDelta) > aMaxError)
1387 return Standard_False;
1388
1389 anError += aDelta*aDelta;
1390 aU2Prev += aDelta;
1391
1392 ///
1393 aDelta = aDetVar2/aDetMainSyst-anOtherVar;
1394 anError += aDelta*aDelta;
1395 anOtherVar += aDelta;
b5ef9d91 1396
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 +
1407 theCoeffs.mVecD);
1408 const Standard_Real aSQNorm = aMSum.Norm2();
1409 return (aSQNorm < aTol2);
1410 }
1411 else
1412 {
1413 theMainVariableValue = aMainVarPrev;
1414 }
1415
1416 anErrorPrev = anError;
ecc4f148 1417 }
1418 while(anError > aTol2);
1419
1420 theMainVariableValue = aMainVarPrev;
1421
1422 return Standard_True;
1423}
1424
1425//=======================================================================
1426//function : InscribePoint
e002f1ce 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).
ecc4f148 1431//=======================================================================
b5ef9d91 1432Standard_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)
ecc4f148 1438{
02effd35 1439 if((theUfTarget - theUGiven <= theTol2D) &&
1440 (theUGiven - theUlTarget <= theTol2D))
ecc4f148 1441 {//it has already inscribed
1442
1443 /*
1444 Utf U Utl
1445 + * +
1446 */
1447
02effd35 1448 if(theFlForce)
1449 {
1450 Standard_Real anUtemp = theUGiven + thePeriod;
1451 if((theUfTarget - anUtemp <= theTol2D) &&
1452 (anUtemp - theUlTarget <= theTol2D))
1453 {
1454 theUGiven = anUtemp;
1455 return Standard_True;
1456 }
1457
1458 anUtemp = theUGiven - thePeriod;
1459 if((theUfTarget - anUtemp <= theTol2D) &&
1460 (anUtemp - theUlTarget <= theTol2D))
1461 {
1462 theUGiven = anUtemp;
1463 }
1464 }
1465
ecc4f148 1466 return Standard_True;
1467 }
1468
1469 if(IsEqual(thePeriod, 0.0))
1470 {//it cannot be inscribed
1471 return Standard_False;
1472 }
1473
1474 Standard_Real aDU = theUGiven - theUfTarget;
1475
1476 if(aDU > 0.0)
1477 aDU = -thePeriod;
1478 else
1479 aDU = +thePeriod;
1480
1481 while(((theUGiven - theUfTarget)*aDU < 0.0) && !((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D)))
1482 {
1483 theUGiven += aDU;
1484 }
1485
1486 return ((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D));
1487}
1488
1489//=======================================================================
1490//function : InscribeInterval
1491//purpose : Adjusts theUfGiven and after that fits theUlGiven to result
1492//=======================================================================
1493static 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)
1499{
1500 Standard_Real anUpar = theUfGiven;
1501 const Standard_Real aDelta = theUlGiven - theUfGiven;
02effd35 1502 if(InscribePoint(theUfTarget, theUlTarget, anUpar,
1503 theTol2D, thePeriod, Standard_False))
ecc4f148 1504 {
1505 theUfGiven = anUpar;
1506 theUlGiven = theUfGiven + aDelta;
1507 }
1508 else
1509 {
1510 anUpar = theUlGiven;
02effd35 1511 if(InscribePoint(theUfTarget, theUlTarget, anUpar,
1512 theTol2D, thePeriod, Standard_False))
ecc4f148 1513 {
1514 theUlGiven = anUpar;
1515 theUfGiven = theUlGiven - aDelta;
1516 }
1517 else
7fd59977 1518 {
1519 return Standard_False;
1520 }
1521 }
ecc4f148 1522
1523 return Standard_True;
1524}
1525
1526//=======================================================================
1527//function : InscribeAndSortArray
1528//purpose : Sort from Min to Max value
1529//=======================================================================
1530static 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)
1536{
1537 for(Standard_Integer i = 0; i < theNOfMembers; i++)
1538 {
1539 if(Precision::IsInfinite(theArr[i]))
1540 {
1541 if(theArr[i] < 0.0)
1542 theArr[i] = -theArr[i];
1543
1544 continue;
1545 }
1546
02effd35 1547 InscribePoint(theUf, theUl, theArr[i], theTol2D, thePeriod, Standard_False);
ecc4f148 1548
1549 for(Standard_Integer j = i - 1; j >= 0; j--)
1550 {
1551
1552 if(theArr[j+1] < theArr[j])
1553 {
1554 Standard_Real anUtemp = theArr[j+1];
1555 theArr[j+1] = theArr[j];
1556 theArr[j] = anUtemp;
1557 }
1558 }
1559 }
1560}
1561
1562
1563//=======================================================================
1564//function : AddPointIntoWL
1565//purpose : Surf1 is a surface, whose U-par is variable.
1566//=======================================================================
1567static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
1568 const IntSurf_Quadric& theQuad2,
b5ef9d91 1569 const stCoeffsValue& theCoeffs,
ecc4f148 1570 const Standard_Boolean isTheReverse,
b5ef9d91 1571 const Standard_Boolean isThePrecise,
ecc4f148 1572 const gp_Pnt2d& thePntOnSurf1,
1573 const gp_Pnt2d& thePntOnSurf2,
1574 const Standard_Real theUfSurf1,
1575 const Standard_Real theUlSurf1,
b5ef9d91 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,
ecc4f148 1582 const Standard_Real thePeriodOfSurf1,
1583 const Handle(IntSurf_LineOn2S)& theLine,
b5ef9d91 1584 const Standard_Integer theWLIndex,
e8feb725 1585 const Standard_Real theTol3D,
02effd35 1586 const Standard_Real theTol2D,
b5ef9d91 1587 const Standard_Boolean theFlForce,
1588 const Standard_Boolean theFlBefore = Standard_False)
ecc4f148 1589{
1590 gp_Pnt aPt1(theQuad1.Value(thePntOnSurf1.X(), thePntOnSurf1.Y())),
1591 aPt2(theQuad2.Value(thePntOnSurf2.X(), thePntOnSurf2.Y()));
1592
b5ef9d91 1593 Standard_Real aU1par = thePntOnSurf1.X();
1594 if(!InscribePoint(theUfSurf1, theUlSurf1, aU1par, theTol2D,
02effd35 1595 thePeriodOfSurf1, theFlForce))
ecc4f148 1596 return Standard_False;
1597
b5ef9d91 1598 Standard_Real aU2par = thePntOnSurf2.X();
1599 if(!InscribePoint(theUfSurf2, theUlSurf2, aU2par, theTol2D,
1600 thePeriodOfSurf1, Standard_False))
1601 return Standard_False;
1602
1603 Standard_Real aV1par = thePntOnSurf1.Y();
1604 if((aV1par - theVlSurf1 > theTol2D) || (theVfSurf1 - aV1par > theTol2D))
1605 return Standard_False;
1606
1607 Standard_Real aV2par = thePntOnSurf2.Y();
1608 if((aV2par - theVlSurf2 > theTol2D) || (theVfSurf2 - aV2par > theTol2D))
1609 return Standard_False;
1610
ecc4f148 1611 IntSurf_PntOn2S aPnt;
1612
1613 if(isTheReverse)
1614 {
1615 aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
b5ef9d91 1616 aU2par, aV2par,
1617 aU1par, aV1par);
ecc4f148 1618 }
1619 else
1620 {
1621 aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
b5ef9d91 1622 aU1par, aV1par,
1623 aU2par, aV2par);
ecc4f148 1624 }
1625
b5ef9d91 1626 Standard_Integer aNbPnts = theLine->NbPoints();
e8feb725 1627 if(aNbPnts > 0)
1628 {
1629 Standard_Real aUl = 0.0, aVl = 0.0;
1630 const IntSurf_PntOn2S aPlast = theLine->Value(aNbPnts);
1631 if(isTheReverse)
1632 aPlast.ParametersOnS2(aUl, aVl);
1633 else
1634 aPlast.ParametersOnS1(aUl, aVl);
1635
b5ef9d91 1636 if(!theFlBefore && (aU1par <= aUl))
1637 {//Parameter value must be increased if theFlBefore == FALSE.
e8feb725 1638 return Standard_False;
1639 }
1640
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))
1647 {
1648 theLine->RemovePoint(aNbPnts);
1649 }
1650 }
1651
ecc4f148 1652 theLine->Add(aPnt);
b5ef9d91 1653
1654 if(!isThePrecise)
1655 return Standard_True;
1656
1657 //Try to precise existing WLine
1658 aNbPnts = theLine->NbPoints();
1659 if(aNbPnts >= 3)
1660 {
1661 Standard_Real aU1 = 0.0, aU2 = 0.0, aU3 = 0.0, aV = 0.0;
1662 if(isTheReverse)
1663 {
1664 theLine->Value(aNbPnts).ParametersOnS2(aU3, aV);
1665 theLine->Value(aNbPnts-1).ParametersOnS2(aU2, aV);
1666 theLine->Value(aNbPnts-2).ParametersOnS2(aU1, aV);
1667 }
1668 else
1669 {
1670 theLine->Value(aNbPnts).ParametersOnS1(aU3, aV);
1671 theLine->Value(aNbPnts-1).ParametersOnS1(aU2, aV);
1672 theLine->Value(aNbPnts-2).ParametersOnS1(aU1, aV);
1673 }
1674
1675 const Standard_Real aStepPrev = aU2-aU1;
1676 const Standard_Real aStep = aU3-aU2;
1677
1678 const Standard_Integer aDeltaStep = RealToInt(aStepPrev/aStep);
1679
1680 if((1 < aDeltaStep) && (aDeltaStep < 2000))
1681 {
1682 SeekAdditionalPoints( theQuad1, theQuad2, theLine,
1683 theCoeffs, theWLIndex, aDeltaStep, aNbPnts-2,
1684 aNbPnts-1, theUfSurf2, theUlSurf2,
1685 theTol2D, thePeriodOfSurf1, isTheReverse);
1686 }
1687 }
1688
ecc4f148 1689 return Standard_True;
1690}
1691
1692//=======================================================================
1693//function : AddBoundaryPoint
1694//purpose :
1695//=======================================================================
1696static 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,
e8feb725 1702 const Standard_Real theTol3D,
ecc4f148 1703 const Standard_Real theTol2D,
1704 const Standard_Real thePeriod,
ecc4f148 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,
b5ef9d91 1711 const Standard_Integer theWLIndex,
ecc4f148 1712 const Standard_Boolean isTheReverse,
02effd35 1713 const Standard_Boolean theFlForce,
ecc4f148 1714 Standard_Boolean& isTheFound1,
1715 Standard_Boolean& isTheFound2)
1716{
1717 Standard_Real aUSurf1f = 0.0, //const
1718 aUSurf1l = 0.0,
1719 aVSurf1f = 0.0,
1720 aVSurf1l = 0.0;
1721 Standard_Real aUSurf2f = 0.0, //const
1722 aUSurf2l = 0.0,
1723 aVSurf2f = 0.0,
1724 aVSurf2l = 0.0;
1725
1726 theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
1727 theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
1728
1729 SearchBoundType aTS1 = SearchNONE, aTS2 = SearchNONE;
1730 Standard_Real aV1zad = aVSurf1f, aV2zad = aVSurf2f;
1731
1732 Standard_Real anUpar1 = theU1, anUpar2 = theU1;
1733 Standard_Real aVf = theV1, aVl = theV1Prev;
b5ef9d91 1734
1735 if( (Abs(aVf-aVSurf1f) <= theTol2D) ||
1736 ((aVf-aVSurf1f)*(aVl-aVSurf1f) <= 0.0))
ecc4f148 1737 {
1738 aTS1 = SearchV1;
1739 aV1zad = aVSurf1f;
b5ef9d91 1740 isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1f, theV2, theU2, theU1, anUpar1);
ecc4f148 1741 }
b5ef9d91 1742 else if((Abs(aVf-aVSurf1l) <= theTol2D) ||
1743 ((aVf-aVSurf1l)*(aVl-aVSurf1l) <= 0.0))
ecc4f148 1744 {
1745 aTS1 = SearchV1;
1746 aV1zad = aVSurf1l;
b5ef9d91 1747 isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1l, theV2, theU2, theU1, anUpar1);
ecc4f148 1748 }
1749
1750 aVf = theV2;
1751 aVl = theV2Prev;
ecc4f148 1752
b5ef9d91 1753 if((Abs(aVf-aVSurf2f) <= theTol2D) ||
1754 ((aVf-aVSurf2f)*(aVl-aVSurf2f) <= 0.0))
ecc4f148 1755 {
1756 aTS2 = SearchV2;
1757 aV2zad = aVSurf2f;
b5ef9d91 1758 isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2f, theV1, theU2, theU1, anUpar2);
ecc4f148 1759 }
b5ef9d91 1760 else if((Abs(aVf-aVSurf2l) <= theTol2D) ||
1761 ((aVf-aVSurf2l)*(aVl-aVSurf2l) <= 0.0))
ecc4f148 1762 {
1763 aTS2 = SearchV2;
1764 aV2zad = aVSurf2l;
b5ef9d91 1765 isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2l, theV1, theU2, theU1, anUpar2);
ecc4f148 1766 }
1767
b5ef9d91 1768 if(!isTheFound1 && !isTheFound2)
1769 return Standard_True;
1770
ecc4f148 1771 if(anUpar1 < anUpar2)
1772 {
1773 if(isTheFound1)
1774 {
b5ef9d91 1775 Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
1776 if(!CylCylComputeParameters(anUpar1, theWLIndex, theCoeffs, aU2, aV1, aV2))
ecc4f148 1777 {
b5ef9d91 1778 isTheFound1 = Standard_False;
ecc4f148 1779 }
b5ef9d91 1780
1781 if(aTS1 == SearchV1)
1782 aV1 = aV1zad;
1783
1784 if(aTS1 == SearchV2)
1785 aV2 = aV2zad;
1786
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))
ecc4f148 1793 {
1794 isTheFound1 = Standard_False;
1795 }
1796 }
1797
1798 if(isTheFound2)
1799 {
b5ef9d91 1800 Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
1801 if(!CylCylComputeParameters(anUpar2, theWLIndex, theCoeffs, aU2, aV1, aV2))
ecc4f148 1802 {
b5ef9d91 1803 isTheFound2 = Standard_False;
ecc4f148 1804 }
b5ef9d91 1805
1806 if(aTS2 == SearchV1)
1807 aV1 = aV1zad;
1808
1809 if(aTS2 == SearchV2)
1810 aV2 = aV2zad;
1811
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))
ecc4f148 1818 {
1819 isTheFound2 = Standard_False;
1820 }
1821 }
1822 }
1823 else
1824 {
1825 if(isTheFound2)
1826 {
b5ef9d91 1827 Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
1828 if(!CylCylComputeParameters(anUpar2, theWLIndex, theCoeffs, aU2, aV1, aV2))
ecc4f148 1829 {
b5ef9d91 1830 isTheFound2 = Standard_False;
ecc4f148 1831 }
b5ef9d91 1832
1833 if(aTS2 == SearchV1)
1834 aV1 = aV1zad;
1835
1836 if(aTS2 == SearchV2)
1837 aV2 = aV2zad;
1838
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))
ecc4f148 1845 {
1846 isTheFound2 = Standard_False;
1847 }
1848 }
1849
1850 if(isTheFound1)
1851 {
b5ef9d91 1852 Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
1853 if(!CylCylComputeParameters(anUpar1, theWLIndex, theCoeffs, aU2, aV1, aV2))
ecc4f148 1854 {
b5ef9d91 1855 isTheFound1 = Standard_False;
ecc4f148 1856 }
b5ef9d91 1857
1858 if(aTS1 == SearchV1)
1859 aV1 = aV1zad;
1860
1861 if(aTS1 == SearchV2)
1862 aV2 = aV2zad;
1863
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))
ecc4f148 1870 {
1871 isTheFound1 = Standard_False;
1872 }
1873 }
1874 }
1875
7fd59977 1876 return Standard_True;
1877}
1878
ecc4f148 1879//=======================================================================
1880//function : SeekAdditionalPoints
1881//purpose :
1882//=======================================================================
1883static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
1884 const IntSurf_Quadric& theQuad2,
b5ef9d91 1885 const Handle(IntSurf_LineOn2S)& theLine,
ecc4f148 1886 const stCoeffsValue& theCoeffs,
b5ef9d91 1887 const Standard_Integer theWLIndex,
ecc4f148 1888 const Standard_Integer theMinNbPoints,
b5ef9d91 1889 const Standard_Integer theStartPointOnLine,
1890 const Standard_Integer theEndPointOnLine,
ecc4f148 1891 const Standard_Real theU2f,
1892 const Standard_Real theU2l,
1893 const Standard_Real theTol2D,
1894 const Standard_Real thePeriodOfSurf2,
ecc4f148 1895 const Standard_Boolean isTheReverse)
1896{
b5ef9d91 1897 if(theLine.IsNull())
1898 return;
1899
1900 Standard_Integer aNbPoints = theEndPointOnLine - theStartPointOnLine + 1;
ecc4f148 1901 if(aNbPoints >= theMinNbPoints)
1902 {
1903 return;
1904 }
1905
b5ef9d91 1906 Standard_Real aMinDeltaParam = theTol2D;
1907
1908 {
1909 Standard_Real u1 = 0.0, v1 = 0.0, u2 = 0.0, v2 = 0.0;
1910
1911 if(isTheReverse)
1912 {
1913 theLine->Value(theStartPointOnLine).ParametersOnS2(u1, v1);
1914 theLine->Value(theEndPointOnLine).ParametersOnS2(u2, v2);
1915 }
1916 else
1917 {
1918 theLine->Value(theStartPointOnLine).ParametersOnS1(u1, v1);
1919 theLine->Value(theEndPointOnLine).ParametersOnS1(u2, v2);
1920 }
1921
1922 aMinDeltaParam = Max(Abs(u2 - u1)/IntToReal(theMinNbPoints), aMinDeltaParam);
1923 }
1924
1925 Standard_Integer aLastPointIndex = theEndPointOnLine;
ecc4f148 1926 Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0;
1927
1928 Standard_Integer aNbPointsPrev = 0;
1929 while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
1930 {
1931 aNbPointsPrev = aNbPoints;
b5ef9d91 1932 for(Standard_Integer fp = theStartPointOnLine, lp = 0; fp < aLastPointIndex; fp = lp + 1)
ecc4f148 1933 {
02effd35 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
1936
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
ecc4f148 1939
1940 lp = fp+1;
1941
1942 if(isTheReverse)
1943 {
b5ef9d91 1944 theLine->Value(fp).ParametersOnS2(U1f, V1f);
1945 theLine->Value(lp).ParametersOnS2(U1l, V1l);
02effd35 1946
b5ef9d91 1947 theLine->Value(fp).ParametersOnS1(U2f, V2f);
1948 theLine->Value(lp).ParametersOnS1(U2l, V2l);
ecc4f148 1949 }
1950 else
1951 {
b5ef9d91 1952 theLine->Value(fp).ParametersOnS1(U1f, V1f);
1953 theLine->Value(lp).ParametersOnS1(U1l, V1l);
02effd35 1954
b5ef9d91 1955 theLine->Value(fp).ParametersOnS2(U2f, V2f);
1956 theLine->Value(lp).ParametersOnS2(U2l, V2l);
02effd35 1957 }
1958
b5ef9d91 1959 if(Abs(U1l - U1f) <= aMinDeltaParam)
02effd35 1960 {
1961 //Step is minimal. It is not necessary to divide it.
1962 continue;
ecc4f148 1963 }
1964
1965 U1prec = 0.5*(U1f+U1l);
1966
b5ef9d91 1967 if(!CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec))
1968 continue;
ecc4f148 1969
02effd35 1970 InscribePoint(theU2f, theU2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False);
ecc4f148 1971
b5ef9d91 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()));
ecc4f148 1974
7c32c7c4 1975#ifdef OCCT_DEBUG
1976 //cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << endl;
1977#endif
1978
ecc4f148 1979 IntSurf_PntOn2S anIP;
1980 if(isTheReverse)
1981 {
1982 anIP.SetValue(aPInt, U2prec, V2prec, U1prec, V1prec);
1983 }
1984 else
1985 {
1986 anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec);
1987 }
1988
b5ef9d91 1989 theLine->InsertBefore(lp, anIP);
ecc4f148 1990
b5ef9d91 1991 aNbPoints++;
1992 aLastPointIndex++;
1993 }
1994
1995 if(aNbPoints >= theMinNbPoints)
1996 {
1997 return;
ecc4f148 1998 }
1999 }
2000}
2001
2002//=======================================================================
2003//function : CriticalPointsComputing
2004//purpose :
2005//=======================================================================
2006static 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[])
2015{
2016 theU1crit[0] = 0.0;
2017 theU1crit[1] = thePeriod;
2018 theU1crit[2] = theUSurf1f;
2019 theU1crit[3] = theUSurf1l;
2020
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))
2024 {
2025 Standard_Real anArg = (aCOS - theCoeffs.mC) / theCoeffs.mB;
2026 if(anArg > 1.0)
2027 anArg = 1.0;
2028 if(anArg < -1.0)
2029 anArg = -1.0;
2030
2031 theU1crit[4] = -acos(anArg) + theCoeffs.mFI1;
2032 theU1crit[5] = acos(anArg) + theCoeffs.mFI1;
2033 }
2034
2035 Standard_Real aSf = cos(theUSurf2f - theCoeffs.mFI2);
2036 Standard_Real aSl = cos(theUSurf2l - theCoeffs.mFI2);
2037 MinMax(aSf, aSl);
2038
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();
2043
2044 //preparative treatment of array
2045 InscribeAndSortArray(theU1crit, theNbCritPointsMax, 0.0, thePeriod, theTol2D, thePeriod);
2046 for(Standard_Integer i = 1; i < theNbCritPointsMax; i++)
2047 {
2048 Standard_Real &a = theU1crit[i],
2049 &b = theU1crit[i-1];
e8feb725 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))
ecc4f148 2052 {
2053 a = (a + b)/2.0;
2054 b = Precision::Infinite();
2055 }
2056 }
2057}
2058
2059//=======================================================================
2060//function : IntCyCyTrim
2061//purpose :
2062//=======================================================================
2063Standard_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)
2073{
2074 Standard_Real aUSurf1f = 0.0, //const
2075 aUSurf1l = 0.0,
2076 aVSurf1f = 0.0,
2077 aVSurf1l = 0.0;
2078 Standard_Real aUSurf2f = 0.0, //const
2079 aUSurf2l = 0.0,
2080 aVSurf2f = 0.0,
2081 aVSurf2l = 0.0;
2082
2083 theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
2084 theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
2085
ecc4f148 2086 const gp_Cylinder& aCyl1 = theQuad1.Cylinder(),
2087 aCyl2 = theQuad2.Cylinder();
2088
2089 IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
2090
2091 if (!anInter.IsDone())
2092 {
2093 return Standard_False;
2094 }
2095
2096 IntAna_ResultType aTypInt = anInter.TypeInter();
2097
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.
2101
2102 return Standard_False;
2103 }
2104
2105 theSlin.Clear();
2106 theSPnt.Clear();
7c32c7c4 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);
ecc4f148 2111 const Standard_Real aPeriod = 2.0*M_PI;
2112 const Standard_Real aStepMin = theTol2D,
2113 aStepMax = (aUSurf1l-aUSurf1f)/IntToReal(aNbPoints);
7c32c7c4 2114 const Standard_Integer aNbWLines = 2;
2115
ecc4f148 2116 const stCoeffsValue anEquationCoeffs(aCyl1, aCyl2);
2117
2118 //Boundaries
2119 const Standard_Integer aNbOfBoundaries = 2;
2120 Standard_Real aU1f[aNbOfBoundaries] = {-Precision::Infinite(), -Precision::Infinite()};
2121 Standard_Real aU1l[aNbOfBoundaries] = {Precision::Infinite(), Precision::Infinite()};
2122
2123 if(anEquationCoeffs.mB > 0.0)
2124 {
2125 if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) < -1.0)
2126 {//There is NOT intersection
2127 return Standard_True;
2128 }
2129 else if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
2130 {//U=[0;2*PI]+aFI1
2131 aU1f[0] = anEquationCoeffs.mFI1;
2132 aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
2133 }
2134 else if((1 + anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
2135 (anEquationCoeffs.mB <= 1 - anEquationCoeffs.mC))
2136 {
2137 Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
2138 if(anArg > 1.0)
2139 anArg = 1.0;
2140 if(anArg < -1.0)
2141 anArg = -1.0;
2142
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;
2149 }
2150 else if((1 - anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
2151 (anEquationCoeffs.mB <= 1 + anEquationCoeffs.mC))
2152 {
2153 Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
2154 if(anArg > 1.0)
2155 anArg = 1.0;
2156 if(anArg < -1.0)
2157 anArg = -1.0;
2158
2159 const Standard_Real aDAngle = acos(anArg);
2160 //U=[aDAngle;2*PI-aDAngle]+aFI1
2161
2162 aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
2163 aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
2164 }
2165 else if(anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
2166 {
2167 Standard_Real anArg1 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB,
2168 anArg2 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
2169 if(anArg1 > 1.0)
2170 anArg1 = 1.0;
2171 if(anArg1 < -1.0)
2172 anArg1 = -1.0;
2173
2174 if(anArg2 > 1.0)
2175 anArg2 = 1.0;
2176 if(anArg2 < -1.0)
2177 anArg2 = -1.0;
2178
2179 const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
2180 //(U=[aDAngle1;aDAngle2]+aFI1) ||
2181 //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
2182
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;
2187 }
2188 else
2189 {
2190 Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
2191 }
2192 }
2193 else if(anEquationCoeffs.mB < 0.0)
2194 {
2195 if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) > 1.0)
2196 {//There is NOT intersection
2197 return Standard_True;
2198 }
2199 else if(-anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
2200 {//U=[0;2*PI]+aFI1
2201 aU1f[0] = anEquationCoeffs.mFI1;
2202 aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
2203 }
2204 else if((-anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
2205 ( anEquationCoeffs.mB <= anEquationCoeffs.mC - 1))
2206 {
2207 Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
2208 if(anArg > 1.0)
2209 anArg = 1.0;
2210 if(anArg < -1.0)
2211 anArg = -1.0;
2212
2213 const Standard_Real aDAngle = acos(anArg);
2214 //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
2215
2216 aU1f[0] = anEquationCoeffs.mFI1;
2217 aU1l[0] = aDAngle + anEquationCoeffs.mFI1;
2218 aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
2219 aU1l[1] = aPeriod + anEquationCoeffs.mFI1;
2220 }
2221 else if((anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
2222 (anEquationCoeffs.mB <= -anEquationCoeffs.mB - 1))
2223 {
2224 Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
2225 if(anArg > 1.0)
2226 anArg = 1.0;
2227 if(anArg < -1.0)
2228 anArg = -1.0;
2229
2230 const Standard_Real aDAngle = acos(anArg);
2231 //U=[aDAngle;2*PI-aDAngle]+aFI1
2232
2233 aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
2234 aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
2235 }
2236 else if(-anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
2237 {
2238 Standard_Real anArg1 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB,
2239 anArg2 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
2240 if(anArg1 > 1.0)
2241 anArg1 = 1.0;
2242 if(anArg1 < -1.0)
2243 anArg1 = -1.0;
2244
2245 if(anArg2 > 1.0)
2246 anArg2 = 1.0;
2247 if(anArg2 < -1.0)
2248 anArg2 = -1.0;
2249
2250 const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
2251 //(U=[aDAngle1;aDAngle2]+aFI1) ||
2252 //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
2253
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;
2258 }
2259 else
2260 {
2261 Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
2262 }
2263 }
2264 else
2265 {
2266 Standard_Failure::Raise("Error. Exception. Unhandled case (B-parameter computation)!");
2267 }
2268
2269 for(Standard_Integer i = 0; i < aNbOfBoundaries; i++)
2270 {
2271 if(Precision::IsInfinite(aU1f[i]) && Precision::IsInfinite(aU1l[i]))
2272 continue;
2273
2274 InscribeInterval(aUSurf1f, aUSurf1l, aU1f[i], aU1l[i], theTol2D, aPeriod);
2275 }
2276
2277 if( !Precision::IsInfinite(aU1f[0]) && !Precision::IsInfinite(aU1f[1]) &&
2278 !Precision::IsInfinite(aU1l[0]) && !Precision::IsInfinite(aU1l[1]))
2279 {
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]);
2285
2286 aU1f[1] = -Precision::Infinite();
2287 aU1l[1] = Precision::Infinite();
2288 }
2289 }
2290
2291 //Critical points
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.
02effd35 2297 //[6...9] - in these points an intersection line goes through
ecc4f148 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()};
2310
2311 CriticalPointsComputing(anEquationCoeffs, aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
2312 aPeriod, theTol2D, aNbCritPointsMax, anU1crit);
2313
2314
2315 //Getting Walking-line
2316
b2af2f56 2317 enum WLFStatus
2318 {
2319 WLFStatus_Absent = 0,
2320 WLFStatus_Exist = 1,
2321 WLFStatus_Broken = 2
2322 };
2323
ecc4f148 2324 for(Standard_Integer aCurInterval = 0; aCurInterval < aNbOfBoundaries; aCurInterval++)
2325 {
2326 if(Precision::IsInfinite(aU1f[aCurInterval]) && Precision::IsInfinite(aU1l[aCurInterval]))
2327 continue;
2328
7c32c7c4 2329 Standard_Boolean isAddedIntoWL[aNbWLines];
2330 for(Standard_Integer i = 0; i < aNbWLines; i++)
2331 isAddedIntoWL[i] = Standard_False;
ecc4f148 2332
2333 Standard_Real anUf = aU1f[aCurInterval], anUl = aU1l[aCurInterval];
e8feb725 2334 const Standard_Boolean isDeltaPeriod = IsEqual(anUl-anUf, aPeriod);
ecc4f148 2335
2336 //Inscribe and sort critical points
2337 InscribeAndSortArray(anU1crit, aNbCritPointsMax, anUf, anUl, theTol2D, aPeriod);
2338
2339 while(anUf < anUl)
2340 {
7c32c7c4 2341 Standard_Real aU2[aNbWLines], aV1[aNbWLines], aV2[aNbWLines];
b2af2f56 2342 WLFStatus aWLFindStatus[aNbWLines];
7c32c7c4 2343 Standard_Real aV1Prev[aNbWLines], aV2Prev[aNbWLines];
b5ef9d91 2344 Standard_Real anUexpect[aNbWLines];
7c32c7c4 2345 Standard_Boolean isAddingWLEnabled[aNbWLines];
2346
2347 Handle(IntSurf_LineOn2S) aL2S[aNbWLines];
2348 Handle(IntPatch_WLine) aWLine[aNbWLines];
2349 for(Standard_Integer i = 0; i < aNbWLines; i++)
2350 {
2351 aL2S[i] = new IntSurf_LineOn2S();
2352 aWLine[i] = new IntPatch_WLine(aL2S[i], Standard_False);
b2af2f56 2353 aWLFindStatus[i] = WLFStatus_Absent;
7c32c7c4 2354 isAddingWLEnabled[i] = Standard_True;
2355 aU2[i] = aV1[i] = aV2[i] = 0.0;
2356 aV1Prev[i] = aV2Prev[i] = 0.0;
b5ef9d91 2357 anUexpect[i] = anUf;
7c32c7c4 2358 }
2359
ecc4f148 2360 Standard_Real anU1 = anUf;
2361
2362 Standard_Real aCriticalDelta[aNbCritPointsMax];
2363 for(Standard_Integer i = 0; i < aNbCritPointsMax; i++)
2364 aCriticalDelta[i] = anU1 - anU1crit[i];
2365
ecc4f148 2366 Standard_Boolean isFirst = Standard_True;
2367
2368 while(anU1 <= anUl)
2369 {
2370 for(Standard_Integer i = 0; i < aNbCritPointsMax; i++)
2371 {
2372 if((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0)
2373 {
2374 anU1 = anU1crit[i];
7c32c7c4 2375
2376 for(Standard_Integer i = 0; i < aNbWLines; i++)
b5ef9d91 2377 {
b2af2f56 2378 aWLFindStatus[i] = WLFStatus_Broken;
b5ef9d91 2379 anUexpect[i] = anU1;
2380 }
7c32c7c4 2381
ecc4f148 2382 break;
2383 }
2384 }
2385
b5ef9d91 2386 if(IsEqual(anU1, anUl))
2387 {
2388 for(Standard_Integer i = 0; i < aNbWLines; i++)
2389 {
2390 aWLFindStatus[i] = WLFStatus_Broken;
2391 anUexpect[i] = anU1;
2392
2393 if(isDeltaPeriod)
2394 {
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.
2402
2403 isAddingWLEnabled[i] = (!isAddedIntoWL[i]);
2404 }
2405 else
2406 {
2407 isAddingWLEnabled[i] = ((theTol2D >= (anUexpect[i] - anU1)) ||
2408 (aWLFindStatus[i] == WLFStatus_Absent));
2409 }
2410 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2411 }
2412 else
2413 {
2414 for(Standard_Integer i = 0; i < aNbWLines; i++)
2415 {
2416 isAddingWLEnabled[i] = ((theTol2D >= (anUexpect[i] - anU1)) ||
2417 (aWLFindStatus[i] == WLFStatus_Absent));
2418 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2419 }
ecc4f148 2420
7c32c7c4 2421 for(Standard_Integer i = 0; i < aNbWLines; i++)
2422 {
2423 const Standard_Integer aNbPntsWL = aWLine[i].IsNull() ? 0 :
2424 aWLine[i]->Curve()->NbPoints();
b5ef9d91 2425
2426 CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]);
02effd35 2427
7c32c7c4 2428 InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False);
02effd35 2429
7c32c7c4 2430 if(aNbPntsWL == 0)
2431 {//the line has not contained any points yet
e002f1ce 2432 if(((aUSurf2f + aPeriod - aUSurf2l) <= 2.0*theTol2D) &&
7c32c7c4 2433 ((Abs(aU2[i] - aUSurf2f) < theTol2D) ||
2434 (Abs(aU2[i]-aUSurf2l) < theTol2D)))
02effd35 2435 {
e002f1ce 2436 //In this case aU2[i] can have two values: current aU2[i] or
b5ef9d91 2437 //aU2[i]+aPeriod (aU2[i]-aPeriod). It is necessary to choose
2438 //correct value.
e002f1ce 2439
b5ef9d91 2440 Standard_Boolean isIncreasing = Standard_True;
2441 CylCylMonotonicity(anU1, i, anEquationCoeffs, aPeriod, isIncreasing);
7c32c7c4 2442
e002f1ce 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.
2448 if(isIncreasing)
2449 {
2450 aU2[i] = aUSurf2f;
2451 }
2452 else
2453 {
2454 aU2[i] = aUSurf2l;
7c32c7c4 2455 }
02effd35 2456 }
2457 }
e002f1ce 2458 else
2459 {//aNbPntsWL > 0
7c32c7c4 2460 if(((aUSurf2l - aUSurf2f) >= aPeriod) &&
2461 ((Abs(aU2[i] - aUSurf2f) < theTol2D) ||
2462 (Abs(aU2[i]-aUSurf2l) < theTol2D)))
2463 {//end of the line
2464 Standard_Real aU2prev = 0.0, aV2prev = 0.0;
2465 if(isTheReverse)
2466 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS1(aU2prev, aV2prev);
02effd35 2467 else
7c32c7c4 2468 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS2(aU2prev, aV2prev);
2469
2470 if(2.0*Abs(aU2prev - aU2[i]) > aPeriod)
2471 {
2472 if(aU2prev > aU2[i])
2473 aU2[i] += aPeriod;
2474 else
2475 aU2[i] -= aPeriod;
2476 }
02effd35 2477 }
2478 }
02effd35 2479
e002f1ce 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).
b2af2f56 2484 if(Abs(aU2[i]) <= theTol2D)
2485 aU2[i] = 0.0;
2486 else if(Abs(aU2[i] - aPeriod) <= theTol2D)
2487 aU2[i] = aPeriod;
2488 else if(Abs(aU2[i] - aUSurf2f) <= theTol2D)
2489 aU2[i] = aUSurf2f;
2490 else if(Abs(aU2[i] - aUSurf2l) <= theTol2D)
2491 aU2[i] = aUSurf2l;
2492 }
2493
b5ef9d91 2494 CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs, aV1[i], aV2[i]);
7c32c7c4 2495
2496 if(isFirst)
2497 {
2498 aV1Prev[i] = aV1[i];
2499 aV2Prev[i] = aV2[i];
02effd35 2500 }
7c32c7c4 2501 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
ecc4f148 2502
7c32c7c4 2503 isFirst = Standard_False;
ecc4f148 2504
7c32c7c4 2505 //Looking for points into WLine
2506 Standard_Boolean isBroken = Standard_False;
2507 for(Standard_Integer i = 0; i < aNbWLines; i++)
ecc4f148 2508 {
7c32c7c4 2509 if(!isAddingWLEnabled[i])
ecc4f148 2510 {
b5ef9d91 2511 Standard_Boolean isBoundIntersect = Standard_False;
2512 if( (Abs(aV1[i] - aVSurf1f) <= theTol2D) ||
2513 ((aV1[i]-aVSurf1f)*(aV1Prev[i]-aVSurf1f) < 0.0))
2514 {
2515 isBoundIntersect = Standard_True;
2516 }
2517 else if( (Abs(aV1[i] - aVSurf1l) <= theTol2D) ||
2518 ( (aV1[i]-aVSurf1l)*(aV1Prev[i]-aVSurf1l) < 0.0))
2519 {
2520 isBoundIntersect = Standard_True;
2521 }
2522 else if( (Abs(aV2[i] - aVSurf2f) <= theTol2D) ||
2523 ( (aV2[i]-aVSurf2f)*(aV2Prev[i]-aVSurf2f) < 0.0))
2524 {
2525 isBoundIntersect = Standard_True;
2526 }
2527 else if( (Abs(aV2[i] - aVSurf2l) <= theTol2D) ||
2528 ( (aV2[i]-aVSurf2l)*(aV2Prev[i]-aVSurf2l) < 0.0))
2529 {
2530 isBoundIntersect = Standard_True;
2531 }
2532
7c32c7c4 2533 aV1Prev[i] = aV1[i];
2534 aV2Prev[i] = aV2[i];
e8feb725 2535
b2af2f56 2536 if(aWLFindStatus[i] == WLFStatus_Broken)
7c32c7c4 2537 isBroken = Standard_True;
02effd35 2538
b5ef9d91 2539 if(!isBoundIntersect)
2540 {
2541 continue;
2542 }
2543 else
2544 {
2545 anUexpect[i] = anU1;
2546 }
ecc4f148 2547 }
ecc4f148 2548
b5ef9d91 2549 const Standard_Boolean isInscribe =
2550 ((aUSurf2f-aU2[i]) <= theTol2D) && ((aU2[i]-aUSurf2l) <= theTol2D) &&
7c32c7c4 2551 ((aVSurf1f - aV1[i]) <= theTol2D) && ((aV1[i] - aVSurf1l) <= theTol2D) &&
b5ef9d91 2552 ((aVSurf2f - aV2[i]) <= theTol2D) && ((aV2[i] - aVSurf2l) <= theTol2D);
2553
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()));
2561
2562
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;
2569
2570 if((aWLFindStatus[i] == WLFStatus_Absent))
ecc4f148 2571 {
b5ef9d91 2572 if(((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1-aUSurf1l) < theTol2D))
02effd35 2573 {
b5ef9d91 2574 isForce = Standard_True;
2575 }
2576 }
e8feb725 2577
b5ef9d91 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);
02effd35 2583
b5ef9d91 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));
2589
2590 if((aWLFindStatus[i] == WLFStatus_Exist) && (isFound1 || isFound2) && !isPrevVBound)
2591 {
2592 aWLFindStatus[i] = WLFStatus_Broken; //start a new line
2593 }
2594 else if(isInscribe)
2595 {
2596 if((aWLFindStatus[i] == WLFStatus_Absent) && (isFound1 || isFound2))
2597 {
2598 aWLFindStatus[i] = WLFStatus_Exist;
ecc4f148 2599 }
ecc4f148 2600
b5ef9d91 2601 if(( aWLFindStatus[i] != WLFStatus_Broken) || (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl))
ecc4f148 2602 {
5b9e1842 2603 if(aWLine[i]->NbPnts() > 0)
2604 {
2605 Standard_Real aU2p = 0.0, aV2p = 0.0;
2606 if(isTheReverse)
2607 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
2608 else
2609 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
2610
2611 const Standard_Real aDelta = aU2[i] - aU2p;
2612
2613 if(2*Abs(aDelta) > aPeriod)
2614 {
2615 if(aDelta > 0.0)
2616 {
2617 aU2[i] -= aPeriod;
2618 }
2619 else
2620 {
2621 aU2[i] += aPeriod;
2622 }
2623 }
2624 }
2625
b5ef9d91 2626 if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse, Standard_True,
7c32c7c4 2627 gp_Pnt2d(anU1, aV1[i]), gp_Pnt2d(aU2[i], aV2[i]),
b5ef9d91 2628 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
2629 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod,
2630 aWLine[i]->Curve(), i, theTol3D, theTol2D, isForce))
ecc4f148 2631 {
b2af2f56 2632 if(aWLFindStatus[i] == WLFStatus_Absent)
e8feb725 2633 {
b2af2f56 2634 aWLFindStatus[i] = WLFStatus_Exist;
e8feb725 2635 }
ecc4f148 2636 }
5b9e1842 2637 else if(!isFound1 && !isFound2)
2638 {//We do not add any point while doing this iteration
2639 if(aWLFindStatus[i] == WLFStatus_Exist)
2640 {
2641 aWLFindStatus[i] = WLFStatus_Broken;
2642 }
2643 }
2644 }
2645 }
2646 else
2647 {//We do not add any point while doing this iteration
2648 if(aWLFindStatus[i] == WLFStatus_Exist)
2649 {
2650 aWLFindStatus[i] = WLFStatus_Broken;
ecc4f148 2651 }
2652 }
b5ef9d91 2653
7c32c7c4 2654 aV1Prev[i] = aV1[i];
2655 aV2Prev[i] = aV2[i];
2656
b2af2f56 2657 if(aWLFindStatus[i] == WLFStatus_Broken)
7c32c7c4 2658 isBroken = Standard_True;
2659 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
ecc4f148 2660
7c32c7c4 2661 if(isBroken)
ecc4f148 2662 {//current lines are filled. Go to the next lines
2663 anUf = anU1;
b5ef9d91 2664
5b9e1842 2665 Standard_Boolean isAdded = Standard_True;
2666
b5ef9d91 2667 for(Standard_Integer i = 0; i < aNbWLines; i++)
2668 {
5b9e1842 2669 if(isAddingWLEnabled[i])
2670 {
b5ef9d91 2671 continue;
5b9e1842 2672 }
b5ef9d91 2673
5b9e1842 2674 isAdded = Standard_False;
b5ef9d91 2675
5b9e1842 2676 Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
b5ef9d91 2677
5b9e1842 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);
b5ef9d91 2683
5b9e1842 2684 if(isFound1 || isFound2)
2685 {
2686 isAdded = Standard_True;
2687 }
b5ef9d91 2688
5b9e1842 2689 if(aWLine[i]->NbPnts() > 0)
b5ef9d91 2690 {
5b9e1842 2691 Standard_Real aU2p = 0.0, aV2p = 0.0;
2692 if(isTheReverse)
2693 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
2694 else
2695 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
2696
2697 const Standard_Real aDelta = aU2[i] - aU2p;
2698
2699 if(2*Abs(aDelta) > aPeriod)
2700 {
2701 if(aDelta > 0.0)
2702 {
2703 aU2[i] -= aPeriod;
2704 }
2705 else
2706 {
2707 aU2[i] += aPeriod;
2708 }
2709 }
2710 }
b5ef9d91 2711
5b9e1842 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))
2718 {
2719 isAdded = Standard_True;
b5ef9d91 2720 }
5b9e1842 2721 }
b5ef9d91 2722
5b9e1842 2723 if(!isAdded)
2724 {
2725 Standard_Real anUmaxAdded = RealFirst();
2726 for(Standard_Integer i = 0; i < aNbWLines; i++)
2727 {
2728 Standard_Real aU1c = 0.0, aV1c = 0.0;
2729 if(isTheReverse)
2730 aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS2(aU1c, aV1c);
2731 else
2732 aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS1(aU1c, aV1c);
b5ef9d91 2733
5b9e1842 2734 anUmaxAdded = Max(anUmaxAdded, aU1c);
b5ef9d91 2735 }
2736
5b9e1842 2737 for(Standard_Integer i = 0; i < aNbWLines; i++)
2738 {
2739 if(isAddingWLEnabled[i])
2740 {
2741 continue;
2742 }
b5ef9d91 2743
5b9e1842 2744 CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs, aU2[i], aV1[i], aV2[i]);
b5ef9d91 2745
5b9e1842 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);
2752 }
b5ef9d91 2753 }
2754
ecc4f148 2755 break;
2756 }
2757
7c32c7c4 2758 //Step computing
ecc4f148 2759
7c32c7c4 2760 {
b5ef9d91 2761 const Standard_Real aDeltaV1 = (aVSurf1l - aVSurf1f)/IntToReal(aNbPoints),
2762 aDeltaV2 = (aVSurf2l - aVSurf2f)/IntToReal(aNbPoints);
2763
2764 math_Matrix aMatr(1, 3, 1, 5);
ecc4f148 2765
b5ef9d91 2766 Standard_Real aMinUexp = RealLast();
2767
2768 for(Standard_Integer i = 0; i < aNbWLines; i++)
2769 {
2770 if(theTol2D < (anUexpect[i] - anU1))
2771 {
2772 continue;
2773 }
ecc4f148 2774
b5ef9d91 2775 if(aWLFindStatus[i] == WLFStatus_Absent)
2776 {
2777 anUexpect[i] += aStepMax;
2778 aMinUexp = Min(aMinUexp, anUexpect[i]);
2779 continue;
2780 }
7c32c7c4 2781
b5ef9d91 2782 Standard_Real aStepTmp = aStepMax;
ecc4f148 2783
b5ef9d91 2784 const Standard_Real aSinU1 = sin(anU1),
2785 aCosU1 = cos(anU1),
2786 aSinU2 = sin(aU2[i]),
2787 aCosU2 = cos(aU2[i]);
ecc4f148 2788
b5ef9d91 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);
ecc4f148 2796
b5ef9d91 2797 if(!StepComputing(aMatr, aV1[i], aV2[i], aDeltaV1, aDeltaV2, aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aStepTmp))
2798 {
2799 //To avoid cycling-up
2800 anUexpect[i] += aStepMax;
2801 aMinUexp = Min(aMinUexp, anUexpect[i]);
ecc4f148 2802
b5ef9d91 2803 continue;
2804 }
2805
2806 if(aStepTmp < aStepMin)
2807 aStepTmp = aStepMin;
ecc4f148 2808
b5ef9d91 2809 if(aStepTmp > aStepMax)
2810 aStepTmp = aStepMax;
2811
2812 anUexpect[i] = anU1 + aStepTmp;
2813 aMinUexp = Min(aMinUexp, anUexpect[i]);
2814 }
ecc4f148 2815
b5ef9d91 2816 anU1 = aMinUexp;
2817 }
ecc4f148 2818
b5ef9d91 2819 if(Precision::PConfusion() >= (anUl - anU1))
ecc4f148 2820 anU1 = anUl;
2821
2822 anUf = anU1;
2823
7c32c7c4 2824 for(Standard_Integer i = 0; i < aNbWLines; i++)
baf72cd2 2825 {
7c32c7c4 2826 if(aWLine[i]->NbPnts() != 1)
2827 isAddedIntoWL[i] = Standard_False;
b5ef9d91 2828
2829 if(anU1 == anUl)
2830 {//strictly equal. Tolerance is considered above.
2831 anUexpect[i] = anUl;
2832 }
baf72cd2 2833 }
ecc4f148 2834 }
ecc4f148 2835
7c32c7c4 2836 for(Standard_Integer i = 0; i < aNbWLines; i++)
ecc4f148 2837 {
7c32c7c4 2838 if((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
baf72cd2 2839 {
7c32c7c4 2840 isTheEmpty = Standard_False;
2841 Standard_Real u1, v1, u2, v2;
2842 aWLine[i]->Point(1).Parameters(u1, v1, u2, v2);
2843 IntPatch_Point aP;
2844 aP.SetParameter(u1);
2845 aP.SetParameters(u1, v1, u2, v2);
2846 aP.SetTolerance(theTol3D);
2847 aP.SetValue(aWLine[i]->Point(1).Value());
2848
2849 theSPnt.Append(aP);
baf72cd2 2850 }
7c32c7c4 2851 else if(aWLine[i]->NbPnts() > 1)
baf72cd2 2852 {
7c32c7c4 2853 Standard_Boolean isGood = Standard_True;
ecc4f148 2854
7c32c7c4 2855 if(aWLine[i]->NbPnts() == 2)
2856 {
2857 const IntSurf_PntOn2S& aPf = aWLine[i]->Point(1);
2858 const IntSurf_PntOn2S& aPl = aWLine[i]->Point(2);
2859
2860 if(aPf.IsSame(aPl, Precision::Confusion()))
2861 isGood = Standard_False;
2862 }
ecc4f148 2863
7c32c7c4 2864 if(isGood)
2865 {
2866 isTheEmpty = Standard_False;
2867 isAddedIntoWL[i] = Standard_True;
2868 SeekAdditionalPoints( theQuad1, theQuad2, aWLine[i]->Curve(),
b5ef9d91 2869 anEquationCoeffs, i, aNbPoints, 1,
2870 aWLine[i]->NbPnts(), aUSurf2f, aUSurf2l,
2871 theTol2D, aPeriod, isTheReverse);
7c32c7c4 2872
2873 aWLine[i]->ComputeVertexParameters(theTol3D);
2874 theSlin.Append(aWLine[i]);
2875 }
2876 }
2877 else
2878 {
2879 isAddedIntoWL[i] = Standard_False;
baf72cd2 2880 }
b5ef9d91 2881
2882#ifdef OCCT_DEBUG
2883 //aWLine[i]->Dump();
2884#endif
2885 }
2886 }
2887 }
2888
2889
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++)
2893 {
2894 for(Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); aNbLin++)
2895 {
2896 const Handle(IntPatch_WLine)& aWLine1 =
2897 Handle(IntPatch_WLine)::DownCast(theSlin.Value(aNbLin));
2898
2899 const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1);
2900 const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aWLine1->NbPnts());
2901
2902 const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNbPnt).PntOn2S();
2903 if( aPntCur.IsSame(aPntFWL1, Precision::Confusion()) ||
2904 aPntCur.IsSame(aPntLWL1, Precision::Confusion()))
2905 {
2906 theSPnt.Remove(aNbPnt);
2907 aNbPnt--;
2908 break;
ecc4f148 2909 }
2910 }
2911 }
2912
b5ef9d91 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++)
2916 {
2917 const IntPatch_Point& aPnt2S = theSPnt.Value(aNbPnt);
2918
2919 Standard_Real u1, v1, u2, v2;
2920 aPnt2S.Parameters(u1, v1, u2, v2);
2921
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());
2925
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;
2929 if(isTheReverse)
2930 {
2931 anUf = Max(u2 - aStepMax, aUSurf1f);
2932 anUl = u2;
2933 aCurU2 = u1;
2934 }
2935 else
2936 {
2937 anUf = Max(u1 - aStepMax, aUSurf1f);
2938 anUl = u1;
2939 aCurU2 = u2;
2940 }
2941 Standard_Real aDelta = RealLast();
2942 for (Standard_Integer i = 0; i < aNbWLines; i++)
2943 {
2944 Standard_Real anU2t = 0.0;
2945 if(!CylCylComputeParameters(anUl, i, anEquationCoeffs, anU2t))
2946 continue;
2947
2948 const Standard_Real aDU = Abs(anU2t - aCurU2);
2949 if(aDU < aDelta)
2950 {
2951 aDelta = aDU;
2952 anIndex = i;
2953 }
2954 }
2955
2956 //Try to fill aWLine by additional points
2957 while(anUl - anUf > RealSmall())
2958 {
2959 Standard_Real anU2 = 0.0, anV1 = 0.0, anV2 = 0.0;
2960 Standard_Boolean isDone =
2961 CylCylComputeParameters(anUf, anIndex, anEquationCoeffs,
2962 anU2, anV1, anV2);
2963 anUf += aDU;
5b9e1842 2964
b5ef9d91 2965 if(!isDone)
2966 {
2967 continue;
2968 }
2969
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))
2976 {
2977 break;
2978 }
2979 }
2980
2981 if(aWLine->NbPnts() > 1)
2982 {
2983 SeekAdditionalPoints( theQuad1, theQuad2, aWLine->Curve(),
2984 anEquationCoeffs, anIndex, aNbMinPoints,
2985 1, aWLine->NbPnts(), aUSurf2f, aUSurf2l,
2986 theTol2D, aPeriod, isTheReverse);
2987
2988 aWLine->ComputeVertexParameters(theTol3D);
2989 theSlin.Append(aWLine);
5b9e1842 2990
b5ef9d91 2991 theSPnt.Remove(aNbPnt);
2992 aNbPnt--;
2993 }
2994 }
2995
ecc4f148 2996 return Standard_True;
2997}
7fd59977 2998
2999//=======================================================================
3000//function : IntCySp
3001//purpose :
3002//=======================================================================
3003Standard_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)
3011
3012{
3013 Standard_Integer i;
3014
3015 IntSurf_TypeTrans trans1,trans2;
3016 IntAna_ResultType typint;
3017 IntPatch_Point ptsol;
3018 gp_Circ cirsol;
3019
3020 gp_Cylinder Cy;
3021 gp_Sphere Sp;
3022
3023 if (!Reversed) {
3024 Cy = Quad1.Cylinder();
3025 Sp = Quad2.Sphere();
3026 }
3027 else {
3028 Cy = Quad2.Cylinder();
3029 Sp = Quad1.Sphere();
3030 }
3031 IntAna_QuadQuadGeo inter(Cy,Sp,Tol);
3032
3033 if (!inter.IsDone()) {return Standard_False;}
3034
3035 typint = inter.TypeInter();
3036 Standard_Integer NbSol = inter.NbSolutions();
3037 Empty = Standard_False;
3038
3039 switch (typint) {
3040
3041 case IntAna_Empty :
3042 {
3043 Empty = Standard_True;
3044 }
3045 break;
3046
3047 case IntAna_Point :
3048 {
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);
3055 spnt.Append(ptsol);
3056 }
3057 break;
3058
3059 case IntAna_Circle:
3060 {
3061 cirsol = inter.Circle(1);
3062 gp_Vec Tgt;
3063 gp_Pnt ptref;
3064 ElCLib::D1(0.,cirsol,ptref,Tgt);
3065
3066 if (NbSol == 1) {
3067 gp_Vec TestCurvature(ptref,Sp.Location());
3068 gp_Vec Normsp,Normcyl;
3069 if (!Reversed) {
3070 Normcyl = Quad1.Normale(ptref);
3071 Normsp = Quad2.Normale(ptref);
3072 }
3073 else {
3074 Normcyl = Quad2.Normale(ptref);
3075 Normsp = Quad1.Normale(ptref);
3076 }
3077
3078 IntSurf_Situation situcyl;
3079 IntSurf_Situation situsp;
3080
3081 if (Normcyl.Dot(TestCurvature) > 0.) {
3082 situsp = IntSurf_Outside;
3083 if (Normsp.Dot(Normcyl) > 0.) {
3084 situcyl = IntSurf_Inside;
3085 }
3086 else {
3087 situcyl = IntSurf_Outside;
3088 }
3089 }
3090 else {
3091 situsp = IntSurf_Inside;
3092 if (Normsp.Dot(Normcyl) > 0.) {
3093 situcyl = IntSurf_Outside;
3094 }
3095 else {
3096 situcyl = IntSurf_Inside;
3097 }
3098 }
3099 Handle(IntPatch_GLine) glig;
3100 if (!Reversed) {
3101 glig = new IntPatch_GLine(cirsol, Standard_True, situcyl, situsp);
3102 }
3103 else {
3104 glig = new IntPatch_GLine(cirsol, Standard_True, situsp, situcyl);
3105 }
3106 slin.Append(glig);
3107 }
3108 else {
3109 if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
3110 trans1 = IntSurf_Out;
3111 trans2 = IntSurf_In;
3112 }
3113 else {
3114 trans1 = IntSurf_In;
3115 trans2 = IntSurf_Out;
3116 }
3117 Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3118 slin.Append(glig);
3119
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;
3126 }
3127 else if(qwe<-0.0000001) {
3128 trans1 = IntSurf_In;
3129 trans2 = IntSurf_Out;
3130 }
3131 else {
3132 trans1=trans2=IntSurf_Undecided;
3133 }
3134 glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3135 slin.Append(glig);
3136 }
3137 }
3138 break;
3139
3140 case IntAna_NoGeometricSolution:
3141 {
3142 gp_Pnt psol;
3143 Standard_Real U1,V1,U2,V2;
3144 IntAna_IntQuadQuad anaint(Cy,Sp,Tol);
3145 if (!anaint.IsDone()) {
3146 return Standard_False;
3147 }
3148
3149 if (anaint.NbPnt()==0 && anaint.NbCurve()==0) {
3150 Empty = Standard_True;
3151 }
3152 else {
3153
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);
3161 spnt.Append(ptsol);
3162 }
3163
3164 gp_Pnt ptvalid,ptf,ptl;
3165 gp_Vec tgvalid;
3166 Standard_Real first,last,para;
3167 IntAna_Curve curvsol;
3168 Standard_Boolean tgfound;
3169 Standard_Integer kount;
3170
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);
3177
3178 para = last;
3179 kount = 1;
3180 tgfound = Standard_False;
3181
3182 while (!tgfound) {
3183 para = (1.123*first + para)/2.123;
3184 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
3185 if (!tgfound) {
3186 kount ++;
3187 tgfound = kount > 5;
3188 }
3189 }
3190 Handle(IntPatch_ALine) alig;
3191 if (kount <= 5) {
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;
3197 }
3198 else if(qwe<-0.00000001) {
3199 trans1 = IntSurf_In;
3200 trans2 = IntSurf_Out;
3201 }
3202 else {
3203 trans1=trans2=IntSurf_Undecided;
3204 }
3205 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
3206 }
3207 else {
3208 alig = new IntPatch_ALine(curvsol,Standard_False);
3209 }
3210 Standard_Boolean TempFalse1a = Standard_False;
3211 Standard_Boolean TempFalse2a = Standard_False;
3212
3213 //-- ptf et ptl : points debut et fin de alig
3214
3215 ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1a,ptf,first,
3216 TempFalse2a,ptl,last,Multpoint,Tol);
3217 slin.Append(alig);
3218 } //-- boucle sur les lignes
3219 } //-- solution avec au moins une lihne
3220 }
3221 break;
3222
3223 default:
3224 {
3225 return Standard_False;
3226 }
3227 }
3228 return Standard_True;
3229}
3230//=======================================================================
3231//function : IntCyCo
3232//purpose :
3233//=======================================================================
3234Standard_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)
3242
3243{
3244 IntPatch_Point ptsol;
3245
3246 Standard_Integer i;
3247
3248 IntSurf_TypeTrans trans1,trans2;
3249 IntAna_ResultType typint;
3250 gp_Circ cirsol;
3251
3252 gp_Cylinder Cy;
3253 gp_Cone Co;
3254
3255 if (!Reversed) {
3256 Cy = Quad1.Cylinder();
3257 Co = Quad2.Cone();
3258 }
3259 else {
3260 Cy = Quad2.Cylinder();
3261 Co = Quad1.Cone();
3262 }
3263 IntAna_QuadQuadGeo inter(Cy,Co,Tol);
3264
3265 if (!inter.IsDone()) {return Standard_False;}
3266
3267 typint = inter.TypeInter();
3268 Standard_Integer NbSol = inter.NbSolutions();
3269 Empty = Standard_False;
3270
3271 switch (typint) {
3272
3273 case IntAna_Empty : {
3274 Empty = Standard_True;
3275 }
3276 break;
3277
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);
3285 spnt.Append(ptsol);
3286 }
3287 break;
3288
3289 case IntAna_Circle: {
3290 gp_Vec Tgt;
3291 gp_Pnt ptref;
3292 Standard_Integer j;
3293 Standard_Real qwe;
3294 //
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;
3302 }
3303 else if(qwe<-0.00000001) {
3304 trans1 = IntSurf_In;
3305 trans2 = IntSurf_Out;
3306 }
3307 else {
3308 trans1=trans2=IntSurf_Undecided;
3309 }
3310 Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3311 slin.Append(glig);
3312 }
3313 }
3314 break;
3315
3316 case IntAna_NoGeometricSolution: {
3317 gp_Pnt psol;
3318 Standard_Real U1,V1,U2,V2;
3319 IntAna_IntQuadQuad anaint(Cy,Co,Tol);
3320 if (!anaint.IsDone()) {
3321 return Standard_False;
3322 }
3323
3324 if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) {
3325 Empty = Standard_True;
3326 }
3327 else {
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);
3335 spnt.Append(ptsol);
3336 }
3337
3338 gp_Pnt ptvalid, ptf, ptl;
3339 gp_Vec tgvalid;
3340 //
3341 Standard_Real first,last,para;
3342 Standard_Boolean tgfound,firstp,lastp,kept;
3343 Standard_Integer kount;
3344 //
3345 //
3346 //IntAna_Curve curvsol;
3347 IntAna_Curve aC;
3348 IntAna_ListOfCurve aLC;
3349 IntAna_ListIteratorOfListOfCurve aIt;
7fd59977 3350
3351 //
3352 NbSol = anaint.NbCurve();
3353 for (i = 1; i <= NbSol; ++i) {
3354 kept = Standard_False;
3355 //curvsol = anaint.Curve(i);
3356 aC=anaint.Curve(i);
3357 aLC.Clear();
96a95605 3358 ExploreCurve(Cy, Co, aC, 10.*Tol, aLC);
7fd59977 3359 //
3360 aIt.Initialize(aLC);
3361 for (; aIt.More(); aIt.Next()) {
3362 IntAna_Curve& curvsol=aIt.Value();
3363 //
3364 curvsol.Domain(first, last);
3365 firstp = !curvsol.IsFirstOpen();
3366 lastp = !curvsol.IsLastOpen();
3367 if (firstp) {
3368 ptf = curvsol.Value(first);
3369 }
3370 if (lastp) {
3371 ptl = curvsol.Value(last);
3372 }
3373 para = last;
3374 kount = 1;
3375 tgfound = Standard_False;
3376
3377 while (!tgfound) {
3378 para = (1.123*first + para)/2.123;
3379 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
3380 if (!tgfound) {
3381 kount ++;
3382 tgfound = kount > 5;
3383 }
3384 }
3385 Handle(IntPatch_ALine) alig;
3386 if (kount <= 5) {
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;
3392 }
3393 else if(qwe<-0.00000001) {
3394 trans1 = IntSurf_In;
3395 trans2 = IntSurf_Out;
3396 }
3397 else {
3398 trans1=trans2=IntSurf_Undecided;
3399 }
3400 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
3401 kept = Standard_True;
3402 }
3403 else {
3404 ptvalid = curvsol.Value(para);
3405 alig = new IntPatch_ALine(curvsol,Standard_False);
3406 kept = Standard_True;
3407 //-- cout << "Transition indeterminee" << endl;
3408 }
3409 if (kept) {
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);
3414 slin.Append(alig);
3415 }
3416 } // for (; aIt.More(); aIt.Next())
3417 } // for (i = 1; i <= NbSol; ++i)
3418 }
3419 }
3420 break;
3421
3422 default:
3423 return Standard_False;
3424
3425 } // switch (typint)
3426
3427 return Standard_True;
3428}
3429//=======================================================================
3430//function : ExploreCurve
3431//purpose :
3432//=======================================================================
3433Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy,
3434 const gp_Cone& aCo,
3435 IntAna_Curve& aC,
3436 const Standard_Real aTol,
3437 IntAna_ListOfCurve& aLC)
3438
3439{
3440 Standard_Boolean bFind=Standard_False;
3441 Standard_Real aTheta, aT1, aT2, aDst;
3442 gp_Pnt aPapx, aPx;
3443 //
3444 //aC.Dump();
3445 //
3446 aLC.Clear();
3447 aLC.Append(aC);
3448 //
3449 aPapx=aCo.Apex();
3450 //
3451 aC.Domain(aT1, aT2);
3452 //
3453 aPx=aC.Value(aT1);
3454 aDst=aPx.Distance(aPapx);
3455 if (aDst<aTol) {
3456 return bFind;
3457 }
3458 aPx=aC.Value(aT2);
3459 aDst=aPx.Distance(aPapx);
3460 if (aDst<aTol) {
3461 return bFind;
3462 }
3463 //
3464 bFind=aC.FindParameter(aPapx, aTheta);
3465 if (!bFind){
3466 return bFind;
3467 }
3468 //
3469 aPx=aC.Value(aTheta);
3470 aDst=aPx.Distance(aPapx);
3471 if (aDst>aTol) {
3472 return !bFind;
3473 }
3474 //
3475 // need to be splitted at aTheta
3476 IntAna_Curve aC1, aC2;
3477 //
3478 aC1=aC;
3479 aC1.SetDomain(aT1, aTheta);
3480 aC2=aC;
3481 aC2.SetDomain(aTheta, aT2);
3482 //
3483 aLC.Clear();
3484 aLC.Append(aC1);
3485 aLC.Append(aC2);
3486 //
3487 return bFind;
3488}
3489//=======================================================================
3490//function : IsToReverse
3491//purpose :
3492//=======================================================================
3493Standard_Boolean IsToReverse(const gp_Cylinder& Cy1,
3494 const gp_Cylinder& Cy2,
3495 const Standard_Real Tol)
3496{
3497 Standard_Boolean bRet;
3498 Standard_Real aR1,aR2, dR, aSc1, aSc2;
3499 //
3500 bRet=Standard_False;
3501 //
3502 aR1=Cy1.Radius();
3503 aR2=Cy2.Radius();
3504 dR=aR1-aR2;
3505 if (dR<0.) {
3506 dR=-dR;
3507 }
3508 if (dR>Tol) {
3509 bRet=aR1>aR2;
3510 }
3511 else {
3512 gp_Dir aDZ(0.,0.,1.);
3513 //
3514 const gp_Dir& aDir1=Cy1.Axis().Direction();
3515 aSc1=aDir1*aDZ;
3516 if (aSc1<0.) {
3517 aSc1=-aSc1;
3518 }
3519 //
3520 const gp_Dir& aDir2=Cy2.Axis().Direction();
3521 aSc2=aDir2*aDZ;
3522 if (aSc2<0.) {
3523 aSc2=-aSc2;
3524 }
3525 //
3526 if (aSc2<aSc1) {
3527 bRet=!bRet;
3528 }
3529 }//if (dR<Tol) {
3530 return bRet;
3531}
Open CASCADE Technology repository