src/GProp/GProp_SGProps.gxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and / or modify it
   7 // under the terms of the GNU Lesser General Public version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 #include <Standard_NotImplemented.hxx>
  16 #include <math_Vector.hxx>
  17 #include <math.hxx>
  18 #include <gp_Pnt2d.hxx>
  19 #include <gp_Vec2d.hxx>
  20 #include <gp_Pnt.hxx>
  21 #include <gp_Vec.hxx>
  22
  23 #include <TColStd_Array1OfReal.hxx>
  24 #include <Precision.hxx>
  25
  26 class HMath_Vector{
  27   math_Vector *pvec;
  28   void operator=(const math_Vector&){}
  29  public:
  30   HMath_Vector(){ pvec = 0;}
  31   HMath_Vector(math_Vector* pv){ pvec = pv;}
  32   ~HMath_Vector(){ if(pvec != 0) delete pvec;}
  33   void operator=(math_Vector* pv){ if(pvec != pv && pvec != 0) delete pvec;  pvec = pv;}
  34   Standard_Real& operator()(Standard_Integer i){ return (*pvec).operator()(i);}
  35   const Standard_Real& operator()(Standard_Integer i) const{ return (*pvec).operator()(i);}
  36   const math_Vector* operator->() const{ return pvec;}
  37   math_Vector* operator->(){ return pvec;}
  38   math_Vector* Vector(){ return pvec;}
  39   math_Vector* Init(Standard_Real v, Standard_Integer i = 0, Standard_Integer iEnd = 0){
  40     if(pvec == 0) return pvec;
  41     if(iEnd - i == 0) pvec->Init(v);
  42     else { Standard_Integer End = (iEnd <= pvec->Upper()) ? iEnd : pvec->Upper();
  43            for(; i <= End; i++) pvec->operator()(i) = v; }
  44     return pvec;
  45   }
  46 };
  47
  48 static Standard_Real EPS_PARAM          = 1.e-12;
  49 static Standard_Real EPS_DIM            = 1.e-20;
  50 static Standard_Real ERROR_ALGEBR_RATIO = 2.0/3.0;
  51
  52 static Standard_Integer  GPM        = 61;
  53 static Standard_Integer  SUBS_POWER = 32;
  54 static Standard_Integer  SM         = 1953;
  55
  56 static math_Vector LGaussP0(1,GPM);
  57 static math_Vector LGaussW0(1,GPM);
  58 static math_Vector LGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
  59 static math_Vector LGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
  60
  61 static math_Vector* LGaussP[] = {&LGaussP0,&LGaussP1};
  62 static math_Vector* LGaussW[] = {&LGaussW0,&LGaussW1};
  63
  64 static HMath_Vector L1    = new math_Vector(1,SM,0.0);
  65 static HMath_Vector L2    = new math_Vector(1,SM,0.0);
  66 static HMath_Vector DimL  = new math_Vector(1,SM,0.0);
  67 static HMath_Vector ErrL  = new math_Vector(1,SM,0.0);
  68 static HMath_Vector ErrUL = new math_Vector(1,SM,0.0);
  69 static HMath_Vector IxL   = new math_Vector(1,SM,0.0);
  70 static HMath_Vector IyL   = new math_Vector(1,SM,0.0);
  71 static HMath_Vector IzL   = new math_Vector(1,SM,0.0);
  72 static HMath_Vector IxxL  = new math_Vector(1,SM,0.0);
  73 static HMath_Vector IyyL  = new math_Vector(1,SM,0.0);
  74 static HMath_Vector IzzL  = new math_Vector(1,SM,0.0);
  75 static HMath_Vector IxyL  = new math_Vector(1,SM,0.0);
  76 static HMath_Vector IxzL  = new math_Vector(1,SM,0.0);
  77 static HMath_Vector IyzL  = new math_Vector(1,SM,0.0);
  78
  79 static math_Vector UGaussP0(1,GPM);
  80 static math_Vector UGaussW0(1,GPM);
  81 static math_Vector UGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
  82 static math_Vector UGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
  83
  84 static math_Vector* UGaussP[] = {&UGaussP0,&UGaussP1};
  85 static math_Vector* UGaussW[] = {&UGaussW0,&UGaussW1};
  86
  87 static HMath_Vector U1   = new math_Vector(1,SM,0.0);
  88 static HMath_Vector U2   = new math_Vector(1,SM,0.0);
  89 static HMath_Vector DimU = new math_Vector(1,SM,0.0);
  90 static HMath_Vector ErrU = new math_Vector(1,SM,0.0);
  91 static HMath_Vector IxU  = new math_Vector(1,SM,0.0);
  92 static HMath_Vector IyU  = new math_Vector(1,SM,0.0);
  93 static HMath_Vector IzU  = new math_Vector(1,SM,0.0);
  94 static HMath_Vector IxxU = new math_Vector(1,SM,0.0);
  95 static HMath_Vector IyyU = new math_Vector(1,SM,0.0);
  96 static HMath_Vector IzzU = new math_Vector(1,SM,0.0);
  97 static HMath_Vector IxyU = new math_Vector(1,SM,0.0);
  98 static HMath_Vector IxzU = new math_Vector(1,SM,0.0);
  99 static HMath_Vector IyzU = new math_Vector(1,SM,0.0);
 100
 101 static inline Standard_Real MultiplicationInf(Standard_Real theMA, Standard_Real theMB)
 102 {
 103   if((theMA == 0.0) || (theMB == 0.0))  //strictly zerro (without any tolerances)
 104     return 0.0;
 105
 106   if(Precision::IsPositiveInfinite(theMA))
 107   {
 108     if(theMB < 0.0)
 109       return -Precision::Infinite();
 110     else
 111       return Precision::Infinite();
 112   }
 113
 114   if(Precision::IsPositiveInfinite(theMB))
 115   {
 116     if(theMA < 0.0)
 117       return -Precision::Infinite();
 118     else
 119       return Precision::Infinite();
 120   }
 121
 122   if(Precision::IsNegativeInfinite(theMA))
 123   {
 124     if(theMB < 0.0)
 125       return +Precision::Infinite();
 126     else
 127       return -Precision::Infinite();
 128   }
 129
 130   if(Precision::IsNegativeInfinite(theMB))
 131   {
 132     if(theMA < 0.0)
 133       return +Precision::Infinite();
 134     else
 135       return -Precision::Infinite();
 136   }
 137
 138   return (theMA * theMB);
 139 }
 140
 141 static inline Standard_Real AdditionInf(Standard_Real theMA, Standard_Real theMB)
 142 {
 143   if(Precision::IsPositiveInfinite(theMA))
 144   {
 145     if(Precision::IsNegativeInfinite(theMB))
 146       return 0.0;
 147     else
 148       return Precision::Infinite();
 149   }
 150
 151   if(Precision::IsPositiveInfinite(theMB))
 152   {
 153     if(Precision::IsNegativeInfinite(theMA))
 154       return 0.0;
 155     else
 156       return Precision::Infinite();
 157   }
 158
 159   if(Precision::IsNegativeInfinite(theMA))
 160   {
 161     if(Precision::IsPositiveInfinite(theMB))
 162       return 0.0;
 163     else
 164       return -Precision::Infinite();
 165   }
 166
 167   if(Precision::IsNegativeInfinite(theMB))
 168   {
 169     if(Precision::IsPositiveInfinite(theMA))
 170       return 0.0;
 171     else
 172       return -Precision::Infinite();
 173   }
 174
 175   return (theMA + theMB);
 176 }
 177
 178 static inline Standard_Real Multiplication(Standard_Real theMA, Standard_Real theMB)
 179 {
 180   return (theMA * theMB);
 181 }
 182
 183 static inline Standard_Real Addition(Standard_Real theMA, Standard_Real theMB)
 184 {
 185   return (theMA + theMB);
 186 }
 187
 188 static Standard_Integer FillIntervalBounds(Standard_Real               A,
 189                                            Standard_Real               B,
 190                                            const TColStd_Array1OfReal& Knots,
 191                                            HMath_Vector&               VA,
 192                                            HMath_Vector&               VB)
 193 {
 194   Standard_Integer i = 1, iEnd = Knots.Upper(), j = 1, k = 1;
 195   VA(j++) = A;
 196   for(; i <= iEnd; i++){
 197     Standard_Real kn = Knots(i);
 198     if(A < kn)
 199     {
 200       if(kn < B)
 201       {
 202         VA(j++) = VB(k++) = kn;
 203       }
 204       else
 205       {
 206         break;
 207       }
 208     }
 209   }
 210   VB(k) = B;
 211   return k;
 212 }
 213
 214 static inline Standard_Integer MaxSubs(Standard_Integer n, Standard_Integer coeff = SUBS_POWER){
 215   //  return n = IntegerLast()/coeff < n? IntegerLast(): n*coeff + 1;
 216   return Min((n * coeff + 1),SM);
 217 }
 218
 219 static Standard_Integer LFillIntervalBounds(Standard_Real               A,
 220                                             Standard_Real               B,
 221                                             const TColStd_Array1OfReal& Knots,
 222                                             const Standard_Integer      NumSubs)
 223 {
 224   Standard_Integer iEnd = Knots.Upper(), jEnd = L1->Upper();
 225   if(iEnd - 1 > jEnd){
 226     iEnd = MaxSubs(iEnd-1,NumSubs);
 227     L1    = new math_Vector(1,iEnd);
 228     L2    = new math_Vector(1,iEnd);
 229     DimL  = new math_Vector(1,iEnd);
 230     ErrL  = new math_Vector(1,iEnd,0.0);
 231     ErrUL = new math_Vector(1,iEnd,0.0);
 232     IxL   = new math_Vector(1,iEnd);
 233     IyL   = new math_Vector(1,iEnd);
 234     IzL   = new math_Vector(1,iEnd);
 235     IxxL  = new math_Vector(1,iEnd);
 236     IyyL  = new math_Vector(1,iEnd);
 237     IzzL  = new math_Vector(1,iEnd);
 238     IxyL  = new math_Vector(1,iEnd);
 239     IxzL  = new math_Vector(1,iEnd);
 240     IyzL  = new math_Vector(1,iEnd);
 241   }
 242   return FillIntervalBounds(A, B, Knots, L1, L2);
 243 }
 244
 245 static Standard_Integer UFillIntervalBounds(Standard_Real               A,
 246                                             Standard_Real               B,
 247                                             const TColStd_Array1OfReal& Knots,
 248                                             const Standard_Integer      NumSubs)
 249 {
 250   Standard_Integer iEnd = Knots.Upper(), jEnd = U1->Upper();
 251   if(iEnd - 1 > jEnd){
 252     iEnd = MaxSubs(iEnd-1,NumSubs);
 253     U1   = new math_Vector(1,iEnd);
 254     U2   = new math_Vector(1,iEnd);
 255     DimU = new math_Vector(1,iEnd);
 256     ErrU = new math_Vector(1,iEnd,0.0);
 257     IxU  = new math_Vector(1,iEnd);
 258     IyU  = new math_Vector(1,iEnd);
 259     IzU  = new math_Vector(1,iEnd);
 260     IxxU = new math_Vector(1,iEnd);
 261     IyyU = new math_Vector(1,iEnd);
 262     IzzU = new math_Vector(1,iEnd);
 263     IxyU = new math_Vector(1,iEnd);
 264     IxzU = new math_Vector(1,iEnd);
 265     IyzU = new math_Vector(1,iEnd);
 266   }
 267   return FillIntervalBounds(A, B, Knots, U1, U2);
 268 }
 269
 270 static Standard_Real CCompute(Face&                  S,
 271                               Domain&                D,
 272                               const gp_Pnt&          loc,
 273                               Standard_Real&         Dim,
 274                               gp_Pnt&                g,
 275                               gp_Mat&                inertia,
 276                               const Standard_Real    EpsDim,
 277                               const Standard_Boolean isErrorCalculation,
 278                               const Standard_Boolean isVerifyComputation)
 279 {
 280   Standard_Real  (*FuncAdd)(Standard_Real, Standard_Real);
 281   Standard_Real  (*FuncMul)(Standard_Real, Standard_Real);
 282
 283   FuncAdd = Addition;
 284   FuncMul = Multiplication;
 285
 286   Standard_Boolean isNaturalRestriction = S.NaturalRestriction();
 287
 288   Standard_Integer NumSubs = SUBS_POWER;
 289
 290   Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
 291   Dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
 292   Standard_Real x, y, z;
 293   //boundary curve parametrization
 294   Standard_Real l1, l2, lm, lr, l;
 295   //Face parametrization in U and V direction
 296   Standard_Real BV1, BV2, v;
 297   Standard_Real BU1, BU2, u1, u2, um, ur, u;
 298   S.Bounds (BU1, BU2, BV1, BV2);
 299   u1 = BU1;
 300
 301   if(Precision::IsInfinite(BU1) || Precision::IsInfinite(BU2) ||
 302      Precision::IsInfinite(BV1) || Precision::IsInfinite(BV2))
 303   {
 304     FuncAdd = AdditionInf;
 305     FuncMul = MultiplicationInf;
 306   }
 307
 308
 309   //location point used to compute the inertia
 310   Standard_Real xloc, yloc, zloc;
 311   loc.Coord (xloc, yloc, zloc); // use member of parent class
 312   //Jacobien (x, y, z) -> (u, v) = ||n||
 313   Standard_Real ds;
 314   //On the Face
 315   gp_Pnt Ps;
 316   gp_Vec VNor;
 317   //On the boundary curve u-v
 318   gp_Pnt2d Puv;
 319   gp_Vec2d Vuv;
 320   Standard_Real Dul;  // Dul = Du / Dl
 321   Standard_Real CDim[2], CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz;
 322   Standard_Real LocDim[2], LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy, LocIxz, LocIyz;
 323
 324   Standard_Real ErrorU, ErrorL, ErrorLMax = 0.0, Eps=0.0, EpsL=0.0, EpsU=0.0;
 325
 326   Standard_Integer iD = 0, NbLSubs, iLS, iLSubEnd, iGL, iGLEnd, NbLGaussP[2], LRange[2], iL, kL, kLEnd, IL, JL;
 327   Standard_Integer i, NbUSubs, iUS, iUSubEnd, iGU, iGUEnd, NbUGaussP[2], URange[2], iU, kU, kUEnd, IU, JU;
 328   Standard_Integer UMaxSubs, LMaxSubs;
 329   iGLEnd = isErrorCalculation? 2: 1;
 330   for(i = 0; i < 2; i++) {
 331     LocDim[i] = 0.0;
 332     CDim[i] = 0.0;
 333   }
 334
 335   NbUGaussP[0] = S.SIntOrder(EpsDim);
 336   NbUGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbUGaussP[0]));
 337   math::GaussPoints(NbUGaussP[0],UGaussP0);  math::GaussWeights(NbUGaussP[0],UGaussW0);
 338   math::GaussPoints(NbUGaussP[1],UGaussP1);  math::GaussWeights(NbUGaussP[1],UGaussW1);
 339
 340   NbUSubs = S.SUIntSubs();
 341   TColStd_Array1OfReal UKnots(1,NbUSubs+1);
 342   S.UKnots(UKnots);
 343
 344   while (isNaturalRestriction || D.More())
 345   {
 346     if(isNaturalRestriction)
 347     {
 348       NbLGaussP[0] = Min(2*NbUGaussP[0],math::GaussPointsMax());
 349     }
 350     else
 351     {
 352       S.Load(D.Value());  ++iD;
 353       NbLGaussP[0] = S.LIntOrder(EpsDim);
 354     }
 355
 356     NbLGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbLGaussP[0]));
 357     math::GaussPoints(NbLGaussP[0],LGaussP0);  math::GaussWeights(NbLGaussP[0],LGaussW0);
 358     math::GaussPoints(NbLGaussP[1],LGaussP1);  math::GaussWeights(NbLGaussP[1],LGaussW1);
 359
 360     NbLSubs = isNaturalRestriction? S.SVIntSubs(): S.LIntSubs();
 361
 362     TColStd_Array1OfReal LKnots(1,NbLSubs+1);
 363     if(isNaturalRestriction)
 364     {
 365       S.VKnots(LKnots);
 366       l1 = BV1; l2 = BV2;
 367     }
 368     else
 369     {
 370       S.LKnots(LKnots);
 371       l1 = S.FirstParameter();  l2 = S.LastParameter();
 372     }
 373
 374     ErrorL = 0.0;
 375     kLEnd = 1; JL = 0;
 376     //OCC503(apo): if(Abs(l2-l1) < EPS_PARAM) continue;
 377     if(Abs(l2-l1) > EPS_PARAM)
 378     {
 379       iLSubEnd = LFillIntervalBounds(l1, l2, LKnots, NumSubs);
 380       LMaxSubs = MaxSubs(iLSubEnd);
 381       if(LMaxSubs > DimL.Vector()->Upper())
 382         LMaxSubs = DimL.Vector()->Upper();
 383
 384       DimL.Init(0.0,1,LMaxSubs);  ErrL.Init(0.0,1,LMaxSubs);  ErrUL.Init(0.0,1,LMaxSubs);
 385
 386       do{// while: L
 387         if(++JL > iLSubEnd)
 388         {
 389           LRange[0] = IL = ErrL->Max();
 390           LRange[1] = JL;
 391           L1(JL) = (L1(IL) + L2(IL))/2.0;
 392           L2(JL) = L2(IL);
 393           L2(IL) = L1(JL);
 394         }
 395         else
 396           LRange[0] = IL = JL;
 397
 398         if(JL == LMaxSubs || Abs(L2(JL) - L1(JL)) < EPS_PARAM)
 399           if(kLEnd == 1)
 400           {
 401             DimL(JL) = ErrL(JL) = IxL(JL) = IyL(JL) = IzL(JL) =
 402               IxxL(JL) = IyyL(JL) = IzzL(JL) = IxyL(JL) = IxzL(JL) = IyzL(JL) = 0.0;
 403           }else{
 404             JL--;
 405             EpsL = ErrorL;  Eps = EpsL/0.9;
 406             break;
 407           }
 408         else
 409           for(kL=0; kL < kLEnd; kL++)
 410           {
 411             iLS = LRange[kL];
 412             lm = 0.5*(L2(iLS) + L1(iLS));
 413             lr = 0.5*(L2(iLS) - L1(iLS));
 414             CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
 415
 416             for(iGL=0; iGL < iGLEnd; iGL++)
 417             {
 418               CDim[iGL] = 0.0;
 419               for(iL=1; iL<=NbLGaussP[iGL]; iL++)
 420               {
 421                 l = lm + lr*(*LGaussP[iGL])(iL);
 422                 if(isNaturalRestriction)
 423                 {
 424                   v = l; u2 = BU2; Dul = (*LGaussW[iGL])(iL);
 425                 }
 426                 else
 427                 {
 428                   S.D12d (l, Puv, Vuv);
 429                   Dul = Vuv.Y()*(*LGaussW[iGL])(iL);  // Dul = Du / Dl
 430                   if(Abs(Dul) < EPS_PARAM)
 431                     continue;
 432
 433                   v  = Puv.Y();
 434                   u2 = Puv.X();
 435                   //Check on cause out off bounds of value current parameter
 436                   if(v < BV1)
 437                     v = BV1;
 438                   else if(v > BV2)
 439                     v = BV2;
 440
 441                   if(u2 < BU1)
 442                     u2 = BU1;
 443                   else if(u2 > BU2)
 444                     u2 = BU2;
 445                 }
 446
 447                 ErrUL(iLS) = 0.0;
 448                 kUEnd = 1;
 449                 JU = 0;
 450
 451                 if(Abs(u2-u1) < EPS_PARAM)
 452                   continue;
 453
 454                 iUSubEnd = UFillIntervalBounds(u1, u2, UKnots, NumSubs);
 455                 UMaxSubs = MaxSubs(iUSubEnd);
 456                 if(UMaxSubs > DimU.Vector()->Upper())
 457                   UMaxSubs = DimU.Vector()->Upper();
 458
 459                 DimU.Init(0.0,1,UMaxSubs);  ErrU.Init(0.0,1,UMaxSubs);  ErrorU = 0.0;
 460
 461                 do{//while: U
 462                   if(++JU > iUSubEnd)
 463                   {
 464                     URange[0] = IU = ErrU->Max();
 465                     URange[1] = JU;
 466                     U1(JU) = (U1(IU)+U2(IU))/2.0;
 467                     U2(JU) = U2(IU);
 468                     U2(IU) = U1(JU);
 469                   }
 470                   else
 471                     URange[0] = IU = JU;
 472
 473                   if(JU == UMaxSubs || Abs(U2(JU) - U1(JU)) < EPS_PARAM)
 474                     if(kUEnd == 1)
 475                     {
 476                       DimU(JU) = ErrU(JU) = IxU(JU) = IyU(JU) = IzU(JU) =
 477                         IxxU(JU) = IyyU(JU) = IzzU(JU) = IxyU(JU) = IxzU(JU) = IyzU(JU) = 0.0;
 478                     }else
 479                     {
 480                       JU--;
 481                       EpsU = ErrorU;  Eps = EpsU*Abs((u2-u1)*Dul)/0.1;  EpsL = 0.9*Eps;
 482                       break;
 483                     }
 484                   else
 485                     for(kU=0; kU < kUEnd; kU++)
 486                     {
 487                       iUS = URange[kU];
 488                       um = 0.5*(U2(iUS) + U1(iUS));
 489                       ur = 0.5*(U2(iUS) - U1(iUS));
 490                       LocIx = LocIy = LocIz = LocIxx = LocIyy = LocIzz = LocIxy = LocIxz = LocIyz = 0.0;
 491                       iGUEnd = iGLEnd - iGL;
 492                       for(iGU=0; iGU < iGUEnd; iGU++)
 493                       {
 494                         LocDim[iGU] = 0.0;
 495                         for(iU=1; iU <= NbUGaussP[iGU]; iU++)
 496                         {
 497                           u = um + ur*(*UGaussP[iGU])(iU);
 498                           S.Normal(u, v, Ps, VNor);
 499                           ds = VNor.Magnitude();    //Jacobien(x,y,z) -> (u,v)=||n||
 500                           ds *= (*UGaussW[iGU])(iU);
 501                           LocDim[iGU] += ds;
 502
 503                           if(iGU > 0)
 504                             continue;
 505
 506                           Ps.Coord(x, y, z);
 507                           x = FuncAdd(x, -xloc);
 508                           y = FuncAdd(y, -yloc);
 509                           z = FuncAdd(z, -zloc);
 510
 511                           const Standard_Real XdS = FuncMul(x, ds);
 512                           const Standard_Real YdS = FuncMul(y, ds);
 513                           const Standard_Real ZdS = FuncMul(z, ds);
 514
 515                           LocIx = FuncAdd(LocIx, XdS);
 516                           LocIy = FuncAdd(LocIy, YdS);
 517                           LocIz = FuncAdd(LocIz, ZdS);
 518                           LocIxy = FuncAdd(LocIxy, FuncMul(x, YdS));
 519                           LocIyz = FuncAdd(LocIyz, FuncMul(y, ZdS));
 520                           LocIxz = FuncAdd(LocIxz, FuncMul(x, ZdS));
 521
 522                           const Standard_Real XXdS = FuncMul(x, XdS);
 523                           const Standard_Real YYdS = FuncMul(y, YdS);
 524                           const Standard_Real ZZdS = FuncMul(z, ZdS);
 525
 526                           LocIxx = FuncAdd(LocIxx, FuncAdd(YYdS, ZZdS));
 527                           LocIyy = FuncAdd(LocIyy, FuncAdd(XXdS, ZZdS));
 528                           LocIzz = FuncAdd(LocIzz, FuncAdd(XXdS, YYdS));
 529                         }//for: iU
 530                       }//for: iGU
 531
 532                       DimU(iUS) = FuncMul(LocDim[0],ur);
 533                       if(iGL > 0)
 534                         continue;
 535
 536                       ErrU(iUS) = FuncMul(Abs(LocDim[1]-LocDim[0]), ur);
 537                       IxU(iUS) = FuncMul(LocIx, ur);
 538                       IyU(iUS) = FuncMul(LocIy, ur);
 539                       IzU(iUS) = FuncMul(LocIz, ur);
 540                       IxxU(iUS) = FuncMul(LocIxx, ur);
 541                       IyyU(iUS) = FuncMul(LocIyy, ur);
 542                       IzzU(iUS) = FuncMul(LocIzz, ur);
 543                       IxyU(iUS) = FuncMul(LocIxy, ur);
 544                       IxzU(iUS) = FuncMul(LocIxz, ur);
 545                       IyzU(iUS) = FuncMul(LocIyz, ur);
 546                     }//for: kU (iUS)
 547
 548                     if(JU == iUSubEnd)
 549                       kUEnd = 2;
 550
 551                     if(kUEnd == 2)
 552                       ErrorU = ErrU(ErrU->Max());
 553                 }while((ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1);
 554
 555                 for(i=1; i<=JU; i++)
 556                   CDim[iGL] = FuncAdd(CDim[iGL], FuncMul(DimU(i), Dul));
 557
 558                 if(iGL > 0)
 559                   continue;
 560
 561                 ErrUL(iLS) = ErrorU*Abs((u2-u1)*Dul);
 562                 for(i=1; i<=JU; i++)
 563                 {
 564                   CIx = FuncAdd(CIx, FuncMul(IxU(i), Dul));
 565                   CIy = FuncAdd(CIy, FuncMul(IyU(i), Dul));
 566                   CIz = FuncAdd(CIz, FuncMul(IzU(i), Dul));
 567                   CIxx = FuncAdd(CIxx, FuncMul(IxxU(i), Dul));
 568                   CIyy = FuncAdd(CIyy, FuncMul(IyyU(i), Dul));
 569                   CIzz = FuncAdd(CIzz, FuncMul(IzzU(i), Dul));
 570                   CIxy = FuncAdd(CIxy, FuncMul(IxyU(i), Dul));
 571                   CIxz = FuncAdd(CIxz, FuncMul(IxzU(i), Dul));
 572                   CIyz = FuncAdd(CIyz, FuncMul(IyzU(i), Dul));
 573                 }
 574               }//for: iL
 575             }//for: iGL
 576
 577             DimL(iLS) = FuncMul(CDim[0], lr);
 578             if(iGLEnd == 2)
 579               ErrL(iLS) = FuncAdd(FuncMul(Abs(CDim[1]-CDim[0]),lr), ErrUL(iLS));
 580
 581             IxL(iLS) = FuncMul(CIx, lr);
 582             IyL(iLS) = FuncMul(CIy, lr);
 583             IzL(iLS) = FuncMul(CIz, lr);
 584             IxxL(iLS) = FuncMul(CIxx, lr);
 585             IyyL(iLS) = FuncMul(CIyy, lr);
 586             IzzL(iLS) = FuncMul(CIzz, lr);
 587             IxyL(iLS) = FuncMul(CIxy, lr);
 588             IxzL(iLS) = FuncMul(CIxz, lr);
 589             IyzL(iLS) = FuncMul(CIyz, lr);
 590           }//for: (kL)iLS
 591           //  Calculate/correct epsilon of computation by current value of Dim
 592           //That is need for not spend time for
 593           if(JL == iLSubEnd)
 594           {
 595             kLEnd = 2;
 596             Standard_Real DDim = 0.0;
 597             for(i=1; i<=JL; i++)
 598               DDim += DimL(i);
 599
 600             DDim = Abs(DDim*EpsDim);
 601             if(DDim > Eps)
 602             {
 603               Eps = DDim;
 604               EpsL = 0.9*Eps;
 605             }
 606           }
 607
 608           if(kLEnd == 2)
 609             ErrorL = ErrL(ErrL->Max());
 610         }while((ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1);
 611
 612         for(i=1; i<=JL; i++)
 613         {
 614           Dim = FuncAdd(Dim, DimL(i));
 615           Ix = FuncAdd(Ix, IxL(i));
 616           Iy = FuncAdd(Iy, IyL(i));
 617           Iz = FuncAdd(Iz, IzL(i));
 618           Ixx = FuncAdd(Ixx, IxxL(i));
 619           Iyy = FuncAdd(Iyy, IyyL(i));
 620           Izz = FuncAdd(Izz, IzzL(i));
 621           Ixy = FuncAdd(Ixy, IxyL(i));
 622           Ixz = FuncAdd(Ixz, IxzL(i));
 623           Iyz = FuncAdd(Iyz, IyzL(i));
 624         }
 625
 626         ErrorLMax = Max(ErrorLMax, ErrorL);
 627       }
 628
 629       if(isNaturalRestriction)
 630         break;
 631
 632       D.Next();
 633     }
 634
 635     if(Abs(Dim) >= EPS_DIM)
 636     {
 637       Ix /= Dim;
 638       Iy /= Dim;
 639       Iz /= Dim;
 640       g.SetCoord (Ix, Iy, Iz);
 641     }
 642     else
 643     {
 644       Dim =0.0;
 645       g.SetCoord (0., 0.,0.);
 646     }
 647
 648     inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
 649                       gp_XYZ (-Ixy, Iyy, -Iyz),
 650                       gp_XYZ (-Ixz, -Iyz, Izz));
 651
 652     if(iGLEnd == 2)
 653       Eps = Dim != 0.0? ErrorLMax/Abs(Dim): 0.0;
 654     else
 655       Eps = EpsDim;
 656
 657     return Eps;
 658 }
 659
 660 static Standard_Real Compute(Face& S, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia,
 661                              Standard_Real EpsDim)
 662 {
 663   Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
 664   Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
 665   EpsDim = Abs(EpsDim);
 666   Domain D;
 667   return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
 668 }
 669
 670 static Standard_Real Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia,
 671                              Standard_Real EpsDim)
 672 {
 673   Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
 674   Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
 675   EpsDim = Abs(EpsDim);
 676   return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
 677 }
 678
 679 static void Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& dim, gp_Pnt& g, gp_Mat& inertia)
 680 {
 681   Standard_Real  (*FuncAdd)(Standard_Real, Standard_Real);
 682   Standard_Real  (*FuncMul)(Standard_Real, Standard_Real);
 683
 684   FuncAdd = Addition;
 685   FuncMul = Multiplication;
 686
 687   Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
 688   dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
 689
 690   Standard_Real x, y, z;
 691   Standard_Integer NbCGaussgp_Pnts = 0;
 692
 693   Standard_Real l1, l2, lm, lr, l;   //boundary curve parametrization
 694   Standard_Real v1, v2,         v;   //Face parametrization in v direction
 695   Standard_Real u1, u2, um, ur, u;
 696   Standard_Real ds;                  //Jacobien (x, y, z) -> (u, v) = ||n||
 697
 698   gp_Pnt P;                    //On the Face
 699   gp_Vec VNor;
 700
 701   gp_Pnt2d Puv;                //On the boundary curve u-v
 702   gp_Vec2d Vuv;
 703   Standard_Real Dul;                 // Dul = Du / Dl
 704   Standard_Real CArea, CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz;
 705   Standard_Real LocArea, LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy,
 706     LocIxz, LocIyz;
 707
 708
 709   S.Bounds (u1, u2, v1, v2);
 710
 711   if(Precision::IsInfinite(u1) || Precision::IsInfinite(u2) ||
 712      Precision::IsInfinite(v1) || Precision::IsInfinite(v2))
 713   {
 714     FuncAdd = AdditionInf;
 715     FuncMul = MultiplicationInf;
 716   }
 717
 718
 719   Standard_Integer NbUGaussgp_Pnts = Min(S.UIntegrationOrder (),
 720     math::GaussPointsMax());
 721   Standard_Integer NbVGaussgp_Pnts = Min(S.VIntegrationOrder (),
 722     math::GaussPointsMax());
 723
 724   Standard_Integer NbGaussgp_Pnts = Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts);
 725
 726   //Number of Gauss points for the integration
 727   //on the Face
 728   math_Vector GaussSPV (1, NbGaussgp_Pnts);
 729   math_Vector GaussSWV (1, NbGaussgp_Pnts);
 730   math::GaussPoints  (NbGaussgp_Pnts,GaussSPV);
 731   math::GaussWeights (NbGaussgp_Pnts,GaussSWV);
 732
 733
 734   //location point used to compute the inertia
 735   Standard_Real xloc, yloc, zloc;
 736   loc.Coord (xloc, yloc, zloc);
 737
 738   while (D.More()) {
 739
 740     S.Load(D.Value());
 741     NbCGaussgp_Pnts =  Min(S.IntegrationOrder (), math::GaussPointsMax());
 742
 743     math_Vector GaussCP (1, NbCGaussgp_Pnts);
 744     math_Vector GaussCW (1, NbCGaussgp_Pnts);
 745     math::GaussPoints  (NbCGaussgp_Pnts,GaussCP);
 746     math::GaussWeights (NbCGaussgp_Pnts,GaussCW);
 747
 748     CArea = 0.0;
 749     CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
 750     l1 = S.FirstParameter ();
 751     l2 = S.LastParameter  ();
 752     lm = 0.5 * (l2 + l1);
 753     lr = 0.5 * (l2 - l1);
 754
 755     for (Standard_Integer i = 1; i <= NbCGaussgp_Pnts; i++) {
 756       l = lm + lr * GaussCP (i);
 757       S.D12d(l, Puv, Vuv);
 758       v   = Puv.Y();
 759       u2  = Puv.X();
 760       Dul = Vuv.Y();
 761       Dul *= GaussCW (i);
 762       um  = 0.5 * (u2 + u1);
 763       ur  = 0.5 * (u2 - u1);
 764       LocArea = LocIx  = LocIy  = LocIz = LocIxx = LocIyy = LocIzz =
 765         LocIxy  = LocIxz = LocIyz = 0.0;
 766       for (Standard_Integer j = 1; j <= NbGaussgp_Pnts; j++) {
 767         u = FuncAdd(um, FuncMul(ur, GaussSPV (j)));
 768         S.Normal (u, v, P, VNor);
 769         ds = VNor.Magnitude();    //normal.Magnitude
 770         ds = FuncMul(ds, Dul) * GaussSWV (j);
 771         LocArea = FuncAdd(LocArea, ds);
 772         P.Coord (x, y, z);
 773
 774         x = FuncAdd(x, -xloc);
 775         y = FuncAdd(y, -yloc);
 776         z = FuncAdd(z, -zloc);
 777
 778         const Standard_Real XdS = FuncMul(x, ds);
 779         const Standard_Real YdS = FuncMul(y, ds);
 780         const Standard_Real ZdS = FuncMul(z, ds);
 781
 782         LocIx = FuncAdd(LocIx, XdS);
 783         LocIy = FuncAdd(LocIy, YdS);
 784         LocIz = FuncAdd(LocIz, ZdS);
 785         LocIxy = FuncAdd(LocIxy, FuncMul(x, YdS));
 786         LocIyz = FuncAdd(LocIyz, FuncMul(y, ZdS));
 787         LocIxz = FuncAdd(LocIxz, FuncMul(x, ZdS));
 788
 789         const Standard_Real XXdS = FuncMul(x, XdS);
 790         const Standard_Real YYdS = FuncMul(y, YdS);
 791         const Standard_Real ZZdS = FuncMul(z, ZdS);
 792
 793         LocIxx = FuncAdd(LocIxx, FuncAdd(YYdS, ZZdS));
 794         LocIyy = FuncAdd(LocIyy, FuncAdd(XXdS, ZZdS));
 795         LocIzz = FuncAdd(LocIzz, FuncAdd(XXdS, YYdS));
 796       }
 797
 798       CArea = FuncAdd(CArea, FuncMul(LocArea, ur));
 799       CIx   = FuncAdd(CIx, FuncMul(LocIx, ur));
 800       CIy   = FuncAdd(CIy, FuncMul(LocIy, ur));
 801       CIz   = FuncAdd(CIz, FuncMul(LocIz, ur));
 802       CIxx  = FuncAdd(CIxx, FuncMul(LocIxx, ur));
 803       CIyy  = FuncAdd(CIyy, FuncMul(LocIyy, ur));
 804       CIzz  = FuncAdd(CIzz, FuncMul(LocIzz, ur));
 805       CIxy  = FuncAdd(CIxy, FuncMul(LocIxy, ur));
 806       CIxz  = FuncAdd(CIxz, FuncMul(LocIxz, ur));
 807       CIyz  = FuncAdd(CIyz, FuncMul(LocIyz, ur));
 808     }
 809
 810     dim = FuncAdd(dim, FuncMul(CArea, lr));
 811     Ix  = FuncAdd(Ix,  FuncMul(CIx, lr));
 812     Iy  = FuncAdd(Iy,  FuncMul(CIy, lr));
 813     Iz  = FuncAdd(Iz,  FuncMul(CIz, lr));
 814     Ixx = FuncAdd(Ixx, FuncMul(CIxx, lr));
 815     Iyy = FuncAdd(Iyy, FuncMul(CIyy, lr));
 816     Izz = FuncAdd(Izz, FuncMul(CIzz, lr));
 817     Ixy = FuncAdd(Ixy, FuncMul(CIxy, lr));
 818     Ixz = FuncAdd(Iyz, FuncMul(CIxz, lr));
 819     Iyz = FuncAdd(Ixz, FuncMul(CIyz, lr));
 820     D.Next();
 821   }
 822
 823   if (Abs(dim) >= EPS_DIM) {
 824     Ix /= dim;
 825     Iy /= dim;
 826     Iz /= dim;
 827     g.SetCoord (Ix, Iy, Iz);
 828   }
 829   else {
 830     dim =0.;
 831     g.SetCoord (0., 0.,0.);
 832   }
 833
 834   inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
 835                     gp_XYZ (-Ixy, Iyy, -Iyz),
 836                     gp_XYZ (-Ixz, -Iyz, Izz));
 837 }
 838
 839 static void Compute(const Face& S,
 840                     const gp_Pnt& loc,
 841                     Standard_Real& dim,
 842                     gp_Pnt& g,
 843                     gp_Mat& inertia)
 844 {
 845   Standard_Real  (*FuncAdd)(Standard_Real, Standard_Real);
 846   Standard_Real  (*FuncMul)(Standard_Real, Standard_Real);
 847
 848   FuncAdd = Addition;
 849   FuncMul = Multiplication;
 850
 851   Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
 852   dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
 853
 854   Standard_Real LowerU, UpperU, LowerV, UpperV;
 855   S.Bounds (LowerU, UpperU, LowerV, UpperV);
 856
 857   if(Precision::IsInfinite(LowerU) || Precision::IsInfinite(UpperU) ||
 858      Precision::IsInfinite(LowerV) || Precision::IsInfinite(UpperV))
 859   {
 860     FuncAdd = AdditionInf;
 861     FuncMul = MultiplicationInf;
 862   }
 863
 864   Standard_Integer UOrder = Min(S.UIntegrationOrder (),
 865     math::GaussPointsMax());
 866   Standard_Integer VOrder = Min(S.VIntegrationOrder (),
 867     math::GaussPointsMax());
 868   gp_Pnt P;
 869   gp_Vec VNor;
 870   Standard_Real dsi, ds;
 871   Standard_Real ur, um, u, vr, vm, v;
 872   Standard_Real x, y, z;
 873   Standard_Real Ixi, Iyi, Izi, Ixxi, Iyyi, Izzi, Ixyi, Ixzi, Iyzi;
 874   Standard_Real xloc, yloc, zloc;
 875   loc.Coord (xloc, yloc, zloc);
 876
 877   Standard_Integer i, j;
 878   math_Vector GaussPU (1, UOrder);     //gauss points and weights
 879   math_Vector GaussWU (1, UOrder);
 880   math_Vector GaussPV (1, VOrder);
 881   math_Vector GaussWV (1, VOrder);
 882
 883   //Recuperation des points de Gauss dans le fichier GaussPoints.
 884   math::GaussPoints  (UOrder,GaussPU);
 885   math::GaussWeights (UOrder,GaussWU);
 886   math::GaussPoints  (VOrder,GaussPV);
 887   math::GaussWeights (VOrder,GaussWV);
 888
 889   // Calcul des integrales aux points de gauss :
 890   um = 0.5 * FuncAdd(UpperU,  LowerU);
 891   vm = 0.5 * FuncAdd(UpperV,  LowerV);
 892   ur = 0.5 * FuncAdd(UpperU, -LowerU);
 893   vr = 0.5 * FuncAdd(UpperV, -LowerV);
 894
 895   for (j = 1; j <= VOrder; j++) {
 896     v = FuncAdd(vm, FuncMul(vr, GaussPV(j)));
 897     dsi = Ixi = Iyi = Izi = Ixxi = Iyyi = Izzi = Ixyi = Ixzi = Iyzi = 0.0;
 898
 899     for (i = 1; i <= UOrder; i++) {
 900       u = FuncAdd(um, FuncMul(ur, GaussPU (i)));
 901       S.Normal (u, v, P, VNor);
 902       ds = FuncMul(VNor.Magnitude(), GaussWU (i));
 903       P.Coord (x, y, z);
 904
 905       x = FuncAdd(x, -xloc);
 906       y = FuncAdd(y, -yloc);
 907       z = FuncAdd(z, -zloc);
 908
 909       dsi = FuncAdd(dsi, ds);
 910
 911       const Standard_Real XdS = FuncMul(x, ds);
 912       const Standard_Real YdS = FuncMul(y, ds);
 913       const Standard_Real ZdS = FuncMul(z, ds);
 914
 915       Ixi = FuncAdd(Ixi, XdS);
 916       Iyi = FuncAdd(Iyi, YdS);
 917       Izi = FuncAdd(Izi, ZdS);
 918       Ixyi = FuncAdd(Ixyi, FuncMul(x, YdS));
 919       Iyzi = FuncAdd(Iyzi, FuncMul(y, ZdS));
 920       Ixzi = FuncAdd(Ixzi, FuncMul(x, ZdS));
 921
 922       const Standard_Real XXdS = FuncMul(x, XdS);
 923       const Standard_Real YYdS = FuncMul(y, YdS);
 924       const Standard_Real ZZdS = FuncMul(z, ZdS);
 925
 926       Ixxi = FuncAdd(Ixxi, FuncAdd(YYdS, ZZdS));
 927       Iyyi = FuncAdd(Iyyi, FuncAdd(XXdS, ZZdS));
 928       Izzi = FuncAdd(Izzi, FuncAdd(XXdS, YYdS));
 929     }
 930
 931     dim   = FuncAdd(dim, FuncMul(dsi, GaussWV (j)));
 932     Ix    = FuncAdd(Ix,  FuncMul(Ixi, GaussWV (j)));
 933     Iy    = FuncAdd(Iy,  FuncMul(Iyi, GaussWV (j)));
 934     Iz    = FuncAdd(Iz,  FuncMul(Izi, GaussWV (j)));
 935     Ixx   = FuncAdd(Ixx, FuncMul(Ixxi, GaussWV (j)));
 936     Iyy   = FuncAdd(Iyy, FuncMul(Iyyi, GaussWV (j)));
 937     Izz   = FuncAdd(Izz, FuncMul(Izzi, GaussWV (j)));
 938     Ixy   = FuncAdd(Ixy, FuncMul(Ixyi, GaussWV (j)));
 939     Iyz   = FuncAdd(Iyz, FuncMul(Iyzi, GaussWV (j)));
 940     Ixz   = FuncAdd(Ixz, FuncMul(Ixzi, GaussWV (j)));
 941   }
 942
 943   vr  = FuncMul(vr, ur);
 944   Ixx = FuncMul(vr, Ixx);
 945   Iyy = FuncMul(vr, Iyy);
 946   Izz = FuncMul(vr, Izz);
 947   Ixy = FuncMul(vr, Ixy);
 948   Ixz = FuncMul(vr, Ixz);
 949   Iyz = FuncMul(vr, Iyz);
 950
 951   if (Abs(dim) >= EPS_DIM)
 952   {
 953     Ix    /= dim;
 954     Iy    /= dim;
 955     Iz    /= dim;
 956     dim   *= vr;
 957     g.SetCoord (Ix, Iy, Iz);
 958   }
 959   else
 960   {
 961     dim =0.;
 962     g.SetCoord (0.,0.,0.);
 963   }
 964
 965   inertia = gp_Mat (gp_XYZ ( Ixx, -Ixy, -Ixz),
 966                     gp_XYZ (-Ixy,  Iyy, -Iyz),
 967                     gp_XYZ (-Ixz, -Iyz,  Izz));
 968 }
 969
 970 GProp_SGProps::GProp_SGProps(){}
 971
 972 GProp_SGProps::GProp_SGProps (const Face&   S,
 973                               const gp_Pnt& SLocation
 974                               )
 975 {
 976   SetLocation(SLocation);
 977   Perform(S);
 978 }
 979
 980 GProp_SGProps::GProp_SGProps (Face&   S,
 981                               Domain& D,
 982                               const gp_Pnt& SLocation
 983                               )
 984 {
 985   SetLocation(SLocation);
 986   Perform(S,D);
 987 }
 988
 989 GProp_SGProps::GProp_SGProps(Face& S, const gp_Pnt& SLocation, const Standard_Real Eps){
 990   SetLocation(SLocation);
 991   Perform(S, Eps);
 992 }
 993
 994 GProp_SGProps::GProp_SGProps(Face& S, Domain& D, const gp_Pnt& SLocation, const Standard_Real Eps){
 995   SetLocation(SLocation);
 996   Perform(S, D, Eps);
 997 }
 998
 999 void GProp_SGProps::SetLocation(const gp_Pnt& SLocation){
1000   loc = SLocation;
1001 }
1002
1003 void GProp_SGProps::Perform(const Face& S){
1004   Compute(S,loc,dim,g,inertia);
1005   myEpsilon = 1.0;
1006   return;
1007 }
1008
1009 void GProp_SGProps::Perform(Face& S, Domain& D){
1010   Compute(S,D,loc,dim,g,inertia);
1011   myEpsilon = 1.0;
1012   return;
1013 }
1014
1015 Standard_Real GProp_SGProps::Perform(Face& S, const Standard_Real Eps){
1016   return myEpsilon = Compute(S,loc,dim,g,inertia,Eps);
1017 }
1018
1019 Standard_Real GProp_SGProps::Perform(Face& S, Domain& D, const Standard_Real Eps){
1020   return myEpsilon = Compute(S,D,loc,dim,g,inertia,Eps);
1021 }
1022
1023
1024 Standard_Real GProp_SGProps::GetEpsilon(){
1025   return myEpsilon;
1026 }