src/PLib/PLib_HermitJacobi.cxx

   1 // File:        PLib_HermitJacobi.cxx
   2 // Created:     Wed Oct 22 11:19:36 1997
   3 // Author:      Sergey SOKOLOV
   4 //              <ssv@velox.nnov.matra-dtv.fr>
   5
   6 #include <PLib_HermitJacobi.ixx>
   7 #include <PLib.hxx>
   8 #include <PLib_LocalArray.hxx>
   9 #include <TColStd_HArray1OfReal.hxx>
  10
  11 //=======================================================================
  12 //function : PLib_HermitJacobi
  13 //purpose  :
  14 //=======================================================================
  15
  16 PLib_HermitJacobi::PLib_HermitJacobi(const Standard_Integer WorkDegree,
  17                                      const GeomAbs_Shape ConstraintOrder) :
  18                                      myH(1,2*(PLib::NivConstr(ConstraintOrder)+1),
  19                                          1,2*(PLib::NivConstr(ConstraintOrder)+1)),
  20                                      myWCoeff(1,2*(PLib::NivConstr(ConstraintOrder)+1)+1)
  21 {
  22   Standard_Integer NivConstr = PLib::NivConstr(ConstraintOrder);
  23   PLib::HermiteCoefficients(-1.,1.,NivConstr,NivConstr,myH);
  24
  25   myJacobi = new PLib_JacobiPolynomial (WorkDegree,ConstraintOrder);
  26
  27   myWCoeff.Init(0.);
  28   myWCoeff(1) = 1.;
  29   switch(NivConstr) {
  30     case 0: myWCoeff(3) = -1.; break;
  31     case 1: myWCoeff(3) = -2.; myWCoeff(5) = 1.; break;
  32     case 2: myWCoeff(3) = -3.; myWCoeff(5) = 3.; myWCoeff(7) = -1.; break;
  33   }
  34 }
  35
  36 //=======================================================================
  37 //function : MaxError
  38 //purpose  :
  39 //=======================================================================
  40
  41 Standard_Real PLib_HermitJacobi::MaxError(const Standard_Integer Dimension,
  42                                           Standard_Real& HermJacCoeff,
  43                                           const Standard_Integer NewDegree) const
  44 {
  45   return myJacobi->MaxError(Dimension,HermJacCoeff,NewDegree);
  46 }
  47
  48 //=======================================================================
  49 //function : ReduceDegree
  50 //purpose  :
  51 //=======================================================================
  52
  53 void PLib_HermitJacobi::ReduceDegree(const Standard_Integer Dimension,
  54                                      const Standard_Integer MaxDegree,
  55                                      const Standard_Real Tol,
  56                                      Standard_Real& HermJacCoeff,
  57                                      Standard_Integer& NewDegree,
  58                                      Standard_Real& MaxError) const
  59 {
  60   myJacobi->ReduceDegree(Dimension,MaxDegree,Tol,
  61                          HermJacCoeff,NewDegree,MaxError);
  62 }
  63
  64 //=======================================================================
  65 //function : AverageError
  66 //purpose  :
  67 //=======================================================================
  68
  69 Standard_Real PLib_HermitJacobi::AverageError(const Standard_Integer Dimension,
  70                                               Standard_Real& HermJacCoeff,
  71                                               const Standard_Integer NewDegree) const
  72 {
  73   return myJacobi->AverageError(Dimension,HermJacCoeff,NewDegree);
  74 }
  75
  76 //=======================================================================
  77 //function : ToCoefficients
  78 //purpose  :
  79 //=======================================================================
  80
  81 void PLib_HermitJacobi::ToCoefficients(const Standard_Integer Dimension,
  82                                        const Standard_Integer Degree,
  83                                        const TColStd_Array1OfReal& HermJacCoeff,
  84                                        TColStd_Array1OfReal& Coefficients) const
  85 {
  86   Standard_Integer i,k,idim,i1,i2;
  87   Standard_Real h1, h2;
  88   Standard_Integer NivConstr  = this->NivConstr(),
  89                    DegreeH    = 2*NivConstr+1;
  90   Standard_Integer ibegHJC = HermJacCoeff.Lower(), kdim;
  91
  92   TColStd_Array1OfReal AuxCoeff(0,(Degree+1)*Dimension-1);
  93   AuxCoeff.Init(0.);
  94
  95   for (k=0; k<=DegreeH; k++) {
  96     kdim =  k*Dimension;
  97     for (i=0; i<=NivConstr; i++) {
  98       h1 =  myH(i+1, k+1);
  99       h2 =  myH(i+NivConstr+2,k+1);
 100       i1 =  ibegHJC + i*Dimension;
 101       i2 =  ibegHJC + (i+NivConstr+1)*Dimension;
 102
 103       for (idim=0; idim<Dimension; idim++) {
 104         AuxCoeff(idim + kdim) += HermJacCoeff(i1 + idim) * h1 +
 105                                  HermJacCoeff(i2 + idim) * h2;
 106       }
 107     }
 108   }
 109   kdim = (Degree+1)*Dimension;
 110   for (k=(DegreeH+1)*Dimension; k<kdim; k++) {
 111     AuxCoeff(k) = HermJacCoeff(ibegHJC + k);
 112   }
 113
 114   if(Degree > DegreeH)
 115     myJacobi->ToCoefficients(Dimension,Degree,AuxCoeff,Coefficients);
 116   else {
 117     Standard_Integer ibegC = Coefficients.Lower();
 118     kdim = (Degree+1)*Dimension;
 119     for(k=0; k < kdim; k++)
 120         Coefficients(ibegC+k) = AuxCoeff(k);
 121   }
 122 }
 123
 124 //=======================================================================
 125 //function : D0123
 126 //purpose  : common part of D0,D1,D2,D3 (FORTRAN subroutine MPOBAS)
 127 //=======================================================================
 128
 129 void PLib_HermitJacobi::D0123(const Standard_Integer NDeriv,
 130                               const Standard_Real U,
 131                               TColStd_Array1OfReal& BasisValue,
 132                               TColStd_Array1OfReal& BasisD1,
 133                               TColStd_Array1OfReal& BasisD2,
 134                               TColStd_Array1OfReal& BasisD3)
 135 {
 136   PLib_LocalArray jac0 (4 * 20);
 137   PLib_LocalArray jac1 (4 * 20);
 138   PLib_LocalArray jac2 (4 * 20);
 139   PLib_LocalArray jac3 (4 * 20);
 140   PLib_LocalArray wvalues (4);
 141
 142   Standard_Integer i, j;
 143   Standard_Integer NivConstr  = this->NivConstr(),
 144                    WorkDegree = this->WorkDegree(),
 145                    DegreeH    = 2*NivConstr+1;
 146   Standard_Integer ibeg0 = BasisValue.Lower(),
 147                    ibeg1 = BasisD1.Lower(),
 148                    ibeg2 = BasisD2.Lower(),
 149                    ibeg3 = BasisD3.Lower();
 150   Standard_Integer JacDegree = WorkDegree-DegreeH-1;
 151   Standard_Real W0;
 152
 153   TColStd_Array1OfReal JacValue0(jac0[0], 0, Max(0,JacDegree));
 154   TColStd_Array1OfReal WValues(wvalues[0],0,NDeriv);
 155   WValues.Init(0.);
 156
 157 // Evaluation des polynomes d'hermite
 158   math_Matrix HermitValues(0,DegreeH, 0, NDeriv, 0.);
 159   if(NDeriv == 0)
 160     for (i=0; i<=DegreeH; i++) {
 161       PLib::NoDerivativeEvalPolynomial(U,DegreeH,1, DegreeH,
 162                                        myH(i+1,1), HermitValues(i,0));
 163     }
 164   else
 165     for (i=0; i<=DegreeH; i++) {
 166       PLib::EvalPolynomial(U,NDeriv,DegreeH,1,
 167                            myH(i+1,1), HermitValues(i,0));
 168     }
 169
 170 // Evaluation des polynomes de Jaccobi
 171   if(JacDegree >= 0) {
 172
 173     switch (NDeriv) {
 174     case 0 :
 175       myJacobi->D0(U,JacValue0);
 176       break;
 177     case 1 :
 178       {
 179         TColStd_Array1OfReal JacValue1(jac1[0], 0, JacDegree);
 180         myJacobi->D1(U,JacValue0,JacValue1);
 181         break;
 182       }
 183     case 2 :
 184       {
 185         TColStd_Array1OfReal JacValue1(jac1[0], 0, JacDegree);
 186         TColStd_Array1OfReal JacValue2(jac2[0], 0, JacDegree);
 187         myJacobi->D2(U,JacValue0,JacValue1,JacValue2);
 188         break;
 189       }
 190     case 3 :
 191       {
 192         TColStd_Array1OfReal JacValue1(jac1[0], 0, JacDegree);
 193         TColStd_Array1OfReal JacValue2(jac2[0], 0, JacDegree);
 194         TColStd_Array1OfReal JacValue3(jac3[0], 0, JacDegree);
 195         myJacobi->D3(U,JacValue0,JacValue1,JacValue2,JacValue3);
 196       }
 197     }
 198
 199 // Evaluation de W(t)
 200     if(NDeriv == 0)
 201       PLib::NoDerivativeEvalPolynomial(U,DegreeH+1,1,DegreeH+1,myWCoeff(1), WValues(0));
 202     else
 203       PLib::EvalPolynomial(U,NDeriv,DegreeH+1,1,myWCoeff(1), WValues(0));
 204   }
 205
 206 // Evaluation a l'ordre 0
 207   for (i=0; i<=DegreeH; i++) {
 208     BasisValue(ibeg0+i) = HermitValues(i,0);
 209   }
 210   W0 = WValues(0);
 211   for (i=DegreeH+1, j=0; i<=WorkDegree; i++, j++) {
 212     BasisValue(ibeg0+i) = W0 * jac0[j];
 213   }
 214
 215 // Evaluation a l'ordre 1
 216   if (NDeriv >= 1) {
 217     Standard_Real W1=WValues(1);
 218     for (i=0; i<=DegreeH; i++) {
 219       BasisD1(ibeg1+i) = HermitValues(i,1);
 220     }
 221     for (i=DegreeH+1, j=0; i<=WorkDegree; i++, j++) {
 222       BasisD1(ibeg1+i) = W0 * jac1[j] +
 223                          W1 * jac0[j];
 224     }
 225 // Evaluation a l'ordre 2
 226     if (NDeriv >= 2) {
 227       Standard_Real W2=WValues(2);
 228       for (i=0; i<=DegreeH; i++) {
 229         BasisD2(ibeg2+i) = HermitValues(i,2);
 230       }
 231       for (i=DegreeH+1, j=0; i<=WorkDegree; i++, j++) {
 232         BasisD2(ibeg2+i) =
 233           W0 * jac2[j] + 2 * W1 * jac1[j] + W2 * jac0[j];
 234       }
 235
 236 // Evaluation a l'ordre 3
 237       if (NDeriv == 3) {
 238          Standard_Real W3 = WValues(3);
 239         for (i=0; i<=DegreeH; i++) {
 240           BasisD3(ibeg3+i) = HermitValues(i,3);
 241         }
 242         for (i=DegreeH+1, j=0; i<=WorkDegree; i++,j++) {
 243           BasisD3(ibeg3+i) = W0*jac3[j] + W3*jac0[j]
 244                            + 3*(W1*jac2[j] + W2*jac1[j]);
 245         }
 246       }
 247     }
 248   }
 249 }
 250
 251 //=======================================================================
 252 //function : D0
 253 //purpose  :
 254 //=======================================================================
 255
 256 void PLib_HermitJacobi::D0(const Standard_Real U, TColStd_Array1OfReal& BasisValue)
 257 {
 258   D0123(0,U,BasisValue,BasisValue,BasisValue,BasisValue);
 259 }
 260
 261 //=======================================================================
 262 //function : D1
 263 //purpose  :
 264 //=======================================================================
 265
 266 void PLib_HermitJacobi::D1(const Standard_Real U,
 267                            TColStd_Array1OfReal& BasisValue, TColStd_Array1OfReal& BasisD1)
 268 {
 269   D0123(1,U,BasisValue,BasisD1,BasisD1,BasisD1);
 270 }
 271
 272 //=======================================================================
 273 //function : D2
 274 //purpose  :
 275 //=======================================================================
 276
 277 void PLib_HermitJacobi::D2(const Standard_Real U, TColStd_Array1OfReal& BasisValue,
 278                            TColStd_Array1OfReal& BasisD1,TColStd_Array1OfReal& BasisD2)
 279 {
 280   D0123(2,U,BasisValue,BasisD1,BasisD2,BasisD2);
 281 }
 282
 283 //=======================================================================
 284 //function : D3
 285 //purpose  :
 286 //=======================================================================
 287
 288 void PLib_HermitJacobi::D3(const Standard_Real U,
 289                            TColStd_Array1OfReal& BasisValue, TColStd_Array1OfReal& BasisD1,
 290                            TColStd_Array1OfReal& BasisD2,    TColStd_Array1OfReal& BasisD3)
 291 {
 292   D0123(3,U,BasisValue,BasisD1,BasisD2,BasisD3);
 293 }