src/GeomFill/GeomFill_PolynomialConvertor.cxx

   1 // Created on: 1997-07-18
   2 // Created by: Philippe MANGIN
   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 <Convert_CompPolynomialToPoles.hxx>
  19 #include <GeomFill_PolynomialConvertor.hxx>
  20 #include <gp_Mat.hxx>
  21 #include <gp_Pnt.hxx>
  22 #include <gp_Vec.hxx>
  23 #include <PLib.hxx>
  24 #include <StdFail_NotDone.hxx>
  25 #include <TColStd_Array1OfReal.hxx>
  26 #include <TColStd_HArray1OfInteger.hxx>
  27 #include <TColStd_HArray1OfReal.hxx>
  28 #include <TColStd_HArray2OfReal.hxx>
  29
  30 GeomFill_PolynomialConvertor::GeomFill_PolynomialConvertor():
  31                     Ordre(8),
  32                     myinit(Standard_False),
  33                     BH(1, Ordre, 1, Ordre)
  34 {
  35 }
  36
  37 Standard_Boolean GeomFill_PolynomialConvertor::Initialized() const
  38 {
  39  return myinit;
  40 }
  41
  42 void GeomFill_PolynomialConvertor::Init()
  43 {
  44   if (myinit) return; //On n'initialise qu'une fois
  45   Standard_Integer ii, jj;
  46   Standard_Real terme;
  47   math_Matrix H(1,Ordre, 1,Ordre), B(1,Ordre, 1,Ordre);
  48   Handle(TColStd_HArray1OfReal)
  49     Coeffs = new (TColStd_HArray1OfReal) (1, Ordre*Ordre),
  50     TrueInter  = new  (TColStd_HArray1OfReal) (1,2);
  51
  52   Handle(TColStd_HArray2OfReal)
  53     Poles1d = new (TColStd_HArray2OfReal) (1, Ordre, 1, Ordre),
  54     Inter   = new (TColStd_HArray2OfReal) (1,1,1,2);
  55
  56   //Calcul de B
  57   Inter->SetValue(1, 1, -1);
  58   Inter->SetValue(1, 2,  1);
  59   TrueInter->SetValue(1, -1);
  60   TrueInter->SetValue(2,  1);
  61
  62   Coeffs->Init(0);
  63   for (ii=1; ii<=Ordre; ii++) {  Coeffs->SetValue(ii+(ii-1)*Ordre, 1); }
  64
  65    //Convertion ancienne formules
  66      Handle(TColStd_HArray1OfInteger) Ncf
  67        = new (TColStd_HArray1OfInteger)(1,1);
  68      Ncf->Init(Ordre);
  69
  70      Convert_CompPolynomialToPoles
  71           AConverter(1, 1, 8, 8,
  72                      Ncf,
  73                      Coeffs,
  74                      Inter,
  75                      TrueInter);
  76 /*  Convert_CompPolynomialToPoles
  77          AConverter(8, Ordre-1,  Ordre-1,
  78                     Coeffs,
  79                     Inter,
  80                     TrueInter); En attente du bon Geomlite*/
  81   AConverter.Poles(Poles1d);
  82
  83   for (jj=1; jj<=Ordre; jj++) {
  84       for (ii=1; ii<=Ordre; ii++) {
  85         terme = Poles1d->Value(ii,jj);
  86         if (Abs(terme-1) < 1.e-9) terme = 1 ; //petite retouche
  87         if (Abs(terme+1) < 1.e-9) terme = -1;
  88         B(ii, jj) =  terme;
  89       }
  90   }
  91   //Calcul de H
  92   myinit = PLib::HermiteCoefficients(-1, 1, Ordre/2-1, Ordre/2-1, H);
  93   H.Transpose();
  94
  95   if (!myinit) return;
  96
  97   // reste l'essentiel
  98   BH = B * H;
  99 }
 100
 101 void GeomFill_PolynomialConvertor::Section(const gp_Pnt& FirstPnt,
 102                                        const gp_Pnt& Center,
 103                                        const gp_Vec& Dir,
 104                                        const Standard_Real Angle,
 105                                        TColgp_Array1OfPnt& Poles) const
 106 {
 107   math_Vector Vx(1, Ordre), Vy(1, Ordre);
 108   math_Vector Px(1, Ordre), Py(1, Ordre);
 109   Standard_Integer ii;
 110   Standard_Real Cos_b = Cos(Angle), Sin_b = Sin(Angle);
 111   Standard_Real beta, beta2, beta3;
 112   gp_Vec V1(Center, FirstPnt), V2;
 113   V2 =  Dir^V1;
 114   beta = Angle/2;
 115   beta2 = beta * beta;
 116   beta3 = beta * beta2;
 117
 118   // Calcul de la transformation
 119   gp_Mat M(V1.X(), V2.X(), 0,
 120            V1.Y(), V2.Y(), 0,
 121            V1.Z(), V2.Z(), 0);
 122
 123   // Calcul des contraintes  -----------
 124   Vx(1) = 1;                   Vy(1) = 0;
 125   Vx(2) = 0;                   Vy(2) = beta;
 126   Vx(3) = -beta2;              Vy(3) = 0;
 127   Vx(4) = 0;                   Vy(4) = -beta3;
 128   Vx(5) = Cos_b;               Vy(5) = Sin_b;
 129   Vx(6) = -beta*Sin_b;         Vy(6) = beta*Cos_b;
 130   Vx(7) = -beta2*Cos_b;        Vy(7) = -beta2*Sin_b;
 131   Vx(8) = beta3*Sin_b;         Vy(8) = -beta3*Cos_b;
 132
 133   // Calcul des poles
 134   Px = BH * Vx;
 135   Py = BH * Vy;
 136   gp_XYZ pnt;
 137   for (ii=1; ii<=Ordre; ii++) {
 138      pnt.SetCoord(Px(ii), Py(ii), 0);
 139      pnt *= M;
 140      pnt += Center.XYZ();
 141      Poles(ii).ChangeCoord() = pnt;
 142   }
 143 }
 144
 145 void GeomFill_PolynomialConvertor::Section(const gp_Pnt& FirstPnt,
 146                                        const gp_Vec& DFirstPnt,
 147                                        const gp_Pnt& Center,
 148                                        const gp_Vec& DCenter,
 149                                        const gp_Vec& Dir,
 150                                        const gp_Vec& DDir,
 151                                        const Standard_Real Angle,
 152                                        const Standard_Real DAngle,
 153                                        TColgp_Array1OfPnt& Poles,
 154                                        TColgp_Array1OfVec& DPoles) const
 155 {
 156   math_Vector Vx(1, Ordre), Vy(1, Ordre),
 157               DVx(1, Ordre), DVy(1, Ordre);
 158   math_Vector Px(1, Ordre), Py(1, Ordre),
 159               DPx(1, Ordre), DPy(1, Ordre);
 160   Standard_Integer ii;
 161   Standard_Real Cos_b = Cos(Angle), Sin_b = Sin(Angle);
 162   Standard_Real beta, beta2, beta3, bprim;
 163   gp_Vec V1(Center, FirstPnt), V1Prim, V2;
 164   V2 =  Dir^V1;
 165   beta = Angle/2;
 166   bprim = DAngle/2;
 167   beta2 = beta * beta;
 168   beta3 = beta * beta2;
 169
 170   // Calcul des  transformations
 171   gp_Mat M (V1.X(), V2.X(), 0,
 172             V1.Y(), V2.Y(), 0,
 173             V1.Z(), V2.Z(), 0);
 174   V1Prim = DFirstPnt - DCenter;
 175   V2 = (DDir^V1) + (Dir^V1Prim);
 176   gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
 177                 V1Prim.Y(), V2.Y(), 0,
 178                 V1Prim.Z(), V2.Z(), 0);
 179
 180
 181
 182   // Calcul des contraintes  -----------
 183   Vx(1) = 1;                   Vy(1) = 0;
 184   Vx(2) = 0;                   Vy(2) = beta;
 185   Vx(3) = -beta2;              Vy(3) = 0;
 186   Vx(4) = 0;                   Vy(4) = -beta3;
 187   Vx(5) = Cos_b;               Vy(5) = Sin_b;
 188   Vx(6) = -beta*Sin_b;         Vy(6) = beta*Cos_b;
 189   Vx(7) = -beta2*Cos_b;        Vy(7) = -beta2*Sin_b;
 190   Vx(8) = beta3*Sin_b;         Vy(8) = -beta3*Cos_b;
 191
 192   Standard_Real b_bprim = bprim*beta,
 193                 b2_bprim = bprim*beta2;
 194   DVx(1) = 0;                  DVy(1) = 0;
 195   DVx(2) = 0;                  DVy(2) = bprim;
 196   DVx(3) = -2*b_bprim;         DVy(3) = 0;
 197   DVx(4) = 0;                  DVy(4) = -3*b2_bprim;
 198   DVx(5) = -2*bprim*Sin_b;     DVy(5) = 2*bprim*Cos_b;
 199   DVx(6) = -bprim*Sin_b - 2*b_bprim*Cos_b;
 200                   DVy(6) = bprim*Cos_b - 2*b_bprim*Sin_b;
 201   DVx(7) = 2*b_bprim*(-Cos_b + beta*Sin_b);
 202                   DVy(7) = -2*b_bprim*(Sin_b+beta*Cos_b);
 203   DVx(8) = b2_bprim*(3*Sin_b + 2*beta*Cos_b);
 204                   DVy(8) = b2_bprim*(2*beta*Sin_b - 3*Cos_b);
 205
 206   // Calcul des poles
 207   Px = BH * Vx;
 208   Py = BH * Vy;
 209   DPx = BH * DVx;
 210   DPy = BH * DVy;
 211   gp_XYZ P, DP, aux;
 212
 213   for (ii=1; ii<=Ordre; ii++) {
 214      P.SetCoord(Px(ii), Py(ii), 0);
 215      Poles(ii).ChangeCoord() = M*P + Center.XYZ();
 216      P *= MPrim;
 217      DP.SetCoord(DPx(ii), DPy(ii), 0);
 218      DP *= M;
 219      aux.SetLinearForm(1, P, 1, DP, DCenter.XYZ());
 220      DPoles(ii).SetXYZ(aux);
 221   }
 222 }
 223
 224 void GeomFill_PolynomialConvertor::Section(const gp_Pnt& FirstPnt,
 225                                        const gp_Vec& DFirstPnt,
 226                                        const gp_Vec& D2FirstPnt,
 227                                        const gp_Pnt& Center,
 228                                        const gp_Vec& DCenter,
 229                                        const gp_Vec& D2Center,
 230                                        const gp_Vec& Dir,
 231                                        const gp_Vec& DDir,
 232                                        const gp_Vec& D2Dir,
 233                                        const Standard_Real Angle,
 234                                        const Standard_Real DAngle,
 235                                        const Standard_Real D2Angle,
 236                                        TColgp_Array1OfPnt& Poles,
 237                                        TColgp_Array1OfVec& DPoles,
 238                                        TColgp_Array1OfVec& D2Poles) const
 239 {
 240   math_Vector Vx(1, Ordre), Vy(1, Ordre),
 241               DVx(1, Ordre), DVy(1, Ordre),
 242               D2Vx(1, Ordre), D2Vy(1, Ordre);
 243   math_Vector Px(1, Ordre), Py(1, Ordre),
 244               DPx(1, Ordre), DPy(1, Ordre),
 245               D2Px(1, Ordre), D2Py(1, Ordre);
 246
 247   Standard_Integer ii;
 248   Standard_Real aux, Cos_b = Cos(Angle), Sin_b = Sin(Angle);
 249   Standard_Real beta, beta2, beta3, bprim, bprim2, bsecn;
 250   gp_Vec V1(Center, FirstPnt), V1Prim, V1Secn, V2;
 251   V2 =  Dir^V1;
 252   beta = Angle/2;
 253   bprim = DAngle/2;
 254   bsecn = D2Angle/2;
 255   bsecn = D2Angle/2;
 256   beta2 = beta * beta;
 257   beta3 = beta * beta2;
 258   bprim2 =  bprim*bprim;
 259
 260   // Calcul des  transformations
 261   gp_Mat M (V1.X(), V2.X(), 0,
 262             V1.Y(), V2.Y(), 0,
 263             V1.Z(), V2.Z(), 0);
 264   V1Prim = DFirstPnt - DCenter;
 265   V2 = (DDir^V1) + (Dir^V1Prim);
 266
 267   gp_Mat MPrim (V1Prim.X(), V2.X(), 0,
 268                 V1Prim.Y(), V2.Y(), 0,
 269                 V1Prim.Z(), V2.Z(), 0);
 270   V1Secn = D2FirstPnt - D2Center;
 271   V2 = DDir^V1Prim;
 272   V2 *= 2;
 273   V2 +=  (D2Dir^V1) + (Dir^V1Secn);
 274
 275
 276   gp_Mat MSecn (V1Secn.X(), V2.X(), 0,
 277                 V1Secn.Y(), V2.Y(), 0,
 278                 V1Secn.Z(), V2.Z(), 0);
 279
 280
 281   // Calcul des contraintes  -----------
 282   Vx(1) = 1;                   Vy(1) = 0;
 283   Vx(2) = 0;                   Vy(2) = beta;
 284   Vx(3) = -beta2;              Vy(3) = 0;
 285   Vx(4) = 0;                   Vy(4) = -beta3;
 286   Vx(5) = Cos_b;               Vy(5) = Sin_b;
 287   Vx(6) = -beta*Sin_b;         Vy(6) = beta*Cos_b;
 288   Vx(7) = -beta2*Cos_b;        Vy(7) = -beta2*Sin_b;
 289   Vx(8) = beta3*Sin_b;         Vy(8) = -beta3*Cos_b;
 290
 291   Standard_Real b_bprim = bprim*beta,
 292                 b2_bprim = bprim*beta2,
 293                 b_bsecn  = bsecn*beta;
 294   DVx(1) = 0;                  DVy(1) = 0;
 295   DVx(2) = 0;                  DVy(2) = bprim;
 296   DVx(3) = -2*b_bprim;         DVy(3) = 0;
 297   DVx(4) = 0;                  DVy(4) = -3*b2_bprim;
 298   DVx(5) = -2*bprim*Sin_b;     DVy(5) = 2*bprim*Cos_b;
 299   DVx(6) = -bprim*Sin_b - 2*b_bprim*Cos_b;
 300                   DVy(6) = bprim*Cos_b - 2*b_bprim*Sin_b;
 301   DVx(7) = 2*b_bprim*(-Cos_b + beta*Sin_b);
 302                   DVy(7) = -2*b_bprim*(Sin_b + beta*Cos_b);
 303   DVx(8) = b2_bprim*(3*Sin_b + 2*beta*Cos_b);
 304                   DVy(8) = b2_bprim*(2*beta*Sin_b - 3*Cos_b);
 305
 306
 307   D2Vx(1) = 0;                  D2Vy(1) = 0;
 308   D2Vx(2) = 0;                  D2Vy(2) = bsecn;
 309   D2Vx(3) = -2*(bprim2+b_bsecn); D2Vy(3) = 0;
 310   D2Vx(4) = 0;    D2Vy(4) = -3*beta*(2*bprim2+b_bsecn);
 311   D2Vx(5) = -2*(bsecn*Sin_b + 2*bprim2*Cos_b);
 312                   D2Vy(5) = 2*(bsecn*Cos_b - 2*bprim2*Sin_b);
 313   D2Vx(6) = (4*beta*bprim2-bsecn)*Sin_b - 2*(2*bprim2+b_bsecn)*Cos_b;
 314
 315                   D2Vy(6) = (bsecn - 4*beta*bprim2)*Cos_b
 316                           - 2*(b_bsecn + 2*bprim2)*Sin_b;
 317
 318   aux =  2*(bprim2+b_bsecn);
 319   D2Vx(7) = aux*(-Cos_b + beta*Sin_b)
 320           + 2*beta*bprim2*(2*beta*Cos_b + 3*Sin_b);
 321
 322                   D2Vy(7) = -aux*(Sin_b + beta*Cos_b)
 323                           - 2*beta*bprim2*(3*Cos_b - 2*beta*Sin_b);
 324
 325   aux = beta*(2*bprim2+b_bsecn);
 326   D2Vx(8) = aux * (3*Sin_b + 2*beta*Cos_b)
 327           + 4*beta2*bprim2 * (2*Cos_b - beta*Sin_b);
 328
 329                   D2Vy(8)= aux * (2*beta*Sin_b - 3*Cos_b)
 330                          + 4*beta2*bprim2 * (2*Sin_b + beta*Cos_b);
 331
 332
 333   // Calcul des poles
 334   Px = BH * Vx;
 335   Py = BH * Vy;
 336   DPx = BH * DVx;
 337   DPy = BH * DVy;
 338   D2Px = BH * D2Vx;
 339   D2Py = BH * D2Vy;
 340
 341   gp_XYZ P, DP, D2P, auxyz;
 342   for (ii=1; ii<=Ordre; ii++) {
 343      P.SetCoord(Px(ii), Py(ii), 0);
 344      DP.SetCoord(DPx(ii), DPy(ii), 0);
 345      D2P.SetCoord(D2Px(ii), D2Py(ii), 0);
 346
 347      Poles(ii).ChangeCoord() = M*P + Center.XYZ();
 348      auxyz.SetLinearForm(1, MPrim*P,
 349                          1, M*DP,
 350                          DCenter.XYZ());
 351      DPoles(ii).SetXYZ(auxyz);
 352      P  *= MSecn;
 353      DP *= MPrim;
 354      D2P*= M;
 355      auxyz.SetLinearForm(1, P,
 356                          2, DP,
 357                          1, D2P,
 358                          D2Center.XYZ());
 359      D2Poles(ii).SetXYZ(auxyz);
 360   }
 361 }
 362