projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0025307: Problems resulted from re-triangulation
[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//----------------------------------------------------------------------
22#ifndef DEB
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
ecc4f148 980 P2.SetCoord(P2t.X() - ProjP2OnDirCyl1*DirCyl.X(),
981 P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y(),
982 P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
7fd59977 983 //--
984 Standard_Real R1pR2=R1+R2;
ecc4f148 985 if(DistA1A2>(R1pR2+Tol))
986 {
7eed5d29 987 typeres=IntAna_Empty;
988 nbint=0;
7fd59977 989 }
ecc4f148 990 else if(DistA1A2>(R1pR2))
991 {
7eed5d29 992 //-- 1 Tangent line -------------------------------------OK
993 typeres=IntAna_Line;
994
995 nbint=1;
996 dir1=DirCyl;
997 Standard_Real R1_R1pR2=R1/R1pR2;
ecc4f148 998 pt1.SetCoord( P1.X() + R1_R1pR2 * (P2.X()-P1.X()),
999 P1.Y() + R1_R1pR2 * (P2.Y()-P1.Y()),
1000 P1.Z() + R1_R1pR2 * (P2.Z()-P1.Z()));
7fd59977 1001 }
ecc4f148 1002 else if(DistA1A2>RmR)
1003 {
7eed5d29 1004 //-- 2 lines ---------------------------------------------OK
1005 typeres=IntAna_Line;
1006 nbint=2;
1007 dir1=DirCyl;
1008 gp_Vec P1P2(P1,P2);
1009 gp_Dir DirA1A2=gp_Dir(P1P2);
1010 gp_Dir Ortho_dir1_P1P2 = dir1.Crossed(DirA1A2);
1011 dir2=dir1;
1012 Standard_Real Alpha=0.5*(R1*R1-R2*R2+DistA1A2*DistA1A2)/(DistA1A2);
1013
ecc4f148 1014 //Standard_Real Beta = Sqrt(R1*R1-Alpha*Alpha);
7eed5d29 1015 Standard_Real anSqrtArg = R1*R1-Alpha*Alpha;
1016 Standard_Real Beta = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
1017
ecc4f148 1018 if((Beta+Beta)<Tol)
1019 {
7eed5d29 1020 nbint=1;
1021 pt1.SetCoord( P1.X() + Alpha*DirA1A2.X()
1022 ,P1.Y() + Alpha*DirA1A2.Y()
1023 ,P1.Z() + Alpha*DirA1A2.Z());
1024 }
ecc4f148 1025 else
1026 {
1027 pt1.SetCoord( P1.X() + Alpha*DirA1A2.X() + Beta*Ortho_dir1_P1P2.X(),
1028 P1.Y() + Alpha*DirA1A2.Y() + Beta*Ortho_dir1_P1P2.Y(),
1029 P1.Z() + Alpha*DirA1A2.Z() + Beta*Ortho_dir1_P1P2.Z());
7eed5d29 1030
ecc4f148 1031 pt2.SetCoord( P1.X() + Alpha*DirA1A2.X() - Beta*Ortho_dir1_P1P2.X(),
1032 P1.Y() + Alpha*DirA1A2.Y() - Beta*Ortho_dir1_P1P2.Y(),
1033 P1.Z() + Alpha*DirA1A2.Z() - Beta*Ortho_dir1_P1P2.Z());
7eed5d29 1034 }
7fd59977 1035 }
ecc4f148 1036 else if(DistA1A2>(RmR-Tol))
1037 {
7eed5d29 1038 //-- 1 Tangent ------------------------------------------OK
1039 typeres=IntAna_Line;
1040 nbint=1;
1041 dir1=DirCyl;
1042 Standard_Real R1_RmR=R1/RmR;
7fd59977 1043
ecc4f148 1044 if(R1 < R2)
1045 R1_RmR = -R1_RmR;
7fd59977 1046
ecc4f148 1047 pt1.SetCoord( P1.X() + R1_RmR * (P2.X()-P1.X()),
1048 P1.Y() + R1_RmR * (P2.Y()-P1.Y()),
1049 P1.Z() + R1_RmR * (P2.Z()-P1.Z()));
7fd59977 1050 }
1051 else {
7eed5d29 1052 nbint=0;
1053 typeres=IntAna_Empty;
7fd59977 1054 }
1055 }
1056 }
1057 else { //-- No Parallel Axis ---------------------------------OK
1058 if((RmR_Relative<=myEPSILON_CYLINDER_DELTA_RADIUS)
ecc4f148 1059 && (DistA1A2 <= myEPSILON_CYLINDER_DELTA_DISTANCE))
1060 {
7fd59977 1061 //-- PI/2 between the two axis and Intersection
1062 //-- and identical radius
1063 typeres=IntAna_Ellipse;
1064 nbint=2;
1065 gp_Dir DirCyl1=Cyl1.Position().Direction();
1066 gp_Dir DirCyl2=Cyl2.Position().Direction();
1067 pt1=pt2=A1A2.PtIntersect();
1068
1069 Standard_Real A=DirCyl1.Angle(DirCyl2);
1070 Standard_Real B;
c6541a0c 1071 B=Abs(Sin(0.5*(M_PI-A)));
7fd59977 1072 A=Abs(Sin(0.5*A));
1073
ecc4f148 1074 if(A==0.0 || B==0.0)
1075 {
7eed5d29 1076 typeres=IntAna_Same;
1077 return;
7fd59977 1078 }
1079
7fd59977 1080 gp_Vec dircyl1(DirCyl1);gp_Vec dircyl2(DirCyl2);
1081 dir1 = gp_Dir(dircyl1.Added(dircyl2));
1082 dir2 = gp_Dir(dircyl1.Subtracted(dircyl2));
7eed5d29 1083
7fd59977 1084 param2 = Cyl1.Radius() / A;
1085 param1 = Cyl1.Radius() / B;
1086 param2bis= param1bis = Cyl1.Radius();
ecc4f148 1087 if(param1 < param1bis)
1088 {
1089 A=param1;
1090 param1=param1bis;
1091 param1bis=A;
7fd59977 1092 }
ecc4f148 1093
1094 if(param2 < param2bis)
1095 {
1096 A=param2;
1097 param2=param2bis;
1098 param2bis=A;
7fd59977 1099 }
1100 }
ecc4f148 1101 else
1102 {
1103 if(Abs(DistA1A2-Cyl1.Radius()-Cyl2.Radius())<Tol)
1104 {
7eed5d29 1105 typeres = IntAna_Point;
1106 Standard_Real d,p1,p2;
1107
1108 gp_Dir D1 = Cyl1.Axis().Direction();
1109 gp_Dir D2 = Cyl2.Axis().Direction();
1110 A1A2.Distance(d,p1,p2);
1111 gp_Pnt P = Cyl1.Axis().Location();
1112 gp_Pnt P1(P.X() - p1*D1.X(),
1113 P.Y() - p1*D1.Y(),
1114 P.Z() - p1*D1.Z());
ecc4f148 1115
7eed5d29 1116 P = Cyl2.Axis().Location();
ecc4f148 1117
7eed5d29 1118 gp_Pnt P2(P.X() - p2*D2.X(),
1119 P.Y() - p2*D2.Y(),
1120 P.Z() - p2*D2.Z());
ecc4f148 1121
7eed5d29 1122 gp_Vec P1P2(P1,P2);
1123 D1=gp_Dir(P1P2);
1124 p1=Cyl1.Radius();
ecc4f148 1125
7eed5d29 1126 pt1.SetCoord(P1.X() + p1*D1.X(),
1127 P1.Y() + p1*D1.Y(),
1128 P1.Z() + p1*D1.Z());
1129 nbint = 1;
7fd59977 1130 }
ecc4f148 1131 else
1132 {
7eed5d29 1133 typeres=IntAna_NoGeometricSolution;
7fd59977 1134 }
1135 }
1136 }
1137}
1138//=======================================================================
1139//function : IntAna_QuadQuadGeo
1140//purpose : Cylinder - Cone
1141//=======================================================================
a34f083b 1142IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1143 const gp_Cone& Con,
1144 const Standard_Real Tol)
7fd59977 1145: done(Standard_False),
1146 nbint(0),
1147 typeres(IntAna_Empty),
1148 pt1(0,0,0),
1149 pt2(0,0,0),
7eed5d29 1150 pt3(0,0,0),
1151 pt4(0,0,0),
7fd59977 1152 param1(0),
1153 param2(0),
7eed5d29 1154 param3(0),
1155 param4(0),
7fd59977 1156 param1bis(0),
1157 param2bis(0),
1158 myCommonGen(Standard_False),
1159 myPChar(0,0,0)
1160{
1161 InitTolerances();
1162 Perform(Cyl,Con,Tol);
1163}
1164//=======================================================================
1165//function : Perform
1166//purpose :
1167//=======================================================================
1168 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
7eed5d29 1169 const gp_Cone& Con,
1170 const Standard_Real )
7fd59977 1171{
1172 done=Standard_True;
1173 AxeOperator A1A2(Cyl.Axis(),Con.Axis());
1174 if(A1A2.Same()) {
1175 gp_Pnt Pt=Con.Apex();
1176 Standard_Real dist=Cyl.Radius()/(Tan(Con.SemiAngle()));
1177 gp_Dir dir=Cyl.Position().Direction();
1178 pt1.SetCoord( Pt.X() + dist*dir.X()
7eed5d29 1179 ,Pt.Y() + dist*dir.Y()
1180 ,Pt.Z() + dist*dir.Z());
7fd59977 1181 pt2.SetCoord( Pt.X() - dist*dir.X()
7eed5d29 1182 ,Pt.Y() - dist*dir.Y()
1183 ,Pt.Z() - dist*dir.Z());
7fd59977 1184 dir1=dir2=dir;
1185 param1=param2=Cyl.Radius();
1186 nbint=2;
1187 typeres=IntAna_Circle;
1188
1189 }
1190 else {
1191 typeres=IntAna_NoGeometricSolution;
1192 }
1193}
1194//=======================================================================
1195//function :
1196//purpose : Cylinder - Sphere
1197//=======================================================================
1198 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
7eed5d29 1199 const gp_Sphere& Sph,
1200 const Standard_Real Tol)
7fd59977 1201: done(Standard_False),
1202 nbint(0),
1203 typeres(IntAna_Empty),
1204 pt1(0,0,0),
1205 pt2(0,0,0),
7eed5d29 1206 pt3(0,0,0),
1207 pt4(0,0,0),
7fd59977 1208 param1(0),
1209 param2(0),
7eed5d29 1210 param3(0),
1211 param4(0),
7fd59977 1212 param1bis(0),
1213 param2bis(0),
1214 myCommonGen(Standard_False),
1215 myPChar(0,0,0)
1216{
1217 InitTolerances();
1218 Perform(Cyl,Sph,Tol);
1219}
1220//=======================================================================
1221//function : Perform
1222//purpose :
1223//=======================================================================
1224 void IntAna_QuadQuadGeo::Perform( const gp_Cylinder& Cyl
7eed5d29 1225 ,const gp_Sphere& Sph
1226 ,const Standard_Real)
7fd59977 1227{
1228 done=Standard_True;
1229 gp_Pnt Pt=Sph.Location();
1230 AxeOperator A1A2(Cyl.Axis(),Sph.Position().Axis());
1231 if((A1A2.Intersect() && Pt.Distance(A1A2.PtIntersect())==0.0 )
1232 || (A1A2.Same())) {
1233 if(Sph.Radius() < Cyl.Radius()) {
1234 typeres = IntAna_Empty;
1235 }
1236 else {
1237 Standard_Real dist=Sqrt( Sph.Radius() * Sph.Radius() - Cyl.Radius() * Cyl.Radius() );
1238 gp_Dir dir=Cyl.Position().Direction();
1239 dir1 = dir2 = dir;
1240 typeres=IntAna_Circle;
1241 pt1.SetCoord( Pt.X() + dist*dir.X()
7eed5d29 1242 ,Pt.Y() + dist*dir.Y()
1243 ,Pt.Z() + dist*dir.Z());
7fd59977 1244 nbint=1;
1245 param1 = Cyl.Radius();
1246 if(dist>RealEpsilon()) {
7eed5d29 1247 pt2.SetCoord( Pt.X() - dist*dir.X()
1248 ,Pt.Y() - dist*dir.Y()
1249 ,Pt.Z() - dist*dir.Z());
1250 param2=Cyl.Radius();
1251 nbint=2;
7fd59977 1252 }
1253 }
1254 }
1255 else {
1256 typeres=IntAna_NoGeometricSolution;
1257 }
1258}
1259
1260//=======================================================================
1261//function : IntAna_QuadQuadGeo
1262//purpose : Cone - Cone
1263//=======================================================================
1264 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con1,
7eed5d29 1265 const gp_Cone& Con2,
1266 const Standard_Real Tol)
7fd59977 1267: done(Standard_False),
1268 nbint(0),
1269 typeres(IntAna_Empty),
1270 pt1(0,0,0),
1271 pt2(0,0,0),
7eed5d29 1272 pt3(0,0,0),
1273 pt4(0,0,0),
7fd59977 1274 param1(0),
1275 param2(0),
7eed5d29 1276 param3(0),
1277 param4(0),
7fd59977 1278 param1bis(0),
1279 param2bis(0),
1280 myCommonGen(Standard_False),
1281 myPChar(0,0,0)
1282{
1283 InitTolerances();
1284 Perform(Con1,Con2,Tol);
1285}
1286//
1287//=======================================================================
1288//function : Perform
1289//purpose :
1290//=======================================================================
1291 void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con1,
7eed5d29 1292 const gp_Cone& Con2,
1293 const Standard_Real Tol)
7fd59977 1294{
1295 done=Standard_True;
1296 //
1297 Standard_Real tg1, tg2, aDA1A2, aTol2;
1298 gp_Pnt aPApex1, aPApex2;
4bd102b8 1299
1300 Standard_Real TOL_APEX_CONF = 1.e-10;
1301
7fd59977 1302 //
1303 tg1=Tan(Con1.SemiAngle());
1304 tg2=Tan(Con2.SemiAngle());
1305
1306 if((tg1 * tg2) < 0.) {
1307 tg2 = -tg2;
1308 }
1309 //
7fd59977 1310 aTol2=Tol*Tol;
1311 aPApex1=Con1.Apex();
1312 aPApex2=Con2.Apex();
1313 aDA1A2=aPApex1.SquareDistance(aPApex2);
7fd59977 1314 //
1315 AxeOperator A1A2(Con1.Axis(),Con2.Axis());
1316 //
1317 // 1
1318 if(A1A2.Same()) {
1319 //-- two circles
1320 Standard_Real x;
1321 gp_Pnt P=Con1.Apex();
1322 gp_Dir D=Con1.Position().Direction();
1323 Standard_Real d=gp_Vec(D).Dot(gp_Vec(P,Con2.Apex()));
1324
1325 if(Abs(tg1-tg2)>myEPSILON_ANGLE_CONE) {
4bd102b8 1326 if (fabs(d) < TOL_APEX_CONF) {
7eed5d29 1327 typeres = IntAna_Point;
1328 nbint = 1;
1329 pt1 = P;
1330 return;
4bd102b8 1331 }
7fd59977 1332 x=(d*tg2)/(tg1+tg2);
1333 pt1.SetCoord( P.X() + x*D.X()
7eed5d29 1334 ,P.Y() + x*D.Y()
1335 ,P.Z() + x*D.Z());
7fd59977 1336 param1=Abs(x*tg1);
1337
1338 x=(d*tg2)/(tg2-tg1);
1339 pt2.SetCoord( P.X() + x*D.X()
7eed5d29 1340 ,P.Y() + x*D.Y()
1341 ,P.Z() + x*D.Z());
7fd59977 1342 param2=Abs(x*tg1);
1343 dir1 = dir2 = D;
1344 nbint=2;
1345 typeres=IntAna_Circle;
1346 }
1347 else {
4bd102b8 1348 if (fabs(d) < TOL_APEX_CONF) {
7eed5d29 1349 typeres=IntAna_Same;
7fd59977 1350 }
1351 else {
7eed5d29 1352 typeres=IntAna_Circle;
1353 nbint=1;
1354 x=d*0.5;
1355 pt1.SetCoord( P.X() + x*D.X()
1356 ,P.Y() + x*D.Y()
1357 ,P.Z() + x*D.Z());
1358 param1 = Abs(x * tg1);
1359 dir1 = D;
7fd59977 1360 }
1361 }
1362 } //-- fin A1A2.Same
1363 // 2
1364 else if((Abs(tg1-tg2)<myEPSILON_ANGLE_CONE) && (A1A2.Parallel())) {
1365 //-- voir AnVer12mai98
1366 Standard_Real DistA1A2=A1A2.Distance();
1367 gp_Dir DA1=Con1.Position().Direction();
1368 gp_Vec O1O2(Con1.Apex(),Con2.Apex());
b045e6a4 1369 gp_Dir O1O2n(O1O2); // normalization of the vector before projection
1370 Standard_Real O1O2_DA1=gp_Vec(DA1).Dot(gp_Vec(O1O2n));
1371
1372 gp_Vec O1_Proj_A2(O1O2n.X()-O1O2_DA1*DA1.X(),
7eed5d29 1373 O1O2n.Y()-O1O2_DA1*DA1.Y(),
1374 O1O2n.Z()-O1O2_DA1*DA1.Z());
7fd59977 1375 gp_Dir DB1=gp_Dir(O1_Proj_A2);
b045e6a4 1376
7fd59977 1377 Standard_Real yO1O2=O1O2.Dot(gp_Vec(DA1));
1378 Standard_Real ABSTG1 = Abs(tg1);
1379 Standard_Real X2 = (DistA1A2/ABSTG1 - yO1O2)*0.5;
1380 Standard_Real X1 = X2+yO1O2;
1381
1382 gp_Pnt P1(Con1.Apex().X() + X1*( DA1.X() + ABSTG1*DB1.X()),
7eed5d29 1383 Con1.Apex().Y() + X1*( DA1.Y() + ABSTG1*DB1.Y()),
1384 Con1.Apex().Z() + X1*( DA1.Z() + ABSTG1*DB1.Z()));
7fd59977 1385
1386 gp_Pnt MO1O2(0.5*(Con1.Apex().X()+Con2.Apex().X()),
7eed5d29 1387 0.5*(Con1.Apex().Y()+Con2.Apex().Y()),
1388 0.5*(Con1.Apex().Z()+Con2.Apex().Z()));
7fd59977 1389 gp_Vec P1MO1O2(P1,MO1O2);
1390
1391 gp_Dir DA1_X_DB1=DA1.Crossed(DB1);
1392 gp_Dir OrthoPln = DA1_X_DB1.Crossed(gp_Dir(P1MO1O2));
1393
1394 IntAna_QuadQuadGeo INTER_QUAD_PLN(gp_Pln(P1,OrthoPln),Con1,Tol,Tol);
1395 if(INTER_QUAD_PLN.IsDone()) {
1396 switch(INTER_QUAD_PLN.TypeInter()) {
7eed5d29 1397 case IntAna_Ellipse: {
1398 typeres=IntAna_Ellipse;
1399 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1400 pt1 = E.Location();
1401 dir1 = E.Position().Direction();
1402 dir2 = E.Position().XDirection();
1403 param1 = E.MajorRadius();
1404 param1bis = E.MinorRadius();
1405 nbint = 1;
1406 break;
7fd59977 1407 }
1408 case IntAna_Circle: {
7eed5d29 1409 typeres=IntAna_Circle;
1410 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1411 pt1 = C.Location();
1412 dir1 = C.Position().XDirection();
1413 dir2 = C.Position().YDirection();
1414 param1 = C.Radius();
1415 nbint = 1;
1416 break;
7fd59977 1417 }
1418 case IntAna_Hyperbola: {
7eed5d29 1419 typeres=IntAna_Hyperbola;
1420 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1421 pt1 = pt2 = H.Location();
1422 dir1 = H.Position().Direction();
1423 dir2 = H.Position().XDirection();
1424 param1 = param2 = H.MajorRadius();
1425 param1bis = param2bis = H.MinorRadius();
1426 nbint = 2;
1427 break;
7fd59977 1428 }
1429 case IntAna_Line: {
7eed5d29 1430 typeres=IntAna_Line;
1431 gp_Lin H=INTER_QUAD_PLN.Line(1);
1432 pt1 = pt2 = H.Location();
1433 dir1 = dir2 = H.Position().Direction();
1434 param1 = param2 = 0.0;
1435 param1bis = param2bis = 0.0;
1436 nbint = 2;
1437 break;
7fd59977 1438 }
1439 default:
7eed5d29 1440 typeres=IntAna_NoGeometricSolution;
7fd59977 1441 }
1442 }
1443 }// else if((Abs(tg1-tg2)<EPSILON_ANGLE_CONE) && (A1A2.Parallel()))
7fd59977 1444 // 3
1445 else if (aDA1A2<aTol2) {
7fd59977 1446 //
1447 // When apices are coinsided there can be 3 possible cases
1448 // 3.1 - empty solution (iRet=0)
1449 // 3.2 - one line when cone1 touches cone2 (iRet=1)
1450 // 3.3 - two lines when cone1 intersects cone2 (iRet=2)
1451 //
1452 Standard_Integer iRet;
1453 Standard_Real aGamma, aBeta1, aBeta2;
1454 Standard_Real aD1, aR1, aTgBeta1, aTgBeta2, aHalfPI;
1455 Standard_Real aCosGamma, aSinGamma, aDx, aR2, aRD2, aD2;
1456 gp_Pnt2d aP0, aPA1, aP1, aPA2;
1457 gp_Vec2d aVAx2;
1458 gp_Ax1 aAx1, aAx2;
1459 //
1460 // Preliminary analysis. Determination of iRet
1461 //
1462 iRet=0;
c6541a0c 1463 aHalfPI=0.5*M_PI;
7fd59977 1464 aD1=1.;
1465 aPA1.SetCoord(aD1, 0.);
1466 aP0.SetCoord(0., 0.);
1467 //
1468 aAx1=Con1.Axis();
1469 aAx2=Con2.Axis();
1470 aGamma=aAx1.Angle(aAx2);
1471 if (aGamma>aHalfPI){
c6541a0c 1472 aGamma=M_PI-aGamma;
7fd59977 1473 }
1474 aCosGamma=Cos(aGamma);
1475 aSinGamma=Sin(aGamma);
1476 //
1477 aBeta1=Con1.SemiAngle();
1478 aTgBeta1=Tan(aBeta1);
1479 aTgBeta1=Abs(aTgBeta1);
1480 //
1481 aBeta2=Con2.SemiAngle();
1482 aTgBeta2=Tan(aBeta2);
1483 aTgBeta2=Abs(aTgBeta2);
1484 //
1485 aR1=aD1*aTgBeta1;
1486 aP1.SetCoord(aD1, aR1);
1487 //
1488 // PA2
1489 aVAx2.SetCoord(aCosGamma, aSinGamma);
1490 gp_Dir2d aDAx2(aVAx2);
1491 gp_Lin2d aLAx2(aP0, aDAx2);
1492 //
1493 gp_Vec2d aV(aP0, aP1);
1494 aDx=aVAx2.Dot(aV);
1495 aPA2=aP0.Translated(aDx*aDAx2);
1496 //
1497 // aR2
1498 aDx=aPA2.Distance(aP0);
1499 aR2=aDx*aTgBeta2;
1500 //
1501 // aRD2
1502 aRD2=aPA2.Distance(aP1);
1503 //
1504 if (aRD2>(aR2+Tol)) {
1505 iRet=0;
7fd59977 1506 typeres=IntAna_Empty; //nothing
4101383e 1507 return;
7fd59977 1508 }
1509 //
1510 iRet=1; //touch case => 1 line
1511 if (aRD2<(aR2-Tol)) {
1512 iRet=2;//intersection => couple of lines
1513 }
1514 //
1515 // Finding the solution in 3D
1516 //
1517 Standard_Real aDa;
1518 gp_Pnt aQApex1, aQA1, aQA2, aQX, aQX1, aQX2;
1519 gp_Dir aD3Ax1, aD3Ax2;
1520 gp_Lin aLin;
1521 IntAna_QuadQuadGeo aIntr;
1522 //
1523 aQApex1=Con1.Apex();
1524 aD3Ax1=aAx1.Direction();
1525 aQA1.SetCoord(aQApex1.X()+aD1*aD3Ax1.X(),
7eed5d29 1526 aQApex1.Y()+aD1*aD3Ax1.Y(),
1527 aQApex1.Z()+aD1*aD3Ax1.Z());
7fd59977 1528 //
1529 aDx=aD3Ax1.Dot(aAx2.Direction());
1530 if (aDx<0.) {
1531 aAx2.Reverse();
1532 }
1533 aD3Ax2=aAx2.Direction();
1534 //
1535 aD2=aD1*sqrt((1.+aTgBeta1*aTgBeta1)/(1.+aTgBeta2*aTgBeta2));
1536 //
1537 aQA2.SetCoord(aQApex1.X()+aD2*aD3Ax2.X(),
7eed5d29 1538 aQApex1.Y()+aD2*aD3Ax2.Y(),
1539 aQApex1.Z()+aD2*aD3Ax2.Z());
7fd59977 1540 //
1541 gp_Pln aPln1(aQA1, aD3Ax1);
1542 gp_Pln aPln2(aQA2, aD3Ax2);
1543 //
1544 aIntr.Perform(aPln1, aPln2, Tol, Tol);
a060129f 1545 if (!aIntr.IsDone() || 0 == aIntr.NbSolutions()) {
7fd59977 1546 iRet=-1; // just in case. it must not be so
1547 typeres=IntAna_NoGeometricSolution;
1548 return;
1549 }
1550 //
1551 aLin=aIntr.Line(1);
1552 const gp_Dir& aDLin=aLin.Direction();
1553 gp_Vec aVLin(aDLin);
1554 gp_Pnt aOrig=aLin.Location();
1555 gp_Vec aVr(aQA1, aOrig);
1556 aDx=aVLin.Dot(aVr);
1557 aQX=aOrig.Translated(aDx*aVLin);
1558 //
1559 // Final part
1560 //
1561 typeres=IntAna_Line;
1562 //
1563 param1=0.;
1564 param2 =0.;
1565 param1bis=0.;
1566 param2bis=0.;
1567 //
1568 if (iRet==1) {
1569 // one line
1570 nbint=1;
1571 pt1=aQApex1;
1572 gp_Vec aVX(aQApex1, aQX);
1573 dir1=gp_Dir(aVX);
7fd59977 1574 }
1575
1576 else {//iRet=2
1577 // two lines
1578 nbint=2;
1579 aDa=aQA1.Distance(aQX);
1580 aDx=sqrt(aR1*aR1-aDa*aDa);
1581 aQX1=aQX.Translated(aDx*aVLin);
1582 aQX2=aQX.Translated(-aDx*aVLin);
1583 //
1584 pt1=aQApex1;
1585 pt2=aQApex1;
1586 gp_Vec aVX1(aQApex1, aQX1);
1587 dir1=gp_Dir(aVX1);
1588 gp_Vec aVX2(aQApex1, aQX2);
1589 dir2=gp_Dir(aVX2);
7fd59977 1590 }
1591 } //else if (aDA1A2<aTol2) {
7fd59977 1592 //Case when cones have common generatrix
1593 else if(A1A2.Intersect()) {
1594 //Check if apex of one cone belongs another one
1595 Standard_Real u, v, tol2 = Tol*Tol;
1596 ElSLib::Parameters(Con2, aPApex1, u, v);
1597 gp_Pnt p = ElSLib::Value(u, v, Con2);
1598 if(aPApex1.SquareDistance(p) > tol2) {
1599 typeres=IntAna_NoGeometricSolution;
1600 return;
1601 }
1602 //
1603 ElSLib::Parameters(Con1, aPApex2, u, v);
1604 p = ElSLib::Value(u, v, Con1);
1605 if(aPApex2.SquareDistance(p) > tol2) {
1606 typeres=IntAna_NoGeometricSolution;
1607 return;
1608 }
1609
1610 //Cones have a common generatrix passing through apexes
1611 myCommonGen = Standard_True;
1612
1613 //common generatrix of cones
1614 gp_Lin aGen(aPApex1, gp_Dir(gp_Vec(aPApex1, aPApex2)));
1615
1616 //Intersection point of axes
1617 gp_Pnt aPAxeInt = A1A2.PtIntersect();
1618
1619 //Characteristic point of intersection curve
1620 u = ElCLib::Parameter(aGen, aPAxeInt);
1621 myPChar = ElCLib::Value(u, aGen);
1622
1623
1624 //Other generatrixes of cones laying in maximal plane
c6541a0c
D		 1625 gp_Lin aGen1 = aGen.Rotated(Con1.Axis(), M_PI);
1626 gp_Lin aGen2 = aGen.Rotated(Con2.Axis(), M_PI);
7fd59977 1627 //
1628 //Intersection point of generatrixes
1629 gp_Dir aN; //solution plane normal
1630 gp_Dir aD1 = aGen1.Direction();
1631
1632 gp_Dir aD2(aD1.Crossed(aGen.Direction()));
1633
1634 if(aD1.IsParallel(aGen2.Direction(), Precision::Angular())) {
1635 aN = aD1.Crossed(aD2);
1636 }
1637 else if(aGen1.SquareDistance(aGen2) > tol2) {
1638 //Something wrong ???
1639 typeres=IntAna_NoGeometricSolution;
1640 return;
1641 }
1642 else {
1643 gp_Dir D1 = aGen1.Position().Direction();
1644 gp_Dir D2 = aGen2.Position().Direction();
1645 gp_Pnt O1 = aGen1.Location();
1646 gp_Pnt O2 = aGen2.Location();
1647 Standard_Real D1DotD2 = D1.Dot(D2);
1648 Standard_Real aSin = 1.-D1DotD2*D1DotD2;
1649 gp_Vec O1O2 (O1,O2);
1650 Standard_Real U2 = (D1.XYZ()*(O1O2.Dot(D1))-(O1O2.XYZ())).Dot(D2.XYZ());
1651 U2 /= aSin;
1652 gp_Pnt aPGint(ElCLib::Value(U2, aGen2));
1653
1654 aD1 = gp_Dir(gp_Vec(aPGint, myPChar));
1655 aN = aD1.Crossed(aD2);
1656 }
1657 //Plane that must contain intersection curves
1658 gp_Pln anIntPln(myPChar, aN);
1659
1660 IntAna_QuadQuadGeo INTER_QUAD_PLN(anIntPln,Con1,Tol,Tol);
1661
1662 if(INTER_QUAD_PLN.IsDone()) {
1663 switch(INTER_QUAD_PLN.TypeInter()) {
7eed5d29 1664 case IntAna_Ellipse: {
1665 typeres=IntAna_Ellipse;
1666 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1667 pt1 = E.Location();
1668 dir1 = E.Position().Direction();
1669 dir2 = E.Position().XDirection();
1670 param1 = E.MajorRadius();
1671 param1bis = E.MinorRadius();
1672 nbint = 1;
1673 break;
7fd59977 1674 }
1675 case IntAna_Circle: {
7eed5d29 1676 typeres=IntAna_Circle;
1677 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1678 pt1 = C.Location();
1679 dir1 = C.Position().XDirection();
1680 dir2 = C.Position().YDirection();
1681 param1 = C.Radius();
1682 nbint = 1;
1683 break;
7fd59977 1684 }
1685 case IntAna_Parabola: {
7eed5d29 1686 typeres=IntAna_Parabola;
1687 gp_Parab Prb=INTER_QUAD_PLN.Parabola(1);
1688 pt1 = Prb.Location();
1689 dir1 = Prb.Position().Direction();
1690 dir2 = Prb.Position().XDirection();
1691 param1 = Prb.Focal();
1692 nbint = 1;
1693 break;
7fd59977 1694 }
1695 case IntAna_Hyperbola: {
7eed5d29 1696 typeres=IntAna_Hyperbola;
1697 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1698 pt1 = pt2 = H.Location();
1699 dir1 = H.Position().Direction();
1700 dir2 = H.Position().XDirection();
1701 param1 = param2 = H.MajorRadius();
1702 param1bis = param2bis = H.MinorRadius();
1703 nbint = 2;
1704 break;
7fd59977 1705 }
1706 default:
7eed5d29 1707 typeres=IntAna_NoGeometricSolution;
7fd59977 1708 }
1709 }
1710 }
4101383e 1711
7fd59977 1712 else {
1713 typeres=IntAna_NoGeometricSolution;
1714 }
1715}
1716//=======================================================================
1717//function : IntAna_QuadQuadGeo
1718//purpose : Sphere - Cone
1719//=======================================================================
1720 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
7eed5d29 1721 const gp_Cone& Con,
1722 const Standard_Real Tol)
7fd59977 1723: done(Standard_False),
1724 nbint(0),
1725 typeres(IntAna_Empty),
1726 pt1(0,0,0),
1727 pt2(0,0,0),
7eed5d29 1728 pt3(0,0,0),
1729 pt4(0,0,0),
7fd59977 1730 param1(0),
1731 param2(0),
7eed5d29 1732 param3(0),
1733 param4(0),
7fd59977 1734 param1bis(0),
1735 param2bis(0),
1736 myCommonGen(Standard_False),
1737 myPChar(0,0,0)
1738{
1739 InitTolerances();
1740 Perform(Sph,Con,Tol);
1741}
1742//=======================================================================
1743//function : Perform
1744//purpose :
1745//=======================================================================
1746 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
7eed5d29 1747 const gp_Cone& Con,
1748 const Standard_Real)
7fd59977 1749{
77088633 1750
1751 //
7fd59977 1752 done=Standard_True;
77088633 1753 //
7fd59977 1754 AxeOperator A1A2(Con.Axis(),Sph.Position().Axis());
1755 gp_Pnt Pt=Sph.Location();
77088633 1756 //
7fd59977 1757 if((A1A2.Intersect() && (Pt.Distance(A1A2.PtIntersect())==0.0))
1758 || A1A2.Same()) {
1759 gp_Pnt ConApex= Con.Apex();
1760 Standard_Real dApexSphCenter=Pt.Distance(ConApex);
1761 gp_Dir ConDir;
1762 if(dApexSphCenter>RealEpsilon()) {
1763 ConDir = gp_Dir(gp_Vec(ConApex,Pt));
1764 }
1765 else {
1766 ConDir = Con.Position().Direction();
1767 }
1768
1769 Standard_Real Rad=Sph.Radius();
1770 Standard_Real tga=Tan(Con.SemiAngle());
1771
1772
1773 //-- 2 circles
1774 //-- x: Roots of (x**2 + y**2 = Rad**2)
1775 //-- tga = y / (x+dApexSphCenter)
1776 Standard_Real tgatga = tga * tga;
1777 math_DirectPolynomialRoots Eq( 1.0+tgatga
7eed5d29 1778 ,2.0*tgatga*dApexSphCenter
1779 ,-Rad*Rad + dApexSphCenter*dApexSphCenter*tgatga);
7fd59977 1780 if(Eq.IsDone()) {
1781 Standard_Integer nbsol=Eq.NbSolutions();
1782 if(nbsol==0) {
7eed5d29 1783 typeres=IntAna_Empty;
7fd59977 1784 }
1785 else {
7eed5d29 1786 typeres=IntAna_Circle;
1787 if(nbsol>=1) {
1788 Standard_Real x = Eq.Value(1);
1789 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1790 nbint=1;
1791 pt1.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1792 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1793 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1794 param1 = tga * dApexSphCenterpx;
1795 param1 = Abs(param1);
1796 dir1 = ConDir;
1797 if(param1<=myEPSILON_MINI_CIRCLE_RADIUS) {
1798 typeres=IntAna_PointAndCircle;
1799 param1=0.0;
1800 }
1801 }
1802 if(nbsol>=2) {
1803 Standard_Real x=Eq.Value(2);
1804 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1805 nbint=2;
1806 pt2.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1807 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1808 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1809 param2 = tga * dApexSphCenterpx;
1810 param2 = Abs(param2);
1811 dir2=ConDir;
1812 if(param2<=myEPSILON_MINI_CIRCLE_RADIUS) {
1813 typeres=IntAna_PointAndCircle;
1814 param2=0.0;
1815 }
1816 }
7fd59977 1817 }
1818 }
1819 else {
1820 done=Standard_False;
1821 }
1822 }
1823 else {
1824 typeres=IntAna_NoGeometricSolution;
1825 }
1826}
1827
1828//=======================================================================
1829//function : IntAna_QuadQuadGeo
1830//purpose : Sphere - Sphere
1831//=======================================================================
1832 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Sphere& Sph1
7eed5d29 1833 ,const gp_Sphere& Sph2
1834 ,const Standard_Real Tol)
7fd59977 1835: done(Standard_False),
1836 nbint(0),
1837 typeres(IntAna_Empty),
1838 pt1(0,0,0),
1839 pt2(0,0,0),
7eed5d29 1840 pt3(0,0,0),
1841 pt4(0,0,0),
7fd59977 1842 param1(0),
1843 param2(0),
7eed5d29 1844 param3(0),
1845 param4(0),
7fd59977 1846 param1bis(0),
1847 param2bis(0),
1848 myCommonGen(Standard_False),
1849 myPChar(0,0,0)
1850{
1851 InitTolerances();
1852 Perform(Sph1,Sph2,Tol);
1853}
1854//=======================================================================
1855//function : Perform
1856//purpose :
1857//=======================================================================
1858 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph1,
7eed5d29 1859 const gp_Sphere& Sph2,
1860 const Standard_Real Tol)
7fd59977 1861{
1862 done=Standard_True;
1863 gp_Pnt O1=Sph1.Location();
1864 gp_Pnt O2=Sph2.Location();
1865 Standard_Real dO1O2=O1.Distance(O2);
1866 Standard_Real R1=Sph1.Radius();
1867 Standard_Real R2=Sph2.Radius();
1868 Standard_Real Rmin,Rmax;
1869 typeres=IntAna_Empty;
1870 param2bis=0.0; //-- pour eviter param2bis not used ....
1871
1872 if(R1>R2) { Rmin=R2; Rmax=R1; } else { Rmin=R1; Rmax=R2; }
1873
1874 if(dO1O2<=Tol && (Abs(R1-R2) <= Tol)) {
1875 typeres = IntAna_Same;
1876 }
1877 else {
1878 if(dO1O2<=Tol) { return; }
1879 gp_Dir Dir=gp_Dir(gp_Vec(O1,O2));
1880 Standard_Real t = Rmax - dO1O2 - Rmin;
1881
1882 //----------------------------------------------------------------------
1883 //-- |----------------- R1 --------------------|
1884 //-- |----dO1O2-----|-----------R2----------|
1885 //-- --->--<-- t
1886 //--
1887 //-- |------ R1 ------|---------dO1O2----------|
1888 //-- |-------------------R2-----------------------|
1889 //-- --->--<-- t
1890 //----------------------------------------------------------------------
1891 if(t >= 0.0 && t <=Tol) {
1892 typeres = IntAna_Point;
1893 nbint = 1;
1894 Standard_Real t2;
1895 if(R1==Rmax) t2=(R1 + (R2 + dO1O2)) * 0.5;
1896 else t2=(-R1+(dO1O2-R2))*0.5;
7eed5d29 1897
7fd59977 1898 pt1.SetCoord( O1.X() + t2*Dir.X()
7eed5d29 1899 ,O1.Y() + t2*Dir.Y()
1900 ,O1.Z() + t2*Dir.Z());
7fd59977 1901 }
1902 else {
1903 //-----------------------------------------------------------------
1904 //-- |----------------- dO1O2 --------------------|
1905 //-- |----R1-----|-----------R2----------|-Tol-|
1906 //--
1907 //-- |----------------- Rmax --------------------|
1908 //-- |----Rmin----|-------dO1O2-------|-Tol-|
1909 //--
1910 //-----------------------------------------------------------------
1911 if((dO1O2 > (R1+R2+Tol)) || (Rmax > (dO1O2+Rmin+Tol))) {
7eed5d29 1912 typeres=IntAna_Empty;
7fd59977 1913 }
1914 else {
7eed5d29 1915 //---------------------------------------------------------------
1916 //--
1917 //--
1918 //---------------------------------------------------------------
1919 Standard_Real Alpha=0.5*(R1*R1-R2*R2+dO1O2*dO1O2)/(dO1O2);
1920 Standard_Real Beta = R1*R1-Alpha*Alpha;
1921 Beta = (Beta>0.0)? Sqrt(Beta) : 0.0;
1922
1923 if(Beta<= myEPSILON_MINI_CIRCLE_RADIUS) {
1924 typeres = IntAna_Point;
1925 Alpha = (R1 + (dO1O2 - R2)) * 0.5;
1926 }
1927 else {
1928 typeres = IntAna_Circle;
1929 dir1 = Dir;
1930 param1 = Beta;
1931 }
1932 pt1.SetCoord( O1.X() + Alpha*Dir.X()
1933 ,O1.Y() + Alpha*Dir.Y()
1934 ,O1.Z() + Alpha*Dir.Z());
1935
1936 nbint=1;
7fd59977 1937 }
1938 }
1939 }
1940}
7eed5d29 1941
1942//=======================================================================
1943//function : IntAna_QuadQuadGeo
1944//purpose : Plane - Torus
1945//=======================================================================
1946IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& Pln,
1947 const gp_Torus& Tor,
1948 const Standard_Real Tol)
1949: done(Standard_False),
1950 nbint(0),
1951 typeres(IntAna_Empty),
1952 pt1(0,0,0),
1953 pt2(0,0,0),
1954 pt3(0,0,0),
1955 pt4(0,0,0),
1956 param1(0),
1957 param2(0),
1958 param3(0),
1959 param4(0),
1960 param1bis(0),
1961 param2bis(0),
1962 myCommonGen(Standard_False),
1963 myPChar(0,0,0)
1964{
1965 InitTolerances();
1966 Perform(Pln,Tor,Tol);
1967}
1968//=======================================================================
1969//function : Perform
1970//purpose :
1971//=======================================================================
1972void IntAna_QuadQuadGeo::Perform(const gp_Pln& Pln,
1973 const gp_Torus& Tor,
1974 const Standard_Real Tol)
1975{
1976 done = Standard_True;
1977 //
1978 Standard_Real aRMin, aRMaj;
1979 //
1980 aRMin = Tor.MinorRadius();
1981 aRMaj = Tor.MajorRadius();
1982 if (aRMin >= aRMaj) {
1983 typeres = IntAna_NoGeometricSolution;
1984 return;
1985 }
1986 //
1987 const gp_Ax1 aPlnAx = Pln.Axis();
1988 const gp_Ax1 aTorAx = Tor.Axis();
1989 //
1990 Standard_Boolean bParallel, bNormal;
1991 //
1992 bParallel = aTorAx.IsParallel(aPlnAx, myEPSILON_AXES_PARA);
1993 bNormal = !bParallel ? aTorAx.IsNormal(aPlnAx, myEPSILON_AXES_PARA) : Standard_False;
1994 if (!bNormal && !bParallel) {
1995 typeres = IntAna_NoGeometricSolution;
1996 return;
1997 }
1998 //
1999 Standard_Real aDist;
2000 //
2001 gp_Pnt aTorLoc = aTorAx.Location();
2002 if (bParallel) {
2003 Standard_Real aDt, X, Y, Z, A, B, C, D;
2004 //
2005 Pln.Coefficients(A,B,C,D);
2006 aTorLoc.Coord(X,Y,Z);
2007 aDist = A*X + B*Y + C*Z + D;
2008 //
2009 if ((Abs(aDist) - aRMin) > Tol) {
2010 typeres=IntAna_Empty;
2011 return;
2012 }
2013 //
2014 typeres = IntAna_Circle;
2015 //
2016 pt1.SetCoord(X - aDist*A, Y - aDist*B, Z - aDist*C);
2017 aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
2018 param1 = aRMaj + aDt;
2019 dir1 = aTorAx.Direction();
2020 nbint = 1;
2021 if ((Abs(aDist) < aRMin) && (aDt > Tol)) {
2022 pt2 = pt1;
2023 param2 = aRMaj - aDt;
2024 dir2 = dir1;
2025 nbint = 2;
2026 }
2027 }
2028 //
2029 else {
2030 aDist = Pln.Distance(aTorLoc);
2031 if (aDist > myEPSILON_DISTANCE) {
2032 typeres = IntAna_NoGeometricSolution;
2033 return;
2034 }
2035 //
2036 typeres = IntAna_Circle;
2037 param2 = param1 = aRMin;
2038 dir2 = dir1 = aPlnAx.Direction();
2039 nbint = 2;
2040 //
2041 gp_Dir aDir = aTorAx.Direction()^dir1;
2042 pt1.SetXYZ(aTorLoc.XYZ() + aRMaj*aDir.XYZ());
2043 pt2.SetXYZ(aTorLoc.XYZ() - aRMaj*aDir.XYZ());
2044 }
2045}
2046
2047//=======================================================================
2048//function : IntAna_QuadQuadGeo
2049//purpose : Cylinder - Torus
2050//=======================================================================
2051IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
2052 const gp_Torus& Tor,
2053 const Standard_Real Tol)
2054: done(Standard_False),
2055 nbint(0),
2056 typeres(IntAna_Empty),
2057 pt1(0,0,0),
2058 pt2(0,0,0),
2059 pt3(0,0,0),
2060 pt4(0,0,0),
2061 param1(0),
2062 param2(0),
2063 param3(0),
2064 param4(0),
2065 param1bis(0),
2066 param2bis(0),
2067 myCommonGen(Standard_False),
2068 myPChar(0,0,0)
2069{
2070 InitTolerances();
2071 Perform(Cyl,Tor,Tol);
2072}
2073//=======================================================================
2074//function : Perform
2075//purpose :
2076//=======================================================================
2077void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
2078 const gp_Torus& Tor,
2079 const Standard_Real Tol)
2080{
2081 done = Standard_True;
2082 //
2083 Standard_Real aRMin, aRMaj;
2084 //
2085 aRMin = Tor.MinorRadius();
2086 aRMaj = Tor.MajorRadius();
2087 if (aRMin >= aRMaj) {
2088 typeres = IntAna_NoGeometricSolution;
2089 return;
2090 }
2091 //
2092 const gp_Ax1 aCylAx = Cyl.Axis();
2093 const gp_Ax1 aTorAx = Tor.Axis();
2094 //
2095 const gp_Lin aLin(aTorAx);
2096 const gp_Pnt aLocCyl = Cyl.Location();
2097 //
2098 if (!aTorAx.IsParallel(aCylAx, myEPSILON_AXES_PARA) ||
2099 (aLin.Distance(aLocCyl) > myEPSILON_DISTANCE)) {
2100 typeres = IntAna_NoGeometricSolution;
2101 return;
2102 }
2103 //
2104 Standard_Real aRCyl;
2105 //
2106 aRCyl = Cyl.Radius();
2107 if (((aRCyl + Tol) < (aRMaj - aRMin)) || ((aRCyl - Tol) > (aRMaj + aRMin))) {
2108 typeres = IntAna_Empty;
2109 return;
2110 }
2111 //
2112 typeres = IntAna_Circle;
2113 //
2114 Standard_Real aDist = Sqrt(Abs(aRMin*aRMin - (aRCyl-aRMaj)*(aRCyl-aRMaj)));
2115 gp_XYZ aTorLoc = aTorAx.Location().XYZ();
2116 //
2117 dir1 = aTorAx.Direction();
2118 pt1.SetXYZ(aTorLoc + aDist*dir1.XYZ());
2119 param1 = aRCyl;
2120 nbint = 1;
2121 if ((aDist > Tol) && (aRCyl > (aRMaj - aRMin)) &&
2122 (aRCyl < (aRMaj + aRMin))) {
2123 dir2 = dir1;
2124 pt2.SetXYZ(aTorLoc - aDist*dir2.XYZ());
2125 param2 = param1;
2126 nbint = 2;
2127 }
2128}
2129
2130//=======================================================================
2131//function : IntAna_QuadQuadGeo
2132//purpose : Cone - Torus
2133//=======================================================================
2134IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con,
2135 const gp_Torus& Tor,
2136 const Standard_Real Tol)
2137: done(Standard_False),
2138 nbint(0),
2139 typeres(IntAna_Empty),
2140 pt1(0,0,0),
2141 pt2(0,0,0),
2142 pt3(0,0,0),
2143 pt4(0,0,0),
2144 param1(0),
2145 param2(0),
2146 param3(0),
2147 param4(0),
2148 param1bis(0),
2149 param2bis(0),
2150 myCommonGen(Standard_False),
2151 myPChar(0,0,0)
2152{
2153 InitTolerances();
2154 Perform(Con,Tor,Tol);
2155}
2156//=======================================================================
2157//function : Perform
2158//purpose :
2159//=======================================================================
2160void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con,
2161 const gp_Torus& Tor,
2162 const Standard_Real Tol)
2163{
2164 done = Standard_True;
2165 //
2166 Standard_Real aRMin, aRMaj;
2167 //
2168 aRMin = Tor.MinorRadius();
2169 aRMaj = Tor.MajorRadius();
2170 if (aRMin >= aRMaj) {
2171 typeres = IntAna_NoGeometricSolution;
2172 return;
2173 }
2174 //
2175 const gp_Ax1 aConAx = Con.Axis();
2176 const gp_Ax1 aTorAx = Tor.Axis();
2177 //
2178 const gp_Lin aLin(aTorAx);
2179 const gp_Pnt aConApex = Con.Apex();
2180 //
2181 if (!aTorAx.IsParallel(aConAx, myEPSILON_AXES_PARA) ||
2182 (aLin.Distance(aConApex) > myEPSILON_DISTANCE)) {
2183 typeres = IntAna_NoGeometricSolution;
2184 return;
2185 }
2186 //
6092c0c8 2187 Standard_Real anAngle, aDist, aParam[4], aDt;
7eed5d29 2188 Standard_Integer i;
2189 gp_Pnt aTorLoc, aPCT, aPN, aPt[4];
2190 gp_Dir aDir[4];
2191 //
2192 anAngle = Con.SemiAngle();
2193 aTorLoc = aTorAx.Location();
2194 //
2195 aPN.SetXYZ(aTorLoc.XYZ() + aRMaj*Tor.YAxis().Direction().XYZ());
2196 gp_Dir aDN (gp_Vec(aTorLoc, aPN));
2197 gp_Ax1 anAxCLRot(aConApex, aDN);
2198 gp_Lin aConL = aLin.Rotated(anAxCLRot, anAngle);
2199 gp_Dir aDL = aConL.Position().Direction();
2200 gp_Dir aXDir = Tor.XAxis().Direction();
2201 //
2202 typeres = IntAna_Empty;
2203 //
2204 for (i = 0; i < 2; ++i) {
2205 if (i) {
2206 aXDir.Reverse();
2207 }
2208 aPCT.SetXYZ(aTorLoc.XYZ() + aRMaj*aXDir.XYZ());
2209 //
2210 aDist = aConL.Distance(aPCT);
2211 if (aDist > aRMin+Tol) {
2212 continue;
2213 }
2214 //
2215 typeres = IntAna_Circle;
2216 //
2217 gp_XYZ aPh = aPCT.XYZ() - aDist*aConL.Normal(aPCT).Direction().XYZ();
6092c0c8 2218 aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
7eed5d29 2219 //
2220 gp_Pnt aP;
6092c0c8 2221 gp_XYZ aDVal = aDt*aDL.XYZ();
7eed5d29 2222 aP.SetXYZ(aPh + aDVal);
2223 aParam[nbint] = aLin.Distance(aP);
2224 aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2225 aDir[nbint] = aTorAx.Direction();
2226 ++nbint;
6092c0c8 2227 if ((aDist < aRMin) && (aDt > Tol)) {
7eed5d29 2228 aP.SetXYZ(aPh - aDVal);
2229 aParam[nbint] = aLin.Distance(aP);
2230 aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2231 aDir[nbint] = aDir[nbint-1];
2232 ++nbint;
2233 }
2234 }
2235 //
2236 for (i = 0; i < nbint; ++i) {
2237 switch (i) {
2238 case 0:{
2239 pt1 = aPt[i];
2240 param1 = aParam[i];
2241 dir1 = aDir[i];
2242 break;
2243 }
2244 case 1:{
2245 pt2 = aPt[i];
2246 param2 = aParam[i];
2247 dir2 = aDir[i];
2248 break;
2249 }
2250 case 2:{
2251 pt3 = aPt[i];
2252 param3 = aParam[i];
2253 dir3 = aDir[i];
2254 break;
2255 }
2256 case 3:{
2257 pt4 = aPt[i];
2258 param4 = aParam[i];
2259 dir4 = aDir[i];
2260 break;
2261 }
2262 default:
2263 break;
2264 }
2265 }
2266}
2267
2268//=======================================================================
2269//function : IntAna_QuadQuadGeo
2270//purpose : Sphere - Torus
2271//=======================================================================
2272IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
2273 const gp_Torus& Tor,
2274 const Standard_Real Tol)
2275: done(Standard_False),
2276 nbint(0),
2277 typeres(IntAna_Empty),
2278 pt1(0,0,0),
2279 pt2(0,0,0),
2280 pt3(0,0,0),
2281 pt4(0,0,0),
2282 param1(0),
2283 param2(0),
2284 param3(0),
2285 param4(0),
2286 param1bis(0),
2287 param2bis(0),
2288 myCommonGen(Standard_False),
2289 myPChar(0,0,0)
2290{
2291 InitTolerances();
2292 Perform(Sph,Tor,Tol);
2293}
2294//=======================================================================
2295//function : Perform
2296//purpose :
2297//=======================================================================
2298void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
2299 const gp_Torus& Tor,
2300 const Standard_Real Tol)
2301{
2302 done = Standard_True;
2303 //
2304 Standard_Real aRMin, aRMaj;
2305 //
2306 aRMin = Tor.MinorRadius();
2307 aRMaj = Tor.MajorRadius();
2308 if (aRMin >= aRMaj) {
2309 typeres = IntAna_NoGeometricSolution;
2310 return;
2311 }
2312 //
2313 const gp_Ax1 aTorAx = Tor.Axis();
2314 const gp_Lin aLin(aTorAx);
2315 const gp_Pnt aSphLoc = Sph.Location();
2316 //
2317 if (aLin.Distance(aSphLoc) > myEPSILON_DISTANCE) {
2318 typeres = IntAna_NoGeometricSolution;
2319 return;
2320 }
2321 //
2322 Standard_Real aRSph, aDist;
2323 gp_Pnt aTorLoc;
2324 //
2325 gp_Dir aXDir = Tor.XAxis().Direction();
2326 aTorLoc.SetXYZ(aTorAx.Location().XYZ() + aRMaj*aXDir.XYZ());
2327 aRSph = Sph.Radius();
2328 //
2329 gp_Vec aVec12(aTorLoc, aSphLoc);
2330 aDist = aVec12.Magnitude();
2331 if (((aDist - Tol) > (aRMin + aRSph)) ||
2332 ((aDist + Tol) < Abs(aRMin - aRSph))) {
2333 typeres = IntAna_Empty;
2334 return;
2335 }
2336 //
2337 typeres = IntAna_Circle;
2338 //
2339 Standard_Real anAlpha, aBeta;
2340 //
2341 anAlpha = 0.5*(aRMin*aRMin - aRSph*aRSph + aDist*aDist ) / aDist;
2342 aBeta = Sqrt(Abs(aRMin*aRMin - anAlpha*anAlpha));
2343 //
2344 gp_Dir aDir12(aVec12);
2345 gp_XYZ aPh = aTorLoc.XYZ() + anAlpha*aDir12.XYZ();
2346 gp_Dir aDC = Tor.YAxis().Direction()^aDir12;
2347 //
2348 gp_Pnt aP;
2349 gp_XYZ aDVal = aBeta*aDC.XYZ();
2350 aP.SetXYZ(aPh + aDVal);
2351 param1 = aLin.Distance(aP);
2352 pt1.SetXYZ(aP.XYZ() - param1*aXDir.XYZ());
2353 dir1 = aTorAx.Direction();
2354 nbint = 1;
2355 if ((aDist < (aRSph + aRMin)) && (aDist > Abs(aRSph - aRMin)) &&
2356 (aDVal.Modulus() > Tol)) {
2357 aP.SetXYZ(aPh - aDVal);
2358 param2 = aLin.Distance(aP);
2359 pt2.SetXYZ(aP.XYZ() - param2*aXDir.XYZ());
2360 dir2 = dir1;
2361 nbint = 2;
2362 }
2363}
2364
2365//=======================================================================
2366//function : IntAna_QuadQuadGeo
2367//purpose : Torus - Torus
2368//=======================================================================
2369IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Torus& Tor1,
2370 const gp_Torus& Tor2,
2371 const Standard_Real Tol)
2372: done(Standard_False),
2373 nbint(0),
2374 typeres(IntAna_Empty),
2375 pt1(0,0,0),
2376 pt2(0,0,0),
2377 pt3(0,0,0),
2378 pt4(0,0,0),
2379 param1(0),
2380 param2(0),
2381 param3(0),
2382 param4(0),
2383 param1bis(0),
2384 param2bis(0),
2385 myCommonGen(Standard_False),
2386 myPChar(0,0,0)
2387{
2388 InitTolerances();
2389 Perform(Tor1,Tor2,Tol);
2390}
2391//=======================================================================
2392//function : Perform
2393//purpose :
2394//=======================================================================
2395void IntAna_QuadQuadGeo::Perform(const gp_Torus& Tor1,
2396 const gp_Torus& Tor2,
2397 const Standard_Real Tol)
2398{
2399 done = Standard_True;
2400 //
2401 Standard_Real aRMin1, aRMin2, aRMaj1, aRMaj2;
2402 //
2403 aRMin1 = Tor1.MinorRadius();
2404 aRMaj1 = Tor1.MajorRadius();
2405 aRMin2 = Tor2.MinorRadius();
2406 aRMaj2 = Tor2.MajorRadius();
2407 if (aRMin1 >= aRMaj1 || aRMin2 >= aRMaj2) {
2408 typeres = IntAna_NoGeometricSolution;
2409 return;
2410 }
2411 //
2412 const gp_Ax1 anAx1 = Tor1.Axis();
2413 const gp_Ax1 anAx2 = Tor2.Axis();
2414 //
2415 gp_Lin aL1(anAx1);
2416 if (!anAx1.IsParallel(anAx2, myEPSILON_AXES_PARA) ||
2417 (aL1.Distance(anAx2.Location()) > myEPSILON_DISTANCE)) {
2418 typeres = IntAna_NoGeometricSolution;
2419 return;
2420 }
2421 //
2422 gp_Pnt aLoc1, aLoc2;
2423 //
2424 aLoc1 = anAx1.Location();
2425 aLoc2 = anAx2.Location();
2426 //
2427 if (aLoc1.IsEqual(aLoc2, Tol) &&
2428 (Abs(aRMin1 - aRMin2) <= Tol) &&
2429 (Abs(aRMaj1 - aRMaj2) <= Tol)) {
2430 typeres = IntAna_Same;
2431 return;
2432 }
2433 //
2434 Standard_Real aDist;
2435 gp_Pnt aP1, aP2;
2436 //
2437 gp_Dir aXDir1 = Tor1.XAxis().Direction();
2438 aP1.SetXYZ(aLoc1.XYZ() + aRMaj1*aXDir1.XYZ());
2439 aP2.SetXYZ(aLoc2.XYZ() + aRMaj2*aXDir1.XYZ());
2440 //
2441 gp_Vec aV12(aP1, aP2);
2442 aDist = aV12.Magnitude();
2443 if (((aDist - Tol) > (aRMin1 + aRMin2)) ||
2444 ((aDist + Tol) < Abs(aRMin1 - aRMin2))) {
2445 typeres = IntAna_Empty;
2446 return;
2447 }
2448 //
2449 typeres = IntAna_Circle;
2450 //
2451 Standard_Real anAlpha, aBeta;
2452 //
2453 anAlpha = 0.5*(aRMin1*aRMin1 - aRMin2*aRMin2 + aDist*aDist ) / aDist;
2454 aBeta = Sqrt(Abs(aRMin1*aRMin1 - anAlpha*anAlpha));
2455 //
2456 gp_Dir aDir12(aV12);
2457 gp_XYZ aPh = aP1.XYZ() + anAlpha*aDir12.XYZ();
2458 gp_Dir aDC = Tor1.YAxis().Direction()^aDir12;
2459 //
2460 gp_Pnt aP;
2461 gp_XYZ aDVal = aBeta*aDC.XYZ();
2462 aP.SetXYZ(aPh + aDVal);
2463 param1 = aL1.Distance(aP);
2464 pt1.SetXYZ(aP.XYZ() - param1*aXDir1.XYZ());
2465 dir1 = anAx1.Direction();
2466 nbint = 1;
2467 if ((aDist < (aRMin1 + aRMin2)) && (aDist > Abs(aRMin1 - aRMin2)) &&
2468 aDVal.Modulus() > Tol) {
2469 aP.SetXYZ(aPh - aDVal);
2470 param2 = aL1.Distance(aP);
2471 pt2.SetXYZ(aP.XYZ() - param2*aXDir1.XYZ());
2472 dir2 = dir1;
2473 nbint = 2;
2474 }
2475}
2476
7fd59977 2477//=======================================================================
2478//function : Point
2479//purpose : Returns a Point
2480//=======================================================================
2481 gp_Pnt IntAna_QuadQuadGeo::Point(const Standard_Integer n) const
2482{
2483 if(!done) { StdFail_NotDone::Raise(); }
2484 if(n>nbint || n<1) { Standard_DomainError::Raise(); }
2485 if(typeres==IntAna_PointAndCircle) {
2486 if(n!=1) { Standard_DomainError::Raise(); }
2487 if(param1==0.0) return(pt1);
2488 return(pt2);
2489 }
2490 else if(typeres==IntAna_Point) {
2491 if(n==1) return(pt1);
2492 return(pt2);
2493 }
2494
2495 // WNT (what can you expect from MicroSoft ?)
2496 return gp_Pnt(0,0,0);
2497}
2498//=======================================================================
2499//function : Line
2500//purpose : Returns a Line
2501//=======================================================================
2502 gp_Lin IntAna_QuadQuadGeo::Line(const Standard_Integer n) const
2503{
2504 if(!done) { StdFail_NotDone::Raise(); }
2505 if((n>nbint) || (n<1) || (typeres!=IntAna_Line)) {
2506 Standard_DomainError::Raise();
2507 }
2508 if(n==1) { return(gp_Lin(pt1,dir1)); }
2509 else { return(gp_Lin(pt2,dir2)); }
2510}
2511//=======================================================================
2512//function : Circle
2513//purpose : Returns a Circle
2514//=======================================================================
2515 gp_Circ IntAna_QuadQuadGeo::Circle(const Standard_Integer n) const
2516{
2517 if(!done) { StdFail_NotDone::Raise(); }
2518 if(typeres==IntAna_PointAndCircle) {
2519 if(n!=1) { Standard_DomainError::Raise(); }
2520 if(param2==0.0) return(gp_Circ(DirToAx2(pt1,dir1),param1));
2521 return(gp_Circ(DirToAx2(pt2,dir2),param2));
2522 }
2523 else if((n>nbint) || (n<1) || (typeres!=IntAna_Circle)) {
2524 Standard_DomainError::Raise();
2525 }
7eed5d29 2526 if (n==1) { return(gp_Circ(DirToAx2(pt1,dir1),param1));}
2527 else if (n==2) { return(gp_Circ(DirToAx2(pt2,dir2),param2));}
2528 else if (n==3) { return(gp_Circ(DirToAx2(pt3,dir3),param3));}
2529 else { return(gp_Circ(DirToAx2(pt4,dir4),param4));}
7fd59977 2530}
2531
2532//=======================================================================
2533//function : Ellipse
2534//purpose : Returns a Elips
2535//=======================================================================
2536 gp_Elips IntAna_QuadQuadGeo::Ellipse(const Standard_Integer n) const
2537{
2538 if(!done) { StdFail_NotDone::Raise(); }
2539 if((n>nbint) || (n<1) || (typeres!=IntAna_Ellipse)) {
2540 Standard_DomainError::Raise();
2541 }
2542
2543 if(n==1) {
2544 Standard_Real R1=param1, R2=param1bis, aTmp;
2545 if (R1<R2) {
2546 aTmp=R1; R1=R2; R2=aTmp;
2547 }
2548 gp_Ax2 anAx2(pt1, dir1 ,dir2);
2549 gp_Elips anElips (anAx2, R1, R2);
2550 return anElips;
2551 }
2552 else {
2553 Standard_Real R1=param2, R2=param2bis, aTmp;
2554 if (R1<R2) {
2555 aTmp=R1; R1=R2; R2=aTmp;
2556 }
2557 gp_Ax2 anAx2(pt2, dir2 ,dir1);
2558 gp_Elips anElips (anAx2, R1, R2);
2559 return anElips;
2560 }
2561}
2562//=======================================================================
2563//function : Parabola
2564//purpose : Returns a Parabola
2565//=======================================================================
2566 gp_Parab IntAna_QuadQuadGeo::Parabola(const Standard_Integer n) const
2567{
2568 if(!done) {
2569 StdFail_NotDone::Raise();
2570 }
2571 if (typeres!=IntAna_Parabola) {
2572 Standard_DomainError::Raise();
2573 }
2574 if((n>nbint) || (n!=1)) {
2575 Standard_OutOfRange::Raise();
2576 }
2577 return(gp_Parab(gp_Ax2( pt1
7eed5d29 2578 ,dir1
2579 ,dir2)
2580 ,param1));
7fd59977 2581}
2582//=======================================================================
2583//function : Hyperbola
2584//purpose : Returns a Hyperbola
2585//=======================================================================
2586 gp_Hypr IntAna_QuadQuadGeo::Hyperbola(const Standard_Integer n) const
2587{
2588 if(!done) {
2589 StdFail_NotDone::Raise();
2590 }
2591 if((n>nbint) || (n<1) || (typeres!=IntAna_Hyperbola)) {
2592 Standard_DomainError::Raise();
2593 }
2594 if(n==1) {
2595 return(gp_Hypr(gp_Ax2( pt1
7eed5d29 2596 ,dir1
2597 ,dir2)
2598 ,param1,param1bis));
7fd59977 2599 }
2600 else {
2601 return(gp_Hypr(gp_Ax2( pt2
7eed5d29 2602 ,dir1
2603 ,dir2.Reversed())
2604 ,param2,param2bis));
7fd59977 2605 }
2606}
7fd59977 2607//=======================================================================
2608//function : HasCommonGen
2609//purpose :
2610//=======================================================================
7fd59977 2611Standard_Boolean IntAna_QuadQuadGeo::HasCommonGen() const
2612{
2613 return myCommonGen;
2614}
7fd59977 2615//=======================================================================
2616//function : PChar
2617//purpose :
2618//=======================================================================
7fd59977 2619const gp_Pnt& IntAna_QuadQuadGeo::PChar() const
2620{
2621 return myPChar;
2622}
77088633 2623//=======================================================================
2624//function : RefineDir
2625//purpose :
2626//=======================================================================
2627void RefineDir(gp_Dir& aDir)
2628{
2629 Standard_Integer k, m, n;
2630 Standard_Real aC[3];
2631 //
2632 aDir.Coord(aC[0], aC[1], aC[2]);
2633 //
2634 m=0;
2635 n=0;
2636 for (k=0; k<3; ++k) {
2637 if (aC[k]==1. || aC[k]==-1.) {
2638 ++m;
2639 }
2640 else if (aC[k]!=0.) {
2641 ++n;
2642 }
2643 }
2644 //
2645 if (m && n) {
2646 Standard_Real aEps, aR1, aR2, aNum;
2647 //
2648 aEps=RealEpsilon();
2649 aR1=1.-aEps;
2650 aR2=1.+aEps;
2651 //
2652 for (k=0; k<3; ++k) {
2653 m=(aC[k]>0.);
2654 aNum=(m)? aC[k] : -aC[k];
2655 if (aNum>aR1 && aNum<aR2) {
7eed5d29 2656 if (m) {
2657 aC[k]=1.;
2658 }
2659 else {
2660 aC[k]=-1.;
2661 }
2662 //
2663 aC[(k+1)%3]=0.;
2664 aC[(k+2)%3]=0.;
2665 break;
77088633 2666 }
2667 }
2668 aDir.SetCoord(aC[0], aC[1], aC[2]);
2669 }
2670}
7fd59977 2671
2672
2673
Open CASCADE Technology repository