0023830: BRepExtrema_DistShapeShape does not find intersection of face with edge
[occt.git] / src / Extrema / Extrema_GenExtCS.cxx
1 // Created on: 1995-07-18
2 // Created by: Modelistation
3 // Copyright (c) 1995-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
23 #include <Extrema_GenExtCS.ixx>
24
25 #include <math_Vector.hxx>
26 #include <math_FunctionSetRoot.hxx>
27 #include <Precision.hxx>
28
29 #include <Geom_Line.hxx>
30 #include <Adaptor3d_HCurve.hxx>
31 #include <GeomAdaptor_Curve.hxx>
32
33 #include <Extrema_ExtCC.hxx>
34 #include <Extrema_POnCurv.hxx>
35 #include <Extrema_ExtPS.hxx>
36
37 //=======================================================================
38 //function : Extrema_GenExtCS
39 //purpose  : 
40 //=======================================================================
41
42 Extrema_GenExtCS::Extrema_GenExtCS()
43 {
44   myDone = Standard_False;
45   myInit = Standard_False;
46 }
47
48 //=======================================================================
49 //function : Extrema_GenExtCS
50 //purpose  : 
51 //=======================================================================
52
53  Extrema_GenExtCS::Extrema_GenExtCS(const Adaptor3d_Curve& C, 
54                                     const Adaptor3d_Surface& S, 
55                                     const Standard_Integer NbT, 
56                                     const Standard_Integer NbU, 
57                                     const Standard_Integer NbV, 
58                                     const Standard_Real Tol1, 
59                                     const Standard_Real Tol2)
60 {
61   Initialize(S, NbU, NbV, Tol2);
62   Perform(C, NbT, Tol1);
63 }
64
65 //=======================================================================
66 //function : Extrema_GenExtCS
67 //purpose  : 
68 //=======================================================================
69
70  Extrema_GenExtCS::Extrema_GenExtCS(const Adaptor3d_Curve& C, 
71                                     const Adaptor3d_Surface& S, 
72                                     const Standard_Integer NbT, 
73                                     const Standard_Integer NbU, 
74                                     const Standard_Integer NbV, 
75                                     const Standard_Real tmin, 
76                                     const Standard_Real tsup, 
77                                     const Standard_Real Umin, 
78                                     const Standard_Real Usup,
79                                     const Standard_Real Vmin, 
80                                     const Standard_Real Vsup, 
81                                     const Standard_Real Tol1, 
82                                     const Standard_Real Tol2)
83 {
84   Initialize(S, NbU, NbV, Umin,Usup,Vmin,Vsup,Tol2);
85   Perform(C, NbT, tmin, tsup, Tol1);
86 }
87
88 //=======================================================================
89 //function : Initialize
90 //purpose  : 
91 //=======================================================================
92
93 void Extrema_GenExtCS::Initialize(const Adaptor3d_Surface& S, 
94                                   const Standard_Integer NbU, 
95                                   const Standard_Integer NbV, 
96                                   const Standard_Real Tol2)
97 {
98   myumin = S.FirstUParameter();
99   myusup = S.LastUParameter();
100   myvmin = S.FirstVParameter();
101   myvsup = S.LastVParameter();
102   Initialize(S,NbU,NbV,myumin,myusup,myvmin,myvsup,Tol2);
103 }
104
105 //=======================================================================
106 //function : Initialize
107 //purpose  : 
108 //=======================================================================
109
110 void Extrema_GenExtCS::Initialize(const Adaptor3d_Surface& S, 
111                                   const Standard_Integer NbU, 
112                                   const Standard_Integer NbV, 
113                                   const Standard_Real Umin, 
114                                   const Standard_Real Usup, 
115                                   const Standard_Real Vmin, 
116                                   const Standard_Real Vsup, 
117                                   const Standard_Real Tol2)
118 {
119   myS = (Adaptor3d_SurfacePtr)&S;
120   myusample = NbU;
121   myvsample = NbV;
122   myumin = Umin;
123   myusup = Usup;
124   myvmin = Vmin;
125   myvsup = Vsup;
126   mytol2 = Tol2;
127 }
128
129 //=======================================================================
130 //function : Perform
131 //purpose  : 
132 //=======================================================================
133
134 void Extrema_GenExtCS::Perform(const Adaptor3d_Curve& C, 
135                                const Standard_Integer NbT,
136                                const Standard_Real Tol1)
137 {
138   mytmin = C.FirstParameter();
139   mytsup = C.LastParameter();
140   Perform(C, NbT, mytmin, mytsup,Tol1);
141 }
142
143 //=======================================================================
144 //function : Perform
145 //purpose  : 
146 //=======================================================================
147
148 void Extrema_GenExtCS::Perform(const Adaptor3d_Curve& C, 
149                                const Standard_Integer NbT,
150                                const Standard_Real tmin, 
151                                const Standard_Real tsup, 
152                                const Standard_Real Tol1)
153 {
154   myDone = Standard_False;
155   myF.Initialize(C,*myS);
156   mytmin = tmin;
157   mytsup = tsup;
158   mytol1 = Tol1;
159   mytsample = NbT;
160   // Modif de lvt pour trimer la surface non pas aux infinis mais  a +/- 10000
161
162   Standard_Real trimusup = myusup, trimumin = myumin,trimvsup = myvsup,trimvmin = myvmin;
163   if (Precision::IsInfinite(trimusup)){
164     trimusup = 1.0e+10;
165   }
166   if (Precision::IsInfinite(trimvsup)){
167     trimvsup = 1.0e+10;
168   }
169   if (Precision::IsInfinite(trimumin)){
170     trimumin = -1.0e+10;
171   }
172   if (Precision::IsInfinite(trimvmin)){
173     trimvmin = -1.0e+10;
174   }
175   //
176   math_Vector Tol(1, 3);
177   Tol(1) = mytol1;
178   Tol(2) = mytol2;
179   Tol(3) = mytol2;
180   math_Vector UV(1,3), UVinf(1,3), UVsup(1,3);
181   UVinf(1) = mytmin;
182   UVinf(2) = trimumin;
183   UVinf(3) = trimvmin;
184
185   UVsup(1) = mytsup;
186   UVsup(2) = trimusup;
187   UVsup(3) = trimvsup;
188   // 18/02/02 akm vvv : (OCC163) bad extremas - extrusion surface versus the line.
189   //                    Try to preset the initial solution as extrema between
190   //                    extrusion direction and the curve.
191   if (myS->GetType() == GeomAbs_SurfaceOfExtrusion)
192   {
193     gp_Dir aDir = myS->Direction();
194     Handle(Adaptor3d_HCurve) aCurve = myS->BasisCurve();
195     Standard_Real dfUFirst = aCurve->FirstParameter();
196     // Create iso line of U=U0
197     GeomAdaptor_Curve anAx(new Geom_Line(aCurve->Value(dfUFirst), aDir),
198                            trimvmin, trimvsup);
199     Extrema_ExtCC aLocator(C, anAx);
200     if (aLocator.IsDone() && aLocator.NbExt()>0)
201     {
202       Standard_Integer iExt;
203       // Try to find all extremas
204       Extrema_POnCurv aP1, aP2;
205       for (iExt=1; iExt<=aLocator.NbExt(); iExt++)
206       {
207   aLocator.Points (iExt, aP1, aP2);
208   // Parameter on curve
209   UV(1) = aP1.Parameter();
210   // To find parameters on surf, try ExtPS
211   Extrema_ExtPS aPreciser (aP1.Value(), *myS, mytol2, mytol2);
212   if (aPreciser.IsDone())
213   {
214     // Managed to find extremas between point and surface
215     Standard_Integer iPExt;
216     for (iPExt=1; iPExt<=aPreciser.NbExt(); iPExt++)
217     {
218       aPreciser.Point(iPExt).Parameter(UV(2),UV(3));
219       math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
220     }
221   }
222   else
223   {
224     // Failed... try the point on iso line
225     UV(2) = dfUFirst;
226     UV(3) = aP2.Parameter();
227     math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
228   }
229       } // for (iExt=1; iExt<=aLocator.NbExt(); iExt++)
230     } // if (aLocator.IsDone() && aLocator.NbExt()>0)
231   } // if (myS.Type() == GeomAbs_ExtrusionSurface)
232   else
233   {
234     Standard_Real aCUAdd = (mytsup - mytmin) / mytsample;
235     Standard_Real aSUAdd = (myusup - myumin) / myusample;
236     Standard_Real aSVAdd = (myvsup - myvmin) / myvsample;
237     TColgp_HArray1OfPnt aCPs(1, mytsample);
238     TColgp_HArray2OfPnt aSPs(1, myusample, 1, myvsample);
239     Standard_Integer aRestIterCount = 3;
240       // The value is calculated by the bug CR23830.
241     Standard_Integer aCUDen = 2, aSUDen = 2, aSVDen = 2;
242     Standard_Boolean anAreAvSqsInited = Standard_False;
243     Standard_Real aCUSq = 0, aSUSq = 0, aSVSq = 0;
244     while (aRestIterCount--)
245     {
246       Standard_Real aMinCU, aMinSU, aMinSV, aMaxCU, aMaxSU, aMaxSV;
247       Standard_Real aMinSqDist = DBL_MAX, aMaxSqDist = -DBL_MAX;
248       for (Standard_Integer aSUNom = 1; aSUNom < aSUDen; aSUNom += 2)
249       {
250         Standard_Real aSU0 = myumin + (aSUNom * aSUAdd) / aSUDen;
251         for (Standard_Integer aSVNom = 1; aSVNom < aSVDen; aSVNom += 2)
252         {
253           Standard_Real aSV0 = myvmin + (aSVNom * aSVAdd) / aSVDen;
254           for (Standard_Integer aCUNom = 1; aCUNom < aCUDen; aCUNom += 2)
255           {
256             Standard_Real aCU0 = mytmin + (aCUNom * aCUAdd) / aCUDen;
257             Standard_Real aCU = aCU0;
258             for (Standard_Integer aCUI = 1; aCUI <= mytsample;
259               aCUI++, aCU += aCUAdd)
260             {
261               aCPs.SetValue(aCUI, C.Value(aCU));
262             }
263             //
264             aCU = aCU0;
265             Standard_Real aSU = aSU0;
266             for (Standard_Integer aSUI = 1; aSUI <= myusample;
267               aSUI++, aSU += aSUAdd)
268             {
269               Standard_Real aSV = aSV0;
270               for (Standard_Integer aSVI = 1; aSVI <= myvsample;
271                 aSVI++, aSV += aSVAdd)
272               {
273                 gp_Pnt aSP = myS->Value(aSU, aSV);
274                 aSPs.ChangeValue(aSUI, aSVI) = aSP;
275                 Standard_Real aCU = aCU0;
276                 for (Standard_Integer aCUI = 1; aCUI <= mytsample;
277                   aCUI++, aCU += aCUAdd)
278                 {
279                   Standard_Real aSqDist = aSP.SquareDistance(aCPs.Value(aCUI));
280                   if (aSqDist < aMinSqDist)
281                   {
282                     aMinSqDist = aSqDist;
283                     aMinCU = aCU;
284                     aMinSU = aSU;
285                     aMinSV = aSV;
286                   }
287                   if (aMaxSqDist < aSqDist)
288                   {
289                     aMaxSqDist = aSqDist;
290                     aMaxCU = aCU;
291                     aMaxSU = aSU;
292                     aMaxSV = aSV;
293                   }
294                 }
295               }
296             }
297           }
298         }
299       }
300       // Find min approximation.
301       UV(1) = aMinCU;
302       UV(2) = aMinSU;
303       UV(3) = aMinSV;
304       math_FunctionSetRoot anA(myF, UV, Tol, UVinf, UVsup, 100, aRestIterCount);
305       // Find max approximation.
306       if (!anA.IsDivergent() || !aRestIterCount)
307       {
308         UV(1) = aMaxCU;
309         UV(2) = aMaxSU;
310         UV(3) = aMaxSV;
311         math_FunctionSetRoot(myF, UV, Tol, UVinf, UVsup);
312         break;
313       }
314       //
315       if (!anAreAvSqsInited)
316       {
317         anAreAvSqsInited = Standard_True;
318         //
319         const gp_Pnt * aCP1 = &aCPs.Value(1);
320         for (Standard_Integer aCUI = 2; aCUI <= mytsample; aCUI++)
321         {
322           const gp_Pnt & aCP2 = aCPs.Value(aCUI);
323           aCUSq += aCP1->SquareDistance(aCP2);
324           aCP1 = &aCP2;
325         }
326         aCUSq /= mytsample - 1;
327         //
328         for (Standard_Integer aSVI = 1; aSVI <= myvsample; aSVI++)
329         {
330           const gp_Pnt * aSP1 = &aSPs.Value(1, aSVI);
331           for (Standard_Integer aSUI = 2; aSUI <= myusample; aSUI++)
332           {
333             const gp_Pnt & aSP2 = aSPs.Value(aSUI, aSVI);
334             aSUSq += aSP1->SquareDistance(aSP2);
335             aSP1 = &aSP2;
336           }
337         }
338         aSUSq /= myvsample * (myusample - 1);
339         //
340         for (Standard_Integer aSUI = 1; aSUI <= myusample; aSUI++)
341         {
342           const gp_Pnt * aSP1 = &aSPs.Value(aSUI, 1);
343           for (Standard_Integer aSVI = 2; aSVI <= myvsample; aSVI++)
344           {
345             const gp_Pnt & aSP2 = aSPs.Value(aSUI, aSVI);
346             aSVSq += aSP1->SquareDistance(aSP2);
347             aSP1 = &aSP2;
348           }
349         }
350         aSVSq /= myusample * (myvsample - 1);
351       }
352       //
353       if ((aSUSq <= aCUSq) && (aSVSq <= aCUSq))
354       {
355         aCUDen += aCUDen;
356         aCUSq /= 4;
357       }
358       else if ((aCUSq <= aSUSq) && (aSVSq <= aSUSq))
359       {
360         aSUDen += aSUDen;
361         aSUSq /= 4;
362       }
363       else
364       {
365         aSVDen += aSVDen;
366         aSVSq /= 4;
367       }
368     }
369   }
370
371   myDone = Standard_True;
372 }
373
374 //=======================================================================
375 //function : IsDone
376 //purpose  : 
377 //=======================================================================
378
379 Standard_Boolean Extrema_GenExtCS::IsDone() const 
380 {
381   return myDone;
382 }
383
384 //=======================================================================
385 //function : NbExt
386 //purpose  : 
387 //=======================================================================
388
389 Standard_Integer Extrema_GenExtCS::NbExt() const 
390 {
391   if (!IsDone()) { StdFail_NotDone::Raise(); }
392   return myF.NbExt();
393 }
394
395 //=======================================================================
396 //function : Value
397 //purpose  : 
398 //=======================================================================
399
400 Standard_Real Extrema_GenExtCS::SquareDistance(const Standard_Integer N) const 
401 {
402   if (!IsDone()) { StdFail_NotDone::Raise(); }
403   return myF.SquareDistance(N);
404 }
405
406 //=======================================================================
407 //function : PointOnCurve
408 //purpose  : 
409 //=======================================================================
410
411 const Extrema_POnCurv& Extrema_GenExtCS::PointOnCurve(const Standard_Integer N) const 
412 {
413   if (!IsDone()) { StdFail_NotDone::Raise(); }
414   return myF.PointOnCurve(N);
415 }
416
417 //=======================================================================
418 //function : PointOnSurface
419 //purpose  : 
420 //=======================================================================
421
422 const Extrema_POnSurf& Extrema_GenExtCS::PointOnSurface(const Standard_Integer N) const 
423 {
424   if (!IsDone()) { StdFail_NotDone::Raise(); }
425   return myF.PointOnSurface(N);
426 }
427
428 //=======================================================================
429 //function : BidonSurface
430 //purpose  : 
431 //=======================================================================
432
433 Adaptor3d_SurfacePtr Extrema_GenExtCS::BidonSurface() const 
434 {
435   return (Adaptor3d_SurfacePtr)0L;
436 }
437
438 //=======================================================================
439 //function : BidonCurve
440 //purpose  : 
441 //=======================================================================
442
443 Adaptor3d_CurvePtr Extrema_GenExtCS::BidonCurve() const 
444 {
445   return (Adaptor3d_CurvePtr)0L;
446 }
447