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