src/GProp/GProp_VGProps.gxx

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