src/BRepGProp/BRepGProp_Vinert.cxx

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