0023985: There is no section between attached faces.
[occt.git] / src / IntAna / IntAna_IntConicQuad.cxx
1 // Created on: 1992-07-27
2 // Created by: Laurent BUCHARD
3 // Copyright (c) 1992-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
5 //
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
10 //
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13 //
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
20
21
22 #ifndef DEB
23 #define No_Standard_RangeError
24 #define No_Standard_OutOfRange
25 #endif
26
27 #define CREATE  IntAna_IntConicQuad::IntAna_IntConicQuad
28 #define PERFORM void IntAna_IntConicQuad::Perform
29
30
31
32 #include <IntAna_IntConicQuad.ixx>
33
34 #include <IntAna_QuadQuadGeo.hxx>
35
36 #include <IntAna2d_AnaIntersection.hxx>
37 #include <IntAna2d_IntPoint.hxx>
38 #include <IntAna_ResultType.hxx>
39
40 #include <gp_Vec.hxx>
41 #include <gp_Lin2d.hxx>
42 #include <gp_Circ2d.hxx>
43 #include <gp_Ax3.hxx>
44
45 #include <math_DirectPolynomialRoots.hxx>
46 #include <math_TrigonometricFunctionRoots.hxx>
47 #include <ElCLib.hxx>
48
49
50 static Standard_Real PIpPI = M_PI + M_PI;
51 //=============================================================================
52 //==                                          E m p t y   C o n s t r u c t o r
53 //== 
54 CREATE(void) {
55   done=Standard_False;
56 }
57 //=============================================================================
58 //==                                                 L i n e  -   Q u a d r i c  
59 //==
60 CREATE(const gp_Lin& L,const IntAna_Quadric& Quad) {
61   Perform(L,Quad);
62 }
63
64 PERFORM(const gp_Lin& L,const IntAna_Quadric& Quad) {
65
66   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte;
67   done=inquadric=parallel=Standard_False;
68
69   //----------------------------------------------------------------------
70   //-- Substitution de x=t Lx + Lx0       ( exprime dans                 )
71   //--                 y=t Ly + Ly0       (  le systeme de coordonnees   )
72   //--                 z=t Lz + Lz0       (  canonique                   )
73   //--
74   //-- Dans     Qxx x**2 + Qyy y**2 + Qzz z**2
75   //--          + 2 ( Qxy x y  + Qxz x z  + Qyz y z  )
76   //--          + 2 ( Qx x + Qy y + Qz z )
77   //--          + QCte
78   //--
79   //-- Done un polynome en t : A2 t**2 + A1 t + A0 avec :
80   //----------------------------------------------------------------------
81
82
83   Standard_Real Lx0,Ly0,Lz0,Lx,Ly,Lz;
84
85
86   nbpts=0;
87
88   L.Direction().Coord(Lx,Ly,Lz);
89   L.Location().Coord(Lx0,Ly0,Lz0);
90  
91   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte);
92   
93   Standard_Real A0=(QCte + Qxx*Lx0*Lx0 + Qyy*Ly0*Ly0  + Qzz*Lz0*Lz0
94           + 2.0* (  Lx0*( Qx + Qxy*Ly0 + Qxz*Lz0) 
95                   + Ly0*( Qy + Qyz*Lz0 )
96                   + Qz*Lz0 )
97           );
98           
99   
100   Standard_Real A1=2.0*( Lx*( Qx + Qxx*Lx0 + Qxy*Ly0 + Qxz*Lz0 )
101               +Ly*( Qy + Qxy*Lx0 + Qyy*Ly0 + Qyz*Lz0 )
102               +Lz*( Qz + Qxz*Lx0 + Qyz*Ly0 + Qzz*Lz0 ));
103
104   Standard_Real A2=(Qxx*Lx*Lx + Qyy*Ly*Ly + Qzz*Lz*Lz
105            + 2.0*( Lx*( Qxy*Ly + Qxz*Lz )
106                   + Qyz*Ly*Lz));
107
108   math_DirectPolynomialRoots LinQuadPol(A2,A1,A0);
109   
110   if( LinQuadPol.IsDone()) {
111     done=Standard_True;
112     if(LinQuadPol.InfiniteRoots()) {
113       inquadric=Standard_True;
114     }
115     else {
116       nbpts= LinQuadPol.NbSolutions();
117       
118       for(Standard_Integer i=1 ; i<=nbpts; i++) {
119         Standard_Real t= LinQuadPol.Value(i);
120         paramonc[i-1] = t;
121         pnts[i-1]=gp_Pnt( Lx0+Lx*t
122                          ,Ly0+Ly*t
123                          ,Lz0+Lz*t);
124       }
125     }
126   }
127 }
128
129 //=============================================================================
130 //==                                            C i r c l e   -   Q u a d r i c  
131 //==
132 CREATE(const gp_Circ& C,const IntAna_Quadric& Quad) {
133   Perform(C,Quad); 
134 }
135
136 PERFORM(const gp_Circ& C,const IntAna_Quadric& Quad) {
137   
138   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte;
139   
140   //----------------------------------------------------------------------
141   //-- Dans le repere liee a C.Position() : 
142   //-- xC = R * Cos[t]
143   //-- yC = R * Sin[t]
144   //-- zC = 0
145   //--
146   //-- On exprime la quadrique dans ce repere et on substitue 
147   //-- xC,yC et zC    a    x,y et z 
148   //-- 
149   //-- On Obtient un polynome en Cos[t] et Sin[t] de degre 2
150   //----------------------------------------------------------------------
151   done=inquadric=parallel=Standard_False;  
152
153   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte);
154   Quad.NewCoefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte,C.Position());
155   
156   Standard_Real R=C.Radius();
157   Standard_Real RR=R*R;
158   
159   Standard_Real P_CosCos = RR * Qxx;    //-- Cos Cos
160   Standard_Real P_SinSin = RR * Qyy;    //-- Sin Sin
161   Standard_Real P_Sin    = R  * Qy;     //-- 2 Sin
162   Standard_Real P_Cos    = R  * Qx;     //-- 2 Cos
163   Standard_Real P_CosSin = RR * Qxy;    //-- 2 Cos Sin
164   Standard_Real P_Cte    = QCte;        //-- 1
165   
166   math_TrigonometricFunctionRoots CircQuadPol( P_CosCos-P_SinSin
167                                               ,P_CosSin
168                                               ,P_Cos+P_Cos
169                                               ,P_Sin+P_Sin
170                                               ,P_Cte+P_SinSin
171                                               ,0.0,PIpPI);
172   
173   if(CircQuadPol.IsDone()) {
174     done=Standard_True;
175     if(CircQuadPol.InfiniteRoots()) {
176       inquadric=Standard_True;
177     }
178     else {
179       nbpts= CircQuadPol.NbSolutions();
180       
181       for(Standard_Integer i=1 ; i<=nbpts; i++) {
182         Standard_Real t= CircQuadPol.Value(i);
183         paramonc[i-1] = t;
184         pnts[i-1]     = ElCLib::CircleValue(t,C.Position(),R);
185       }
186     }
187   }
188 }
189
190
191 //=============================================================================
192 //==                                                  E l i p s - Q u a d r i c 
193 //==
194 CREATE(const gp_Elips& E,const IntAna_Quadric& Quad) {
195   Perform(E,Quad);
196 }
197
198 PERFORM(const gp_Elips& E,const IntAna_Quadric& Quad) {
199   
200   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte;
201   
202   done=inquadric=parallel=Standard_False;  
203   
204   //----------------------------------------------------------------------
205   //-- Dans le repere liee a E.Position() : 
206   //-- xE = R * Cos[t]
207   //-- yE = r * Sin[t]
208   //-- zE = 0
209   //--
210   //-- On exprime la quadrique dans ce repere et on substitue 
211   //-- xE,yE et zE    a    x,y et z 
212   //-- 
213   //-- On Obtient un polynome en Cos[t] et Sin[t] de degre 2
214   //----------------------------------------------------------------------
215   
216   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte);
217   Quad.NewCoefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte,E.Position());
218   
219   Standard_Real R=E.MajorRadius();  
220   Standard_Real r=E.MinorRadius();   
221
222   
223   Standard_Real P_CosCos = R*R * Qxx;    //-- Cos Cos
224   Standard_Real P_SinSin = r*r * Qyy;    //-- Sin Sin
225   Standard_Real P_Sin    = r   * Qy;     //-- 2 Sin
226   Standard_Real P_Cos    = R   * Qx;     //-- 2 Cos
227   Standard_Real P_CosSin = R*r * Qxy;    //-- 2 Cos Sin
228   Standard_Real P_Cte    = QCte;         //-- 1
229   
230   math_TrigonometricFunctionRoots ElipsQuadPol( P_CosCos-P_SinSin
231                                               ,P_CosSin
232                                               ,P_Cos+P_Cos
233                                               ,P_Sin+P_Sin
234                                               ,P_Cte+P_SinSin
235                                               ,0.0,PIpPI);
236   
237   if(ElipsQuadPol.IsDone()) {
238     done=Standard_True;
239     if(ElipsQuadPol.InfiniteRoots()) {
240       inquadric=Standard_True;
241     }
242     else {
243       nbpts= ElipsQuadPol.NbSolutions();
244       for(Standard_Integer i=1 ; i<=nbpts; i++) {
245         Standard_Real t= ElipsQuadPol.Value(i);
246         paramonc[i-1] = t;
247         pnts[i-1]     = ElCLib::EllipseValue(t,E.Position(),R,r);
248       }
249     }
250   }
251 }
252
253 //=============================================================================
254 //==                                                  P a r a b - Q u a d r i c 
255 //==
256 CREATE(const gp_Parab& P,const IntAna_Quadric& Quad) {
257   Perform(P,Quad);
258 }
259
260 PERFORM(const gp_Parab& P,const IntAna_Quadric& Quad) {
261   
262   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte;
263   
264   done=inquadric=parallel=Standard_False;
265
266   //----------------------------------------------------------------------
267   //-- Dans le repere liee a P.Position() : 
268   //-- xP = y*y / (2 p)
269   //-- yP = y
270   //-- zP = 0
271   //--
272   //-- On exprime la quadrique dans ce repere et on substitue 
273   //-- xP,yP et zP    a    x,y et z 
274   //-- 
275   //-- On Obtient un polynome en y de degre 4
276   //----------------------------------------------------------------------
277   
278   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte);
279   Quad.NewCoefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte,P.Position());
280   
281   Standard_Real f=P.Focal();
282   Standard_Real Un_Sur_2p = 0.25 / f;
283
284   Standard_Real A4 = Qxx * Un_Sur_2p * Un_Sur_2p;
285   Standard_Real A3 = (Qxy+Qxy) * Un_Sur_2p;
286   Standard_Real A2 = Qyy + (Qx+Qx) * Un_Sur_2p;
287   Standard_Real A1 = Qy+Qy;
288   Standard_Real A0 = QCte;
289   
290   math_DirectPolynomialRoots ParabQuadPol(A4,A3,A2,A1,A0);
291   
292   if( ParabQuadPol.IsDone()) {
293     done=Standard_True;
294     if(ParabQuadPol.InfiniteRoots()) {
295       inquadric=Standard_True;
296     }
297     else {
298       nbpts= ParabQuadPol.NbSolutions();
299       
300       for(Standard_Integer i=1 ; i<=nbpts; i++) {
301         Standard_Real t= ParabQuadPol.Value(i);
302         paramonc[i-1] = t;
303         pnts[i-1]     = ElCLib::ParabolaValue(t,P.Position(),f);
304       }
305     }
306   }  
307 }
308
309 //=============================================================================
310 //==                                                    H y p r - Q u a d r i c 
311 //==
312 CREATE(const gp_Hypr& H,const IntAna_Quadric& Quad) {
313   Perform(H,Quad);
314 }
315
316 PERFORM(const gp_Hypr& H,const IntAna_Quadric& Quad) {
317   
318   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte;
319
320   done=inquadric=parallel=Standard_False;  
321
322   //----------------------------------------------------------------------
323   //-- Dans le repere liee a P.Position() : 
324   //-- xH = R Ch[t]
325   //-- yH = r Sh[t]
326   //-- zH = 0
327   //--
328   //-- On exprime la quadrique dans ce repere et on substitue 
329   //-- xP,yP et zP    a    x,y et z 
330   //-- 
331   //-- On Obtient un polynome en Exp[t] de degre 4
332   //----------------------------------------------------------------------
333   
334   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte);
335   Quad.NewCoefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,QCte,H.Position());
336   
337   Standard_Real R=H.MajorRadius();
338   Standard_Real r=H.MinorRadius();
339
340   Standard_Real RR=R*R;
341   Standard_Real rr=r*r;
342   Standard_Real Rr=R*r;
343
344   Standard_Real A4 = RR * Qxx + Rr* ( Qxy + Qxy ) + rr * Qyy;
345   Standard_Real A3 = 4.0* ( R*Qx + r*Qy );
346   Standard_Real A2 = 2.0* ( (QCte+QCte) + Qxx*RR - Qyy*rr );
347   Standard_Real A1 = 4.0* ( R*Qx - r*Qy);
348   Standard_Real A0 = Qxx*RR - Rr*(Qxy+Qxy) + Qyy*rr; 
349   
350   math_DirectPolynomialRoots HyperQuadPol(A4,A3,A2,A1,A0);
351   
352   if( HyperQuadPol.IsDone()) {
353     done=Standard_True;
354     if(HyperQuadPol.InfiniteRoots()) {
355       inquadric=Standard_True;
356     }
357     else {
358       nbpts= HyperQuadPol.NbSolutions();
359       Standard_Integer bonnessolutions = 0;
360       for(Standard_Integer i=1 ; i<=nbpts; i++) {
361         Standard_Real t= HyperQuadPol.Value(i);
362         if(t>=RealEpsilon()) {
363           Standard_Real Lnt = Log(t);
364           paramonc[bonnessolutions] = Lnt;
365           pnts[bonnessolutions]     = ElCLib::HyperbolaValue(Lnt,H.Position(),R,r);
366           bonnessolutions++;
367         }
368       }
369       nbpts=bonnessolutions;
370     }
371   } 
372 }
373 //=============================================================================
374
375
376
377
378 IntAna_IntConicQuad::IntAna_IntConicQuad (const gp_Lin& L, const gp_Pln& P,
379                                           const Standard_Real Tolang,
380                                           const Standard_Real Tol,
381                                           const Standard_Real Len) {
382   Perform(L,P,Tolang,Tol,Len);
383 }
384
385
386 IntAna_IntConicQuad::IntAna_IntConicQuad (const gp_Circ& C, const gp_Pln& P,
387                                           const Standard_Real Tolang,
388                                           const Standard_Real Tol) {
389   Perform(C,P,Tolang,Tol);
390 }
391
392
393 IntAna_IntConicQuad::IntAna_IntConicQuad (const gp_Elips& E, const gp_Pln& P,
394                                           const Standard_Real Tolang,
395                                           const Standard_Real Tol) {
396   Perform(E,P,Tolang,Tol);
397 }
398
399
400 IntAna_IntConicQuad::IntAna_IntConicQuad (const gp_Parab& Pb, const gp_Pln& P,
401                                           const Standard_Real Tolang){
402   Perform(Pb,P,Tolang);
403 }
404
405
406 IntAna_IntConicQuad::IntAna_IntConicQuad (const gp_Hypr& H, const gp_Pln& P,
407                                           const Standard_Real Tolang){
408   Perform(H,P,Tolang);
409 }
410
411
412 void IntAna_IntConicQuad::Perform (const gp_Lin& L, const gp_Pln& P,
413                                    const Standard_Real Tolang,
414                                    const Standard_Real Tol,
415                                    const Standard_Real Len) {
416   
417   // Tolang represente la tolerance angulaire a partir de laquelle on considere
418   // que l angle entre 2 vecteurs est nul. On raisonnera sur le cosinus de cet
419   // angle, (on a Cos(t) equivalent a t au voisinage de Pi/2).
420   
421   done=Standard_False;
422   
423   Standard_Real A,B,C,D;
424   Standard_Real Al,Bl,Cl;
425   Standard_Real Dis,Direc;
426   
427   P.Coefficients(A,B,C,D);
428   gp_Pnt Orig(L.Location());
429   L.Direction().Coord(Al,Bl,Cl);
430
431   Direc=A*Al+B*Bl+C*Cl;
432   Dis = A*Orig.X() + B*Orig.Y() + C*Orig.Z() + D;
433   //
434   parallel=Standard_False;
435   if (Abs(Direc) < Tolang) {
436     parallel=Standard_True;
437     if (Len!=0 && Direc!=0) {
438       //check the distance from bounding point of the line to the plane
439       gp_Pnt aP1, aP2;
440       //
441       aP1.SetCoord(Orig.X()-Dis*A, Orig.Y()-Dis*B, Orig.Z()-Dis*C);
442       aP2.SetCoord(aP1.X()+Len*Al, aP1.Y()+Len*Bl, aP1.Z()+Len*Cl);
443       if (P.Distance(aP2) > Tol) {
444         parallel=Standard_False;
445       } 
446     }
447   }
448   if (parallel) {
449     if (Abs(Dis) < Tolang) {
450       inquadric=Standard_True;
451     }
452     else {
453       inquadric=Standard_False;
454     }
455   }
456   else {
457     parallel=Standard_False;
458     inquadric=Standard_False;
459     nbpts = 1;
460     paramonc [0] = - Dis/Direc;
461     pnts[0].SetCoord(Orig.X()+paramonc[0]*Al,
462                      Orig.Y()+paramonc[0]*Bl,
463                      Orig.Z()+paramonc[0]*Cl);
464   }
465   done=Standard_True;
466 }
467
468
469 void IntAna_IntConicQuad::Perform (const gp_Circ& C, const gp_Pln& P,
470                                   const Standard_Real Tolang,
471                                   const Standard_Real Tol)
472 {
473   
474   done=Standard_False;
475   
476   gp_Pln Plconic(gp_Ax3(C.Position()));
477   IntAna_QuadQuadGeo IntP(Plconic,P,Tolang,Tol);
478   if (!IntP.IsDone()) {return;}
479   if (IntP.TypeInter() == IntAna_Empty) {
480     parallel=Standard_True;
481     Standard_Real distmax = P.Distance(C.Location()) + C.Radius()*Tolang;
482     if (distmax < Tol) {
483       inquadric = Standard_True;
484     }
485     else {
486       inquadric = Standard_False;
487     }
488     done=Standard_True;
489   }
490   else     if(IntP.TypeInter() == IntAna_Same) { 
491     inquadric = Standard_True;
492     done = Standard_True;
493   }
494   else {
495     inquadric=Standard_False;
496     parallel=Standard_False;
497     gp_Lin Ligsol(IntP.Line(1));
498     
499     gp_Vec V0(Plconic.Location(),Ligsol.Location());
500     gp_Vec Axex(Plconic.Position().XDirection());
501     gp_Vec Axey(Plconic.Position().YDirection());
502     
503     gp_Pnt2d Orig(Axex.Dot(V0),Axey.Dot(V0));
504     gp_Vec2d Dire(Axex.Dot(Ligsol.Direction()),
505                   Axey.Dot(Ligsol.Direction()));
506     
507     gp_Lin2d Ligs(Orig,Dire);
508     gp_Pnt2d Pnt2dBid(0.0,0.0);
509     gp_Dir2d Dir2dBid(1.0,0.0);
510     gp_Ax2d Ax2dBid(Pnt2dBid,Dir2dBid);
511     gp_Circ2d Cir(Ax2dBid,C.Radius());
512     
513     IntAna2d_AnaIntersection Int2d(Ligs,Cir);
514     
515     if (!Int2d.IsDone()) {return;}
516     
517     nbpts=Int2d.NbPoints();
518     for (Standard_Integer i=1; i<=nbpts; i++) {
519       
520       gp_Pnt2d resul(Int2d.Point(i).Value());
521       Standard_Real X= resul.X();
522       Standard_Real Y= resul.Y();
523       pnts[i-1].SetCoord(Plconic.Location().X() + X*Axex.X() + Y*Axey.X(),
524                          Plconic.Location().Y() + X*Axex.Y() + Y*Axey.Y(),
525                          Plconic.Location().Z() + X*Axex.Z() + Y*Axey.Z());
526       paramonc[i-1]=Int2d.Point(i).ParamOnSecond();
527     }
528     done=Standard_True;
529   }
530 }
531
532
533 void IntAna_IntConicQuad::Perform (const gp_Elips& E, const gp_Pln& Pln,
534                                    const Standard_Real,
535                                    const Standard_Real){
536   Perform(E,Pln);
537 }
538
539
540 void IntAna_IntConicQuad::Perform (const gp_Parab& P, const gp_Pln& Pln,
541                                    const Standard_Real){
542   Perform(P,Pln);
543 }
544
545
546 void IntAna_IntConicQuad::Perform (const gp_Hypr& H, const gp_Pln& Pln,
547                                    const Standard_Real){
548   Perform(H,Pln);
549 }
550
551