projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0025782: The result of intersection between two cylinders is incorrect
[occt.git] / src / IntAna / IntAna_QuadQuadGeo.cxx
CommitLineData
b311480e 1// Created on: 1992-08-06
2// Created by: Laurent BUCHARD
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.
b311480e 16
7fd59977 17//----------------------------------------------------------------------
18//-- Purposse: Geometric Intersection between two Natural Quadric
19//-- If the intersection is not a conic,
20//-- analytical methods must be called.
21//----------------------------------------------------------------------
0797d9d3 22#ifndef OCCT_DEBUG
7fd59977 23#define No_Standard_RangeError
24#define No_Standard_OutOfRange
25#endif
26
27#include <IntAna_QuadQuadGeo.ixx>
28
29#include <IntAna_IntConicQuad.hxx>
30#include <StdFail_NotDone.hxx>
31#include <Standard_DomainError.hxx>
32#include <Standard_OutOfRange.hxx>
33#include <math_DirectPolynomialRoots.hxx>
34
35#include <gp.hxx>
36#include <gp_Pln.hxx>
37#include <gp_Vec.hxx>
38#include <ElSLib.hxx>
39#include <ElCLib.hxx>
40
41#include <gp_Dir.hxx>
42#include <gp_XYZ.hxx>
43#include <gp_Pnt2d.hxx>
44#include <gp_Vec2d.hxx>
45#include <gp_Dir2d.hxx>
46
47
48static
49 gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D);
77088633 50static
51 void RefineDir(gp_Dir& aDir);
7fd59977 52
53//=======================================================================
a34f083b 54//class : AxeOperator
7fd59977 55//purpose : O p e r a t i o n s D i v e r s e s s u r d e s A x 1
56//=======================================================================
57class AxeOperator {
58 public:
59 AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2);
60
61 void Distance(Standard_Real& dist,
7eed5d29 62 Standard_Real& Param1,
63 Standard_Real& Param2);
7fd59977 64
65 gp_Pnt PtIntersect() {
66 return ptintersect;
67 }
68 Standard_Boolean Coplanar(void) {
69 return thecoplanar;
70 }
71 Standard_Boolean Same(void) {
72 return theparallel && (thedistance<myEPSILON_DISTANCE);
73 }
74 Standard_Real Distance(void) {
75 return thedistance ;
76 }
77 Standard_Boolean Intersect(void) {
78 return (thecoplanar && (!theparallel));
79 }
80 Standard_Boolean Parallel(void) {
81 return theparallel;
82 }
83 Standard_Boolean Normal(void) {
84 return thenormal;
85 }
86
87 protected:
88 Standard_Real Det33(const Standard_Real a11,
7eed5d29 89 const Standard_Real a12,
90 const Standard_Real a13,
91 const Standard_Real a21,
92 const Standard_Real a22,
93 const Standard_Real a23,
94 const Standard_Real a31,
95 const Standard_Real a32,
96 const Standard_Real a33) {
7fd59977 97 Standard_Real theReturn =
98 a11*(a22*a33-a32*a23) - a21*(a12*a33-a32*a13) + a31*(a12*a23-a22*a13) ;
99 return theReturn ;
100 }
101
102 private:
103 gp_Pnt ptintersect;
104 gp_Ax1 Axe1;
105 gp_Ax1 Axe2;
106 Standard_Real thedistance;
107 Standard_Boolean theparallel;
108 Standard_Boolean thecoplanar;
109 Standard_Boolean thenormal;
110 //
111 Standard_Real myEPSILON_DISTANCE;
112 Standard_Real myEPSILON_AXES_PARA;
113};
114
115//=======================================================================
116//function : AxeOperator::AxeOperator
117//purpose :
118//=======================================================================
a34f083b 119AxeOperator::AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2)
7fd59977 120{
121 myEPSILON_DISTANCE=0.00000000000001;
122 myEPSILON_AXES_PARA=0.000000000001;
123 Axe1=A1;
124 Axe2=A2;
125 //---------------------------------------------------------------------
126 gp_Dir V1=Axe1.Direction();
127 gp_Dir V2=Axe2.Direction();
128 gp_Pnt P1=Axe1.Location();
129 gp_Pnt P2=Axe2.Location();
77088633 130 //
131 RefineDir(V1);
132 RefineDir(V2);
7fd59977 133 thecoplanar= Standard_False;
134 thenormal = Standard_False;
135
136 //--- check if the two axis are parallel
137 theparallel=V1.IsParallel(V2, myEPSILON_AXES_PARA);
138 //--- Distance between the two axis
139 gp_XYZ perp(A1.Direction().XYZ().Crossed(A2.Direction().XYZ()));
140 if (theparallel) {
141 gp_Lin L1(A1);
142 thedistance = L1.Distance(A2.Location());
143 }
144 else {
145 thedistance = Abs(gp_Vec(perp.Normalized()).Dot(gp_Vec(Axe1.Location(),
7eed5d29 146 Axe2.Location())));
7fd59977 147 }
148 //--- check if Axis are Coplanar
149 Standard_Real D33;
150 if(thedistance<myEPSILON_DISTANCE) {
151 D33=Det33(V1.X(),V1.Y(),V1.Z()
7eed5d29 152 ,V2.X(),V2.Y(),V2.Z()
153 ,P1.X()-P2.X(),P1.Y()-P2.Y(),P1.Z()-P2.Z());
7fd59977 154 if(Abs(D33)<=myEPSILON_DISTANCE) {
155 thecoplanar=Standard_True;
156 }
157 }
158 else {
159 thecoplanar=Standard_True;
160 thenormal=(V1.Dot(V2)==0.0)? Standard_True : Standard_False;
161 }
162 //--- check if the two axis are concurrent
163 if(thecoplanar && (!theparallel)) {
164 Standard_Real smx=P2.X() - P1.X();
165 Standard_Real smy=P2.Y() - P1.Y();
166 Standard_Real smz=P2.Z() - P1.Z();
167 Standard_Real Det1,Det2,Det3,A;
168 Det1=V1.Y() * V2.X() - V1.X() * V2.Y();
169 Det2=V1.Z() * V2.Y() - V1.Y() * V2.Z();
170 Det3=V1.Z() * V2.X() - V1.X() * V2.Z();
171
172 if((Det1!=0.0) && ((Abs(Det1) >= Abs(Det2))&&(Abs(Det1) >= Abs(Det3)))) {
173 A=(smy * V2.X() - smx * V2.Y())/Det1;
174 }
175 else if((Det2!=0.0)
7eed5d29 176 && ((Abs(Det2) >= Abs(Det1))
177 &&(Abs(Det2) >= Abs(Det3)))) {
7fd59977 178 A=(smz * V2.Y() - smy * V2.Z())/Det2;
179 }
180 else {
181 A=(smz * V2.X() - smx * V2.Z())/Det3;
182 }
183 ptintersect.SetCoord( P1.X() + A * V1.X()
7eed5d29 184 ,P1.Y() + A * V1.Y()
185 ,P1.Z() + A * V1.Z());
7fd59977 186 }
187 else {
188 ptintersect.SetCoord(0,0,0); //-- Pour eviter des FPE
189 }
190}
191//=======================================================================
192//function : Distance
193//purpose :
194//=======================================================================
a34f083b 195void AxeOperator::Distance(Standard_Real& dist,
196 Standard_Real& Param1,
197 Standard_Real& Param2)
7fd59977 198 {
a34f083b 199 gp_Vec O1O2(Axe1.Location(),Axe2.Location());
7fd59977 200 gp_Dir U1 = Axe1.Direction(); //-- juste pour voir.
201 gp_Dir U2 = Axe2.Direction();
202
203 gp_Dir N = U1.Crossed(U2);
204 Standard_Real D = Det33(U1.X(),U2.X(),N.X(),
7eed5d29 205 U1.Y(),U2.Y(),N.Y(),
206 U1.Z(),U2.Z(),N.Z());
7fd59977 207 if(D) {
208 dist = Det33(U1.X(),U2.X(),O1O2.X(),
7eed5d29 209 U1.Y(),U2.Y(),O1O2.Y(),
210 U1.Z(),U2.Z(),O1O2.Z()) / D;
7fd59977 211 Param1 = Det33(O1O2.X(),U2.X(),N.X(),
7eed5d29 212 O1O2.Y(),U2.Y(),N.Y(),
213 O1O2.Z(),U2.Z(),N.Z()) / (-D);
7fd59977 214 //------------------------------------------------------------
215 //-- On resout P1 * Dir1 + P2 * Dir2 + d * N = O1O2
216 //-- soit : Segment perpendiculaire : O1+P1 D1
217 //-- O2-P2 D2
218 Param2 = Det33(U1.X(),O1O2.X(),N.X(),
7eed5d29 219 U1.Y(),O1O2.Y(),N.Y(),
220 U1.Z(),O1O2.Z(),N.Z()) / (D);
7fd59977 221 }
222}
223//=======================================================================
224//function : DirToAx2
225//purpose : returns a gp_Ax2 where D is the main direction
226//=======================================================================
227gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D)
228{
229 Standard_Real x=D.X(); Standard_Real ax=Abs(x);
230 Standard_Real y=D.Y(); Standard_Real ay=Abs(y);
231 Standard_Real z=D.Z(); Standard_Real az=Abs(z);
232 if( (ax==0.0) || ((ax<ay) && (ax<az)) ) {
233 return(gp_Ax2(P,D,gp_Dir(gp_Vec(0.0,-z,y))));
234 }
235 else if( (ay==0.0) || ((ay<ax) && (ay<az)) ) {
236 return(gp_Ax2(P,D,gp_Dir(gp_Vec(-z,0.0,x))));
237 }
238 else {
239 return(gp_Ax2(P,D,gp_Dir(gp_Vec(-y,x,0.0))));
240 }
241}
242//=======================================================================
243//function : IntAna_QuadQuadGeo
244//purpose : Empty constructor
245//=======================================================================
a34f083b 246IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(void)
7fd59977 247 : done(Standard_False),
248 nbint(0),
249 typeres(IntAna_Empty),
250 pt1(0,0,0),
251 pt2(0,0,0),
7eed5d29 252 pt3(0,0,0),
253 pt4(0,0,0),
7fd59977 254 param1(0),
255 param2(0),
7eed5d29 256 param3(0),
257 param4(0),
7fd59977 258 param1bis(0),
259 param2bis(0),
260 myCommonGen(Standard_False),
261 myPChar(0,0,0)
262{
263 InitTolerances();
264}
265//=======================================================================
266//function : InitTolerances
267//purpose :
268//=======================================================================
a34f083b 269void IntAna_QuadQuadGeo::InitTolerances()
7fd59977 270{
271 myEPSILON_DISTANCE = 0.00000000000001;
272 myEPSILON_ANGLE_CONE = 0.000000000001;
273 myEPSILON_MINI_CIRCLE_RADIUS = 0.000000001;
274 myEPSILON_CYLINDER_DELTA_RADIUS = 0.0000000000001;
275 myEPSILON_CYLINDER_DELTA_DISTANCE= 0.0000001;
276 myEPSILON_AXES_PARA = 0.000000000001;
277}
278//=======================================================================
279//function : IntAna_QuadQuadGeo
280//purpose : Pln Pln
281//=======================================================================
a34f083b 282IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P1,
283 const gp_Pln& P2,
284 const Standard_Real TolAng,
285 const Standard_Real Tol)
7fd59977 286: done(Standard_False),
287 nbint(0),
288 typeres(IntAna_Empty),
289 pt1(0,0,0),
290 pt2(0,0,0),
7eed5d29 291 pt3(0,0,0),
292 pt4(0,0,0),
7fd59977 293 param1(0),
294 param2(0),
7eed5d29 295 param3(0),
296 param4(0),
7fd59977 297 param1bis(0),
298 param2bis(0),
299 myCommonGen(Standard_False),
300 myPChar(0,0,0)
301{
302 InitTolerances();
303 Perform(P1,P2,TolAng,Tol);
304}
305//=======================================================================
306//function : Perform
307//purpose :
308//=======================================================================
a34f083b 309void IntAna_QuadQuadGeo::Perform (const gp_Pln& P1,
310 const gp_Pln& P2,
311 const Standard_Real TolAng,
312 const Standard_Real Tol)
7fd59977 313{
a34f083b 314 Standard_Real A1, B1, C1, D1, A2, B2, C2, D2, dist1, dist2, aMVD;
315 //
7fd59977 316 done=Standard_False;
a34f083b 317 param2bis=0.;
7fd59977 318 //
7fd59977 319 P1.Coefficients(A1,B1,C1,D1);
320 P2.Coefficients(A2,B2,C2,D2);
a34f083b 321 //
322 gp_Vec aVN1(A1,B1,C1);
323 gp_Vec aVN2(A2,B2,C2);
324 gp_Vec vd(aVN1.Crossed(aVN2));
325 //
326 const gp_Pnt& aLocP1=P1.Location();
327 const gp_Pnt& aLocP2=P2.Location();
328 //
329 dist1=A2*aLocP1.X() + B2*aLocP1.Y() + C2*aLocP1.Z() + D2;
330 dist2=A1*aLocP2.X() + B1*aLocP2.Y() + C1*aLocP2.Z() + D1;
331 //
332 aMVD=vd.Magnitude();
333 if(aMVD <=TolAng) {
7fd59977 334 // normalles are collinear - planes are same or parallel
a34f083b 335 typeres = (Abs(dist1) <= Tol && Abs(dist2) <= Tol) ? IntAna_Same
336 : IntAna_Empty;
7fd59977 337 }
338 else {
a34f083b 339 Standard_Real denom, denom2, ddenom, par1, par2;
340 Standard_Real X1, Y1, Z1, X2, Y2, Z2, aEps;
341 //
342 aEps=1.e-16;
343 denom=A1*A2 + B1*B2 + C1*C2;
344 denom2 = denom*denom;
345 ddenom = 1. - denom2;
346
347 denom = ( Abs(ddenom) <= aEps ) ? aEps : ddenom;
7fd59977 348
a34f083b 349 par1 = dist1/denom;
350 par2 = -dist2/denom;
7fd59977 351
a34f083b 352 gp_Vec inter1(aVN1.Crossed(vd));
353 gp_Vec inter2(aVN2.Crossed(vd));
7fd59977 354
a34f083b 355 X1=aLocP1.X() + par1*inter1.X();
356 Y1=aLocP1.Y() + par1*inter1.Y();
357 Z1=aLocP1.Z() + par1*inter1.Z();
358 X2=aLocP2.X() + par2*inter2.X();
359 Y2=aLocP2.Y() + par2*inter2.Y();
360 Z2=aLocP2.Z() + par2*inter2.Z();
7fd59977 361
362 pt1=gp_Pnt((X1+X2)*0.5, (Y1+Y2)*0.5, (Z1+Z2)*0.5);
363 dir1 = gp_Dir(vd);
364 typeres = IntAna_Line;
365 nbint = 1;
a34f083b 366 //
367 //-------------------------------------------------------
368 // When the value of the angle between the planes is small
369 // the origin of intersection line is computed with error
370 // [ ~0.0001 ] that can not br considered as small one
371 // e.g.
372 // for {A~=2.e-6, dist1=4.2e-5, dist2==1.e-4} =>
373 // {denom=3.4e-12, par1=12550297.6, par2=32605552.9, etc}
374 // So,
375 // the origin should be refined if it is possible
376 //
377 Standard_Real aTreshAng, aTreshDist;
378 //
379 aTreshAng=2.e-6; // 1.e-4 deg
380 aTreshDist=1.e-12;
381 //
382 if (aMVD < aTreshAng) {
383 Standard_Real aDist1, aDist2;
384 //
385 aDist1=A1*pt1.X() + B1*pt1.Y() + C1*pt1.Z() + D1;
386 aDist2=A2*pt1.X() + B2*pt1.Y() + C2*pt1.Z() + D2;
387 //
388 if (fabs(aDist1)>aTreshDist || fabs(aDist2)>aTreshDist) {
389 Standard_Boolean bIsDone, bIsParallel;
390 IntAna_IntConicQuad aICQ;
391 //
392 // 1.
393 gp_Dir aDN1(aVN1);
394 gp_Lin aL1(pt1, aDN1);
395 //
396 aICQ.Perform(aL1, P1, TolAng, Tol);
397 bIsDone=aICQ.IsDone();
398 if (!bIsDone) {
399 return;
400 }
401 //
402 const gp_Pnt& aPnt1=aICQ.Point(1);
403 //----------------------------------
404 // 2.
405 gp_Dir aDL2(dir1.Crossed(aDN1));
406 gp_Lin aL2(aPnt1, aDL2);
407 //
408 aICQ.Perform(aL2, P2, TolAng, Tol);
409 bIsDone=aICQ.IsDone();
410 if (!bIsDone) {
411 return;
412 }
413 //
414 bIsParallel=aICQ.IsParallel();
415 if (bIsParallel) {
416 return;
417 }
418 //
419 const gp_Pnt& aPnt2=aICQ.Point(1);
420 //
421 pt1=aPnt2;
422 }
423 }
7fd59977 424 }
425 done=Standard_True;
426}
427//=======================================================================
428//function : IntAna_QuadQuadGeo
429//purpose : Pln Cylinder
430//=======================================================================
a34f083b 431IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Pln& P
432 ,const gp_Cylinder& Cl
433 ,const Standard_Real Tolang
434 ,const Standard_Real Tol
435 ,const Standard_Real H)
436 : done(Standard_False),
437 nbint(0),
438 typeres(IntAna_Empty),
439 pt1(0,0,0),
440 pt2(0,0,0),
441 pt3(0,0,0),
442 pt4(0,0,0),
443 param1(0),
444 param2(0),
445 param3(0),
446 param4(0),
447 param1bis(0),
448 param2bis(0),
449 myCommonGen(Standard_False),
450 myPChar(0,0,0)
7fd59977 451{
452 InitTolerances();
04cbc9d3 453 Perform(P,Cl,Tolang,Tol,H);
7fd59977 454}
455//=======================================================================
456//function : Perform
457//purpose :
458//=======================================================================
459 void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
04cbc9d3 460 ,const gp_Cylinder& Cl
461 ,const Standard_Real Tolang
462 ,const Standard_Real Tol
463 ,const Standard_Real H)
7fd59977 464{
465 done = Standard_False;
466 Standard_Real dist,radius;
467 Standard_Real A,B,C,D;
468 Standard_Real X,Y,Z;
469 Standard_Real sint,cost,h;
470 gp_XYZ axex,axey,omega;
471
472
473 param2bis=0.0;
474 radius = Cl.Radius();
475
476 gp_Lin axec(Cl.Axis());
477 gp_XYZ normp(P.Axis().Direction().XYZ());
478
479 P.Coefficients(A,B,C,D);
480 axec.Location().Coord(X,Y,Z);
a34f083b 481 // la distance axe/plan est evaluee a l origine de l axe.
482 dist = A*X + B*Y + C*Z + D;
7fd59977 483
484 Standard_Real tolang = Tolang;
485 Standard_Boolean newparams = Standard_False;
486
487 gp_Vec ldv( axec.Direction() );
488 gp_Vec npv( normp );
489 Standard_Real dA = Abs( ldv.Angle( npv ) );
c6541a0c 490 if( dA > (M_PI/4.) )
7fd59977 491 {
c6541a0c 492 Standard_Real dang = Abs( ldv.Angle( npv ) ) - M_PI/2.;
7fd59977 493 Standard_Real dangle = Abs( dang );
494 if( dangle > Tolang )
7eed5d29 495 {
496 Standard_Real sinda = Abs( Sin( dangle ) );
497 Standard_Real dif = Abs( sinda - Tol );
498 if( dif < Tol )
499 {
500 tolang = sinda * 2.;
501 newparams = Standard_True;
502 }
503 }
7fd59977 504 }
505
506 nbint = 0;
04cbc9d3 507 IntAna_IntConicQuad inter(axec,P,tolang,Tol,H);
7fd59977 508
509 if (inter.IsParallel()) {
510 // Le resultat de l intersection Plan-Cylindre est de type droite.
511 // il y a 1 ou 2 droites
512
513 typeres = IntAna_Line;
514 omega.SetCoord(X-dist*A,Y-dist*B,Z-dist*C);
515
516 if (Abs(Abs(dist)-radius) < Tol)
517 {
7eed5d29 518 nbint = 1;
519 pt1.SetXYZ(omega);
520
521 if( newparams )
522 {
523 gp_XYZ omegaXYZ(X,Y,Z);
524 gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
525 Standard_Real Xt,Yt,Zt,distt;
526 omegaXYZtrnsl.Coord(Xt,Yt,Zt);
527 distt = A*Xt + B*Yt + C*Zt + D;
a34f083b 528 gp_XYZ omega1(omegaXYZtrnsl.X()-distt*A,
529 omegaXYZtrnsl.Y()-distt*B,
530 omegaXYZtrnsl.Z()-distt*C );
7eed5d29 531 gp_Pnt ppt1;
532 ppt1.SetXYZ( omega1 );
533 gp_Vec vv1(pt1,ppt1);
534 gp_Dir dd1( vv1 );
535 dir1 = dd1;
536 }
537 else
538 dir1 = axec.Direction();
7fd59977 539 }
540 else if (Abs(dist) < radius)
541 {
7eed5d29 542 nbint = 2;
543 h = Sqrt(radius*radius - dist*dist);
544 axey = axec.Direction().XYZ().Crossed(normp); // axey est normalise
545
546 pt1.SetXYZ(omega - h*axey);
547 pt2.SetXYZ(omega + h*axey);
548
549 if( newparams )
550 {
551 gp_XYZ omegaXYZ(X,Y,Z);
552 gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
553 Standard_Real Xt,Yt,Zt,distt,ht;
554 omegaXYZtrnsl.Coord(Xt,Yt,Zt);
555 distt = A*Xt + B*Yt + C*Zt + D;
556 // ht = Sqrt(radius*radius - distt*distt);
557 Standard_Real anSqrtArg = radius*radius - distt*distt;
558 ht = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
559
a34f083b 560 gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A,
561 omegaXYZtrnsl.Y()-distt*B,
562 omegaXYZtrnsl.Z()-distt*C );
7eed5d29 563 gp_Pnt ppt1,ppt2;
564 ppt1.SetXYZ( omega1 - ht*axey);
565 ppt2.SetXYZ( omega1 + ht*axey);
566 gp_Vec vv1(pt1,ppt1);
567 gp_Vec vv2(pt2,ppt2);
568 gp_Dir dd1( vv1 );
569 gp_Dir dd2( vv2 );
570 dir1 = dd1;
571 dir2 = dd2;
572 }
573 else
574 {
575 dir1 = axec.Direction();
576 dir2 = axec.Direction();
577 }
7fd59977 578 }
579 // else nbint = 0
580
581 // debug JAG : le nbint = 0 doit etre remplace par typeres = IntAna_Empty
582 // et ne pas etre seulement supprime...
583
584 else {
585 typeres = IntAna_Empty;
586 }
587 }
588 else { // Il y a un point d intersection. C est le centre du cercle
589 // ou de l ellipse solution.
590
591 nbint = 1;
592 axey = normp.Crossed(axec.Direction().XYZ());
593 sint = axey.Modulus();
594
595 pt1 = inter.Point(1);
596
597 if (sint < Tol/radius) {
598
599 // on construit un cercle avec comme axes X et Y ceux du cylindre
600 typeres = IntAna_Circle;
601
602 dir1 = axec.Direction(); // axe Z
603 dir2 = Cl.Position().XDirection();
604 param1 = radius;
605 }
606 else {
607
608 // on construit un ellipse
609 typeres = IntAna_Ellipse;
610 cost = Abs(axec.Direction().XYZ().Dot(normp));
611 axex = axey.Crossed(normp);
612
613 dir1.SetXYZ(normp); //Modif ds ce bloc
614 dir2.SetXYZ(axex);
615
616 param1 = radius/cost;
617 param1bis = radius;
618 }
619 }
788cbaf4 620
7fd59977 621 done = Standard_True;
622}
623//=======================================================================
624//function : IntAna_QuadQuadGeo
625//purpose : Pln Cone
626//=======================================================================
627 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
7eed5d29 628 const gp_Cone& Co,
629 const Standard_Real Tolang,
630 const Standard_Real Tol)
7fd59977 631: done(Standard_False),
632 nbint(0),
633 typeres(IntAna_Empty),
634 pt1(0,0,0),
635 pt2(0,0,0),
7eed5d29 636 pt3(0,0,0),
637 pt4(0,0,0),
7fd59977 638 param1(0),
639 param2(0),
7eed5d29 640 param3(0),
641 param4(0),
7fd59977 642 param1bis(0),
643 param2bis(0),
644 myCommonGen(Standard_False),
645 myPChar(0,0,0)
646{
647 InitTolerances();
648 Perform(P,Co,Tolang,Tol);
649}
650//=======================================================================
651//function : Perform
652//purpose :
653//=======================================================================
654 void IntAna_QuadQuadGeo::Perform(const gp_Pln& P,
7eed5d29 655 const gp_Cone& Co,
656 const Standard_Real Tolang,
657 const Standard_Real Tol)
7fd59977 658{
659
660 done = Standard_False;
661 nbint = 0;
662
663 Standard_Real A,B,C,D;
664 Standard_Real X,Y,Z;
665 Standard_Real dist,sint,cost,sina,cosa,angl,costa;
666 Standard_Real dh;
667 gp_XYZ axex,axey;
668
669 gp_Lin axec(Co.Axis());
670 P.Coefficients(A,B,C,D);
671 gp_Pnt apex(Co.Apex());
672
673 apex.Coord(X,Y,Z);
674 dist = A*X + B*Y + C*Z + D; // distance signee sommet du cone/ Plan
675
676 gp_XYZ normp = P.Axis().Direction().XYZ();
677 if(P.Direct()==Standard_False) { //-- lbr le 14 jan 97
678 normp.Reverse();
679 }
680
681 axey = normp.Crossed(Co.Axis().Direction().XYZ());
682 axex = axey.Crossed(normp);
683
684
685 angl = Co.SemiAngle();
686
687 cosa = Cos(angl);
688 sina = Abs(Sin(angl));
689
690
691 // Angle entre la normale au plan et l axe du cone, ramene entre 0. et PI/2.
692
693 sint = axey.Modulus();
694 cost = Abs(Co.Axis().Direction().XYZ().Dot(normp));
695
696 // Le calcul de costa permet de determiner si le plan contient
697 // un generatrice du cone : on calcul Sin((PI/2. - t) - angl)
698
699 costa = cost*cosa - sint*sina; // sin((PI/2 -t)-angl)=cos(t+angl)
700 // avec t ramene entre 0 et pi/2.
701
702 if (Abs(dist) < Tol) {
703 // on considere que le plan contient le sommet du cone.
704 // les solutions possibles sont donc : 1 point, 1 droite, 2 droites
705 // selon l inclinaison du plan.
706
707 if (Abs(costa) < Tolang) { // plan parallele a la generatrice
708 typeres = IntAna_Line;
709 nbint = 1;
710 gp_XYZ ptonaxe(apex.XYZ() + 10.*(Co.Axis().Direction().XYZ()));
711 // point sur l axe du cone cote z positif
712
713 dist = A*ptonaxe.X() + B*ptonaxe.Y() + C*ptonaxe.Z() + D;
714 ptonaxe = ptonaxe - dist*normp;
715 pt1 = apex;
716 dir1.SetXYZ(ptonaxe - pt1.XYZ());
717 }
718 else if (cost < sina) { // plan "interieur" au cone
719 typeres = IntAna_Line;
720 nbint = 2;
721 pt1 = apex;
722 pt2 = apex;
723 dh = Sqrt(sina*sina-cost*cost)/cosa;
724 dir1.SetXYZ(axex + dh*axey);
725 dir2.SetXYZ(axex - dh*axey);
726 }
727 else { // plan "exterieur" au cone
728 typeres = IntAna_Point;
729 nbint = 1;
730 pt1 = apex;
731 }
732 }
733 else {
734 // Solutions possibles : cercle, ellipse, parabole, hyperbole selon
735 // l inclinaison du plan.
736 Standard_Real deltacenter, distance;
737
738 if (cost < Tolang) {
739 // Le plan contient la direction de l axe du cone. La solution est
740 // l hyperbole
741 typeres = IntAna_Hyperbola;
742 nbint = 2;
743 pt1.SetXYZ(apex.XYZ()-dist*normp);
744 pt2 = pt1;
745 dir1=normp;
746 dir2.SetXYZ(axex);
747 param1 = param2 = Abs(dist/Tan(angl));
748 param1bis = param2bis = Abs(dist);
749 }
750 else {
751
752 IntAna_IntConicQuad inter(axec,P,Tolang); // on a necessairement 1 point.
753
754 gp_Pnt center(inter.Point(1));
755
756 // En fonction de la position de l intersection par rapport au sommet
757 // du cone, on change l axe x en -x et y en -y. Le parametre du sommet
758 // sur axec est negatif (voir definition du cone)
759
760 distance = apex.Distance(center);
761
762 if (inter.ParamOnConic(1) + Co.RefRadius()/Tan(angl) < 0.) {
7eed5d29 763 axex.Reverse();
764 axey.Reverse();
7fd59977 765 }
766
767 if (Abs(costa) < Tolang) { // plan parallele a une generatrice
7eed5d29 768 typeres = IntAna_Parabola;
769 nbint = 1;
770 deltacenter = distance/2./cosa;
771 axex.Normalize();
772 pt1.SetXYZ(center.XYZ()-deltacenter*axex);
773 dir1 = normp;
774 dir2.SetXYZ(axex);
775 param1 = deltacenter*sina*sina;
7fd59977 776 }
777 else if (sint < Tolang) { // plan perpendiculaire a l axe
7eed5d29 778 typeres = IntAna_Circle;
779 nbint = 1;
780 pt1 = center;
781 dir1 = Co.Position().Direction();
782 dir2 = Co.Position().XDirection();
783 param1 = apex.Distance(center)*Abs(Tan(angl));
7fd59977 784 }
785 else if (cost < sina ) {
7eed5d29 786 typeres = IntAna_Hyperbola;
787 nbint = 2;
788 axex.Normalize();
789
790 deltacenter = sint*sina*sina*distance/(sina*sina - cost*cost);
791 pt1.SetXYZ(center.XYZ() - deltacenter*axex);
792 pt2 = pt1;
793 dir1 = normp;
794 dir2.SetXYZ(axex);
795 param1 = param2 = cost*sina*cosa*distance /(sina*sina-cost*cost);
796 param1bis = param2bis = cost*sina*distance / Sqrt(sina*sina-cost*cost);
7fd59977 797
798 }
799 else { // on a alors cost > sina
7eed5d29 800 typeres = IntAna_Ellipse;
801 nbint = 1;
802 Standard_Real radius = cost*sina*cosa*distance/(cost*cost-sina*sina);
803 deltacenter = sint*sina*sina*distance/(cost*cost-sina*sina);
804 axex.Normalize();
805 pt1.SetXYZ(center.XYZ() + deltacenter*axex);
806 dir1 = normp;
807 dir2.SetXYZ(axex);
808 param1 = radius;
809 param1bis = cost*sina*distance/ Sqrt(cost*cost - sina*sina);
7fd59977 810 }
811 }
812 }
813
814 //-- On a du mal a gerer plus loin (Value ProjLib, Params ... )
815 //-- des hyperboles trop bizarres
816 //-- On retourne False -> Traitement par biparametree
817 static Standard_Real EllipseLimit = 1.0E+9; //OCC513(apo) 1000000
818 static Standard_Real HyperbolaLimit = 2.0E+6; //OCC537(apo) 50000
819 if(typeres==IntAna_Ellipse && nbint>=1) {
820 if(Abs(param1) > EllipseLimit || Abs(param1bis) > EllipseLimit) {
821 done=Standard_False;
822 return;
823 }
824 }
825 if(typeres==IntAna_Hyperbola && nbint>=2) {
826 if(Abs(param2) > HyperbolaLimit || Abs(param2bis) > HyperbolaLimit) {
827 done = Standard_False;
828 return;
829 }
830 }
831 if(typeres==IntAna_Hyperbola && nbint>=1) {
832 if(Abs(param1) > HyperbolaLimit || Abs(param1bis) > HyperbolaLimit) {
833 done=Standard_False;
834 return;
835 }
836 }
837
838 done = Standard_True;
839}
840
841//=======================================================================
842//function : IntAna_QuadQuadGeo
843//purpose : Pln Sphere
844//=======================================================================
a34f083b 845IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
846 const gp_Sphere& S)
7fd59977 847: done(Standard_False),
848 nbint(0),
849 typeres(IntAna_Empty),
850 pt1(0,0,0),
851 pt2(0,0,0),
7eed5d29 852 pt3(0,0,0),
853 pt4(0,0,0),
7fd59977 854 param1(0),
855 param2(0),
7eed5d29 856 param3(0),
857 param4(0),
7fd59977 858 param1bis(0),
859 param2bis(0),
860 myCommonGen(Standard_False),
861 myPChar(0,0,0)
862{
863 InitTolerances();
864 Perform(P,S);
865}
866//=======================================================================
867//function : Perform
868//purpose :
869//=======================================================================
a34f083b 870void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
871 ,const gp_Sphere& S)
7fd59977 872{
873
874 done = Standard_False;
875 Standard_Real A,B,C,D,dist, radius;
876 Standard_Real X,Y,Z;
877
878 nbint = 0;
879// debug JAG : on met typeres = IntAna_Empty par defaut...
880 typeres = IntAna_Empty;
881
882 P.Coefficients(A,B,C,D);
883 S.Location().Coord(X,Y,Z);
884 radius = S.Radius();
885
886 dist = A * X + B * Y + C * Z + D;
887
888 if (Abs( Abs(dist) - radius) < Epsilon(radius)) {
889 // on a une seule solution : le point projection du centre de la sphere
890 // sur le plan
891 nbint = 1;
892 typeres = IntAna_Point;
893 pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
894 }
895 else if (Abs(dist) < radius) {
896 // on a un cercle solution
897 nbint = 1;
898 typeres = IntAna_Circle;
899 pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
900 dir1 = P.Axis().Direction();
901 if(P.Direct()==Standard_False) dir1.Reverse();
902 dir2 = P.Position().XDirection();
903 param1 = Sqrt(radius*radius - dist*dist);
904 }
905 param2bis=0.0; //-- pour eviter param2bis not used ....
906 done = Standard_True;
907}
908
909//=======================================================================
910//function : IntAna_QuadQuadGeo
911//purpose : Cylinder - Cylinder
912//=======================================================================
a34f083b 913IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl1,
914 const gp_Cylinder& Cyl2,
915 const Standard_Real Tol)
7fd59977 916: done(Standard_False),
917 nbint(0),
918 typeres(IntAna_Empty),
919 pt1(0,0,0),
920 pt2(0,0,0),
7eed5d29 921 pt3(0,0,0),
922 pt4(0,0,0),
7fd59977 923 param1(0),
924 param2(0),
7eed5d29 925 param3(0),
926 param4(0),
7fd59977 927 param1bis(0),
928 param2bis(0),
929 myCommonGen(Standard_False),
930 myPChar(0,0,0)
931{
932 InitTolerances();
933 Perform(Cyl1,Cyl2,Tol);
934}
935//=======================================================================
936//function : Perform
937//purpose :
938//=======================================================================
a34f083b 939void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
940 const gp_Cylinder& Cyl2,
941 const Standard_Real Tol)
7fd59977 942{
943 done=Standard_True;
944 //---------------------------- Parallel axes -------------------------
945 AxeOperator A1A2(Cyl1.Axis(),Cyl2.Axis());
946 Standard_Real R1=Cyl1.Radius();
947 Standard_Real R2=Cyl2.Radius();
948 Standard_Real RmR, RmR_Relative;
949 RmR=(R1>R2)? (R1-R2) : (R2-R1);
950 {
96a95605 951 Standard_Real Rmax;
7fd59977 952 Rmax=(R1>R2)? R1 : R2;
7fd59977 953 RmR_Relative=RmR/Rmax;
954 }
955
956 Standard_Real DistA1A2=A1A2.Distance();
957
ecc4f148 958 if(A1A2.Parallel())
959 {
960 if(DistA1A2<=Tol)
961 {
962 if(RmR<=Tol)
963 {
7eed5d29 964 typeres=IntAna_Same;
7fd59977 965 }
ecc4f148 966 else
967 {
7eed5d29 968 typeres=IntAna_Empty;
7fd59977 969 }
970 }
ecc4f148 971 else
972 { //-- DistA1A2 > Tol
7fd59977 973 gp_Pnt P1=Cyl1.Location();
974 gp_Pnt P2t=Cyl2.Location();
975 gp_Pnt P2;
976 //-- P2t is projected on the plane (P1,DirCylX,DirCylY)
977 gp_Dir DirCyl = Cyl1.Position().Direction();
978 Standard_Real ProjP2OnDirCyl1=gp_Vec(DirCyl).Dot(gp_Vec(P1,P2t));
979
b70d2b09 980 //P2 is a projection the location of the 2nd cylinder on the base
981 //of the 1st cylinder
ecc4f148 982 P2.SetCoord(P2t.X() - ProjP2OnDirCyl1*DirCyl.X(),
983 P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y(),
984 P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
7fd59977 985 //--
986 Standard_Real R1pR2=R1+R2;
ecc4f148 987 if(DistA1A2>(R1pR2+Tol))
988 {
7eed5d29 989 typeres=IntAna_Empty;
990 nbint=0;
7fd59977 991 }
b70d2b09 992 else if((R1pR2 - DistA1A2) <= RealSmall())
ecc4f148 993 {
7eed5d29 994 //-- 1 Tangent line -------------------------------------OK
995 typeres=IntAna_Line;
996
997 nbint=1;
998 dir1=DirCyl;
999 Standard_Real R1_R1pR2=R1/R1pR2;
ecc4f148 1000 pt1.SetCoord( P1.X() + R1_R1pR2 * (P2.X()-P1.X()),
1001 P1.Y() + R1_R1pR2 * (P2.Y()-P1.Y()),
1002 P1.Z() + R1_R1pR2 * (P2.Z()-P1.Z()));
7fd59977 1003 }
ecc4f148 1004 else if(DistA1A2>RmR)
1005 {
7eed5d29 1006 //-- 2 lines ---------------------------------------------OK
1007 typeres=IntAna_Line;
1008 nbint=2;
1009 dir1=DirCyl;
7eed5d29 1010 dir2=dir1;
7eed5d29 1011
b70d2b09 1012 const Standard_Real aR1R1 = R1*R1;
1013
1014 /*
1015 P1
1016 o
1017 * | *
1018 * O1| *
1019 A o-----o-----o B
1020 * | *
1021 * | *
1022 o
1023 P2
1024
1025 Two cylinders have axes collinear. Therefore, problem can be reformulated as
1026 to find intersection point of two circles (the bases of the cylinders) on
1027 the plane: 1st circle has center P1 and radius R1 (the radius of the
1028 1st cylinder) and 2nd circle has center P2 and radius R2 (the radius of the
1029 2nd cylinder). The plane is the base of the 1st cylinder. Points A and B
1030 are intersection point of these circles. Distance P1P2 is equal to DistA1A2.
1031 O1 is the intersection point of P1P2 and AB segments.
1032
1033 At that, if distance AB < Tol we consider that the circles are tangent and
1034 has only one intersection point.
1035
1036 AB = 2*R1*sin(angle AP1P2).
1037 Accordingly,
1038 AB^2 < Tol^2 => 4*R1*R1*sin(angle AP1P2)^2 < Tol*Tol.
1039 */
1040
7eed5d29 1041
b70d2b09 1042 //Cosine and Square of Sine of the A-P1-P2 angle
1043 const Standard_Real aCos = 0.5*(aR1R1-R2*R2+DistA1A2*DistA1A2)/(R1*DistA1A2);
1044 const Standard_Real aSin2 = 1-aCos*aCos;
1045
1046 const Standard_Boolean isTangent =((4.0*aR1R1*aSin2) < Tol*Tol);
1047
1048 //Normalized vector P1P2
1049 const gp_Vec DirA1A2((P2.XYZ() - P1.XYZ())/DistA1A2);
1050
1051 if(isTangent)
ecc4f148 1052 {
b70d2b09 1053 //Intercept the segment from P1 point along P1P2 direction
1054 //and having |P1O1| length
7eed5d29 1055 nbint=1;
b70d2b09 1056 pt1.SetXYZ(P1.XYZ() + DirA1A2.XYZ()*R1*aCos);
7eed5d29 1057 }
ecc4f148 1058 else
b70d2b09 1059 {
1060 //Sine of the A-P1-P2 angle (if aSin2 < 0 then isTangent == TRUE =>
1061 //go to another branch)
1062 const Standard_Real aSin = sqrt(aSin2);
1063
1064 //1. Rotate P1P2 to the angle A-P1-P2 relative to P1
1065 //(clockwise and anticlockwise for getting
1066 //two intersection points).
1067 //2. Intercept the segment from P1 along direction,
1068 //determined in the preview paragraph and having R1 length
1069 const gp_Dir &aXDir = Cyl1.Position().XDirection(),
1070 &aYDir = Cyl1.Position().YDirection();
1071 const gp_Vec aR1Xdir = R1*aXDir.XYZ(),
1072 aR1Ydir = R1*aYDir.XYZ();
1073
1074 //Source 2D-coordinates of the P1P2 vector normalized
1075 //in coordinate system, based on the X- and Y-directions
1076 //of the 1st cylinder in the plane of the 1st cylinder base
1077 //(P1 is the origin of the coordinate system).
1078 const Standard_Real aDx = DirA1A2.Dot(aXDir),
1079 aDy = DirA1A2.Dot(aYDir);
1080
1081 //New coordinate (after rotation) of the P1P2 vector normalized.
1082 Standard_Real aNewDx = aDx*aCos - aDy*aSin,
1083 aNewDy = aDy*aCos + aDx*aSin;
1084 pt1.SetXYZ(P1.XYZ() + aNewDx*aR1Xdir.XYZ() + aNewDy*aR1Ydir.XYZ());
1085
1086 aNewDx = aDx*aCos + aDy*aSin;
1087 aNewDy = aDy*aCos - aDx*aSin;
1088 pt2.SetXYZ(P1.XYZ() + aNewDx*aR1Xdir.XYZ() + aNewDy*aR1Ydir.XYZ());
7eed5d29 1089 }
7fd59977 1090 }
ecc4f148 1091 else if(DistA1A2>(RmR-Tol))
1092 {
7eed5d29 1093 //-- 1 Tangent ------------------------------------------OK
1094 typeres=IntAna_Line;
1095 nbint=1;
1096 dir1=DirCyl;
1097 Standard_Real R1_RmR=R1/RmR;
7fd59977 1098
ecc4f148 1099 if(R1 < R2)
1100 R1_RmR = -R1_RmR;
7fd59977 1101
ecc4f148 1102 pt1.SetCoord( P1.X() + R1_RmR * (P2.X()-P1.X()),
1103 P1.Y() + R1_RmR * (P2.Y()-P1.Y()),
1104 P1.Z() + R1_RmR * (P2.Z()-P1.Z()));
7fd59977 1105 }
1106 else {
7eed5d29 1107 nbint=0;
1108 typeres=IntAna_Empty;
7fd59977 1109 }
1110 }
1111 }
1112 else { //-- No Parallel Axis ---------------------------------OK
1113 if((RmR_Relative<=myEPSILON_CYLINDER_DELTA_RADIUS)
ecc4f148 1114 && (DistA1A2 <= myEPSILON_CYLINDER_DELTA_DISTANCE))
1115 {
7fd59977 1116 //-- PI/2 between the two axis and Intersection
1117 //-- and identical radius
1118 typeres=IntAna_Ellipse;
1119 nbint=2;
1120 gp_Dir DirCyl1=Cyl1.Position().Direction();
1121 gp_Dir DirCyl2=Cyl2.Position().Direction();
1122 pt1=pt2=A1A2.PtIntersect();
1123
1124 Standard_Real A=DirCyl1.Angle(DirCyl2);
1125 Standard_Real B;
c6541a0c 1126 B=Abs(Sin(0.5*(M_PI-A)));
7fd59977 1127 A=Abs(Sin(0.5*A));
1128
ecc4f148 1129 if(A==0.0 || B==0.0)
1130 {
7eed5d29 1131 typeres=IntAna_Same;
1132 return;
7fd59977 1133 }
1134
7fd59977 1135 gp_Vec dircyl1(DirCyl1);gp_Vec dircyl2(DirCyl2);
1136 dir1 = gp_Dir(dircyl1.Added(dircyl2));
1137 dir2 = gp_Dir(dircyl1.Subtracted(dircyl2));
7eed5d29 1138
7fd59977 1139 param2 = Cyl1.Radius() / A;
1140 param1 = Cyl1.Radius() / B;
1141 param2bis= param1bis = Cyl1.Radius();
ecc4f148 1142 if(param1 < param1bis)
1143 {
1144 A=param1;
1145 param1=param1bis;
1146 param1bis=A;
7fd59977 1147 }
ecc4f148 1148
1149 if(param2 < param2bis)
1150 {
1151 A=param2;
1152 param2=param2bis;
1153 param2bis=A;
7fd59977 1154 }
1155 }
ecc4f148 1156 else
1157 {
1158 if(Abs(DistA1A2-Cyl1.Radius()-Cyl2.Radius())<Tol)
1159 {
7eed5d29 1160 typeres = IntAna_Point;
1161 Standard_Real d,p1,p2;
1162
1163 gp_Dir D1 = Cyl1.Axis().Direction();
1164 gp_Dir D2 = Cyl2.Axis().Direction();
1165 A1A2.Distance(d,p1,p2);
1166 gp_Pnt P = Cyl1.Axis().Location();
1167 gp_Pnt P1(P.X() - p1*D1.X(),
1168 P.Y() - p1*D1.Y(),
1169 P.Z() - p1*D1.Z());
ecc4f148 1170
7eed5d29 1171 P = Cyl2.Axis().Location();
ecc4f148 1172
7eed5d29 1173 gp_Pnt P2(P.X() - p2*D2.X(),
1174 P.Y() - p2*D2.Y(),
1175 P.Z() - p2*D2.Z());
ecc4f148 1176
7eed5d29 1177 gp_Vec P1P2(P1,P2);
1178 D1=gp_Dir(P1P2);
1179 p1=Cyl1.Radius();
ecc4f148 1180
7eed5d29 1181 pt1.SetCoord(P1.X() + p1*D1.X(),
1182 P1.Y() + p1*D1.Y(),
1183 P1.Z() + p1*D1.Z());
1184 nbint = 1;
7fd59977 1185 }
ecc4f148 1186 else
1187 {
7eed5d29 1188 typeres=IntAna_NoGeometricSolution;
7fd59977 1189 }
1190 }
1191 }
1192}
1193//=======================================================================
1194//function : IntAna_QuadQuadGeo
1195//purpose : Cylinder - Cone
1196//=======================================================================
a34f083b 1197IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1198 const gp_Cone& Con,
1199 const Standard_Real Tol)
7fd59977 1200: done(Standard_False),
1201 nbint(0),
1202 typeres(IntAna_Empty),
1203 pt1(0,0,0),
1204 pt2(0,0,0),
7eed5d29 1205 pt3(0,0,0),
1206 pt4(0,0,0),
7fd59977 1207 param1(0),
1208 param2(0),
7eed5d29 1209 param3(0),
1210 param4(0),
7fd59977 1211 param1bis(0),
1212 param2bis(0),
1213 myCommonGen(Standard_False),
1214 myPChar(0,0,0)
1215{
1216 InitTolerances();
1217 Perform(Cyl,Con,Tol);
1218}
1219//=======================================================================
1220//function : Perform
1221//purpose :
1222//=======================================================================
1223 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
7eed5d29 1224 const gp_Cone& Con,
1225 const Standard_Real )
7fd59977 1226{
1227 done=Standard_True;
1228 AxeOperator A1A2(Cyl.Axis(),Con.Axis());
1229 if(A1A2.Same()) {
1230 gp_Pnt Pt=Con.Apex();
1231 Standard_Real dist=Cyl.Radius()/(Tan(Con.SemiAngle()));
1232 gp_Dir dir=Cyl.Position().Direction();
1233 pt1.SetCoord( Pt.X() + dist*dir.X()
7eed5d29 1234 ,Pt.Y() + dist*dir.Y()
1235 ,Pt.Z() + dist*dir.Z());
7fd59977 1236 pt2.SetCoord( Pt.X() - dist*dir.X()
7eed5d29 1237 ,Pt.Y() - dist*dir.Y()
1238 ,Pt.Z() - dist*dir.Z());
7fd59977 1239 dir1=dir2=dir;
1240 param1=param2=Cyl.Radius();
1241 nbint=2;
1242 typeres=IntAna_Circle;
1243
1244 }
1245 else {
1246 typeres=IntAna_NoGeometricSolution;
1247 }
1248}
1249//=======================================================================
1250//function :
1251//purpose : Cylinder - Sphere
1252//=======================================================================
1253 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
7eed5d29 1254 const gp_Sphere& Sph,
1255 const Standard_Real Tol)
7fd59977 1256: done(Standard_False),
1257 nbint(0),
1258 typeres(IntAna_Empty),
1259 pt1(0,0,0),
1260 pt2(0,0,0),
7eed5d29 1261 pt3(0,0,0),
1262 pt4(0,0,0),
7fd59977 1263 param1(0),
1264 param2(0),
7eed5d29 1265 param3(0),
1266 param4(0),
7fd59977 1267 param1bis(0),
1268 param2bis(0),
1269 myCommonGen(Standard_False),
1270 myPChar(0,0,0)
1271{
1272 InitTolerances();
1273 Perform(Cyl,Sph,Tol);
1274}
1275//=======================================================================
1276//function : Perform
1277//purpose :
1278//=======================================================================
1279 void IntAna_QuadQuadGeo::Perform( const gp_Cylinder& Cyl
7eed5d29 1280 ,const gp_Sphere& Sph
1281 ,const Standard_Real)
7fd59977 1282{
1283 done=Standard_True;
1284 gp_Pnt Pt=Sph.Location();
1285 AxeOperator A1A2(Cyl.Axis(),Sph.Position().Axis());
1286 if((A1A2.Intersect() && Pt.Distance(A1A2.PtIntersect())==0.0 )
1287 || (A1A2.Same())) {
1288 if(Sph.Radius() < Cyl.Radius()) {
1289 typeres = IntAna_Empty;
1290 }
1291 else {
1292 Standard_Real dist=Sqrt( Sph.Radius() * Sph.Radius() - Cyl.Radius() * Cyl.Radius() );
1293 gp_Dir dir=Cyl.Position().Direction();
1294 dir1 = dir2 = dir;
1295 typeres=IntAna_Circle;
1296 pt1.SetCoord( Pt.X() + dist*dir.X()
7eed5d29 1297 ,Pt.Y() + dist*dir.Y()
1298 ,Pt.Z() + dist*dir.Z());
7fd59977 1299 nbint=1;
1300 param1 = Cyl.Radius();
1301 if(dist>RealEpsilon()) {
7eed5d29 1302 pt2.SetCoord( Pt.X() - dist*dir.X()
1303 ,Pt.Y() - dist*dir.Y()
1304 ,Pt.Z() - dist*dir.Z());
1305 param2=Cyl.Radius();
1306 nbint=2;
7fd59977 1307 }
1308 }
1309 }
1310 else {
1311 typeres=IntAna_NoGeometricSolution;
1312 }
1313}
1314
1315//=======================================================================
1316//function : IntAna_QuadQuadGeo
1317//purpose : Cone - Cone
1318//=======================================================================
1319 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con1,
7eed5d29 1320 const gp_Cone& Con2,
1321 const Standard_Real Tol)
7fd59977 1322: done(Standard_False),
1323 nbint(0),
1324 typeres(IntAna_Empty),
1325 pt1(0,0,0),
1326 pt2(0,0,0),
7eed5d29 1327 pt3(0,0,0),
1328 pt4(0,0,0),
7fd59977 1329 param1(0),
1330 param2(0),
7eed5d29 1331 param3(0),
1332 param4(0),
7fd59977 1333 param1bis(0),
1334 param2bis(0),
1335 myCommonGen(Standard_False),
1336 myPChar(0,0,0)
1337{
1338 InitTolerances();
1339 Perform(Con1,Con2,Tol);
1340}
1341//
1342//=======================================================================
1343//function : Perform
1344//purpose :
1345//=======================================================================
1346 void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con1,
7eed5d29 1347 const gp_Cone& Con2,
1348 const Standard_Real Tol)
7fd59977 1349{
1350 done=Standard_True;
1351 //
1352 Standard_Real tg1, tg2, aDA1A2, aTol2;
1353 gp_Pnt aPApex1, aPApex2;
4bd102b8 1354
1355 Standard_Real TOL_APEX_CONF = 1.e-10;
1356
7fd59977 1357 //
1358 tg1=Tan(Con1.SemiAngle());
1359 tg2=Tan(Con2.SemiAngle());
1360
1361 if((tg1 * tg2) < 0.) {
1362 tg2 = -tg2;
1363 }
1364 //
7fd59977 1365 aTol2=Tol*Tol;
1366 aPApex1=Con1.Apex();
1367 aPApex2=Con2.Apex();
1368 aDA1A2=aPApex1.SquareDistance(aPApex2);
7fd59977 1369 //
1370 AxeOperator A1A2(Con1.Axis(),Con2.Axis());
1371 //
1372 // 1
1373 if(A1A2.Same()) {
1374 //-- two circles
1375 Standard_Real x;
1376 gp_Pnt P=Con1.Apex();
1377 gp_Dir D=Con1.Position().Direction();
1378 Standard_Real d=gp_Vec(D).Dot(gp_Vec(P,Con2.Apex()));
1379
1380 if(Abs(tg1-tg2)>myEPSILON_ANGLE_CONE) {
4bd102b8 1381 if (fabs(d) < TOL_APEX_CONF) {
7eed5d29 1382 typeres = IntAna_Point;
1383 nbint = 1;
1384 pt1 = P;
1385 return;
4bd102b8 1386 }
7fd59977 1387 x=(d*tg2)/(tg1+tg2);
1388 pt1.SetCoord( P.X() + x*D.X()
7eed5d29 1389 ,P.Y() + x*D.Y()
1390 ,P.Z() + x*D.Z());
7fd59977 1391 param1=Abs(x*tg1);
1392
1393 x=(d*tg2)/(tg2-tg1);
1394 pt2.SetCoord( P.X() + x*D.X()
7eed5d29 1395 ,P.Y() + x*D.Y()
1396 ,P.Z() + x*D.Z());
7fd59977 1397 param2=Abs(x*tg1);
1398 dir1 = dir2 = D;
1399 nbint=2;
1400 typeres=IntAna_Circle;
1401 }
1402 else {
4bd102b8 1403 if (fabs(d) < TOL_APEX_CONF) {
7eed5d29 1404 typeres=IntAna_Same;
7fd59977 1405 }
1406 else {
7eed5d29 1407 typeres=IntAna_Circle;
1408 nbint=1;
1409 x=d*0.5;
1410 pt1.SetCoord( P.X() + x*D.X()
1411 ,P.Y() + x*D.Y()
1412 ,P.Z() + x*D.Z());
1413 param1 = Abs(x * tg1);
1414 dir1 = D;
7fd59977 1415 }
1416 }
1417 } //-- fin A1A2.Same
1418 // 2
1419 else if((Abs(tg1-tg2)<myEPSILON_ANGLE_CONE) && (A1A2.Parallel())) {
1420 //-- voir AnVer12mai98
1421 Standard_Real DistA1A2=A1A2.Distance();
1422 gp_Dir DA1=Con1.Position().Direction();
1423 gp_Vec O1O2(Con1.Apex(),Con2.Apex());
b045e6a4 1424 gp_Dir O1O2n(O1O2); // normalization of the vector before projection
1425 Standard_Real O1O2_DA1=gp_Vec(DA1).Dot(gp_Vec(O1O2n));
1426
1427 gp_Vec O1_Proj_A2(O1O2n.X()-O1O2_DA1*DA1.X(),
7eed5d29 1428 O1O2n.Y()-O1O2_DA1*DA1.Y(),
1429 O1O2n.Z()-O1O2_DA1*DA1.Z());
7fd59977 1430 gp_Dir DB1=gp_Dir(O1_Proj_A2);
b045e6a4 1431
7fd59977 1432 Standard_Real yO1O2=O1O2.Dot(gp_Vec(DA1));
1433 Standard_Real ABSTG1 = Abs(tg1);
1434 Standard_Real X2 = (DistA1A2/ABSTG1 - yO1O2)*0.5;
1435 Standard_Real X1 = X2+yO1O2;
1436
1437 gp_Pnt P1(Con1.Apex().X() + X1*( DA1.X() + ABSTG1*DB1.X()),
7eed5d29 1438 Con1.Apex().Y() + X1*( DA1.Y() + ABSTG1*DB1.Y()),
1439 Con1.Apex().Z() + X1*( DA1.Z() + ABSTG1*DB1.Z()));
7fd59977 1440
1441 gp_Pnt MO1O2(0.5*(Con1.Apex().X()+Con2.Apex().X()),
7eed5d29 1442 0.5*(Con1.Apex().Y()+Con2.Apex().Y()),
1443 0.5*(Con1.Apex().Z()+Con2.Apex().Z()));
7fd59977 1444 gp_Vec P1MO1O2(P1,MO1O2);
1445
1446 gp_Dir DA1_X_DB1=DA1.Crossed(DB1);
1447 gp_Dir OrthoPln = DA1_X_DB1.Crossed(gp_Dir(P1MO1O2));
1448
1449 IntAna_QuadQuadGeo INTER_QUAD_PLN(gp_Pln(P1,OrthoPln),Con1,Tol,Tol);
1450 if(INTER_QUAD_PLN.IsDone()) {
1451 switch(INTER_QUAD_PLN.TypeInter()) {
7eed5d29 1452 case IntAna_Ellipse: {
1453 typeres=IntAna_Ellipse;
1454 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1455 pt1 = E.Location();
1456 dir1 = E.Position().Direction();
1457 dir2 = E.Position().XDirection();
1458 param1 = E.MajorRadius();
1459 param1bis = E.MinorRadius();
1460 nbint = 1;
1461 break;
7fd59977 1462 }
1463 case IntAna_Circle: {
7eed5d29 1464 typeres=IntAna_Circle;
1465 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1466 pt1 = C.Location();
1467 dir1 = C.Position().XDirection();
1468 dir2 = C.Position().YDirection();
1469 param1 = C.Radius();
1470 nbint = 1;
1471 break;
7fd59977 1472 }
1473 case IntAna_Hyperbola: {
7eed5d29 1474 typeres=IntAna_Hyperbola;
1475 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1476 pt1 = pt2 = H.Location();
1477 dir1 = H.Position().Direction();
1478 dir2 = H.Position().XDirection();
1479 param1 = param2 = H.MajorRadius();
1480 param1bis = param2bis = H.MinorRadius();
1481 nbint = 2;
1482 break;
7fd59977 1483 }
1484 case IntAna_Line: {
7eed5d29 1485 typeres=IntAna_Line;
1486 gp_Lin H=INTER_QUAD_PLN.Line(1);
1487 pt1 = pt2 = H.Location();
1488 dir1 = dir2 = H.Position().Direction();
1489 param1 = param2 = 0.0;
1490 param1bis = param2bis = 0.0;
1491 nbint = 2;
1492 break;
7fd59977 1493 }
1494 default:
7eed5d29 1495 typeres=IntAna_NoGeometricSolution;
7fd59977 1496 }
1497 }
1498 }// else if((Abs(tg1-tg2)<EPSILON_ANGLE_CONE) && (A1A2.Parallel()))
7fd59977 1499 // 3
1500 else if (aDA1A2<aTol2) {
7fd59977 1501 //
1502 // When apices are coinsided there can be 3 possible cases
1503 // 3.1 - empty solution (iRet=0)
1504 // 3.2 - one line when cone1 touches cone2 (iRet=1)
1505 // 3.3 - two lines when cone1 intersects cone2 (iRet=2)
1506 //
1507 Standard_Integer iRet;
1508 Standard_Real aGamma, aBeta1, aBeta2;
1509 Standard_Real aD1, aR1, aTgBeta1, aTgBeta2, aHalfPI;
1510 Standard_Real aCosGamma, aSinGamma, aDx, aR2, aRD2, aD2;
1511 gp_Pnt2d aP0, aPA1, aP1, aPA2;
1512 gp_Vec2d aVAx2;
1513 gp_Ax1 aAx1, aAx2;
1514 //
1515 // Preliminary analysis. Determination of iRet
1516 //
1517 iRet=0;
c6541a0c 1518 aHalfPI=0.5*M_PI;
7fd59977 1519 aD1=1.;
1520 aPA1.SetCoord(aD1, 0.);
1521 aP0.SetCoord(0., 0.);
1522 //
1523 aAx1=Con1.Axis();
1524 aAx2=Con2.Axis();
1525 aGamma=aAx1.Angle(aAx2);
1526 if (aGamma>aHalfPI){
c6541a0c 1527 aGamma=M_PI-aGamma;
7fd59977 1528 }
1529 aCosGamma=Cos(aGamma);
1530 aSinGamma=Sin(aGamma);
1531 //
1532 aBeta1=Con1.SemiAngle();
1533 aTgBeta1=Tan(aBeta1);
1534 aTgBeta1=Abs(aTgBeta1);
1535 //
1536 aBeta2=Con2.SemiAngle();
1537 aTgBeta2=Tan(aBeta2);
1538 aTgBeta2=Abs(aTgBeta2);
1539 //
1540 aR1=aD1*aTgBeta1;
1541 aP1.SetCoord(aD1, aR1);
1542 //
1543 // PA2
1544 aVAx2.SetCoord(aCosGamma, aSinGamma);
1545 gp_Dir2d aDAx2(aVAx2);
1546 gp_Lin2d aLAx2(aP0, aDAx2);
1547 //
1548 gp_Vec2d aV(aP0, aP1);
1549 aDx=aVAx2.Dot(aV);
1550 aPA2=aP0.Translated(aDx*aDAx2);
1551 //
1552 // aR2
1553 aDx=aPA2.Distance(aP0);
1554 aR2=aDx*aTgBeta2;
1555 //
1556 // aRD2
1557 aRD2=aPA2.Distance(aP1);
1558 //
1559 if (aRD2>(aR2+Tol)) {
1560 iRet=0;
7fd59977 1561 typeres=IntAna_Empty; //nothing
4101383e 1562 return;
7fd59977 1563 }
1564 //
1565 iRet=1; //touch case => 1 line
1566 if (aRD2<(aR2-Tol)) {
1567 iRet=2;//intersection => couple of lines
1568 }
1569 //
1570 // Finding the solution in 3D
1571 //
1572 Standard_Real aDa;
1573 gp_Pnt aQApex1, aQA1, aQA2, aQX, aQX1, aQX2;
1574 gp_Dir aD3Ax1, aD3Ax2;
1575 gp_Lin aLin;
1576 IntAna_QuadQuadGeo aIntr;
1577 //
1578 aQApex1=Con1.Apex();
1579 aD3Ax1=aAx1.Direction();
1580 aQA1.SetCoord(aQApex1.X()+aD1*aD3Ax1.X(),
7eed5d29 1581 aQApex1.Y()+aD1*aD3Ax1.Y(),
1582 aQApex1.Z()+aD1*aD3Ax1.Z());
7fd59977 1583 //
1584 aDx=aD3Ax1.Dot(aAx2.Direction());
1585 if (aDx<0.) {
1586 aAx2.Reverse();
1587 }
1588 aD3Ax2=aAx2.Direction();
1589 //
1590 aD2=aD1*sqrt((1.+aTgBeta1*aTgBeta1)/(1.+aTgBeta2*aTgBeta2));
1591 //
1592 aQA2.SetCoord(aQApex1.X()+aD2*aD3Ax2.X(),
7eed5d29 1593 aQApex1.Y()+aD2*aD3Ax2.Y(),
1594 aQApex1.Z()+aD2*aD3Ax2.Z());
7fd59977 1595 //
1596 gp_Pln aPln1(aQA1, aD3Ax1);
1597 gp_Pln aPln2(aQA2, aD3Ax2);
1598 //
1599 aIntr.Perform(aPln1, aPln2, Tol, Tol);
a060129f 1600 if (!aIntr.IsDone() || 0 == aIntr.NbSolutions()) {
7fd59977 1601 iRet=-1; // just in case. it must not be so
1602 typeres=IntAna_NoGeometricSolution;
1603 return;
1604 }
1605 //
1606 aLin=aIntr.Line(1);
1607 const gp_Dir& aDLin=aLin.Direction();
1608 gp_Vec aVLin(aDLin);
1609 gp_Pnt aOrig=aLin.Location();
1610 gp_Vec aVr(aQA1, aOrig);
1611 aDx=aVLin.Dot(aVr);
1612 aQX=aOrig.Translated(aDx*aVLin);
1613 //
1614 // Final part
1615 //
1616 typeres=IntAna_Line;
1617 //
1618 param1=0.;
1619 param2 =0.;
1620 param1bis=0.;
1621 param2bis=0.;
1622 //
1623 if (iRet==1) {
1624 // one line
1625 nbint=1;
1626 pt1=aQApex1;
1627 gp_Vec aVX(aQApex1, aQX);
1628 dir1=gp_Dir(aVX);
7fd59977 1629 }
1630
1631 else {//iRet=2
1632 // two lines
1633 nbint=2;
1634 aDa=aQA1.Distance(aQX);
1635 aDx=sqrt(aR1*aR1-aDa*aDa);
1636 aQX1=aQX.Translated(aDx*aVLin);
1637 aQX2=aQX.Translated(-aDx*aVLin);
1638 //
1639 pt1=aQApex1;
1640 pt2=aQApex1;
1641 gp_Vec aVX1(aQApex1, aQX1);
1642 dir1=gp_Dir(aVX1);
1643 gp_Vec aVX2(aQApex1, aQX2);
1644 dir2=gp_Dir(aVX2);
7fd59977 1645 }
1646 } //else if (aDA1A2<aTol2) {
7fd59977 1647 //Case when cones have common generatrix
1648 else if(A1A2.Intersect()) {
1649 //Check if apex of one cone belongs another one
1650 Standard_Real u, v, tol2 = Tol*Tol;
1651 ElSLib::Parameters(Con2, aPApex1, u, v);
1652 gp_Pnt p = ElSLib::Value(u, v, Con2);
1653 if(aPApex1.SquareDistance(p) > tol2) {
1654 typeres=IntAna_NoGeometricSolution;
1655 return;
1656 }
1657 //
1658 ElSLib::Parameters(Con1, aPApex2, u, v);
1659 p = ElSLib::Value(u, v, Con1);
1660 if(aPApex2.SquareDistance(p) > tol2) {
1661 typeres=IntAna_NoGeometricSolution;
1662 return;
1663 }
1664
1665 //Cones have a common generatrix passing through apexes
1666 myCommonGen = Standard_True;
1667
1668 //common generatrix of cones
1669 gp_Lin aGen(aPApex1, gp_Dir(gp_Vec(aPApex1, aPApex2)));
1670
1671 //Intersection point of axes
1672 gp_Pnt aPAxeInt = A1A2.PtIntersect();
1673
1674 //Characteristic point of intersection curve
1675 u = ElCLib::Parameter(aGen, aPAxeInt);
1676 myPChar = ElCLib::Value(u, aGen);
1677
1678
1679 //Other generatrixes of cones laying in maximal plane
c6541a0c
D		 1680 gp_Lin aGen1 = aGen.Rotated(Con1.Axis(), M_PI);
1681 gp_Lin aGen2 = aGen.Rotated(Con2.Axis(), M_PI);
7fd59977 1682 //
1683 //Intersection point of generatrixes
1684 gp_Dir aN; //solution plane normal
1685 gp_Dir aD1 = aGen1.Direction();
1686
1687 gp_Dir aD2(aD1.Crossed(aGen.Direction()));
1688
1689 if(aD1.IsParallel(aGen2.Direction(), Precision::Angular())) {
1690 aN = aD1.Crossed(aD2);
1691 }
1692 else if(aGen1.SquareDistance(aGen2) > tol2) {
1693 //Something wrong ???
1694 typeres=IntAna_NoGeometricSolution;
1695 return;
1696 }
1697 else {
1698 gp_Dir D1 = aGen1.Position().Direction();
1699 gp_Dir D2 = aGen2.Position().Direction();
1700 gp_Pnt O1 = aGen1.Location();
1701 gp_Pnt O2 = aGen2.Location();
1702 Standard_Real D1DotD2 = D1.Dot(D2);
1703 Standard_Real aSin = 1.-D1DotD2*D1DotD2;
1704 gp_Vec O1O2 (O1,O2);
1705 Standard_Real U2 = (D1.XYZ()*(O1O2.Dot(D1))-(O1O2.XYZ())).Dot(D2.XYZ());
1706 U2 /= aSin;
1707 gp_Pnt aPGint(ElCLib::Value(U2, aGen2));
1708
1709 aD1 = gp_Dir(gp_Vec(aPGint, myPChar));
1710 aN = aD1.Crossed(aD2);
1711 }
1712 //Plane that must contain intersection curves
1713 gp_Pln anIntPln(myPChar, aN);
1714
1715 IntAna_QuadQuadGeo INTER_QUAD_PLN(anIntPln,Con1,Tol,Tol);
1716
1717 if(INTER_QUAD_PLN.IsDone()) {
1718 switch(INTER_QUAD_PLN.TypeInter()) {
7eed5d29 1719 case IntAna_Ellipse: {
1720 typeres=IntAna_Ellipse;
1721 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1722 pt1 = E.Location();
1723 dir1 = E.Position().Direction();
1724 dir2 = E.Position().XDirection();
1725 param1 = E.MajorRadius();
1726 param1bis = E.MinorRadius();
1727 nbint = 1;
1728 break;
7fd59977 1729 }
1730 case IntAna_Circle: {
7eed5d29 1731 typeres=IntAna_Circle;
1732 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1733 pt1 = C.Location();
1734 dir1 = C.Position().XDirection();
1735 dir2 = C.Position().YDirection();
1736 param1 = C.Radius();
1737 nbint = 1;
1738 break;
7fd59977 1739 }
1740 case IntAna_Parabola: {
7eed5d29 1741 typeres=IntAna_Parabola;
1742 gp_Parab Prb=INTER_QUAD_PLN.Parabola(1);
1743 pt1 = Prb.Location();
1744 dir1 = Prb.Position().Direction();
1745 dir2 = Prb.Position().XDirection();
1746 param1 = Prb.Focal();
1747 nbint = 1;
1748 break;
7fd59977 1749 }
1750 case IntAna_Hyperbola: {
7eed5d29 1751 typeres=IntAna_Hyperbola;
1752 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1753 pt1 = pt2 = H.Location();
1754 dir1 = H.Position().Direction();
1755 dir2 = H.Position().XDirection();
1756 param1 = param2 = H.MajorRadius();
1757 param1bis = param2bis = H.MinorRadius();
1758 nbint = 2;
1759 break;
7fd59977 1760 }
1761 default:
7eed5d29 1762 typeres=IntAna_NoGeometricSolution;
7fd59977 1763 }
1764 }
1765 }
4101383e 1766
7fd59977 1767 else {
1768 typeres=IntAna_NoGeometricSolution;
1769 }
1770}
1771//=======================================================================
1772//function : IntAna_QuadQuadGeo
1773//purpose : Sphere - Cone
1774//=======================================================================
1775 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
7eed5d29 1776 const gp_Cone& Con,
1777 const Standard_Real Tol)
7fd59977 1778: done(Standard_False),
1779 nbint(0),
1780 typeres(IntAna_Empty),
1781 pt1(0,0,0),
1782 pt2(0,0,0),
7eed5d29 1783 pt3(0,0,0),
1784 pt4(0,0,0),
7fd59977 1785 param1(0),
1786 param2(0),
7eed5d29 1787 param3(0),
1788 param4(0),
7fd59977 1789 param1bis(0),
1790 param2bis(0),
1791 myCommonGen(Standard_False),
1792 myPChar(0,0,0)
1793{
1794 InitTolerances();
1795 Perform(Sph,Con,Tol);
1796}
1797//=======================================================================
1798//function : Perform
1799//purpose :
1800//=======================================================================
1801 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
7eed5d29 1802 const gp_Cone& Con,
1803 const Standard_Real)
7fd59977 1804{
77088633 1805
1806 //
7fd59977 1807 done=Standard_True;
77088633 1808 //
7fd59977 1809 AxeOperator A1A2(Con.Axis(),Sph.Position().Axis());
1810 gp_Pnt Pt=Sph.Location();
77088633 1811 //
7fd59977 1812 if((A1A2.Intersect() && (Pt.Distance(A1A2.PtIntersect())==0.0))
1813 || A1A2.Same()) {
1814 gp_Pnt ConApex= Con.Apex();
1815 Standard_Real dApexSphCenter=Pt.Distance(ConApex);
1816 gp_Dir ConDir;
1817 if(dApexSphCenter>RealEpsilon()) {
1818 ConDir = gp_Dir(gp_Vec(ConApex,Pt));
1819 }
1820 else {
1821 ConDir = Con.Position().Direction();
1822 }
1823
1824 Standard_Real Rad=Sph.Radius();
1825 Standard_Real tga=Tan(Con.SemiAngle());
1826
1827
1828 //-- 2 circles
1829 //-- x: Roots of (x**2 + y**2 = Rad**2)
1830 //-- tga = y / (x+dApexSphCenter)
1831 Standard_Real tgatga = tga * tga;
1832 math_DirectPolynomialRoots Eq( 1.0+tgatga
7eed5d29 1833 ,2.0*tgatga*dApexSphCenter
1834 ,-Rad*Rad + dApexSphCenter*dApexSphCenter*tgatga);
7fd59977 1835 if(Eq.IsDone()) {
1836 Standard_Integer nbsol=Eq.NbSolutions();
1837 if(nbsol==0) {
7eed5d29 1838 typeres=IntAna_Empty;
7fd59977 1839 }
1840 else {
7eed5d29 1841 typeres=IntAna_Circle;
1842 if(nbsol>=1) {
1843 Standard_Real x = Eq.Value(1);
1844 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1845 nbint=1;
1846 pt1.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1847 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1848 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1849 param1 = tga * dApexSphCenterpx;
1850 param1 = Abs(param1);
1851 dir1 = ConDir;
1852 if(param1<=myEPSILON_MINI_CIRCLE_RADIUS) {
1853 typeres=IntAna_PointAndCircle;
1854 param1=0.0;
1855 }
1856 }
1857 if(nbsol>=2) {
1858 Standard_Real x=Eq.Value(2);
1859 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1860 nbint=2;
1861 pt2.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1862 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1863 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1864 param2 = tga * dApexSphCenterpx;
1865 param2 = Abs(param2);
1866 dir2=ConDir;
1867 if(param2<=myEPSILON_MINI_CIRCLE_RADIUS) {
1868 typeres=IntAna_PointAndCircle;
1869 param2=0.0;
1870 }
1871 }
7fd59977 1872 }
1873 }
1874 else {
1875 done=Standard_False;
1876 }
1877 }
1878 else {
1879 typeres=IntAna_NoGeometricSolution;
1880 }
1881}
1882
1883//=======================================================================
1884//function : IntAna_QuadQuadGeo
1885//purpose : Sphere - Sphere
1886//=======================================================================
1887 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Sphere& Sph1
7eed5d29 1888 ,const gp_Sphere& Sph2
1889 ,const Standard_Real Tol)
7fd59977 1890: done(Standard_False),
1891 nbint(0),
1892 typeres(IntAna_Empty),
1893 pt1(0,0,0),
1894 pt2(0,0,0),
7eed5d29 1895 pt3(0,0,0),
1896 pt4(0,0,0),
7fd59977 1897 param1(0),
1898 param2(0),
7eed5d29 1899 param3(0),
1900 param4(0),
7fd59977 1901 param1bis(0),
1902 param2bis(0),
1903 myCommonGen(Standard_False),
1904 myPChar(0,0,0)
1905{
1906 InitTolerances();
1907 Perform(Sph1,Sph2,Tol);
1908}
1909//=======================================================================
1910//function : Perform
1911//purpose :
1912//=======================================================================
1913 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph1,
7eed5d29 1914 const gp_Sphere& Sph2,
1915 const Standard_Real Tol)
7fd59977 1916{
1917 done=Standard_True;
1918 gp_Pnt O1=Sph1.Location();
1919 gp_Pnt O2=Sph2.Location();
1920 Standard_Real dO1O2=O1.Distance(O2);
1921 Standard_Real R1=Sph1.Radius();
1922 Standard_Real R2=Sph2.Radius();
1923 Standard_Real Rmin,Rmax;
1924 typeres=IntAna_Empty;
1925 param2bis=0.0; //-- pour eviter param2bis not used ....
1926
1927 if(R1>R2) { Rmin=R2; Rmax=R1; } else { Rmin=R1; Rmax=R2; }
1928
1929 if(dO1O2<=Tol && (Abs(R1-R2) <= Tol)) {
1930 typeres = IntAna_Same;
1931 }
1932 else {
1933 if(dO1O2<=Tol) { return; }
1934 gp_Dir Dir=gp_Dir(gp_Vec(O1,O2));
1935 Standard_Real t = Rmax - dO1O2 - Rmin;
1936
1937 //----------------------------------------------------------------------
1938 //-- |----------------- R1 --------------------|
1939 //-- |----dO1O2-----|-----------R2----------|
1940 //-- --->--<-- t
1941 //--
1942 //-- |------ R1 ------|---------dO1O2----------|
1943 //-- |-------------------R2-----------------------|
1944 //-- --->--<-- t
1945 //----------------------------------------------------------------------
1946 if(t >= 0.0 && t <=Tol) {
1947 typeres = IntAna_Point;
1948 nbint = 1;
1949 Standard_Real t2;
1950 if(R1==Rmax) t2=(R1 + (R2 + dO1O2)) * 0.5;
1951 else t2=(-R1+(dO1O2-R2))*0.5;
7eed5d29 1952
7fd59977 1953 pt1.SetCoord( O1.X() + t2*Dir.X()
7eed5d29 1954 ,O1.Y() + t2*Dir.Y()
1955 ,O1.Z() + t2*Dir.Z());
7fd59977 1956 }
1957 else {
1958 //-----------------------------------------------------------------
1959 //-- |----------------- dO1O2 --------------------|
1960 //-- |----R1-----|-----------R2----------|-Tol-|
1961 //--
1962 //-- |----------------- Rmax --------------------|
1963 //-- |----Rmin----|-------dO1O2-------|-Tol-|
1964 //--
1965 //-----------------------------------------------------------------
1966 if((dO1O2 > (R1+R2+Tol)) || (Rmax > (dO1O2+Rmin+Tol))) {
7eed5d29 1967 typeres=IntAna_Empty;
7fd59977 1968 }
1969 else {
7eed5d29 1970 //---------------------------------------------------------------
1971 //--
1972 //--
1973 //---------------------------------------------------------------
1974 Standard_Real Alpha=0.5*(R1*R1-R2*R2+dO1O2*dO1O2)/(dO1O2);
1975 Standard_Real Beta = R1*R1-Alpha*Alpha;
1976 Beta = (Beta>0.0)? Sqrt(Beta) : 0.0;
1977
1978 if(Beta<= myEPSILON_MINI_CIRCLE_RADIUS) {
1979 typeres = IntAna_Point;
1980 Alpha = (R1 + (dO1O2 - R2)) * 0.5;
1981 }
1982 else {
1983 typeres = IntAna_Circle;
1984 dir1 = Dir;
1985 param1 = Beta;
1986 }
1987 pt1.SetCoord( O1.X() + Alpha*Dir.X()
1988 ,O1.Y() + Alpha*Dir.Y()
1989 ,O1.Z() + Alpha*Dir.Z());
1990
1991 nbint=1;
7fd59977 1992 }
1993 }
1994 }
1995}
7eed5d29 1996
1997//=======================================================================
1998//function : IntAna_QuadQuadGeo
1999//purpose : Plane - Torus
2000//=======================================================================
2001IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& Pln,
2002 const gp_Torus& Tor,
2003 const Standard_Real Tol)
2004: done(Standard_False),
2005 nbint(0),
2006 typeres(IntAna_Empty),
2007 pt1(0,0,0),
2008 pt2(0,0,0),
2009 pt3(0,0,0),
2010 pt4(0,0,0),
2011 param1(0),
2012 param2(0),
2013 param3(0),
2014 param4(0),
2015 param1bis(0),
2016 param2bis(0),
2017 myCommonGen(Standard_False),
2018 myPChar(0,0,0)
2019{
2020 InitTolerances();
2021 Perform(Pln,Tor,Tol);
2022}
2023//=======================================================================
2024//function : Perform
2025//purpose :
2026//=======================================================================
2027void IntAna_QuadQuadGeo::Perform(const gp_Pln& Pln,
2028 const gp_Torus& Tor,
2029 const Standard_Real Tol)
2030{
2031 done = Standard_True;
2032 //
2033 Standard_Real aRMin, aRMaj;
2034 //
2035 aRMin = Tor.MinorRadius();
2036 aRMaj = Tor.MajorRadius();
2037 if (aRMin >= aRMaj) {
2038 typeres = IntAna_NoGeometricSolution;
2039 return;
2040 }
2041 //
2042 const gp_Ax1 aPlnAx = Pln.Axis();
2043 const gp_Ax1 aTorAx = Tor.Axis();
2044 //
2045 Standard_Boolean bParallel, bNormal;
2046 //
2047 bParallel = aTorAx.IsParallel(aPlnAx, myEPSILON_AXES_PARA);
2048 bNormal = !bParallel ? aTorAx.IsNormal(aPlnAx, myEPSILON_AXES_PARA) : Standard_False;
2049 if (!bNormal && !bParallel) {
2050 typeres = IntAna_NoGeometricSolution;
2051 return;
2052 }
2053 //
2054 Standard_Real aDist;
2055 //
2056 gp_Pnt aTorLoc = aTorAx.Location();
2057 if (bParallel) {
2058 Standard_Real aDt, X, Y, Z, A, B, C, D;
2059 //
2060 Pln.Coefficients(A,B,C,D);
2061 aTorLoc.Coord(X,Y,Z);
2062 aDist = A*X + B*Y + C*Z + D;
2063 //
2064 if ((Abs(aDist) - aRMin) > Tol) {
2065 typeres=IntAna_Empty;
2066 return;
2067 }
2068 //
2069 typeres = IntAna_Circle;
2070 //
2071 pt1.SetCoord(X - aDist*A, Y - aDist*B, Z - aDist*C);
2072 aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
2073 param1 = aRMaj + aDt;
2074 dir1 = aTorAx.Direction();
2075 nbint = 1;
2076 if ((Abs(aDist) < aRMin) && (aDt > Tol)) {
2077 pt2 = pt1;
2078 param2 = aRMaj - aDt;
2079 dir2 = dir1;
2080 nbint = 2;
2081 }
2082 }
2083 //
2084 else {
2085 aDist = Pln.Distance(aTorLoc);
2086 if (aDist > myEPSILON_DISTANCE) {
2087 typeres = IntAna_NoGeometricSolution;
2088 return;
2089 }
2090 //
2091 typeres = IntAna_Circle;
2092 param2 = param1 = aRMin;
2093 dir2 = dir1 = aPlnAx.Direction();
2094 nbint = 2;
2095 //
2096 gp_Dir aDir = aTorAx.Direction()^dir1;
2097 pt1.SetXYZ(aTorLoc.XYZ() + aRMaj*aDir.XYZ());
2098 pt2.SetXYZ(aTorLoc.XYZ() - aRMaj*aDir.XYZ());
2099 }
2100}
2101
2102//=======================================================================
2103//function : IntAna_QuadQuadGeo
2104//purpose : Cylinder - Torus
2105//=======================================================================
2106IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
2107 const gp_Torus& Tor,
2108 const Standard_Real Tol)
2109: done(Standard_False),
2110 nbint(0),
2111 typeres(IntAna_Empty),
2112 pt1(0,0,0),
2113 pt2(0,0,0),
2114 pt3(0,0,0),
2115 pt4(0,0,0),
2116 param1(0),
2117 param2(0),
2118 param3(0),
2119 param4(0),
2120 param1bis(0),
2121 param2bis(0),
2122 myCommonGen(Standard_False),
2123 myPChar(0,0,0)
2124{
2125 InitTolerances();
2126 Perform(Cyl,Tor,Tol);
2127}
2128//=======================================================================
2129//function : Perform
2130//purpose :
2131//=======================================================================
2132void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
2133 const gp_Torus& Tor,
2134 const Standard_Real Tol)
2135{
2136 done = Standard_True;
2137 //
2138 Standard_Real aRMin, aRMaj;
2139 //
2140 aRMin = Tor.MinorRadius();
2141 aRMaj = Tor.MajorRadius();
2142 if (aRMin >= aRMaj) {
2143 typeres = IntAna_NoGeometricSolution;
2144 return;
2145 }
2146 //
2147 const gp_Ax1 aCylAx = Cyl.Axis();
2148 const gp_Ax1 aTorAx = Tor.Axis();
2149 //
2150 const gp_Lin aLin(aTorAx);
2151 const gp_Pnt aLocCyl = Cyl.Location();
2152 //
2153 if (!aTorAx.IsParallel(aCylAx, myEPSILON_AXES_PARA) ||
2154 (aLin.Distance(aLocCyl) > myEPSILON_DISTANCE)) {
2155 typeres = IntAna_NoGeometricSolution;
2156 return;
2157 }
2158 //
2159 Standard_Real aRCyl;
2160 //
2161 aRCyl = Cyl.Radius();
2162 if (((aRCyl + Tol) < (aRMaj - aRMin)) || ((aRCyl - Tol) > (aRMaj + aRMin))) {
2163 typeres = IntAna_Empty;
2164 return;
2165 }
2166 //
2167 typeres = IntAna_Circle;
2168 //
2169 Standard_Real aDist = Sqrt(Abs(aRMin*aRMin - (aRCyl-aRMaj)*(aRCyl-aRMaj)));
2170 gp_XYZ aTorLoc = aTorAx.Location().XYZ();
2171 //
2172 dir1 = aTorAx.Direction();
2173 pt1.SetXYZ(aTorLoc + aDist*dir1.XYZ());
2174 param1 = aRCyl;
2175 nbint = 1;
2176 if ((aDist > Tol) && (aRCyl > (aRMaj - aRMin)) &&
2177 (aRCyl < (aRMaj + aRMin))) {
2178 dir2 = dir1;
2179 pt2.SetXYZ(aTorLoc - aDist*dir2.XYZ());
2180 param2 = param1;
2181 nbint = 2;
2182 }
2183}
2184
2185//=======================================================================
2186//function : IntAna_QuadQuadGeo
2187//purpose : Cone - Torus
2188//=======================================================================
2189IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con,
2190 const gp_Torus& Tor,
2191 const Standard_Real Tol)
2192: done(Standard_False),
2193 nbint(0),
2194 typeres(IntAna_Empty),
2195 pt1(0,0,0),
2196 pt2(0,0,0),
2197 pt3(0,0,0),
2198 pt4(0,0,0),
2199 param1(0),
2200 param2(0),
2201 param3(0),
2202 param4(0),
2203 param1bis(0),
2204 param2bis(0),
2205 myCommonGen(Standard_False),
2206 myPChar(0,0,0)
2207{
2208 InitTolerances();
2209 Perform(Con,Tor,Tol);
2210}
2211//=======================================================================
2212//function : Perform
2213//purpose :
2214//=======================================================================
2215void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con,
2216 const gp_Torus& Tor,
2217 const Standard_Real Tol)
2218{
2219 done = Standard_True;
2220 //
2221 Standard_Real aRMin, aRMaj;
2222 //
2223 aRMin = Tor.MinorRadius();
2224 aRMaj = Tor.MajorRadius();
2225 if (aRMin >= aRMaj) {
2226 typeres = IntAna_NoGeometricSolution;
2227 return;
2228 }
2229 //
2230 const gp_Ax1 aConAx = Con.Axis();
2231 const gp_Ax1 aTorAx = Tor.Axis();
2232 //
2233 const gp_Lin aLin(aTorAx);
2234 const gp_Pnt aConApex = Con.Apex();
2235 //
2236 if (!aTorAx.IsParallel(aConAx, myEPSILON_AXES_PARA) ||
2237 (aLin.Distance(aConApex) > myEPSILON_DISTANCE)) {
2238 typeres = IntAna_NoGeometricSolution;
2239 return;
2240 }
2241 //
6092c0c8 2242 Standard_Real anAngle, aDist, aParam[4], aDt;
7eed5d29 2243 Standard_Integer i;
2244 gp_Pnt aTorLoc, aPCT, aPN, aPt[4];
2245 gp_Dir aDir[4];
2246 //
2247 anAngle = Con.SemiAngle();
2248 aTorLoc = aTorAx.Location();
2249 //
2250 aPN.SetXYZ(aTorLoc.XYZ() + aRMaj*Tor.YAxis().Direction().XYZ());
2251 gp_Dir aDN (gp_Vec(aTorLoc, aPN));
2252 gp_Ax1 anAxCLRot(aConApex, aDN);
2253 gp_Lin aConL = aLin.Rotated(anAxCLRot, anAngle);
2254 gp_Dir aDL = aConL.Position().Direction();
2255 gp_Dir aXDir = Tor.XAxis().Direction();
2256 //
2257 typeres = IntAna_Empty;
2258 //
2259 for (i = 0; i < 2; ++i) {
2260 if (i) {
2261 aXDir.Reverse();
2262 }
2263 aPCT.SetXYZ(aTorLoc.XYZ() + aRMaj*aXDir.XYZ());
2264 //
2265 aDist = aConL.Distance(aPCT);
2266 if (aDist > aRMin+Tol) {
2267 continue;
2268 }
2269 //
2270 typeres = IntAna_Circle;
2271 //
2272 gp_XYZ aPh = aPCT.XYZ() - aDist*aConL.Normal(aPCT).Direction().XYZ();
6092c0c8 2273 aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
7eed5d29 2274 //
2275 gp_Pnt aP;
6092c0c8 2276 gp_XYZ aDVal = aDt*aDL.XYZ();
7eed5d29 2277 aP.SetXYZ(aPh + aDVal);
2278 aParam[nbint] = aLin.Distance(aP);
2279 aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2280 aDir[nbint] = aTorAx.Direction();
2281 ++nbint;
6092c0c8 2282 if ((aDist < aRMin) && (aDt > Tol)) {
7eed5d29 2283 aP.SetXYZ(aPh - aDVal);
2284 aParam[nbint] = aLin.Distance(aP);
2285 aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2286 aDir[nbint] = aDir[nbint-1];
2287 ++nbint;
2288 }
2289 }
2290 //
2291 for (i = 0; i < nbint; ++i) {
2292 switch (i) {
2293 case 0:{
2294 pt1 = aPt[i];
2295 param1 = aParam[i];
2296 dir1 = aDir[i];
2297 break;
2298 }
2299 case 1:{
2300 pt2 = aPt[i];
2301 param2 = aParam[i];
2302 dir2 = aDir[i];
2303 break;
2304 }
2305 case 2:{
2306 pt3 = aPt[i];
2307 param3 = aParam[i];
2308 dir3 = aDir[i];
2309 break;
2310 }
2311 case 3:{
2312 pt4 = aPt[i];
2313 param4 = aParam[i];
2314 dir4 = aDir[i];
2315 break;
2316 }
2317 default:
2318 break;
2319 }
2320 }
2321}
2322
2323//=======================================================================
2324//function : IntAna_QuadQuadGeo
2325//purpose : Sphere - Torus
2326//=======================================================================
2327IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
2328 const gp_Torus& Tor,
2329 const Standard_Real Tol)
2330: done(Standard_False),
2331 nbint(0),
2332 typeres(IntAna_Empty),
2333 pt1(0,0,0),
2334 pt2(0,0,0),
2335 pt3(0,0,0),
2336 pt4(0,0,0),
2337 param1(0),
2338 param2(0),
2339 param3(0),
2340 param4(0),
2341 param1bis(0),
2342 param2bis(0),
2343 myCommonGen(Standard_False),
2344 myPChar(0,0,0)
2345{
2346 InitTolerances();
2347 Perform(Sph,Tor,Tol);
2348}
2349//=======================================================================
2350//function : Perform
2351//purpose :
2352//=======================================================================
2353void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
2354 const gp_Torus& Tor,
2355 const Standard_Real Tol)
2356{
2357 done = Standard_True;
2358 //
2359 Standard_Real aRMin, aRMaj;
2360 //
2361 aRMin = Tor.MinorRadius();
2362 aRMaj = Tor.MajorRadius();
2363 if (aRMin >= aRMaj) {
2364 typeres = IntAna_NoGeometricSolution;
2365 return;
2366 }
2367 //
2368 const gp_Ax1 aTorAx = Tor.Axis();
2369 const gp_Lin aLin(aTorAx);
2370 const gp_Pnt aSphLoc = Sph.Location();
2371 //
2372 if (aLin.Distance(aSphLoc) > myEPSILON_DISTANCE) {
2373 typeres = IntAna_NoGeometricSolution;
2374 return;
2375 }
2376 //
2377 Standard_Real aRSph, aDist;
2378 gp_Pnt aTorLoc;
2379 //
2380 gp_Dir aXDir = Tor.XAxis().Direction();
2381 aTorLoc.SetXYZ(aTorAx.Location().XYZ() + aRMaj*aXDir.XYZ());
2382 aRSph = Sph.Radius();
2383 //
2384 gp_Vec aVec12(aTorLoc, aSphLoc);
2385 aDist = aVec12.Magnitude();
2386 if (((aDist - Tol) > (aRMin + aRSph)) ||
2387 ((aDist + Tol) < Abs(aRMin - aRSph))) {
2388 typeres = IntAna_Empty;
2389 return;
2390 }
2391 //
2392 typeres = IntAna_Circle;
2393 //
2394 Standard_Real anAlpha, aBeta;
2395 //
2396 anAlpha = 0.5*(aRMin*aRMin - aRSph*aRSph + aDist*aDist ) / aDist;
2397 aBeta = Sqrt(Abs(aRMin*aRMin - anAlpha*anAlpha));
2398 //
2399 gp_Dir aDir12(aVec12);
2400 gp_XYZ aPh = aTorLoc.XYZ() + anAlpha*aDir12.XYZ();
2401 gp_Dir aDC = Tor.YAxis().Direction()^aDir12;
2402 //
2403 gp_Pnt aP;
2404 gp_XYZ aDVal = aBeta*aDC.XYZ();
2405 aP.SetXYZ(aPh + aDVal);
2406 param1 = aLin.Distance(aP);
2407 pt1.SetXYZ(aP.XYZ() - param1*aXDir.XYZ());
2408 dir1 = aTorAx.Direction();
2409 nbint = 1;
2410 if ((aDist < (aRSph + aRMin)) && (aDist > Abs(aRSph - aRMin)) &&
2411 (aDVal.Modulus() > Tol)) {
2412 aP.SetXYZ(aPh - aDVal);
2413 param2 = aLin.Distance(aP);
2414 pt2.SetXYZ(aP.XYZ() - param2*aXDir.XYZ());
2415 dir2 = dir1;
2416 nbint = 2;
2417 }
2418}
2419
2420//=======================================================================
2421//function : IntAna_QuadQuadGeo
2422//purpose : Torus - Torus
2423//=======================================================================
2424IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Torus& Tor1,
2425 const gp_Torus& Tor2,
2426 const Standard_Real Tol)
2427: done(Standard_False),
2428 nbint(0),
2429 typeres(IntAna_Empty),
2430 pt1(0,0,0),
2431 pt2(0,0,0),
2432 pt3(0,0,0),
2433 pt4(0,0,0),
2434 param1(0),
2435 param2(0),
2436 param3(0),
2437 param4(0),
2438 param1bis(0),
2439 param2bis(0),
2440 myCommonGen(Standard_False),
2441 myPChar(0,0,0)
2442{
2443 InitTolerances();
2444 Perform(Tor1,Tor2,Tol);
2445}
2446//=======================================================================
2447//function : Perform
2448//purpose :
2449//=======================================================================
2450void IntAna_QuadQuadGeo::Perform(const gp_Torus& Tor1,
2451 const gp_Torus& Tor2,
2452 const Standard_Real Tol)
2453{
2454 done = Standard_True;
2455 //
2456 Standard_Real aRMin1, aRMin2, aRMaj1, aRMaj2;
2457 //
2458 aRMin1 = Tor1.MinorRadius();
2459 aRMaj1 = Tor1.MajorRadius();
2460 aRMin2 = Tor2.MinorRadius();
2461 aRMaj2 = Tor2.MajorRadius();
2462 if (aRMin1 >= aRMaj1 || aRMin2 >= aRMaj2) {
2463 typeres = IntAna_NoGeometricSolution;
2464 return;
2465 }
2466 //
2467 const gp_Ax1 anAx1 = Tor1.Axis();
2468 const gp_Ax1 anAx2 = Tor2.Axis();
2469 //
2470 gp_Lin aL1(anAx1);
2471 if (!anAx1.IsParallel(anAx2, myEPSILON_AXES_PARA) ||
2472 (aL1.Distance(anAx2.Location()) > myEPSILON_DISTANCE)) {
2473 typeres = IntAna_NoGeometricSolution;
2474 return;
2475 }
2476 //
2477 gp_Pnt aLoc1, aLoc2;
2478 //
2479 aLoc1 = anAx1.Location();
2480 aLoc2 = anAx2.Location();
2481 //
2482 if (aLoc1.IsEqual(aLoc2, Tol) &&
2483 (Abs(aRMin1 - aRMin2) <= Tol) &&
2484 (Abs(aRMaj1 - aRMaj2) <= Tol)) {
2485 typeres = IntAna_Same;
2486 return;
2487 }
2488 //
2489 Standard_Real aDist;
2490 gp_Pnt aP1, aP2;
2491 //
2492 gp_Dir aXDir1 = Tor1.XAxis().Direction();
2493 aP1.SetXYZ(aLoc1.XYZ() + aRMaj1*aXDir1.XYZ());
2494 aP2.SetXYZ(aLoc2.XYZ() + aRMaj2*aXDir1.XYZ());
2495 //
2496 gp_Vec aV12(aP1, aP2);
2497 aDist = aV12.Magnitude();
2498 if (((aDist - Tol) > (aRMin1 + aRMin2)) ||
2499 ((aDist + Tol) < Abs(aRMin1 - aRMin2))) {
2500 typeres = IntAna_Empty;
2501 return;
2502 }
2503 //
2504 typeres = IntAna_Circle;
2505 //
2506 Standard_Real anAlpha, aBeta;
2507 //
2508 anAlpha = 0.5*(aRMin1*aRMin1 - aRMin2*aRMin2 + aDist*aDist ) / aDist;
2509 aBeta = Sqrt(Abs(aRMin1*aRMin1 - anAlpha*anAlpha));
2510 //
2511 gp_Dir aDir12(aV12);
2512 gp_XYZ aPh = aP1.XYZ() + anAlpha*aDir12.XYZ();
2513 gp_Dir aDC = Tor1.YAxis().Direction()^aDir12;
2514 //
2515 gp_Pnt aP;
2516 gp_XYZ aDVal = aBeta*aDC.XYZ();
2517 aP.SetXYZ(aPh + aDVal);
2518 param1 = aL1.Distance(aP);
2519 pt1.SetXYZ(aP.XYZ() - param1*aXDir1.XYZ());
2520 dir1 = anAx1.Direction();
2521 nbint = 1;
2522 if ((aDist < (aRMin1 + aRMin2)) && (aDist > Abs(aRMin1 - aRMin2)) &&
2523 aDVal.Modulus() > Tol) {
2524 aP.SetXYZ(aPh - aDVal);
2525 param2 = aL1.Distance(aP);
2526 pt2.SetXYZ(aP.XYZ() - param2*aXDir1.XYZ());
2527 dir2 = dir1;
2528 nbint = 2;
2529 }
2530}
2531
7fd59977 2532//=======================================================================
2533//function : Point
2534//purpose : Returns a Point
2535//=======================================================================
2536 gp_Pnt IntAna_QuadQuadGeo::Point(const Standard_Integer n) const
2537{
2538 if(!done) { StdFail_NotDone::Raise(); }
2539 if(n>nbint || n<1) { Standard_DomainError::Raise(); }
2540 if(typeres==IntAna_PointAndCircle) {
2541 if(n!=1) { Standard_DomainError::Raise(); }
2542 if(param1==0.0) return(pt1);
2543 return(pt2);
2544 }
2545 else if(typeres==IntAna_Point) {
2546 if(n==1) return(pt1);
2547 return(pt2);
2548 }
2549
2550 // WNT (what can you expect from MicroSoft ?)
2551 return gp_Pnt(0,0,0);
2552}
2553//=======================================================================
2554//function : Line
2555//purpose : Returns a Line
2556//=======================================================================
2557 gp_Lin IntAna_QuadQuadGeo::Line(const Standard_Integer n) const
2558{
2559 if(!done) { StdFail_NotDone::Raise(); }
2560 if((n>nbint) || (n<1) || (typeres!=IntAna_Line)) {
2561 Standard_DomainError::Raise();
2562 }
2563 if(n==1) { return(gp_Lin(pt1,dir1)); }
2564 else { return(gp_Lin(pt2,dir2)); }
2565}
2566//=======================================================================
2567//function : Circle
2568//purpose : Returns a Circle
2569//=======================================================================
2570 gp_Circ IntAna_QuadQuadGeo::Circle(const Standard_Integer n) const
2571{
2572 if(!done) { StdFail_NotDone::Raise(); }
2573 if(typeres==IntAna_PointAndCircle) {
2574 if(n!=1) { Standard_DomainError::Raise(); }
2575 if(param2==0.0) return(gp_Circ(DirToAx2(pt1,dir1),param1));
2576 return(gp_Circ(DirToAx2(pt2,dir2),param2));
2577 }
2578 else if((n>nbint) || (n<1) || (typeres!=IntAna_Circle)) {
2579 Standard_DomainError::Raise();
2580 }
7eed5d29 2581 if (n==1) { return(gp_Circ(DirToAx2(pt1,dir1),param1));}
2582 else if (n==2) { return(gp_Circ(DirToAx2(pt2,dir2),param2));}
2583 else if (n==3) { return(gp_Circ(DirToAx2(pt3,dir3),param3));}
2584 else { return(gp_Circ(DirToAx2(pt4,dir4),param4));}
7fd59977 2585}
2586
2587//=======================================================================
2588//function : Ellipse
2589//purpose : Returns a Elips
2590//=======================================================================
2591 gp_Elips IntAna_QuadQuadGeo::Ellipse(const Standard_Integer n) const
2592{
2593 if(!done) { StdFail_NotDone::Raise(); }
2594 if((n>nbint) || (n<1) || (typeres!=IntAna_Ellipse)) {
2595 Standard_DomainError::Raise();
2596 }
2597
2598 if(n==1) {
2599 Standard_Real R1=param1, R2=param1bis, aTmp;
2600 if (R1<R2) {
2601 aTmp=R1; R1=R2; R2=aTmp;
2602 }
2603 gp_Ax2 anAx2(pt1, dir1 ,dir2);
2604 gp_Elips anElips (anAx2, R1, R2);
2605 return anElips;
2606 }
2607 else {
2608 Standard_Real R1=param2, R2=param2bis, aTmp;
2609 if (R1<R2) {
2610 aTmp=R1; R1=R2; R2=aTmp;
2611 }
2612 gp_Ax2 anAx2(pt2, dir2 ,dir1);
2613 gp_Elips anElips (anAx2, R1, R2);
2614 return anElips;
2615 }
2616}
2617//=======================================================================
2618//function : Parabola
2619//purpose : Returns a Parabola
2620//=======================================================================
2621 gp_Parab IntAna_QuadQuadGeo::Parabola(const Standard_Integer n) const
2622{
2623 if(!done) {
2624 StdFail_NotDone::Raise();
2625 }
2626 if (typeres!=IntAna_Parabola) {
2627 Standard_DomainError::Raise();
2628 }
2629 if((n>nbint) || (n!=1)) {
2630 Standard_OutOfRange::Raise();
2631 }
2632 return(gp_Parab(gp_Ax2( pt1
7eed5d29 2633 ,dir1
2634 ,dir2)
2635 ,param1));
7fd59977 2636}
2637//=======================================================================
2638//function : Hyperbola
2639//purpose : Returns a Hyperbola
2640//=======================================================================
2641 gp_Hypr IntAna_QuadQuadGeo::Hyperbola(const Standard_Integer n) const
2642{
2643 if(!done) {
2644 StdFail_NotDone::Raise();
2645 }
2646 if((n>nbint) || (n<1) || (typeres!=IntAna_Hyperbola)) {
2647 Standard_DomainError::Raise();
2648 }
2649 if(n==1) {
2650 return(gp_Hypr(gp_Ax2( pt1
7eed5d29 2651 ,dir1
2652 ,dir2)
2653 ,param1,param1bis));
7fd59977 2654 }
2655 else {
2656 return(gp_Hypr(gp_Ax2( pt2
7eed5d29 2657 ,dir1
2658 ,dir2.Reversed())
2659 ,param2,param2bis));
7fd59977 2660 }
2661}
7fd59977 2662//=======================================================================
2663//function : HasCommonGen
2664//purpose :
2665//=======================================================================
7fd59977 2666Standard_Boolean IntAna_QuadQuadGeo::HasCommonGen() const
2667{
2668 return myCommonGen;
2669}
7fd59977 2670//=======================================================================
2671//function : PChar
2672//purpose :
2673//=======================================================================
7fd59977 2674const gp_Pnt& IntAna_QuadQuadGeo::PChar() const
2675{
2676 return myPChar;
2677}
77088633 2678//=======================================================================
2679//function : RefineDir
2680//purpose :
2681//=======================================================================
2682void RefineDir(gp_Dir& aDir)
2683{
2684 Standard_Integer k, m, n;
2685 Standard_Real aC[3];
2686 //
2687 aDir.Coord(aC[0], aC[1], aC[2]);
2688 //
2689 m=0;
2690 n=0;
2691 for (k=0; k<3; ++k) {
2692 if (aC[k]==1. || aC[k]==-1.) {
2693 ++m;
2694 }
2695 else if (aC[k]!=0.) {
2696 ++n;
2697 }
2698 }
2699 //
2700 if (m && n) {
2701 Standard_Real aEps, aR1, aR2, aNum;
2702 //
2703 aEps=RealEpsilon();
2704 aR1=1.-aEps;
2705 aR2=1.+aEps;
2706 //
2707 for (k=0; k<3; ++k) {
2708 m=(aC[k]>0.);
2709 aNum=(m)? aC[k] : -aC[k];
2710 if (aNum>aR1 && aNum<aR2) {
7eed5d29 2711 if (m) {
2712 aC[k]=1.;
2713 }
2714 else {
2715 aC[k]=-1.;
2716 }
2717 //
2718 aC[(k+1)%3]=0.;
2719 aC[(k+2)%3]=0.;
2720 break;
77088633 2721 }
2722 }
2723 aDir.SetCoord(aC[0], aC[1], aC[2]);
2724 }
2725}
7fd59977 2726
2727
2728
Open CASCADE Technology repository