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