projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0023178: Intersection of cylinders fails to produce results
[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
7365fad6 17#include <algorithm>
d30895f5 18#include <Bnd_Range.hxx>
7fd59977 19#include <IntAna_ListIteratorOfListOfCurve.hxx>
d30895f5 20#include <IntAna_ListOfCurve.hxx>
ecc4f148 21#include <IntPatch_WLine.hxx>
d30895f5 22#include <Standard_DivideByZero.hxx>
b5ef9d91 23#include <math_Matrix.hxx>
24
25//If Abs(a) <= aNulValue then it is considered that a = 0.
26static const Standard_Real aNulValue = 1.0e-11;
27
d30895f5 28static void ShortCosForm( const Standard_Real theCosFactor,
29 const Standard_Real theSinFactor,
30 Standard_Real& theCoeff,
31 Standard_Real& theAngle);
fa0291ff 32//
7fd59977 33static
34 Standard_Boolean ExploreCurve(const gp_Cylinder& aCy,
35 const gp_Cone& aCo,
36 IntAna_Curve& aC,
37 const Standard_Real aTol,
38 IntAna_ListOfCurve& aLC);
7fd59977 39
b5ef9d91 40static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
41 const Standard_Real theUlTarget,
42 Standard_Real& theUGiven,
43 const Standard_Real theTol2D,
44 const Standard_Real thePeriod,
45 const Standard_Boolean theFlForce);
46
d30895f5 47
48class ComputationMethods
49{
50public:
51 //Stores equations coefficients
52 struct stCoeffsValue
53 {
54 stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&);
55
56 math_Vector mVecA1;
57 math_Vector mVecA2;
58 math_Vector mVecB1;
59 math_Vector mVecB2;
60 math_Vector mVecC1;
61 math_Vector mVecC2;
62 math_Vector mVecD;
63
64 Standard_Real mK21; //sinU2
65 Standard_Real mK11; //sinU1
66 Standard_Real mL21; //cosU2
67 Standard_Real mL11; //cosU1
68 Standard_Real mM1; //Free member
69
70 Standard_Real mK22; //sinU2
71 Standard_Real mK12; //sinU1
72 Standard_Real mL22; //cosU2
73 Standard_Real mL12; //cosU1
74 Standard_Real mM2; //Free member
75
76 Standard_Real mK1;
77 Standard_Real mL1;
78 Standard_Real mK2;
79 Standard_Real mL2;
80
81 Standard_Real mFIV1;
82 Standard_Real mPSIV1;
83 Standard_Real mFIV2;
84 Standard_Real mPSIV2;
85
86 Standard_Real mB;
87 Standard_Real mC;
88 Standard_Real mFI1;
89 Standard_Real mFI2;
90 };
91
92
93 //! Determines, if U2(U1) function is increasing.
94 static Standard_Boolean CylCylMonotonicity(const Standard_Real theU1par,
95 const Standard_Integer theWLIndex,
96 const stCoeffsValue& theCoeffs,
97 const Standard_Real thePeriod,
98 Standard_Boolean& theIsIncreasing);
99
100 //! Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0,
101 //! esimates the tolerance of U2-computing (estimation result is
102 //! assigned to *theDelta value).
103 static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
104 const Standard_Integer theWLIndex,
105 const stCoeffsValue& theCoeffs,
106 Standard_Real& theU2,
107 Standard_Real* const theDelta = 0);
108
109 static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1,
110 const Standard_Real theU2,
111 const stCoeffsValue& theCoeffs,
112 Standard_Real& theV1,
113 Standard_Real& theV2);
114
115 static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
116 const Standard_Integer theWLIndex,
117 const stCoeffsValue& theCoeffs,
118 Standard_Real& theU2,
119 Standard_Real& theV1,
120 Standard_Real& theV2);
121
122};
123
124ComputationMethods::stCoeffsValue::stCoeffsValue(const gp_Cylinder& theCyl1,
125 const gp_Cylinder& theCyl2):
126 mVecA1(-theCyl1.Radius()*theCyl1.XAxis().Direction().XYZ()),
127 mVecA2(theCyl2.Radius()*theCyl2.XAxis().Direction().XYZ()),
128 mVecB1(-theCyl1.Radius()*theCyl1.YAxis().Direction().XYZ()),
129 mVecB2(theCyl2.Radius()*theCyl2.YAxis().Direction().XYZ()),
130 mVecC1(theCyl1.Axis().Direction().XYZ()),
131 mVecC2(theCyl2.Axis().Direction().XYZ().Reversed()),
132 mVecD(theCyl2.Location().XYZ() - theCyl1.Location().XYZ())
133{
134 enum CoupleOfEquation
135 {
136 COENONE = 0,
137 COE12 = 1,
138 COE23 = 2,
139 COE13 = 3
140 }aFoundCouple = COENONE;
141
142
143 Standard_Real aDetV1V2 = 0.0;
144
145 const Standard_Real aDelta1 = mVecC1(1)*mVecC2(2)-mVecC1(2)*mVecC2(1); //1-2
146 const Standard_Real aDelta2 = mVecC1(2)*mVecC2(3)-mVecC1(3)*mVecC2(2); //2-3
147 const Standard_Real aDelta3 = mVecC1(1)*mVecC2(3)-mVecC1(3)*mVecC2(1); //1-3
148 const Standard_Real anAbsD1 = Abs(aDelta1); //1-2
149 const Standard_Real anAbsD2 = Abs(aDelta2); //2-3
150 const Standard_Real anAbsD3 = Abs(aDelta3); //1-3
151
152 if(anAbsD1 >= anAbsD2)
153 {
154 if(anAbsD3 > anAbsD1)
155 {
156 aFoundCouple = COE13;
157 aDetV1V2 = aDelta3;
158 }
159 else
160 {
161 aFoundCouple = COE12;
162 aDetV1V2 = aDelta1;
163 }
164 }
165 else
166 {
167 if(anAbsD3 > anAbsD2)
168 {
169 aFoundCouple = COE13;
170 aDetV1V2 = aDelta3;
171 }
172 else
173 {
174 aFoundCouple = COE23;
175 aDetV1V2 = aDelta2;
176 }
177 }
178
179 // In point of fact, every determinant (aDelta1, aDelta2 and aDelta3) is
180 // cross-product between directions (i.e. sine of angle).
181 // If sine is too small then sine is (approx.) equal to angle itself.
182 // Therefore, in this case we should compare sine with angular tolerance.
183 // This constant is used for check if axes are parallel (see constructor
184 // AxeOperator::AxeOperator(...) in IntAna_QuadQuadGeo.cxx file).
185 if(Abs(aDetV1V2) < Precision::Angular())
186 {
187 Standard_Failure::Raise("Error. Exception in divide by zerro (IntCyCyTrim)!!!!");
188 }
189
190 switch(aFoundCouple)
191 {
192 case COE12:
193 break;
194 case COE23:
195 {
196 math_Vector aVTemp(mVecA1);
197 mVecA1(1) = aVTemp(2);
198 mVecA1(2) = aVTemp(3);
199 mVecA1(3) = aVTemp(1);
200
201 aVTemp = mVecA2;
202 mVecA2(1) = aVTemp(2);
203 mVecA2(2) = aVTemp(3);
204 mVecA2(3) = aVTemp(1);
205
206 aVTemp = mVecB1;
207 mVecB1(1) = aVTemp(2);
208 mVecB1(2) = aVTemp(3);
209 mVecB1(3) = aVTemp(1);
210
211 aVTemp = mVecB2;
212 mVecB2(1) = aVTemp(2);
213 mVecB2(2) = aVTemp(3);
214 mVecB2(3) = aVTemp(1);
215
216 aVTemp = mVecC1;
217 mVecC1(1) = aVTemp(2);
218 mVecC1(2) = aVTemp(3);
219 mVecC1(3) = aVTemp(1);
220
221 aVTemp = mVecC2;
222 mVecC2(1) = aVTemp(2);
223 mVecC2(2) = aVTemp(3);
224 mVecC2(3) = aVTemp(1);
225
226 aVTemp = mVecD;
227 mVecD(1) = aVTemp(2);
228 mVecD(2) = aVTemp(3);
229 mVecD(3) = aVTemp(1);
230
231 break;
232 }
233 case COE13:
234 {
235 math_Vector aVTemp = mVecA1;
236 mVecA1(2) = aVTemp(3);
237 mVecA1(3) = aVTemp(2);
238
239 aVTemp = mVecA2;
240 mVecA2(2) = aVTemp(3);
241 mVecA2(3) = aVTemp(2);
242
243 aVTemp = mVecB1;
244 mVecB1(2) = aVTemp(3);
245 mVecB1(3) = aVTemp(2);
246
247 aVTemp = mVecB2;
248 mVecB2(2) = aVTemp(3);
249 mVecB2(3) = aVTemp(2);
250
251 aVTemp = mVecC1;
252 mVecC1(2) = aVTemp(3);
253 mVecC1(3) = aVTemp(2);
254
255 aVTemp = mVecC2;
256 mVecC2(2) = aVTemp(3);
257 mVecC2(3) = aVTemp(2);
258
259 aVTemp = mVecD;
260 mVecD(2) = aVTemp(3);
261 mVecD(3) = aVTemp(2);
262
263 break;
264 }
265 default:
266 break;
267 }
268
269 //------- For V1 (begin)
270 //sinU2
271 mK21 = (mVecC2(2)*mVecB2(1)-mVecC2(1)*mVecB2(2))/aDetV1V2;
272 //sinU1
273 mK11 = (mVecC2(2)*mVecB1(1)-mVecC2(1)*mVecB1(2))/aDetV1V2;
274 //cosU2
275 mL21 = (mVecC2(2)*mVecA2(1)-mVecC2(1)*mVecA2(2))/aDetV1V2;
276 //cosU1
277 mL11 = (mVecC2(2)*mVecA1(1)-mVecC2(1)*mVecA1(2))/aDetV1V2;
278 //Free member
279 mM1 = (mVecC2(2)*mVecD(1)-mVecC2(1)*mVecD(2))/aDetV1V2;
280 //------- For V1 (end)
281
282 //------- For V2 (begin)
283 //sinU2
284 mK22 = (mVecC1(1)*mVecB2(2)-mVecC1(2)*mVecB2(1))/aDetV1V2;
285 //sinU1
286 mK12 = (mVecC1(1)*mVecB1(2)-mVecC1(2)*mVecB1(1))/aDetV1V2;
287 //cosU2
288 mL22 = (mVecC1(1)*mVecA2(2)-mVecC1(2)*mVecA2(1))/aDetV1V2;
289 //cosU1
290 mL12 = (mVecC1(1)*mVecA1(2)-mVecC1(2)*mVecA1(1))/aDetV1V2;
291 //Free member
292 mM2 = (mVecC1(1)*mVecD(2)-mVecC1(2)*mVecD(1))/aDetV1V2;
293 //------- For V1 (end)
294
295 ShortCosForm(mL11, mK11, mK1, mFIV1);
296 ShortCosForm(mL21, mK21, mL1, mPSIV1);
297 ShortCosForm(mL12, mK12, mK2, mFIV2);
298 ShortCosForm(mL22, mK22, mL2, mPSIV2);
299
300 const Standard_Real aA1=mVecC1(3)*mK21+mVecC2(3)*mK22-mVecB2(3), //sinU2
301 aA2=mVecC1(3)*mL21+mVecC2(3)*mL22-mVecA2(3), //cosU2
302 aB1=mVecB1(3)-mVecC1(3)*mK11-mVecC2(3)*mK12, //sinU1
303 aB2=mVecA1(3)-mVecC1(3)*mL11-mVecC2(3)*mL12; //cosU1
304
305 mC =mVecD(3) - mVecC1(3)*mM1 -mVecC2(3)*mM2; //Free
306
307 Standard_Real aA = 0.0;
308
309 ShortCosForm(aB2,aB1,mB,mFI1);
310 ShortCosForm(aA2,aA1,aA,mFI2);
311
312 mB /= aA;
313 mC /= aA;
314}
315
316class WorkWithBoundaries
317{
318public:
319 enum SearchBoundType
320 {
321 SearchNONE = 0,
322 SearchV1 = 1,
323 SearchV2 = 2
324 };
325
326 struct StPInfo
327 {
328 StPInfo()
329 {
330 mySurfID = 0;
331 myU1 = RealLast();
332 myV1 = RealLast();
333 myU2 = RealLast();
334 myV2 = RealLast();
335 }
336
337 //Equal to 0 for 1st surface non-zerro for 2nd one.
338 Standard_Integer mySurfID;
339
340 Standard_Real myU1;
341 Standard_Real myV1;
342 Standard_Real myU2;
343 Standard_Real myV2;
344
345 bool operator>(const StPInfo& theOther) const
346 {
347 return myU1 > theOther.myU1;
348 }
349
350 bool operator<(const StPInfo& theOther) const
351 {
352 return myU1 < theOther.myU1;
353 }
354
355 bool operator==(const StPInfo& theOther) const
356 {
357 return myU1 == theOther.myU1;
358 }
359 };
360
361 WorkWithBoundaries(const IntSurf_Quadric& theQuad1,
362 const IntSurf_Quadric& theQuad2,
363 const ComputationMethods::stCoeffsValue& theCoeffs,
364 const Bnd_Box2d& theUVSurf1,
365 const Bnd_Box2d& theUVSurf2,
366 const Standard_Integer theNbWLines,
367 const Standard_Real thePeriod,
368 const Standard_Real theTol3D,
369 const Standard_Real theTol2D,
370 const Standard_Boolean isTheReverse) :
371 myQuad1(theQuad1), myQuad2(theQuad2), myCoeffs(theCoeffs),
372 myUVSurf1(theUVSurf1), myUVSurf2(theUVSurf2), myNbWLines(theNbWLines),
373 myPeriod(thePeriod), myTol3D(theTol3D), myTol2D(theTol2D),
374 myIsReverse(isTheReverse)
375 {
376 };
377
378 void AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL,
379 const Standard_Real theU1,
380 const Standard_Real theU2,
381 const Standard_Real theV1,
382 const Standard_Real theV1Prev,
383 const Standard_Real theV2,
384 const Standard_Real theV2Prev,
385 const Standard_Integer theWLIndex,
386 const Standard_Boolean theFlForce,
387 Standard_Boolean& isTheFound1,
388 Standard_Boolean& isTheFound2) const;
389
390 Standard_Boolean BoundariesComputing(Standard_Real theU1f[],
391 Standard_Real theU1l[]) const;
392
393 void BoundaryEstimation(const gp_Cylinder& theCy1,
394 const gp_Cylinder& theCy2,
395 Bnd_Range& theOutBoxS1,
396 Bnd_Range& theOutBoxS2) const;
397
398protected:
399
400 Standard_Boolean
401 SearchOnVBounds(const SearchBoundType theSBType,
402 const Standard_Real theVzad,
403 const Standard_Real theVInit,
404 const Standard_Real theInitU2,
405 const Standard_Real theInitMainVar,
406 Standard_Real& theMainVariableValue) const;
407
408 const WorkWithBoundaries& operator=(const WorkWithBoundaries&);
409
410private:
411 friend class ComputationMethods;
412
413 const IntSurf_Quadric& myQuad1;
414 const IntSurf_Quadric& myQuad2;
415 const ComputationMethods::stCoeffsValue& myCoeffs;
416 const Bnd_Box2d& myUVSurf1;
417 const Bnd_Box2d& myUVSurf2;
418 const Standard_Integer myNbWLines;
419 const Standard_Real myPeriod;
420 const Standard_Real myTol3D;
421 const Standard_Real myTol2D;
422 const Standard_Boolean myIsReverse;
423};
424
425static
426 Standard_Boolean ExploreCurve(const gp_Cylinder& aCy,
427 const gp_Cone& aCo,
428 IntAna_Curve& aC,
429 const Standard_Real aTol,
430 IntAna_ListOfCurve& aLC);
431
b5ef9d91 432static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
433 const IntSurf_Quadric& theQuad2,
434 const Handle(IntSurf_LineOn2S)& theLine,
d30895f5 435 const ComputationMethods::stCoeffsValue& theCoeffs,
b5ef9d91 436 const Standard_Integer theWLIndex,
437 const Standard_Integer theMinNbPoints,
438 const Standard_Integer theStartPointOnLine,
439 const Standard_Integer theEndPointOnLine,
b5ef9d91 440 const Standard_Real theTol2D,
441 const Standard_Real thePeriodOfSurf2,
442 const Standard_Boolean isTheReverse);
443
444//=======================================================================
445//function : MinMax
446//purpose : Replaces theParMIN = MIN(theParMIN, theParMAX),
9e20ed57 447// theParMAX = MAX(theParMIN, theParMAX).
b5ef9d91 448//=======================================================================
9e20ed57 449static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX)
450{
451 if(theParMIN > theParMAX)
452 {
453 const Standard_Real aux = theParMAX;
454 theParMAX = theParMIN;
455 theParMIN = aux;
456 }
457}
458
d30895f5 459//=======================================================================
460//function : ExtremaLineLine
461//purpose : Computes extrema between the given lines. Returns parameters
462// on correspond curve.
463//=======================================================================
464static inline void ExtremaLineLine(const gp_Ax1& theC1,
465 const gp_Ax1& theC2,
466 const Standard_Real theCosA,
467 const Standard_Real theSqSinA,
468 Standard_Real& thePar1,
469 Standard_Real& thePar2)
470{
471 const gp_Dir &aD1 = theC1.Direction(),
472 &aD2 = theC2.Direction();
473
474 const gp_XYZ aL1L2 = theC2.Location().XYZ() - theC1.Location().XYZ();
475 const Standard_Real aD1L = aD1.XYZ().Dot(aL1L2),
476 aD2L = aD2.XYZ().Dot(aL1L2);
477
478 thePar1 = (aD1L - theCosA * aD2L) / theSqSinA;
479 thePar2 = (theCosA * aD1L - aD2L) / theSqSinA;
480}
481
b5ef9d91 482//=======================================================================
483//function : VBoundaryPrecise
484//purpose : By default, we shall consider, that V1 and V2 will increase
485// if U1 increases. But if it is not, new V1set and/or V2set
486// must be computed as [V1current - DeltaV1] (analogically
487// for V2). This function processes this case.
488//=======================================================================
489static void VBoundaryPrecise( const math_Matrix& theMatr,
490 const Standard_Real theV1AfterDecrByDelta,
491 const Standard_Real theV2AfterDecrByDelta,
b5ef9d91 492 Standard_Real& theV1Set,
493 Standard_Real& theV2Set)
494{
495 //Now we are going to define if V1 (and V2) increases
496 //(or decreases) when U1 will increase.
497 const Standard_Integer aNbDim = 3;
498 math_Matrix aSyst(1, aNbDim, 1, aNbDim);
499
500 aSyst.SetCol(1, theMatr.Col(1));
501 aSyst.SetCol(2, theMatr.Col(2));
502 aSyst.SetCol(3, theMatr.Col(4));
503
504 const Standard_Real aDet = aSyst.Determinant();
505
506 aSyst.SetCol(1, theMatr.Col(3));
507 const Standard_Real aDet1 = aSyst.Determinant();
508
509 aSyst.SetCol(1, theMatr.Col(1));
510 aSyst.SetCol(2, theMatr.Col(3));
511
512 const Standard_Real aDet2 = aSyst.Determinant();
513
514 if(aDet*aDet1 > 0.0)
515 {
7365fad6 516 theV1Set = theV1AfterDecrByDelta;
b5ef9d91 517 }
518
519 if(aDet*aDet2 > 0.0)
520 {
7365fad6 521 theV2Set = theV2AfterDecrByDelta;
b5ef9d91 522 }
523}
524
525//=======================================================================
526//function : DeltaU1Computing
527//purpose : Computes new step for U1 parameter.
528//=======================================================================
529static inline
530 Standard_Boolean DeltaU1Computing(const math_Matrix& theSyst,
531 const math_Vector& theFree,
532 Standard_Real& theDeltaU1Found)
533{
534 Standard_Real aDet = theSyst.Determinant();
535
536 if(Abs(aDet) > aNulValue)
537 {
538 math_Matrix aSyst1(theSyst);
539 aSyst1.SetCol(2, theFree);
540
541 theDeltaU1Found = Abs(aSyst1.Determinant()/aDet);
542 return Standard_True;
543 }
544
545 return Standard_False;
546}
547
548//=======================================================================
549//function : StepComputing
550//purpose :
551//
552//Attention!!!:
553// theMatr must have 3*5-dimension strictly.
554// For system
555// {a11*V1+a12*V2+a13*dU1+a14*dU2=b1;
556// {a21*V1+a22*V2+a23*dU1+a24*dU2=b2;
557// {a31*V1+a32*V2+a33*dU1+a34*dU2=b3;
558// theMatr must be following:
559// (a11 a12 a13 a14 b1)
560// (a21 a22 a23 a24 b2)
561// (a31 a32 a33 a34 b3)
562//=======================================================================
563static Standard_Boolean StepComputing(const math_Matrix& theMatr,
564 const Standard_Real theV1Cur,
565 const Standard_Real theV2Cur,
566 const Standard_Real theDeltaV1,
567 const Standard_Real theDeltaV2,
b5ef9d91 568 Standard_Real& theDeltaU1Found/*,
569 Standard_Real& theDeltaU2Found,
570 Standard_Real& theV1Found,
571 Standard_Real& theV2Found*/)
572{
d30895f5 573#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
7365fad6 574 bool flShow = false;
575
576 if(flShow)
577 {
578 printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
579 theMatr(1,1), theMatr(1,2), theMatr(1,3), theMatr(1,4), theMatr(1,5));
580 printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
581 theMatr(2,1), theMatr(2,2), theMatr(2,3), theMatr(2,4), theMatr(2,5));
582 printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
583 theMatr(3,1), theMatr(3,2), theMatr(3,3), theMatr(3,4), theMatr(3,5));
584 }
585#endif
586
b5ef9d91 587 Standard_Boolean isSuccess = Standard_False;
588 theDeltaU1Found/* = theDeltaU2Found*/ = RealLast();
589 //theV1Found = theV1set;
590 //theV2Found = theV2Set;
591 const Standard_Integer aNbDim = 3;
592
593 math_Matrix aSyst(1, aNbDim, 1, aNbDim);
594 math_Vector aFree(1, aNbDim);
595
596 //By default, increasing V1(U1) and V2(U1) functions is
597 //considered
7365fad6 598 Standard_Real aV1Set = theV1Cur + theDeltaV1,
599 aV2Set = theV2Cur + theDeltaV2;
b5ef9d91 600
7365fad6 601 //However, what is indeed?
602 VBoundaryPrecise( theMatr, theV1Cur - theDeltaV1,
603 theV2Cur - theDeltaV2, aV1Set, aV2Set);
b5ef9d91 604
605 aSyst.SetCol(2, theMatr.Col(3));
606 aSyst.SetCol(3, theMatr.Col(4));
607
608 for(Standard_Integer i = 0; i < 2; i++)
609 {
610 if(i == 0)
611 {//V1 is known
612 aSyst.SetCol(1, theMatr.Col(2));
613 aFree.Set(1, aNbDim, theMatr.Col(5)-aV1Set*theMatr.Col(1));
614 }
615 else
616 {//i==1 => V2 is known
617 aSyst.SetCol(1, theMatr.Col(1));
618 aFree.Set(1, aNbDim, theMatr.Col(5)-aV2Set*theMatr.Col(2));
619 }
620
621 Standard_Real aNewDU = theDeltaU1Found;
622 if(DeltaU1Computing(aSyst, aFree, aNewDU))
623 {
624 isSuccess = Standard_True;
625 if(aNewDU < theDeltaU1Found)
626 {
627 theDeltaU1Found = aNewDU;
628 }
629 }
630 }
631
632 if(!isSuccess)
633 {
634 aFree = theMatr.Col(5) - aV1Set*theMatr.Col(1) - aV2Set*theMatr.Col(2);
635 math_Matrix aSyst1(1, aNbDim, 1, 2);
636 aSyst1.SetCol(1, aSyst.Col(2));
637 aSyst1.SetCol(2, aSyst.Col(3));
638
639 //Now we have overdetermined system.
640
641 const Standard_Real aDet1 = theMatr(1,3)*theMatr(2,4) - theMatr(2,3)*theMatr(1,4);
642 const Standard_Real aDet2 = theMatr(1,3)*theMatr(3,4) - theMatr(3,3)*theMatr(1,4);
643 const Standard_Real aDet3 = theMatr(2,3)*theMatr(3,4) - theMatr(3,3)*theMatr(2,4);
644 const Standard_Real anAbsD1 = Abs(aDet1);
645 const Standard_Real anAbsD2 = Abs(aDet2);
646 const Standard_Real anAbsD3 = Abs(aDet3);
647
648 if(anAbsD1 >= anAbsD2)
649 {
650 if(anAbsD1 >= anAbsD3)
651 {
652 //Det1
653 if(anAbsD1 <= aNulValue)
654 return isSuccess;
655
656 theDeltaU1Found = Abs(aFree(1)*theMatr(2,4) - aFree(2)*theMatr(1,4))/anAbsD1;
657 isSuccess = Standard_True;
658 }
659 else
660 {
661 //Det3
662 if(anAbsD3 <= aNulValue)
663 return isSuccess;
664
665 theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
666 isSuccess = Standard_True;
667 }
668 }
669 else
670 {
671 if(anAbsD2 >= anAbsD3)
672 {
673 //Det2
674 if(anAbsD2 <= aNulValue)
675 return isSuccess;
676
677 theDeltaU1Found = Abs(aFree(1)*theMatr(3,4) - aFree(3)*theMatr(1,4))/anAbsD2;
678 isSuccess = Standard_True;
679 }
680 else
681 {
682 //Det3
683 if(anAbsD3 <= aNulValue)
684 return isSuccess;
685
686 theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
687 isSuccess = Standard_True;
688 }
689 }
690 }
691
692 return isSuccess;
693}
9e20ed57 694
7fd59977 695//=======================================================================
696//function : ProcessBounds
697//purpose :
698//=======================================================================
699void ProcessBounds(const Handle(IntPatch_ALine)& alig, //-- ligne courante
700 const IntPatch_SequenceOfLine& slin,
701 const IntSurf_Quadric& Quad1,
702 const IntSurf_Quadric& Quad2,
703 Standard_Boolean& procf,
704 const gp_Pnt& ptf, //-- Debut Ligne Courante
705 const Standard_Real first, //-- Paramf
706 Standard_Boolean& procl,
707 const gp_Pnt& ptl, //-- Fin Ligne courante
708 const Standard_Real last, //-- Paraml
709 Standard_Boolean& Multpoint,
710 const Standard_Real Tol)
711{
712 Standard_Integer j,k;
713 Standard_Real U1,V1,U2,V2;
714 IntPatch_Point ptsol;
715 Standard_Real d;
716
717 if (procf && procl) {
718 j = slin.Length() + 1;
719 }
720 else {
721 j = 1;
722 }
723
724
725 //-- On prend les lignes deja enregistrees
726
727 while (j <= slin.Length()) {
728 if(slin.Value(j)->ArcType() == IntPatch_Analytic) {
729 const Handle(IntPatch_ALine)& aligold = *((Handle(IntPatch_ALine)*)&slin.Value(j));
730 k = 1;
731
732 //-- On prend les vertex des lignes deja enregistrees
733
734 while (k <= aligold->NbVertex()) {
735 ptsol = aligold->Vertex(k);
736 if (!procf) {
737 d=ptf.Distance(ptsol.Value());
738 if (d <= Tol) {
739 if (!ptsol.IsMultiple()) {
740 //-- le point ptsol (de aligold) est declare multiple sur aligold
741 Multpoint = Standard_True;
742 ptsol.SetMultiple(Standard_True);
743 aligold->Replace(k,ptsol);
744 }
745 ptsol.SetParameter(first);
746 alig->AddVertex(ptsol);
747 alig->SetFirstPoint(alig->NbVertex());
748 procf = Standard_True;
749
750 //-- On restore le point avec son parametre sur aligold
751 ptsol = aligold->Vertex(k);
752
753 }
754 }
755 if (!procl) {
756 if (ptl.Distance(ptsol.Value()) <= Tol) {
757 if (!ptsol.IsMultiple()) {
758 Multpoint = Standard_True;
759 ptsol.SetMultiple(Standard_True);
760 aligold->Replace(k,ptsol);
761 }
762 ptsol.SetParameter(last);
763 alig->AddVertex(ptsol);
764 alig->SetLastPoint(alig->NbVertex());
765 procl = Standard_True;
766
767 //-- On restore le point avec son parametre sur aligold
768 ptsol = aligold->Vertex(k);
769
770 }
771 }
772 if (procf && procl) {
773 k = aligold->NbVertex()+1;
774 }
775 else {
776 k = k+1;
777 }
778 }
779 if (procf && procl) {
780 j = slin.Length()+1;
781 }
782 else {
783 j = j+1;
784 }
785 }
786 }
787 if (!procf && !procl) {
788 Quad1.Parameters(ptf,U1,V1);
789 Quad2.Parameters(ptf,U2,V2);
790 ptsol.SetValue(ptf,Tol,Standard_False);
791 ptsol.SetParameters(U1,V1,U2,V2);
792 ptsol.SetParameter(first);
793 if (ptf.Distance(ptl) <= Tol) {
794 ptsol.SetMultiple(Standard_True); // a voir
795 Multpoint = Standard_True; // a voir de meme
796 alig->AddVertex(ptsol);
797 alig->SetFirstPoint(alig->NbVertex());
798
799 ptsol.SetParameter(last);
800 alig->AddVertex(ptsol);
801 alig->SetLastPoint(alig->NbVertex());
802 }
803 else {
804 alig->AddVertex(ptsol);
805 alig->SetFirstPoint(alig->NbVertex());
806 Quad1.Parameters(ptl,U1,V1);
807 Quad2.Parameters(ptl,U2,V2);
808 ptsol.SetValue(ptl,Tol,Standard_False);
809 ptsol.SetParameters(U1,V1,U2,V2);
810 ptsol.SetParameter(last);
811 alig->AddVertex(ptsol);
812 alig->SetLastPoint(alig->NbVertex());
813 }
814 }
815 else if (!procf) {
816 Quad1.Parameters(ptf,U1,V1);
817 Quad2.Parameters(ptf,U2,V2);
818 ptsol.SetValue(ptf,Tol,Standard_False);
819 ptsol.SetParameters(U1,V1,U2,V2);
820 ptsol.SetParameter(first);
821 alig->AddVertex(ptsol);
822 alig->SetFirstPoint(alig->NbVertex());
823 }
824 else if (!procl) {
825 Quad1.Parameters(ptl,U1,V1);
826 Quad2.Parameters(ptl,U2,V2);
827 ptsol.SetValue(ptl,Tol,Standard_False);
828 ptsol.SetParameters(U1,V1,U2,V2);
829 ptsol.SetParameter(last);
830 alig->AddVertex(ptsol);
831 alig->SetLastPoint(alig->NbVertex());
832 }
833}
d30895f5 834
7fd59977 835//=======================================================================
d30895f5 836//function : CyCyAnalyticalIntersect
7fd59977 837//purpose :
838//=======================================================================
d30895f5 839Standard_Boolean CyCyAnalyticalIntersect( const IntSurf_Quadric& Quad1,
840 const IntSurf_Quadric& Quad2,
841 const IntAna_QuadQuadGeo& theInter,
842 const Standard_Real Tol,
843 const Standard_Boolean isTheReverse,
844 Standard_Boolean& Empty,
845 Standard_Boolean& Same,
846 Standard_Boolean& Multpoint,
847 IntPatch_SequenceOfLine& slin,
848 IntPatch_SequenceOfPoint& spnt)
7fd59977 849{
850 IntPatch_Point ptsol;
d30895f5 851
7fd59977 852 Standard_Integer i;
853
854 IntSurf_TypeTrans trans1,trans2;
855 IntAna_ResultType typint;
856
857 gp_Elips elipsol;
858 gp_Lin linsol;
859
860 gp_Cylinder Cy1(Quad1.Cylinder());
861 gp_Cylinder Cy2(Quad2.Cylinder());
862
d30895f5 863 typint = theInter.TypeInter();
864 Standard_Integer NbSol = theInter.NbSolutions();
7fd59977 865 Empty = Standard_False;
866 Same = Standard_False;
867
ecc4f148 868 switch (typint)
d30895f5 869 {
ecc4f148 870 case IntAna_Empty:
7fd59977 871 {
872 Empty = Standard_True;
d30895f5 873 }
7fd59977 874 break;
875
876 case IntAna_Same:
d30895f5 877 {
7fd59977 878 Same = Standard_True;
d30895f5 879 }
7fd59977 880 break;
d30895f5 881
ecc4f148 882 case IntAna_Point:
7fd59977 883 {
d30895f5 884 gp_Pnt psol(theInter.Point(1));
7fd59977 885 ptsol.SetValue(psol,Tol,Standard_True);
d30895f5 886
887 Standard_Real U1,V1,U2,V2;
888
889 if(isTheReverse)
890 {
891 Quad2.Parameters(psol, U1, V1);
892 Quad1.Parameters(psol, U2, V2);
893 }
894 else
895 {
896 Quad1.Parameters(psol, U1, V1);
897 Quad2.Parameters(psol, U2, V2);
898 }
899
7fd59977 900 ptsol.SetParameters(U1,V1,U2,V2);
901 spnt.Append(ptsol);
902 }
903 break;
904
905 case IntAna_Line:
906 {
907 gp_Pnt ptref;
ecc4f148 908 if (NbSol == 1)
909 { // Cylinders are tangent to each other by line
d30895f5 910 linsol = theInter.Line(1);
ecc4f148 911 ptref = linsol.Location();
d30895f5 912
913 //Radius-vectors
ecc4f148 914 gp_Dir crb1(gp_Vec(ptref,Cy1.Location()));
915 gp_Dir crb2(gp_Vec(ptref,Cy2.Location()));
d30895f5 916
917 //outer normal lines
ecc4f148 918 gp_Vec norm1(Quad1.Normale(ptref));
919 gp_Vec norm2(Quad2.Normale(ptref));
920 IntSurf_Situation situcyl1;
921 IntSurf_Situation situcyl2;
922
923 if (crb1.Dot(crb2) < 0.)
924 { // centre de courbures "opposes"
d30895f5 925 //ATTENTION!!!
926 // Normal and Radius-vector of the 1st(!) cylinder
927 // is used for judging what the situation of the 2nd(!)
928 // cylinder is.
929
ecc4f148 930 if (norm1.Dot(crb1) > 0.)
931 {
932 situcyl2 = IntSurf_Inside;
933 }
934 else
935 {
936 situcyl2 = IntSurf_Outside;
937 }
938
939 if (norm2.Dot(crb2) > 0.)
940 {
941 situcyl1 = IntSurf_Inside;
942 }
943 else
944 {
945 situcyl1 = IntSurf_Outside;
946 }
947 }
948 else
d30895f5 949 {
ecc4f148 950 if (Cy1.Radius() < Cy2.Radius())
951 {
952 if (norm1.Dot(crb1) > 0.)
953 {
954 situcyl2 = IntSurf_Inside;
955 }
956 else
957 {
958 situcyl2 = IntSurf_Outside;
959 }
960
961 if (norm2.Dot(crb2) > 0.)
962 {
963 situcyl1 = IntSurf_Outside;
964 }
965 else
966 {
967 situcyl1 = IntSurf_Inside;
968 }
969 }
970 else
971 {
972 if (norm1.Dot(crb1) > 0.)
973 {
974 situcyl2 = IntSurf_Outside;
975 }
976 else
977 {
978 situcyl2 = IntSurf_Inside;
979 }
980
981 if (norm2.Dot(crb2) > 0.)
982 {
983 situcyl1 = IntSurf_Inside;
984 }
985 else
986 {
987 situcyl1 = IntSurf_Outside;
988 }
989 }
990 }
991
992 Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_True, situcyl1, situcyl2);
993 slin.Append(glig);
7fd59977 994 }
ecc4f148 995 else
996 {
997 for (i=1; i <= NbSol; i++)
998 {
d30895f5 999 linsol = theInter.Line(i);
ecc4f148 1000 ptref = linsol.Location();
1001 gp_Vec lsd = linsol.Direction();
d30895f5 1002
1003 //Theoretically, qwe = +/- 1.0.
ecc4f148 1004 Standard_Real qwe = lsd.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
1005 if (qwe >0.00000001)
1006 {
1007 trans1 = IntSurf_Out;
1008 trans2 = IntSurf_In;
1009 }
1010 else if (qwe <-0.00000001)
1011 {
1012 trans1 = IntSurf_In;
1013 trans2 = IntSurf_Out;
1014 }
1015 else
1016 {
1017 trans1=trans2=IntSurf_Undecided;
1018 }
1019
1020 Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
1021 slin.Append(glig);
1022 }
7fd59977 1023 }
1024 }
1025 break;
1026
1027 case IntAna_Ellipse:
1028 {
1029 gp_Vec Tgt;
1030 gp_Pnt ptref;
fa0291ff 1031 IntPatch_Point pmult1, pmult2;
d30895f5 1032
1033 elipsol = theInter.Ellipse(1);
1034
fa0291ff 1035 gp_Pnt pttang1(ElCLib::Value(0.5*M_PI, elipsol));
1036 gp_Pnt pttang2(ElCLib::Value(1.5*M_PI, elipsol));
d30895f5 1037
7fd59977 1038 Multpoint = Standard_True;
1039 pmult1.SetValue(pttang1,Tol,Standard_True);
1040 pmult2.SetValue(pttang2,Tol,Standard_True);
1041 pmult1.SetMultiple(Standard_True);
1042 pmult2.SetMultiple(Standard_True);
d30895f5 1043
7fd59977 1044 Standard_Real oU1,oV1,oU2,oV2;
d30895f5 1045
1046 if(isTheReverse)
1047 {
1048 Quad2.Parameters(pttang1,oU1,oV1);
1049 Quad1.Parameters(pttang1,oU2,oV2);
1050 }
1051 else
1052 {
1053 Quad1.Parameters(pttang1,oU1,oV1);
1054 Quad2.Parameters(pttang1,oU2,oV2);
1055 }
1056
7fd59977 1057 pmult1.SetParameters(oU1,oV1,oU2,oV2);
1058
d30895f5 1059 if(isTheReverse)
1060 {
1061 Quad2.Parameters(pttang2,oU1,oV1);
1062 Quad1.Parameters(pttang2,oU2,oV2);
1063 }
1064 else
1065 {
1066 Quad1.Parameters(pttang2,oU1,oV1);
1067 Quad2.Parameters(pttang2,oU2,oV2);
1068 }
1069
7fd59977 1070 pmult2.SetParameters(oU1,oV1,oU2,oV2);
7fd59977 1071
d30895f5 1072 // on traite la premiere ellipse
1073
7fd59977 1074 //-- Calcul de la Transition de la ligne
1075 ElCLib::D1(0.,elipsol,ptref,Tgt);
d30895f5 1076
1077 //Theoretically, qwe = +/- |Tgt|.
7fd59977 1078 Standard_Real qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
ecc4f148 1079 if (qwe>0.00000001)
1080 {
1081 trans1 = IntSurf_Out;
1082 trans2 = IntSurf_In;
7fd59977 1083 }
ecc4f148 1084 else if (qwe<-0.00000001)
1085 {
1086 trans1 = IntSurf_In;
1087 trans2 = IntSurf_Out;
7fd59977 1088 }
ecc4f148 1089 else
d30895f5 1090 {
ecc4f148 1091 trans1=trans2=IntSurf_Undecided;
7fd59977 1092 }
ecc4f148 1093
7fd59977 1094 //-- Transition calculee au point 0 -> Trans2 , Trans1
1095 //-- car ici, on devarit calculer en PI
1096 Handle(IntPatch_GLine) glig = new IntPatch_GLine(elipsol,Standard_False,trans2,trans1);
fa0291ff 1097 //
1098 {
ecc4f148 1099 Standard_Real aU1, aV1, aU2, aV2;
1100 IntPatch_Point aIP;
1101 gp_Pnt aP (ElCLib::Value(0., elipsol));
1102 //
1103 aIP.SetValue(aP,Tol,Standard_False);
1104 aIP.SetMultiple(Standard_False);
1105 //
d30895f5 1106
1107 if(isTheReverse)
1108 {
1109 Quad2.Parameters(aP, aU1, aV1);
1110 Quad1.Parameters(aP, aU2, aV2);
1111 }
1112 else
1113 {
1114 Quad1.Parameters(aP, aU1, aV1);
1115 Quad2.Parameters(aP, aU2, aV2);
1116 }
1117
ecc4f148 1118 aIP.SetParameters(aU1, aV1, aU2, aV2);
1119 //
1120 aIP.SetParameter(0.);
1121 glig->AddVertex(aIP);
1122 glig->SetFirstPoint(1);
1123 //
1124 aIP.SetParameter(2.*M_PI);
1125 glig->AddVertex(aIP);
1126 glig->SetLastPoint(2);
fa0291ff 1127 }
1128 //
1129 pmult1.SetParameter(0.5*M_PI);
7fd59977 1130 glig->AddVertex(pmult1);
fa0291ff 1131 //
c6541a0c 1132 pmult2.SetParameter(1.5*M_PI);
7fd59977 1133 glig->AddVertex(pmult2);
fa0291ff 1134
1135 //
7fd59977 1136 slin.Append(glig);
1137
1138 //-- Transitions calculee au point 0 OK
fa0291ff 1139 //
7fd59977 1140 // on traite la deuxieme ellipse
d30895f5 1141 elipsol = theInter.Ellipse(2);
7fd59977 1142
1143 Standard_Real param1 = ElCLib::Parameter(elipsol,pttang1);
1144 Standard_Real param2 = ElCLib::Parameter(elipsol,pttang2);
ecc4f148 1145 Standard_Real parampourtransition = 0.0;
1146 if (param1 < param2)
1147 {
1148 pmult1.SetParameter(0.5*M_PI);
1149 pmult2.SetParameter(1.5*M_PI);
1150 parampourtransition = M_PI;
7fd59977 1151 }
1152 else {
ecc4f148 1153 pmult1.SetParameter(1.5*M_PI);
1154 pmult2.SetParameter(0.5*M_PI);
1155 parampourtransition = 0.0;
7fd59977 1156 }
1157
1158 //-- Calcul des transitions de ligne pour la premiere ligne
1159 ElCLib::D1(parampourtransition,elipsol,ptref,Tgt);
d30895f5 1160
1161 //Theoretically, qwe = +/- |Tgt|.
7fd59977 1162 qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
ecc4f148 1163 if(qwe> 0.00000001)
1164 {
1165 trans1 = IntSurf_Out;
1166 trans2 = IntSurf_In;
7fd59977 1167 }
ecc4f148 1168 else if(qwe< -0.00000001)
1169 {
1170 trans1 = IntSurf_In;
1171 trans2 = IntSurf_Out;
7fd59977 1172 }
ecc4f148 1173 else
1174 {
1175 trans1=trans2=IntSurf_Undecided;
7fd59977 1176 }
ecc4f148 1177
7fd59977 1178 //-- La transition a ete calculee sur un point de cette ligne
1179 glig = new IntPatch_GLine(elipsol,Standard_False,trans1,trans2);
fa0291ff 1180 //
1181 {
ecc4f148 1182 Standard_Real aU1, aV1, aU2, aV2;
1183 IntPatch_Point aIP;
1184 gp_Pnt aP (ElCLib::Value(0., elipsol));
1185 //
1186 aIP.SetValue(aP,Tol,Standard_False);
1187 aIP.SetMultiple(Standard_False);
1188 //
7fd59977 1189
d30895f5 1190 if(isTheReverse)
ecc4f148 1191 {
d30895f5 1192 Quad2.Parameters(aP, aU1, aV1);
1193 Quad1.Parameters(aP, aU2, aV2);
ecc4f148 1194 }
d30895f5 1195 else
ecc4f148 1196 {
d30895f5 1197 Quad1.Parameters(aP, aU1, aV1);
1198 Quad2.Parameters(aP, aU2, aV2);
ecc4f148 1199 }
ecc4f148 1200
d30895f5 1201 aIP.SetParameters(aU1, aV1, aU2, aV2);
1202 //
1203 aIP.SetParameter(0.);
1204 glig->AddVertex(aIP);
1205 glig->SetFirstPoint(1);
1206 //
1207 aIP.SetParameter(2.*M_PI);
1208 glig->AddVertex(aIP);
1209 glig->SetLastPoint(2);
1210 }
1211 //
1212 glig->AddVertex(pmult1);
1213 glig->AddVertex(pmult2);
1214 //
1215 slin.Append(glig);
ecc4f148 1216 }
ecc4f148 1217 break;
ecc4f148 1218
d30895f5 1219 case IntAna_Parabola:
1220 case IntAna_Hyperbola:
1221 Standard_Failure::Raise("IntCyCy(): Wrong intersection type!");
ecc4f148 1222
d30895f5 1223 case IntAna_Circle:
1224 // Circle is useful when we will work with trimmed surfaces
1225 // (two cylinders can be tangent by their basises, e.g. circle)
1226 case IntAna_NoGeometricSolution:
1227 default:
1228 return Standard_False;
1229 }
ecc4f148 1230
d30895f5 1231 return Standard_True;
1232}
ecc4f148 1233
d30895f5 1234//=======================================================================
1235//function : ShortCosForm
1236//purpose : Represents theCosFactor*cosA+theSinFactor*sinA as
1237// theCoeff*cos(A-theAngle) if it is possibly (all angles are
1238// in radians).
1239//=======================================================================
1240static void ShortCosForm( const Standard_Real theCosFactor,
1241 const Standard_Real theSinFactor,
1242 Standard_Real& theCoeff,
1243 Standard_Real& theAngle)
1244{
1245 theCoeff = sqrt(theCosFactor*theCosFactor+theSinFactor*theSinFactor);
1246 theAngle = 0.0;
1247 if(IsEqual(theCoeff, 0.0))
1248 {
1249 theAngle = 0.0;
1250 return;
1251 }
ecc4f148 1252
d30895f5 1253 theAngle = acos(Abs(theCosFactor/theCoeff));
ecc4f148 1254
d30895f5 1255 if(theSinFactor > 0.0)
1256 {
1257 if(IsEqual(theCosFactor, 0.0))
1258 {
1259 theAngle = M_PI/2.0;
1260 }
1261 else if(theCosFactor < 0.0)
1262 {
1263 theAngle = M_PI-theAngle;
1264 }
1265 }
1266 else if(IsEqual(theSinFactor, 0.0))
1267 {
1268 if(theCosFactor < 0.0)
1269 {
1270 theAngle = M_PI;
1271 }
1272 }
1273 if(theSinFactor < 0.0)
1274 {
1275 if(theCosFactor > 0.0)
1276 {
1277 theAngle = 2.0*M_PI-theAngle;
1278 }
1279 else if(IsEqual(theCosFactor, 0.0))
1280 {
1281 theAngle = 3.0*M_PI/2.0;
1282 }
1283 else if(theCosFactor < 0.0)
1284 {
1285 theAngle = M_PI+theAngle;
1286 }
1287 }
ecc4f148 1288}
1289
b5ef9d91 1290//=======================================================================
1291//function : CylCylMonotonicity
1292//purpose : Determines, if U2(U1) function is increasing.
1293//=======================================================================
d30895f5 1294Standard_Boolean ComputationMethods::CylCylMonotonicity(const Standard_Real theU1par,
1295 const Standard_Integer theWLIndex,
1296 const stCoeffsValue& theCoeffs,
1297 const Standard_Real thePeriod,
1298 Standard_Boolean& theIsIncreasing)
b5ef9d91 1299{
1300 // U2(U1) = FI2 + (+/-)acos(B*cos(U1 - FI1) + C);
1301 //Accordingly,
1302 //Func. U2(X1) = FI2 + X1;
1303 //Func. X1(X2) = anArccosFactor*X2;
1304 //Func. X2(X3) = acos(X3);
1305 //Func. X3(X4) = B*X4 + C;
1306 //Func. X4(U1) = cos(U1 - FI1).
1307 //
1308 //Consequently,
1309 //U2(X1) is always increasing.
1310 //X1(X2) is increasing if anArccosFactor > 0.0 and
1311 //is decreasing otherwise.
1312 //X2(X3) is always decreasing.
1313 //Therefore, U2(X3) is decreasing if anArccosFactor > 0.0 and
1314 //is increasing otherwise.
1315 //X3(X4) is increasing if B > 0 and is decreasing otherwise.
1316 //X4(U1) is increasing if U1 - FI1 in [PI, 2*PI) and
1317 //is decreasing U1 - FI1 in [0, PI) or if (U1 - FI1 == 2PI).
1318 //After that, we can predict behaviour of U2(U1) function:
1319 //if it is increasing or decreasing.
1320
1321 //For "+/-" sign. If isPlus == TRUE, "+" is chosen, otherwise, we choose "-".
1322 Standard_Boolean isPlus = Standard_False;
1323
1324 switch(theWLIndex)
1325 {
1326 case 0:
1327 isPlus = Standard_True;
1328 break;
1329 case 1:
1330 isPlus = Standard_False;
1331 break;
1332 default:
1333 //Standard_Failure::Raise("Error. Range Error!!!!");
1334 return Standard_False;
1335 }
1336
1337 Standard_Real aU1Temp = theU1par - theCoeffs.mFI1;
1338 InscribePoint(0, thePeriod, aU1Temp, 0.0, thePeriod, Standard_False);
1339
1340 theIsIncreasing = Standard_True;
1341
1342 if(((M_PI - aU1Temp) < RealSmall()) && (aU1Temp < thePeriod))
1343 {
1344 theIsIncreasing = Standard_False;
1345 }
1346
1347 if(theCoeffs.mB < 0.0)
1348 {
1349 theIsIncreasing = !theIsIncreasing;
1350 }
1351
1352 if(!isPlus)
1353 {
1354 theIsIncreasing = !theIsIncreasing;
1355 }
1356
1357 return Standard_True;
1358}
1359
7365fad6 1360//=======================================================================
1361//function : CylCylComputeParameters
1362//purpose : Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0,
1363// esimates the tolerance of U2-computing (estimation result is
1364// assigned to *theDelta value).
1365//=======================================================================
d30895f5 1366Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
b5ef9d91 1367 const Standard_Integer theWLIndex,
1368 const stCoeffsValue& theCoeffs,
7365fad6 1369 Standard_Real& theU2,
d30895f5 1370 Standard_Real* const theDelta)
b5ef9d91 1371{
7365fad6 1372 //This formula is got from some experience and can be changed.
1373 const Standard_Real aTol0 = Min(10.0*Epsilon(1.0)*theCoeffs.mB, aNulValue);
1374 const Standard_Real aTol = 1.0 - aTol0;
b5ef9d91 1375
7365fad6 1376 if(theWLIndex < 0 || theWLIndex > 1)
b5ef9d91 1377 return Standard_False;
b5ef9d91 1378
7365fad6 1379 const Standard_Real aSign = theWLIndex ? -1.0 : 1.0;
1380
1381 Standard_Real anArg = cos(theU1par - theCoeffs.mFI1);
1382 anArg = theCoeffs.mB*anArg + theCoeffs.mC;
1383
1384 if(anArg > aTol)
1385 {
1386 if(theDelta)
1387 *theDelta = 0.0;
b5ef9d91 1388
b5ef9d91 1389 anArg = 1.0;
7365fad6 1390 }
1391 else if(anArg < -aTol)
1392 {
1393 if(theDelta)
1394 *theDelta = 0.0;
b5ef9d91 1395
b5ef9d91 1396 anArg = -1.0;
7365fad6 1397 }
1398 else if(theDelta)
1399 {
1400 //There is a case, when
1401 // const double aPar = 0.99999999999721167;
1402 // const double aFI2 = 2.3593296083566181e-006;
1403
1404 //Then
1405 // aPar - cos(aFI2) == -5.10703e-015 ==> cos(aFI2) == aPar.
1406 //Theoreticaly, in this case
1407 // aFI2 == +/- acos(aPar).
1408 //However,
1409 // acos(aPar) - aFI2 == 2.16362e-009.
1410 //Error is quite big.
1411
1412 //This error should be estimated. Let use following way, which is based
1413 //on expanding to Taylor series.
1414
1415 // acos(p)-acos(p+x) = x/sqrt(1-p*p).
1416
1417 //If p == (1-d) (when p > 0) or p == (-1+d) (when p < 0) then
1418 // acos(p)-acos(p+x) = x/sqrt(d*(2-d)).
1419
1420 //Here always aTol0 <= d <= 1. Max(x) is considered (!) to be equal to aTol0.
1421 //In this case
1422 // 8*aTol0 <= acos(p)-acos(p+x) <= sqrt(2/(2-aTol0)-1),
1423 // because 0 < aTol0 < 1.
1424 //E.g. when aTol0 = 1.0e-11,
1425 // 8.0e-11 <= acos(p)-acos(p+x) < 2.24e-6.
1426
1427 const Standard_Real aDelta = Min(1.0-anArg, 1.0+anArg);
1428 Standard_DivideByZero_Raise_if((aDelta*aDelta < RealSmall()) || (aDelta >= 2.0),
1429 "IntPatch_ImpImpIntersection_4.gxx, CylCylComputeParameters()");
1430 *theDelta = aTol0/sqrt(aDelta*(2.0-aDelta));
1431 }
b5ef9d91 1432
7365fad6 1433 theU2 = acos(anArg);
1434 theU2 = theCoeffs.mFI2 + aSign*theU2;
b5ef9d91 1435
1436 return Standard_True;
1437}
1438
d30895f5 1439Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1,
b5ef9d91 1440 const Standard_Real theU2,
1441 const stCoeffsValue& theCoeffs,
1442 Standard_Real& theV1,
1443 Standard_Real& theV2)
1444{
1445 theV1 = theCoeffs.mK21 * sin(theU2) +
1446 theCoeffs.mK11 * sin(theU1) +
1447 theCoeffs.mL21 * cos(theU2) +
1448 theCoeffs.mL11 * cos(theU1) + theCoeffs.mM1;
1449
1450 theV2 = theCoeffs.mK22 * sin(theU2) +
1451 theCoeffs.mK12 * sin(theU1) +
1452 theCoeffs.mL22 * cos(theU2) +
1453 theCoeffs.mL12 * cos(theU1) + theCoeffs.mM2;
1454
1455 return Standard_True;
1456}
1457
1458
d30895f5 1459Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
b5ef9d91 1460 const Standard_Integer theWLIndex,
1461 const stCoeffsValue& theCoeffs,
1462 Standard_Real& theU2,
1463 Standard_Real& theV1,
1464 Standard_Real& theV2)
1465{
1466 if(!CylCylComputeParameters(theU1par, theWLIndex, theCoeffs, theU2))
1467 return Standard_False;
1468
1469 if(!CylCylComputeParameters(theU1par, theU2, theCoeffs, theV1, theV2))
1470 return Standard_False;
1471
1472 return Standard_True;
1473}
1474
ecc4f148 1475//=======================================================================
1476//function : SearchOnVBounds
1477//purpose :
1478//=======================================================================
d30895f5 1479Standard_Boolean WorkWithBoundaries::
1480 SearchOnVBounds(const SearchBoundType theSBType,
ecc4f148 1481 const Standard_Real theVzad,
b5ef9d91 1482 const Standard_Real theVInit,
ecc4f148 1483 const Standard_Real theInitU2,
1484 const Standard_Real theInitMainVar,
d30895f5 1485 Standard_Real& theMainVariableValue) const
ecc4f148 1486{
b5ef9d91 1487 const Standard_Integer aNbDim = 3;
ecc4f148 1488 const Standard_Real aMaxError = 4.0*M_PI; // two periods
1489
b5ef9d91 1490 theMainVariableValue = theInitMainVar;
1491 const Standard_Real aTol2 = 1.0e-18;
1492 Standard_Real aMainVarPrev = theInitMainVar, aU2Prev = theInitU2, anOtherVar = theVInit;
ecc4f148 1493
b5ef9d91 1494 //Structure of aMatr:
1495 // C_{1}*U_{1} & C_{2}*U_{2} & C_{3}*V_{*},
1496 //where C_{1}, C_{2} and C_{3} are math_Vector.
1497 math_Matrix aMatr(1, aNbDim, 1, aNbDim);
1498
ecc4f148 1499 Standard_Real anError = RealLast();
b5ef9d91 1500 Standard_Real anErrorPrev = anError;
ecc4f148 1501 Standard_Integer aNbIter = 0;
1502 do
1503 {
1504 if(++aNbIter > 1000)
1505 return Standard_False;
1506
b5ef9d91 1507 const Standard_Real aSinU1 = sin(aMainVarPrev),
1508 aCosU1 = cos(aMainVarPrev),
1509 aSinU2 = sin(aU2Prev),
1510 aCosU2 = cos(aU2Prev);
ecc4f148 1511
d30895f5 1512 math_Vector aVecFreeMem = (myCoeffs.mVecA2 * aU2Prev +
1513 myCoeffs.mVecB2) * aSinU2 -
1514 (myCoeffs.mVecB2 * aU2Prev -
1515 myCoeffs.mVecA2) * aCosU2 +
1516 (myCoeffs.mVecA1 * aMainVarPrev +
1517 myCoeffs.mVecB1) * aSinU1 -
1518 (myCoeffs.mVecB1 * aMainVarPrev -
1519 myCoeffs.mVecA1) * aCosU1 +
1520 myCoeffs.mVecD;
b5ef9d91 1521
1522 math_Vector aMSum(1, 3);
ecc4f148 1523
1524 switch(theSBType)
1525 {
1526 case SearchV1:
d30895f5 1527 aMatr.SetCol(3, myCoeffs.mVecC2);
1528 aMSum = myCoeffs.mVecC1 * theVzad;
b5ef9d91 1529 aVecFreeMem -= aMSum;
d30895f5 1530 aMSum += myCoeffs.mVecC2*anOtherVar;
ecc4f148 1531 break;
1532
1533 case SearchV2:
d30895f5 1534 aMatr.SetCol(3, myCoeffs.mVecC1);
1535 aMSum = myCoeffs.mVecC2 * theVzad;
b5ef9d91 1536 aVecFreeMem -= aMSum;
d30895f5 1537 aMSum += myCoeffs.mVecC1*anOtherVar;
ecc4f148 1538 break;
1539
1540 default:
1541 return Standard_False;
1542 }
1543
d30895f5 1544 aMatr.SetCol(1, myCoeffs.mVecA1 * aSinU1 - myCoeffs.mVecB1 * aCosU1);
1545 aMatr.SetCol(2, myCoeffs.mVecA2 * aSinU2 - myCoeffs.mVecB2 * aCosU2);
b5ef9d91 1546
1547 Standard_Real aDetMainSyst = aMatr.Determinant();
ecc4f148 1548
b5ef9d91 1549 if(Abs(aDetMainSyst) < aNulValue)
ecc4f148 1550 {
1551 return Standard_False;
1552 }
1553
b5ef9d91 1554 math_Matrix aM1(aMatr), aM2(aMatr), aM3(aMatr);
1555 aM1.SetCol(1, aVecFreeMem);
1556 aM2.SetCol(2, aVecFreeMem);
1557 aM3.SetCol(3, aVecFreeMem);
ecc4f148 1558
b5ef9d91 1559 const Standard_Real aDetMainVar = aM1.Determinant();
1560 const Standard_Real aDetVar1 = aM2.Determinant();
1561 const Standard_Real aDetVar2 = aM3.Determinant();
ecc4f148 1562
1563 Standard_Real aDelta = aDetMainVar/aDetMainSyst-aMainVarPrev;
1564
1565 if(Abs(aDelta) > aMaxError)
1566 return Standard_False;
1567
1568 anError = aDelta*aDelta;
1569 aMainVarPrev += aDelta;
1570
1571 ///
1572 aDelta = aDetVar1/aDetMainSyst-aU2Prev;
1573
1574 if(Abs(aDelta) > aMaxError)
1575 return Standard_False;
1576
1577 anError += aDelta*aDelta;
1578 aU2Prev += aDelta;
1579
1580 ///
1581 aDelta = aDetVar2/aDetMainSyst-anOtherVar;
1582 anError += aDelta*aDelta;
1583 anOtherVar += aDelta;
b5ef9d91 1584
1585 if(anError > anErrorPrev)
1586 {//Method diverges. Keep the best result
51740958 1587 const Standard_Real aSinU1Last = sin(aMainVarPrev),
1588 aCosU1Last = cos(aMainVarPrev),
1589 aSinU2Last = sin(aU2Prev),
1590 aCosU2Last = cos(aU2Prev);
d30895f5 1591 aMSum -= (myCoeffs.mVecA1*aCosU1Last +
1592 myCoeffs.mVecB1*aSinU1Last +
1593 myCoeffs.mVecA2*aCosU2Last +
1594 myCoeffs.mVecB2*aSinU2Last +
1595 myCoeffs.mVecD);
b5ef9d91 1596 const Standard_Real aSQNorm = aMSum.Norm2();
1597 return (aSQNorm < aTol2);
1598 }
1599 else
1600 {
1601 theMainVariableValue = aMainVarPrev;
1602 }
1603
1604 anErrorPrev = anError;
ecc4f148 1605 }
1606 while(anError > aTol2);
1607
1608 theMainVariableValue = aMainVarPrev;
1609
1610 return Standard_True;
1611}
1612
1613//=======================================================================
1614//function : InscribePoint
e002f1ce 1615//purpose : If theFlForce==TRUE theUGiven will be changed forcefully
1616// even if theUGiven is already inscibed in the boundary
1617// (if it is possible; i.e. if new theUGiven is inscribed
1618// in the boundary, too).
ecc4f148 1619//=======================================================================
b5ef9d91 1620Standard_Boolean InscribePoint( const Standard_Real theUfTarget,
1621 const Standard_Real theUlTarget,
1622 Standard_Real& theUGiven,
1623 const Standard_Real theTol2D,
1624 const Standard_Real thePeriod,
1625 const Standard_Boolean theFlForce)
ecc4f148 1626{
7365fad6 1627 if(Precision::IsInfinite(theUGiven))
1628 {
1629 return Standard_False;
1630 }
1631
02effd35 1632 if((theUfTarget - theUGiven <= theTol2D) &&
1633 (theUGiven - theUlTarget <= theTol2D))
ecc4f148 1634 {//it has already inscribed
1635
1636 /*
1637 Utf U Utl
1638 + * +
1639 */
1640
02effd35 1641 if(theFlForce)
1642 {
1643 Standard_Real anUtemp = theUGiven + thePeriod;
1644 if((theUfTarget - anUtemp <= theTol2D) &&
1645 (anUtemp - theUlTarget <= theTol2D))
1646 {
1647 theUGiven = anUtemp;
1648 return Standard_True;
1649 }
1650
1651 anUtemp = theUGiven - thePeriod;
1652 if((theUfTarget - anUtemp <= theTol2D) &&
1653 (anUtemp - theUlTarget <= theTol2D))
1654 {
1655 theUGiven = anUtemp;
1656 }
1657 }
1658
ecc4f148 1659 return Standard_True;
1660 }
1661
d30895f5 1662 const Standard_Real aUf = theUfTarget - theTol2D;
1663 const Standard_Real aUl = aUf + thePeriod;
ecc4f148 1664
d30895f5 1665 theUGiven = ElCLib::InPeriod(theUGiven, aUf, aUl);
ecc4f148 1666
7365fad6 1667 return ((theUfTarget - theUGiven <= theTol2D) &&
1668 (theUGiven - theUlTarget <= theTol2D));
ecc4f148 1669}
1670
1671//=======================================================================
1672//function : InscribeInterval
1673//purpose : Adjusts theUfGiven and after that fits theUlGiven to result
1674//=======================================================================
1675static Standard_Boolean InscribeInterval(const Standard_Real theUfTarget,
1676 const Standard_Real theUlTarget,
1677 Standard_Real& theUfGiven,
1678 Standard_Real& theUlGiven,
1679 const Standard_Real theTol2D,
1680 const Standard_Real thePeriod)
1681{
1682 Standard_Real anUpar = theUfGiven;
1683 const Standard_Real aDelta = theUlGiven - theUfGiven;
02effd35 1684 if(InscribePoint(theUfTarget, theUlTarget, anUpar,
1685 theTol2D, thePeriod, Standard_False))
ecc4f148 1686 {
1687 theUfGiven = anUpar;
1688 theUlGiven = theUfGiven + aDelta;
1689 }
1690 else
1691 {
1692 anUpar = theUlGiven;
02effd35 1693 if(InscribePoint(theUfTarget, theUlTarget, anUpar,
1694 theTol2D, thePeriod, Standard_False))
ecc4f148 1695 {
1696 theUlGiven = anUpar;
1697 theUfGiven = theUlGiven - aDelta;
1698 }
1699 else
7fd59977 1700 {
1701 return Standard_False;
1702 }
1703 }
ecc4f148 1704
1705 return Standard_True;
1706}
1707
1708//=======================================================================
7365fad6 1709//function : ExcludeNearElements
1710//purpose : Checks if theArr contains two almost equal elements.
1711// If it is true then one of equal elements will be excluded
1712// (made infinite).
1713// Returns TRUE if at least one element of theArr has been changed.
1714//ATTENTION!!!
1715// 1. Every not infinite element of theArr is considered to be
1716// in [0, T] interval (where T is the period);
1717// 2. theArr must be sorted in ascending order.
ecc4f148 1718//=======================================================================
7365fad6 1719static Standard_Boolean ExcludeNearElements(Standard_Real theArr[],
1720 const Standard_Integer theNOfMembers,
79997052 1721 const Standard_Real theUSurf1f,
1722 const Standard_Real theUSurf1l,
7365fad6 1723 const Standard_Real theTol)
ecc4f148 1724{
7365fad6 1725 Standard_Boolean aRetVal = Standard_False;
1726 for(Standard_Integer i = 1; i < theNOfMembers; i++)
ecc4f148 1727 {
7365fad6 1728 Standard_Real &anA = theArr[i],
1729 &anB = theArr[i-1];
ecc4f148 1730
7365fad6 1731 //Here, anA >= anB
ecc4f148 1732
7365fad6 1733 if(Precision::IsInfinite(anA))
1734 break;
ecc4f148 1735
7365fad6 1736 if((anA-anB) < theTol)
ecc4f148 1737 {
79997052 1738 if((anB != 0.0) && (anB != theUSurf1f) && (anB != theUSurf1l))
d30895f5 1739 anA = (anA + anB)/2.0;
79997052 1740 else
1741 anA = anB;
ecc4f148 1742
7365fad6 1743 //Make this element infinite an forget it
1744 //(we will not use it in next iterations).
1745 anB = Precision::Infinite();
1746 aRetVal = Standard_True;
ecc4f148 1747 }
1748 }
ecc4f148 1749
7365fad6 1750 return aRetVal;
1751}
ecc4f148 1752
1753//=======================================================================
1754//function : AddPointIntoWL
1755//purpose : Surf1 is a surface, whose U-par is variable.
1756//=======================================================================
1757static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
1758 const IntSurf_Quadric& theQuad2,
d30895f5 1759 const ComputationMethods::stCoeffsValue& theCoeffs,
ecc4f148 1760 const Standard_Boolean isTheReverse,
b5ef9d91 1761 const Standard_Boolean isThePrecise,
ecc4f148 1762 const gp_Pnt2d& thePntOnSurf1,
1763 const gp_Pnt2d& thePntOnSurf2,
1764 const Standard_Real theUfSurf1,
1765 const Standard_Real theUlSurf1,
b5ef9d91 1766 const Standard_Real theUfSurf2,
1767 const Standard_Real theUlSurf2,
1768 const Standard_Real theVfSurf1,
1769 const Standard_Real theVlSurf1,
1770 const Standard_Real theVfSurf2,
1771 const Standard_Real theVlSurf2,
ecc4f148 1772 const Standard_Real thePeriodOfSurf1,
1773 const Handle(IntSurf_LineOn2S)& theLine,
b5ef9d91 1774 const Standard_Integer theWLIndex,
e8feb725 1775 const Standard_Real theTol3D,
02effd35 1776 const Standard_Real theTol2D,
b5ef9d91 1777 const Standard_Boolean theFlForce,
1778 const Standard_Boolean theFlBefore = Standard_False)
ecc4f148 1779{
1780 gp_Pnt aPt1(theQuad1.Value(thePntOnSurf1.X(), thePntOnSurf1.Y())),
1781 aPt2(theQuad2.Value(thePntOnSurf2.X(), thePntOnSurf2.Y()));
1782
b5ef9d91 1783 Standard_Real aU1par = thePntOnSurf1.X();
1784 if(!InscribePoint(theUfSurf1, theUlSurf1, aU1par, theTol2D,
02effd35 1785 thePeriodOfSurf1, theFlForce))
ecc4f148 1786 return Standard_False;
1787
b5ef9d91 1788 Standard_Real aU2par = thePntOnSurf2.X();
1789 if(!InscribePoint(theUfSurf2, theUlSurf2, aU2par, theTol2D,
1790 thePeriodOfSurf1, Standard_False))
1791 return Standard_False;
1792
1793 Standard_Real aV1par = thePntOnSurf1.Y();
1794 if((aV1par - theVlSurf1 > theTol2D) || (theVfSurf1 - aV1par > theTol2D))
1795 return Standard_False;
1796
1797 Standard_Real aV2par = thePntOnSurf2.Y();
1798 if((aV2par - theVlSurf2 > theTol2D) || (theVfSurf2 - aV2par > theTol2D))
1799 return Standard_False;
1800
ecc4f148 1801 IntSurf_PntOn2S aPnt;
1802
1803 if(isTheReverse)
1804 {
1805 aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
b5ef9d91 1806 aU2par, aV2par,
1807 aU1par, aV1par);
ecc4f148 1808 }
1809 else
1810 {
1811 aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
b5ef9d91 1812 aU1par, aV1par,
1813 aU2par, aV2par);
ecc4f148 1814 }
1815
b5ef9d91 1816 Standard_Integer aNbPnts = theLine->NbPoints();
e8feb725 1817 if(aNbPnts > 0)
1818 {
1819 Standard_Real aUl = 0.0, aVl = 0.0;
1820 const IntSurf_PntOn2S aPlast = theLine->Value(aNbPnts);
1821 if(isTheReverse)
1822 aPlast.ParametersOnS2(aUl, aVl);
1823 else
1824 aPlast.ParametersOnS1(aUl, aVl);
1825
b5ef9d91 1826 if(!theFlBefore && (aU1par <= aUl))
d30895f5 1827 {
1828 //Parameter value must be increased if theFlBefore == FALSE.
1829
1830 aU1par += thePeriodOfSurf1;
1831
1832 //The condition is as same as in
1833 //InscribePoint(...) function
1834 if((theUfSurf1 - aU1par > theTol2D) ||
1835 (aU1par - theUlSurf1 > theTol2D))
1836 {
1837 //New aU1par is out of target interval.
1838 //Go back to old value.
1839
1840 return Standard_False;
1841 }
e8feb725 1842 }
1843
1844 //theTol2D is minimal step along parameter changed.
1845 //Therefore, if we apply this minimal step two
1846 //neighbour points will be always "same". Consequently,
1847 //we should reduce tolerance for IsSame checking.
1848 const Standard_Real aDTol = 1.0-Epsilon(1.0);
1849 if(aPnt.IsSame(aPlast, theTol3D*aDTol, theTol2D*aDTol))
1850 {
1851 theLine->RemovePoint(aNbPnts);
1852 }
1853 }
1854
ecc4f148 1855 theLine->Add(aPnt);
b5ef9d91 1856
1857 if(!isThePrecise)
1858 return Standard_True;
1859
1860 //Try to precise existing WLine
1861 aNbPnts = theLine->NbPoints();
1862 if(aNbPnts >= 3)
1863 {
1864 Standard_Real aU1 = 0.0, aU2 = 0.0, aU3 = 0.0, aV = 0.0;
1865 if(isTheReverse)
1866 {
1867 theLine->Value(aNbPnts).ParametersOnS2(aU3, aV);
1868 theLine->Value(aNbPnts-1).ParametersOnS2(aU2, aV);
1869 theLine->Value(aNbPnts-2).ParametersOnS2(aU1, aV);
1870 }
1871 else
1872 {
1873 theLine->Value(aNbPnts).ParametersOnS1(aU3, aV);
1874 theLine->Value(aNbPnts-1).ParametersOnS1(aU2, aV);
1875 theLine->Value(aNbPnts-2).ParametersOnS1(aU1, aV);
1876 }
1877
1878 const Standard_Real aStepPrev = aU2-aU1;
1879 const Standard_Real aStep = aU3-aU2;
1880
1881 const Standard_Integer aDeltaStep = RealToInt(aStepPrev/aStep);
1882
1883 if((1 < aDeltaStep) && (aDeltaStep < 2000))
1884 {
1885 SeekAdditionalPoints( theQuad1, theQuad2, theLine,
1886 theCoeffs, theWLIndex, aDeltaStep, aNbPnts-2,
d30895f5 1887 aNbPnts-1, theTol2D, thePeriodOfSurf1, isTheReverse);
b5ef9d91 1888 }
1889 }
1890
ecc4f148 1891 return Standard_True;
1892}
1893
1894//=======================================================================
1895//function : AddBoundaryPoint
1896//purpose :
1897//=======================================================================
d30895f5 1898void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL,
ecc4f148 1899 const Standard_Real theU1,
1900 const Standard_Real theU2,
1901 const Standard_Real theV1,
1902 const Standard_Real theV1Prev,
1903 const Standard_Real theV2,
1904 const Standard_Real theV2Prev,
b5ef9d91 1905 const Standard_Integer theWLIndex,
02effd35 1906 const Standard_Boolean theFlForce,
ecc4f148 1907 Standard_Boolean& isTheFound1,
d30895f5 1908 Standard_Boolean& isTheFound2) const
ecc4f148 1909{
1910 Standard_Real aUSurf1f = 0.0, //const
1911 aUSurf1l = 0.0,
1912 aVSurf1f = 0.0,
1913 aVSurf1l = 0.0;
1914 Standard_Real aUSurf2f = 0.0, //const
1915 aUSurf2l = 0.0,
1916 aVSurf2f = 0.0,
1917 aVSurf2l = 0.0;
1918
d30895f5 1919 myUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
1920 myUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
1921
1922 const Standard_Integer aSize = 4;
1923 const Standard_Real anArrVzad[aSize] = {aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l};
ecc4f148 1924
d30895f5 1925 StPInfo aUVPoint[aSize];
b5ef9d91 1926
d30895f5 1927 for(Standard_Integer anIDSurf = 0; anIDSurf < 4; anIDSurf+=2)
ecc4f148 1928 {
d30895f5 1929 const Standard_Real aVf = (anIDSurf == 0) ? theV1 : theV2,
1930 aVl = (anIDSurf == 0) ? theV1Prev : theV2Prev;
ecc4f148 1931
d30895f5 1932 const SearchBoundType aTS = (anIDSurf == 0) ? SearchV1 : SearchV2;
ecc4f148 1933
d30895f5 1934 for(Standard_Integer anIDBound = 0; anIDBound < 2; anIDBound++)
1935 {
ecc4f148 1936
d30895f5 1937 const Standard_Integer anIndex = anIDSurf+anIDBound;
b5ef9d91 1938
d30895f5 1939 aUVPoint[anIndex].mySurfID = anIDSurf;
7365fad6 1940
d30895f5 1941 if((Abs(aVf-anArrVzad[anIndex]) > myTol2D) &&
1942 ((aVf-anArrVzad[anIndex])*(aVl-anArrVzad[anIndex]) > 0.0))
ecc4f148 1943 {
d30895f5 1944 continue;
ecc4f148 1945 }
ecc4f148 1946
d30895f5 1947 const Standard_Boolean aRes = SearchOnVBounds(aTS, anArrVzad[anIndex],
1948 (anIDSurf == 0) ? theV2 : theV1,
1949 theU2, theU1,
1950 aUVPoint[anIndex].myU1);
1951
1952 if(!aRes || aUVPoint[anIndex].myU1 >= theU1)
ecc4f148 1953 {
d30895f5 1954 aUVPoint[anIndex].myU1 = RealLast();
1955 continue;
ecc4f148 1956 }
d30895f5 1957 else
ecc4f148 1958 {
d30895f5 1959 Standard_Real &aU1 = aUVPoint[anIndex].myU1,
1960 &aU2 = aUVPoint[anIndex].myU2,
1961 &aV1 = aUVPoint[anIndex].myV1,
1962 &aV2 = aUVPoint[anIndex].myV2;
1963 aU2 = theU2;
1964 aV1 = theV1;
1965 aV2 = theV2;
1966
1967 if(!ComputationMethods::
1968 CylCylComputeParameters(aU1, theWLIndex, myCoeffs, aU2, aV1, aV2))
1969 {
1970 aU1 = RealLast();
1971 continue;
1972 }
1973
1974 if(aTS == SearchV1)
1975 aV1 = anArrVzad[anIndex];
1976 else //if(aTS[anIndex] == SearchV2)
1977 aV2 = anArrVzad[anIndex];
ecc4f148 1978 }
1979 }
1980 }
b5ef9d91 1981
d30895f5 1982 std::sort(aUVPoint, aUVPoint+aSize);
1983
1984 isTheFound1 = isTheFound2 = Standard_False;
b5ef9d91 1985
d30895f5 1986 for(Standard_Integer i = 0; i < aSize; i++)
1987 {
1988 if(aUVPoint[i].myU1 == RealLast())
1989 break;
b5ef9d91 1990
d30895f5 1991 if(!AddPointIntoWL(myQuad1, myQuad2, myCoeffs, myIsReverse, Standard_False,
1992 gp_Pnt2d(aUVPoint[i].myU1, aUVPoint[i].myV1),
1993 gp_Pnt2d(aUVPoint[i].myU2, aUVPoint[i].myV2),
1994 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
1995 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, myPeriod,
1996 theWL->Curve(), theWLIndex, myTol3D, myTol2D, theFlForce))
7365fad6 1997 {
d30895f5 1998 continue;
7365fad6 1999 }
ecc4f148 2000
d30895f5 2001 if(aUVPoint[i].mySurfID == 0)
ecc4f148 2002 {
d30895f5 2003 isTheFound1 = Standard_True;
ecc4f148 2004 }
7365fad6 2005 else
2006 {
d30895f5 2007 isTheFound2 = Standard_True;
7365fad6 2008 }
ecc4f148 2009 }
7fd59977 2010}
2011
ecc4f148 2012//=======================================================================
2013//function : SeekAdditionalPoints
2014//purpose :
2015//=======================================================================
2016static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
2017 const IntSurf_Quadric& theQuad2,
b5ef9d91 2018 const Handle(IntSurf_LineOn2S)& theLine,
d30895f5 2019 const ComputationMethods::stCoeffsValue& theCoeffs,
b5ef9d91 2020 const Standard_Integer theWLIndex,
ecc4f148 2021 const Standard_Integer theMinNbPoints,
b5ef9d91 2022 const Standard_Integer theStartPointOnLine,
2023 const Standard_Integer theEndPointOnLine,
ecc4f148 2024 const Standard_Real theTol2D,
2025 const Standard_Real thePeriodOfSurf2,
ecc4f148 2026 const Standard_Boolean isTheReverse)
2027{
b5ef9d91 2028 if(theLine.IsNull())
2029 return;
2030
2031 Standard_Integer aNbPoints = theEndPointOnLine - theStartPointOnLine + 1;
ecc4f148 2032 if(aNbPoints >= theMinNbPoints)
2033 {
2034 return;
2035 }
2036
b5ef9d91 2037 Standard_Real aMinDeltaParam = theTol2D;
2038
2039 {
2040 Standard_Real u1 = 0.0, v1 = 0.0, u2 = 0.0, v2 = 0.0;
2041
2042 if(isTheReverse)
2043 {
2044 theLine->Value(theStartPointOnLine).ParametersOnS2(u1, v1);
2045 theLine->Value(theEndPointOnLine).ParametersOnS2(u2, v2);
2046 }
2047 else
2048 {
2049 theLine->Value(theStartPointOnLine).ParametersOnS1(u1, v1);
2050 theLine->Value(theEndPointOnLine).ParametersOnS1(u2, v2);
2051 }
2052
2053 aMinDeltaParam = Max(Abs(u2 - u1)/IntToReal(theMinNbPoints), aMinDeltaParam);
2054 }
2055
2056 Standard_Integer aLastPointIndex = theEndPointOnLine;
ecc4f148 2057 Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0;
2058
2059 Standard_Integer aNbPointsPrev = 0;
2060 while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
2061 {
2062 aNbPointsPrev = aNbPoints;
b5ef9d91 2063 for(Standard_Integer fp = theStartPointOnLine, lp = 0; fp < aLastPointIndex; fp = lp + 1)
ecc4f148 2064 {
02effd35 2065 Standard_Real U1f = 0.0, V1f = 0.0; //first point in 1st suraface
2066 Standard_Real U1l = 0.0, V1l = 0.0; //last point in 1st suraface
2067
2068 Standard_Real U2f = 0.0, V2f = 0.0; //first point in 2nd suraface
2069 Standard_Real U2l = 0.0, V2l = 0.0; //last point in 2nd suraface
ecc4f148 2070
2071 lp = fp+1;
2072
2073 if(isTheReverse)
2074 {
b5ef9d91 2075 theLine->Value(fp).ParametersOnS2(U1f, V1f);
2076 theLine->Value(lp).ParametersOnS2(U1l, V1l);
02effd35 2077
b5ef9d91 2078 theLine->Value(fp).ParametersOnS1(U2f, V2f);
2079 theLine->Value(lp).ParametersOnS1(U2l, V2l);
ecc4f148 2080 }
2081 else
2082 {
b5ef9d91 2083 theLine->Value(fp).ParametersOnS1(U1f, V1f);
2084 theLine->Value(lp).ParametersOnS1(U1l, V1l);
02effd35 2085
b5ef9d91 2086 theLine->Value(fp).ParametersOnS2(U2f, V2f);
2087 theLine->Value(lp).ParametersOnS2(U2l, V2l);
02effd35 2088 }
2089
b5ef9d91 2090 if(Abs(U1l - U1f) <= aMinDeltaParam)
02effd35 2091 {
2092 //Step is minimal. It is not necessary to divide it.
2093 continue;
ecc4f148 2094 }
2095
2096 U1prec = 0.5*(U1f+U1l);
2097
d30895f5 2098 if(!ComputationMethods::
2099 CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec))
2100 {
b5ef9d91 2101 continue;
d30895f5 2102 }
ecc4f148 2103
d30895f5 2104 MinMax(U2f, U2l);
2105 if(!InscribePoint(U2f, U2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False))
2106 {
2107 continue;
2108 }
ecc4f148 2109
b5ef9d91 2110 const gp_Pnt aP1(theQuad1.Value(U1prec, V1prec)), aP2(theQuad2.Value(U2prec, V2prec));
2111 const gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
ecc4f148 2112
d30895f5 2113#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
2114 std::cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << std::endl;
7c32c7c4 2115#endif
2116
ecc4f148 2117 IntSurf_PntOn2S anIP;
2118 if(isTheReverse)
2119 {
2120 anIP.SetValue(aPInt, U2prec, V2prec, U1prec, V1prec);
2121 }
2122 else
2123 {
2124 anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec);
2125 }
2126
b5ef9d91 2127 theLine->InsertBefore(lp, anIP);
ecc4f148 2128
b5ef9d91 2129 aNbPoints++;
2130 aLastPointIndex++;
2131 }
2132
2133 if(aNbPoints >= theMinNbPoints)
2134 {
2135 return;
ecc4f148 2136 }
2137 }
2138}
2139
79997052 2140//=======================================================================
2141//function : BoundariesComputing
2142//purpose : Computes true domain of future intersection curve (allows
2143// avoiding excess iterations)
2144//=======================================================================
d30895f5 2145Standard_Boolean WorkWithBoundaries::BoundariesComputing(Standard_Real theU1f[],
2146 Standard_Real theU1l[]) const
79997052 2147{
d30895f5 2148 if(myCoeffs.mB > 0.0)
79997052 2149 {
d30895f5 2150 if(myCoeffs.mB + Abs(myCoeffs.mC) < -1.0)
79997052 2151 {//There is NOT intersection
2152 return Standard_False;
2153 }
d30895f5 2154 else if(myCoeffs.mB + Abs(myCoeffs.mC) <= 1.0)
79997052 2155 {//U=[0;2*PI]+aFI1
d30895f5 2156 theU1f[0] = myCoeffs.mFI1;
2157 theU1l[0] = myPeriod + myCoeffs.mFI1;
79997052 2158 }
d30895f5 2159 else if((1 + myCoeffs.mC <= myCoeffs.mB) &&
2160 (myCoeffs.mB <= 1 - myCoeffs.mC))
79997052 2161 {
d30895f5 2162 Standard_Real anArg = -(myCoeffs.mC + 1) / myCoeffs.mB;
79997052 2163 if(anArg > 1.0)
2164 anArg = 1.0;
2165 if(anArg < -1.0)
2166 anArg = -1.0;
2167
2168 const Standard_Real aDAngle = acos(anArg);
2169 //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
d30895f5 2170 theU1f[0] = myCoeffs.mFI1;
2171 theU1l[0] = aDAngle + myCoeffs.mFI1;
2172 theU1f[1] = myPeriod - aDAngle + myCoeffs.mFI1;
2173 theU1l[1] = myPeriod + myCoeffs.mFI1;
79997052 2174 }
d30895f5 2175 else if((1 - myCoeffs.mC <= myCoeffs.mB) &&
2176 (myCoeffs.mB <= 1 + myCoeffs.mC))
79997052 2177 {
d30895f5 2178 Standard_Real anArg = (1 - myCoeffs.mC) / myCoeffs.mB;
79997052 2179 if(anArg > 1.0)
2180 anArg = 1.0;
2181 if(anArg < -1.0)
2182 anArg = -1.0;
2183
2184 const Standard_Real aDAngle = acos(anArg);
2185 //U=[aDAngle;2*PI-aDAngle]+aFI1
2186
d30895f5 2187 theU1f[0] = aDAngle + myCoeffs.mFI1;
2188 theU1l[0] = myPeriod - aDAngle + myCoeffs.mFI1;
79997052 2189 }
d30895f5 2190 else if(myCoeffs.mB - Abs(myCoeffs.mC) >= 1.0)
79997052 2191 {
d30895f5 2192 Standard_Real anArg1 = (1 - myCoeffs.mC) / myCoeffs.mB,
2193 anArg2 = -(myCoeffs.mC + 1) / myCoeffs.mB;
79997052 2194 if(anArg1 > 1.0)
2195 anArg1 = 1.0;
2196 if(anArg1 < -1.0)
2197 anArg1 = -1.0;
2198
2199 if(anArg2 > 1.0)
2200 anArg2 = 1.0;
2201 if(anArg2 < -1.0)
2202 anArg2 = -1.0;
2203
2204 const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
2205 //(U=[aDAngle1;aDAngle2]+aFI1) ||
2206 //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
2207
d30895f5 2208 theU1f[0] = aDAngle1 + myCoeffs.mFI1;
2209 theU1l[0] = aDAngle2 + myCoeffs.mFI1;
2210 theU1f[1] = myPeriod - aDAngle2 + myCoeffs.mFI1;
2211 theU1l[1] = myPeriod - aDAngle1 + myCoeffs.mFI1;
79997052 2212 }
2213 else
2214 {
2215 return Standard_False;
2216 }
2217 }
d30895f5 2218 else if(myCoeffs.mB < 0.0)
79997052 2219 {
d30895f5 2220 if(myCoeffs.mB + Abs(myCoeffs.mC) > 1.0)
79997052 2221 {//There is NOT intersection
2222 return Standard_False;
2223 }
d30895f5 2224 else if(-myCoeffs.mB + Abs(myCoeffs.mC) <= 1.0)
79997052 2225 {//U=[0;2*PI]+aFI1
d30895f5 2226 theU1f[0] = myCoeffs.mFI1;
2227 theU1l[0] = myPeriod + myCoeffs.mFI1;
79997052 2228 }
d30895f5 2229 else if((-myCoeffs.mC - 1 <= myCoeffs.mB) &&
2230 ( myCoeffs.mB <= myCoeffs.mC - 1))
79997052 2231 {
d30895f5 2232 Standard_Real anArg = (1 - myCoeffs.mC) / myCoeffs.mB;
79997052 2233 if(anArg > 1.0)
2234 anArg = 1.0;
2235 if(anArg < -1.0)
2236 anArg = -1.0;
2237
2238 const Standard_Real aDAngle = acos(anArg);
2239 //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
2240
d30895f5 2241 theU1f[0] = myCoeffs.mFI1;
2242 theU1l[0] = aDAngle + myCoeffs.mFI1;
2243 theU1f[1] = myPeriod - aDAngle + myCoeffs.mFI1;
2244 theU1l[1] = myPeriod + myCoeffs.mFI1;
79997052 2245 }
d30895f5 2246 else if((myCoeffs.mC - 1 <= myCoeffs.mB) &&
2247 (myCoeffs.mB <= -myCoeffs.mB - 1))
79997052 2248 {
d30895f5 2249 Standard_Real anArg = -(myCoeffs.mC + 1) / myCoeffs.mB;
79997052 2250 if(anArg > 1.0)
2251 anArg = 1.0;
2252 if(anArg < -1.0)
2253 anArg = -1.0;
2254
2255 const Standard_Real aDAngle = acos(anArg);
2256 //U=[aDAngle;2*PI-aDAngle]+aFI1
2257
d30895f5 2258 theU1f[0] = aDAngle + myCoeffs.mFI1;
2259 theU1l[0] = myPeriod - aDAngle + myCoeffs.mFI1;
79997052 2260 }
d30895f5 2261 else if(-myCoeffs.mB - Abs(myCoeffs.mC) >= 1.0)
79997052 2262 {
d30895f5 2263 Standard_Real anArg1 = -(myCoeffs.mC + 1) / myCoeffs.mB,
2264 anArg2 = (1 - myCoeffs.mC) / myCoeffs.mB;
79997052 2265 if(anArg1 > 1.0)
2266 anArg1 = 1.0;
2267 if(anArg1 < -1.0)
2268 anArg1 = -1.0;
2269
2270 if(anArg2 > 1.0)
2271 anArg2 = 1.0;
2272 if(anArg2 < -1.0)
2273 anArg2 = -1.0;
2274
2275 const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
2276 //(U=[aDAngle1;aDAngle2]+aFI1) ||
2277 //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
2278
d30895f5 2279 theU1f[0] = aDAngle1 + myCoeffs.mFI1;
2280 theU1l[0] = aDAngle2 + myCoeffs.mFI1;
2281 theU1f[1] = myPeriod - aDAngle2 + myCoeffs.mFI1;
2282 theU1l[1] = myPeriod - aDAngle1 + myCoeffs.mFI1;
79997052 2283 }
2284 else
2285 {
2286 return Standard_False;
2287 }
2288 }
2289 else
2290 {
2291 return Standard_False;
2292 }
2293
2294 return Standard_True;
2295}
2296
ecc4f148 2297//=======================================================================
2298//function : CriticalPointsComputing
79997052 2299//purpose : theNbCritPointsMax contains true number of critical points.
2300// It must be initialized correctly before function calling
ecc4f148 2301//=======================================================================
d30895f5 2302static void CriticalPointsComputing(const ComputationMethods::stCoeffsValue& theCoeffs,
ecc4f148 2303 const Standard_Real theUSurf1f,
2304 const Standard_Real theUSurf1l,
2305 const Standard_Real theUSurf2f,
2306 const Standard_Real theUSurf2l,
2307 const Standard_Real thePeriod,
2308 const Standard_Real theTol2D,
7365fad6 2309 Standard_Integer& theNbCritPointsMax,
ecc4f148 2310 Standard_Real theU1crit[])
2311{
d30895f5 2312 //[0...1] - in these points parameter U1 goes through
2313 // the seam-edge of the first cylinder.
2314 //[2...3] - First and last U1 parameter.
2315 //[4...5] - in these points parameter U2 goes through
2316 // the seam-edge of the second cylinder.
2317 //[6...9] - in these points an intersection line goes through
2318 // U-boundaries of the second surface.
79997052 2319 //[10...11] - Boundary of monotonicity interval of U2(U1) function
2320 // (see CylCylMonotonicity() function)
7365fad6 2321
ecc4f148 2322 theU1crit[0] = 0.0;
2323 theU1crit[1] = thePeriod;
2324 theU1crit[2] = theUSurf1f;
2325 theU1crit[3] = theUSurf1l;
2326
2327 const Standard_Real aCOS = cos(theCoeffs.mFI2);
2328 const Standard_Real aBSB = Abs(theCoeffs.mB);
2329 if((theCoeffs.mC - aBSB <= aCOS) && (aCOS <= theCoeffs.mC + aBSB))
2330 {
2331 Standard_Real anArg = (aCOS - theCoeffs.mC) / theCoeffs.mB;
2332 if(anArg > 1.0)
2333 anArg = 1.0;
2334 if(anArg < -1.0)
2335 anArg = -1.0;
2336
2337 theU1crit[4] = -acos(anArg) + theCoeffs.mFI1;
7365fad6 2338 theU1crit[5] = acos(anArg) + theCoeffs.mFI1;
ecc4f148 2339 }
2340
2341 Standard_Real aSf = cos(theUSurf2f - theCoeffs.mFI2);
2342 Standard_Real aSl = cos(theUSurf2l - theCoeffs.mFI2);
2343 MinMax(aSf, aSl);
2344
7365fad6 2345 //In accorance with pure mathematic, theU1crit[6] and [8]
2346 //must be -Precision::Infinite() instead of used +Precision::Infinite()
2347 theU1crit[6] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
2348 -acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
2349 Precision::Infinite();
2350 theU1crit[7] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
2351 -acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
2352 Precision::Infinite();
2353 theU1crit[8] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
2354 acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
2355 Precision::Infinite();
2356 theU1crit[9] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
2357 acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
2358 Precision::Infinite();
2359
79997052 2360 theU1crit[10] = theCoeffs.mFI1;
2361 theU1crit[11] = M_PI+theCoeffs.mFI1;
2362
7365fad6 2363 //preparative treatment of array. This array must have faled to contain negative
2364 //infinity number
2365
2366 for(Standard_Integer i = 0; i < theNbCritPointsMax; i++)
2367 {
2368 if(Precision::IsInfinite(theU1crit[i]))
2369 {
2370 continue;
2371 }
2372
2373 theU1crit[i] = fmod(theU1crit[i], thePeriod);
2374 if(theU1crit[i] < 0.0)
2375 theU1crit[i] += thePeriod;
2376 }
2377
2378 //Here all not infinite elements of theU1crit are in [0, thePeriod) range
ecc4f148 2379
7365fad6 2380 do
ecc4f148 2381 {
7365fad6 2382 std::sort(theU1crit, theU1crit + theNbCritPointsMax);
2383 }
79997052 2384 while(ExcludeNearElements(theU1crit, theNbCritPointsMax,
2385 theUSurf1f, theUSurf1l, theTol2D));
7365fad6 2386
2387 //Here all not infinite elements in theU1crit are different and sorted.
2388
2389 while(theNbCritPointsMax > 0)
2390 {
2391 Standard_Real &anB = theU1crit[theNbCritPointsMax-1];
2392 if(Precision::IsInfinite(anB))
ecc4f148 2393 {
7365fad6 2394 theNbCritPointsMax--;
2395 continue;
ecc4f148 2396 }
7365fad6 2397
2398 //1st not infinte element is found
2399
2400 if(theNbCritPointsMax == 1)
2401 break;
2402
2403 //Here theNbCritPointsMax > 1
2404
2405 Standard_Real &anA = theU1crit[0];
2406
2407 //Compare 1st and last significant elements of theU1crit
2408 //They may still differs by period.
2409
2410 if (Abs(anB - anA - thePeriod) < theTol2D)
2411 {//E.g. anA == 2.0e-17, anB == (thePeriod-1.0e-18)
2412 anA = (anA + anB - thePeriod)/2.0;
2413 anB = Precision::Infinite();
2414 theNbCritPointsMax--;
2415 }
2416
2417 //Out of "while(theNbCritPointsMax > 0)" cycle.
2418 break;
ecc4f148 2419 }
7365fad6 2420
2421 //Attention! Here theU1crit may be unsorted.
ecc4f148 2422}
2423
2424//=======================================================================
d30895f5 2425//function : BoundaryEstimation
2426//purpose : Rough estimation of the parameter range.
2427//=======================================================================
2428void WorkWithBoundaries::BoundaryEstimation(const gp_Cylinder& theCy1,
2429 const gp_Cylinder& theCy2,
2430 Bnd_Range& theOutBoxS1,
2431 Bnd_Range& theOutBoxS2) const
2432{
2433 const gp_Dir &aD1 = theCy1.Axis().Direction(),
2434 &aD2 = theCy2.Axis().Direction();
2435 const Standard_Real aR1 = theCy1.Radius(),
2436 aR2 = theCy2.Radius();
2437
2438 //Let consider a parallelogram. Its edges are parallel to aD1 and aD2.
2439 //Its altitudes are equal to 2*aR1 and 2*aR2 (diameters of the cylinders).
2440 //In fact, this parallelogram is a projection of the cylinders to the plane
2441 //created by the intersected axes aD1 and aD2 (if the axes are skewed then
2442 //one axis can be translated by parallel shifting till intersection).
2443
2444 const Standard_Real aCosA = aD1.Dot(aD2);
2445 const Standard_Real aSqSinA = 1.0 - aCosA * aCosA;
2446
2447 //If sine is small then it can be compared with angle.
2448 if (aSqSinA < Precision::Angular()*Precision::Angular())
2449 return;
2450
2451 //Half of delta V. Delta V is a distance between
2452 //projections of two opposite parallelogram vertices
2453 //(joined by the maximal diagonal) to the cylinder axis.
2454 const Standard_Real aSinA = sqrt(aSqSinA);
2455 const Standard_Real aHDV1 = (aR1 * aCosA + aR2)/aSinA,
2456 aHDV2 = (aR2 * aCosA + aR1)/aSinA;
2457
2458#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
2459 //The code in this block is created for test only.It is stupidly to create
2460 //OCCT-test for the method, which will be changed possibly never.
2461 std::cout << "Reference values: aHDV1 = " << aHDV1 << "; aHDV2 = " << aHDV2 << std::endl;
2462#endif
2463
2464 //V-parameters of intersection point of the axes (in case of skewed axes,
2465 //see comment above).
2466 Standard_Real aV01 = 0.0, aV02 = 0.0;
2467 ExtremaLineLine(theCy1.Axis(), theCy2.Axis(), aCosA, aSqSinA, aV01, aV02);
2468
2469 theOutBoxS1.Add(aV01 - aHDV1);
2470 theOutBoxS1.Add(aV01 + aHDV1);
2471
2472 theOutBoxS2.Add(aV02 - aHDV2);
2473 theOutBoxS2.Add(aV02 + aHDV2);
2474
2475 theOutBoxS1.Enlarge(Precision::Confusion());
2476 theOutBoxS2.Enlarge(Precision::Confusion());
2477
2478 Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0;
2479 myUVSurf1.Get(aU1, aV1, aU2, aV2);
2480 theOutBoxS1.Common(Bnd_Range(aV1, aV2));
2481
2482 myUVSurf2.Get(aU1, aV1, aU2, aV2);
2483 theOutBoxS2.Common(Bnd_Range(aV1, aV2));
2484}
2485
2486//=======================================================================
2487//function : IntCyCy
ecc4f148 2488//purpose :
2489//=======================================================================
d30895f5 2490Standard_Boolean IntCyCy( const IntSurf_Quadric& theQuad1,
2491 const IntSurf_Quadric& theQuad2,
2492 const Standard_Real theTol3D,
2493 const Standard_Real theTol2D,
2494 const Bnd_Box2d& theUVSurf1,
2495 const Bnd_Box2d& theUVSurf2,
2496 const Standard_Boolean isTheReverse,
2497 Standard_Boolean& isTheEmpty,
2498 Standard_Boolean& isTheSameSurface,
2499 Standard_Boolean& isTheMultiplePoint,
2500 IntPatch_SequenceOfLine& theSlin,
2501 IntPatch_SequenceOfPoint& theSPnt)
ecc4f148 2502{
d30895f5 2503 isTheEmpty = Standard_True;
2504 isTheSameSurface = Standard_False;
2505 isTheMultiplePoint = Standard_False;
2506 theSlin.Clear();
2507 theSPnt.Clear();
2508
2509 const gp_Cylinder &aCyl1 = theQuad1.Cylinder(),
2510 &aCyl2 = theQuad2.Cylinder();
2511
2512 IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
2513
2514 if (!anInter.IsDone())
2515 {
2516 return Standard_False;
2517 }
2518
2519 if(anInter.TypeInter() != IntAna_NoGeometricSolution)
2520 {
2521 return CyCyAnalyticalIntersect( theQuad1, theQuad2, anInter,
2522 theTol3D, isTheReverse, isTheEmpty,
2523 isTheSameSurface, isTheMultiplePoint,
2524 theSlin, theSPnt);
2525 }
2526
ecc4f148 2527 Standard_Real aUSurf1f = 0.0, //const
2528 aUSurf1l = 0.0,
2529 aVSurf1f = 0.0,
2530 aVSurf1l = 0.0;
2531 Standard_Real aUSurf2f = 0.0, //const
2532 aUSurf2l = 0.0,
2533 aVSurf2f = 0.0,
2534 aVSurf2l = 0.0;
2535
2536 theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
2537 theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
2538
7c32c7c4 2539 const Standard_Integer aNbMaxPoints = 2000;
2540 const Standard_Integer aNbMinPoints = 200;
2541 const Standard_Integer aNbPoints = Min(Max(aNbMinPoints,
2542 RealToInt(20.0*aCyl1.Radius())), aNbMaxPoints);
ecc4f148 2543 const Standard_Real aPeriod = 2.0*M_PI;
2544 const Standard_Real aStepMin = theTol2D,
7365fad6 2545 aStepMax = (aUSurf1l-aUSurf1f > M_PI/100.0) ?
2546 (aUSurf1l-aUSurf1f)/IntToReal(aNbPoints) :
2547 aUSurf1l-aUSurf1f;
2548
7c32c7c4 2549 const Standard_Integer aNbWLines = 2;
2550
d30895f5 2551 const ComputationMethods::stCoeffsValue anEquationCoeffs(aCyl1, aCyl2);
2552
2553 const WorkWithBoundaries aBoundWork(theQuad1, theQuad2, anEquationCoeffs,
2554 theUVSurf1, theUVSurf2, aNbWLines,
2555 aPeriod, theTol3D, theTol2D, isTheReverse);
2556
2557 Bnd_Range aRangeS1, aRangeS2;
2558 aBoundWork.BoundaryEstimation(aCyl1, aCyl2, aRangeS1, aRangeS2);
2559 if (aRangeS1.IsVoid() || aRangeS2.IsVoid())
2560 return Standard_True;
ecc4f148 2561
2562 //Boundaries
2563 const Standard_Integer aNbOfBoundaries = 2;
2564 Standard_Real aU1f[aNbOfBoundaries] = {-Precision::Infinite(), -Precision::Infinite()};
2565 Standard_Real aU1l[aNbOfBoundaries] = {Precision::Infinite(), Precision::Infinite()};
2566
d30895f5 2567 if(!aBoundWork.BoundariesComputing(aU1f, aU1l))
79997052 2568 return Standard_True;
ecc4f148 2569
2570 for(Standard_Integer i = 0; i < aNbOfBoundaries; i++)
2571 {
2572 if(Precision::IsInfinite(aU1f[i]) && Precision::IsInfinite(aU1l[i]))
2573 continue;
2574
2575 InscribeInterval(aUSurf1f, aUSurf1l, aU1f[i], aU1l[i], theTol2D, aPeriod);
2576 }
2577
2578 if( !Precision::IsInfinite(aU1f[0]) && !Precision::IsInfinite(aU1f[1]) &&
2579 !Precision::IsInfinite(aU1l[0]) && !Precision::IsInfinite(aU1l[1]))
2580 {
2581 if( ((aU1f[1] <= aU1l[0]) || (aU1l[1] <= aU1l[0])) &&
2582 ((aU1f[0] <= aU1l[1]) || (aU1l[0] <= aU1l[1])))
2583 {//Join all intervals to one
2584 aU1f[0] = Min(aU1f[0], aU1f[1]);
2585 aU1l[0] = Max(aU1l[0], aU1l[1]);
2586
2587 aU1f[1] = -Precision::Infinite();
2588 aU1l[1] = Precision::Infinite();
2589 }
2590 }
2591
2592 //Critical points
79997052 2593 const Standard_Integer aNbCritPointsMax = 12;
ecc4f148 2594 Standard_Real anU1crit[aNbCritPointsMax] = {Precision::Infinite(),
2595 Precision::Infinite(),
2596 Precision::Infinite(),
2597 Precision::Infinite(),
2598 Precision::Infinite(),
2599 Precision::Infinite(),
2600 Precision::Infinite(),
2601 Precision::Infinite(),
2602 Precision::Infinite(),
79997052 2603 Precision::Infinite(),
2604 Precision::Infinite(),
ecc4f148 2605 Precision::Infinite()};
2606
7365fad6 2607 Standard_Integer aNbCritPoints = aNbCritPointsMax;
ecc4f148 2608 CriticalPointsComputing(anEquationCoeffs, aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
7365fad6 2609 aPeriod, theTol2D, aNbCritPoints, anU1crit);
ecc4f148 2610
2611 //Getting Walking-line
2612
b2af2f56 2613 enum WLFStatus
2614 {
2615 WLFStatus_Absent = 0,
2616 WLFStatus_Exist = 1,
2617 WLFStatus_Broken = 2
2618 };
2619
ecc4f148 2620 for(Standard_Integer aCurInterval = 0; aCurInterval < aNbOfBoundaries; aCurInterval++)
2621 {
2622 if(Precision::IsInfinite(aU1f[aCurInterval]) && Precision::IsInfinite(aU1l[aCurInterval]))
2623 continue;
2624
7c32c7c4 2625 Standard_Boolean isAddedIntoWL[aNbWLines];
2626 for(Standard_Integer i = 0; i < aNbWLines; i++)
2627 isAddedIntoWL[i] = Standard_False;
ecc4f148 2628
2629 Standard_Real anUf = aU1f[aCurInterval], anUl = aU1l[aCurInterval];
e8feb725 2630 const Standard_Boolean isDeltaPeriod = IsEqual(anUl-anUf, aPeriod);
ecc4f148 2631
2632 //Inscribe and sort critical points
7365fad6 2633 for(Standard_Integer i = 0; i < aNbCritPoints; i++)
2634 {
2635 InscribePoint(anUf, anUl, anU1crit[i], theTol2D, aPeriod, Standard_False);
2636 }
2637
2638 std::sort(anU1crit, anU1crit + aNbCritPoints);
ecc4f148 2639
2640 while(anUf < anUl)
2641 {
7c32c7c4 2642 Standard_Real aU2[aNbWLines], aV1[aNbWLines], aV2[aNbWLines];
b2af2f56 2643 WLFStatus aWLFindStatus[aNbWLines];
7c32c7c4 2644 Standard_Real aV1Prev[aNbWLines], aV2Prev[aNbWLines];
b5ef9d91 2645 Standard_Real anUexpect[aNbWLines];
7c32c7c4 2646 Standard_Boolean isAddingWLEnabled[aNbWLines];
2647
2648 Handle(IntSurf_LineOn2S) aL2S[aNbWLines];
2649 Handle(IntPatch_WLine) aWLine[aNbWLines];
2650 for(Standard_Integer i = 0; i < aNbWLines; i++)
2651 {
2652 aL2S[i] = new IntSurf_LineOn2S();
2653 aWLine[i] = new IntPatch_WLine(aL2S[i], Standard_False);
b2af2f56 2654 aWLFindStatus[i] = WLFStatus_Absent;
7c32c7c4 2655 isAddingWLEnabled[i] = Standard_True;
2656 aU2[i] = aV1[i] = aV2[i] = 0.0;
2657 aV1Prev[i] = aV2Prev[i] = 0.0;
b5ef9d91 2658 anUexpect[i] = anUf;
7c32c7c4 2659 }
2660
04e93070 2661 Standard_Real aCriticalDelta[aNbCritPointsMax] = {0};
7365fad6 2662 for(Standard_Integer aCritPID = 0; aCritPID < aNbCritPoints; aCritPID++)
2663 { //We are not intersted in elements of aCriticalDelta array
2664 //if their index is greater than or equal to aNbCritPoints
2665
2666 aCriticalDelta[aCritPID] = anUf - anU1crit[aCritPID];
2667 }
ecc4f148 2668
7365fad6 2669 Standard_Real anU1 = anUf;
ecc4f148 2670 Standard_Boolean isFirst = Standard_True;
2671
2672 while(anU1 <= anUl)
2673 {
7365fad6 2674 for(Standard_Integer i = 0; i < aNbCritPoints; i++)
ecc4f148 2675 {
2676 if((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0)
2677 {
2678 anU1 = anU1crit[i];
7c32c7c4 2679
51740958 2680 for(Standard_Integer j = 0; j < aNbWLines; j++)
b5ef9d91 2681 {
51740958 2682 aWLFindStatus[j] = WLFStatus_Broken;
2683 anUexpect[j] = anU1;
b5ef9d91 2684 }
7c32c7c4 2685
ecc4f148 2686 break;
2687 }
2688 }
2689
b5ef9d91 2690 if(IsEqual(anU1, anUl))
2691 {
2692 for(Standard_Integer i = 0; i < aNbWLines; i++)
2693 {
2694 aWLFindStatus[i] = WLFStatus_Broken;
2695 anUexpect[i] = anU1;
2696
2697 if(isDeltaPeriod)
2698 {
7365fad6 2699 //if isAddedIntoWL[i] == TRUE WLine contains only one point
b5ef9d91 2700 //(which was end point of previous WLine). If we will
2701 //add point found on the current step WLine will contain only
2702 //two points. At that both these points will be equal to the
2703 //points found earlier. Therefore, new WLine will repeat
2704 //already existing WLine. Consequently, it is necessary
2705 //to forbid building new line in this case.
2706
2707 isAddingWLEnabled[i] = (!isAddedIntoWL[i]);
2708 }
2709 else
2710 {
2711 isAddingWLEnabled[i] = ((theTol2D >= (anUexpect[i] - anU1)) ||
2712 (aWLFindStatus[i] == WLFStatus_Absent));
2713 }
2714 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2715 }
2716 else
2717 {
2718 for(Standard_Integer i = 0; i < aNbWLines; i++)
2719 {
2720 isAddingWLEnabled[i] = ((theTol2D >= (anUexpect[i] - anU1)) ||
2721 (aWLFindStatus[i] == WLFStatus_Absent));
2722 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2723 }
ecc4f148 2724
7c32c7c4 2725 for(Standard_Integer i = 0; i < aNbWLines; i++)
2726 {
2727 const Standard_Integer aNbPntsWL = aWLine[i].IsNull() ? 0 :
2728 aWLine[i]->Curve()->NbPoints();
7365fad6 2729
2730 if( (aWLFindStatus[i] == WLFStatus_Broken) ||
2731 (aWLFindStatus[i] == WLFStatus_Absent))
2732 {//Begin and end of WLine must be on boundary point
2733 //or on seam-edge strictly (if it is possible).
b5ef9d91 2734
7365fad6 2735 Standard_Real aTol = theTol2D;
d30895f5 2736 ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs,
2737 aU2[i], &aTol);
7365fad6 2738 InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False);
02effd35 2739
7365fad6 2740 aTol = Max(aTol, theTol2D);
02effd35 2741
7365fad6 2742 if(Abs(aU2[i]) <= aTol)
2743 aU2[i] = 0.0;
2744 else if(Abs(aU2[i] - aPeriod) <= aTol)
2745 aU2[i] = aPeriod;
2746 else if(Abs(aU2[i] - aUSurf2f) <= aTol)
2747 aU2[i] = aUSurf2f;
2748 else if(Abs(aU2[i] - aUSurf2l) <= aTol)
2749 aU2[i] = aUSurf2l;
2750 }
2751 else
2752 {
d30895f5 2753 ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]);
7365fad6 2754 InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False);
2755 }
2756
7c32c7c4 2757 if(aNbPntsWL == 0)
2758 {//the line has not contained any points yet
e002f1ce 2759 if(((aUSurf2f + aPeriod - aUSurf2l) <= 2.0*theTol2D) &&
7c32c7c4 2760 ((Abs(aU2[i] - aUSurf2f) < theTol2D) ||
2761 (Abs(aU2[i]-aUSurf2l) < theTol2D)))
02effd35 2762 {
e002f1ce 2763 //In this case aU2[i] can have two values: current aU2[i] or
b5ef9d91 2764 //aU2[i]+aPeriod (aU2[i]-aPeriod). It is necessary to choose
2765 //correct value.
e002f1ce 2766
b5ef9d91 2767 Standard_Boolean isIncreasing = Standard_True;
d30895f5 2768 ComputationMethods::CylCylMonotonicity(anU1+aStepMin, i, anEquationCoeffs,
2769 aPeriod, isIncreasing);
7c32c7c4 2770
e002f1ce 2771 //If U2(U1) is increasing and U2 is considered to be equal aUSurf2l
2772 //then after the next step (when U1 will be increased) U2 will be
2773 //increased too. And we will go out of surface boundary.
2774 //Therefore, If U2(U1) is increasing then U2 must be equal aUSurf2f.
2775 //Analogically, if U2(U1) is decreasing.
2776 if(isIncreasing)
2777 {
2778 aU2[i] = aUSurf2f;
2779 }
2780 else
2781 {
2782 aU2[i] = aUSurf2l;
7c32c7c4 2783 }
02effd35 2784 }
2785 }
e002f1ce 2786 else
2787 {//aNbPntsWL > 0
7c32c7c4 2788 if(((aUSurf2l - aUSurf2f) >= aPeriod) &&
2789 ((Abs(aU2[i] - aUSurf2f) < theTol2D) ||
2790 (Abs(aU2[i]-aUSurf2l) < theTol2D)))
2791 {//end of the line
2792 Standard_Real aU2prev = 0.0, aV2prev = 0.0;
2793 if(isTheReverse)
2794 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS1(aU2prev, aV2prev);
02effd35 2795 else
7c32c7c4 2796 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS2(aU2prev, aV2prev);
2797
2798 if(2.0*Abs(aU2prev - aU2[i]) > aPeriod)
2799 {
2800 if(aU2prev > aU2[i])
2801 aU2[i] += aPeriod;
2802 else
2803 aU2[i] -= aPeriod;
2804 }
02effd35 2805 }
2806 }
02effd35 2807
d30895f5 2808 ComputationMethods::CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs,
2809 aV1[i], aV2[i]);
7c32c7c4 2810
2811 if(isFirst)
2812 {
2813 aV1Prev[i] = aV1[i];
2814 aV2Prev[i] = aV2[i];
02effd35 2815 }
7c32c7c4 2816 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
ecc4f148 2817
7c32c7c4 2818 isFirst = Standard_False;
ecc4f148 2819
7c32c7c4 2820 //Looking for points into WLine
2821 Standard_Boolean isBroken = Standard_False;
2822 for(Standard_Integer i = 0; i < aNbWLines; i++)
ecc4f148 2823 {
7c32c7c4 2824 if(!isAddingWLEnabled[i])
ecc4f148 2825 {
b5ef9d91 2826 Standard_Boolean isBoundIntersect = Standard_False;
2827 if( (Abs(aV1[i] - aVSurf1f) <= theTol2D) ||
2828 ((aV1[i]-aVSurf1f)*(aV1Prev[i]-aVSurf1f) < 0.0))
2829 {
2830 isBoundIntersect = Standard_True;
2831 }
2832 else if( (Abs(aV1[i] - aVSurf1l) <= theTol2D) ||
2833 ( (aV1[i]-aVSurf1l)*(aV1Prev[i]-aVSurf1l) < 0.0))
2834 {
2835 isBoundIntersect = Standard_True;
2836 }
2837 else if( (Abs(aV2[i] - aVSurf2f) <= theTol2D) ||
2838 ( (aV2[i]-aVSurf2f)*(aV2Prev[i]-aVSurf2f) < 0.0))
2839 {
2840 isBoundIntersect = Standard_True;
2841 }
2842 else if( (Abs(aV2[i] - aVSurf2l) <= theTol2D) ||
2843 ( (aV2[i]-aVSurf2l)*(aV2Prev[i]-aVSurf2l) < 0.0))
2844 {
2845 isBoundIntersect = Standard_True;
2846 }
2847
b2af2f56 2848 if(aWLFindStatus[i] == WLFStatus_Broken)
7c32c7c4 2849 isBroken = Standard_True;
02effd35 2850
b5ef9d91 2851 if(!isBoundIntersect)
2852 {
2853 continue;
2854 }
2855 else
2856 {
2857 anUexpect[i] = anU1;
2858 }
ecc4f148 2859 }
ecc4f148 2860
b5ef9d91 2861 const Standard_Boolean isInscribe =
2862 ((aUSurf2f-aU2[i]) <= theTol2D) && ((aU2[i]-aUSurf2l) <= theTol2D) &&
7c32c7c4 2863 ((aVSurf1f - aV1[i]) <= theTol2D) && ((aV1[i] - aVSurf1l) <= theTol2D) &&
b5ef9d91 2864 ((aVSurf2f - aV2[i]) <= theTol2D) && ((aV2[i] - aVSurf2l) <= theTol2D);
2865
2866 //isVIntersect == TRUE if intersection line intersects two (!)
2867 //V-bounds of cylinder (1st or 2nd - no matter)
2868 const Standard_Boolean isVIntersect =
2869 ( ((aVSurf1f-aV1[i])*(aVSurf1f-aV1Prev[i]) < RealSmall()) &&
2870 ((aVSurf1l-aV1[i])*(aVSurf1l-aV1Prev[i]) < RealSmall())) ||
2871 ( ((aVSurf2f-aV2[i])*(aVSurf2f-aV2Prev[i]) < RealSmall()) &&
2872 ((aVSurf2l-aV2[i])*(aVSurf2l-aV2Prev[i]) < RealSmall()));
2873
2874
2875 //isFound1 == TRUE if intersection line intersects V-bounds
2876 // (First or Last - no matter) of the 1st cylynder
2877 //isFound2 == TRUE if intersection line intersects V-bounds
2878 // (First or Last - no matter) of the 2nd cylynder
2879 Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
2880 Standard_Boolean isForce = Standard_False;
2881
5b111128 2882 if (aWLFindStatus[i] == WLFStatus_Absent)
ecc4f148 2883 {
b5ef9d91 2884 if(((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1-aUSurf1l) < theTol2D))
02effd35 2885 {
b5ef9d91 2886 isForce = Standard_True;
2887 }
2888 }
e8feb725 2889
d30895f5 2890 aBoundWork.AddBoundaryPoint(aWLine[i], anU1, aU2[i], aV1[i], aV1Prev[i],
2891 aV2[i], aV2Prev[i], i, isForce,
2892 isFound1, isFound2);
02effd35 2893
b5ef9d91 2894 const Standard_Boolean isPrevVBound = !isVIntersect &&
2895 ((Abs(aV1Prev[i] - aVSurf1f) <= theTol2D) ||
2896 (Abs(aV1Prev[i] - aVSurf1l) <= theTol2D) ||
2897 (Abs(aV2Prev[i] - aVSurf2f) <= theTol2D) ||
2898 (Abs(aV2Prev[i] - aVSurf2l) <= theTol2D));
2899
7365fad6 2900 aV1Prev[i] = aV1[i];
2901 aV2Prev[i] = aV2[i];
2902
b5ef9d91 2903 if((aWLFindStatus[i] == WLFStatus_Exist) && (isFound1 || isFound2) && !isPrevVBound)
2904 {
2905 aWLFindStatus[i] = WLFStatus_Broken; //start a new line
2906 }
2907 else if(isInscribe)
2908 {
2909 if((aWLFindStatus[i] == WLFStatus_Absent) && (isFound1 || isFound2))
2910 {
2911 aWLFindStatus[i] = WLFStatus_Exist;
ecc4f148 2912 }
ecc4f148 2913
d30895f5 2914 if( (aWLFindStatus[i] != WLFStatus_Broken) ||
2915 (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl))
ecc4f148 2916 {
5b9e1842 2917 if(aWLine[i]->NbPnts() > 0)
2918 {
2919 Standard_Real aU2p = 0.0, aV2p = 0.0;
2920 if(isTheReverse)
2921 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
2922 else
2923 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
2924
2925 const Standard_Real aDelta = aU2[i] - aU2p;
2926
2927 if(2*Abs(aDelta) > aPeriod)
2928 {
2929 if(aDelta > 0.0)
2930 {
2931 aU2[i] -= aPeriod;
2932 }
2933 else
2934 {
2935 aU2[i] += aPeriod;
2936 }
2937 }
2938 }
2939
b5ef9d91 2940 if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse, Standard_True,
7c32c7c4 2941 gp_Pnt2d(anU1, aV1[i]), gp_Pnt2d(aU2[i], aV2[i]),
b5ef9d91 2942 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
2943 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod,
2944 aWLine[i]->Curve(), i, theTol3D, theTol2D, isForce))
ecc4f148 2945 {
b2af2f56 2946 if(aWLFindStatus[i] == WLFStatus_Absent)
e8feb725 2947 {
b2af2f56 2948 aWLFindStatus[i] = WLFStatus_Exist;
e8feb725 2949 }
ecc4f148 2950 }
5b9e1842 2951 else if(!isFound1 && !isFound2)
2952 {//We do not add any point while doing this iteration
2953 if(aWLFindStatus[i] == WLFStatus_Exist)
2954 {
2955 aWLFindStatus[i] = WLFStatus_Broken;
2956 }
2957 }
2958 }
2959 }
2960 else
2961 {//We do not add any point while doing this iteration
2962 if(aWLFindStatus[i] == WLFStatus_Exist)
2963 {
2964 aWLFindStatus[i] = WLFStatus_Broken;
ecc4f148 2965 }
2966 }
b5ef9d91 2967
b2af2f56 2968 if(aWLFindStatus[i] == WLFStatus_Broken)
7c32c7c4 2969 isBroken = Standard_True;
2970 }//for(Standard_Integer i = 0; i < aNbWLines; i++)
ecc4f148 2971
7c32c7c4 2972 if(isBroken)
ecc4f148 2973 {//current lines are filled. Go to the next lines
2974 anUf = anU1;
b5ef9d91 2975
5b9e1842 2976 Standard_Boolean isAdded = Standard_True;
2977
b5ef9d91 2978 for(Standard_Integer i = 0; i < aNbWLines; i++)
2979 {
5b9e1842 2980 if(isAddingWLEnabled[i])
2981 {
b5ef9d91 2982 continue;
5b9e1842 2983 }
b5ef9d91 2984
5b9e1842 2985 isAdded = Standard_False;
b5ef9d91 2986
5b9e1842 2987 Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
b5ef9d91 2988
d30895f5 2989 aBoundWork.AddBoundaryPoint(aWLine[i], anU1, aU2[i], aV1[i], aV1Prev[i],
2990 aV2[i], aV2Prev[i], i,
2991 Standard_False, isFound1, isFound2);
b5ef9d91 2992
5b9e1842 2993 if(isFound1 || isFound2)
2994 {
2995 isAdded = Standard_True;
2996 }
b5ef9d91 2997
5b9e1842 2998 if(aWLine[i]->NbPnts() > 0)
b5ef9d91 2999 {
5b9e1842 3000 Standard_Real aU2p = 0.0, aV2p = 0.0;
3001 if(isTheReverse)
3002 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
3003 else
3004 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
3005
3006 const Standard_Real aDelta = aU2[i] - aU2p;
3007
3008 if(2*Abs(aDelta) > aPeriod)
3009 {
3010 if(aDelta > 0.0)
3011 {
3012 aU2[i] -= aPeriod;
3013 }
3014 else
3015 {
3016 aU2[i] += aPeriod;
3017 }
3018 }
3019 }
b5ef9d91 3020
5b9e1842 3021 if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse,
3022 Standard_True, gp_Pnt2d(anU1, aV1[i]),
3023 gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
3024 aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
3025 aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
3026 i, theTol3D, theTol2D, Standard_False))
3027 {
3028 isAdded = Standard_True;
b5ef9d91 3029 }
5b9e1842 3030 }
b5ef9d91 3031
5b9e1842 3032 if(!isAdded)
3033 {
3034 Standard_Real anUmaxAdded = RealFirst();
7365fad6 3035
5b9e1842 3036 {
7365fad6 3037 Standard_Boolean isChanged = Standard_False;
3038 for(Standard_Integer i = 0; i < aNbWLines; i++)
3039 {
3040 if(aWLFindStatus[i] == WLFStatus_Absent)
3041 continue;
b5ef9d91 3042
7365fad6 3043 Standard_Real aU1c = 0.0, aV1c = 0.0;
3044 if(isTheReverse)
3045 aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS2(aU1c, aV1c);
3046 else
3047 aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS1(aU1c, aV1c);
3048
3049 anUmaxAdded = Max(anUmaxAdded, aU1c);
3050 isChanged = Standard_True;
3051 }
3052
3053 if(!isChanged)
3054 { //If anUmaxAdded were not changed in previous cycle then
3055 //we would break existing WLines.
3056 break;
3057 }
b5ef9d91 3058 }
3059
5b9e1842 3060 for(Standard_Integer i = 0; i < aNbWLines; i++)
3061 {
3062 if(isAddingWLEnabled[i])
3063 {
3064 continue;
3065 }
b5ef9d91 3066
d30895f5 3067 ComputationMethods::CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs,
3068 aU2[i], aV1[i], aV2[i]);
b5ef9d91 3069
5b9e1842 3070 AddPointIntoWL( theQuad1, theQuad2, anEquationCoeffs, isTheReverse,
3071 Standard_True, gp_Pnt2d(anUmaxAdded, aV1[i]),
3072 gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
3073 aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
3074 aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
3075 i, theTol3D, theTol2D, Standard_False);
3076 }
b5ef9d91 3077 }
3078
ecc4f148 3079 break;
3080 }
3081
7c32c7c4 3082 //Step computing
ecc4f148 3083
7c32c7c4 3084 {
d30895f5 3085 const Standard_Real aDeltaV1 = aRangeS1.Delta()/IntToReal(aNbPoints),
3086 aDeltaV2 = aRangeS2.Delta()/IntToReal(aNbPoints);
b5ef9d91 3087
3088 math_Matrix aMatr(1, 3, 1, 5);
ecc4f148 3089
b5ef9d91 3090 Standard_Real aMinUexp = RealLast();
3091
3092 for(Standard_Integer i = 0; i < aNbWLines; i++)
3093 {
3094 if(theTol2D < (anUexpect[i] - anU1))
3095 {
3096 continue;
3097 }
ecc4f148 3098
b5ef9d91 3099 if(aWLFindStatus[i] == WLFStatus_Absent)
3100 {
3101 anUexpect[i] += aStepMax;
3102 aMinUexp = Min(aMinUexp, anUexpect[i]);
3103 continue;
3104 }
7c32c7c4 3105
b5ef9d91 3106 Standard_Real aStepTmp = aStepMax;
ecc4f148 3107
b5ef9d91 3108 const Standard_Real aSinU1 = sin(anU1),
3109 aCosU1 = cos(anU1),
3110 aSinU2 = sin(aU2[i]),
3111 aCosU2 = cos(aU2[i]);
ecc4f148 3112
b5ef9d91 3113 aMatr.SetCol(1, anEquationCoeffs.mVecC1);
3114 aMatr.SetCol(2, anEquationCoeffs.mVecC2);
3115 aMatr.SetCol(3, anEquationCoeffs.mVecA1*aSinU1 - anEquationCoeffs.mVecB1*aCosU1);
3116 aMatr.SetCol(4, anEquationCoeffs.mVecA2*aSinU2 - anEquationCoeffs.mVecB2*aCosU2);
3117 aMatr.SetCol(5, anEquationCoeffs.mVecA1*aCosU1 + anEquationCoeffs.mVecB1*aSinU1 +
3118 anEquationCoeffs.mVecA2*aCosU2 + anEquationCoeffs.mVecB2*aSinU2 +
3119 anEquationCoeffs.mVecD);
ecc4f148 3120
7365fad6 3121 if(!StepComputing(aMatr, aV1[i], aV2[i], aDeltaV1, aDeltaV2, aStepTmp))
b5ef9d91 3122 {
3123 //To avoid cycling-up
3124 anUexpect[i] += aStepMax;
3125 aMinUexp = Min(aMinUexp, anUexpect[i]);
ecc4f148 3126
b5ef9d91 3127 continue;
3128 }
3129
3130 if(aStepTmp < aStepMin)
3131 aStepTmp = aStepMin;
ecc4f148 3132
b5ef9d91 3133 if(aStepTmp > aStepMax)
3134 aStepTmp = aStepMax;
3135
3136 anUexpect[i] = anU1 + aStepTmp;
3137 aMinUexp = Min(aMinUexp, anUexpect[i]);
3138 }
ecc4f148 3139
b5ef9d91 3140 anU1 = aMinUexp;
3141 }
ecc4f148 3142
b5ef9d91 3143 if(Precision::PConfusion() >= (anUl - anU1))
ecc4f148 3144 anU1 = anUl;
3145
3146 anUf = anU1;
3147
7c32c7c4 3148 for(Standard_Integer i = 0; i < aNbWLines; i++)
baf72cd2 3149 {
7c32c7c4 3150 if(aWLine[i]->NbPnts() != 1)
3151 isAddedIntoWL[i] = Standard_False;
b5ef9d91 3152
3153 if(anU1 == anUl)
3154 {//strictly equal. Tolerance is considered above.
3155 anUexpect[i] = anUl;
3156 }
baf72cd2 3157 }
ecc4f148 3158 }
ecc4f148 3159
7c32c7c4 3160 for(Standard_Integer i = 0; i < aNbWLines; i++)
ecc4f148 3161 {
7c32c7c4 3162 if((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
baf72cd2 3163 {
7c32c7c4 3164 isTheEmpty = Standard_False;
3165 Standard_Real u1, v1, u2, v2;
3166 aWLine[i]->Point(1).Parameters(u1, v1, u2, v2);
3167 IntPatch_Point aP;
3168 aP.SetParameter(u1);
3169 aP.SetParameters(u1, v1, u2, v2);
3170 aP.SetTolerance(theTol3D);
3171 aP.SetValue(aWLine[i]->Point(1).Value());
3172
3173 theSPnt.Append(aP);
baf72cd2 3174 }
7c32c7c4 3175 else if(aWLine[i]->NbPnts() > 1)
baf72cd2 3176 {
7c32c7c4 3177 Standard_Boolean isGood = Standard_True;
ecc4f148 3178
7c32c7c4 3179 if(aWLine[i]->NbPnts() == 2)
3180 {
3181 const IntSurf_PntOn2S& aPf = aWLine[i]->Point(1);
3182 const IntSurf_PntOn2S& aPl = aWLine[i]->Point(2);
3183
3184 if(aPf.IsSame(aPl, Precision::Confusion()))
3185 isGood = Standard_False;
3186 }
ecc4f148 3187
7c32c7c4 3188 if(isGood)
3189 {
3190 isTheEmpty = Standard_False;
3191 isAddedIntoWL[i] = Standard_True;
3192 SeekAdditionalPoints( theQuad1, theQuad2, aWLine[i]->Curve(),
b5ef9d91 3193 anEquationCoeffs, i, aNbPoints, 1,
d30895f5 3194 aWLine[i]->NbPnts(), theTol2D, aPeriod,
3195 isTheReverse);
7c32c7c4 3196
3197 aWLine[i]->ComputeVertexParameters(theTol3D);
3198 theSlin.Append(aWLine[i]);
3199 }
3200 }
3201 else
3202 {
3203 isAddedIntoWL[i] = Standard_False;
baf72cd2 3204 }
b5ef9d91 3205
d30895f5 3206#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
b5ef9d91 3207 //aWLine[i]->Dump();
3208#endif
3209 }
3210 }
3211 }
3212
3213
7365fad6 3214 //Delete the points in theSPnt, which
3215 //lie at least in one of the line in theSlin.
b5ef9d91 3216 for(Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
3217 {
3218 for(Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); aNbLin++)
3219 {
d30895f5 3220 Handle(IntPatch_WLine) aWLine1 (Handle(IntPatch_WLine)::
3221 DownCast(theSlin.Value(aNbLin)));
b5ef9d91 3222
3223 const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1);
3224 const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aWLine1->NbPnts());
3225
3226 const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNbPnt).PntOn2S();
3227 if( aPntCur.IsSame(aPntFWL1, Precision::Confusion()) ||
3228 aPntCur.IsSame(aPntLWL1, Precision::Confusion()))
3229 {
3230 theSPnt.Remove(aNbPnt);
3231 aNbPnt--;
3232 break;
ecc4f148 3233 }
3234 }
3235 }
3236
b5ef9d91 3237 const Standard_Real aDU = aStepMin + Epsilon(aStepMin);
3238 //Try to add new points in the neighbourhood of existing point
3239 for(Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
3240 {
3241 const IntPatch_Point& aPnt2S = theSPnt.Value(aNbPnt);
3242
3243 Standard_Real u1, v1, u2, v2;
3244 aPnt2S.Parameters(u1, v1, u2, v2);
3245
3246 Handle(IntSurf_LineOn2S) aL2S = new IntSurf_LineOn2S();
3247 Handle(IntPatch_WLine) aWLine = new IntPatch_WLine(aL2S, Standard_False);
3248 aWLine->Curve()->Add(aPnt2S.PntOn2S());
3249
3250 //Define the index of WLine, which lies the point aPnt2S in.
3251 Standard_Real anUf = 0.0, anUl = 0.0, aCurU2 = 0.0;
3252 Standard_Integer anIndex = 0;
3253 if(isTheReverse)
3254 {
3255 anUf = Max(u2 - aStepMax, aUSurf1f);
3256 anUl = u2;
3257 aCurU2 = u1;
3258 }
3259 else
3260 {
3261 anUf = Max(u1 - aStepMax, aUSurf1f);
3262 anUl = u1;
3263 aCurU2 = u2;
3264 }
3265 Standard_Real aDelta = RealLast();
3266 for (Standard_Integer i = 0; i < aNbWLines; i++)
3267 {
3268 Standard_Real anU2t = 0.0;
d30895f5 3269 if(!ComputationMethods::CylCylComputeParameters(anUl, i, anEquationCoeffs, anU2t))
b5ef9d91 3270 continue;
3271
51740958 3272 const Standard_Real aDU2 = Abs(anU2t - aCurU2);
3273 if(aDU2 < aDelta)
b5ef9d91 3274 {
51740958 3275 aDelta = aDU2;
b5ef9d91 3276 anIndex = i;
3277 }
3278 }
3279
3280 //Try to fill aWLine by additional points
3281 while(anUl - anUf > RealSmall())
3282 {
3283 Standard_Real anU2 = 0.0, anV1 = 0.0, anV2 = 0.0;
3284 Standard_Boolean isDone =
d30895f5 3285 ComputationMethods::CylCylComputeParameters(anUf, anIndex, anEquationCoeffs,
b5ef9d91 3286 anU2, anV1, anV2);
3287 anUf += aDU;
5b9e1842 3288
b5ef9d91 3289 if(!isDone)
3290 {
3291 continue;
3292 }
3293
d30895f5 3294 if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse,
3295 Standard_True, gp_Pnt2d(anUf, anV1), gp_Pnt2d(anU2, anV2),
b5ef9d91 3296 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
3297 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l,
3298 aPeriod, aWLine->Curve(), anIndex, theTol3D,
3299 theTol2D, Standard_False, Standard_True))
3300 {
3301 break;
3302 }
3303 }
3304
3305 if(aWLine->NbPnts() > 1)
3306 {
3307 SeekAdditionalPoints( theQuad1, theQuad2, aWLine->Curve(),
3308 anEquationCoeffs, anIndex, aNbMinPoints,
d30895f5 3309 1, aWLine->NbPnts(), theTol2D, aPeriod,
3310 isTheReverse);
b5ef9d91 3311
3312 aWLine->ComputeVertexParameters(theTol3D);
3313 theSlin.Append(aWLine);
5b9e1842 3314
b5ef9d91 3315 theSPnt.Remove(aNbPnt);
3316 aNbPnt--;
3317 }
3318 }
3319
ecc4f148 3320 return Standard_True;
3321}
7fd59977 3322
3323//=======================================================================
3324//function : IntCySp
3325//purpose :
3326//=======================================================================
3327Standard_Boolean IntCySp(const IntSurf_Quadric& Quad1,
3328 const IntSurf_Quadric& Quad2,
3329 const Standard_Real Tol,
3330 const Standard_Boolean Reversed,
3331 Standard_Boolean& Empty,
3332 Standard_Boolean& Multpoint,
3333 IntPatch_SequenceOfLine& slin,
3334 IntPatch_SequenceOfPoint& spnt)
3335
3336{
3337 Standard_Integer i;
3338
3339 IntSurf_TypeTrans trans1,trans2;
3340 IntAna_ResultType typint;
3341 IntPatch_Point ptsol;
3342 gp_Circ cirsol;
3343
3344 gp_Cylinder Cy;
3345 gp_Sphere Sp;
3346
3347 if (!Reversed) {
3348 Cy = Quad1.Cylinder();
3349 Sp = Quad2.Sphere();
3350 }
3351 else {
3352 Cy = Quad2.Cylinder();
3353 Sp = Quad1.Sphere();
3354 }
3355 IntAna_QuadQuadGeo inter(Cy,Sp,Tol);
3356
3357 if (!inter.IsDone()) {return Standard_False;}
3358
3359 typint = inter.TypeInter();
3360 Standard_Integer NbSol = inter.NbSolutions();
3361 Empty = Standard_False;
3362
3363 switch (typint) {
3364
3365 case IntAna_Empty :
3366 {
3367 Empty = Standard_True;
3368 }
3369 break;
3370
3371 case IntAna_Point :
3372 {
3373 gp_Pnt psol(inter.Point(1));
3374 Standard_Real U1,V1,U2,V2;
3375 Quad1.Parameters(psol,U1,V1);
3376 Quad2.Parameters(psol,U2,V2);
3377 ptsol.SetValue(psol,Tol,Standard_True);
3378 ptsol.SetParameters(U1,V1,U2,V2);
3379 spnt.Append(ptsol);
3380 }
3381 break;
3382
3383 case IntAna_Circle:
3384 {
3385 cirsol = inter.Circle(1);
3386 gp_Vec Tgt;
3387 gp_Pnt ptref;
3388 ElCLib::D1(0.,cirsol,ptref,Tgt);
3389
3390 if (NbSol == 1) {
3391 gp_Vec TestCurvature(ptref,Sp.Location());
3392 gp_Vec Normsp,Normcyl;
3393 if (!Reversed) {
3394 Normcyl = Quad1.Normale(ptref);
3395 Normsp = Quad2.Normale(ptref);
3396 }
3397 else {
3398 Normcyl = Quad2.Normale(ptref);
3399 Normsp = Quad1.Normale(ptref);
3400 }
3401
3402 IntSurf_Situation situcyl;
3403 IntSurf_Situation situsp;
3404
3405 if (Normcyl.Dot(TestCurvature) > 0.) {
3406 situsp = IntSurf_Outside;
3407 if (Normsp.Dot(Normcyl) > 0.) {
3408 situcyl = IntSurf_Inside;
3409 }
3410 else {
3411 situcyl = IntSurf_Outside;
3412 }
3413 }
3414 else {
3415 situsp = IntSurf_Inside;
3416 if (Normsp.Dot(Normcyl) > 0.) {
3417 situcyl = IntSurf_Outside;
3418 }
3419 else {
3420 situcyl = IntSurf_Inside;
3421 }
3422 }
3423 Handle(IntPatch_GLine) glig;
3424 if (!Reversed) {
3425 glig = new IntPatch_GLine(cirsol, Standard_True, situcyl, situsp);
3426 }
3427 else {
3428 glig = new IntPatch_GLine(cirsol, Standard_True, situsp, situcyl);
3429 }
3430 slin.Append(glig);
3431 }
3432 else {
3433 if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
3434 trans1 = IntSurf_Out;
3435 trans2 = IntSurf_In;
3436 }
3437 else {
3438 trans1 = IntSurf_In;
3439 trans2 = IntSurf_Out;
3440 }
3441 Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3442 slin.Append(glig);
3443
3444 cirsol = inter.Circle(2);
3445 ElCLib::D1(0.,cirsol,ptref,Tgt);
3446 Standard_Real qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
3447 if(qwe> 0.0000001) {
3448 trans1 = IntSurf_Out;
3449 trans2 = IntSurf_In;
3450 }
3451 else if(qwe<-0.0000001) {
3452 trans1 = IntSurf_In;
3453 trans2 = IntSurf_Out;
3454 }
3455 else {
3456 trans1=trans2=IntSurf_Undecided;
3457 }
3458 glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3459 slin.Append(glig);
3460 }
3461 }
3462 break;
3463
3464 case IntAna_NoGeometricSolution:
3465 {
3466 gp_Pnt psol;
3467 Standard_Real U1,V1,U2,V2;
3468 IntAna_IntQuadQuad anaint(Cy,Sp,Tol);
3469 if (!anaint.IsDone()) {
3470 return Standard_False;
3471 }
3472
3473 if (anaint.NbPnt()==0 && anaint.NbCurve()==0) {
3474 Empty = Standard_True;
3475 }
3476 else {
3477
3478 NbSol = anaint.NbPnt();
3479 for (i = 1; i <= NbSol; i++) {
3480 psol = anaint.Point(i);
3481 Quad1.Parameters(psol,U1,V1);
3482 Quad2.Parameters(psol,U2,V2);
3483 ptsol.SetValue(psol,Tol,Standard_True);
3484 ptsol.SetParameters(U1,V1,U2,V2);
3485 spnt.Append(ptsol);
3486 }
3487
3488 gp_Pnt ptvalid,ptf,ptl;
3489 gp_Vec tgvalid;
3490 Standard_Real first,last,para;
3491 IntAna_Curve curvsol;
3492 Standard_Boolean tgfound;
3493 Standard_Integer kount;
3494
3495 NbSol = anaint.NbCurve();
3496 for (i = 1; i <= NbSol; i++) {
3497 curvsol = anaint.Curve(i);
3498 curvsol.Domain(first,last);
3499 ptf = curvsol.Value(first);
3500 ptl = curvsol.Value(last);
3501
3502 para = last;
3503 kount = 1;
3504 tgfound = Standard_False;
3505
3506 while (!tgfound) {
3507 para = (1.123*first + para)/2.123;
3508 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
3509 if (!tgfound) {
3510 kount ++;
3511 tgfound = kount > 5;
3512 }
3513 }
3514 Handle(IntPatch_ALine) alig;
3515 if (kount <= 5) {
3516 Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
3517 Quad1.Normale(ptvalid));
3518 if(qwe> 0.00000001) {
3519 trans1 = IntSurf_Out;
3520 trans2 = IntSurf_In;
3521 }
3522 else if(qwe<-0.00000001) {
3523 trans1 = IntSurf_In;
3524 trans2 = IntSurf_Out;
3525 }
3526 else {
3527 trans1=trans2=IntSurf_Undecided;
3528 }
3529 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
3530 }
3531 else {
3532 alig = new IntPatch_ALine(curvsol,Standard_False);
3533 }
3534 Standard_Boolean TempFalse1a = Standard_False;
3535 Standard_Boolean TempFalse2a = Standard_False;
3536
3537 //-- ptf et ptl : points debut et fin de alig
3538
3539 ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1a,ptf,first,
3540 TempFalse2a,ptl,last,Multpoint,Tol);
3541 slin.Append(alig);
3542 } //-- boucle sur les lignes
3543 } //-- solution avec au moins une lihne
3544 }
3545 break;
3546
3547 default:
3548 {
3549 return Standard_False;
3550 }
3551 }
3552 return Standard_True;
3553}
3554//=======================================================================
3555//function : IntCyCo
3556//purpose :
3557//=======================================================================
3558Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
3559 const IntSurf_Quadric& Quad2,
3560 const Standard_Real Tol,
3561 const Standard_Boolean Reversed,
3562 Standard_Boolean& Empty,
3563 Standard_Boolean& Multpoint,
3564 IntPatch_SequenceOfLine& slin,
3565 IntPatch_SequenceOfPoint& spnt)
3566
3567{
3568 IntPatch_Point ptsol;
3569
3570 Standard_Integer i;
3571
3572 IntSurf_TypeTrans trans1,trans2;
3573 IntAna_ResultType typint;
3574 gp_Circ cirsol;
3575
3576 gp_Cylinder Cy;
3577 gp_Cone Co;
3578
3579 if (!Reversed) {
3580 Cy = Quad1.Cylinder();
3581 Co = Quad2.Cone();
3582 }
3583 else {
3584 Cy = Quad2.Cylinder();
3585 Co = Quad1.Cone();
3586 }
3587 IntAna_QuadQuadGeo inter(Cy,Co,Tol);
3588
3589 if (!inter.IsDone()) {return Standard_False;}
3590
3591 typint = inter.TypeInter();
3592 Standard_Integer NbSol = inter.NbSolutions();
3593 Empty = Standard_False;
3594
3595 switch (typint) {
3596
3597 case IntAna_Empty : {
3598 Empty = Standard_True;
3599 }
3600 break;
3601
3602 case IntAna_Point :{
3603 gp_Pnt psol(inter.Point(1));
3604 Standard_Real U1,V1,U2,V2;
3605 Quad1.Parameters(psol,U1,V1);
3606 Quad1.Parameters(psol,U2,V2);
3607 ptsol.SetValue(psol,Tol,Standard_True);
3608 ptsol.SetParameters(U1,V1,U2,V2);
3609 spnt.Append(ptsol);
3610 }
3611 break;
3612
3613 case IntAna_Circle: {
3614 gp_Vec Tgt;
3615 gp_Pnt ptref;
3616 Standard_Integer j;
3617 Standard_Real qwe;
3618 //
3619 for(j=1; j<=2; ++j) {
3620 cirsol = inter.Circle(j);
3621 ElCLib::D1(0.,cirsol,ptref,Tgt);
3622 qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
3623 if(qwe> 0.00000001) {
3624 trans1 = IntSurf_Out;
3625 trans2 = IntSurf_In;
3626 }
3627 else if(qwe<-0.00000001) {
3628 trans1 = IntSurf_In;
3629 trans2 = IntSurf_Out;
3630 }
3631 else {
3632 trans1=trans2=IntSurf_Undecided;
3633 }
3634 Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3635 slin.Append(glig);
3636 }
3637 }
3638 break;
3639
3640 case IntAna_NoGeometricSolution: {
3641 gp_Pnt psol;
3642 Standard_Real U1,V1,U2,V2;
3643 IntAna_IntQuadQuad anaint(Cy,Co,Tol);
3644 if (!anaint.IsDone()) {
3645 return Standard_False;
3646 }
3647
3648 if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) {
3649 Empty = Standard_True;
3650 }
3651 else {
3652 NbSol = anaint.NbPnt();
3653 for (i = 1; i <= NbSol; i++) {
3654 psol = anaint.Point(i);
3655 Quad1.Parameters(psol,U1,V1);
3656 Quad2.Parameters(psol,U2,V2);
3657 ptsol.SetValue(psol,Tol,Standard_True);
3658 ptsol.SetParameters(U1,V1,U2,V2);
3659 spnt.Append(ptsol);
3660 }
3661
3662 gp_Pnt ptvalid, ptf, ptl;
3663 gp_Vec tgvalid;
3664 //
3665 Standard_Real first,last,para;
3666 Standard_Boolean tgfound,firstp,lastp,kept;
3667 Standard_Integer kount;
3668 //
3669 //
3670 //IntAna_Curve curvsol;
3671 IntAna_Curve aC;
3672 IntAna_ListOfCurve aLC;
3673 IntAna_ListIteratorOfListOfCurve aIt;
7fd59977 3674
3675 //
3676 NbSol = anaint.NbCurve();
3677 for (i = 1; i <= NbSol; ++i) {
3678 kept = Standard_False;
3679 //curvsol = anaint.Curve(i);
3680 aC=anaint.Curve(i);
3681 aLC.Clear();
96a95605 3682 ExploreCurve(Cy, Co, aC, 10.*Tol, aLC);
7fd59977 3683 //
3684 aIt.Initialize(aLC);
3685 for (; aIt.More(); aIt.Next()) {
3686 IntAna_Curve& curvsol=aIt.Value();
3687 //
3688 curvsol.Domain(first, last);
3689 firstp = !curvsol.IsFirstOpen();
3690 lastp = !curvsol.IsLastOpen();
3691 if (firstp) {
3692 ptf = curvsol.Value(first);
3693 }
3694 if (lastp) {
3695 ptl = curvsol.Value(last);
3696 }
3697 para = last;
3698 kount = 1;
3699 tgfound = Standard_False;
3700
3701 while (!tgfound) {
3702 para = (1.123*first + para)/2.123;
3703 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
3704 if (!tgfound) {
3705 kount ++;
3706 tgfound = kount > 5;
3707 }
3708 }
3709 Handle(IntPatch_ALine) alig;
3710 if (kount <= 5) {
3711 Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
3712 Quad1.Normale(ptvalid));
3713 if(qwe> 0.00000001) {
3714 trans1 = IntSurf_Out;
3715 trans2 = IntSurf_In;
3716 }
3717 else if(qwe<-0.00000001) {
3718 trans1 = IntSurf_In;
3719 trans2 = IntSurf_Out;
3720 }
3721 else {
3722 trans1=trans2=IntSurf_Undecided;
3723 }
3724 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
3725 kept = Standard_True;
3726 }
3727 else {
3728 ptvalid = curvsol.Value(para);
3729 alig = new IntPatch_ALine(curvsol,Standard_False);
3730 kept = Standard_True;
d30895f5 3731 //-- std::cout << "Transition indeterminee" << std::endl;
7fd59977 3732 }
3733 if (kept) {
3734 Standard_Boolean Nfirstp = !firstp;
3735 Standard_Boolean Nlastp = !lastp;
3736 ProcessBounds(alig,slin,Quad1,Quad2,Nfirstp,ptf,first,
3737 Nlastp,ptl,last,Multpoint,Tol);
3738 slin.Append(alig);
3739 }
3740 } // for (; aIt.More(); aIt.Next())
3741 } // for (i = 1; i <= NbSol; ++i)
3742 }
3743 }
3744 break;
3745
3746 default:
3747 return Standard_False;
3748
3749 } // switch (typint)
3750
3751 return Standard_True;
3752}
3753//=======================================================================
3754//function : ExploreCurve
3755//purpose :
3756//=======================================================================
3757Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy,
3758 const gp_Cone& aCo,
3759 IntAna_Curve& aC,
3760 const Standard_Real aTol,
3761 IntAna_ListOfCurve& aLC)
3762
3763{
3764 Standard_Boolean bFind=Standard_False;
3765 Standard_Real aTheta, aT1, aT2, aDst;
3766 gp_Pnt aPapx, aPx;
3767 //
3768 //aC.Dump();
3769 //
3770 aLC.Clear();
3771 aLC.Append(aC);
3772 //
3773 aPapx=aCo.Apex();
3774 //
3775 aC.Domain(aT1, aT2);
3776 //
3777 aPx=aC.Value(aT1);
3778 aDst=aPx.Distance(aPapx);
3779 if (aDst<aTol) {
3780 return bFind;
3781 }
3782 aPx=aC.Value(aT2);
3783 aDst=aPx.Distance(aPapx);
3784 if (aDst<aTol) {
3785 return bFind;
3786 }
3787 //
3788 bFind=aC.FindParameter(aPapx, aTheta);
3789 if (!bFind){
3790 return bFind;
3791 }
3792 //
3793 aPx=aC.Value(aTheta);
3794 aDst=aPx.Distance(aPapx);
3795 if (aDst>aTol) {
3796 return !bFind;
3797 }
3798 //
3799 // need to be splitted at aTheta
3800 IntAna_Curve aC1, aC2;
3801 //
3802 aC1=aC;
3803 aC1.SetDomain(aT1, aTheta);
3804 aC2=aC;
3805 aC2.SetDomain(aTheta, aT2);
3806 //
3807 aLC.Clear();
3808 aLC.Append(aC1);
3809 aLC.Append(aC2);
3810 //
3811 return bFind;
3812}
Open CASCADE Technology repository