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