projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0028389: Data Exchange - Import of STEP Saved Views and Clipping Planes
[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
d30895f5 501//=======================================================================
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
b5ef9d91 524//=======================================================================
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
b5ef9d91 1348//=======================================================================
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
ecc4f148 1533//=======================================================================
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
79997052 2216//=======================================================================
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
ecc4f148 2407//=======================================================================
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);
3215 else
3216 aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS1(aU1c, aV1c);
3217
3218 anUmaxAdded = Max(anUmaxAdded, aU1c);
3219 isChanged = Standard_True;
3220 }
3221
3222 if(!isChanged)
3223 { //If anUmaxAdded were not changed in previous cycle then
3224 //we would break existing WLines.
3225 break;
3226 }
b5ef9d91 3227 }
3228
5b9e1842 3229 for(Standard_Integer i = 0; i < aNbWLines; i++)
3230 {
3231 if(isAddingWLEnabled[i])
3232 {
3233 continue;
3234 }
b5ef9d91 3235
d30895f5 3236 ComputationMethods::CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs,
3237 aU2[i], aV1[i], aV2[i]);
b5ef9d91 3238
5b9e1842 3239 AddPointIntoWL( theQuad1, theQuad2, anEquationCoeffs, isTheReverse,
3240 Standard_True, gp_Pnt2d(anUmaxAdded, aV1[i]),
3241 gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
3242 aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
3243 aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
3244 i, theTol3D, theTol2D, Standard_False);
3245 }
b5ef9d91 3246 }
3247
ecc4f148 3248 break;
3249 }
3250
7c32c7c4 3251 //Step computing
ecc4f148 3252
7c32c7c4 3253 {
4dba155d 3254 //Step of aU1-parameter is computed adaptively. The algorithm
3255 //aims to provide given aDeltaV1 and aDeltaV2 values (if it is
3256 //possible because the intersection line can go along V-isoline)
3257 //in every iteration. It allows avoiding "flying" intersection
3258 //points too far each from other (see issue #24915).
3259
d30895f5 3260 const Standard_Real aDeltaV1 = aRangeS1.Delta()/IntToReal(aNbPoints),
3261 aDeltaV2 = aRangeS2.Delta()/IntToReal(aNbPoints);
b5ef9d91 3262
3263 math_Matrix aMatr(1, 3, 1, 5);
ecc4f148 3264
b5ef9d91 3265 Standard_Real aMinUexp = RealLast();
3266
3267 for(Standard_Integer i = 0; i < aNbWLines; i++)
3268 {
3269 if(theTol2D < (anUexpect[i] - anU1))
3270 {
3271 continue;
3272 }
ecc4f148 3273
b5ef9d91 3274 if(aWLFindStatus[i] == WLFStatus_Absent)
3275 {
3276 anUexpect[i] += aStepMax;
3277 aMinUexp = Min(aMinUexp, anUexpect[i]);
3278 continue;
3279 }
7c32c7c4 3280
b5ef9d91 3281 Standard_Real aStepTmp = aStepMax;
ecc4f148 3282
b5ef9d91 3283 const Standard_Real aSinU1 = sin(anU1),
3284 aCosU1 = cos(anU1),
3285 aSinU2 = sin(aU2[i]),
3286 aCosU2 = cos(aU2[i]);
ecc4f148 3287
b5ef9d91 3288 aMatr.SetCol(1, anEquationCoeffs.mVecC1);
3289 aMatr.SetCol(2, anEquationCoeffs.mVecC2);
3290 aMatr.SetCol(3, anEquationCoeffs.mVecA1*aSinU1 - anEquationCoeffs.mVecB1*aCosU1);
3291 aMatr.SetCol(4, anEquationCoeffs.mVecA2*aSinU2 - anEquationCoeffs.mVecB2*aCosU2);
3292 aMatr.SetCol(5, anEquationCoeffs.mVecA1*aCosU1 + anEquationCoeffs.mVecB1*aSinU1 +
3293 anEquationCoeffs.mVecA2*aCosU2 + anEquationCoeffs.mVecB2*aSinU2 +
3294 anEquationCoeffs.mVecD);
ecc4f148 3295
4dba155d 3296 //The main idea is in solving of linearized system (2)
3297 //(see description to ComputationMethods class) in order to find new U1-value
3298 //to provide new value V1 or V2, which differs from current one by aDeltaV1 or
3299 //aDeltaV2 respectively.
3300
3301 //While linearizing, following Taylor formulas are used:
3302 // cos(x0+dx) = cos(x0) - sin(x0)*dx
3303 // sin(x0+dx) = sin(x0) + cos(x0)*dx
3304
3305 //Consequently cos(U1), cos(U2), sin(U1) and sin(U2) in the system (2)
3306 //must be substituted by corresponding values.
3307
3308 //ATTENTION!!!
3309 //The solution is approximate. More over, all requirements to the
3310 //linearization must be satisfied in order to obtain quality result.
3311
7365fad6 3312 if(!StepComputing(aMatr, aV1[i], aV2[i], aDeltaV1, aDeltaV2, aStepTmp))
b5ef9d91 3313 {
3314 //To avoid cycling-up
3315 anUexpect[i] += aStepMax;
3316 aMinUexp = Min(aMinUexp, anUexpect[i]);
ecc4f148 3317
b5ef9d91 3318 continue;
3319 }
3320
3321 if(aStepTmp < aStepMin)
3322 aStepTmp = aStepMin;
ecc4f148 3323
b5ef9d91 3324 if(aStepTmp > aStepMax)
3325 aStepTmp = aStepMax;
3326
3327 anUexpect[i] = anU1 + aStepTmp;
3328 aMinUexp = Min(aMinUexp, anUexpect[i]);
3329 }
ecc4f148 3330
b5ef9d91 3331 anU1 = aMinUexp;
3332 }
ecc4f148 3333
b5ef9d91 3334 if(Precision::PConfusion() >= (anUl - anU1))
ecc4f148 3335 anU1 = anUl;
3336
3337 anUf = anU1;
3338
7c32c7c4 3339 for(Standard_Integer i = 0; i < aNbWLines; i++)
baf72cd2 3340 {
7c32c7c4 3341 if(aWLine[i]->NbPnts() != 1)
3342 isAddedIntoWL[i] = Standard_False;
b5ef9d91 3343
3344 if(anU1 == anUl)
3345 {//strictly equal. Tolerance is considered above.
3346 anUexpect[i] = anUl;
3347 }
baf72cd2 3348 }
ecc4f148 3349 }
ecc4f148 3350
7c32c7c4 3351 for(Standard_Integer i = 0; i < aNbWLines; i++)
ecc4f148 3352 {
7c32c7c4 3353 if((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
baf72cd2 3354 {
7c32c7c4 3355 isTheEmpty = Standard_False;
3356 Standard_Real u1, v1, u2, v2;
3357 aWLine[i]->Point(1).Parameters(u1, v1, u2, v2);
3358 IntPatch_Point aP;
3359 aP.SetParameter(u1);
3360 aP.SetParameters(u1, v1, u2, v2);
3361 aP.SetTolerance(theTol3D);
3362 aP.SetValue(aWLine[i]->Point(1).Value());
3363
3364 theSPnt.Append(aP);
baf72cd2 3365 }
7c32c7c4 3366 else if(aWLine[i]->NbPnts() > 1)
baf72cd2 3367 {
7c32c7c4 3368 Standard_Boolean isGood = Standard_True;
ecc4f148 3369
7c32c7c4 3370 if(aWLine[i]->NbPnts() == 2)
3371 {
3372 const IntSurf_PntOn2S& aPf = aWLine[i]->Point(1);
3373 const IntSurf_PntOn2S& aPl = aWLine[i]->Point(2);
3374
3375 if(aPf.IsSame(aPl, Precision::Confusion()))
3376 isGood = Standard_False;
3377 }
ecc4f148 3378
7c32c7c4 3379 if(isGood)
3380 {
3381 isTheEmpty = Standard_False;
3382 isAddedIntoWL[i] = Standard_True;
3383 SeekAdditionalPoints( theQuad1, theQuad2, aWLine[i]->Curve(),
b5ef9d91 3384 anEquationCoeffs, i, aNbPoints, 1,
d30895f5 3385 aWLine[i]->NbPnts(), theTol2D, aPeriod,
3386 isTheReverse);
7c32c7c4 3387
3388 aWLine[i]->ComputeVertexParameters(theTol3D);
3389 theSlin.Append(aWLine[i]);
3390 }
3391 }
3392 else
3393 {
3394 isAddedIntoWL[i] = Standard_False;
baf72cd2 3395 }
b5ef9d91 3396
d30895f5 3397#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
b5ef9d91 3398 //aWLine[i]->Dump();
3399#endif
3400 }
3401 }
3402 }
3403
3404
7365fad6 3405 //Delete the points in theSPnt, which
3406 //lie at least in one of the line in theSlin.
b5ef9d91 3407 for(Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
3408 {
3409 for(Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); aNbLin++)
3410 {
d30895f5 3411 Handle(IntPatch_WLine) aWLine1 (Handle(IntPatch_WLine)::
3412 DownCast(theSlin.Value(aNbLin)));
b5ef9d91 3413
3414 const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1);
3415 const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aWLine1->NbPnts());
3416
3417 const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNbPnt).PntOn2S();
3418 if( aPntCur.IsSame(aPntFWL1, Precision::Confusion()) ||
3419 aPntCur.IsSame(aPntLWL1, Precision::Confusion()))
3420 {
3421 theSPnt.Remove(aNbPnt);
3422 aNbPnt--;
3423 break;
ecc4f148 3424 }
3425 }
3426 }
3427
b5ef9d91 3428 const Standard_Real aDU = aStepMin + Epsilon(aStepMin);
3429 //Try to add new points in the neighbourhood of existing point
3430 for(Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
3431 {
3432 const IntPatch_Point& aPnt2S = theSPnt.Value(aNbPnt);
3433
3434 Standard_Real u1, v1, u2, v2;
3435 aPnt2S.Parameters(u1, v1, u2, v2);
3436
3437 Handle(IntSurf_LineOn2S) aL2S = new IntSurf_LineOn2S();
3438 Handle(IntPatch_WLine) aWLine = new IntPatch_WLine(aL2S, Standard_False);
3439 aWLine->Curve()->Add(aPnt2S.PntOn2S());
3440
3441 //Define the index of WLine, which lies the point aPnt2S in.
3442 Standard_Real anUf = 0.0, anUl = 0.0, aCurU2 = 0.0;
3443 Standard_Integer anIndex = 0;
3444 if(isTheReverse)
3445 {
3446 anUf = Max(u2 - aStepMax, aUSurf1f);
3447 anUl = u2;
3448 aCurU2 = u1;
3449 }
3450 else
3451 {
3452 anUf = Max(u1 - aStepMax, aUSurf1f);
3453 anUl = u1;
3454 aCurU2 = u2;
3455 }
3456 Standard_Real aDelta = RealLast();
3457 for (Standard_Integer i = 0; i < aNbWLines; i++)
3458 {
3459 Standard_Real anU2t = 0.0;
d30895f5 3460 if(!ComputationMethods::CylCylComputeParameters(anUl, i, anEquationCoeffs, anU2t))
b5ef9d91 3461 continue;
3462
51740958 3463 const Standard_Real aDU2 = Abs(anU2t - aCurU2);
3464 if(aDU2 < aDelta)
b5ef9d91 3465 {
51740958 3466 aDelta = aDU2;
b5ef9d91 3467 anIndex = i;
3468 }
3469 }
3470
3471 //Try to fill aWLine by additional points
3472 while(anUl - anUf > RealSmall())
3473 {
3474 Standard_Real anU2 = 0.0, anV1 = 0.0, anV2 = 0.0;
3475 Standard_Boolean isDone =
d30895f5 3476 ComputationMethods::CylCylComputeParameters(anUf, anIndex, anEquationCoeffs,
b5ef9d91 3477 anU2, anV1, anV2);
3478 anUf += aDU;
5b9e1842 3479
b5ef9d91 3480 if(!isDone)
3481 {
3482 continue;
3483 }
3484
d30895f5 3485 if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse,
3486 Standard_True, gp_Pnt2d(anUf, anV1), gp_Pnt2d(anU2, anV2),
b5ef9d91 3487 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
3488 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l,
3489 aPeriod, aWLine->Curve(), anIndex, theTol3D,
3490 theTol2D, Standard_False, Standard_True))
3491 {
3492 break;
3493 }
3494 }
3495
3496 if(aWLine->NbPnts() > 1)
3497 {
3498 SeekAdditionalPoints( theQuad1, theQuad2, aWLine->Curve(),
3499 anEquationCoeffs, anIndex, aNbMinPoints,
d30895f5 3500 1, aWLine->NbPnts(), theTol2D, aPeriod,
3501 isTheReverse);
b5ef9d91 3502
3503 aWLine->ComputeVertexParameters(theTol3D);
3504 theSlin.Append(aWLine);
5b9e1842 3505
b5ef9d91 3506 theSPnt.Remove(aNbPnt);
3507 aNbPnt--;
3508 }
3509 }
3510
e146ef9a 3511 return IntPatch_ImpImpIntersection::IntStatus_OK;
ecc4f148 3512}
7fd59977 3513
3514//=======================================================================
3515//function : IntCySp
3516//purpose :
3517//=======================================================================
3518Standard_Boolean IntCySp(const IntSurf_Quadric& Quad1,
3519 const IntSurf_Quadric& Quad2,
3520 const Standard_Real Tol,
3521 const Standard_Boolean Reversed,
3522 Standard_Boolean& Empty,
3523 Standard_Boolean& Multpoint,
3524 IntPatch_SequenceOfLine& slin,
3525 IntPatch_SequenceOfPoint& spnt)
3526
3527{
3528 Standard_Integer i;
3529
3530 IntSurf_TypeTrans trans1,trans2;
3531 IntAna_ResultType typint;
3532 IntPatch_Point ptsol;
3533 gp_Circ cirsol;
3534
3535 gp_Cylinder Cy;
3536 gp_Sphere Sp;
3537
3538 if (!Reversed) {
3539 Cy = Quad1.Cylinder();
3540 Sp = Quad2.Sphere();
3541 }
3542 else {
3543 Cy = Quad2.Cylinder();
3544 Sp = Quad1.Sphere();
3545 }
3546 IntAna_QuadQuadGeo inter(Cy,Sp,Tol);
3547
3548 if (!inter.IsDone()) {return Standard_False;}
3549
3550 typint = inter.TypeInter();
3551 Standard_Integer NbSol = inter.NbSolutions();
3552 Empty = Standard_False;
3553
3554 switch (typint) {
3555
3556 case IntAna_Empty :
3557 {
3558 Empty = Standard_True;
3559 }
3560 break;
3561
3562 case IntAna_Point :
3563 {
3564 gp_Pnt psol(inter.Point(1));
3565 Standard_Real U1,V1,U2,V2;
3566 Quad1.Parameters(psol,U1,V1);
3567 Quad2.Parameters(psol,U2,V2);
3568 ptsol.SetValue(psol,Tol,Standard_True);
3569 ptsol.SetParameters(U1,V1,U2,V2);
3570 spnt.Append(ptsol);
3571 }
3572 break;
3573
3574 case IntAna_Circle:
3575 {
3576 cirsol = inter.Circle(1);
3577 gp_Vec Tgt;
3578 gp_Pnt ptref;
3579 ElCLib::D1(0.,cirsol,ptref,Tgt);
3580
3581 if (NbSol == 1) {
3582 gp_Vec TestCurvature(ptref,Sp.Location());
3583 gp_Vec Normsp,Normcyl;
3584 if (!Reversed) {
3585 Normcyl = Quad1.Normale(ptref);
3586 Normsp = Quad2.Normale(ptref);
3587 }
3588 else {
3589 Normcyl = Quad2.Normale(ptref);
3590 Normsp = Quad1.Normale(ptref);
3591 }
3592
3593 IntSurf_Situation situcyl;
3594 IntSurf_Situation situsp;
3595
3596 if (Normcyl.Dot(TestCurvature) > 0.) {
3597 situsp = IntSurf_Outside;
3598 if (Normsp.Dot(Normcyl) > 0.) {
3599 situcyl = IntSurf_Inside;
3600 }
3601 else {
3602 situcyl = IntSurf_Outside;
3603 }
3604 }
3605 else {
3606 situsp = IntSurf_Inside;
3607 if (Normsp.Dot(Normcyl) > 0.) {
3608 situcyl = IntSurf_Outside;
3609 }
3610 else {
3611 situcyl = IntSurf_Inside;
3612 }
3613 }
3614 Handle(IntPatch_GLine) glig;
3615 if (!Reversed) {
3616 glig = new IntPatch_GLine(cirsol, Standard_True, situcyl, situsp);
3617 }
3618 else {
3619 glig = new IntPatch_GLine(cirsol, Standard_True, situsp, situcyl);
3620 }
3621 slin.Append(glig);
3622 }
3623 else {
3624 if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
3625 trans1 = IntSurf_Out;
3626 trans2 = IntSurf_In;
3627 }
3628 else {
3629 trans1 = IntSurf_In;
3630 trans2 = IntSurf_Out;
3631 }
3632 Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3633 slin.Append(glig);
3634
3635 cirsol = inter.Circle(2);
3636 ElCLib::D1(0.,cirsol,ptref,Tgt);
3637 Standard_Real qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
3638 if(qwe> 0.0000001) {
3639 trans1 = IntSurf_Out;
3640 trans2 = IntSurf_In;
3641 }
3642 else if(qwe<-0.0000001) {
3643 trans1 = IntSurf_In;
3644 trans2 = IntSurf_Out;
3645 }
3646 else {
3647 trans1=trans2=IntSurf_Undecided;
3648 }
3649 glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3650 slin.Append(glig);
3651 }
3652 }
3653 break;
3654
3655 case IntAna_NoGeometricSolution:
3656 {
3657 gp_Pnt psol;
3658 Standard_Real U1,V1,U2,V2;
3659 IntAna_IntQuadQuad anaint(Cy,Sp,Tol);
3660 if (!anaint.IsDone()) {
3661 return Standard_False;
3662 }
3663
3664 if (anaint.NbPnt()==0 && anaint.NbCurve()==0) {
3665 Empty = Standard_True;
3666 }
3667 else {
3668
3669 NbSol = anaint.NbPnt();
3670 for (i = 1; i <= NbSol; i++) {
3671 psol = anaint.Point(i);
3672 Quad1.Parameters(psol,U1,V1);
3673 Quad2.Parameters(psol,U2,V2);
3674 ptsol.SetValue(psol,Tol,Standard_True);
3675 ptsol.SetParameters(U1,V1,U2,V2);
3676 spnt.Append(ptsol);
3677 }
3678
3679 gp_Pnt ptvalid,ptf,ptl;
3680 gp_Vec tgvalid;
3681 Standard_Real first,last,para;
3682 IntAna_Curve curvsol;
3683 Standard_Boolean tgfound;
3684 Standard_Integer kount;
3685
3686 NbSol = anaint.NbCurve();
3687 for (i = 1; i <= NbSol; i++) {
3688 curvsol = anaint.Curve(i);
3689 curvsol.Domain(first,last);
3690 ptf = curvsol.Value(first);
3691 ptl = curvsol.Value(last);
3692
3693 para = last;
3694 kount = 1;
3695 tgfound = Standard_False;
3696
3697 while (!tgfound) {
3698 para = (1.123*first + para)/2.123;
3699 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
3700 if (!tgfound) {
3701 kount ++;
3702 tgfound = kount > 5;
3703 }
3704 }
3705 Handle(IntPatch_ALine) alig;
3706 if (kount <= 5) {
3707 Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
3708 Quad1.Normale(ptvalid));
3709 if(qwe> 0.00000001) {
3710 trans1 = IntSurf_Out;
3711 trans2 = IntSurf_In;
3712 }
3713 else if(qwe<-0.00000001) {
3714 trans1 = IntSurf_In;
3715 trans2 = IntSurf_Out;
3716 }
3717 else {
3718 trans1=trans2=IntSurf_Undecided;
3719 }
3720 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
3721 }
3722 else {
3723 alig = new IntPatch_ALine(curvsol,Standard_False);
3724 }
3725 Standard_Boolean TempFalse1a = Standard_False;
3726 Standard_Boolean TempFalse2a = Standard_False;
3727
3728 //-- ptf et ptl : points debut et fin de alig
3729
3730 ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1a,ptf,first,
3731 TempFalse2a,ptl,last,Multpoint,Tol);
3732 slin.Append(alig);
3733 } //-- boucle sur les lignes
3734 } //-- solution avec au moins une lihne
3735 }
3736 break;
3737
3738 default:
3739 {
3740 return Standard_False;
3741 }
3742 }
3743 return Standard_True;
3744}
3745//=======================================================================
3746//function : IntCyCo
3747//purpose :
3748//=======================================================================
3749Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
3750 const IntSurf_Quadric& Quad2,
3751 const Standard_Real Tol,
3752 const Standard_Boolean Reversed,
3753 Standard_Boolean& Empty,
3754 Standard_Boolean& Multpoint,
3755 IntPatch_SequenceOfLine& slin,
3756 IntPatch_SequenceOfPoint& spnt)
3757
3758{
3759 IntPatch_Point ptsol;
3760
3761 Standard_Integer i;
3762
3763 IntSurf_TypeTrans trans1,trans2;
3764 IntAna_ResultType typint;
3765 gp_Circ cirsol;
3766
3767 gp_Cylinder Cy;
3768 gp_Cone Co;
3769
3770 if (!Reversed) {
3771 Cy = Quad1.Cylinder();
3772 Co = Quad2.Cone();
3773 }
3774 else {
3775 Cy = Quad2.Cylinder();
3776 Co = Quad1.Cone();
3777 }
3778 IntAna_QuadQuadGeo inter(Cy,Co,Tol);
3779
3780 if (!inter.IsDone()) {return Standard_False;}
3781
3782 typint = inter.TypeInter();
3783 Standard_Integer NbSol = inter.NbSolutions();
3784 Empty = Standard_False;
3785
3786 switch (typint) {
3787
3788 case IntAna_Empty : {
3789 Empty = Standard_True;
3790 }
3791 break;
3792
3793 case IntAna_Point :{
3794 gp_Pnt psol(inter.Point(1));
3795 Standard_Real U1,V1,U2,V2;
3796 Quad1.Parameters(psol,U1,V1);
3797 Quad1.Parameters(psol,U2,V2);
3798 ptsol.SetValue(psol,Tol,Standard_True);
3799 ptsol.SetParameters(U1,V1,U2,V2);
3800 spnt.Append(ptsol);
3801 }
3802 break;
3803
3804 case IntAna_Circle: {
3805 gp_Vec Tgt;
3806 gp_Pnt ptref;
3807 Standard_Integer j;
3808 Standard_Real qwe;
3809 //
3810 for(j=1; j<=2; ++j) {
3811 cirsol = inter.Circle(j);
3812 ElCLib::D1(0.,cirsol,ptref,Tgt);
3813 qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
3814 if(qwe> 0.00000001) {
3815 trans1 = IntSurf_Out;
3816 trans2 = IntSurf_In;
3817 }
3818 else if(qwe<-0.00000001) {
3819 trans1 = IntSurf_In;
3820 trans2 = IntSurf_Out;
3821 }
3822 else {
3823 trans1=trans2=IntSurf_Undecided;
3824 }
3825 Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3826 slin.Append(glig);
3827 }
3828 }
3829 break;
3830
3831 case IntAna_NoGeometricSolution: {
3832 gp_Pnt psol;
3833 Standard_Real U1,V1,U2,V2;
3834 IntAna_IntQuadQuad anaint(Cy,Co,Tol);
3835 if (!anaint.IsDone()) {
3836 return Standard_False;
3837 }
3838
3839 if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) {
3840 Empty = Standard_True;
3841 }
3842 else {
3843 NbSol = anaint.NbPnt();
3844 for (i = 1; i <= NbSol; i++) {
3845 psol = anaint.Point(i);
3846 Quad1.Parameters(psol,U1,V1);
3847 Quad2.Parameters(psol,U2,V2);
3848 ptsol.SetValue(psol,Tol,Standard_True);
3849 ptsol.SetParameters(U1,V1,U2,V2);
3850 spnt.Append(ptsol);
3851 }
3852
3853 gp_Pnt ptvalid, ptf, ptl;
3854 gp_Vec tgvalid;
3855 //
3856 Standard_Real first,last,para;
3857 Standard_Boolean tgfound,firstp,lastp,kept;
3858 Standard_Integer kount;
3859 //
3860 //
3861 //IntAna_Curve curvsol;
3862 IntAna_Curve aC;
3863 IntAna_ListOfCurve aLC;
3864 IntAna_ListIteratorOfListOfCurve aIt;
7fd59977 3865
3866 //
3867 NbSol = anaint.NbCurve();
3868 for (i = 1; i <= NbSol; ++i) {
3869 kept = Standard_False;
3870 //curvsol = anaint.Curve(i);
3871 aC=anaint.Curve(i);
3872 aLC.Clear();
96a95605 3873 ExploreCurve(Cy, Co, aC, 10.*Tol, aLC);
7fd59977 3874 //
3875 aIt.Initialize(aLC);
3876 for (; aIt.More(); aIt.Next()) {
3877 IntAna_Curve& curvsol=aIt.Value();
3878 //
3879 curvsol.Domain(first, last);
3880 firstp = !curvsol.IsFirstOpen();
3881 lastp = !curvsol.IsLastOpen();
3882 if (firstp) {
3883 ptf = curvsol.Value(first);
3884 }
3885 if (lastp) {
3886 ptl = curvsol.Value(last);
3887 }
3888 para = last;
3889 kount = 1;
3890 tgfound = Standard_False;
3891
3892 while (!tgfound) {
3893 para = (1.123*first + para)/2.123;
3894 tgfound = curvsol.D1u(para,ptvalid,tgvalid);
3895 if (!tgfound) {
3896 kount ++;
3897 tgfound = kount > 5;
3898 }
3899 }
3900 Handle(IntPatch_ALine) alig;
3901 if (kount <= 5) {
3902 Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
3903 Quad1.Normale(ptvalid));
3904 if(qwe> 0.00000001) {
3905 trans1 = IntSurf_Out;
3906 trans2 = IntSurf_In;
3907 }
3908 else if(qwe<-0.00000001) {
3909 trans1 = IntSurf_In;
3910 trans2 = IntSurf_Out;
3911 }
3912 else {
3913 trans1=trans2=IntSurf_Undecided;
3914 }
3915 alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
3916 kept = Standard_True;
3917 }
3918 else {
3919 ptvalid = curvsol.Value(para);
3920 alig = new IntPatch_ALine(curvsol,Standard_False);
3921 kept = Standard_True;
d30895f5 3922 //-- std::cout << "Transition indeterminee" << std::endl;
7fd59977 3923 }
3924 if (kept) {
3925 Standard_Boolean Nfirstp = !firstp;
3926 Standard_Boolean Nlastp = !lastp;
3927 ProcessBounds(alig,slin,Quad1,Quad2,Nfirstp,ptf,first,
3928 Nlastp,ptl,last,Multpoint,Tol);
3929 slin.Append(alig);
3930 }
3931 } // for (; aIt.More(); aIt.Next())
3932 } // for (i = 1; i <= NbSol; ++i)
3933 }
3934 }
3935 break;
3936
3937 default:
3938 return Standard_False;
3939
3940 } // switch (typint)
3941
3942 return Standard_True;
3943}
3944//=======================================================================
3945//function : ExploreCurve
3946//purpose :
3947//=======================================================================
3948Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy,
3949 const gp_Cone& aCo,
3950 IntAna_Curve& aC,
3951 const Standard_Real aTol,
3952 IntAna_ListOfCurve& aLC)
3953
3954{
3955 Standard_Boolean bFind=Standard_False;
3956 Standard_Real aTheta, aT1, aT2, aDst;
3957 gp_Pnt aPapx, aPx;
3958 //
3959 //aC.Dump();
3960 //
3961 aLC.Clear();
3962 aLC.Append(aC);
3963 //
3964 aPapx=aCo.Apex();
3965 //
3966 aC.Domain(aT1, aT2);
3967 //
3968 aPx=aC.Value(aT1);
3969 aDst=aPx.Distance(aPapx);
3970 if (aDst<aTol) {
3971 return bFind;
3972 }
3973 aPx=aC.Value(aT2);
3974 aDst=aPx.Distance(aPapx);
3975 if (aDst<aTol) {
3976 return bFind;
3977 }
3978 //
3979 bFind=aC.FindParameter(aPapx, aTheta);
3980 if (!bFind){
3981 return bFind;
3982 }
3983 //
3984 aPx=aC.Value(aTheta);
3985 aDst=aPx.Distance(aPapx);
3986 if (aDst>aTol) {
3987 return !bFind;
3988 }
3989 //
3990 // need to be splitted at aTheta
3991 IntAna_Curve aC1, aC2;
3992 //
3993 aC1=aC;
3994 aC1.SetDomain(aT1, aTheta);
3995 aC2=aC;
3996 aC2.SetDomain(aTheta, aT2);
3997 //
3998 aLC.Clear();
3999 aLC.Append(aC1);
4000 aLC.Append(aC2);
4001 //
4002 return bFind;
4003}
Open CASCADE Technology repository