src/GProp/GProp_SGProps.gxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
   3 //
   4 // The content of this file is subject to the Open CASCADE Technology Public
   5 // License Version 6.5 (the "License"). You may not use the content of this file
   6 // except in compliance with the License. Please obtain a copy of the License
   7 // at http://www.opencascade.org and read it completely before using this file.
   8 //
   9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  11 //
  12 // The Original Code and all software distributed under the License is
  13 // distributed on an "AS IS" basis, without warranty of any kind, and the
  14 // Initial Developer hereby disclaims all such warranties, including without
  15 // limitation, any warranties of merchantability, fitness for a particular
  16 // purpose or non-infringement. Please see the License for the specific terms
  17 // and conditions governing the rights and limitations under the License.
  18
  19 #include <Standard_NotImplemented.hxx>
  20 #include <math_Vector.hxx>
  21 #include <math.hxx>
  22 #include <gp_Pnt2d.hxx>
  23 #include <gp_Vec2d.hxx>
  24 #include <gp_Pnt.hxx>
  25 #include <gp_Vec.hxx>
  26
  27 #include <TColStd_Array1OfReal.hxx>
  28 #include <Precision.hxx>
  29
  30 class HMath_Vector{
  31   math_Vector *pvec;
  32   void operator=(const math_Vector&){}
  33  public:
  34   HMath_Vector(){ pvec = 0;}
  35   HMath_Vector(math_Vector* pv){ pvec = pv;}
  36   ~HMath_Vector(){ if(pvec != 0) delete pvec;}
  37   void operator=(math_Vector* pv){ if(pvec != pv && pvec != 0) delete pvec;  pvec = pv;}
  38   Standard_Real& operator()(Standard_Integer i){ return (*pvec).operator()(i);}
  39   const Standard_Real& operator()(Standard_Integer i) const{ return (*pvec).operator()(i);}
  40   const math_Vector* operator->() const{ return pvec;}
  41   math_Vector* operator->(){ return pvec;}
  42   math_Vector* Vector(){ return pvec;}
  43   math_Vector* Init(Standard_Real v, Standard_Integer i = 0, Standard_Integer iEnd = 0){
  44     if(pvec == 0) return pvec;
  45     if(iEnd - i == 0) pvec->Init(v);
  46     else { Standard_Integer End = (iEnd <= pvec->Upper()) ? iEnd : pvec->Upper();
  47            for(; i <= End; i++) pvec->operator()(i) = v; }
  48     return pvec;
  49   }
  50 };
  51
  52 static Standard_Real EPS_PARAM          = 1.e-12;
  53 static Standard_Real EPS_DIM            = 1.e-20;
  54 static Standard_Real ERROR_ALGEBR_RATIO = 2.0/3.0;
  55
  56 static Standard_Integer  GPM        = 61;
  57 static Standard_Integer  SUBS_POWER = 32;
  58 static Standard_Integer  SM         = 1953;
  59
  60 static math_Vector LGaussP0(1,GPM);
  61 static math_Vector LGaussW0(1,GPM);
  62 static math_Vector LGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
  63 static math_Vector LGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
  64
  65 static math_Vector* LGaussP[] = {&LGaussP0,&LGaussP1};
  66 static math_Vector* LGaussW[] = {&LGaussW0,&LGaussW1};
  67
  68 static HMath_Vector L1    = new math_Vector(1,SM,0.0);
  69 static HMath_Vector L2    = new math_Vector(1,SM,0.0);
  70 static HMath_Vector DimL  = new math_Vector(1,SM,0.0);
  71 static HMath_Vector ErrL  = new math_Vector(1,SM,0.0);
  72 static HMath_Vector ErrUL = new math_Vector(1,SM,0.0);
  73 static HMath_Vector IxL   = new math_Vector(1,SM,0.0);
  74 static HMath_Vector IyL   = new math_Vector(1,SM,0.0);
  75 static HMath_Vector IzL   = new math_Vector(1,SM,0.0);
  76 static HMath_Vector IxxL  = new math_Vector(1,SM,0.0);
  77 static HMath_Vector IyyL  = new math_Vector(1,SM,0.0);
  78 static HMath_Vector IzzL  = new math_Vector(1,SM,0.0);
  79 static HMath_Vector IxyL  = new math_Vector(1,SM,0.0);
  80 static HMath_Vector IxzL  = new math_Vector(1,SM,0.0);
  81 static HMath_Vector IyzL  = new math_Vector(1,SM,0.0);
  82
  83 static math_Vector UGaussP0(1,GPM);
  84 static math_Vector UGaussW0(1,GPM);
  85 static math_Vector UGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
  86 static math_Vector UGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
  87
  88 static math_Vector* UGaussP[] = {&UGaussP0,&UGaussP1};
  89 static math_Vector* UGaussW[] = {&UGaussW0,&UGaussW1};
  90
  91 static HMath_Vector U1   = new math_Vector(1,SM,0.0);
  92 static HMath_Vector U2   = new math_Vector(1,SM,0.0);
  93 static HMath_Vector DimU = new math_Vector(1,SM,0.0);
  94 static HMath_Vector ErrU = new math_Vector(1,SM,0.0);
  95 static HMath_Vector IxU  = new math_Vector(1,SM,0.0);
  96 static HMath_Vector IyU  = new math_Vector(1,SM,0.0);
  97 static HMath_Vector IzU  = new math_Vector(1,SM,0.0);
  98 static HMath_Vector IxxU = new math_Vector(1,SM,0.0);
  99 static HMath_Vector IyyU = new math_Vector(1,SM,0.0);
 100 static HMath_Vector IzzU = new math_Vector(1,SM,0.0);
 101 static HMath_Vector IxyU = new math_Vector(1,SM,0.0);
 102 static HMath_Vector IxzU = new math_Vector(1,SM,0.0);
 103 static HMath_Vector IyzU = new math_Vector(1,SM,0.0);
 104
 105 static Standard_Integer FillIntervalBounds(Standard_Real               A,
 106                                            Standard_Real               B,
 107                                            const TColStd_Array1OfReal& Knots,
 108                                            HMath_Vector&               VA,
 109                                            HMath_Vector&               VB)
 110 {
 111   Standard_Integer i = 1, iEnd = Knots.Upper(), j = 1, k = 1;
 112   VA(j++) = A;
 113   for(; i <= iEnd; i++){
 114     Standard_Real kn = Knots(i);
 115     if(A < kn)
 116     {
 117       if(kn < B)
 118       {
 119         VA(j++) = VB(k++) = kn;
 120       }
 121       else
 122       {
 123         break;
 124       }
 125     }
 126   }
 127   VB(k) = B;
 128   return k;
 129 }
 130
 131 static inline Standard_Integer MaxSubs(Standard_Integer n, Standard_Integer coeff = SUBS_POWER){
 132 //  return n = IntegerLast()/coeff < n? IntegerLast(): n*coeff + 1;
 133   return Min((n * coeff + 1),SM);
 134 }
 135
 136 static Standard_Integer LFillIntervalBounds(Standard_Real               A,
 137                                             Standard_Real               B,
 138                                             const TColStd_Array1OfReal& Knots,
 139                                             const Standard_Integer      NumSubs)
 140 {
 141   Standard_Integer iEnd = Knots.Upper(), jEnd = L1->Upper();
 142   if(iEnd - 1 > jEnd){
 143     iEnd = MaxSubs(iEnd-1,NumSubs);
 144     L1    = new math_Vector(1,iEnd);
 145     L2    = new math_Vector(1,iEnd);
 146     DimL  = new math_Vector(1,iEnd);
 147     ErrL  = new math_Vector(1,iEnd,0.0);
 148     ErrUL = new math_Vector(1,iEnd,0.0);
 149     IxL   = new math_Vector(1,iEnd);
 150     IyL   = new math_Vector(1,iEnd);
 151     IzL   = new math_Vector(1,iEnd);
 152     IxxL  = new math_Vector(1,iEnd);
 153     IyyL  = new math_Vector(1,iEnd);
 154     IzzL  = new math_Vector(1,iEnd);
 155     IxyL  = new math_Vector(1,iEnd);
 156     IxzL  = new math_Vector(1,iEnd);
 157     IyzL  = new math_Vector(1,iEnd);
 158   }
 159   return FillIntervalBounds(A, B, Knots, L1, L2);
 160 }
 161
 162 static Standard_Integer UFillIntervalBounds(Standard_Real               A,
 163                                             Standard_Real               B,
 164                                             const TColStd_Array1OfReal& Knots,
 165                                             const Standard_Integer      NumSubs)
 166 {
 167   Standard_Integer iEnd = Knots.Upper(), jEnd = U1->Upper();
 168   if(iEnd - 1 > jEnd){
 169     iEnd = MaxSubs(iEnd-1,NumSubs);
 170     U1   = new math_Vector(1,iEnd);
 171     U2   = new math_Vector(1,iEnd);
 172     DimU = new math_Vector(1,iEnd);
 173     ErrU = new math_Vector(1,iEnd,0.0);
 174     IxU  = new math_Vector(1,iEnd);
 175     IyU  = new math_Vector(1,iEnd);
 176     IzU  = new math_Vector(1,iEnd);
 177     IxxU = new math_Vector(1,iEnd);
 178     IyyU = new math_Vector(1,iEnd);
 179     IzzU = new math_Vector(1,iEnd);
 180     IxyU = new math_Vector(1,iEnd);
 181     IxzU = new math_Vector(1,iEnd);
 182     IyzU = new math_Vector(1,iEnd);
 183   }
 184   return FillIntervalBounds(A, B, Knots, U1, U2);
 185 }
 186
 187 static Standard_Real CCompute(Face&                  S,
 188                               Domain&                D,
 189                               const gp_Pnt&          loc,
 190                               Standard_Real&         Dim,
 191                               gp_Pnt&                g,
 192                               gp_Mat&                inertia,
 193                               const Standard_Real    EpsDim,
 194                               const Standard_Boolean isErrorCalculation,
 195                               const Standard_Boolean isVerifyComputation)
 196 {
 197   Standard_Boolean isNaturalRestriction = S.NaturalRestriction();
 198
 199   Standard_Integer NumSubs = SUBS_POWER;
 200
 201   Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
 202   Dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
 203   Standard_Real x, y, z;
 204   //boundary curve parametrization
 205   Standard_Real l1, l2, lm, lr, l;
 206   //Face parametrization in U and V direction
 207   Standard_Real BV1, BV2, v;
 208   Standard_Real BU1, BU2, u1, u2, um, ur, u;
 209   S.Bounds (BU1, BU2, BV1, BV2);  u1 = BU1;
 210   //location point used to compute the inertia
 211   Standard_Real xloc, yloc, zloc;
 212   loc.Coord (xloc, yloc, zloc); // use member of parent class
 213   //Jacobien (x, y, z) -> (u, v) = ||n||
 214   Standard_Real ds;
 215   //On the Face
 216   gp_Pnt Ps;
 217   gp_Vec VNor;
 218   //On the boundary curve u-v
 219   gp_Pnt2d Puv;
 220   gp_Vec2d Vuv;
 221   Standard_Real Dul;  // Dul = Du / Dl
 222   Standard_Real CDim[2], CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz;
 223   Standard_Real LocDim[2], LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy, LocIxz, LocIyz;
 224
 225   Standard_Real ErrorU, ErrorL, ErrorLMax = 0.0, Eps=0.0, EpsL=0.0, EpsU=0.0;
 226
 227   Standard_Integer iD = 0, NbLSubs, iLS, iLSubEnd, iGL, iGLEnd, NbLGaussP[2], LRange[2], iL, kL, kLEnd, IL, JL;
 228   Standard_Integer i, NbUSubs, iUS, iUSubEnd, iGU, iGUEnd, NbUGaussP[2], URange[2], iU, kU, kUEnd, IU, JU;
 229   Standard_Integer UMaxSubs, LMaxSubs;
 230   iGLEnd = isErrorCalculation? 2: 1;
 231   for(i = 0; i < 2; i++) {
 232     LocDim[i] = 0.0;
 233     CDim[i] = 0.0;
 234   }
 235
 236   NbUGaussP[0] = S.SIntOrder(EpsDim);
 237   NbUGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbUGaussP[0]));
 238   math::GaussPoints(NbUGaussP[0],UGaussP0);  math::GaussWeights(NbUGaussP[0],UGaussW0);
 239   math::GaussPoints(NbUGaussP[1],UGaussP1);  math::GaussWeights(NbUGaussP[1],UGaussW1);
 240
 241   NbUSubs = S.SUIntSubs();
 242   TColStd_Array1OfReal UKnots(1,NbUSubs+1);
 243   S.UKnots(UKnots);
 244
 245
 246   while (isNaturalRestriction || D.More()) {
 247     if(isNaturalRestriction){
 248       NbLGaussP[0] = Min(2*NbUGaussP[0],math::GaussPointsMax());
 249     }else{
 250       S.Load(D.Value());  ++iD;
 251       NbLGaussP[0] = S.LIntOrder(EpsDim);
 252     }
 253
 254
 255     NbLGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbLGaussP[0]));
 256     math::GaussPoints(NbLGaussP[0],LGaussP0);  math::GaussWeights(NbLGaussP[0],LGaussW0);
 257     math::GaussPoints(NbLGaussP[1],LGaussP1);  math::GaussWeights(NbLGaussP[1],LGaussW1);
 258
 259     NbLSubs = isNaturalRestriction? S.SVIntSubs(): S.LIntSubs();
 260
 261     TColStd_Array1OfReal LKnots(1,NbLSubs+1);
 262     if(isNaturalRestriction){
 263       S.VKnots(LKnots);
 264       l1 = BV1; l2 = BV2;
 265     }else{
 266       S.LKnots(LKnots);
 267       l1 = S.FirstParameter();  l2 = S.LastParameter();
 268     }
 269     ErrorL = 0.0;
 270     kLEnd = 1; JL = 0;
 271     //OCC503(apo): if(Abs(l2-l1) < EPS_PARAM) continue;
 272     if(Abs(l2-l1) > EPS_PARAM) {
 273       iLSubEnd = LFillIntervalBounds(l1, l2, LKnots, NumSubs);
 274       LMaxSubs = MaxSubs(iLSubEnd);
 275       if(LMaxSubs > DimL.Vector()->Upper()) LMaxSubs = DimL.Vector()->Upper();
 276       DimL.Init(0.0,1,LMaxSubs);  ErrL.Init(0.0,1,LMaxSubs);  ErrUL.Init(0.0,1,LMaxSubs);
 277       do{// while: L
 278         if(++JL > iLSubEnd){
 279           LRange[0] = IL = ErrL->Max();  LRange[1] = JL;
 280           L1(JL) = (L1(IL) + L2(IL))/2.0;  L2(JL) = L2(IL);  L2(IL) = L1(JL);
 281         }else  LRange[0] = IL = JL;
 282         if(JL == LMaxSubs || Abs(L2(JL) - L1(JL)) < EPS_PARAM)
 283         if(kLEnd == 1){
 284           DimL(JL) = ErrL(JL) = IxL(JL) = IyL(JL) = IzL(JL) =
 285             IxxL(JL) = IyyL(JL) = IzzL(JL) = IxyL(JL) = IxzL(JL) = IyzL(JL) = 0.0;
 286         }else{
 287           JL--;
 288           EpsL = ErrorL;  Eps = EpsL/0.9;
 289           break;
 290         }
 291         else
 292           for(kL=0; kL < kLEnd; kL++){
 293             iLS = LRange[kL];
 294             lm = 0.5*(L2(iLS) + L1(iLS));
 295             lr = 0.5*(L2(iLS) - L1(iLS));
 296             CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
 297             for(iGL=0; iGL < iGLEnd; iGL++){//
 298               CDim[iGL] = 0.0;
 299               for(iL=1; iL<=NbLGaussP[iGL]; iL++){
 300                 l = lm + lr*(*LGaussP[iGL])(iL);
 301                 if(isNaturalRestriction){
 302                   v = l; u2 = BU2; Dul = (*LGaussW[iGL])(iL);
 303                 }else{
 304                   S.D12d (l, Puv, Vuv);
 305                   Dul = Vuv.Y()*(*LGaussW[iGL])(iL);  // Dul = Du / Dl
 306                   if(Abs(Dul) < EPS_PARAM) continue;
 307                   v  = Puv.Y();  u2 = Puv.X();
 308                   //Check on cause out off bounds of value current parameter
 309                   if(v < BV1) v = BV1; else if(v > BV2) v = BV2;
 310                   if(u2 < BU1) u2 = BU1; else if(u2 > BU2) u2 = BU2;
 311                 }
 312                 ErrUL(iLS) = 0.0;
 313                 kUEnd = 1; JU = 0;
 314                 if(Abs(u2-u1) < EPS_PARAM) continue;
 315                 iUSubEnd = UFillIntervalBounds(u1, u2, UKnots, NumSubs);
 316                 UMaxSubs = MaxSubs(iUSubEnd);
 317                 if(UMaxSubs > DimU.Vector()->Upper()) UMaxSubs = DimU.Vector()->Upper();
 318                 DimU.Init(0.0,1,UMaxSubs);  ErrU.Init(0.0,1,UMaxSubs);  ErrorU = 0.0;
 319                 do{//while: U
 320                   if(++JU > iUSubEnd){
 321                     URange[0] = IU = ErrU->Max();  URange[1] = JU;
 322                     U1(JU) = (U1(IU)+U2(IU))/2.0;  U2(JU) = U2(IU);  U2(IU) = U1(JU);
 323                   }else  URange[0] = IU = JU;
 324                   if(JU == UMaxSubs || Abs(U2(JU) - U1(JU)) < EPS_PARAM)
 325                     if(kUEnd == 1){
 326                       DimU(JU) = ErrU(JU) = IxU(JU) = IyU(JU) = IzU(JU) =
 327                         IxxU(JU) = IyyU(JU) = IzzU(JU) = IxyU(JU) = IxzU(JU) = IyzU(JU) = 0.0;
 328                     }else{
 329                       JU--;
 330                       EpsU = ErrorU;  Eps = EpsU*Abs((u2-u1)*Dul)/0.1;  EpsL = 0.9*Eps;
 331                       break;
 332                     }
 333                   else
 334                     for(kU=0; kU < kUEnd; kU++){
 335                       iUS = URange[kU];
 336                       um = 0.5*(U2(iUS) + U1(iUS));
 337                       ur = 0.5*(U2(iUS) - U1(iUS));
 338                       LocIx = LocIy = LocIz = LocIxx = LocIyy = LocIzz = LocIxy = LocIxz = LocIyz = 0.0;
 339                       iGUEnd = iGLEnd - iGL;
 340                       for(iGU=0; iGU < iGUEnd; iGU++){//
 341                         LocDim[iGU] = 0.0;
 342                         for(iU=1; iU<=NbUGaussP[iGU]; iU++){
 343                           u = um + ur*(*UGaussP[iGU])(iU);
 344                           S.Normal(u, v, Ps, VNor);
 345                           ds = VNor.Magnitude();    //Jacobien(x,y,z) -> (u,v)=||n||
 346                           ds *= (*UGaussW[iGU])(iU);
 347                           LocDim[iGU] += ds;
 348                           if(iGU > 0) continue;
 349                           Ps.Coord(x, y, z);
 350                           x -= xloc;  y -= yloc;  z -= zloc;
 351                           LocIx += x*ds;  LocIy += y*ds;  LocIz += z*ds;
 352                           LocIxy += x*y*ds;  LocIyz += y*z*ds;  LocIxz += x*z*ds;
 353                           x *= x;  y *= y;  z *= z;
 354                           LocIxx += (y+z)*ds;  LocIyy += (x+z)*ds;  LocIzz += (x+y)*ds;
 355                         }//for: iU
 356                       }//for: iGU
 357                       DimU(iUS) = LocDim[0]*ur;
 358                       if(iGL > 0) continue;
 359                       ErrU(iUS) = Abs(LocDim[1]-LocDim[0])*ur;
 360                       IxU(iUS) = LocIx*ur; IyU(iUS) = LocIy*ur; IzU(iUS) = LocIz*ur;
 361                       IxxU(iUS) = LocIxx*ur; IyyU(iUS) = LocIyy*ur; IzzU(iUS) = LocIzz*ur;
 362                       IxyU(iUS) = LocIxy*ur; IxzU(iUS) = LocIxz*ur; IyzU(iUS) = LocIyz*ur;
 363                     }//for: kU (iUS)
 364                   if(JU == iUSubEnd)  kUEnd = 2;
 365                   if(kUEnd == 2)  ErrorU = ErrU(ErrU->Max());
 366                 }while((ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1);
 367                 for(i=1; i<=JU; i++)  CDim[iGL] += DimU(i)*Dul;
 368                 if(iGL > 0) continue;
 369                 ErrUL(iLS) = ErrorU*Abs((u2-u1)*Dul);
 370                 for(i=1; i<=JU; i++){
 371                   CIx += IxU(i)*Dul; CIy += IyU(i)*Dul; CIz += IzU(i)*Dul;
 372                   CIxx += IxxU(i)*Dul; CIyy += IyyU(i)*Dul; CIzz += IzzU(i)*Dul;
 373                   CIxy += IxyU(i)*Dul; CIxz += IxzU(i)*Dul; CIyz += IyzU(i)*Dul;
 374                 }
 375               }//for: iL
 376             }//for: iGL
 377             DimL(iLS) = CDim[0]*lr;
 378             if(iGLEnd == 2) ErrL(iLS) = Abs(CDim[1]-CDim[0])*lr + ErrUL(iLS);
 379             IxL(iLS) = CIx*lr; IyL(iLS) = CIy*lr; IzL(iLS) = CIz*lr;
 380             IxxL(iLS) = CIxx*lr; IyyL(iLS) = CIyy*lr; IzzL(iLS) = CIzz*lr;
 381             IxyL(iLS) = CIxy*lr; IxzL(iLS) = CIxz*lr; IyzL(iLS) = CIyz*lr;
 382           }//for: (kL)iLS
 383         //  Calculate/correct epsilon of computation by current value of Dim
 384         //That is need for not spend time for
 385         if(JL == iLSubEnd){
 386           kLEnd = 2;
 387           Standard_Real DDim = 0.0;
 388           for(i=1; i<=JL; i++) DDim += DimL(i);
 389           DDim = Abs(DDim*EpsDim);
 390           if(DDim > Eps) {
 391             Eps = DDim;  EpsL = 0.9*Eps;
 392           }
 393         }
 394         if(kLEnd == 2) ErrorL = ErrL(ErrL->Max());
 395       }while((ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1);
 396       for(i=1; i<=JL; i++){
 397         Dim += DimL(i);
 398         Ix += IxL(i); Iy += IyL(i); Iz += IzL(i);
 399         Ixx += IxxL(i); Iyy += IyyL(i); Izz += IzzL(i);
 400         Ixy += IxyL(i); Ixz += IxzL(i); Iyz += IyzL(i);
 401       }
 402       ErrorLMax = Max(ErrorLMax, ErrorL);
 403     }
 404     if(isNaturalRestriction) break;
 405     D.Next();
 406   }
 407   if(Abs(Dim) >= EPS_DIM){
 408     Ix /= Dim;  Iy /= Dim;  Iz /= Dim;
 409     g.SetCoord (Ix, Iy, Iz);
 410   }else{
 411     Dim =0.0;
 412     g.SetCoord (0., 0.,0.);
 413   }
 414   inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
 415                     gp_XYZ (-Ixy, Iyy, -Iyz),
 416                     gp_XYZ (-Ixz, -Iyz, Izz));
 417
 418   if(iGLEnd == 2) Eps = Dim != 0.0? ErrorLMax/Abs(Dim): 0.0;
 419   else Eps = EpsDim;
 420   return Eps;
 421 }
 422
 423 static Standard_Real Compute(Face& S, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia,
 424                              Standard_Real EpsDim)
 425 {
 426   Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
 427   Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
 428   EpsDim = Abs(EpsDim);
 429   Domain D;
 430   return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
 431 }
 432
 433 static Standard_Real Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia,
 434                              Standard_Real EpsDim)
 435 {
 436   Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
 437   Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
 438   EpsDim = Abs(EpsDim);
 439   return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
 440 }
 441
 442 static void Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& dim, gp_Pnt& g, gp_Mat& inertia){
 443    Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
 444    dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
 445
 446    Standard_Real x, y, z;
 447    Standard_Integer NbCGaussgp_Pnts = 0;
 448
 449    Standard_Real l1, l2, lm, lr, l;   //boundary curve parametrization
 450    Standard_Real v1, v2, vm, vr, v;   //Face parametrization in v direction
 451    Standard_Real u1, u2, um, ur, u;
 452    Standard_Real ds;                  //Jacobien (x, y, z) -> (u, v) = ||n||
 453
 454    gp_Pnt P;                    //On the Face
 455    gp_Vec VNor;
 456
 457    gp_Pnt2d Puv;                //On the boundary curve u-v
 458    gp_Vec2d Vuv;
 459    Standard_Real Dul;                 // Dul = Du / Dl
 460    Standard_Real CArea, CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz;
 461    Standard_Real LocArea, LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy,
 462         LocIxz, LocIyz;
 463
 464
 465    S.Bounds (u1, u2, v1, v2);
 466
 467    Standard_Integer NbUGaussgp_Pnts = Min(S.UIntegrationOrder (),
 468                                           math::GaussPointsMax());
 469    Standard_Integer NbVGaussgp_Pnts = Min(S.VIntegrationOrder (),
 470                                           math::GaussPointsMax());
 471
 472    Standard_Integer NbGaussgp_Pnts = Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts);
 473
 474    //Number of Gauss points for the integration
 475    //on the Face
 476    math_Vector GaussSPV (1, NbGaussgp_Pnts);
 477    math_Vector GaussSWV (1, NbGaussgp_Pnts);
 478    math::GaussPoints  (NbGaussgp_Pnts,GaussSPV);
 479    math::GaussWeights (NbGaussgp_Pnts,GaussSWV);
 480
 481
 482    //location point used to compute the inertia
 483    Standard_Real xloc, yloc, zloc;
 484    loc.Coord (xloc, yloc, zloc);
 485
 486    while (D.More()) {
 487
 488       S.Load(D.Value());
 489       NbCGaussgp_Pnts =  Min(S.IntegrationOrder (), math::GaussPointsMax());
 490
 491       math_Vector GaussCP (1, NbCGaussgp_Pnts);
 492       math_Vector GaussCW (1, NbCGaussgp_Pnts);
 493       math::GaussPoints  (NbCGaussgp_Pnts,GaussCP);
 494       math::GaussWeights (NbCGaussgp_Pnts,GaussCW);
 495
 496       CArea = 0.0;
 497       CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
 498       l1 = S.FirstParameter ();
 499       l2 = S.LastParameter  ();
 500       lm = 0.5 * (l2 + l1);
 501       lr = 0.5 * (l2 - l1);
 502
 503       Puv = S.Value2d (lm);
 504       vm  = Puv.Y();
 505       Puv = S.Value2d (lr);
 506       vr  = Puv.Y();
 507
 508       for (Standard_Integer i = 1; i <= NbCGaussgp_Pnts; i++) {
 509         l = lm + lr * GaussCP (i);
 510         S.D12d(l, Puv, Vuv);
 511         v   = Puv.Y();
 512         u2  = Puv.X();
 513         Dul = Vuv.Y();
 514         Dul *= GaussCW (i);
 515         um  = 0.5 * (u2 + u1);
 516         ur  = 0.5 * (u2 - u1);
 517         LocArea = LocIx  = LocIy  = LocIz = LocIxx = LocIyy = LocIzz =
 518         LocIxy  = LocIxz = LocIyz = 0.0;
 519         for (Standard_Integer j = 1; j <= NbGaussgp_Pnts; j++) {
 520           u = um + ur * GaussSPV (j);
 521           S.Normal (u, v, P, VNor);
 522           ds = VNor.Magnitude();    //normal.Magnitude
 523           ds = ds * Dul * GaussSWV (j);
 524           LocArea +=  ds;
 525           P.Coord (x, y, z);
 526           x      -= xloc;
 527           y      -= yloc;
 528           z      -= zloc;
 529           LocIx  += x * ds;
 530           LocIy  += y * ds;
 531           LocIz  += z * ds;
 532           LocIxy += x * y * ds;
 533           LocIyz += y * z * ds;
 534           LocIxz += x * z * ds;
 535           x *= x;
 536           y *= y;
 537           z *= z;
 538           LocIxx += (y + z) * ds;
 539           LocIyy += (x + z) * ds;
 540           LocIzz += (x + y) * ds;
 541         }
 542         CArea += LocArea * ur;
 543         CIx   += LocIx * ur;
 544         CIy   += LocIy * ur;
 545         CIz   += LocIz * ur;
 546         CIxx  += LocIxx * ur;
 547         CIyy  += LocIyy * ur;
 548         CIzz  += LocIzz * ur;
 549         CIxy  += LocIxy * ur;
 550         CIxz  += LocIxz * ur;
 551         CIyz  += LocIyz * ur;
 552       }
 553       dim += CArea * lr;
 554       Ix  += CIx * lr;
 555       Iy  += CIy * lr;
 556       Iz  += CIz * lr;
 557       Ixx += CIxx * lr;
 558       Iyy += CIyy * lr;
 559       Izz += CIzz * lr;
 560       Ixy += CIxy * lr;
 561       Ixz += CIxz * lr;
 562       Iyz += CIyz * lr;
 563       D.Next();
 564    }
 565    if (Abs(dim) >= EPS_DIM) {
 566      Ix /= dim;
 567      Iy /= dim;
 568      Iz /= dim;
 569      g.SetCoord (Ix, Iy, Iz);
 570    }
 571    else {
 572      dim =0.;
 573      g.SetCoord (0., 0.,0.);
 574    }
 575    inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
 576                      gp_XYZ (-Ixy, Iyy, -Iyz),
 577                      gp_XYZ (-Ixz, -Iyz, Izz));
 578 }
 579
 580
 581
 582 static void Compute(const Face& S, const gp_Pnt& loc, Standard_Real& dim, gp_Pnt& g, gp_Mat& inertia){
 583    Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
 584    dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
 585
 586    Standard_Real LowerU, UpperU, LowerV, UpperV;
 587    S.Bounds (LowerU, UpperU, LowerV, UpperV);
 588    Standard_Integer UOrder = Min(S.UIntegrationOrder (),
 589                                  math::GaussPointsMax());
 590    Standard_Integer VOrder = Min(S.VIntegrationOrder (),
 591                                  math::GaussPointsMax());
 592    gp_Pnt P;
 593    gp_Vec VNor;
 594    Standard_Real dsi, ds;
 595    Standard_Real ur, um, u, vr, vm, v;
 596    Standard_Real x, y, z;
 597    Standard_Real Ixi, Iyi, Izi, Ixxi, Iyyi, Izzi, Ixyi, Ixzi, Iyzi;
 598    Standard_Real xloc, yloc, zloc;
 599    loc.Coord (xloc, yloc, zloc);
 600
 601    Standard_Integer i, j;
 602    math_Vector GaussPU (1, UOrder);     //gauss points and weights
 603    math_Vector GaussWU (1, UOrder);
 604    math_Vector GaussPV (1, VOrder);
 605    math_Vector GaussWV (1, VOrder);
 606
 607    //Recuperation des points de Gauss dans le fichier GaussPoints.
 608    math::GaussPoints  (UOrder,GaussPU);
 609    math::GaussWeights (UOrder,GaussWU);
 610    math::GaussPoints  (VOrder,GaussPV);
 611    math::GaussWeights (VOrder,GaussWV);
 612
 613    // Calcul des integrales aux points de gauss :
 614    um = 0.5 * (UpperU + LowerU);
 615    vm = 0.5 * (UpperV + LowerV);
 616    ur = 0.5 * (UpperU - LowerU);
 617    vr = 0.5 * (UpperV - LowerV);
 618
 619    for (j = 1; j <= VOrder; j++) {
 620      v = vm + vr * GaussPV (j);
 621      dsi = Ixi = Iyi = Izi = Ixxi = Iyyi = Izzi = Ixyi = Ixzi = Iyzi = 0.0;
 622
 623      for (i = 1; i <= UOrder; i++) {
 624        u = um + ur * GaussPU (i);
 625        S.Normal (u, v, P, VNor);
 626        ds = VNor.Magnitude() * GaussWU (i);
 627        P.Coord (x, y, z);
 628        x    -=  xloc;
 629        y    -=  yloc;
 630        z    -=  zloc;
 631        dsi  +=  ds;
 632        Ixi  += x * ds;
 633        Iyi  += y * ds;
 634        Izi  += z * ds;
 635        Ixyi += x * y * ds;
 636        Iyzi += y * z * ds;
 637        Ixzi += x * z * ds;
 638        x    *= x;
 639        y    *= y;
 640        z    *= z;
 641        Ixxi += (y + z) * ds;
 642        Iyyi += (x + z) * ds;
 643        Izzi += (x + y) * ds;
 644      }
 645      dim  += dsi  * GaussWV (j);
 646      Ix    += Ixi  * GaussWV (j);
 647      Iy    += Iyi  * GaussWV (j);
 648      Iz    += Izi  * GaussWV (j);
 649      Ixx   += Ixxi * GaussWV (j);
 650      Iyy   += Iyyi * GaussWV (j);
 651      Izz   += Izzi * GaussWV (j);
 652      Ixy   += Ixyi * GaussWV (j);
 653      Iyz   += Iyzi * GaussWV (j);
 654      Ixz   += Ixzi * GaussWV (j);
 655    }
 656    vr    *= ur;
 657    Ixx   *= vr;
 658    Iyy   *= vr;
 659    Izz   *= vr;
 660    Ixy   *= vr;
 661    Ixz   *= vr;
 662    Iyz   *= vr;
 663    if (Abs(dim) >= EPS_DIM) {
 664      Ix    /= dim;
 665      Iy    /= dim;
 666      Iz    /= dim;
 667      dim   *= vr;
 668      g.SetCoord (Ix, Iy, Iz);
 669    }
 670    else {
 671      dim =0.;
 672      g.SetCoord (0.,0.,0.);
 673    }
 674    inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
 675                      gp_XYZ (-Ixy, Iyy, -Iyz),
 676                      gp_XYZ (-Ixz, -Iyz, Izz));
 677 }
 678
 679 GProp_SGProps::GProp_SGProps(){}
 680
 681 GProp_SGProps::GProp_SGProps (const Face&   S,
 682                               const gp_Pnt& SLocation
 683                               )
 684 {
 685    SetLocation(SLocation);
 686    Perform(S);
 687 }
 688
 689 GProp_SGProps::GProp_SGProps (Face&   S,
 690                               Domain& D,
 691                               const gp_Pnt& SLocation
 692                               )
 693 {
 694    SetLocation(SLocation);
 695    Perform(S,D);
 696 }
 697
 698 GProp_SGProps::GProp_SGProps(Face& S, const gp_Pnt& SLocation, const Standard_Real Eps){
 699   SetLocation(SLocation);
 700   Perform(S, Eps);
 701 }
 702
 703 GProp_SGProps::GProp_SGProps(Face& S, Domain& D, const gp_Pnt& SLocation, const Standard_Real Eps){
 704   SetLocation(SLocation);
 705   Perform(S, D, Eps);
 706 }
 707
 708 void GProp_SGProps::SetLocation(const gp_Pnt& SLocation){
 709   loc = SLocation;
 710 }
 711
 712 void GProp_SGProps::Perform(const Face& S){
 713   Compute(S,loc,dim,g,inertia);
 714   myEpsilon = 1.0;
 715   return;
 716 }
 717
 718 void GProp_SGProps::Perform(Face& S, Domain& D){
 719   Compute(S,D,loc,dim,g,inertia);
 720   myEpsilon = 1.0;
 721   return;
 722 }
 723
 724 Standard_Real GProp_SGProps::Perform(Face& S, const Standard_Real Eps){
 725   return myEpsilon = Compute(S,loc,dim,g,inertia,Eps);
 726 }
 727
 728 Standard_Real GProp_SGProps::Perform(Face& S, Domain& D, const Standard_Real Eps){
 729   return myEpsilon = Compute(S,D,loc,dim,g,inertia,Eps);
 730 }
 731
 732
 733 Standard_Real GProp_SGProps::GetEpsilon(){
 734   return myEpsilon;
 735 }