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