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