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