0023158: ApproxInt_PrmPrmSvSurfaces raises FPE (division by zero) signal
[occt.git] / src / ApproxInt / ApproxInt_PrmPrmSvSurfaces.gxx
1 // Created on: 1993-03-17
2 // Created by: Laurent BUCHARD
3 // Copyright (c) 1993-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 #define TOLTANGENCY 0.0000000001
23
24
25 #include <TColStd_Array1OfReal.hxx>
26 #include <math_FunctionSetRoot.hxx>
27 #include <Precision.hxx>
28
29 #define Debug(expr)  cout<<" expr :"<<expr;
30 #define MySurf1 MyIntersectionOn2S.Function().AuxillarSurface1()
31 #define MySurf2 MyIntersectionOn2S.Function().AuxillarSurface2()
32
33
34 //--------------------------------------------------------------------------------
35 ApproxInt_PrmPrmSvSurfaces::ApproxInt_PrmPrmSvSurfaces( const ThePSurface& Surf1
36                                                        ,const ThePSurface& Surf2):
37        MyHasBeenComputed(Standard_False),
38        MyHasBeenComputedbis(Standard_False),
39        MyIntersectionOn2S(Surf1,Surf2,TOLTANGENCY)
40
41 }
42 //--------------------------------------------------------------------------------
43 Standard_Boolean ApproxInt_PrmPrmSvSurfaces::Compute( Standard_Real& u1
44                                                      ,Standard_Real& v1
45                                                      ,Standard_Real& u2
46                                                      ,Standard_Real& v2
47                                                      ,gp_Pnt& P
48                                                      ,gp_Vec& Tg
49                                                      ,gp_Vec2d& Tguv1
50                                                      ,gp_Vec2d& Tguv2) { 
51   
52   Standard_Real tu1=u1;
53   Standard_Real tu2=u2;
54   Standard_Real tv1=v1;
55   Standard_Real tv2=v2;
56   
57   if(MyHasBeenComputed) { 
58     if(  (MyParOnS1.X()==u1)&&(MyParOnS1.Y()==v1)
59        &&(MyParOnS2.X()==u2)&&(MyParOnS2.Y()==v2)) {
60       return(MyIsTangent);
61     }
62     else if(MyHasBeenComputedbis == Standard_False) { 
63       MyTgbis         = MyTg;
64       MyTguv1bis      = MyTguv1;
65       MyTguv2bis      = MyTguv2;
66       MyPntbis        = MyPnt;
67       MyParOnS1bis    = MyParOnS1;
68       MyParOnS2bis    = MyParOnS2;
69       MyIsTangentbis  = MyIsTangent;
70       MyHasBeenComputedbis = MyHasBeenComputed; 
71     }
72   }
73   if(MyHasBeenComputedbis) { 
74     if(  (MyParOnS1bis.X()==u1)&&(MyParOnS1bis.Y()==v1)
75        &&(MyParOnS2bis.X()==u2)&&(MyParOnS2bis.Y()==v2)) {
76
77       gp_Vec            TV(MyTg);
78       gp_Vec2d          TV1(MyTguv1);
79       gp_Vec2d          TV2(MyTguv2);
80       gp_Pnt            TP(MyPnt);
81       gp_Pnt2d          TP1(MyParOnS1);
82       gp_Pnt2d          TP2(MyParOnS2);
83       Standard_Boolean  TB=MyIsTangent;
84
85       MyTg        = MyTgbis;
86       MyTguv1     = MyTguv1bis;
87       MyTguv2     = MyTguv2bis;
88       MyPnt       = MyPntbis;
89       MyParOnS1   = MyParOnS1bis;
90       MyParOnS2   = MyParOnS2bis;
91       MyIsTangent = MyIsTangentbis;
92
93       MyTgbis         = TV;
94       MyTguv1bis      = TV1;
95       MyTguv2bis      = TV2;
96       MyPntbis        = TP;
97       MyParOnS1bis    = TP1;
98       MyParOnS2bis    = TP2;
99       MyIsTangentbis  = TB;
100
101       return(MyIsTangent);
102     }
103   }
104
105
106   MyIsTangent = Standard_True;
107
108   static TColStd_Array1OfReal Param(1,4);
109   Param(1) = u1; Param(2) = v1;
110   Param(3) = u2; Param(4) = v2;
111   math_FunctionSetRoot  Rsnld(MyIntersectionOn2S.Function());
112   MyIntersectionOn2S.Perform(Param,Rsnld);
113   if (!MyIntersectionOn2S.IsDone())  { 
114     MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
115     return(Standard_False);
116   }
117   if (MyIntersectionOn2S.IsEmpty()) {
118     MyIsTangent=Standard_False;
119     //cout<<"\n----- Parametree Parametree : IsEmpty ds Compute "<<endl;
120     //Debug(u1); Debug(u2); Debug(v1); Debug(v2);   cout<<endl;
121     MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
122     return(Standard_False);
123   }
124   MyHasBeenComputed = Standard_True;
125   MyPnt = P = MyIntersectionOn2S.Point().Value();
126   
127   MyIntersectionOn2S.Point().Parameters(u1,v1,u2,v2);
128   MyParOnS1.SetCoord(tu1,tv1);
129   MyParOnS2.SetCoord(tu2,tv2);
130   
131   if(MyIntersectionOn2S.IsTangent()) { 
132     MyIsTangent=Standard_False;
133     MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
134     return(Standard_False); 
135   }
136   MyTg    = Tg    = MyIntersectionOn2S.Direction();
137   MyTguv1 = Tguv1 = MyIntersectionOn2S.DirectionOnS1();
138   MyTguv2 = Tguv2 = MyIntersectionOn2S.DirectionOnS2();
139
140   //----------------------------------------------------------------------
141   //-- Si ( Tg )    TU et TV sont normes 
142   //-- 
143   //-- On a    Tg   =  DeltaU  *  TU    +   DeltaV  *  TV
144   //-- 
145   //-- soit :  Tg.TU  =  DeltaU  TU.TU  +   DeltaV  TU.TV
146   //--         Tg.TV  =  DeltaU  TV.TU  +   DeltaV  TV.TV 
147   //-- 
148   //-- Donc : 
149   //--
150   //--               Tg.TU TV.TV  - Tg.TV * TU.TV
151   //--   DeltaU = -------------------------------
152   //--               TU.TU TV.TV  - (TU.TV)**2
153   //-- 
154   //--               Tg.TV TU.TU  - Tg.TU * TU.TV
155   //--   DeltaV = -------------------------------
156   //--               TU.TU TV.TV  - (TU.TV)**2
157   //--
158   //--
159
160   Tg.Normalize();    MyTg = Tg; 
161
162   Standard_Real DeltaU,DeltaV;
163   gp_Vec TU,TV;
164   gp_Pnt Pbid;
165   Standard_Real TUTV,TgTU,TgTV,TUTU,TVTV,DIS;
166   //------------------------------------------------------------
167   //-- Calcul de Tguv1
168   //--
169   ThePSurfaceTool::D1(MySurf1,u1,v1,Pbid,TU,TV);
170   
171   TUTU = TU.Dot(TU);
172   TVTV = TV.Dot(TV);
173   TUTV = TU.Dot(TV);
174   TgTU = Tg.Dot(TU);
175   TgTV = Tg.Dot(TV);
176   DIS  = TUTU * TVTV - TUTV * TUTV;
177   if(fabs(DIS)<Precision::Angular()) { 
178     MyIsTangent=Standard_False;
179     MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
180     return(Standard_False); 
181   }
182   
183   DeltaU = (TgTU * TVTV - TgTV * TUTV ) / DIS ; 
184   DeltaV = (TgTV * TUTU - TgTU * TUTV ) / DIS ;
185
186   Tguv1.SetCoord(DeltaU,DeltaV);  MyTguv1 = Tguv1;
187
188   //------------------------------------------------------------
189   //-- Calcul de Tguv2
190   //--  
191   ThePSurfaceTool::D1(MySurf2,u2,v2,Pbid,TU,TV);
192
193   TUTU = TU.Dot(TU);
194   TVTV = TV.Dot(TV);
195   TUTV = TU.Dot(TV);
196   TgTU = Tg.Dot(TU);
197   TgTV = Tg.Dot(TV);
198   DIS  = TUTU * TVTV - TUTV * TUTV;
199   if(fabs(DIS)<Precision::Angular()) { 
200     MyIsTangent=Standard_False;
201     MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
202     return(Standard_False); 
203   }
204
205   DeltaU = (TgTU * TVTV - TgTV * TUTV ) / DIS ; 
206   DeltaV = (TgTV * TUTU - TgTU * TUTV ) / DIS ;
207   
208   Tguv2.SetCoord(DeltaU,DeltaV);  MyTguv2 = Tguv2;
209
210   return(Standard_True);
211 }
212 //--------------------------------------------------------------------------------
213 void ApproxInt_PrmPrmSvSurfaces::Pnt(const Standard_Real u1,
214                                      const Standard_Real v1,
215                                      const Standard_Real u2,
216                                      const Standard_Real v2,
217                                      gp_Pnt& P) { 
218   gp_Pnt aP;
219   gp_Vec aT;
220   gp_Vec2d aTS1,aTS2;
221   Standard_Real tu1=u1;
222   Standard_Real tu2=u2;
223   Standard_Real tv1=v1;
224   Standard_Real tv2=v2;
225 #ifdef DEB
226   Standard_Boolean t=this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
227 #else
228   this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
229 #endif
230   P=MyPnt;
231 }
232 //--------------------------------------------------------------------------------
233 Standard_Boolean ApproxInt_PrmPrmSvSurfaces::Tangency(const Standard_Real u1,
234                                                       const Standard_Real v1,
235                                                       const Standard_Real u2,
236                                                       const Standard_Real v2,
237                                                       gp_Vec& T) { 
238   gp_Pnt aP;
239   gp_Vec aT;
240   gp_Vec2d aTS1,aTS2;
241   Standard_Real tu1=u1;
242   Standard_Real tu2=u2;
243   Standard_Real tv1=v1;
244   Standard_Real tv2=v2;
245   Standard_Boolean t=this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
246   T=MyTg;
247   return(t);
248 }
249 //--------------------------------------------------------------------------------
250 Standard_Boolean ApproxInt_PrmPrmSvSurfaces::TangencyOnSurf1(const Standard_Real u1,
251                                                              const Standard_Real v1,
252                                                              const Standard_Real u2,
253                                                              const Standard_Real v2,
254                                                              gp_Vec2d& T) { 
255   gp_Pnt aP;
256   gp_Vec aT;
257   gp_Vec2d aTS1,aTS2;
258   Standard_Real tu1=u1;
259   Standard_Real tu2=u2;
260   Standard_Real tv1=v1;
261   Standard_Real tv2=v2;
262   Standard_Boolean t=this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
263   T=MyTguv1;
264   return(t);
265 }
266 //--------------------------------------------------------------------------------
267 Standard_Boolean ApproxInt_PrmPrmSvSurfaces::TangencyOnSurf2(const Standard_Real u1,
268                                                              const Standard_Real v1,
269                                                              const Standard_Real u2,
270                                                              const Standard_Real v2,
271                                                              gp_Vec2d& T) { 
272   gp_Pnt aP;
273   gp_Vec aT;
274   gp_Vec2d aTS1,aTS2;
275   Standard_Real tu1=u1;
276   Standard_Real tu2=u2;
277   Standard_Real tv1=v1;
278   Standard_Real tv2=v2;
279   Standard_Boolean t=this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
280   T=MyTguv2;
281   return(t);
282 }
283 //--------------------------------------------------------------------------------
284
285
286
287
288 #if 0 
289   //------------------------------------------------------------
290   //-- Calcul de Tguv1
291   //--
292   ThePSurfaceTool::D1(MySurf1,u1,v1,P,TU,TV);
293   
294   TUTV = TU.Dot(TV);
295   TgTU = Tg.Dot(TU);
296   TgTV = Tg.Dot(TV);
297   UmTUTV2 = 1.0 - TUTV * TUTV;
298   
299   DeltaU = (TgTU - TgTV * TUTV ) / UmTUTV2 ; 
300   DeltaV = (TgTV - TgTU * TUTV ) / UmTUTV2 ;
301
302   Delta = 1.0 / Sqrt(DeltaU * DeltaU + DeltaV * DeltaV);
303   
304   Tguv1.Multiplied(Delta);  MyTguv1 = Tguv1;
305
306   //------------------------------------------------------------
307   //-- Calcul de Tguv2
308   //--  
309   ThePSurfaceTool::D1(MySurf2,u2,v2,P,TU,TV);
310
311   TUTV = TU.Dot(TV);
312   TgTU = Tg.Dot(TU);
313   TgTV = Tg.Dot(TV);
314   UmTUTV2 = 1.0 - TUTV * TUTV;
315   
316   DeltaU = (TgTU - TgTV * TUTV ) / UmTUTV2 ; 
317   DeltaV = (TgTV - TgTU * TUTV ) / UmTUTV2 ;
318   
319   Delta = 1.0 / Sqrt(DeltaU * DeltaU + DeltaV * DeltaV);
320   
321   Tguv2.Multiplied(Delta);  MyTguv2 = Tguv2;
322
323   return(Standard_True);
324 }
325 #endif
326
327
328
329