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. |
26 | static const Standard_Real aNulValue = 1.0e-11; |
27 | |
d30895f5 |
28 | static void ShortCosForm( const Standard_Real theCosFactor, |
29 | const Standard_Real theSinFactor, |
30 | Standard_Real& theCoeff, |
31 | Standard_Real& theAngle); |
fa0291ff |
32 | // |
7fd59977 |
33 | static |
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 |
40 | static 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 | |
48 | class 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 |
88 | public: |
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 | |
162 | ComputationMethods::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 | |
354 | class WorkWithBoundaries |
355 | { |
356 | public: |
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 | |
436 | protected: |
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 | |
452 | private: |
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 | |
467 | static |
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 |
474 | static 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 |
491 | static 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 | //======================================================================= |
506 | static 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 | //======================================================================= |
531 | static 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 | //======================================================================= |
586 | static 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 | //======================================================================= |
620 | static 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 | //======================================================================= |
756 | void 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 |
897 | Standard_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 | //======================================================================= |
1298 | static 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 |
1352 | Standard_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 |
1424 | Standard_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 |
1497 | Standard_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 |
1517 | Standard_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 |
1537 | Standard_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 |
1678 | Standard_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 | //======================================================================= |
1733 | static 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 |
1777 | static 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 | //======================================================================= |
1815 | static 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 |
1960 | void 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 | //======================================================================= |
2092 | static 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 |
2221 | Standard_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 |
2412 | static 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 | //======================================================================= |
2538 | void 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 |
2601 | IntPatch_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 | |
4dba155d |
2818 | Standard_Real aCriticalDelta[aNbCritPointsMax]; |
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 | //======================================================================= |
3518 | Standard_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 | //======================================================================= |
3749 | Standard_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 | //======================================================================= |
3948 | Standard_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 | } |