0026937: Eliminate NO_CXX_EXCEPTION macro support
[occt.git] / src / Extrema / Extrema_GenExtSS.cxx
1 // Created on: 1995-07-18
2 // Created by: Modelistation
3 // Copyright (c) 1995-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
18 #include <Adaptor3d_Surface.hxx>
19 #include <Extrema_GenExtSS.hxx>
20 #include <Extrema_POnSurf.hxx>
21 #include <math_BFGS.hxx>
22 #include <math_FunctionSetRoot.hxx>
23 #include <math_MultipleVarFunctionWithGradient.hxx>
24 #include <math_Vector.hxx>
25 #include <Standard_OutOfRange.hxx>
26 #include <Standard_TypeMismatch.hxx>
27 #include <StdFail_NotDone.hxx>
28
29 //! This class represents distance objective function for surface / surface.
30 class Extrema_FuncDistSS  : public math_MultipleVarFunctionWithGradient
31 {
32 public:
33     DEFINE_STANDARD_ALLOC
34
35     Standard_EXPORT Extrema_FuncDistSS(const Adaptor3d_Surface& S1,
36                                        const Adaptor3d_Surface& S2)
37     : myS1(&S1),
38       myS2(&S2)
39     {
40     }
41
42   Standard_EXPORT Standard_Integer NbVariables() const
43   {
44     return 4;
45   }
46
47   Standard_EXPORT virtual Standard_Boolean Value(const math_Vector& X,Standard_Real& F)
48   {
49     F = myS1->Value(X(1), X(2)).SquareDistance(myS2->Value(X(3), X(4)));
50     return true;
51   }
52
53   Standard_EXPORT Standard_Boolean Gradient(const math_Vector& X,math_Vector& G)
54   {
55     gp_Pnt P1, P2;
56     gp_Vec Du1s1, Dv1s1;
57     gp_Vec Du2s2, Dv2s2;
58     myS1->D1(X(1),X(2),P1,Du1s1,Dv1s1);
59     myS2->D1(X(3),X(4),P2,Du2s2,Dv2s2);
60
61     gp_Vec P1P2 (P2,P1);
62
63     G(1) = P1P2.Dot(Du1s1);
64     G(2) = P1P2.Dot(Dv1s1);
65     G(3) = -P1P2.Dot(Du2s2);
66     G(4) = -P1P2.Dot(Dv2s2);
67
68     return true;
69   }
70
71   Standard_EXPORT virtual  Standard_Boolean Values(const math_Vector& X,Standard_Real& F,math_Vector& G)
72   {
73     F = myS1->Value(X(1), X(2)).SquareDistance(myS2->Value(X(3), X(4)));
74
75     gp_Pnt P1, P2;
76     gp_Vec Du1s1, Dv1s1;
77     gp_Vec Du2s2, Dv2s2;
78     myS1->D1(X(1),X(2),P1,Du1s1,Dv1s1);
79     myS2->D1(X(3),X(4),P2,Du2s2,Dv2s2);
80
81     gp_Vec P1P2 (P2,P1);
82
83     G(1) = P1P2.Dot(Du1s1);
84     G(2) = P1P2.Dot(Dv1s1);
85     G(3) = -P1P2.Dot(Du2s2);
86     G(4) = -P1P2.Dot(Dv2s2);
87
88     return true;
89   }
90
91 protected:
92
93 private:
94
95   const Adaptor3d_Surface *myS1;
96   const Adaptor3d_Surface *myS2;
97 };
98
99 //=======================================================================
100 //function : Extrema_GenExtSS
101 //purpose  : 
102 //=======================================================================
103 Extrema_GenExtSS::Extrema_GenExtSS()
104 {
105   myDone = Standard_False;
106   myInit = Standard_False;
107 }
108
109 //=======================================================================
110 //function : Extrema_GenExtSS
111 //purpose  : 
112 //=======================================================================
113
114 Extrema_GenExtSS::Extrema_GenExtSS(const Adaptor3d_Surface& S1, 
115                                    const Adaptor3d_Surface& S2, 
116                                    const Standard_Integer NbU, 
117                                    const Standard_Integer NbV, 
118                                    const Standard_Real    Tol1, 
119                                    const Standard_Real    Tol2) : myF(S1,S2)
120 {
121   Initialize(S2, NbU, NbV, Tol2);
122   Perform(S1, Tol1);
123 }
124
125 //=======================================================================
126 //function : Extrema_GenExtSS
127 //purpose  : 
128 //=======================================================================
129
130 Extrema_GenExtSS::Extrema_GenExtSS(const Adaptor3d_Surface& S1, 
131                                    const Adaptor3d_Surface& S2, 
132                                    const Standard_Integer NbU, 
133                                    const Standard_Integer NbV, 
134                                    const Standard_Real U1min, 
135                                    const Standard_Real U1sup, 
136                                    const Standard_Real V1min, 
137                                    const Standard_Real V1sup, 
138                                    const Standard_Real U2min, 
139                                    const Standard_Real U2sup, 
140                                    const Standard_Real V2min, 
141                                    const Standard_Real V2sup, 
142                                    const Standard_Real Tol1, 
143                                    const Standard_Real Tol2): myF(S1, S2)
144 {
145   Initialize(S2, NbU, NbV, U2min,U2sup,V2min,V2sup,Tol2);
146   Perform(S1, U1min,U1sup,V1min,V1sup,Tol1);
147 }
148
149 //=======================================================================
150 //function : Initialize
151 //purpose  : 
152 //=======================================================================
153
154 void Extrema_GenExtSS::Initialize(const Adaptor3d_Surface& S2, 
155                                   const Standard_Integer NbU, 
156                                   const Standard_Integer NbV, 
157                                   const Standard_Real Tol2)
158 {
159   myu2min = S2.FirstUParameter();
160   myu2sup = S2.LastUParameter();
161   myv2min = S2.FirstVParameter();
162   myv2sup = S2.LastVParameter();
163   Initialize(S2,NbU,NbV,myu2min,myu2sup,myv2min,myv2sup,Tol2);
164 }
165
166 //=======================================================================
167 //function : Initialize
168 //purpose  : 
169 //=======================================================================
170
171 void Extrema_GenExtSS::Initialize(const Adaptor3d_Surface& S2, 
172                                   const Standard_Integer NbU, 
173                                   const Standard_Integer NbV, 
174                                   const Standard_Real U2min, 
175                                   const Standard_Real U2sup, 
176                                   const Standard_Real V2min, 
177                                   const Standard_Real V2sup, 
178                                   const Standard_Real Tol2)
179 {
180   myS2 = (Adaptor3d_SurfacePtr)&S2;
181   mypoints1 = new TColgp_HArray2OfPnt(0,NbU+1,0,NbV+1);
182   mypoints2 = new TColgp_HArray2OfPnt(0,NbU+1,0,NbV+1);
183   myusample = NbU;
184   myvsample = NbV;
185   myu2min = U2min;
186   myu2sup = U2sup;
187   myv2min = V2min;
188   myv2sup = V2sup;
189   mytol2 = Tol2;
190
191 // Parametrage de l echantillon sur S2
192
193   Standard_Real PasU = myu2sup - myu2min;
194   Standard_Real PasV = myv2sup - myv2min;
195   Standard_Real U0 = PasU / myusample / 100.;
196   Standard_Real V0 = PasV / myvsample / 100.;
197   gp_Pnt P1;
198   PasU = (PasU - U0) / (myusample - 1);
199   PasV = (PasV - V0) / (myvsample - 1);
200   U0 = myu2min + U0/2.;
201   V0 = myv2min + V0/2.;
202
203 // Calcul des distances
204
205   Standard_Integer NoU, NoV;
206   Standard_Real U, V;
207   for ( NoU = 1, U = U0; NoU <= myusample; NoU++, U += PasU) {
208     for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) {
209       P1 = myS2->Value(U, V);
210       mypoints2->SetValue(NoU,NoV,P1);
211     }
212   }
213 }
214
215 //=======================================================================
216 //function : Perform
217 //purpose  : 
218 //=======================================================================
219
220 void Extrema_GenExtSS::Perform(const Adaptor3d_Surface& S1,
221                                const Standard_Real    Tol1)
222 {
223   myu1min = S1.FirstUParameter();
224   myu1sup = S1.LastUParameter();
225   myv1min = S1.FirstVParameter();
226   myv1sup = S1.LastVParameter();
227   Perform(S1, myu1min, myu1sup,myv1min,myv1sup,Tol1);
228 }
229
230 //=======================================================================
231 //function : Perform
232 //purpose  : 
233 //=======================================================================
234
235 void Extrema_GenExtSS::Perform(const Adaptor3d_Surface& S1,
236                                const Standard_Real U1min, 
237                                const Standard_Real U1sup, 
238                                const Standard_Real V1min, 
239                                const Standard_Real V1sup, 
240                                const Standard_Real Tol1)
241 {
242   myF.Initialize(S1,*myS2);
243   myu1min = U1min;
244   myu1sup = U1sup;
245   myv1min = V1min;
246   myv1sup = V1sup;
247   mytol1 = Tol1;
248
249   Standard_Real U1, V1, U2, V2;
250   Standard_Integer NoU1, NoV1, NoU2, NoV2;
251   gp_Pnt P1, P2;
252
253 // Parametrage de l echantillon sur S1
254
255   Standard_Real PasU1 = myu1sup - myu1min;
256   Standard_Real PasV1 = myv1sup - myv1min;
257   Standard_Real U10 = PasU1 / myusample / 100.;
258   Standard_Real V10 = PasV1 / myvsample / 100.;
259   PasU1 = (PasU1 - U10) / (myusample - 1);
260   PasV1 = (PasV1 - V10) / (myvsample - 1);
261   U10 = myu1min + U10/2.;
262   V10 = myv1min + V10/2.;
263
264   Standard_Real PasU2 = myu2sup - myu2min;
265   Standard_Real PasV2 = myv2sup - myv2min;
266   Standard_Real U20 = PasU2 / myusample / 100.;
267   Standard_Real V20 = PasV2 / myvsample / 100.;
268   PasU2 = (PasU2 - U20) / (myusample - 1);
269   PasV2 = (PasV2 - V20) / (myvsample - 1);
270   U20 = myu2min + U20/2.;
271   V20 = myv2min + V20/2.;
272
273 // Calcul des distances
274
275   for ( NoU1 = 1, U1 = U10; NoU1 <= myusample; NoU1++, U1 += PasU1) {
276     for ( NoV1 = 1, V1 = V10; NoV1 <= myvsample; NoV1++, V1 += PasV1) {
277       P1 = S1.Value(U1, V1);
278       mypoints1->SetValue(NoU1,NoV1,P1);
279     }
280   }
281   
282 /*
283 b- Calcul des minima:
284    -----------------
285    b.a) Initialisations:
286 */
287
288   math_Vector Tol(1, 4);
289   Tol(1) = mytol1;
290   Tol(2) = mytol1;
291   Tol(3) = mytol2;
292   Tol(4) = mytol2;
293   math_Vector UV(1,4), UVinf(1,4), UVsup(1,4);
294   UVinf(1) = myu1min;
295   UVinf(2) = myv1min;
296   UVinf(3) = myu2min;
297   UVinf(4) = myv2min;
298   UVsup(1) = myu1sup;
299   UVsup(2) = myv1sup;
300   UVsup(3) = myu2sup;
301   UVsup(4) = myv2sup;
302   
303
304   Standard_Real distmin = RealLast(), distmax = 0.0, TheDist;
305
306   Standard_Integer N1Umin=0,N1Vmin=0,N2Umin=0,N2Vmin=0;
307   gp_Pnt PP1min, PP2min;
308   Standard_Integer N1Umax=0,N1Vmax=0,N2Umax=0,N2Vmax=0;
309   gp_Pnt PP1max, PP2max;
310
311   for ( NoU1 = 1, U1 = U10; NoU1 <= myusample; NoU1++, U1 += PasU1) {
312     for ( NoV1 = 1, V1 = V10; NoV1 <= myvsample; NoV1++, V1 += PasV1) {
313       P1 = mypoints1->Value(NoU1, NoV1);
314       for ( NoU2 = 1, U2 = U20; NoU2 <= myusample; NoU2++, U2 += PasU2) {
315         for ( NoV2 = 1, V2 = V20; NoV2 <= myvsample; NoV2++, V2 += PasV2) {
316           P2 = mypoints2->Value(NoU2, NoV2);
317           TheDist = P1.SquareDistance(P2);
318           if (TheDist < distmin) {
319             distmin = TheDist;
320             N1Umin = NoU1;
321             N1Vmin = NoV1;
322             N2Umin = NoU2;
323             N2Vmin = NoV2;
324             PP1min = P1;
325             PP2min = P2;
326           }
327           if (TheDist > distmax) {
328             distmax = TheDist;
329             N1Umax = NoU1;
330             N1Vmax = NoV1;
331             N2Umax = NoU2;
332             N2Vmax = NoV2;
333             PP1max = P1;
334             PP2max = P2;
335           }
336         }
337       }
338     }
339   }
340   
341   UV(1) = U10 + (N1Umin - 1) * PasU1;
342   UV(2) = V10 + (N1Vmin - 1) * PasV1;
343   UV(3) = U20 + (N2Umin - 1) * PasU2;
344   UV(4) = V20 + (N2Vmin - 1) * PasV2;
345
346   Extrema_FuncDistSS aGFSS(S1, *myS2);
347   math_BFGS aBFGSSolver(4);
348   aBFGSSolver.Perform(aGFSS, UV);
349   if (aBFGSSolver.IsDone())
350   {
351     aBFGSSolver.Location(UV);
352
353     //  Store result in myF.
354     myF.Value(UV , UV);
355     myF.GetStateNumber();
356   }
357   else
358   {
359     // If optimum is not computed successfully then compute by old approach.
360
361     // Restore initial point.
362     UV(1) = U10 + (N1Umin - 1) * PasU1;
363     UV(2) = V10 + (N1Vmin - 1) * PasV1;
364     UV(3) = U20 + (N2Umin - 1) * PasU2;
365     UV(4) = V20 + (N2Vmin - 1) * PasV2;
366
367     math_FunctionSetRoot SR1(myF, Tol);
368     SR1.Perform(myF, UV, UVinf, UVsup);
369   }
370
371   //math_FunctionSetRoot SR1(myF, Tol);
372   //SR1.Perform(myF, UV, UVinf, UVsup);
373
374   UV(1) = U10 + (N1Umax - 1) * PasU1;
375   UV(2) = V10 + (N1Vmax - 1) * PasV1;
376   UV(3) = U20 + (N2Umax - 1) * PasU2;
377   UV(4) = V20 + (N2Vmax - 1) * PasV2;
378
379   // It is impossible to compute max distance in the same manner,
380   // since for the distance functional for max have bad definition.
381   // So, for max computation old approach is used.
382   math_FunctionSetRoot SR2(myF, Tol);
383   SR2.Perform(myF, UV, UVinf, UVsup);
384
385   myDone = Standard_True;
386 }
387
388 //=======================================================================
389 //function : IsDone
390 //purpose  : 
391 //=======================================================================
392
393 Standard_Boolean Extrema_GenExtSS::IsDone() const 
394 {
395   return myDone;
396 }
397
398 //=======================================================================
399 //function : NbExt
400 //purpose  : 
401 //=======================================================================
402
403 Standard_Integer Extrema_GenExtSS::NbExt() const 
404 {
405   if (!IsDone()) { throw StdFail_NotDone(); }
406   return myF.NbExt();
407
408 }
409
410 //=======================================================================
411 //function : Value
412 //purpose  : 
413 //=======================================================================
414
415 Standard_Real Extrema_GenExtSS::SquareDistance(const Standard_Integer N) const 
416 {  
417   if (!IsDone()) { throw StdFail_NotDone(); }
418   return myF.SquareDistance(N);
419 }
420
421 //=======================================================================
422 //function : PointOnS1
423 //purpose  : 
424 //=======================================================================
425
426 const Extrema_POnSurf& Extrema_GenExtSS::PointOnS1(const Standard_Integer N) const 
427 {
428   if (!IsDone()) { throw StdFail_NotDone(); }
429   return myF.PointOnS1(N);
430 }
431
432 //=======================================================================
433 //function : PointOnS2
434 //purpose  : 
435 //=======================================================================
436
437 const Extrema_POnSurf& Extrema_GenExtSS::PointOnS2(const Standard_Integer N) const 
438 {
439   if (!IsDone()) { throw StdFail_NotDone(); }
440   return myF.PointOnS2(N);
441 }
442
443 //=======================================================================
444 //function : Bidon
445 //purpose  : 
446 //=======================================================================
447
448 Adaptor3d_SurfacePtr Extrema_GenExtSS::Bidon() const 
449 {
450   return (Adaptor3d_SurfacePtr)0L;
451 }
452