projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0024915: Wrong intersection curves between two cylinders
[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
fa0291ff 21//
7fd59977 22static
23 Standard_Boolean ExploreCurve(const gp_Cylinder& aCy,
24 const gp_Cone& aCo,
25 IntAna_Curve& aC,
26 const Standard_Real aTol,
27 IntAna_ListOfCurve& aLC);
28static
29 Standard_Boolean IsToReverse(const gp_Cylinder& Cy1,
30 const gp_Cylinder& Cy2,
31 const Standard_Real Tol);
32
33//=======================================================================
34//function : ProcessBounds
35//purpose :
36//=======================================================================
37void ProcessBounds(const Handle(IntPatch_ALine)& alig, //-- ligne courante
38 const IntPatch_SequenceOfLine& slin,
39 const IntSurf_Quadric& Quad1,
40 const IntSurf_Quadric& Quad2,
41 Standard_Boolean& procf,
42 const gp_Pnt& ptf, //-- Debut Ligne Courante
43 const Standard_Real first, //-- Paramf
44 Standard_Boolean& procl,
45 const gp_Pnt& ptl, //-- Fin Ligne courante
46 const Standard_Real last, //-- Paraml
47 Standard_Boolean& Multpoint,
48 const Standard_Real Tol)
49{
50 Standard_Integer j,k;
51 Standard_Real U1,V1,U2,V2;
52 IntPatch_Point ptsol;
53 Standard_Real d;
54
55 if (procf && procl) {
56 j = slin.Length() + 1;
57 }
58 else {
59 j = 1;
60 }
61
62
63 //-- On prend les lignes deja enregistrees
64
65 while (j <= slin.Length()) {
66 if(slin.Value(j)->ArcType() == IntPatch_Analytic) {
67 const Handle(IntPatch_ALine)& aligold = *((Handle(IntPatch_ALine)*)&slin.Value(j));
68 k = 1;
69
70 //-- On prend les vertex des lignes deja enregistrees
71
72 while (k <= aligold->NbVertex()) {
73 ptsol = aligold->Vertex(k);
74 if (!procf) {
75 d=ptf.Distance(ptsol.Value());
76 if (d <= Tol) {
77 if (!ptsol.IsMultiple()) {
78 //-- le point ptsol (de aligold) est declare multiple sur aligold
79 Multpoint = Standard_True;
80 ptsol.SetMultiple(Standard_True);
81 aligold->Replace(k,ptsol);
82 }
83 ptsol.SetParameter(first);
84 alig->AddVertex(ptsol);
85 alig->SetFirstPoint(alig->NbVertex());
86 procf = Standard_True;
87
88 //-- On restore le point avec son parametre sur aligold
89 ptsol = aligold->Vertex(k);
90
91 }
92 }
93 if (!procl) {
94 if (ptl.Distance(ptsol.Value()) <= Tol) {
95 if (!ptsol.IsMultiple()) {
96 Multpoint = Standard_True;
97 ptsol.SetMultiple(Standard_True);
98 aligold->Replace(k,ptsol);
99 }
100 ptsol.SetParameter(last);
101 alig->AddVertex(ptsol);
102 alig->SetLastPoint(alig->NbVertex());
103 procl = Standard_True;
104
105 //-- On restore le point avec son parametre sur aligold
106 ptsol = aligold->Vertex(k);
107
108 }
109 }
110 if (procf && procl) {
111 k = aligold->NbVertex()+1;
112 }
113 else {
114 k = k+1;
115 }
116 }
117 if (procf && procl) {
118 j = slin.Length()+1;
119 }
120 else {
121 j = j+1;
122 }
123 }
124 }
125 if (!procf && !procl) {
126 Quad1.Parameters(ptf,U1,V1);
127 Quad2.Parameters(ptf,U2,V2);
128 ptsol.SetValue(ptf,Tol,Standard_False);
129 ptsol.SetParameters(U1,V1,U2,V2);
130 ptsol.SetParameter(first);
131 if (ptf.Distance(ptl) <= Tol) {
132 ptsol.SetMultiple(Standard_True); // a voir
133 Multpoint = Standard_True; // a voir de meme
134 alig->AddVertex(ptsol);
135 alig->SetFirstPoint(alig->NbVertex());
136
137 ptsol.SetParameter(last);
138 alig->AddVertex(ptsol);
139 alig->SetLastPoint(alig->NbVertex());
140 }
141 else {
142 alig->AddVertex(ptsol);
143 alig->SetFirstPoint(alig->NbVertex());
144 Quad1.Parameters(ptl,U1,V1);
145 Quad2.Parameters(ptl,U2,V2);
146 ptsol.SetValue(ptl,Tol,Standard_False);
147 ptsol.SetParameters(U1,V1,U2,V2);
148 ptsol.SetParameter(last);
149 alig->AddVertex(ptsol);
150 alig->SetLastPoint(alig->NbVertex());
151 }
152 }
153 else if (!procf) {
154 Quad1.Parameters(ptf,U1,V1);
155 Quad2.Parameters(ptf,U2,V2);
156 ptsol.SetValue(ptf,Tol,Standard_False);
157 ptsol.SetParameters(U1,V1,U2,V2);
158 ptsol.SetParameter(first);
159 alig->AddVertex(ptsol);
160 alig->SetFirstPoint(alig->NbVertex());
161 }
162 else if (!procl) {
163 Quad1.Parameters(ptl,U1,V1);
164 Quad2.Parameters(ptl,U2,V2);
165 ptsol.SetValue(ptl,Tol,Standard_False);
166 ptsol.SetParameters(U1,V1,U2,V2);
167 ptsol.SetParameter(last);
168 alig->AddVertex(ptsol);
169 alig->SetLastPoint(alig->NbVertex());
170 }
171}
172//=======================================================================
173//function : IntCyCy
174//purpose :
175//=======================================================================
176Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
177 const IntSurf_Quadric& Quad2,
178 const Standard_Real Tol,
179 Standard_Boolean& Empty,
180 Standard_Boolean& Same,
181 Standard_Boolean& Multpoint,
182 IntPatch_SequenceOfLine& slin,
183 IntPatch_SequenceOfPoint& spnt)
184
185{
186 IntPatch_Point ptsol;
187
188 Standard_Integer i;
189
190 IntSurf_TypeTrans trans1,trans2;
191 IntAna_ResultType typint;
192
193 gp_Elips elipsol;
194 gp_Lin linsol;
195
196 gp_Cylinder Cy1(Quad1.Cylinder());
197 gp_Cylinder Cy2(Quad2.Cylinder());
198
199 IntAna_QuadQuadGeo inter(Cy1,Cy2,Tol);
200
ecc4f148 201 if (!inter.IsDone())
202 {
203 return Standard_False;
204 }
7fd59977 205
206 typint = inter.TypeInter();
207 Standard_Integer NbSol = inter.NbSolutions();
208 Empty = Standard_False;
209 Same = Standard_False;
210
ecc4f148 211 switch (typint)
212 {
213 case IntAna_Empty:
7fd59977 214 {
215 Empty = Standard_True;
216 }
217 break;
218
219 case IntAna_Same:
220 {
221 Same = Standard_True;
222 }
223 break;
ecc4f148 224
225 case IntAna_Point:
7fd59977 226 {
227 gp_Pnt psol(inter.Point(1));
228 Standard_Real U1,V1,U2,V2;
229 Quad1.Parameters(psol,U1,V1);
230 Quad2.Parameters(psol,U2,V2);
231 ptsol.SetValue(psol,Tol,Standard_True);
232 ptsol.SetParameters(U1,V1,U2,V2);
233 spnt.Append(ptsol);
234 }
235 break;
236
237 case IntAna_Line:
238 {
239 gp_Pnt ptref;
ecc4f148 240 if (NbSol == 1)
241 { // Cylinders are tangent to each other by line
242 linsol = inter.Line(1);
243 ptref = linsol.Location();
244 gp_Dir crb1(gp_Vec(ptref,Cy1.Location()));
245 gp_Dir crb2(gp_Vec(ptref,Cy2.Location()));
246 gp_Vec norm1(Quad1.Normale(ptref));
247 gp_Vec norm2(Quad2.Normale(ptref));
248 IntSurf_Situation situcyl1;
249 IntSurf_Situation situcyl2;
250
251 if (crb1.Dot(crb2) < 0.)
252 { // centre de courbures "opposes"
253 if (norm1.Dot(crb1) > 0.)
254 {
255 situcyl2 = IntSurf_Inside;
256 }
257 else
258 {
259 situcyl2 = IntSurf_Outside;
260 }
261
262 if (norm2.Dot(crb2) > 0.)
263 {
264 situcyl1 = IntSurf_Inside;
265 }
266 else
267 {
268 situcyl1 = IntSurf_Outside;
269 }
270 }
271 else
272 {
273 if (Cy1.Radius() < Cy2.Radius())
274 {
275 if (norm1.Dot(crb1) > 0.)
276 {
277 situcyl2 = IntSurf_Inside;
278 }
279 else
280 {
281 situcyl2 = IntSurf_Outside;
282 }
283
284 if (norm2.Dot(crb2) > 0.)
285 {
286 situcyl1 = IntSurf_Outside;
287 }
288 else
289 {
290 situcyl1 = IntSurf_Inside;
291 }
292 }
293 else
294 {
295 if (norm1.Dot(crb1) > 0.)
296 {
297 situcyl2 = IntSurf_Outside;
298 }
299 else
300 {
301 situcyl2 = IntSurf_Inside;
302 }
303
304 if (norm2.Dot(crb2) > 0.)
305 {
306 situcyl1 = IntSurf_Inside;
307 }
308 else
309 {
310 situcyl1 = IntSurf_Outside;
311 }
312 }
313 }
314
315 Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_True, situcyl1, situcyl2);
316 slin.Append(glig);
7fd59977 317 }
ecc4f148 318 else
319 {
320 for (i=1; i <= NbSol; i++)
321 {
322 linsol = inter.Line(i);
323 ptref = linsol.Location();
324 gp_Vec lsd = linsol.Direction();
325 Standard_Real qwe = lsd.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
326 if (qwe >0.00000001)
327 {
328 trans1 = IntSurf_Out;
329 trans2 = IntSurf_In;
330 }
331 else if (qwe <-0.00000001)
332 {
333 trans1 = IntSurf_In;
334 trans2 = IntSurf_Out;
335 }
336 else
337 {
338 trans1=trans2=IntSurf_Undecided;
339 }
340
341 Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
342 slin.Append(glig);
343 }
7fd59977 344 }
345 }
346 break;
347
348 case IntAna_Ellipse:
349 {
350 gp_Vec Tgt;
351 gp_Pnt ptref;
fa0291ff 352 IntPatch_Point pmult1, pmult2;
353
7fd59977 354 elipsol = inter.Ellipse(1);
fa0291ff 355
356 gp_Pnt pttang1(ElCLib::Value(0.5*M_PI, elipsol));
357 gp_Pnt pttang2(ElCLib::Value(1.5*M_PI, elipsol));
7fd59977 358
359 Multpoint = Standard_True;
360 pmult1.SetValue(pttang1,Tol,Standard_True);
361 pmult2.SetValue(pttang2,Tol,Standard_True);
362 pmult1.SetMultiple(Standard_True);
363 pmult2.SetMultiple(Standard_True);
364
365 Standard_Real oU1,oV1,oU2,oV2;
366 Quad1.Parameters(pttang1,oU1,oV1);
367 Quad2.Parameters(pttang1,oU2,oV2);
368 pmult1.SetParameters(oU1,oV1,oU2,oV2);
369
370 Quad1.Parameters(pttang2,oU1,oV1);
371 Quad2.Parameters(pttang2,oU2,oV2);
372 pmult2.SetParameters(oU1,oV1,oU2,oV2);
fa0291ff 373
7fd59977 374 // on traite la premiere ellipse
375
376 //-- Calcul de la Transition de la ligne
377 ElCLib::D1(0.,elipsol,ptref,Tgt);
378 Standard_Real qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
ecc4f148 379 if (qwe>0.00000001)
380 {
381 trans1 = IntSurf_Out;
382 trans2 = IntSurf_In;
7fd59977 383 }
ecc4f148 384 else if (qwe<-0.00000001)
385 {
386 trans1 = IntSurf_In;
387 trans2 = IntSurf_Out;
7fd59977 388 }
ecc4f148 389 else
390 {
391 trans1=trans2=IntSurf_Undecided;
7fd59977 392 }
ecc4f148 393
7fd59977 394 //-- Transition calculee au point 0 -> Trans2 , Trans1
395 //-- car ici, on devarit calculer en PI
396 Handle(IntPatch_GLine) glig = new IntPatch_GLine(elipsol,Standard_False,trans2,trans1);
fa0291ff 397 //
398 {
ecc4f148 399 Standard_Real aU1, aV1, aU2, aV2;
400 IntPatch_Point aIP;
401 gp_Pnt aP (ElCLib::Value(0., elipsol));
402 //
403 aIP.SetValue(aP,Tol,Standard_False);
404 aIP.SetMultiple(Standard_False);
405 //
406 Quad1.Parameters(aP, aU1, aV1);
407 Quad2.Parameters(aP, aU2, aV2);
408 aIP.SetParameters(aU1, aV1, aU2, aV2);
409 //
410 aIP.SetParameter(0.);
411 glig->AddVertex(aIP);
412 glig->SetFirstPoint(1);
413 //
414 aIP.SetParameter(2.*M_PI);
415 glig->AddVertex(aIP);
416 glig->SetLastPoint(2);
fa0291ff 417 }
418 //
419 pmult1.SetParameter(0.5*M_PI);
7fd59977 420 glig->AddVertex(pmult1);
fa0291ff 421 //
c6541a0c 422 pmult2.SetParameter(1.5*M_PI);
7fd59977 423 glig->AddVertex(pmult2);
fa0291ff 424
425 //
7fd59977 426 slin.Append(glig);
427
428 //-- Transitions calculee au point 0 OK
fa0291ff 429 //
7fd59977 430 // on traite la deuxieme ellipse
7fd59977 431 elipsol = inter.Ellipse(2);
432
433 Standard_Real param1 = ElCLib::Parameter(elipsol,pttang1);
434 Standard_Real param2 = ElCLib::Parameter(elipsol,pttang2);
ecc4f148 435 Standard_Real parampourtransition = 0.0;
436 if (param1 < param2)
437 {
438 pmult1.SetParameter(0.5*M_PI);
439 pmult2.SetParameter(1.5*M_PI);
440 parampourtransition = M_PI;
7fd59977 441 }
442 else {
ecc4f148 443 pmult1.SetParameter(1.5*M_PI);
444 pmult2.SetParameter(0.5*M_PI);
445 parampourtransition = 0.0;
7fd59977 446 }
447
448 //-- Calcul des transitions de ligne pour la premiere ligne
449 ElCLib::D1(parampourtransition,elipsol,ptref,Tgt);
450 qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
ecc4f148 451 if(qwe> 0.00000001)
452 {
453 trans1 = IntSurf_Out;
454 trans2 = IntSurf_In;
7fd59977 455 }
ecc4f148 456 else if(qwe< -0.00000001)
457 {
458 trans1 = IntSurf_In;
459 trans2 = IntSurf_Out;
7fd59977 460 }
ecc4f148 461 else
462 {
463 trans1=trans2=IntSurf_Undecided;
7fd59977 464 }
ecc4f148 465
7fd59977 466 //-- La transition a ete calculee sur un point de cette ligne
467 glig = new IntPatch_GLine(elipsol,Standard_False,trans1,trans2);
fa0291ff 468 //
469 {
ecc4f148 470 Standard_Real aU1, aV1, aU2, aV2;
471 IntPatch_Point aIP;
472 gp_Pnt aP (ElCLib::Value(0., elipsol));
473 //
474 aIP.SetValue(aP,Tol,Standard_False);
475 aIP.SetMultiple(Standard_False);
476 //
477 Quad1.Parameters(aP, aU1, aV1);
478 Quad2.Parameters(aP, aU2, aV2);
479 aIP.SetParameters(aU1, aV1, aU2, aV2);
480 //
481 aIP.SetParameter(0.);
482 glig->AddVertex(aIP);
483 glig->SetFirstPoint(1);
484 //
485 aIP.SetParameter(2.*M_PI);
486 glig->AddVertex(aIP);
487 glig->SetLastPoint(2);
7fd59977 488 }
fa0291ff 489 //
7fd59977 490 glig->AddVertex(pmult1);
fa0291ff 491 glig->AddVertex(pmult2);
492 //
7fd59977 493 slin.Append(glig);
494 }
495 break;
7fd59977 496
497 case IntAna_NoGeometricSolution:
498 {
499 Standard_Boolean bReverse;
500 Standard_Real U1,V1,U2,V2;
501 gp_Pnt psol;
502 //
503 bReverse=IsToReverse(Cy1, Cy2, Tol);
ecc4f148 504 if (bReverse)
505 {
506 Cy2=Quad1.Cylinder();
507 Cy1=Quad2.Cylinder();
7fd59977 508 }
509 //
510 IntAna_IntQuadQuad anaint(Cy1,Cy2,Tol);
ecc4f148 511 if (!anaint.IsDone())
512 {
513 return Standard_False;
7fd59977 514 }
515
ecc4f148 516 if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0)
517 {
518 Empty = Standard_True;
7fd59977 519 }
ecc4f148 520 else
521 {
522 NbSol = anaint.NbPnt();
523 for (i = 1; i <= NbSol; i++)
524 {
525 psol = anaint.Point(i);
526 Quad1.Parameters(psol,U1,V1);
527 Quad2.Parameters(psol,U2,V2);
528 ptsol.SetValue(psol,Tol,Standard_True);
529 ptsol.SetParameters(U1,V1,U2,V2);
530 spnt.Append(ptsol);
531 }
532
533 gp_Pnt ptvalid, ptf, ptl;
534 gp_Vec tgvalid;
535
536 Standard_Real first,last,para;
537 IntAna_Curve curvsol;
538 Standard_Boolean tgfound;
539 Standard_Integer kount;
540
541 NbSol = anaint.NbCurve();
542 for (i = 1; i <= NbSol; i++)
543 {
544 curvsol = anaint.Curve(i);
545 curvsol.Domain(first,last);
546 ptf = curvsol.Value(first);
547 ptl = curvsol.Value(last);
548
549 para = last;
550 kount = 1;
551 tgfound = Standard_False;
552
553 while (!tgfound)
554 {
555 para = (1.123*first + para)/2.123;
556 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
557 if (!tgfound)
558 {
559 kount ++;
560 tgfound = kount > 5;
561 }
562 }
563
564 Handle(IntPatch_ALine) alig;
565 if (kount <= 5)
566 {
567 Standard_Real qwe = tgvalid.DotCross( Quad2.Normale(ptvalid),
568 Quad1.Normale(ptvalid));
569 if(qwe>0.00000001)
570 {
571 trans1 = IntSurf_Out;
572 trans2 = IntSurf_In;
573 }
574 else if(qwe<-0.00000001)
575 {
576 trans1 = IntSurf_In;
577 trans2 = IntSurf_Out;
578 }
579 else
580 {
581 trans1=trans2=IntSurf_Undecided;
582 }
583 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
584 }
585 else
586 {
587 alig = new IntPatch_ALine(curvsol,Standard_False);
588 //-- cout << "Transition indeterminee" << endl;
589 }
590
591 Standard_Boolean TempFalse1 = Standard_False;
592 Standard_Boolean TempFalse2 = Standard_False;
593
594 ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1,ptf,first,
595 TempFalse2,ptl,last,Multpoint,Tol);
596 slin.Append(alig);
597 }
7fd59977 598 }
599 }
600 break;
601
602 default:
ecc4f148 603 return Standard_False;
604 }
605
606 return Standard_True;
607}
608
609//=======================================================================
610//function : ShortCosForm
611//purpose : Represents theCosFactor*cosA+theSinFactor*sinA as
612// theCoeff*cos(A-theAngle) if it is possibly (all angles are
613// in radians).
614//=======================================================================
615static void ShortCosForm( const Standard_Real theCosFactor,
616 const Standard_Real theSinFactor,
617 Standard_Real& theCoeff,
618 Standard_Real& theAngle)
619{
620 theCoeff = sqrt(theCosFactor*theCosFactor+theSinFactor*theSinFactor);
621 theAngle = 0.0;
622 if(IsEqual(theCoeff, 0.0))
623 {
624 theAngle = 0.0;
625 return;
626 }
627
628 theAngle = acos(Abs(theCosFactor/theCoeff));
629
630 if(theSinFactor > 0.0)
631 {
632 if(IsEqual(theCosFactor, 0.0))
633 {
634 theAngle = M_PI/2.0;
635 }
636 else if(theCosFactor < 0.0)
637 {
638 theAngle = M_PI-theAngle;
639 }
640 }
641 else if(IsEqual(theSinFactor, 0.0))
642 {
643 if(theCosFactor < 0.0)
644 {
645 theAngle = M_PI;
646 }
647 }
648 if(theSinFactor < 0.0)
649 {
650 if(theCosFactor > 0.0)
651 {
652 theAngle = 2.0*M_PI-theAngle;
653 }
654 else if(IsEqual(theCosFactor, 0.0))
655 {
656 theAngle = 3.0*M_PI/2.0;
657 }
658 else if(theCosFactor < 0.0)
659 {
660 theAngle = M_PI+theAngle;
661 }
662 }
663}
664
665enum SearchBoundType
666{
667 SearchNONE = 0,
668 SearchV1 = 1,
669 SearchV2 = 2
670};
671
672//Stores equations coefficients
673struct stCoeffsValue
674{
675 stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&);
676
677 gp_Vec mVecA1;
678 gp_Vec mVecA2;
679 gp_Vec mVecB1;
680 gp_Vec mVecB2;
681 gp_Vec mVecC1;
682 gp_Vec mVecC2;
683 gp_Vec mVecD;
684
685 Standard_Real mK21; //sinU2
686 Standard_Real mK11; //sinU1
687 Standard_Real mL21; //cosU2
688 Standard_Real mL11; //cosU1
689 Standard_Real mM1; //Free member
690
691 Standard_Real mK22; //sinU2
692 Standard_Real mK12; //sinU1
693 Standard_Real mL22; //cosU2
694 Standard_Real mL12; //cosU1
695 Standard_Real mM2; //Free member
696
697 Standard_Real mK1;
698 Standard_Real mL1;
699 Standard_Real mK2;
700 Standard_Real mL2;
701
702 Standard_Real mFIV1;
703 Standard_Real mPSIV1;
704 Standard_Real mFIV2;
705 Standard_Real mPSIV2;
706
707 Standard_Real mB;
708 Standard_Real mC;
709 Standard_Real mFI1;
710 Standard_Real mFI2;
711};
712
713stCoeffsValue::stCoeffsValue( const gp_Cylinder& theCyl1,
714 const gp_Cylinder& theCyl2)
715{
716 const Standard_Real aNulValue = 0.01*Precision::PConfusion();
717
718 mVecA1 = -theCyl1.Radius()*theCyl1.XAxis().Direction();
719 mVecA2 = theCyl2.Radius()*theCyl2.XAxis().Direction();
720
721 mVecB1 = -theCyl1.Radius()*theCyl1.YAxis().Direction();
722 mVecB2 = theCyl2.Radius()*theCyl2.YAxis().Direction();
723
724 mVecC1 = theCyl1.Axis().Direction();
725 mVecC2 = -(theCyl2.Axis().Direction());
726
727 mVecD = theCyl2.Location().XYZ() - theCyl1.Location().XYZ();
728
729 enum CoupleOfEquation
730 {
731 COENONE = 0,
732 COE12 = 1,
733 COE23 = 2,
734 COE13 = 3
735 }aFoundCouple = COENONE;
736
737
738 Standard_Real aDetV1V2 = mVecC1.X()*mVecC2.Y()-mVecC1.Y()*mVecC2.X();
739
740 if(Abs(aDetV1V2) < aNulValue)
741 {
742 aDetV1V2 = mVecC1.Y()*mVecC2.Z()-mVecC1.Z()*mVecC2.Y();
743 if(Abs(aDetV1V2) < aNulValue)
744 {
745 aDetV1V2 = mVecC1.X()*mVecC2.Z()-mVecC1.Z()*mVecC2.X();
746 if(Abs(aDetV1V2) < aNulValue)
747 {
748 Standard_Failure::Raise("Error. Exception in divide by zerro (IntCyCyTrim)!!!!");
749 }
750 else
751 {
752 aFoundCouple = COE13;
753 }
754 }
755 else
756 {
757 aFoundCouple = COE23;
758 }
759 }
760 else
761 {
762 aFoundCouple = COE12;
763 }
764
765 switch(aFoundCouple)
766 {
767 case COE12:
768 break;
769 case COE23:
770 {
771 gp_Vec aVTemp = mVecA1;
772 mVecA1.SetX(aVTemp.Y());
773 mVecA1.SetY(aVTemp.Z());
774 mVecA1.SetZ(aVTemp.X());
775
776 aVTemp = mVecA2;
777 mVecA2.SetX(aVTemp.Y());
778 mVecA2.SetY(aVTemp.Z());
779 mVecA2.SetZ(aVTemp.X());
780
781 aVTemp = mVecB1;
782 mVecB1.SetX(aVTemp.Y());
783 mVecB1.SetY(aVTemp.Z());
784 mVecB1.SetZ(aVTemp.X());
785
786 aVTemp = mVecB2;
787 mVecB2.SetX(aVTemp.Y());
788 mVecB2.SetY(aVTemp.Z());
789 mVecB2.SetZ(aVTemp.X());
790
791 aVTemp = mVecC1;
792 mVecC1.SetX(aVTemp.Y());
793 mVecC1.SetY(aVTemp.Z());
794 mVecC1.SetZ(aVTemp.X());
795
796 aVTemp = mVecC2;
797 mVecC2.SetX(aVTemp.Y());
798 mVecC2.SetY(aVTemp.Z());
799 mVecC2.SetZ(aVTemp.X());
800
801 aVTemp = mVecD;
802 mVecD.SetX(aVTemp.Y());
803 mVecD.SetY(aVTemp.Z());
804 mVecD.SetZ(aVTemp.X());
805
806 break;
807 }
808 case COE13:
809 {
810 gp_Vec aVTemp = mVecA1;
811 mVecA1.SetY(aVTemp.Z());
812 mVecA1.SetZ(aVTemp.Y());
813
814 aVTemp = mVecA2;
815 mVecA2.SetY(aVTemp.Z());
816 mVecA2.SetZ(aVTemp.Y());
817
818 aVTemp = mVecB1;
819 mVecB1.SetY(aVTemp.Z());
820 mVecB1.SetZ(aVTemp.Y());
821
822 aVTemp = mVecB2;
823 mVecB2.SetY(aVTemp.Z());
824 mVecB2.SetZ(aVTemp.Y());
825
826 aVTemp = mVecC1;
827 mVecC1.SetY(aVTemp.Z());
828 mVecC1.SetZ(aVTemp.Y());
829
830 aVTemp = mVecC2;
831 mVecC2.SetY(aVTemp.Z());
832 mVecC2.SetZ(aVTemp.Y());
833
834 aVTemp = mVecD;
835 mVecD.SetY(aVTemp.Z());
836 mVecD.SetZ(aVTemp.Y());
837
838 break;
839 }
840 default:
841 break;
842 }
843
844 //------- For V1 (begin)
845 //sinU2
846 mK21 = (mVecC2.Y()*mVecB2.X()-mVecC2.X()*mVecB2.Y())/aDetV1V2;
847 //sinU1
848 mK11 = (mVecC2.Y()*mVecB1.X()-mVecC2.X()*mVecB1.Y())/aDetV1V2;
849 //cosU2
850 mL21 = (mVecC2.Y()*mVecA2.X()-mVecC2.X()*mVecA2.Y())/aDetV1V2;
851 //cosU1
852 mL11 = (mVecC2.Y()*mVecA1.X()-mVecC2.X()*mVecA1.Y())/aDetV1V2;
853 //Free member
854 mM1 = (mVecC2.Y()*mVecD.X()-mVecC2.X()*mVecD.Y())/aDetV1V2;
855 //------- For V1 (end)
856
857 //------- For V2 (begin)
858 //sinU2
859 mK22 = (mVecC1.X()*mVecB2.Y()-mVecC1.Y()*mVecB2.X())/aDetV1V2;
860 //sinU1
861 mK12 = (mVecC1.X()*mVecB1.Y()-mVecC1.Y()*mVecB1.X())/aDetV1V2;
862 //cosU2
863 mL22 = (mVecC1.X()*mVecA2.Y()-mVecC1.Y()*mVecA2.X())/aDetV1V2;
864 //cosU1
865 mL12 = (mVecC1.X()*mVecA1.Y()-mVecC1.Y()*mVecA1.X())/aDetV1V2;
866 //Free member
867 mM2 = (mVecC1.X()*mVecD.Y()-mVecC1.Y()*mVecD.X())/aDetV1V2;
868 //------- For V1 (end)
869
870 ShortCosForm(mL11, mK11, mK1, mFIV1);
871 ShortCosForm(mL21, mK21, mL1, mPSIV1);
872 ShortCosForm(mL12, mK12, mK2, mFIV2);
873 ShortCosForm(mL22, mK22, mL2, mPSIV2);
874
875 const Standard_Real aA1=mVecC1.Z()*mK21+mVecC2.Z()*mK22-mVecB2.Z(), //sinU2
876 aA2=mVecC1.Z()*mL21+mVecC2.Z()*mL22-mVecA2.Z(), //cosU2
877 aB1=mVecB1.Z()-mVecC1.Z()*mK11-mVecC2.Z()*mK12, //sinU1
878 aB2=mVecA1.Z()-mVecC1.Z()*mL11-mVecC2.Z()*mL12; //cosU1
879
880 mC =mVecD.Z() -mVecC1.Z()*mM1 -mVecC2.Z()*mM2; //Free
881
882 Standard_Real aA = 0.0;
883
884 ShortCosForm(aB2,aB1,mB,mFI1);
885 ShortCosForm(aA2,aA1,aA,mFI2);
886
887 mB /= aA;
888 mC /= aA;
889}
890
891//=======================================================================
892//function : SearchOnVBounds
893//purpose :
894//=======================================================================
895static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType,
896 const stCoeffsValue& theCoeffs,
897 const Standard_Real theVzad,
898 const Standard_Real theInitU2,
899 const Standard_Real theInitMainVar,
900 Standard_Real& theMainVariableValue)
901{
902 const Standard_Real aMaxError = 4.0*M_PI; // two periods
903
904 theMainVariableValue = RealLast();
905 const Standard_Real aTol2 = Precision::PConfusion()*Precision::PConfusion();
906 Standard_Real aMainVarPrev = theInitMainVar, aU2Prev = theInitU2, anOtherVar = theVzad;
907
908 Standard_Real anError = RealLast();
909 Standard_Integer aNbIter = 0;
910 do
911 {
912 if(++aNbIter > 1000)
913 return Standard_False;
914
915 gp_Vec aVecMainVar = theCoeffs.mVecA1 * sin(aMainVarPrev) - theCoeffs.mVecB1 * cos(aMainVarPrev);
916 gp_Vec aVecFreeMem = (theCoeffs.mVecA2 * aU2Prev + theCoeffs.mVecB2) * sin(aU2Prev) -
917 (theCoeffs.mVecB2 * aU2Prev - theCoeffs.mVecA2) * cos(aU2Prev) +
918 (theCoeffs.mVecA1 * aMainVarPrev + theCoeffs.mVecB1) * sin(aMainVarPrev) -
919 (theCoeffs.mVecB1 * aMainVarPrev - theCoeffs.mVecA1) * cos(aMainVarPrev) + theCoeffs.mVecD;
920
921 gp_Vec aVecVar1 = theCoeffs.mVecA2 * sin(aU2Prev) - theCoeffs.mVecB2 * cos(aU2Prev);
922 gp_Vec aVecVar2;
923
924 switch(theSBType)
925 {
926 case SearchV1:
927 aVecVar2 = theCoeffs.mVecC2;
928 aVecFreeMem -= theCoeffs.mVecC1 * theVzad;
929 break;
930
931 case SearchV2:
932 aVecVar2 = theCoeffs.mVecC1;
933 aVecFreeMem -= theCoeffs.mVecC2 * theVzad;
934 break;
935
936 default:
937 return Standard_False;
938 }
939
940 Standard_Real aDetMainSyst = aVecMainVar.X()*(aVecVar1.Y()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.Y())-
941 aVecMainVar.Y()*(aVecVar1.X()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.X())+
942 aVecMainVar.Z()*(aVecVar1.X()*aVecVar2.Y()-aVecVar1.Y()*aVecVar2.X());
943
944 if(IsEqual(aDetMainSyst, 0.0))
945 {
946 return Standard_False;
947 }
948
949
950 Standard_Real aDetMainVar = aVecFreeMem.X()*(aVecVar1.Y()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.Y())-
951 aVecFreeMem.Y()*(aVecVar1.X()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.X())+
952 aVecFreeMem.Z()*(aVecVar1.X()*aVecVar2.Y()-aVecVar1.Y()*aVecVar2.X());
953
954 Standard_Real aDetVar1 = aVecMainVar.X()*(aVecFreeMem.Y()*aVecVar2.Z()-aVecFreeMem.Z()*aVecVar2.Y())-
955 aVecMainVar.Y()*(aVecFreeMem.X()*aVecVar2.Z()-aVecFreeMem.Z()*aVecVar2.X())+
956 aVecMainVar.Z()*(aVecFreeMem.X()*aVecVar2.Y()-aVecFreeMem.Y()*aVecVar2.X());
957
958 Standard_Real aDetVar2 = aVecMainVar.X()*(aVecVar1.Y()*aVecFreeMem.Z()-aVecVar1.Z()*aVecFreeMem.Y())-
959 aVecMainVar.Y()*(aVecVar1.X()*aVecFreeMem.Z()-aVecVar1.Z()*aVecFreeMem.X())+
960 aVecMainVar.Z()*(aVecVar1.X()*aVecFreeMem.Y()-aVecVar1.Y()*aVecFreeMem.X());
961
962 Standard_Real aDelta = aDetMainVar/aDetMainSyst-aMainVarPrev;
963
964 if(Abs(aDelta) > aMaxError)
965 return Standard_False;
966
967 anError = aDelta*aDelta;
968 aMainVarPrev += aDelta;
969
970 ///
971 aDelta = aDetVar1/aDetMainSyst-aU2Prev;
972
973 if(Abs(aDelta) > aMaxError)
974 return Standard_False;
975
976 anError += aDelta*aDelta;
977 aU2Prev += aDelta;
978
979 ///
980 aDelta = aDetVar2/aDetMainSyst-anOtherVar;
981 anError += aDelta*aDelta;
982 anOtherVar += aDelta;
983 }
984 while(anError > aTol2);
985
986 theMainVariableValue = aMainVarPrev;
987
988 return Standard_True;
989}
990
991//=======================================================================
992//function : InscribePoint
993//purpose :
994//=======================================================================
995static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
996 const Standard_Real theUlTarget,
997 Standard_Real& theUGiven,
998 const Standard_Real theTol2D,
999 const Standard_Real thePeriod)
1000{
1001 if((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D))
1002 {//it has already inscribed
1003
1004 /*
1005 Utf U Utl
1006 + * +
1007 */
1008
1009 return Standard_True;
1010 }
1011
1012 if(IsEqual(thePeriod, 0.0))
1013 {//it cannot be inscribed
1014 return Standard_False;
1015 }
1016
1017 Standard_Real aDU = theUGiven - theUfTarget;
1018
1019 if(aDU > 0.0)
1020 aDU = -thePeriod;
1021 else
1022 aDU = +thePeriod;
1023
1024 while(((theUGiven - theUfTarget)*aDU < 0.0) && !((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D)))
1025 {
1026 theUGiven += aDU;
1027 }
1028
1029 return ((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D));
1030}
1031
1032//=======================================================================
1033//function : InscribeInterval
1034//purpose : Adjusts theUfGiven and after that fits theUlGiven to result
1035//=======================================================================
1036static Standard_Boolean InscribeInterval(const Standard_Real theUfTarget,
1037 const Standard_Real theUlTarget,
1038 Standard_Real& theUfGiven,
1039 Standard_Real& theUlGiven,
1040 const Standard_Real theTol2D,
1041 const Standard_Real thePeriod)
1042{
1043 Standard_Real anUpar = theUfGiven;
1044 const Standard_Real aDelta = theUlGiven - theUfGiven;
1045 if(InscribePoint(theUfTarget, theUlTarget, anUpar, theTol2D, thePeriod))
1046 {
1047 theUfGiven = anUpar;
1048 theUlGiven = theUfGiven + aDelta;
1049 }
1050 else
1051 {
1052 anUpar = theUlGiven;
1053 if(InscribePoint(theUfTarget, theUlTarget, anUpar, theTol2D, thePeriod))
1054 {
1055 theUlGiven = anUpar;
1056 theUfGiven = theUlGiven - aDelta;
1057 }
1058 else
7fd59977 1059 {
1060 return Standard_False;
1061 }
1062 }
ecc4f148 1063
1064 return Standard_True;
1065}
1066
1067//=======================================================================
1068//function : InscribeAndSortArray
1069//purpose : Sort from Min to Max value
1070//=======================================================================
1071static void InscribeAndSortArray( Standard_Real theArr[],
1072 const Standard_Integer theNOfMembers,
1073 const Standard_Real theUf,
1074 const Standard_Real theUl,
1075 const Standard_Real theTol2D,
1076 const Standard_Real thePeriod)
1077{
1078 for(Standard_Integer i = 0; i < theNOfMembers; i++)
1079 {
1080 if(Precision::IsInfinite(theArr[i]))
1081 {
1082 if(theArr[i] < 0.0)
1083 theArr[i] = -theArr[i];
1084
1085 continue;
1086 }
1087
1088 InscribePoint(theUf, theUl, theArr[i], theTol2D, thePeriod);
1089
1090 for(Standard_Integer j = i - 1; j >= 0; j--)
1091 {
1092
1093 if(theArr[j+1] < theArr[j])
1094 {
1095 Standard_Real anUtemp = theArr[j+1];
1096 theArr[j+1] = theArr[j];
1097 theArr[j] = anUtemp;
1098 }
1099 }
1100 }
1101}
1102
1103
1104//=======================================================================
1105//function : AddPointIntoWL
1106//purpose : Surf1 is a surface, whose U-par is variable.
1107//=======================================================================
1108static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
1109 const IntSurf_Quadric& theQuad2,
1110 const Standard_Boolean isTheReverse,
1111 const gp_Pnt2d& thePntOnSurf1,
1112 const gp_Pnt2d& thePntOnSurf2,
1113 const Standard_Real theUfSurf1,
1114 const Standard_Real theUlSurf1,
1115 const Standard_Real thePeriodOfSurf1,
1116 const Handle(IntSurf_LineOn2S)& theLine,
1117 const Standard_Real theTol2D)
1118{
1119 gp_Pnt aPt1(theQuad1.Value(thePntOnSurf1.X(), thePntOnSurf1.Y())),
1120 aPt2(theQuad2.Value(thePntOnSurf2.X(), thePntOnSurf2.Y()));
1121
1122 Standard_Real anUpar = thePntOnSurf1.X();
1123 if(!InscribePoint(theUfSurf1, theUlSurf1, anUpar, theTol2D, thePeriodOfSurf1))
1124 return Standard_False;
1125
1126 IntSurf_PntOn2S aPnt;
1127
1128 if(isTheReverse)
1129 {
1130 aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
1131 thePntOnSurf2.X(), thePntOnSurf2.Y(),
1132 thePntOnSurf1.X(), thePntOnSurf1.Y());
1133 }
1134 else
1135 {
1136 aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
1137 thePntOnSurf1.X(), thePntOnSurf1.Y(),
1138 thePntOnSurf2.X(), thePntOnSurf2.Y());
1139 }
1140
1141 theLine->Add(aPnt);
1142 return Standard_True;
1143}
1144
1145//=======================================================================
1146//function : AddBoundaryPoint
1147//purpose :
1148//=======================================================================
1149static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1,
1150 const IntSurf_Quadric& theQuad2,
1151 const Handle(IntPatch_WLine)& theWL,
1152 const stCoeffsValue& theCoeffs,
1153 const Bnd_Box2d& theUVSurf1,
1154 const Bnd_Box2d& theUVSurf2,
1155 const Standard_Real theTol2D,
1156 const Standard_Real thePeriod,
1157 const Standard_Real theNulValue,
1158 const Standard_Real theU1,
1159 const Standard_Real theU2,
1160 const Standard_Real theV1,
1161 const Standard_Real theV1Prev,
1162 const Standard_Real theV2,
1163 const Standard_Real theV2Prev,
1164 const Standard_Boolean isTheReverse,
1165 const Standard_Real theArccosFactor,
1166 Standard_Boolean& isTheFound1,
1167 Standard_Boolean& isTheFound2)
1168{
1169 Standard_Real aUSurf1f = 0.0, //const
1170 aUSurf1l = 0.0,
1171 aVSurf1f = 0.0,
1172 aVSurf1l = 0.0;
1173 Standard_Real aUSurf2f = 0.0, //const
1174 aUSurf2l = 0.0,
1175 aVSurf2f = 0.0,
1176 aVSurf2l = 0.0;
1177
1178 theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
1179 theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
1180
1181 SearchBoundType aTS1 = SearchNONE, aTS2 = SearchNONE;
1182 Standard_Real aV1zad = aVSurf1f, aV2zad = aVSurf2f;
1183
1184 Standard_Real anUpar1 = theU1, anUpar2 = theU1;
1185 Standard_Real aVf = theV1, aVl = theV1Prev;
1186 MinMax(aVf, aVl);
1187 if((aVf <= aVSurf1f) && (aVSurf1f <= aVl))
1188 {
1189 aTS1 = SearchV1;
1190 aV1zad = aVSurf1f;
1191 isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1f, theU2, theU1, anUpar1);
1192 }
1193 else if((aVf <= aVSurf1l) && (aVSurf1l <= aVl))
1194 {
1195 aTS1 = SearchV1;
1196 aV1zad = aVSurf1l;
1197 isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1l, theU2, theU1, anUpar1);
1198 }
1199
1200 aVf = theV2;
1201 aVl = theV2Prev;
1202 MinMax(aVf, aVl);
1203
1204 if((aVf <= aVSurf2f) && (aVSurf2f <= aVl))
1205 {
1206 aTS2 = SearchV2;
1207 aV2zad = aVSurf2f;
1208 isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2f, theU2, theU1, anUpar2);
1209 }
1210 else if((aVf <= aVSurf2l) && (aVSurf2l <= aVl))
1211 {
1212 aTS2 = SearchV2;
1213 aV2zad = aVSurf2l;
1214 isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2l, theU2, theU1, anUpar2);
1215 }
1216
1217 if(anUpar1 < anUpar2)
1218 {
1219 if(isTheFound1)
1220 {
1221 Standard_Real anArg = theCoeffs.mB * cos(anUpar1 - theCoeffs.mFI1) + theCoeffs.mC;
1222
1223 if(theNulValue > 1.0 - anArg)
1224 anArg = 1.0;
1225 if(anArg + 1.0 < theNulValue)
1226 anArg = -1.0;
1227
1228 Standard_Real aU2 = theCoeffs.mFI2 + theArccosFactor * acos(anArg);
1229
1230 if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod))
1231 {
1232 const Standard_Real aV1 = (aTS1 == SearchV1) ? aV1zad :
1233 theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar1) +
1234 theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar1) + theCoeffs.mM1;
1235 const Standard_Real aV2 = (aTS1 == SearchV2) ? aV2zad :
1236 theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar1) +
1237 theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar1) + theCoeffs.mM2;
1238
1239 AddPointIntoWL(theQuad1, theQuad2, isTheReverse, gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2), aUSurf1f, aUSurf1l, thePeriod, theWL->Curve(), theTol2D);
1240 }
1241 else
1242 {
1243 isTheFound1 = Standard_False;
1244 }
1245 }
1246
1247 if(isTheFound2)
1248 {
1249 Standard_Real anArg = theCoeffs.mB * cos(anUpar2 - theCoeffs.mFI1) + theCoeffs.mC;
1250
1251 if(theNulValue > 1.0 - anArg)
1252 anArg = 1.0;
1253 if(anArg + 1.0 < theNulValue)
1254 anArg = -1.0;
1255
1256 Standard_Real aU2 = theCoeffs.mFI2 + theArccosFactor * acos(anArg);
1257 if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod))
1258 {
1259 const Standard_Real aV1 = (aTS2 == SearchV1) ? aV1zad :
1260 theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar2) +
1261 theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar2) + theCoeffs.mM1;
1262 const Standard_Real aV2 = (aTS2 == SearchV2) ? aV2zad :
1263 theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar2) +
1264 theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar2) + theCoeffs.mM2;
1265
1266 AddPointIntoWL(theQuad1, theQuad2, isTheReverse, gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2), aUSurf1f, aUSurf1l, thePeriod, theWL->Curve(), theTol2D);
1267 }
1268 else
1269 {
1270 isTheFound2 = Standard_False;
1271 }
1272 }
1273 }
1274 else
1275 {
1276 if(isTheFound2)
1277 {
1278 Standard_Real anArg = theCoeffs.mB * cos(anUpar2 - theCoeffs.mFI1) + theCoeffs.mC;
1279
1280 if(theNulValue > 1.0 - anArg)
1281 anArg = 1.0;
1282 if(anArg + 1.0 < theNulValue)
1283 anArg = -1.0;
1284
1285 Standard_Real aU2 = theCoeffs.mFI2 + theArccosFactor * acos(anArg);
1286
1287 if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod))
1288 {
1289 const Standard_Real aV1 = (aTS2 == SearchV1) ? aV1zad :
1290 theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar2) +
1291 theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar2) + theCoeffs.mM1;
1292 const Standard_Real aV2 = (aTS2 == SearchV2) ? aV2zad :
1293 theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar2) +
1294 theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar2) + theCoeffs.mM2;
1295
1296 AddPointIntoWL(theQuad1, theQuad2, isTheReverse, gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2), aUSurf1f, aUSurf1l, thePeriod, theWL->Curve(), theTol2D);
1297 }
1298 else
1299 {
1300 isTheFound2 = Standard_False;
1301 }
1302 }
1303
1304 if(isTheFound1)
1305 {
1306 Standard_Real anArg = theCoeffs.mB*cos(anUpar1-theCoeffs.mFI1)+theCoeffs.mC;
1307
1308 if(theNulValue > 1.0 - anArg)
1309 anArg = 1.0;
1310 if(anArg + 1.0 < theNulValue)
1311 anArg = -1.0;
1312
1313 Standard_Real aU2 = theCoeffs.mFI2+theArccosFactor*acos(anArg);
1314 if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod))
1315 {
1316 const Standard_Real aV1 = (aTS1 == SearchV1) ? aV1zad :
1317 theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar1) +
1318 theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar1) + theCoeffs.mM1;
1319 const Standard_Real aV2 = (aTS1 == SearchV2) ? aV2zad :
1320 theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar1) +
1321 theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar1) + theCoeffs.mM2;
1322
1323 AddPointIntoWL(theQuad1, theQuad2, isTheReverse, gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2), aUSurf1f, aUSurf1l, thePeriod, theWL->Curve(), theTol2D);
1324 }
1325 else
1326 {
1327 isTheFound1 = Standard_False;
1328 }
1329 }
1330 }
1331
7fd59977 1332 return Standard_True;
1333}
1334
ecc4f148 1335//=======================================================================
1336//function : SeekAdditionalPoints
1337//purpose :
1338//=======================================================================
1339static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
1340 const IntSurf_Quadric& theQuad2,
1341 const Handle(IntSurf_LineOn2S)& theLile,
1342 const stCoeffsValue& theCoeffs,
1343 const Standard_Integer theMinNbPoints,
1344 const Standard_Real theU2f,
1345 const Standard_Real theU2l,
1346 const Standard_Real theTol2D,
1347 const Standard_Real thePeriodOfSurf2,
1348 const Standard_Real theArccosFactor,
1349 const Standard_Boolean isTheReverse)
1350{
1351 Standard_Integer aNbPoints = theLile->NbPoints();
1352 if(aNbPoints >= theMinNbPoints)
1353 {
1354 return;
1355 }
1356
1357 Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0;
1358
1359 Standard_Integer aNbPointsPrev = 0;
1360 while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
1361 {
1362 aNbPointsPrev = aNbPoints;
1363 for(Standard_Integer fp = 1, lp = 2; fp < aNbPoints; fp = lp + 1)
1364 {
1365 Standard_Real U1f, V1f; //first point in 1st suraface
1366 Standard_Real U1l, V1l; //last point in 1st suraface
1367
1368 lp = fp+1;
1369
1370 if(isTheReverse)
1371 {
1372 theLile->Value(fp).ParametersOnS2(U1f, V1f);
1373 theLile->Value(lp).ParametersOnS2(U1l, V1l);
1374 }
1375 else
1376 {
1377 theLile->Value(fp).ParametersOnS1(U1f, V1f);
1378 theLile->Value(lp).ParametersOnS1(U1l, V1l);
1379 }
1380
1381 U1prec = 0.5*(U1f+U1l);
1382
1383 Standard_Real anArg = theCoeffs.mB * cos(U1prec - theCoeffs.mFI1) + theCoeffs.mC;
1384 if(anArg > 1.0)
1385 anArg = 1.0;
1386 if(anArg < -1.0)
1387 anArg = -1.0;
1388
1389 U2prec = theCoeffs.mFI2 + theArccosFactor*acos(anArg);
1390 InscribePoint(theU2f, theU2l, U2prec, theTol2D, thePeriodOfSurf2);
1391
1392 gp_Pnt aP1, aP2;
1393 gp_Vec aVec1, aVec2;
1394
1395 if(isTheReverse)
1396 {
1397 V1prec = theCoeffs.mK21 * sin(U2prec) + theCoeffs.mK11 * sin(U1prec) + theCoeffs.mL21 * cos(U2prec) + theCoeffs.mL11 * cos(U1prec) + theCoeffs.mM1;
1398 V2prec = theCoeffs.mK22 * sin(U2prec) + theCoeffs.mK12 * sin(U1prec) + theCoeffs.mL22 * cos(U2prec) + theCoeffs.mL12 * cos(U1prec) + theCoeffs.mM2;
1399
1400 gp_Pnt aP3, aP4;
1401 theQuad1.D1(U2prec, V2prec, aP3, aVec1, aVec2);
1402 theQuad2.D1(U1prec, V1prec, aP4, aVec1, aVec2);
1403
1404 theQuad1.D1(U1prec, V1prec, aP1, aVec1, aVec2);
1405 theQuad2.D1(U2prec, V2prec, aP2, aVec1, aVec2);
1406 }
1407 else
1408 {
1409 V1prec = theCoeffs.mK21 * sin(U2prec) + theCoeffs.mK11 * sin(U1prec) + theCoeffs.mL21 * cos(U2prec) + theCoeffs.mL11 * cos(U1prec) + theCoeffs.mM1;
1410 V2prec = theCoeffs.mK22 * sin(U2prec) + theCoeffs.mK12 * sin(U1prec) + theCoeffs.mL22 * cos(U2prec) + theCoeffs.mL12 * cos(U1prec) + theCoeffs.mM2;
1411
1412 theQuad1.D1(U1prec, V1prec, aP1, aVec1, aVec2);
1413 theQuad2.D1(U2prec, V2prec, aP2, aVec1, aVec2);
1414 }
1415
1416 gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
1417
1418 IntSurf_PntOn2S anIP;
1419 if(isTheReverse)
1420 {
1421 anIP.SetValue(aPInt, U2prec, V2prec, U1prec, V1prec);
1422 }
1423 else
1424 {
1425 anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec);
1426 }
1427
1428 theLile->InsertBefore(lp, anIP);
1429
1430 aNbPoints = theLile->NbPoints();
1431 if(aNbPoints >= theMinNbPoints)
1432 {
1433 return;
1434 }
1435 }
1436 }
1437}
1438
1439//=======================================================================
1440//function : CriticalPointsComputing
1441//purpose :
1442//=======================================================================
1443static void CriticalPointsComputing(const stCoeffsValue& theCoeffs,
1444 const Standard_Real theUSurf1f,
1445 const Standard_Real theUSurf1l,
1446 const Standard_Real theUSurf2f,
1447 const Standard_Real theUSurf2l,
1448 const Standard_Real thePeriod,
1449 const Standard_Real theTol2D,
1450 const Standard_Integer theNbCritPointsMax,
1451 Standard_Real theU1crit[])
1452{
1453 theU1crit[0] = 0.0;
1454 theU1crit[1] = thePeriod;
1455 theU1crit[2] = theUSurf1f;
1456 theU1crit[3] = theUSurf1l;
1457
1458 const Standard_Real aCOS = cos(theCoeffs.mFI2);
1459 const Standard_Real aBSB = Abs(theCoeffs.mB);
1460 if((theCoeffs.mC - aBSB <= aCOS) && (aCOS <= theCoeffs.mC + aBSB))
1461 {
1462 Standard_Real anArg = (aCOS - theCoeffs.mC) / theCoeffs.mB;
1463 if(anArg > 1.0)
1464 anArg = 1.0;
1465 if(anArg < -1.0)
1466 anArg = -1.0;
1467
1468 theU1crit[4] = -acos(anArg) + theCoeffs.mFI1;
1469 theU1crit[5] = acos(anArg) + theCoeffs.mFI1;
1470 }
1471
1472 Standard_Real aSf = cos(theUSurf2f - theCoeffs.mFI2);
1473 Standard_Real aSl = cos(theUSurf2l - theCoeffs.mFI2);
1474 MinMax(aSf, aSl);
1475
1476 theU1crit[6] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? -acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : -Precision::Infinite();
1477 theU1crit[7] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? -acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : Precision::Infinite();
1478 theU1crit[8] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : -Precision::Infinite();
1479 theU1crit[9] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : Precision::Infinite();
1480
1481 //preparative treatment of array
1482 InscribeAndSortArray(theU1crit, theNbCritPointsMax, 0.0, thePeriod, theTol2D, thePeriod);
1483 for(Standard_Integer i = 1; i < theNbCritPointsMax; i++)
1484 {
1485 Standard_Real &a = theU1crit[i],
1486 &b = theU1crit[i-1];
1487 if(Abs(a - b) < theTol2D)
1488 {
1489 a = (a + b)/2.0;
1490 b = Precision::Infinite();
1491 }
1492 }
1493}
1494
1495//=======================================================================
1496//function : IntCyCyTrim
1497//purpose :
1498//=======================================================================
1499Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
1500 const IntSurf_Quadric& theQuad2,
1501 const Standard_Real theTol3D,
1502 const Standard_Real theTol2D,
1503 const Bnd_Box2d& theUVSurf1,
1504 const Bnd_Box2d& theUVSurf2,
1505 const Standard_Boolean isTheReverse,
1506 Standard_Boolean& isTheEmpty,
1507 IntPatch_SequenceOfLine& theSlin,
1508 IntPatch_SequenceOfPoint& theSPnt)
1509{
1510 Standard_Real aUSurf1f = 0.0, //const
1511 aUSurf1l = 0.0,
1512 aVSurf1f = 0.0,
1513 aVSurf1l = 0.0;
1514 Standard_Real aUSurf2f = 0.0, //const
1515 aUSurf2l = 0.0,
1516 aVSurf2f = 0.0,
1517 aVSurf2l = 0.0;
1518
1519 theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
1520 theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
1521
1522 const Standard_Real aNulValue = 0.01*Precision::PConfusion();
1523
1524 const gp_Cylinder& aCyl1 = theQuad1.Cylinder(),
1525 aCyl2 = theQuad2.Cylinder();
1526
1527 IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
1528
1529 if (!anInter.IsDone())
1530 {
1531 return Standard_False;
1532 }
1533
1534 IntAna_ResultType aTypInt = anInter.TypeInter();
1535
1536 if(aTypInt != IntAna_NoGeometricSolution)
1537 { //It is not necessary (because result is an analytic curve) or
1538 //it is impossible to make Walking-line.
1539
1540 return Standard_False;
1541 }
1542
1543 theSlin.Clear();
1544 theSPnt.Clear();
1545 const Standard_Integer aNbPoints = Min(Max(200, RealToInt(20.0*aCyl1.Radius())), 2000);
1546 const Standard_Real aPeriod = 2.0*M_PI;
1547 const Standard_Real aStepMin = theTol2D,
1548 aStepMax = (aUSurf1l-aUSurf1f)/IntToReal(aNbPoints);
1549
1550 const stCoeffsValue anEquationCoeffs(aCyl1, aCyl2);
1551
1552 //Boundaries
1553 const Standard_Integer aNbOfBoundaries = 2;
1554 Standard_Real aU1f[aNbOfBoundaries] = {-Precision::Infinite(), -Precision::Infinite()};
1555 Standard_Real aU1l[aNbOfBoundaries] = {Precision::Infinite(), Precision::Infinite()};
1556
1557 if(anEquationCoeffs.mB > 0.0)
1558 {
1559 if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) < -1.0)
1560 {//There is NOT intersection
1561 return Standard_True;
1562 }
1563 else if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
1564 {//U=[0;2*PI]+aFI1
1565 aU1f[0] = anEquationCoeffs.mFI1;
1566 aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
1567 }
1568 else if((1 + anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
1569 (anEquationCoeffs.mB <= 1 - anEquationCoeffs.mC))
1570 {
1571 Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
1572 if(anArg > 1.0)
1573 anArg = 1.0;
1574 if(anArg < -1.0)
1575 anArg = -1.0;
1576
1577 const Standard_Real aDAngle = acos(anArg);
1578 //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
1579 aU1f[0] = anEquationCoeffs.mFI1;
1580 aU1l[0] = aDAngle + anEquationCoeffs.mFI1;
1581 aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
1582 aU1l[1] = aPeriod + anEquationCoeffs.mFI1;
1583 }
1584 else if((1 - anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
1585 (anEquationCoeffs.mB <= 1 + anEquationCoeffs.mC))
1586 {
1587 Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
1588 if(anArg > 1.0)
1589 anArg = 1.0;
1590 if(anArg < -1.0)
1591 anArg = -1.0;
1592
1593 const Standard_Real aDAngle = acos(anArg);
1594 //U=[aDAngle;2*PI-aDAngle]+aFI1
1595
1596 aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
1597 aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
1598 }
1599 else if(anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
1600 {
1601 Standard_Real anArg1 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB,
1602 anArg2 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
1603 if(anArg1 > 1.0)
1604 anArg1 = 1.0;
1605 if(anArg1 < -1.0)
1606 anArg1 = -1.0;
1607
1608 if(anArg2 > 1.0)
1609 anArg2 = 1.0;
1610 if(anArg2 < -1.0)
1611 anArg2 = -1.0;
1612
1613 const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
1614 //(U=[aDAngle1;aDAngle2]+aFI1) ||
1615 //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
1616
1617 aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1;
1618 aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1;
1619 aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1;
1620 aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1;
1621 }
1622 else
1623 {
1624 Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
1625 }
1626 }
1627 else if(anEquationCoeffs.mB < 0.0)
1628 {
1629 if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) > 1.0)
1630 {//There is NOT intersection
1631 return Standard_True;
1632 }
1633 else if(-anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
1634 {//U=[0;2*PI]+aFI1
1635 aU1f[0] = anEquationCoeffs.mFI1;
1636 aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
1637 }
1638 else if((-anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
1639 ( anEquationCoeffs.mB <= anEquationCoeffs.mC - 1))
1640 {
1641 Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
1642 if(anArg > 1.0)
1643 anArg = 1.0;
1644 if(anArg < -1.0)
1645 anArg = -1.0;
1646
1647 const Standard_Real aDAngle = acos(anArg);
1648 //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
1649
1650 aU1f[0] = anEquationCoeffs.mFI1;
1651 aU1l[0] = aDAngle + anEquationCoeffs.mFI1;
1652 aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
1653 aU1l[1] = aPeriod + anEquationCoeffs.mFI1;
1654 }
1655 else if((anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
1656 (anEquationCoeffs.mB <= -anEquationCoeffs.mB - 1))
1657 {
1658 Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
1659 if(anArg > 1.0)
1660 anArg = 1.0;
1661 if(anArg < -1.0)
1662 anArg = -1.0;
1663
1664 const Standard_Real aDAngle = acos(anArg);
1665 //U=[aDAngle;2*PI-aDAngle]+aFI1
1666
1667 aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
1668 aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
1669 }
1670 else if(-anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
1671 {
1672 Standard_Real anArg1 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB,
1673 anArg2 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
1674 if(anArg1 > 1.0)
1675 anArg1 = 1.0;
1676 if(anArg1 < -1.0)
1677 anArg1 = -1.0;
1678
1679 if(anArg2 > 1.0)
1680 anArg2 = 1.0;
1681 if(anArg2 < -1.0)
1682 anArg2 = -1.0;
1683
1684 const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
1685 //(U=[aDAngle1;aDAngle2]+aFI1) ||
1686 //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
1687
1688 aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1;
1689 aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1;
1690 aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1;
1691 aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1;
1692 }
1693 else
1694 {
1695 Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
1696 }
1697 }
1698 else
1699 {
1700 Standard_Failure::Raise("Error. Exception. Unhandled case (B-parameter computation)!");
1701 }
1702
1703 for(Standard_Integer i = 0; i < aNbOfBoundaries; i++)
1704 {
1705 if(Precision::IsInfinite(aU1f[i]) && Precision::IsInfinite(aU1l[i]))
1706 continue;
1707
1708 InscribeInterval(aUSurf1f, aUSurf1l, aU1f[i], aU1l[i], theTol2D, aPeriod);
1709 }
1710
1711 if( !Precision::IsInfinite(aU1f[0]) && !Precision::IsInfinite(aU1f[1]) &&
1712 !Precision::IsInfinite(aU1l[0]) && !Precision::IsInfinite(aU1l[1]))
1713 {
1714 if( ((aU1f[1] <= aU1l[0]) || (aU1l[1] <= aU1l[0])) &&
1715 ((aU1f[0] <= aU1l[1]) || (aU1l[0] <= aU1l[1])))
1716 {//Join all intervals to one
1717 aU1f[0] = Min(aU1f[0], aU1f[1]);
1718 aU1l[0] = Max(aU1l[0], aU1l[1]);
1719
1720 aU1f[1] = -Precision::Infinite();
1721 aU1l[1] = Precision::Infinite();
1722 }
1723 }
1724
1725 //Critical points
1726 //[0...1] - in these points parameter U1 goes through
1727 // the seam-edge of the first cylinder.
1728 //[2...3] - First and last U1 parameter.
1729 //[4...5] - in these points parameter U2 goes through
1730 // the seam-edge of the second cylinder.
1731 //[6...9] - in these points an intersetion line goes through
1732 // U-boundaries of the second surface.
1733 const Standard_Integer aNbCritPointsMax = 10;
1734 Standard_Real anU1crit[aNbCritPointsMax] = {Precision::Infinite(),
1735 Precision::Infinite(),
1736 Precision::Infinite(),
1737 Precision::Infinite(),
1738 Precision::Infinite(),
1739 Precision::Infinite(),
1740 Precision::Infinite(),
1741 Precision::Infinite(),
1742 Precision::Infinite(),
1743 Precision::Infinite()};
1744
1745 CriticalPointsComputing(anEquationCoeffs, aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
1746 aPeriod, theTol2D, aNbCritPointsMax, anU1crit);
1747
1748
1749 //Getting Walking-line
1750
1751 for(Standard_Integer aCurInterval = 0; aCurInterval < aNbOfBoundaries; aCurInterval++)
1752 {
1753 if(Precision::IsInfinite(aU1f[aCurInterval]) && Precision::IsInfinite(aU1l[aCurInterval]))
1754 continue;
1755
1756 Standard_Boolean isAddedIntoWL1 = Standard_False, isAddedIntoWL2 = Standard_False;
1757
1758 Standard_Real anUf = aU1f[aCurInterval], anUl = aU1l[aCurInterval];
1759
1760 //Inscribe and sort critical points
1761 InscribeAndSortArray(anU1crit, aNbCritPointsMax, anUf, anUl, theTol2D, aPeriod);
1762
1763 while(anUf < anUl)
1764 {
1765 Handle(IntSurf_LineOn2S) aL2S1 = new IntSurf_LineOn2S();
1766 Handle(IntSurf_LineOn2S) aL2S2 = new IntSurf_LineOn2S();
1767
1768 Handle(IntPatch_WLine) aWLine1 = new IntPatch_WLine(aL2S1, Standard_False);
1769 Handle(IntPatch_WLine) aWLine2 = new IntPatch_WLine(aL2S2, Standard_False);
1770
1771 Standard_Integer aWL1FindStatus = 0, aWL2FindStatus = 0;
1772
1773 Standard_Real anU1 = anUf;
1774
1775 Standard_Real aCriticalDelta[aNbCritPointsMax];
1776 for(Standard_Integer i = 0; i < aNbCritPointsMax; i++)
1777 aCriticalDelta[i] = anU1 - anU1crit[i];
1778
1779 Standard_Real aV11Prev = 0.0,
1780 aV12Prev = 0.0,
1781 aV21Prev = 0.0,
1782 aV22Prev = 0.0;
1783 Standard_Boolean isFirst = Standard_True;
1784
1785 while(anU1 <= anUl)
1786 {
1787 for(Standard_Integer i = 0; i < aNbCritPointsMax; i++)
1788 {
1789 if((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0)
1790 {
1791 anU1 = anU1crit[i];
1792 aWL1FindStatus = 2;
1793 aWL2FindStatus = 2;
1794 break;
1795 }
1796 }
1797
1798 Standard_Real anArg = anEquationCoeffs.mB * cos(anU1 - anEquationCoeffs.mFI1) + anEquationCoeffs.mC;
1799
1800 if(aNulValue > 1.0 - anArg)
1801 anArg = 1.0;
1802 if(anArg + 1.0 < aNulValue)
1803 anArg = -1.0;
1804
1805 Standard_Real aU21 = anEquationCoeffs.mFI2 + acos(anArg);
1806 InscribePoint(aUSurf2f, aUSurf2l, aU21, theTol2D, aPeriod);
1807
1808 Standard_Real aU22 = anEquationCoeffs.mFI2 - acos(anArg);
1809 InscribePoint(aUSurf2f, aUSurf2l, aU22, theTol2D, aPeriod);
1810
1811 const Standard_Real aV11 = anEquationCoeffs.mK21 * sin(aU21) +
1812 anEquationCoeffs.mK11 * sin(anU1) +
1813 anEquationCoeffs.mL21 * cos(aU21) +
1814 anEquationCoeffs.mL11 * cos(anU1) + anEquationCoeffs.mM1;
1815 const Standard_Real aV12 = anEquationCoeffs.mK21 * sin(aU22) +
1816 anEquationCoeffs.mK11 * sin(anU1) +
1817 anEquationCoeffs.mL21 * cos(aU22) +
1818 anEquationCoeffs.mL11 * cos(anU1) + anEquationCoeffs.mM1;
1819 const Standard_Real aV21 = anEquationCoeffs.mK22 * sin(aU21) +
1820 anEquationCoeffs.mK12 * sin(anU1) +
1821 anEquationCoeffs.mL22 * cos(aU21) +
1822 anEquationCoeffs.mL12 * cos(anU1) + anEquationCoeffs.mM2;
1823 const Standard_Real aV22 = anEquationCoeffs.mK22 * sin(aU22) +
1824 anEquationCoeffs.mK12 * sin(anU1) +
1825 anEquationCoeffs.mL22 * cos(aU22) +
1826 anEquationCoeffs.mL12 * cos(anU1) + anEquationCoeffs.mM2;
1827
1828 if(isFirst)
1829 {
1830 aV11Prev = aV11;
1831 aV12Prev = aV12;
1832 aV21Prev = aV21;
1833 aV22Prev = aV22;
1834 isFirst = Standard_False;
1835 }
1836
1837 if(((aUSurf2f-aU21) <= theTol2D) && ((aU21-aUSurf2l) <= theTol2D) && (aVSurf1f <= aV11) && (aV11 <= aVSurf1l) && (aVSurf2f <= aV21) && (aV21 <= aVSurf2l))
1838 {
1839 if(!aWL1FindStatus)
1840 {
1841 Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
1842
1843 AddBoundaryPoint(theQuad1, theQuad2, aWLine1, anEquationCoeffs, theUVSurf1, theUVSurf2,
1844 theTol2D, aPeriod, aNulValue, anU1, aU21, aV11, aV11Prev, aV21, aV21Prev, isTheReverse, 1.0, isFound1, isFound2);
1845
1846 if(isFound1 || isFound2)
1847 {
1848 aWL1FindStatus = 1;
1849 }
1850 }
1851
1852 if((aWL1FindStatus != 2) || (aWLine1->NbPnts() >= 1))
1853 {
1854 if(AddPointIntoWL(theQuad1, theQuad2, isTheReverse, gp_Pnt2d(anU1, aV11), gp_Pnt2d(aU21, aV21), aUSurf1f, aUSurf1l, aPeriod, aWLine1->Curve(), theTol2D))
1855 {
1856 if(!aWL1FindStatus)
1857 {
1858 aWL1FindStatus = 1;
1859 }
1860 }
1861 }
1862 }
1863 else
1864 {
1865 if(aWL1FindStatus == 1)
1866 {
1867 Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
1868
1869 AddBoundaryPoint(theQuad1, theQuad2, aWLine1, anEquationCoeffs, theUVSurf1, theUVSurf2,
1870 theTol2D, aPeriod, aNulValue, anU1, aU21, aV11, aV11Prev, aV21, aV21Prev, isTheReverse, 1.0, isFound1, isFound2);
1871
1872 if(isFound1 || isFound2)
1873 aWL1FindStatus = 2; //start a new line
1874 }
1875 }
1876
1877 if(((aUSurf2f-aU22) <= theTol2D) && ((aU22-aUSurf2l) <= theTol2D) && (aVSurf1f <= aV12) && (aV12 <= aVSurf1l)&& (aVSurf2f <= aV22) && (aV22 <= aVSurf2l))
1878 {
1879 if(!aWL2FindStatus)
1880 {
1881 Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
1882
1883 AddBoundaryPoint(theQuad1, theQuad2, aWLine2, anEquationCoeffs, theUVSurf1, theUVSurf2,
1884 theTol2D, aPeriod, aNulValue, anU1, aU22, aV12, aV12Prev, aV22, aV22Prev, isTheReverse, -1.0, isFound1, isFound2);
1885
1886 if(isFound1 || isFound2)
1887 {
1888 aWL2FindStatus = 1;
1889 }
1890 }
1891
1892 if((aWL2FindStatus != 2) || (aWLine2->NbPnts() >= 1))
1893 {
1894 if(AddPointIntoWL(theQuad1, theQuad2, isTheReverse, gp_Pnt2d(anU1, aV12), gp_Pnt2d(aU22, aV22), aUSurf1f, aUSurf1l, aPeriod, aWLine2->Curve(), theTol2D))
1895 {
1896 if(!aWL2FindStatus)
1897 {
1898 aWL2FindStatus = 1;
1899 }
1900 }
1901 }
1902 }
1903 else
1904 {
1905 if(aWL2FindStatus == 1)
1906 {
1907 Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
1908
1909 AddBoundaryPoint(theQuad1, theQuad2, aWLine2, anEquationCoeffs, theUVSurf1, theUVSurf2,
1910 theTol2D, aPeriod, aNulValue, anU1, aU22, aV12, aV12Prev, aV22, aV22Prev, isTheReverse, -1.0, isFound1, isFound2);
1911
1912 if(isFound1 || isFound2)
1913 aWL2FindStatus = 2; //start a new line
1914 }
1915 }
1916
1917 aV11Prev = aV11;
1918 aV12Prev = aV12;
1919 aV21Prev = aV21;
1920 aV22Prev = aV22;
1921
1922 if((aWL1FindStatus == 2) || (aWL2FindStatus == 2))
1923 {//current lines are filled. Go to the next lines
1924 anUf = anU1;
1925 break;
1926 }
1927
1928 Standard_Real aFact1 = !IsEqual(sin(aU21 - anEquationCoeffs.mFI2), 0.0) ?
1929 anEquationCoeffs.mK1 * sin(anU1 - anEquationCoeffs.mFIV1) +
1930 anEquationCoeffs.mL1 * anEquationCoeffs.mB * sin(aU21 - anEquationCoeffs.mPSIV1) *
1931 sin(anU1 - anEquationCoeffs.mFI1)/sin(aU21-anEquationCoeffs.mFI2) : 0.0,
1932 aFact2 = !IsEqual(sin(aU22-anEquationCoeffs.mFI2), 0.0) ?
1933 anEquationCoeffs.mK1 * sin(anU1 - anEquationCoeffs.mFIV1) +
1934 anEquationCoeffs.mL1 * anEquationCoeffs.mB * sin(aU22 - anEquationCoeffs.mPSIV1) *
1935 sin(anU1 - anEquationCoeffs.mFI1)/sin(aU22-anEquationCoeffs.mFI2) : 0.0;
1936
1937 Standard_Real aDeltaV1 = (aVSurf1l - aVSurf1f)/IntToReal(aNbPoints);
1938
1939 if((aV11 < aVSurf1f) && (aFact1 < 0.0))
1940 {//Make close to aVSurf1f by increasing anU1 (for the 1st line)
1941 aDeltaV1 = Min(aDeltaV1, Abs(aV11-aVSurf1f));
1942 }
1943
1944 if((aV12 < aVSurf1f) && (aFact2 < 0.0))
1945 {//Make close to aVSurf1f by increasing anU1 (for the 2nd line)
1946 aDeltaV1 = Min(aDeltaV1, Abs(aV12-aVSurf1f));
1947 }
1948
1949 if((aV11 > aVSurf1l) && (aFact1 > 0.0))
1950 {//Make close to aVSurf1l by increasing anU1 (for the 1st line)
1951 aDeltaV1 = Min(aDeltaV1, Abs(aV11-aVSurf1l));
1952 }
1953
1954 if((aV12 > aVSurf1l) && (aFact2 > 0.0))
1955 {//Make close to aVSurf1l by increasing anU1 (for the 1st line)
1956 aDeltaV1 = Min(aDeltaV1, Abs(aV12-aVSurf1l));
1957 }
1958
1959 Standard_Real aDeltaU1L1 = !IsEqual(aFact1,0.0)? Abs(aDeltaV1/aFact1) : aStepMax,
1960 aDeltaU1L2 = !IsEqual(aFact2,0.0)? Abs(aDeltaV1/aFact2) : aStepMax;
1961
1962 const Standard_Real aDeltaU1V1 = Min(aDeltaU1L1, aDeltaU1L2);
1963
1964 ///////////////////////////
1965 aFact1 = !IsEqual(sin(aU21-anEquationCoeffs.mFI2), 0.0) ?
1966 anEquationCoeffs.mK2 * sin(anU1 - anEquationCoeffs.mFIV2) +
1967 anEquationCoeffs.mL2 * anEquationCoeffs.mB * sin(aU21 - anEquationCoeffs.mPSIV2) *
1968 sin(anU1 - anEquationCoeffs.mFI1)/sin(aU21 - anEquationCoeffs.mFI2) : 0.0;
1969 aFact2 = !IsEqual(sin(aU22-anEquationCoeffs.mFI2), 0.0) ?
1970 anEquationCoeffs.mK2 * sin(anU1 - anEquationCoeffs.mFIV2) +
1971 anEquationCoeffs.mL2 * anEquationCoeffs.mB * sin(aU22 - anEquationCoeffs.mPSIV2) *
1972 sin(anU1 - anEquationCoeffs.mFI1)/sin(aU22 - anEquationCoeffs.mFI2) : 0.0;
1973
1974 Standard_Real aDeltaV2 = (aVSurf2l - aVSurf2f)/IntToReal(aNbPoints);
1975
1976 if((aV21 < aVSurf2f) && (aFact1 < 0.0))
1977 {//Make close to aVSurf2f by increasing anU1 (for the 1st line)
1978 aDeltaV2 = Min(aDeltaV2, Abs(aV21-aVSurf2f));
1979 }
1980
1981 if((aV22 < aVSurf2f) && (aFact2 < 0.0))
1982 {//Make close to aVSurf1f by increasing anU1 (for the 2nd line)
1983 aDeltaV2 = Min(aDeltaV2, Abs(aV22-aVSurf2f));
1984 }
1985
1986 if((aV21 > aVSurf2l) && (aFact1 > 0.0))
1987 {//Make close to aVSurf1l by increasing anU1 (for the 1st line)
1988 aDeltaV2 = Min(aDeltaV2, Abs(aV21-aVSurf2l));
1989 }
1990
1991 if((aV22 > aVSurf2l) && (aFact2 > 0.0))
1992 {//Make close to aVSurf1l by increasing anU1 (for the 1st line)
1993 aDeltaV2 = Min(aDeltaV2, Abs(aV22-aVSurf1l));
1994 }
1995
1996 aDeltaU1L1 = !IsEqual(aFact1,0.0)? Abs(aDeltaV2/aFact1) : aStepMax;
1997 aDeltaU1L2 = !IsEqual(aFact2,0.0)? Abs(aDeltaV2/aFact2) : aStepMax;
1998
1999 const Standard_Real aDeltaU1V2 = Min(aDeltaU1L1, aDeltaU1L2);
2000
2001 Standard_Real aDeltaU1 = Min(aDeltaU1V1, aDeltaU1V2);
2002
2003 if(aDeltaU1 < aStepMin)
2004 aDeltaU1 = aStepMin;
2005
2006 if(aDeltaU1 > aStepMax)
2007 aDeltaU1 = aStepMax;
2008
2009 anU1 += aDeltaU1;
2010
2011 const Standard_Real aDiff = anU1 - anUl;
2012 if((0.0 < aDiff) && (aDiff < aDeltaU1-Precision::PConfusion()))
2013 anU1 = anUl;
2014
2015 anUf = anU1;
2016
2017 if(aWLine1->NbPnts() != 1)
2018 isAddedIntoWL1 = Standard_False;
2019
2020 if(aWLine2->NbPnts() != 1)
2021 isAddedIntoWL2 = Standard_False;
2022 }
2023
2024 if((aWLine1->NbPnts() == 1) && (!isAddedIntoWL1))
2025 {
2026 isTheEmpty = Standard_False;
2027 Standard_Real u1, v1, u2, v2;
2028 aWLine1->Point(1).Parameters(u1, v1, u2, v2);
2029 IntPatch_Point aP;
2030 aP.SetParameter(u1);
2031 aP.SetParameters(u1, v1, u2, v2);
2032 aP.SetTolerance(theTol3D);
2033 aP.SetValue(aWLine1->Point(1).Value());
2034
2035 theSPnt.Append(aP);
2036 }
2037 else if(aWLine1->NbPnts() > 1)
2038 {
2039 isTheEmpty = Standard_False;
2040 isAddedIntoWL1 = Standard_True;
2041
2042 SeekAdditionalPoints(theQuad1, theQuad2, aWLine1->Curve(), anEquationCoeffs, aNbPoints, aUSurf2f, aUSurf2l, theTol2D, aPeriod, 1.0, isTheReverse);
2043
2044 aWLine1->ComputeVertexParameters(theTol3D);
2045 theSlin.Append(aWLine1);
2046 }
2047 else
2048 {
2049 isAddedIntoWL1 = Standard_False;
2050 }
2051
2052 if((aWLine2->NbPnts() == 1) && (!isAddedIntoWL2))
2053 {
2054 isTheEmpty = Standard_False;
2055 Standard_Real u1, v1, u2, v2;
2056 aWLine2->Point(1).Parameters(u1, v1, u2, v2);
2057 IntPatch_Point aP;
2058 aP.SetParameter(u1);
2059 aP.SetParameters(u1, v1, u2, v2);
2060 aP.SetTolerance(theTol3D);
2061 aP.SetValue(aWLine2->Point(1).Value());
2062
2063 theSPnt.Append(aP);
2064 }
2065 else if(aWLine2->NbPnts() > 1)
2066 {
2067 isTheEmpty = Standard_False;
2068 isAddedIntoWL2 = Standard_True;
2069
2070 SeekAdditionalPoints(theQuad1, theQuad2, aWLine2->Curve(), anEquationCoeffs, aNbPoints, aUSurf2f, aUSurf2l, theTol2D, aPeriod, -1.0, isTheReverse);
2071
2072 aWLine2->ComputeVertexParameters(theTol3D);
2073 theSlin.Append(aWLine2);
2074 }
2075 else
2076 {
2077 isAddedIntoWL2 = Standard_False;
2078 }
2079 }
2080 }
2081
2082 return Standard_True;
2083}
7fd59977 2084
2085//=======================================================================
2086//function : IntCySp
2087//purpose :
2088//=======================================================================
2089Standard_Boolean IntCySp(const IntSurf_Quadric& Quad1,
2090 const IntSurf_Quadric& Quad2,
2091 const Standard_Real Tol,
2092 const Standard_Boolean Reversed,
2093 Standard_Boolean& Empty,
2094 Standard_Boolean& Multpoint,
2095 IntPatch_SequenceOfLine& slin,
2096 IntPatch_SequenceOfPoint& spnt)
2097
2098{
2099 Standard_Integer i;
2100
2101 IntSurf_TypeTrans trans1,trans2;
2102 IntAna_ResultType typint;
2103 IntPatch_Point ptsol;
2104 gp_Circ cirsol;
2105
2106 gp_Cylinder Cy;
2107 gp_Sphere Sp;
2108
2109 if (!Reversed) {
2110 Cy = Quad1.Cylinder();
2111 Sp = Quad2.Sphere();
2112 }
2113 else {
2114 Cy = Quad2.Cylinder();
2115 Sp = Quad1.Sphere();
2116 }
2117 IntAna_QuadQuadGeo inter(Cy,Sp,Tol);
2118
2119 if (!inter.IsDone()) {return Standard_False;}
2120
2121 typint = inter.TypeInter();
2122 Standard_Integer NbSol = inter.NbSolutions();
2123 Empty = Standard_False;
2124
2125 switch (typint) {
2126
2127 case IntAna_Empty :
2128 {
2129 Empty = Standard_True;
2130 }
2131 break;
2132
2133 case IntAna_Point :
2134 {
2135 gp_Pnt psol(inter.Point(1));
2136 Standard_Real U1,V1,U2,V2;
2137 Quad1.Parameters(psol,U1,V1);
2138 Quad2.Parameters(psol,U2,V2);
2139 ptsol.SetValue(psol,Tol,Standard_True);
2140 ptsol.SetParameters(U1,V1,U2,V2);
2141 spnt.Append(ptsol);
2142 }
2143 break;
2144
2145 case IntAna_Circle:
2146 {
2147 cirsol = inter.Circle(1);
2148 gp_Vec Tgt;
2149 gp_Pnt ptref;
2150 ElCLib::D1(0.,cirsol,ptref,Tgt);
2151
2152 if (NbSol == 1) {
2153 gp_Vec TestCurvature(ptref,Sp.Location());
2154 gp_Vec Normsp,Normcyl;
2155 if (!Reversed) {
2156 Normcyl = Quad1.Normale(ptref);
2157 Normsp = Quad2.Normale(ptref);
2158 }
2159 else {
2160 Normcyl = Quad2.Normale(ptref);
2161 Normsp = Quad1.Normale(ptref);
2162 }
2163
2164 IntSurf_Situation situcyl;
2165 IntSurf_Situation situsp;
2166
2167 if (Normcyl.Dot(TestCurvature) > 0.) {
2168 situsp = IntSurf_Outside;
2169 if (Normsp.Dot(Normcyl) > 0.) {
2170 situcyl = IntSurf_Inside;
2171 }
2172 else {
2173 situcyl = IntSurf_Outside;
2174 }
2175 }
2176 else {
2177 situsp = IntSurf_Inside;
2178 if (Normsp.Dot(Normcyl) > 0.) {
2179 situcyl = IntSurf_Outside;
2180 }
2181 else {
2182 situcyl = IntSurf_Inside;
2183 }
2184 }
2185 Handle(IntPatch_GLine) glig;
2186 if (!Reversed) {
2187 glig = new IntPatch_GLine(cirsol, Standard_True, situcyl, situsp);
2188 }
2189 else {
2190 glig = new IntPatch_GLine(cirsol, Standard_True, situsp, situcyl);
2191 }
2192 slin.Append(glig);
2193 }
2194 else {
2195 if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
2196 trans1 = IntSurf_Out;
2197 trans2 = IntSurf_In;
2198 }
2199 else {
2200 trans1 = IntSurf_In;
2201 trans2 = IntSurf_Out;
2202 }
2203 Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
2204 slin.Append(glig);
2205
2206 cirsol = inter.Circle(2);
2207 ElCLib::D1(0.,cirsol,ptref,Tgt);
2208 Standard_Real qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
2209 if(qwe> 0.0000001) {
2210 trans1 = IntSurf_Out;
2211 trans2 = IntSurf_In;
2212 }
2213 else if(qwe<-0.0000001) {
2214 trans1 = IntSurf_In;
2215 trans2 = IntSurf_Out;
2216 }
2217 else {
2218 trans1=trans2=IntSurf_Undecided;
2219 }
2220 glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
2221 slin.Append(glig);
2222 }
2223 }
2224 break;
2225
2226 case IntAna_NoGeometricSolution:
2227 {
2228 gp_Pnt psol;
2229 Standard_Real U1,V1,U2,V2;
2230 IntAna_IntQuadQuad anaint(Cy,Sp,Tol);
2231 if (!anaint.IsDone()) {
2232 return Standard_False;
2233 }
2234
2235 if (anaint.NbPnt()==0 && anaint.NbCurve()==0) {
2236 Empty = Standard_True;
2237 }
2238 else {
2239
2240 NbSol = anaint.NbPnt();
2241 for (i = 1; i <= NbSol; i++) {
2242 psol = anaint.Point(i);
2243 Quad1.Parameters(psol,U1,V1);
2244 Quad2.Parameters(psol,U2,V2);
2245 ptsol.SetValue(psol,Tol,Standard_True);
2246 ptsol.SetParameters(U1,V1,U2,V2);
2247 spnt.Append(ptsol);
2248 }
2249
2250 gp_Pnt ptvalid,ptf,ptl;
2251 gp_Vec tgvalid;
2252 Standard_Real first,last,para;
2253 IntAna_Curve curvsol;
2254 Standard_Boolean tgfound;
2255 Standard_Integer kount;
2256
2257 NbSol = anaint.NbCurve();
2258 for (i = 1; i <= NbSol; i++) {
2259 curvsol = anaint.Curve(i);
2260 curvsol.Domain(first,last);
2261 ptf = curvsol.Value(first);
2262 ptl = curvsol.Value(last);
2263
2264 para = last;
2265 kount = 1;
2266 tgfound = Standard_False;
2267
2268 while (!tgfound) {
2269 para = (1.123*first + para)/2.123;
2270 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
2271 if (!tgfound) {
2272 kount ++;
2273 tgfound = kount > 5;
2274 }
2275 }
2276 Handle(IntPatch_ALine) alig;
2277 if (kount <= 5) {
2278 Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
2279 Quad1.Normale(ptvalid));
2280 if(qwe> 0.00000001) {
2281 trans1 = IntSurf_Out;
2282 trans2 = IntSurf_In;
2283 }
2284 else if(qwe<-0.00000001) {
2285 trans1 = IntSurf_In;
2286 trans2 = IntSurf_Out;
2287 }
2288 else {
2289 trans1=trans2=IntSurf_Undecided;
2290 }
2291 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
2292 }
2293 else {
2294 alig = new IntPatch_ALine(curvsol,Standard_False);
2295 }
2296 Standard_Boolean TempFalse1a = Standard_False;
2297 Standard_Boolean TempFalse2a = Standard_False;
2298
2299 //-- ptf et ptl : points debut et fin de alig
2300
2301 ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1a,ptf,first,
2302 TempFalse2a,ptl,last,Multpoint,Tol);
2303 slin.Append(alig);
2304 } //-- boucle sur les lignes
2305 } //-- solution avec au moins une lihne
2306 }
2307 break;
2308
2309 default:
2310 {
2311 return Standard_False;
2312 }
2313 }
2314 return Standard_True;
2315}
2316//=======================================================================
2317//function : IntCyCo
2318//purpose :
2319//=======================================================================
2320Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
2321 const IntSurf_Quadric& Quad2,
2322 const Standard_Real Tol,
2323 const Standard_Boolean Reversed,
2324 Standard_Boolean& Empty,
2325 Standard_Boolean& Multpoint,
2326 IntPatch_SequenceOfLine& slin,
2327 IntPatch_SequenceOfPoint& spnt)
2328
2329{
2330 IntPatch_Point ptsol;
2331
2332 Standard_Integer i;
2333
2334 IntSurf_TypeTrans trans1,trans2;
2335 IntAna_ResultType typint;
2336 gp_Circ cirsol;
2337
2338 gp_Cylinder Cy;
2339 gp_Cone Co;
2340
2341 if (!Reversed) {
2342 Cy = Quad1.Cylinder();
2343 Co = Quad2.Cone();
2344 }
2345 else {
2346 Cy = Quad2.Cylinder();
2347 Co = Quad1.Cone();
2348 }
2349 IntAna_QuadQuadGeo inter(Cy,Co,Tol);
2350
2351 if (!inter.IsDone()) {return Standard_False;}
2352
2353 typint = inter.TypeInter();
2354 Standard_Integer NbSol = inter.NbSolutions();
2355 Empty = Standard_False;
2356
2357 switch (typint) {
2358
2359 case IntAna_Empty : {
2360 Empty = Standard_True;
2361 }
2362 break;
2363
2364 case IntAna_Point :{
2365 gp_Pnt psol(inter.Point(1));
2366 Standard_Real U1,V1,U2,V2;
2367 Quad1.Parameters(psol,U1,V1);
2368 Quad1.Parameters(psol,U2,V2);
2369 ptsol.SetValue(psol,Tol,Standard_True);
2370 ptsol.SetParameters(U1,V1,U2,V2);
2371 spnt.Append(ptsol);
2372 }
2373 break;
2374
2375 case IntAna_Circle: {
2376 gp_Vec Tgt;
2377 gp_Pnt ptref;
2378 Standard_Integer j;
2379 Standard_Real qwe;
2380 //
2381 for(j=1; j<=2; ++j) {
2382 cirsol = inter.Circle(j);
2383 ElCLib::D1(0.,cirsol,ptref,Tgt);
2384 qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
2385 if(qwe> 0.00000001) {
2386 trans1 = IntSurf_Out;
2387 trans2 = IntSurf_In;
2388 }
2389 else if(qwe<-0.00000001) {
2390 trans1 = IntSurf_In;
2391 trans2 = IntSurf_Out;
2392 }
2393 else {
2394 trans1=trans2=IntSurf_Undecided;
2395 }
2396 Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
2397 slin.Append(glig);
2398 }
2399 }
2400 break;
2401
2402 case IntAna_NoGeometricSolution: {
2403 gp_Pnt psol;
2404 Standard_Real U1,V1,U2,V2;
2405 IntAna_IntQuadQuad anaint(Cy,Co,Tol);
2406 if (!anaint.IsDone()) {
2407 return Standard_False;
2408 }
2409
2410 if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) {
2411 Empty = Standard_True;
2412 }
2413 else {
2414 NbSol = anaint.NbPnt();
2415 for (i = 1; i <= NbSol; i++) {
2416 psol = anaint.Point(i);
2417 Quad1.Parameters(psol,U1,V1);
2418 Quad2.Parameters(psol,U2,V2);
2419 ptsol.SetValue(psol,Tol,Standard_True);
2420 ptsol.SetParameters(U1,V1,U2,V2);
2421 spnt.Append(ptsol);
2422 }
2423
2424 gp_Pnt ptvalid, ptf, ptl;
2425 gp_Vec tgvalid;
2426 //
2427 Standard_Real first,last,para;
2428 Standard_Boolean tgfound,firstp,lastp,kept;
2429 Standard_Integer kount;
2430 //
2431 //
2432 //IntAna_Curve curvsol;
2433 IntAna_Curve aC;
2434 IntAna_ListOfCurve aLC;
2435 IntAna_ListIteratorOfListOfCurve aIt;
7fd59977 2436
2437 //
2438 NbSol = anaint.NbCurve();
2439 for (i = 1; i <= NbSol; ++i) {
2440 kept = Standard_False;
2441 //curvsol = anaint.Curve(i);
2442 aC=anaint.Curve(i);
2443 aLC.Clear();
96a95605 2444 ExploreCurve(Cy, Co, aC, 10.*Tol, aLC);
7fd59977 2445 //
2446 aIt.Initialize(aLC);
2447 for (; aIt.More(); aIt.Next()) {
2448 IntAna_Curve& curvsol=aIt.Value();
2449 //
2450 curvsol.Domain(first, last);
2451 firstp = !curvsol.IsFirstOpen();
2452 lastp = !curvsol.IsLastOpen();
2453 if (firstp) {
2454 ptf = curvsol.Value(first);
2455 }
2456 if (lastp) {
2457 ptl = curvsol.Value(last);
2458 }
2459 para = last;
2460 kount = 1;
2461 tgfound = Standard_False;
2462
2463 while (!tgfound) {
2464 para = (1.123*first + para)/2.123;
2465 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
2466 if (!tgfound) {
2467 kount ++;
2468 tgfound = kount > 5;
2469 }
2470 }
2471 Handle(IntPatch_ALine) alig;
2472 if (kount <= 5) {
2473 Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
2474 Quad1.Normale(ptvalid));
2475 if(qwe> 0.00000001) {
2476 trans1 = IntSurf_Out;
2477 trans2 = IntSurf_In;
2478 }
2479 else if(qwe<-0.00000001) {
2480 trans1 = IntSurf_In;
2481 trans2 = IntSurf_Out;
2482 }
2483 else {
2484 trans1=trans2=IntSurf_Undecided;
2485 }
2486 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
2487 kept = Standard_True;
2488 }
2489 else {
2490 ptvalid = curvsol.Value(para);
2491 alig = new IntPatch_ALine(curvsol,Standard_False);
2492 kept = Standard_True;
2493 //-- cout << "Transition indeterminee" << endl;
2494 }
2495 if (kept) {
2496 Standard_Boolean Nfirstp = !firstp;
2497 Standard_Boolean Nlastp = !lastp;
2498 ProcessBounds(alig,slin,Quad1,Quad2,Nfirstp,ptf,first,
2499 Nlastp,ptl,last,Multpoint,Tol);
2500 slin.Append(alig);
2501 }
2502 } // for (; aIt.More(); aIt.Next())
2503 } // for (i = 1; i <= NbSol; ++i)
2504 }
2505 }
2506 break;
2507
2508 default:
2509 return Standard_False;
2510
2511 } // switch (typint)
2512
2513 return Standard_True;
2514}
2515//=======================================================================
2516//function : ExploreCurve
2517//purpose :
2518//=======================================================================
2519Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy,
2520 const gp_Cone& aCo,
2521 IntAna_Curve& aC,
2522 const Standard_Real aTol,
2523 IntAna_ListOfCurve& aLC)
2524
2525{
2526 Standard_Boolean bFind=Standard_False;
2527 Standard_Real aTheta, aT1, aT2, aDst;
2528 gp_Pnt aPapx, aPx;
2529 //
2530 //aC.Dump();
2531 //
2532 aLC.Clear();
2533 aLC.Append(aC);
2534 //
2535 aPapx=aCo.Apex();
2536 //
2537 aC.Domain(aT1, aT2);
2538 //
2539 aPx=aC.Value(aT1);
2540 aDst=aPx.Distance(aPapx);
2541 if (aDst<aTol) {
2542 return bFind;
2543 }
2544 aPx=aC.Value(aT2);
2545 aDst=aPx.Distance(aPapx);
2546 if (aDst<aTol) {
2547 return bFind;
2548 }
2549 //
2550 bFind=aC.FindParameter(aPapx, aTheta);
2551 if (!bFind){
2552 return bFind;
2553 }
2554 //
2555 aPx=aC.Value(aTheta);
2556 aDst=aPx.Distance(aPapx);
2557 if (aDst>aTol) {
2558 return !bFind;
2559 }
2560 //
2561 // need to be splitted at aTheta
2562 IntAna_Curve aC1, aC2;
2563 //
2564 aC1=aC;
2565 aC1.SetDomain(aT1, aTheta);
2566 aC2=aC;
2567 aC2.SetDomain(aTheta, aT2);
2568 //
2569 aLC.Clear();
2570 aLC.Append(aC1);
2571 aLC.Append(aC2);
2572 //
2573 return bFind;
2574}
2575//=======================================================================
2576//function : IsToReverse
2577//purpose :
2578//=======================================================================
2579Standard_Boolean IsToReverse(const gp_Cylinder& Cy1,
2580 const gp_Cylinder& Cy2,
2581 const Standard_Real Tol)
2582{
2583 Standard_Boolean bRet;
2584 Standard_Real aR1,aR2, dR, aSc1, aSc2;
2585 //
2586 bRet=Standard_False;
2587 //
2588 aR1=Cy1.Radius();
2589 aR2=Cy2.Radius();
2590 dR=aR1-aR2;
2591 if (dR<0.) {
2592 dR=-dR;
2593 }
2594 if (dR>Tol) {
2595 bRet=aR1>aR2;
2596 }
2597 else {
2598 gp_Dir aDZ(0.,0.,1.);
2599 //
2600 const gp_Dir& aDir1=Cy1.Axis().Direction();
2601 aSc1=aDir1*aDZ;
2602 if (aSc1<0.) {
2603 aSc1=-aSc1;
2604 }
2605 //
2606 const gp_Dir& aDir2=Cy2.Axis().Direction();
2607 aSc2=aDir2*aDZ;
2608 if (aSc2<0.) {
2609 aSc2=-aSc2;
2610 }
2611 //
2612 if (aSc2<aSc1) {
2613 bRet=!bRet;
2614 }
2615 }//if (dR<Tol) {
2616 return bRet;
2617}
Open CASCADE Technology repository