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