src/AdvApprox/AdvApprox_ApproxAFunction.cxx

   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
  21
  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
  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
  66 static Standard_Boolean AdvApprox_Debug = 0;
  67
  68 //=====================================================
  69 static 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 //===================================================================
 139 static 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 //======================================================================
 398 void 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
 626 AdvApprox_ApproxAFunction::
 627 AdvApprox_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
 655 AdvApprox_ApproxAFunction::
 656 AdvApprox_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 //=======================================================================
 688 void 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,
 790                   TotalNumSS,   LocalDimension,
 791                   myFirst,      myLast,
 792                   *(AdvApprox_EvaluatorFunction*)myEvaluator,
 793                   CutTool, ContinuityOrder, NumMaxCoeffs,
 794                   myMaxSegments, LocalTolerances, code_precis,
 795                   NumCurves,                            // Nombre de courbe en sortie
 796                   NumCoeffPerCurvePtr->ChangeArray1(),  // Nbre de coeff par courbe
 797                   LocalCoefficientsPtr->ChangeArray1(), // Les Coeffs solutions
 798                   IntervalsPtr->ChangeArray1(),         // La Table des decoupes
 799                   ErrorMax,                             // Majoration de l'erreur
 800                   AverageError,                         // Erreur moyenne constatee
 801                   ErrorCode) ;
 802
 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
 977 void 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
 992 Standard_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
1003 void 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
1016 void 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 //=======================================================================
1028 Handle(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 //=======================================================================
1055 Standard_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 //=======================================================================
1068 Handle(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
1096 Standard_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
1106 void  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 }