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