src/PLib/PLib_DoubleJacobiPolynomial.cxx

   1 // Created on: 1997-05-28
   2 // Created by: Sergey SOKOLOV
   3 // Copyright (c) 1997-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
  18 #include <math_Vector.hxx>
  19 #include <PLib_DoubleJacobiPolynomial.hxx>
  20 #include <PLib_JacobiPolynomial.hxx>
  21
  22 //=======================================================================
  23 //function : PLib_DoubleJacobiPolynomial
  24 //purpose  :
  25 //=======================================================================
  26 PLib_DoubleJacobiPolynomial::PLib_DoubleJacobiPolynomial()
  27
  28 {
  29 }
  30
  31 //=======================================================================
  32 //function : PLib_DoubleJacobiPolynomial
  33 //purpose  :
  34 //=======================================================================
  35
  36 PLib_DoubleJacobiPolynomial::PLib_DoubleJacobiPolynomial(const Handle(PLib_JacobiPolynomial)& JacPolU,
  37                                                          const Handle(PLib_JacobiPolynomial)& JacPolV) :
  38                                                          myJacPolU(JacPolU),
  39                                                          myJacPolV(JacPolV)
  40 {
  41   Handle (TColStd_HArray1OfReal) TabMaxU =
  42     new TColStd_HArray1OfReal (0,JacPolU->WorkDegree()-2*(JacPolU->NivConstr()+1));
  43   JacPolU->MaxValue(TabMaxU->ChangeArray1());
  44   myTabMaxU = TabMaxU;
  45
  46   Handle (TColStd_HArray1OfReal) TabMaxV =
  47     new TColStd_HArray1OfReal (0,JacPolV->WorkDegree()-2*(JacPolV->NivConstr()+1));
  48   JacPolV->MaxValue(TabMaxV->ChangeArray1());
  49   myTabMaxV = TabMaxV;
  50 }
  51
  52 //=======================================================================
  53 //function : MaxErrorU
  54 //purpose  :
  55 //=======================================================================
  56
  57 Standard_Real
  58 PLib_DoubleJacobiPolynomial::MaxErrorU(const Standard_Integer Dimension,
  59                                        const Standard_Integer DegreeU,
  60                                        const Standard_Integer DegreeV,
  61                                        const Standard_Integer dJacCoeff,
  62                                        const TColStd_Array1OfReal& JacCoeff) const
  63 {
  64   Standard_Integer ii,idim,dJac,MinU,MinV,WorkDegreeU,WorkDegreeV;
  65   Standard_Real Bid0;
  66
  67   math_Vector MaxErrDim(1,Dimension,0.);
  68
  69   MinU = 2*(myJacPolU->NivConstr()+1);
  70   MinV = 2*(myJacPolV->NivConstr()+1);
  71   WorkDegreeU = myJacPolU->WorkDegree();
  72   WorkDegreeV = myJacPolV->WorkDegree();
  73
  74   Bid0 = myTabMaxV->Value(DegreeV-MinV);
  75   for (idim=1; idim<=Dimension; idim++) {
  76     dJac = dJacCoeff + (idim-1)*(WorkDegreeU+1)*(WorkDegreeV+1);
  77     for (ii=MinU; ii<=DegreeU; ii++) {
  78       MaxErrDim(idim) += (Abs(JacCoeff(ii + DegreeV*(WorkDegreeU+1) + dJac)) *
  79                           myTabMaxU->Value(ii-MinU) * Bid0);
  80     }
  81   }
  82   return (MaxErrDim.Norm());
  83 }
  84
  85 //=======================================================================
  86 //function : MaxErrorV
  87 //purpose  :
  88 //=======================================================================
  89
  90 Standard_Real
  91 PLib_DoubleJacobiPolynomial::MaxErrorV(const Standard_Integer Dimension,
  92                                        const Standard_Integer DegreeU,
  93                                        const Standard_Integer DegreeV,
  94                                        const Standard_Integer dJacCoeff,
  95                                        const TColStd_Array1OfReal& JacCoeff) const
  96 {
  97   Standard_Integer jj,idim,dJac,MinU,MinV,WorkDegreeU,WorkDegreeV;
  98   Standard_Real Bid0;
  99
 100   math_Vector MaxErrDim(1,Dimension,0.);
 101
 102   MinU = 2*(myJacPolU->NivConstr()+1);
 103   MinV = 2*(myJacPolV->NivConstr()+1);
 104   WorkDegreeU = myJacPolU->WorkDegree();
 105   WorkDegreeV = myJacPolV->WorkDegree();
 106
 107   Bid0 = myTabMaxU->Value(DegreeU-MinU);
 108   for (idim=1; idim<=Dimension; idim++) {
 109     dJac = dJacCoeff + (idim-1)*(WorkDegreeU+1)*(WorkDegreeV+1);
 110     for (jj=MinV; jj<=DegreeV; jj++) {
 111       MaxErrDim(idim) += (Abs(JacCoeff(DegreeU + jj*(WorkDegreeU+1) + dJac)) *
 112                           myTabMaxV->Value(jj-MinV) * Bid0);
 113     }
 114   }
 115   return (MaxErrDim.Norm());
 116 }
 117
 118 //=======================================================================
 119 //function : MaxError
 120 //purpose  :
 121 //=======================================================================
 122
 123 Standard_Real
 124 PLib_DoubleJacobiPolynomial::MaxError(const Standard_Integer Dimension,
 125                                       const Standard_Integer MinDegreeU,
 126                                       const Standard_Integer MaxDegreeU,
 127                                       const Standard_Integer MinDegreeV,
 128                                       const Standard_Integer MaxDegreeV,
 129                                       const Standard_Integer dJacCoeff,
 130                                       const TColStd_Array1OfReal& JacCoeff,
 131                                       const Standard_Real Error) const
 132 {
 133   Standard_Integer ii,jj,idim,dJac,MinU,MinV,WorkDegreeU,WorkDegreeV;
 134   Standard_Real Bid0,Bid1;
 135
 136   math_Vector MaxErrDim(1,Dimension,0.);
 137
 138   MinU = 2*(myJacPolU->NivConstr()+1);
 139   MinV = 2*(myJacPolV->NivConstr()+1);
 140   WorkDegreeU = myJacPolU->WorkDegree();
 141   WorkDegreeV = myJacPolV->WorkDegree();
 142
 143 //------------------- Calcul du majorant de l'erreur max ---------------
 144 //----- lorsque sont enleves les coeff. d'indices MinDegreeU a MaxDegreeU ------
 145 //---------------- en U et d'indices MinDegreeV a MaxDegreeV en V --------------
 146
 147   for (idim=1; idim<=Dimension; idim++) {
 148     dJac = dJacCoeff + (idim-1)*(WorkDegreeU+1)*(WorkDegreeV+1);
 149     Bid1 = 0.;
 150     for (jj=MinDegreeV; jj<=MaxDegreeV; jj++) {
 151       Bid0 = 0.;
 152       for (ii=MinDegreeU; ii<=MaxDegreeU; ii++) {
 153         Bid0 += fabs(JacCoeff(ii + jj*(WorkDegreeU+1) + dJac)) * myTabMaxU->Value(ii-MinU);
 154       }
 155       Bid1 += Bid0 * myTabMaxV->Value(jj-MinV);
 156     }
 157     MaxErrDim(idim) = Bid1;
 158   }
 159
 160 //----------------------- Calcul de l' erreur max ----------------------
 161
 162   math_Vector MaxErr2(1,2);
 163   MaxErr2(1) = Error;
 164   MaxErr2(2) = MaxErrDim.Norm();
 165   return (MaxErr2.Norm());
 166 }
 167
 168 //=======================================================================
 169 //function : ReduceDegree
 170 //purpose  :
 171 //=======================================================================
 172
 173 void PLib_DoubleJacobiPolynomial::ReduceDegree(const Standard_Integer Dimension,
 174                                                const Standard_Integer MinDegreeU,
 175                                                const Standard_Integer MaxDegreeU,
 176                                                const Standard_Integer MinDegreeV,
 177                                                const Standard_Integer MaxDegreeV,
 178                                                const Standard_Integer dJacCoeff,
 179                                                const TColStd_Array1OfReal& JacCoeff,
 180                                                const Standard_Real EpmsCut,
 181                                                Standard_Real& MaxError,
 182                                                Standard_Integer& NewDegreeU,
 183                                                Standard_Integer& NewDegreeV) const
 184 {
 185   Standard_Integer NewU,NewV;
 186   Standard_Real ErrU,ErrV;
 187
 188   NewU = MaxDegreeU;
 189   NewV = MaxDegreeV;
 190   math_Vector MaxErr2(1,2);
 191
 192 //**********************************************************************
 193 //-------------------- Coupure des coefficients ------------------------
 194 //**********************************************************************
 195
 196   do {
 197
 198 //------------------- Calcul du majorant de l'erreur max ---------------
 199 //----- lorsque sont enleves les coeff. d'indices MinU a NewU ------
 200 //---------------- en U, le degre en V etant fixe a NewV -----------------
 201     if (NewV > MinDegreeV)
 202       ErrV = MaxErrorU(Dimension,NewU,NewV,dJacCoeff,JacCoeff);
 203     else {
 204       ErrV = 2*EpmsCut;
 205     }
 206
 207 //------------------- Calcul du majorant de l'erreur max ---------------
 208 //----- lorsque sont enleves les coeff. d'indices MinV a NewV ------
 209 //---------------- en V, le degre en U etant fixe a NewU -----------------
 210     if (NewU > MinDegreeU)
 211       ErrU = MaxErrorV(Dimension,NewU,NewV,dJacCoeff,JacCoeff);
 212     else {
 213       ErrU = 2*EpmsCut;
 214     }
 215
 216 //----------------------- Calcul de l' erreur max ----------------------
 217     MaxErr2(1) = MaxError;
 218     MaxErr2(2) = ErrU;
 219     ErrU = MaxErr2.Norm();
 220     MaxErr2(2) = ErrV;
 221     ErrV = MaxErr2.Norm();
 222
 223     if (ErrU > ErrV) {
 224       if (ErrV < EpmsCut) {
 225         MaxError = ErrV;
 226         NewV--;
 227       }
 228     }
 229     else {
 230       if (ErrU < EpmsCut) {
 231         MaxError = ErrU;
 232         NewU--;
 233       }
 234     }
 235   }
 236   while ((ErrU > ErrV && ErrV <= EpmsCut) || (ErrV >= ErrU && ErrU <= EpmsCut));
 237
 238 //-------------------------- Recuperation des degres -------------------
 239
 240   NewDegreeU = Max(NewU,1);
 241   NewDegreeV = Max(NewV,1);
 242 }
 243
 244 //=======================================================================
 245 //function : AverageError
 246 //purpose  :
 247 //=======================================================================
 248
 249 Standard_Real
 250 PLib_DoubleJacobiPolynomial::AverageError(const Standard_Integer Dimension,
 251                                               const Standard_Integer DegreeU,
 252                                               const Standard_Integer DegreeV,
 253                                               const Standard_Integer dJacCoeff,
 254                                               const TColStd_Array1OfReal& JacCoeff) const
 255 {
 256   Standard_Integer ii,jj,idim,dJac,IDebU,IDebV,MinU,MinV,WorkDegreeU,WorkDegreeV;
 257   Standard_Real Bid0,Bid1,AverageErr;
 258
 259 //----------------------------- Initialisations ------------------------
 260
 261   IDebU = 2*(myJacPolU->NivConstr()+1);
 262   IDebV = 2*(myJacPolV->NivConstr()+1);
 263   MinU = Max(IDebU,DegreeU);
 264   MinV = Max(IDebV,DegreeV);
 265   WorkDegreeU = myJacPolU->WorkDegree();
 266   WorkDegreeV = myJacPolV->WorkDegree();
 267   Bid0 = 0.;
 268
 269 //------------------ Calcul du majorant de l'erreur moyenne ------------
 270 //----- lorsque sont enleves les coeff. d'indices DegreeU a WorkDegreeU ------
 271 //---------------- en U et d'indices DegreeV a WorkDegreeV en V --------------
 272
 273   for (idim=1; idim<=Dimension; idim++) {
 274     dJac = dJacCoeff + (idim-1)*(WorkDegreeU+1)*(WorkDegreeV+1);
 275     for (jj=MinV; jj<=WorkDegreeV; jj++) {
 276       for (ii=IDebU; ii<=WorkDegreeU; ii++) {
 277         Bid1 = JacCoeff(ii + jj*(WorkDegreeU+1) + dJac);
 278         Bid0 += Bid1*Bid1;
 279       }
 280     }
 281     for (jj=IDebV; jj<=MinV-1; jj++) {
 282       for (ii=MinU; ii<=WorkDegreeU; ii++) {
 283         Bid1 = JacCoeff(ii + jj*(WorkDegreeU+1) + dJac);
 284         Bid0 += Bid1*Bid1;
 285       }
 286     }
 287   }
 288   AverageErr = sqrt(Bid0/4);
 289   return (AverageErr);
 290 }
 291
 292 //=======================================================================
 293 //function : WDoubleJacobiToCoefficients
 294 //purpose  :
 295 //=======================================================================
 296
 297 void PLib_DoubleJacobiPolynomial::WDoubleJacobiToCoefficients(const Standard_Integer Dimension,
 298                                                               const Standard_Integer DegreeU,
 299                                                               const Standard_Integer DegreeV,
 300                                                               const TColStd_Array1OfReal& JacCoeff,
 301                                                               TColStd_Array1OfReal& Coefficients) const
 302 {
 303   Standard_Integer iu,iv,idim,WorkDegreeU,WorkDegreeV;
 304
 305   Coefficients.Init(0.);
 306
 307   WorkDegreeU = myJacPolU->WorkDegree();
 308   WorkDegreeV = myJacPolV->WorkDegree();
 309
 310   TColStd_Array1OfReal Aux1(0, (DegreeU+1)*(DegreeV+1)*Dimension-1);
 311   TColStd_Array1OfReal Aux2(0, (DegreeU+1)*(DegreeV+1)*Dimension-1);
 312
 313   for (iu=0; iu<=DegreeU; iu++) {
 314     for (iv=0; iv<=DegreeV; iv++) {
 315       for (idim=1; idim<=Dimension; idim++) {
 316         Aux1(idim-1 + iv*Dimension + iu*Dimension*(DegreeV+1)) =
 317           JacCoeff(iu + iv*(WorkDegreeU+1) + (idim-1)*(WorkDegreeU+1)*(WorkDegreeV+1));
 318       }
 319     }
 320   }
 321 //   Passage dans canonique en u.
 322   myJacPolU->ToCoefficients(Dimension*(DegreeV+1),DegreeU,Aux1,Aux2);
 323
 324 //   Permutation des u et des v.
 325   for (iu=0; iu<=DegreeU; iu++) {
 326     for (iv=0; iv<=DegreeV; iv++) {
 327       for (idim=1; idim<=Dimension; idim++) {
 328         Aux1(idim-1 + iu*Dimension + iv*Dimension*(DegreeU+1)) =
 329           Aux2(idim-1 + iv*Dimension + iu*Dimension*(DegreeV+1));
 330       }
 331     }
 332   }
 333 //   Passage dans canonique en v.
 334   myJacPolV->ToCoefficients(Dimension*(DegreeU+1),DegreeV,Aux1,Aux2);
 335
 336 //   Permutation des u et des v.
 337   for (iu=0; iu<=DegreeU; iu++) {
 338     for (iv=0; iv<=DegreeV; iv++) {
 339       for (idim=1; idim<=Dimension; idim++) {
 340         Coefficients(iu + iv*(DegreeU+1) + (idim-1)*(DegreeU+1)*(DegreeV+1)) =
 341           Aux2(idim-1 + iu*Dimension + iv*Dimension*(DegreeU+1));
 342       }
 343     }
 344   }
 345 }
 346