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