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