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