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