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