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