0023024: Update headers of OCCT files
[occt.git] / src / AdvApprox / AdvApprox_ApproxAFunction.cxx
CommitLineData
b311480e 1// Created on: 1995-05-29
2// Created by: Xavier BENVENISTE
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
7fd59977 21
7fd59977 22// Modified PMN 22/05/1996 : Pb de recopie des poles
23// Modified PMN 20/08/1996 : Introduction de la decoupe parametrable.
24// (plus verif des derives en debug)
25// Modified PMN 20/08/1996 : Linearisation de F(t) => Les sous Espaces
26// facilement approchable ont de petites erreurs.
27// Modified PMN 15/04/1997 : Gestion fine de la continuite aux lieux de decoupes
7fd59977 28
29#include <math_Vector.hxx>
30#include <AdvApprox_ApproxAFunction.ixx>
31#include <AdvApprox_EvaluatorFunction.hxx>
32#include <AdvApprox_DichoCutting.hxx>
33#include <AdvApprox_SimpleApprox.hxx>
34
35#include <Precision.hxx>
36#include <GeomAbs_Shape.hxx>
37#include <Convert_CompPolynomialToPoles.hxx>
38#include <BSplCLib.hxx>
39
40#include <TColStd_HArray1OfReal.hxx>
41#include <TColStd_HArray2OfReal.hxx>
42#include <TColStd_HArray1OfInteger.hxx>
43#include <TColStd_Array1OfReal.hxx>
44#include <TColStd_Array1OfInteger.hxx>
45
46#include <gp_Vec.hxx>
47#include <gp_Vec2d.hxx>
48
49#include <TColgp_Array1OfPnt2d.hxx>
50#include <TColgp_Array1OfPnt.hxx>
51#include <TColgp_Array1OfPnt2d.hxx>
52#include <TColgp_Array1OfPnt.hxx>
53#include <TColgp_HArray2OfPnt2d.hxx>
54#include <TColgp_HArray2OfPnt.hxx>
55
56
57#include <Standard_ConstructionError.hxx>
58#include <Standard_OutOfRange.hxx>
59
60#include <PLib.hxx>
61#include <PLib_JacobiPolynomial.hxx>
62
63
64#ifdef DEB
65
66static Standard_Boolean AdvApprox_Debug = 0;
67
68//=====================================================
69static void MAPDBN(const Standard_Integer dimension,
70 const Standard_Real Debut,
71 const Standard_Real Fin,
72 AdvApprox_EvaluatorFunction& Evaluator,
73 const Standard_Integer Iordre)
74// Objet : Controle par difference finis, des derives
75// Warning : En mode Debug, uniquement
76///===================================================
77{
78 Standard_Integer derive, OrdreDer, ier, NDIMEN = dimension;
79 Standard_Real* Ptr;
80 Standard_Real h = 1.e-4*(Fin-Debut+1.e-3), eps = 1.e-3, t, ab[2];
81 math_Vector V1(1,NDIMEN), V2(1,NDIMEN), Diff(1,NDIMEN), Der(1,NDIMEN);
82
83 if (h < 100*RealEpsilon()) {h = 100*RealEpsilon();}
84 ab[0] = Debut;
85 ab[1] = Fin;
86
87 for (OrdreDer=1, derive = 0;
88 OrdreDer <= Iordre; OrdreDer++, derive++) {
89 // Verif au debut
90 Ptr = (Standard_Real*) &V1.Value(1);
91 t = Debut+h;
92 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
93
94 Ptr = (Standard_Real*) &V2.Value(1);
95 t = Debut;
96 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
97
98 Diff = (V1 - V2)/h;
99
100 Ptr = (Standard_Real*) &Der.Value(1);
101 t = Debut;
102 Evaluator(&NDIMEN, ab, &t, &OrdreDer, Ptr, &ier);
103
104
105 if ( (Diff-Der).Norm() > eps * (Der.Norm()+1) ) {
106 cout << " Debug Ft au parametre t+ = " << t << endl;
107 cout << " Positionement sur la derive "<< OrdreDer
108 << " : " << Der << endl;
109 cout << " Erreur estime : " << (Der-Diff) << endl;
110 }
111
112 // Verif a la fin
113 Ptr = (Standard_Real*) &V1.Value(1);
114 t = Fin-h;
115 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
116
117 Ptr = (Standard_Real*) &V2.Value(1);
118 t = Fin;
119 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
120
121 Diff = (V2 - V1)/h;
122
123 Ptr = (Standard_Real*) &Der.Value(1);
124 t = Fin;
125 Evaluator(&NDIMEN, ab, &t, &OrdreDer, Ptr, &ier);
126
127
128 if ( (Diff-Der).Norm() > eps * (Der.Norm()+1) ) {
129 cout << " Debug Ft au parametre t- = " << t << endl;
130 cout << " Positionement sur la derive "<< OrdreDer
131 << " : " << Der << endl;
132 cout << " Erreur estime : " << (Der-Diff) << endl;
133 }
134 }
135}
136#endif
137
138//===================================================================
139static void PrepareConvert(const Standard_Integer NumCurves,
140 const Standard_Integer MaxDegree,
141 const Standard_Integer ContinuityOrder,
142 const Standard_Integer Num1DSS,
143 const Standard_Integer Num2DSS,
144 const Standard_Integer Num3DSS,
145 const TColStd_Array1OfInteger& NumCoeffPerCurve,
146 TColStd_Array1OfReal& Coefficients,
147 const TColStd_Array2OfReal& PolynomialIntervals,
148 const TColStd_Array1OfReal& TrueIntervals,
149 const TColStd_Array1OfReal& LocalTolerance,
150 TColStd_Array1OfReal& ErrorMax,
151 TColStd_Array1OfInteger& Continuity)
152// Pour determiner les continuites locales
153//====================================================================
154
155{
156 // Declaration
157 Standard_Boolean isCi;
158 Standard_Integer icurve, idim, iordre, ii,
159 Dimension=Num1DSS + 2*Num2DSS + 3*Num3DSS,
160 NbSpace = Num1DSS + Num2DSS + Num3DSS;
161 Standard_Real diff, moy, facteur1, facteur2, normal1, normal2, eps;
162 Standard_Real *Res1, *Res2, *Val1, *Val2;
163 Standard_Real *Coef1, *Coef2;
164 Standard_Integer RealDegree = Max(MaxDegree + 1, 2 * ContinuityOrder + 2);
165
166 gp_Vec V1,V2;
167 gp_Vec2d v1, v2;
168
169 TColStd_Array1OfReal Result(1, 2*(ContinuityOrder+1)*Dimension);
170 TColStd_Array1OfReal Prec(1, NbSpace), Suivant(1, NbSpace);
171
172
173 Res1 = (Standard_Real*) &(Result.ChangeValue(1));
174 Res2 = (Standard_Real*) &(Result.ChangeValue((ContinuityOrder+1)*Dimension+1));
175
176
177 //Init
178 Continuity.Init(0);
179 if (ContinuityOrder ==0) return;
180
181 for (icurve=1; icurve<NumCurves; icurve++) {
182 // Init et positionement au noeud
183 isCi = Standard_True;
184 Coef1 = (Standard_Real*) &(Coefficients.Value(
185 (icurve-1) * Dimension * RealDegree + Coefficients.Lower()));
186 Coef2 = Coef1 + Dimension * RealDegree;
187 Standard_Integer Deg1 = NumCoeffPerCurve(NumCoeffPerCurve.Lower()+icurve-1) - 1;
188 PLib::EvalPolynomial(PolynomialIntervals(icurve, 2),
189 ContinuityOrder,
190 Deg1,
191 Dimension,
192 Coef1[0],
193 Res1[0]);
194 Standard_Integer Deg2 = NumCoeffPerCurve(NumCoeffPerCurve.Lower()+icurve) - 1;
195 PLib::EvalPolynomial(PolynomialIntervals(icurve+1, 1),
196 ContinuityOrder,
197 Deg2,
198 Dimension,
199 Coef2[0],
200 Res2[0]);
201
202 // Verif dans chaque sous espaces.
203 for (iordre=1; iordre<=ContinuityOrder && isCi; iordre++) {
204
205 // fixing a bug PRO18577
206
207 Standard_Real Toler = 1.0e-5;
208
209 Standard_Real f1_dividend = PolynomialIntervals(icurve,2)-PolynomialIntervals(icurve,1);
210 Standard_Real f2_dividend = PolynomialIntervals(icurve+1,2)-PolynomialIntervals(icurve+1,1);
211 Standard_Real f1_divizor = TrueIntervals(icurve+1)-TrueIntervals(icurve);
212 Standard_Real f2_divizor = TrueIntervals(icurve+2)-TrueIntervals(icurve+1);
213 Standard_Real fract1, fract2;
214
215 if( Abs(f1_divizor) < Toler ) // this is to avoid divizion by zero
216 //in this case fract1 = 5.14755758946803e-85
217 facteur1 = 0.0;
218 else{ fract1 = f1_dividend/f1_divizor;
219 facteur1 = Pow( fract1, iordre);
220 }
221 if( Abs(f2_divizor) < Toler )
222 //in this case fract2 = 6.77193633669143e-313
223 facteur2 = 0.0;
224 else{ fract2 = f2_dividend/f2_divizor;
225 facteur2 = Pow( fract2, iordre);
226 }
227 normal1 = Pow(f1_divizor, iordre);
228 normal2 = Pow(f2_divizor, iordre);
229
230
231 idim = 1;
232 Val1 = Res1+iordre*Dimension-1;
233 Val2 = Res2+iordre*Dimension-1;
234
235 for (ii=1; ii<=Num1DSS && isCi; ii++, idim++) {
236 eps = LocalTolerance(idim)*0.01;
237 diff = Abs (Val1[ii]*facteur1 - Val2[ii]*facteur2);
238 moy = Abs (Val1[ii]*facteur1 + Val2[ii]*facteur2);
239 // Un premier controle sur la valeur relative
240 if (diff > moy*1.e-9) {
241 Prec(idim) = diff*normal1;
242 Suivant(idim) = diff*normal2;
243 // Et un second sur le majorant d'erreur engendree
244 if (Prec(idim)>eps ||
245 Suivant(idim)>eps) isCi=Standard_False;
246 }
247 else {
248 Prec(idim) = 0;
249 Suivant(idim) = 0;
250 }
251 }
252
253 Val1+=Num1DSS;
254 Val2+=Num1DSS;
255 for (ii=1; ii<=Num2DSS && isCi; ii++, idim++, Val1+=2,Val2+=2) {
256 eps = LocalTolerance(idim)*0.01;
257 v1.SetCoord(Val1[1],Val1[2]);
258 v2.SetCoord(Val2[1],Val2[2]);
259 v1 *= facteur1;
260 v2 *= facteur2;
261 diff = Abs(v1.X() - v2.X()) + Abs(v1.Y() - v2.Y());
262 moy = Abs(v1.X() + v2.X()) + Abs(v1.Y() + v2.Y());
263 if (diff > moy*1.e-9) {
264 Prec(idim) = diff*normal1;
265 Suivant(idim) = diff*normal2;
266 if (Prec(idim)>eps ||
267 Suivant(idim)>eps) isCi=Standard_False;
268 }
269 else {
270 Prec(idim) = 0;
271 Suivant(idim) = 0;
272 }
273 }
274
275 for (ii=1; ii<=Num3DSS && isCi; ii++, idim++, Val1+=3,Val2+=3) {
276 eps = LocalTolerance(idim)*0.01;
277 V1.SetCoord(Val1[1], Val1[2], Val1[3]);
278 V2.SetCoord(Val2[1], Val2[2], Val2[3]);
279 V1 *= facteur1;
280 V2 *= facteur2;
281 diff = Abs(V1.X() - V2.X()) + Abs(V1.Y() - V2.Y()) +
282 Abs(V1.Z() - V2.Z());
283 moy = Abs(V1.X() + V2.X()) + Abs(V1.Y() + V2.Y()) +
284 Abs(V1.Z() + V2.Z());
285 if (diff > moy*1.e-9) {
286 Prec(idim) = diff*normal1;
287 Suivant(idim) = diff*normal2;
288 if (Prec(idim)>eps ||
289 Suivant(idim)>eps) isCi=Standard_False;
290 }
291 else {
292 Prec(idim) = 0;
293 Suivant(idim) = 0;
294 }
295 }
296 // Si c'est bon on update le tout
297 if (isCi) {
298 Continuity(icurve+1) = iordre;
299 Standard_Integer index = (icurve-1)*NbSpace+1;
300 for (ii=index, idim=1; idim<=NbSpace; ii++,idim++) {
301 ErrorMax(ii) += Prec(idim);
302 ErrorMax(ii+NbSpace) += Suivant(idim);
303 }
304 }
305 }
306 }
307}
308
309
310//=======================================================================
311//function : ApproxFunction
312//
313//purpose : Approximation d' UNE fonction non polynomiale (dans l'espace de
314// dimension NDIMEN) par N courbes polynomiales, par la methode
315// des moindres carres. Le parametre de la fonction est conserve.
316//
317// ARGUMENTS D'ENTREE :
318//C ------------------
319//C NDIMEN : Dimension totale de l' espace (somme des dimensions
320//C des sous-espaces).
321//C NBSESP : Nombre de sous-espaces "independants".
322//C NDIMSE : Table des dimensions des sous-espaces.
323//C ABFONC : Bornes de l' intervalle (a,b) de definition de la
324//C fonction a approcher.
325//C FONCNP : Fonction externe de positionnement sur la fonction non
326//C polynomiale a approcher.
327//C IORDRE : Ordre de contrainte aux extremites de chaque courbe
328//C polynomiale cree :
329//C -1 = pas de contraintes,
330//C 0 = contraintes de passage aux bornes (i.e. C0),
331//C 1 = C0 + contraintes de derivees 1eres (i.e. C1),
332//C 2 = C1 + contraintes de derivees 2ndes (i.e. C2).
333//C NBCRMX : Nombre maxi de courbes polynomiales a calculer.
334//C NCFLIM : Nombre LIMITE de coeff des "courbes" polynomiales
335//C d' approximation. Doit etre superieur ou egal a
336//C 2*IORDRE + 2 et inferieur ou egal a 50 et a NCOFMX).
337//C Limite a 14 apres l'appel a VRIRAZ pour eviter le bruit.
338//C EPSAPR : Table des erreurs d' approximations souhaitees
339//C sous-espace par sous-espace.
340//C ICODEO : Code d' init. des parametres de discretisation.
341//C (choix de NBPNT et de NDJAC).
342//C = 1 calcul rapide avec precision moyenne sur les coeff.
343//C = 2 calcul rapide avec meilleure precision " " ".
344//C = 3 calcul un peu plus lent avec bonne precision ".
345//C = 4 calcul lent avec la meilleure precision possible ".
346//C
347//C
348//C ARGUMENTS DE SORTIE :
349//C -------------------
350//C NBCRBE : Nombre de courbes polynomiales creees.
351//C NCFTAB : Table des nombres de coeff. significatifs des NBCRBE
352//C "courbes" calculees.
353//C CRBTAB : Tableau des coeff. des "courbes" polynomiales calculees.
354//C Doit etre dimensionne a CRBTAB(NDIMEN,NCOFMX,NBCRMX).
355//C TABINT : Table des NBCRBE + 1 bornes des intervalles de decoupe de
356//C FONCNP.
357//C ERRMAX : Table des erreurs (sous-espace par sous espace)
358//C MAXIMALES commises dans l' approximation de FONCNP par
359//C par les NBCRBE courbes polynomiales.
360//C ERRMOY : Table des erreurs (sous-espace par sous espace)
361//C MOYENNES commises dans l' approximation de FONCNP par
362//C par les NBCRBE courbes polynomiales.
363//C IERCOD : Code d' erreur :
364//C -1 = ERRMAX > EPSAPR pour au moins un des sous-espace.
365//C On a alors NCRBMX courbes calculees, certaines de
366//C degre mathematique min(NCFLIM,14)-1 ou la precision
367//C demandee n' est pas atteinte.
368//C -2 = Pb dans le mode DEBUG
369//C 0 = Tout est ok.
370//C 1 = Pb d' incoherence des entrees.
371//C 10 = Pb de calcul de l' interpolation des contraintes.
372//C 13 = Pb dans l' allocation dynamique.
373//C 33 = Pb dans la recuperation des donnees du block data
374//C des coeff. d' integration par la methode de GAUSS.
375//C >100 Pb dans l' evaluation de FONCNP, le code d' erreur
376//C renvoye est egal au code d' erreur de FONCNP + 100.
377//C
378//
379// --> La i-eme courbe polynomiale creee correspond a l' approximation
380//C de FONCNP sur l' intervalle (TABINT(i),TABINT(i+1)). TABINT(i)
381// est une suite STRICTEMENT monotone avec TABINT(1)=ABFONC(1) et
382// TABINT(NBCRBE+1)=ABFONC(2).
383//
384// --> Si IERCOD=-1, la precision demandee n' est pas atteinte (ERRMAX
385//C est superieur a EPSAPR sur au moins l' un des sous-espaces), mais
386// on donne le meilleur resultat possible pour min(NCFLIM,14),NBCRMX
387// et EPSAPR choisis par l' utilisateur Dans ce cas (ainsi que pour
388//C IERCOD=0), on a une solution.
389//C
390//C--> ATTENTION : Toute modification du calcul de NDJAC et NBPNT
391//C doit etre reportee dans le programme MAPNBP.
392//
393// HISTORIQUE DES MODIFICATIONS :
394// --------------------------------
395//
396// 20-08-1996 : PMN ; Creation d'apres la routine Fortran MAPFNC
397//======================================================================
398void AdvApprox_ApproxAFunction::Approximation(
399 const Standard_Integer TotalDimension,
400 // Dimension totale de l' espace
401 // (somme des dimensions des sous-espaces).
402 const Standard_Integer TotalNumSS,//Nombre de sous-espaces "independants".
403 const TColStd_Array1OfInteger& LocalDimension,//dimensions des sous-espaces.
404 const Standard_Real First,
405 const Standard_Real Last,
406 // Intervalle (a,b) de definition de la fonction a approcher.
407 AdvApprox_EvaluatorFunction& Evaluator,
408 // Fonction externe de positionnement sur la fonction a approcher.
409 const AdvApprox_Cutting& CutTool,
410 // Maniere de decouper en 2 un intervalle.
411 const Standard_Integer ContinuityOrder,
412 // ordre de continutie a respecter (-1, 0, 1, 2)
413 const Standard_Integer NumMaxCoeffs,
414 // Nombre LIMITE de coeff des "courbes" polynomiales
415 // d' approximation. Doit etre superieur ou egal a
416 // Doit etre superieur ou egal a 2*ContinuityOrder + 2
417 // Limite a 14 pour eviter le bruit.
418 const Standard_Integer MaxSegments,
419 //Nombre maxi de courbes polynomiales a calculer.
420 const TColStd_Array1OfReal& LocalTolerancesArray,
421 //Tolerances d'approximation souhaitees sous-espace par sous-espace.
422 const Standard_Integer code_precis,
423 //Code d' init. des parametres de discretisation.
424 //C (choix de NBPNT et de NDJAC).
425 //C = 1 calcul rapide avec precision moyenne sur les coeff.
426 //C = 2 calcul rapide avec meilleure precision " " ".
427 //C = 3 calcul un peu plus lent avec bonne precision ".
428 //C = 4 calcul lent avec la meilleure precision possible ".
429 Standard_Integer& NumCurves,
430 TColStd_Array1OfInteger& NumCoeffPerCurveArray,
431 TColStd_Array1OfReal& CoefficientArray,
432 TColStd_Array1OfReal& IntervalsArray,
433 TColStd_Array1OfReal& ErrorMaxArray,
434 TColStd_Array1OfReal& AverageErrorArray,
435 Standard_Integer& ErrorCode)
436{
437// Standard_Real EpsPar = Precision::Confusion();
438 Standard_Integer IDIM, NUPIL,TheDeg;
439#ifdef DEB
440 Standard_Integer NDIMEN = TotalDimension;
441#endif
442 Standard_Boolean isCut = Standard_False;
443
444// Defintion des Tableaux C
445 Standard_Real * TABINT = (Standard_Real*) &IntervalsArray(1);
446
447
448 ErrorCode=0;
449 CoefficientArray.Init(0);
450
451//-------------------------- Verification des entrees ------------------
452
453 if ((MaxSegments<1)|| (Abs(Last-First)<1.e-9))
454 {ErrorCode=1;
455 return;}
456
457//--> La dimension totale doit etre la somme des dimensions des
458// sous-espaces
459 IDIM=0;
460 for ( Standard_Integer I=1; I<=TotalNumSS; I++) {IDIM += LocalDimension(I);}
461 if (IDIM != TotalDimension)
462 {ErrorCode=1;
463 return;}
464 GeomAbs_Shape Continuity=GeomAbs_C0;
465 switch (ContinuityOrder) {
466 case 0:
467 Continuity = GeomAbs_C0 ;
468 break ;
469 case 1:
470 Continuity = GeomAbs_C1 ;
471 break ;
472 case 2:
473 Continuity = GeomAbs_C2 ;
474 break ;
475 default:
476 Standard_ConstructionError::Raise();
477 }
478
479//--------------------- Choix du nombre de points ----------------------
480//---------- et du degre de lissage dans la base orthogonale -----------
481//--> NDJAC est le degre de "travail" dans la base orthogonale.
482
483
484 Standard_Integer NbGaussPoints, WorkDegree;
485
486 PLib::JacobiParameters(Continuity, NumMaxCoeffs-1, code_precis,
487 NbGaussPoints, WorkDegree);
488// NDJAC=WorkDegree;
489
490//------------------ Initialisation de la gestion des decoupes ---------
491
492 TABINT[0]=First;
493 TABINT[1]=Last;
494 NUPIL=1;
495 NumCurves=0;
496
497//C**********************************************************************
498//C APPROXIMATION AVEC DECOUPES
499//C**********************************************************************
500 Handle(PLib_JacobiPolynomial) JacobiBase = new (PLib_JacobiPolynomial) (WorkDegree, Continuity);
501//Portage HP le compilateur refuse le debranchement
502 Standard_Integer IS ;
503 Standard_Boolean goto_fin_de_boucle;
504 Standard_Integer MaxDegree = NumMaxCoeffs-1;
505 AdvApprox_SimpleApprox Approx (TotalDimension, TotalNumSS,
506 Continuity,
507 WorkDegree, NbGaussPoints,
508 JacobiBase, Evaluator);
509 while ((NUPIL-NumCurves)!=0) {
510//--> Lorsque l' on a atteint le haut de la pile, c' est fini !
511
512//Portage HP le compilateur refuse le debranchement
513 goto_fin_de_boucle=Standard_False;
514
515//---- Calcul de la courbe d' approximation dans la base de Jacobi -----
516
517 Approx.Perform ( LocalDimension, LocalTolerancesArray,
518 TABINT[NumCurves], TABINT[NumCurves+1],
519 MaxDegree);
520 if (!Approx.IsDone()) {
521 ErrorCode = 1;
522 return;
523 }
524
525
526
527//---------- Calcul du degre de la courbe et de l' erreur max ----------
528
529 IDIM=0;
530 NumCoeffPerCurveArray(NumCurves + 1)=0;
531
532// L'erreur doit etre satisfaite sur tous les sous-espaces sinon, on decoupe
533
534 Standard_Boolean MaxMaxErr=Standard_True;
535 for ( IS=1; IS<=TotalNumSS; IS++)
536 { if (Approx.MaxError(IS) > LocalTolerancesArray(IS))
537 { MaxMaxErr = Standard_False;
538 break;
539 }
540 }
541
542 if (MaxMaxErr == Standard_True)
543 {
544 NumCurves++;
545// NumCoeffPerCurveArray(NumCurves) = Approx.Degree()+1;
546 }
547 else
548 {
549//-> ...sinon on essai de decouper l' intervalle courant en 2...
550 Standard_Real TMIL;
551 Standard_Boolean Large;
552
553 Large = CutTool.Value(TABINT[NumCurves], TABINT[NumCurves+1],
554 TMIL);
555
556 if (NUPIL<MaxSegments && Large) {
557
558// if (IS!=1) {NumCurves--;}
559 isCut = Standard_True; // Ca y est, on le sait !
560 Standard_Real * IDEB = TABINT+NumCurves+1;
561 Standard_Real * IDEB1 = IDEB+1;
562 Standard_Integer ILONG= NUPIL-NumCurves-1;
563 for (Standard_Integer iI=ILONG; iI>=0; iI--) {
564 IDEB1[iI] = IDEB[iI];
565 }
566 IDEB[0] = TMIL;
567 NUPIL++;
568//--> ... et on recommence.
569//Portage HP le compilateur refuse le debranchement
570 goto_fin_de_boucle=Standard_True;
571// break;
572 }
573 else {
574//--> Si la pile est pleine...
575// pour l'instant ErrorCode=-1;
576 NumCurves++;
577// NumCoeffPerCurveArray(NumCurves) = Approx.Degree()+1;
578 }
579 }
580// IDIM += NDSES;
581//Portage HP le compilateur refuse le debranchement
582 if (goto_fin_de_boucle) continue;
583 for (IS=1; IS<=TotalNumSS; IS++) {
584 ErrorMaxArray.SetValue(IS+(NumCurves-1)*TotalNumSS,Approx. MaxError(IS));
585//---> ...et calcul de l' erreur moyenne.
586 AverageErrorArray.SetValue(IS+(NumCurves-1)*TotalNumSS,Approx.AverageError(IS));
587 }
588
589 Handle(TColStd_HArray1OfReal) HJacCoeff = Approx.Coefficients();
590 TheDeg = Approx.Degree();
591 if (isCut && (TheDeg < 2*ContinuityOrder+1) )
592 // Pour ne pas bruiter les derives aux bouts, et garder une continuite
593 // correcte sur la BSpline resultat.
594 TheDeg = 2*ContinuityOrder+1;
595
596 NumCoeffPerCurveArray(NumCurves) = TheDeg+1;
597 TColStd_Array1OfReal Coefficients(0,(TheDeg+1)*TotalDimension-1);
598 JacobiBase->ToCoefficients (TotalDimension, TheDeg,
599 HJacCoeff->Array1(), Coefficients);
600 Standard_Integer i,j, f = (TheDeg+1)*TotalDimension;
601 for (i=0,j=(NumCurves-1)*TotalDimension*NumMaxCoeffs+1;
602 i<f; i++, j++) {
603 CoefficientArray.SetValue(j, Coefficients.Value(i));
604 }
605
606#ifdef DEB
607 // Test des derives par difference finis
608 Standard_Integer IORDRE = ContinuityOrder;
609
610 if (IORDRE>0 && AdvApprox_Debug) {
611 MAPDBN(NDIMEN, TABINT[NumCurves-1],
612 TABINT[NumCurves], Evaluator, IORDRE);
613 }
614#endif
615//Portage HP le compilateur refuse le debranchement
616// fin_de_boucle:
617// {;} fin de la boucle while
618 }
619 return;
620}
621//=======================================================================
622//function : AdvApprox_ApproxAFunction
623//purpose : Constructeur avec Decoupe Dichotomique
624//=======================================================================
625
626AdvApprox_ApproxAFunction::
627AdvApprox_ApproxAFunction(const Standard_Integer Num1DSS,
628 const Standard_Integer Num2DSS,
629 const Standard_Integer Num3DSS,
630 const Handle_TColStd_HArray1OfReal& OneDTol,
631 const Handle_TColStd_HArray1OfReal& TwoDTol,
632 const Handle_TColStd_HArray1OfReal& ThreeDTol,
633 const Standard_Real First,
634 const Standard_Real Last,
635 const GeomAbs_Shape Continuity,
636 const Standard_Integer MaxDeg,
637 const Standard_Integer MaxSeg,
638 const AdvApprox_EvaluatorFunction& Func) :
639 my1DTolerances(OneDTol),
640 my2DTolerances(TwoDTol),
641 my3DTolerances(ThreeDTol),
642 myFirst(First),
643 myLast(Last),
644 myContinuity(Continuity),
645 myMaxDegree(MaxDeg),
646 myMaxSegments(MaxSeg),
647 myDone(Standard_False),
648 myHasResult(Standard_False),
649 myEvaluator((Standard_Address)&Func)
650{
651 AdvApprox_DichoCutting Cut;
652 Perform(Num1DSS, Num2DSS, Num3DSS, Cut);
653}
654
655AdvApprox_ApproxAFunction::
656AdvApprox_ApproxAFunction(const Standard_Integer Num1DSS,
657 const Standard_Integer Num2DSS,
658 const Standard_Integer Num3DSS,
659 const Handle_TColStd_HArray1OfReal& OneDTol,
660 const Handle_TColStd_HArray1OfReal& TwoDTol,
661 const Handle_TColStd_HArray1OfReal& ThreeDTol,
662 const Standard_Real First,
663 const Standard_Real Last,
664 const GeomAbs_Shape Continuity,
665 const Standard_Integer MaxDeg,
666 const Standard_Integer MaxSeg,
667 const AdvApprox_EvaluatorFunction& Func,
668 const AdvApprox_Cutting& CutTool) :
669 my1DTolerances(OneDTol),
670 my2DTolerances(TwoDTol),
671 my3DTolerances(ThreeDTol),
672 myFirst(First),
673 myLast(Last),
674 myContinuity(Continuity),
675 myMaxDegree(MaxDeg),
676 myMaxSegments(MaxSeg),
677 myDone(Standard_False),
678 myHasResult(Standard_False),
679 myEvaluator((Standard_Address)&Func)
680{
681 Perform(Num1DSS, Num2DSS, Num3DSS, CutTool);
682}
683
684//=======================================================================
685//function : AdvApprox_ApproxAFunction::Perform
686//purpose : Make all the Work !!
687//=======================================================================
688void AdvApprox_ApproxAFunction::Perform(const Standard_Integer Num1DSS,
689 const Standard_Integer Num2DSS,
690 const Standard_Integer Num3DSS,
691 const AdvApprox_Cutting& CutTool)
692{
693 if (Num1DSS < 0 ||
694 Num2DSS < 0 ||
695 Num3DSS < 0 ||
696 Num1DSS + Num2DSS + Num3DSS <= 0 ||
697 myLast < myFirst ||
698 myMaxDegree < 1 ||
699 myMaxSegments < 0)
700 Standard_ConstructionError::Raise();
701 if (myMaxDegree > 14) {
702 myMaxDegree = 14 ;
703 }
704 //
705 // Input
706 //
707 myNumSubSpaces[0] = Num1DSS ;
708 myNumSubSpaces[1] = Num2DSS ;
709 myNumSubSpaces[2] = Num3DSS ;
710 Standard_Integer TotalNumSS =
711 Num1DSS + Num2DSS + Num3DSS,
712 ii,
713 jj,
714 kk,
715 index,
716 dim_index,
717 local_index;
718 Standard_Integer TotalDimension =
719 myNumSubSpaces[0] + 2 * myNumSubSpaces[1] + 3 * myNumSubSpaces[2] ;
720 Standard_Real error_value ;
721
722 Standard_Integer ContinuityOrder=0 ;
723 switch (myContinuity) {
724 case GeomAbs_C0:
725 ContinuityOrder = 0 ;
726 break ;
727 case GeomAbs_C1:
728 ContinuityOrder = 1 ;
729 break ;
730 case GeomAbs_C2:
731 ContinuityOrder = 2 ;
732 break ;
733 default:
734 Standard_ConstructionError::Raise();
735 }
736 Standard_Real ApproxStartEnd[2] ;
737 Standard_Integer NumMaxCoeffs = Max(myMaxDegree + 1, 2 * ContinuityOrder + 2);
738 myMaxDegree = NumMaxCoeffs - 1;
739 Standard_Integer code_precis = 1 ;
740 //
741 // WARNING : the polynomial coefficients are the
742 // taylor expansion of the polynomial at 0.0e0 !
743 //
744 ApproxStartEnd[0] = -1.0e0 ;
745 ApproxStartEnd[1] = 1.0e0 ;
746
747 TColStd_Array1OfInteger LocalDimension(1,TotalNumSS) ;
748
749
750 index = 1 ;
751 TColStd_Array1OfReal LocalTolerances(1,TotalNumSS) ;
752
753 for(jj = 1; jj <= myNumSubSpaces[0] ; jj++) {
754 LocalTolerances.SetValue(index,my1DTolerances->Value(jj)) ;
755 LocalDimension.SetValue(index,1) ;
756 index += 1 ;
757 }
758 for(jj = 1 ; jj <= myNumSubSpaces[1] ; jj++) {
759 LocalTolerances.SetValue(index,my2DTolerances->Value(jj)) ;
760 LocalDimension.SetValue(index,2) ;
761 index += 1 ;
762 }
763 for(jj = 1; jj <= myNumSubSpaces[2] ; jj++) {
764 LocalTolerances.SetValue(index,my3DTolerances->Value(jj)) ;
765 LocalDimension.SetValue(index,3) ;
766 index += 1 ;
767 }
768 //
769 // Ouput
770 //
771
772 Standard_Integer ErrorCode = 0,
773 NumCurves,
774 size =
775 myMaxSegments * NumMaxCoeffs * TotalDimension ;
776 Handle(TColStd_HArray1OfInteger) NumCoeffPerCurvePtr =
777 new TColStd_HArray1OfInteger (1,myMaxSegments) ;
778
779 Handle(TColStd_HArray1OfReal) LocalCoefficientsPtr =
780 new TColStd_HArray1OfReal(1,size) ;
781
782 Handle(TColStd_HArray1OfReal) IntervalsPtr =
783 new TColStd_HArray1OfReal (1,myMaxSegments+1) ;
784
785 TColStd_Array1OfReal ErrorMax(1,myMaxSegments * TotalNumSS) ;
786
787 TColStd_Array1OfReal AverageError(1,myMaxSegments * TotalNumSS) ;
788
789 Approximation (TotalDimension, TotalNumSS, LocalDimension,
790 myFirst, myLast,
791 *(AdvApprox_EvaluatorFunction*)myEvaluator,
792 CutTool,
793 ContinuityOrder,
794 NumMaxCoeffs, myMaxSegments,
795 LocalTolerances, code_precis,
796 NumCurves, // Nombre de courbe en sortie
797 NumCoeffPerCurvePtr->ChangeArray1(), // Nbre de coeff par courbe
798 LocalCoefficientsPtr->ChangeArray1(),// Les Coeffs solutions
799 IntervalsPtr->ChangeArray1(), // La Table des decoupes
800 ErrorMax, // Majoration de l'erreur
801 AverageError, // Erreur moyenne constatee
802 ErrorCode) ;
803 if (ErrorCode == 0 || ErrorCode == -1) {
804 //
805 // si tout est OK ou bien on a un resultat dont l une des erreurs max est
806 // plus grande que la tolerance demandee
807
808 TColStd_Array1OfInteger TabContinuity(1, NumCurves);
809 TColStd_Array2OfReal PolynomialIntervalsPtr (1,NumCurves,1,2);
810 for (ii = PolynomialIntervalsPtr.LowerRow() ;
811 ii <= PolynomialIntervalsPtr.UpperRow() ;
812 ii++) {
813 // On force un degre 1 minimum (PRO5474)
814 NumCoeffPerCurvePtr->ChangeValue(ii) =
815 Max(2, NumCoeffPerCurvePtr->Value(ii));
816 PolynomialIntervalsPtr.SetValue(ii,1,ApproxStartEnd[0]) ;
817 PolynomialIntervalsPtr.SetValue(ii,2,ApproxStartEnd[1]) ;
818 }
819
820 PrepareConvert(NumCurves, myMaxDegree, ContinuityOrder,
821 Num1DSS, Num2DSS, Num3DSS,
822 NumCoeffPerCurvePtr->Array1(),
823 LocalCoefficientsPtr->ChangeArray1(),
824 PolynomialIntervalsPtr, IntervalsPtr->Array1(),
825 LocalTolerances, ErrorMax,
826 TabContinuity);
827
828 Convert_CompPolynomialToPoles
829 AConverter(NumCurves,
830 TotalDimension,
831 myMaxDegree,
832 TabContinuity,
833 NumCoeffPerCurvePtr->Array1(),
834 LocalCoefficientsPtr->Array1(),
835 PolynomialIntervalsPtr,
836 IntervalsPtr->Array1());
837
838 if (AConverter.IsDone()) {
839 Handle(TColStd_HArray2OfReal) PolesPtr ;
840 AConverter.Poles(PolesPtr) ;
841 AConverter.Knots(myKnots) ;
842 AConverter.Multiplicities(myMults) ;
843 myDegree = AConverter.Degree() ;
844 index = 0 ;
845 if (myNumSubSpaces[0] > 0) {
846 my1DPoles = new TColStd_HArray2OfReal(1,PolesPtr->ColLength(),
847 1,myNumSubSpaces[0]) ;
848 my1DMaxError = new TColStd_HArray1OfReal(1,myNumSubSpaces[0]) ;
849 my1DAverageError = new TColStd_HArray1OfReal(1,myNumSubSpaces[0]) ;
850 for (ii = 1 ; ii <= PolesPtr->ColLength() ; ii++) {
851 for (jj = 1 ; jj <= myNumSubSpaces[0] ; jj++) {
852 my1DPoles->SetValue(ii,jj, PolesPtr->Value(ii,jj)) ;
853 }
854 }
855
856 for (jj = 1 ; jj <= myNumSubSpaces[0] ; jj++) {
857 error_value = 0.0e0 ;
858 for (ii = 1 ; ii <= NumCurves ; ii++) {
859 local_index = (ii - 1) * TotalNumSS ;
860 error_value = Max(ErrorMax(local_index + jj),error_value) ;
861
862 }
863 my1DMaxError->SetValue(jj, error_value) ;
864 }
865 for (jj = 1 ; jj <= myNumSubSpaces[0] ; jj++) {
866 error_value = 0.0e0 ;
867 for (ii = 1 ; ii <= NumCurves ; ii++) {
868 local_index = (ii - 1) * TotalNumSS ;
869 error_value += AverageError(local_index + jj) ;
870
871 }
872 error_value /= (Standard_Real) NumCurves ;
873 my1DAverageError->SetValue(jj, error_value) ;
874 }
875 index += myNumSubSpaces[0] ;
876 }
877
878 dim_index = index; //Pour le cas ou il n'y a pas de 2D
879
880 if (myNumSubSpaces[1] > 0) {
881 gp_Pnt2d Point2d ;
882 my2DPoles = new TColgp_HArray2OfPnt2d(1,PolesPtr->ColLength(),
883 1,myNumSubSpaces[1]) ;
884 my2DMaxError = new TColStd_HArray1OfReal(1,myNumSubSpaces[1]) ;
885 my2DAverageError = new TColStd_HArray1OfReal(1,myNumSubSpaces[1]) ;
886 for (ii = 1 ; ii <= PolesPtr->ColLength() ; ii++) {
887 for (jj = 1 ; jj <= myNumSubSpaces[1] ; jj++) {
888 local_index = index + (jj-1) * 2 ;
889 for (kk = 1; kk <= 2 ; kk++) {
890 Point2d.SetCoord(kk,
891 PolesPtr->Value(ii,local_index + kk)) ;
892 }
893 my2DPoles->SetValue(ii,jj, Point2d) ;
894 }
895 }
896
897 for (jj = 1 ; jj <= myNumSubSpaces[1] ; jj++) {
898 error_value = 0.0e0 ;
899 for (ii = 1 ; ii <= NumCurves ; ii++) {
900 local_index = (ii - 1) * TotalNumSS + index ;
901 error_value = Max(ErrorMax(local_index + jj),error_value) ;
902
903 }
904 my2DMaxError->SetValue(jj, error_value) ;
905 }
906 for (jj = 1 ; jj <= myNumSubSpaces[1] ; jj++) {
907 error_value = 0.0e0 ;
908 for (ii = 1 ; ii <= NumCurves ; ii++) {
909 local_index = (ii - 1) * TotalNumSS + index ;
910 error_value += AverageError(local_index + jj) ;
911
912 }
913 error_value /= (Standard_Real) NumCurves ;
914 my2DAverageError->SetValue(jj, error_value) ;
915 }
916 index += myNumSubSpaces[1] ;
917 // Pour les poles il faut doubler le decalage :
918 dim_index = index + myNumSubSpaces[1];
919 }
920
921 if (myNumSubSpaces[2] > 0) {
922 gp_Pnt Point ;
923 my3DPoles = new TColgp_HArray2OfPnt(1,PolesPtr->ColLength(),
924 1,myNumSubSpaces[2]) ;
925 my3DMaxError = new TColStd_HArray1OfReal(1,myNumSubSpaces[2]) ;
926 my3DAverageError = new TColStd_HArray1OfReal(1,myNumSubSpaces[2]) ;
927 for (ii = 1 ; ii <= PolesPtr->ColLength() ; ii++) {
928 for (jj = 1 ; jj <= myNumSubSpaces[2] ; jj++) {
929 local_index = dim_index + (jj-1) * 3 ;
930 for (kk = 1; kk <= 3 ; kk++) {
931 Point.SetCoord(kk,
932 PolesPtr->Value(ii,local_index + kk)) ;
933 }
934 my3DPoles->SetValue(ii,jj,Point) ;
935
936 }
937 }
938
939 for (jj = 1 ; jj <= myNumSubSpaces[2] ; jj++) {
940 error_value = 0.0e0 ;
941 for (ii = 1 ; ii <= NumCurves ; ii++) {
942 local_index = (ii - 1) * TotalNumSS + index ;
943 error_value = Max(ErrorMax(local_index + jj),error_value) ;
944
945 }
946 my3DMaxError->SetValue(jj, error_value) ;
947 }
948 for (jj = 1 ; jj <= myNumSubSpaces[2] ; jj++) {
949 error_value = 0.0e0 ;
950 for (ii = 1 ; ii <= NumCurves ; ii++) {
951 local_index = (ii - 1) * TotalNumSS + index ;
952 error_value += AverageError(local_index + jj) ;
953
954 }
955 error_value /= (Standard_Real) NumCurves ;
956 my3DAverageError->SetValue(jj, error_value) ;
957 }
958 }
959 if (ErrorCode == 0) {
960 myDone = Standard_True ;
961 myHasResult = Standard_True ;
962 }
963 else if (ErrorCode == -1) {
964 myHasResult = Standard_True ;
965 }
966
967 }
968 }
969}
970
971
972//=======================================================================
973//function : Poles
974//purpose :
975//=======================================================================
976
977void AdvApprox_ApproxAFunction::Poles(const Standard_Integer Index,
978 TColgp_Array1OfPnt& P) const
979{
980 Standard_Integer ii ;
981 for (ii = P.Lower() ; ii <= P.Upper() ; ii++) {
982 P.SetValue(ii,my3DPoles->Value(ii,Index)) ;
983 }
984}
985
986
987//=======================================================================
988//function : NbPoles
989//purpose :
990//=======================================================================
991
992Standard_Integer AdvApprox_ApproxAFunction::NbPoles() const
993{ if (myDone || myHasResult) return BSplCLib::NbPoles(myDegree,
994 Standard_False,
995 myMults->Array1()) ;
996 return 0 ; }
997
998//=======================================================================
999//function : Poles
1000//purpose :
1001//=======================================================================
1002
1003void AdvApprox_ApproxAFunction::Poles2d(const Standard_Integer Index,
1004 TColgp_Array1OfPnt2d& P) const
1005{
1006 Standard_Integer ii ;
1007 for (ii = P.Lower() ; ii <= P.Upper() ; ii++) {
1008 P.SetValue(ii,my2DPoles->Value(ii,Index)) ;
1009 }
1010}
1011//=======================================================================
1012//function : Poles
1013//purpose :
1014//=======================================================================
1015
1016void AdvApprox_ApproxAFunction::Poles1d(const Standard_Integer Index,
1017 TColStd_Array1OfReal& P) const
1018{
1019 Standard_Integer ii ;
1020 for (ii = P.Lower() ; ii <= P.Upper() ; ii++) {
1021 P.SetValue(ii,my1DPoles->Value(ii,Index)) ;
1022 }
1023}
1024//=======================================================================
1025//function : MaxError
1026//purpose :
1027//=======================================================================
1028Handle(TColStd_HArray1OfReal)
1029 AdvApprox_ApproxAFunction::MaxError(const Standard_Integer D) const
1030
1031{
1032 Handle(TColStd_HArray1OfReal) EPtr ;
1033 if (D <= 0 ||
1034 D > 3) {
1035
1036 Standard_OutOfRange::Raise() ;
1037 }
1038 switch (D) {
1039 case 1:
1040 EPtr = my1DMaxError ;
1041 break ;
1042 case 2:
1043 EPtr = my2DMaxError ;
1044 break ;
1045 case 3:
1046 EPtr = my3DMaxError ;
1047 break ;
1048 }
1049 return EPtr ;
1050}
1051//=======================================================================
1052//function : MaxError
1053//purpose :
1054//=======================================================================
1055Standard_Real AdvApprox_ApproxAFunction::MaxError(const Standard_Integer D,
1056 const Standard_Integer Index) const
1057{
1058 Handle(TColStd_HArray1OfReal) EPtr =
1059 MaxError(D) ;
1060
1061 return EPtr->Value(Index) ;
1062}
1063
1064//=======================================================================
1065//function : AverageError
1066//purpose :
1067//=======================================================================
1068Handle(TColStd_HArray1OfReal)
1069 AdvApprox_ApproxAFunction::AverageError(const Standard_Integer D) const
1070
1071{
1072 Handle(TColStd_HArray1OfReal) EPtr ;
1073 if (D <= 0 ||
1074 D > 3) {
1075
1076 Standard_OutOfRange::Raise() ;
1077 }
1078 switch (D) {
1079 case 1:
1080 EPtr = my1DAverageError ;
1081 break ;
1082 case 2:
1083 EPtr = my2DAverageError ;
1084 break ;
1085 case 3:
1086 EPtr = my3DAverageError ;
1087 break ;
1088 }
1089 return EPtr ;
1090}
1091//=======================================================================
1092//function : AverageError
1093//purpose :
1094//=======================================================================
1095
1096Standard_Real AdvApprox_ApproxAFunction::AverageError(
1097 const Standard_Integer D,
1098 const Standard_Integer Index) const
1099
1100{
1101 Handle(TColStd_HArray1OfReal) EPtr =
1102 AverageError(D) ;
1103 return EPtr->Value(Index) ;
1104}
1105
1106void AdvApprox_ApproxAFunction::Dump(Standard_OStream& o) const
1107{
1108 Standard_Integer ii;
1109 o << "Dump of ApproxAFunction" << endl;
1110 if (myNumSubSpaces[0] > 0) {
1111 o << "Error(s) 1d = " << endl;
1112 for (ii=1; ii <= myNumSubSpaces[0]; ii++) {
1113 o << " " << MaxError(1, ii) << endl;
1114 }
1115 }
1116
1117 if (myNumSubSpaces[1] > 0) {
1118 o << "Error(s) 2d = " << endl;
1119 for (ii=1; ii <= myNumSubSpaces[1]; ii++) {
1120 o << " " << MaxError(2, ii) << endl;
1121 }
1122 }
1123
1124 if (myNumSubSpaces[2] > 0) {
1125 o << "Error(s) 3d = " << endl;
1126 for (ii=1; ii <= myNumSubSpaces[2]; ii++) {
1127 o << " " << MaxError(3, ii) << endl;
1128 }
1129 }
1130}