0024428: Implementation of LGPL license
[occt.git] / src / BRepBlend / BRepBlend_AppFuncRoot.cxx
CommitLineData
b311480e 1// Created on: 1998-05-12
2// Created by: Philippe NOUAILLE
3// Copyright (c) 1998-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
973c2be1 8// This library is free software; you can redistribute it and / or modify it
9// under the terms of the GNU Lesser General Public 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
17#include <BRepBlend_AppFuncRoot.ixx>
18#include <Blend_AppFunction.hxx>
19
20#include <Blend_Point.hxx>
21#include <BRepBlend_Line.hxx>
22
23#include <math_FunctionSetRoot.hxx>
24
25#include <TColgp_HArray1OfPnt.hxx>
26#include <TColgp_HArray1OfPnt2d.hxx>
27#include <TColgp_HArray1OfVec.hxx>
28#include <TColgp_HArray1OfVec2d.hxx>
29
30#include <TColStd_HArray1OfReal.hxx>
31#include <TColStd_HArray1OfInteger.hxx>
32
33BRepBlend_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
81bba717 55 // Tables
7fd59977 56 Func.GetShape( NbPoles, NbKnots, Degree, NbPoles2d);
57
81bba717 58 // Calculation of BaryCentre (rationnal case).
7fd59977 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
81bba717 81// Purpose : Calculation of section for v = Param, if calculation fails
82// Standard_False is raised.
7fd59977 83//================================================================================
84Standard_Boolean BRepBlend_AppFuncRoot::D0(const Standard_Real Param,
35e08fe8 85 const Standard_Real /*First*/,
86 const Standard_Real /*Last*/,
87 TColgp_Array1OfPnt& Poles,
88 TColgp_Array1OfPnt2d& Poles2d,
89 TColStd_Array1OfReal& Weigths)
7fd59977 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,
35e08fe8 96 Poles,
97 Poles2d,
98 Weigths);
7fd59977 99 return Ok;
100}
101
102//================================================================================
103// Function: D1
81bba717 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.
7fd59977 106//================================================================================
107Standard_Boolean BRepBlend_AppFuncRoot::D1(const Standard_Real Param,
35e08fe8 108 const Standard_Real /*First*/,
109 const Standard_Real /*Last*/,
7fd59977 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
81bba717 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.
7fd59977 137//===========================================================================
138Standard_Boolean BRepBlend_AppFuncRoot::D2(const Standard_Real Param,
35e08fe8 139 const Standard_Real /*First*/,
140 const Standard_Real /*Last*/,
7fd59977 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
164Standard_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
172void 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
181void BRepBlend_AppFuncRoot::Knots(TColStd_Array1OfReal& TKnots) const
182{
183 Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
184 Func->Knots(TKnots);
185}
186
187void BRepBlend_AppFuncRoot::Mults(TColStd_Array1OfInteger& TMults) const
188{
189 Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
190 Func->Mults(TMults);
191}
192
193Standard_Boolean BRepBlend_AppFuncRoot::IsRational() const
194{
195 Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
196 return (*Func).IsRational();
197}
198
199Standard_Integer BRepBlend_AppFuncRoot::NbIntervals(const GeomAbs_Shape S) const
200{
201 Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
202 return Func->NbIntervals(S);
203}
204
205void 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
211void 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
217void 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
226void 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
239void 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
250gp_Pnt BRepBlend_AppFuncRoot::BarycentreOfSurf() const
251{
252 return myBary;
253}
254
255Standard_Real BRepBlend_AppFuncRoot::MaximalSection() const
256{
257 Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
258 return Func->GetSectionSize();
259}
260
261void 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//
81bba717 272// Purpose : Find point solution with parameter Param (on 2 Surfaces)
7fd59977 273//
81bba717 274// Algorithm :
275// 1) Approximative solution is found from already calculated Points
276// 2) Convergence is done by a method of type Newton
7fd59977 277//
81bba717 278// Possible causes of fails :
279// - Singularity on surfaces.
280// - no information oin the "line" resulting from processing.
7fd59977 281//
282//================================================================================
283
284Standard_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();
81bba717 290 // (1) Find a point of init
7fd59977 291 Standard_Integer I1=1, I2=myLine->NbPoints(), Index;
292 Standard_Real t1, t2;
293
81bba717 294 // (1.a) It is checked if it is inside
7fd59977 295 if (Param < myLine->Point(I1).Parameter()) {return Standard_False;}
296 if (Param > myLine->Point(I2).Parameter()) {return Standard_False;}
297
81bba717 298 // (1.b) Find the interval
7fd59977 299 Trouve = SearchLocation(Param, I1, I2, Index);
300
81bba717 301 // (1.c) If the point is already calculated it is returned
7fd59977 302 if (Trouve) {
303 Pnt = myLine->Point(Index);
304 Vec(XInit,Pnt);
305 }
306 else {
81bba717 307 // (1.d) Intialisation by linear interpolation
7fd59977 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
81bba717 323 // (2) Calculation of the solution ------------------------
7fd59977 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 DEB
332 cout << "AppFunc : RNLD Not done en t = " << Param << endl;
333# endif
334 return Standard_False;
335 }
336 rsnld.Root(Sol);
337
81bba717 338 // (3) Storage of the point
7fd59977 339 Point(Func,Param,Sol,Pnt);
340
81bba717 341 // (4) Insertion of the point if the calculation seems long.
7fd59977 342 if ((!Trouve)&&(rsnld.NbIterations()>3)) {
343#ifdef DEB
81bba717 344 cout << "Evaluation in t = " << Param << "given" << endl;
7fd59977 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//
81bba717 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)
7fd59977 362//================================================================================
363Standard_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