0024002: Overall code and build procedure refactoring -- automatic
[occt.git] / src / BRepBlend / BRepBlend_AppFuncRoot.cxx
1 // Created on: 1998-05-12
2 // Created by: Philippe NOUAILLE
3 // Copyright (c) 1998-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 <Blend_AppFunction.hxx>
19 #include <Blend_Point.hxx>
20 #include <BRepBlend_AppFuncRoot.hxx>
21 #include <BRepBlend_Line.hxx>
22 #include <gp_Pnt.hxx>
23 #include <math_FunctionSetRoot.hxx>
24 #include <Standard_OutOfRange.hxx>
25 #include <Standard_Type.hxx>
26 #include <TColgp_HArray1OfPnt.hxx>
27 #include <TColgp_HArray1OfPnt2d.hxx>
28 #include <TColgp_HArray1OfVec.hxx>
29 #include <TColgp_HArray1OfVec2d.hxx>
30 #include <TColStd_HArray1OfInteger.hxx>
31 #include <TColStd_HArray1OfReal.hxx>
32
33 BRepBlend_AppFuncRoot::BRepBlend_AppFuncRoot(Handle(BRepBlend_Line)& Line,
34                                              Blend_AppFunction&      Func,
35                                              const Standard_Real     Tol3d,
36                                              const Standard_Real     Tol2d)
37 :myLine(Line),
38  myFunc(&Func),
39  myTolerance(1,Func.NbVariables()),
40  X1(1,Func.NbVariables()),
41  X2(1,Func.NbVariables()),
42  XInit(1,Func.NbVariables()),
43  Sol(1,Func.NbVariables())
44 {
45   Standard_Integer NbPoles, NbKnots, Degree, NbPoles2d;
46   Standard_Integer ii;
47   
48   //  Tolerances    
49   Func.GetTolerance(myTolerance, Tol3d);
50   Standard_Integer dim = Func.NbVariables();
51   for (ii=1; ii<= dim; ii++) {
52     if (myTolerance(ii)>Tol2d) { myTolerance(ii) = Tol2d;}
53   }
54   
55   //  Tables
56   Func.GetShape( NbPoles, NbKnots, Degree, NbPoles2d);
57   
58   // Calculation of BaryCentre (rationnal case).
59   if (Func.IsRational()) {
60     Standard_Real Xmax =-1.e100, Xmin = 1.e100, 
61     Ymax =-1.e100, Ymin = 1.e100, 
62     Zmax =-1.e100, Zmin = 1.e100;
63     Blend_Point P;
64     for (ii=1; ii<=myLine->NbPoints(); ii++) {
65       P = myLine->Point(ii);
66       Xmax = Max ( Max(P.PointOnS1().X(), P.PointOnS2().X()), Xmax);
67       Xmin = Min ( Min(P.PointOnS1().X(), P.PointOnS2().X()), Xmin);
68       Ymax = Max ( Max(P.PointOnS1().Y(), P.PointOnS2().Y()), Ymax);
69       Ymin = Min ( Min(P.PointOnS1().Y(), P.PointOnS2().Y()), Ymin);
70       Zmax = Max ( Max(P.PointOnS1().Z(), P.PointOnS2().Z()), Zmax);
71       Zmin = Min ( Min(P.PointOnS1().Z(), P.PointOnS2().Z()), Zmin);
72       
73       myBary.SetCoord((Xmax+Xmin)/2, (Ymax+Ymin)/2, (Zmax+Zmin)/2);
74     }
75   }
76   else {myBary.SetCoord(0,0,0);}
77 }
78
79 //================================================================================ 
80 // Function: D0
81 // Purpose : Calculation of section for v = Param, if calculation fails
82 //           Standard_False is raised. 
83 //================================================================================
84 Standard_Boolean BRepBlend_AppFuncRoot::D0(const Standard_Real Param,
85                                            const Standard_Real /*First*/,
86                                            const Standard_Real /*Last*/,
87                                                  TColgp_Array1OfPnt& Poles,
88                                                  TColgp_Array1OfPnt2d& Poles2d,
89                                                  TColStd_Array1OfReal& Weigths) 
90 {
91   Standard_Boolean Ok=Standard_True;
92   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
93   Ok = SearchPoint( *Func, Param, myPnt);
94   
95   if (Ok) (*Func).Section(myPnt,
96                           Poles,
97                           Poles2d,
98                           Weigths);
99   return Ok;
100 }
101
102 //================================================================================ 
103 // Function: D1
104 // Purpose : Calculation of the partial derivative of the section corresponding to v
105 //           for v = Param, if the calculation fails Standard_False is raised.
106 //================================================================================ 
107 Standard_Boolean BRepBlend_AppFuncRoot::D1(const Standard_Real Param,
108                                            const Standard_Real /*First*/,
109                                            const Standard_Real /*Last*/,
110                                            TColgp_Array1OfPnt& Poles,
111                                            TColgp_Array1OfVec& DPoles,
112                                            TColgp_Array1OfPnt2d& Poles2d,
113                                            TColgp_Array1OfVec2d& DPoles2d,
114                                            TColStd_Array1OfReal& Weigths,
115                                            TColStd_Array1OfReal& DWeigths) 
116 {
117   Standard_Boolean Ok=Standard_True;
118   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
119   
120   Ok = SearchPoint( *Func, Param, myPnt);
121   
122   if (Ok) {
123     Ok = (*Func).Section(myPnt,
124                          Poles, DPoles,
125                          Poles2d, DPoles2d,
126                          Weigths, DWeigths);
127   }
128   
129   return Ok;
130 }
131
132 //=========================================================================== 
133 // Function: D2
134 // Purpose : Calculation of the derivative and second partial of the 
135 //           section corresponding to v.
136 //           For v = Param, if the calculation fails Standard_False is raised.  
137 //=========================================================================== 
138 Standard_Boolean BRepBlend_AppFuncRoot::D2(const Standard_Real Param,
139                                            const Standard_Real /*First*/,
140                                            const Standard_Real /*Last*/,
141                                            TColgp_Array1OfPnt& Poles,
142                                            TColgp_Array1OfVec& DPoles,
143                                            TColgp_Array1OfVec& D2Poles,
144                                            TColgp_Array1OfPnt2d& Poles2d,
145                                            TColgp_Array1OfVec2d& DPoles2d,
146                                            TColgp_Array1OfVec2d& D2Poles2d,
147                                            TColStd_Array1OfReal& Weigths,
148                                            TColStd_Array1OfReal& DWeigths,
149                                            TColStd_Array1OfReal& D2Weigths) 
150 {
151   Standard_Boolean Ok=Standard_True;
152   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
153   
154   Ok = SearchPoint( *Func, Param, myPnt); 
155   if (Ok) {
156     Ok = (*Func).Section(myPnt,
157                          Poles, DPoles, D2Poles,
158                          Poles2d, DPoles2d, D2Poles2d,
159                          Weigths, DWeigths, D2Weigths);
160   }
161   return Ok;
162 }
163
164 Standard_Integer BRepBlend_AppFuncRoot::Nb2dCurves() const
165 {
166   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
167   Standard_Integer i,j,k,nbpol2d;
168   (*Func).GetShape(i,j,k,nbpol2d);
169   return nbpol2d;
170 }
171
172 void BRepBlend_AppFuncRoot::SectionShape(Standard_Integer& NbPoles,
173                                          Standard_Integer& NbKnots,
174                                          Standard_Integer& Degree) const
175 {
176   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
177   Standard_Integer ii;
178   (*Func).GetShape( NbPoles, NbKnots, Degree, ii);
179 }
180
181 void BRepBlend_AppFuncRoot::Knots(TColStd_Array1OfReal& TKnots) const
182 {
183   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc; 
184   Func->Knots(TKnots);
185 }
186
187 void BRepBlend_AppFuncRoot::Mults(TColStd_Array1OfInteger& TMults) const
188 {  
189   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
190   Func->Mults(TMults);
191 }
192
193 Standard_Boolean BRepBlend_AppFuncRoot::IsRational() const
194 {
195   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
196   return  (*Func).IsRational();
197 }
198
199 Standard_Integer BRepBlend_AppFuncRoot::NbIntervals(const GeomAbs_Shape S) const
200 {
201   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
202   return  Func->NbIntervals(S);
203 }
204
205 void BRepBlend_AppFuncRoot::Intervals(TColStd_Array1OfReal& T,const GeomAbs_Shape S) const
206 {
207   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
208   Func->Intervals(T, S);
209 }
210
211 void BRepBlend_AppFuncRoot::SetInterval(const Standard_Real First,const Standard_Real Last) 
212 {
213   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
214   Func->Set(First, Last);
215 }
216
217 void BRepBlend_AppFuncRoot::Resolution(const Standard_Integer Index,
218                                        const Standard_Real Tol,
219                                        Standard_Real& TolU,
220                                        Standard_Real& TolV) const
221 {
222   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
223   Func->Resolution(Index,Tol,TolU,TolV);
224 }
225
226 void BRepBlend_AppFuncRoot::GetTolerance(const Standard_Real BoundTol,
227                                          const Standard_Real SurfTol,
228                                          const Standard_Real AngleTol,
229                                          TColStd_Array1OfReal& Tol3d) const
230 {
231   Standard_Integer ii;
232   math_Vector V3d(1, Tol3d.Length()), V1d(1, Tol3d.Length());
233   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc; 
234   
235   Func->GetTolerance(BoundTol, SurfTol, AngleTol, V3d, V1d); 
236   for (ii=1; ii<=Tol3d.Length(); ii++) Tol3d(ii) = V3d(ii);
237 }
238
239 void BRepBlend_AppFuncRoot::SetTolerance(const Standard_Real Tol3d, 
240                                          const Standard_Real Tol2d)
241 {
242   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc; 
243   Standard_Integer ii, dim = Func->NbVariables();
244   Func->GetTolerance(myTolerance, Tol3d);
245   for (ii=1; ii<=dim; ii++) {
246     if (myTolerance(ii)>Tol2d) { myTolerance(ii) = Tol2d;}
247   }
248
249
250 gp_Pnt BRepBlend_AppFuncRoot::BarycentreOfSurf() const
251 {
252   return myBary;
253 }
254
255 Standard_Real BRepBlend_AppFuncRoot::MaximalSection() const
256 {
257   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc; 
258   return Func->GetSectionSize(); 
259 }
260
261 void BRepBlend_AppFuncRoot::GetMinimalWeight(TColStd_Array1OfReal& Weigths) const
262 {
263   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
264   Func->GetMinimalWeight(Weigths);
265 }
266
267
268 //================================================================================ 
269 //
270 // Function : SearchPoint
271 //
272 // Purpose : Find point solution with parameter Param (on 2 Surfaces)
273 //
274 // Algorithm : 
275 //     1) Approximative solution is found from already calculated Points
276 //     2) Convergence is done by a method of type Newton
277 // 
278 // Possible causes of fails : 
279 //        - Singularity on surfaces.
280 //        - no information oin the "line" resulting from processing. 
281 //            
282 //================================================================================  
283
284 Standard_Boolean BRepBlend_AppFuncRoot::SearchPoint(Blend_AppFunction& Func,
285                                                     const Standard_Real Param,
286                                                     Blend_Point& Pnt)
287 {
288   Standard_Boolean Trouve;
289   Standard_Integer dim = Func.NbVariables();
290   // (1) Find a point of init
291   Standard_Integer I1=1, I2=myLine->NbPoints(), Index;
292   Standard_Real t1, t2;
293   
294   //  (1.a) It is checked if it is inside
295   if (Param < myLine->Point(I1).Parameter()) {return Standard_False;}
296   if (Param > myLine->Point(I2).Parameter()) {return Standard_False;}
297   
298   //  (1.b) Find the interval
299   Trouve = SearchLocation(Param, I1, I2, Index);
300   
301   //  (1.c) If the point is already calculated it is returned
302   if (Trouve) {
303     Pnt = myLine->Point(Index);
304     Vec(XInit,Pnt);
305   }
306   else {
307     //  (1.d) Intialisation by linear interpolation
308     Pnt = myLine->Point(Index);
309     Vec(X1,Pnt);
310     t1 = Pnt.Parameter();
311
312     Pnt = myLine->Point(Index+1);
313     Vec(X2,Pnt);
314     t2 = Pnt.Parameter();
315
316     Standard_Real Parammt1 = (Param-t1) / (t2-t1);
317     Standard_Real t2mParam = (t2-Param) / (t2-t1);
318     for(Standard_Integer i = 1; i <= dim; i++){
319       XInit(i) = X2(i) * Parammt1 + X1(i) * t2mParam;
320     }
321   }
322
323   // (2) Calculation of the solution ------------------------
324   Func.Set(Param);
325   Func.GetBounds(X1, X2);
326   math_FunctionSetRoot rsnld(Func, myTolerance, 30);
327   
328   rsnld.Perform(Func, XInit, X1, X2);
329   
330   if (!rsnld.IsDone()) {
331 # ifdef BREPBLEND_DEB
332     cout << "AppFunc : RNLD Not done en t = " <<  Param  << endl;
333 # endif 
334     return Standard_False;
335   }
336   rsnld.Root(Sol);
337   
338   // (3) Storage of the point
339   Point(Func,Param,Sol,Pnt);
340
341   // (4) Insertion of the point if the calculation seems long.
342   if ((!Trouve)&&(rsnld.NbIterations()>3)) {
343 #ifdef OCCT_DEBUG
344     cout << "Evaluation in t = " <<  Param << "given" << endl;
345     rsnld.Dump(cout);
346 #endif
347     myLine->InsertBefore(Index+1, Pnt);
348   }
349   return Standard_True;
350 }
351
352
353 //=============================================================================
354 //
355 // Function : SearchLocation
356 //
357 // Purpose : Binary search of the line of the parametric interval containing
358 //           Param in the list of calculated points (myline)
359 //           if the point of parameter Param is already stored in the list
360 //           True is raised and ParamIndex corresponds to line of Point.
361 //           Complexity of this algorithm is log(n)/log(2)
362 //================================================================================ 
363 Standard_Boolean BRepBlend_AppFuncRoot::SearchLocation(const Standard_Real Param,
364                                                        const Standard_Integer FirstIndex,
365                                                        const Standard_Integer LastIndex,
366                                                        Standard_Integer& ParamIndex) const
367 {
368   Standard_Integer Ideb = FirstIndex, Ifin =  LastIndex, Idemi;
369   Standard_Real Valeur;
370   
371   Valeur = myLine->Point(Ideb).Parameter();
372   if (Param == Valeur) {
373     ParamIndex = Ideb;
374     return Standard_True; 
375   }
376   
377   Valeur = myLine->Point(Ifin).Parameter();
378   if (Param == Valeur) {
379     ParamIndex = Ifin;
380     return Standard_True; 
381   } 
382   
383   while ( Ideb+1 != Ifin) {
384     Idemi = (Ideb+Ifin)/2;
385     Valeur = myLine->Point(Idemi).Parameter();
386     if (Valeur < Param) {Ideb = Idemi;}
387     else { 
388       if ( Valeur > Param) { Ifin = Idemi;}
389       else                 { ParamIndex = Idemi;
390                              return Standard_True;}
391     }
392   }
393   
394   ParamIndex = Ideb;
395   return Standard_False;
396 }
397