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