0024428: Implementation of LGPL license
[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
7 // under the terms of the GNU Lesser General Public 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                           x = Precision::IsInfinite(x) ? Precision::Infinite() : x*x;
522                           y = Precision::IsInfinite(y) ? Precision::Infinite() : y*y;
523                           z = Precision::IsInfinite(z) ? Precision::Infinite() : z*z;
524                           LocIxx = FuncAdd(LocIxx, FuncAdd(YdS, ZdS));
525                           LocIyy = FuncAdd(LocIyy, FuncAdd(XdS, ZdS));
526                           LocIzz = FuncAdd(LocIzz, FuncAdd(XdS, YdS));
527                         }//for: iU
528                       }//for: iGU
529
530                       DimU(iUS) = FuncMul(LocDim[0],ur);
531                       if(iGL > 0)
532                         continue;
533
534                       ErrU(iUS) = FuncMul(Abs(LocDim[1]-LocDim[0]), ur);
535                       IxU(iUS) = FuncMul(LocIx, ur);
536                       IyU(iUS) = FuncMul(LocIy, ur);
537                       IzU(iUS) = FuncMul(LocIz, ur);
538                       IxxU(iUS) = FuncMul(LocIxx, ur);
539                       IyyU(iUS) = FuncMul(LocIyy, ur);
540                       IzzU(iUS) = FuncMul(LocIzz, ur);
541                       IxyU(iUS) = FuncMul(LocIxy, ur);
542                       IxzU(iUS) = FuncMul(LocIxz, ur);
543                       IyzU(iUS) = FuncMul(LocIyz, ur);
544                     }//for: kU (iUS)
545
546                     if(JU == iUSubEnd)
547                       kUEnd = 2;
548
549                     if(kUEnd == 2)
550                       ErrorU = ErrU(ErrU->Max());
551                 }while((ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1);
552
553                 for(i=1; i<=JU; i++)
554                   CDim[iGL] = FuncAdd(CDim[iGL], FuncMul(DimU(i), Dul));
555
556                 if(iGL > 0)
557                   continue;
558
559                 ErrUL(iLS) = ErrorU*Abs((u2-u1)*Dul);
560                 for(i=1; i<=JU; i++)
561                 {
562                   CIx = FuncAdd(CIx, FuncMul(IxU(i), Dul));
563                   CIy = FuncAdd(CIy, FuncMul(IyU(i), Dul));
564                   CIz = FuncAdd(CIz, FuncMul(IzU(i), Dul));
565                   CIxx = FuncAdd(CIxx, FuncMul(IxxU(i), Dul));
566                   CIyy = FuncAdd(CIyy, FuncMul(IyyU(i), Dul));
567                   CIzz = FuncAdd(CIzz, FuncMul(IzzU(i), Dul));
568                   CIxy = FuncAdd(CIxy, FuncMul(IxyU(i), Dul));
569                   CIxz = FuncAdd(CIxz, FuncMul(IxzU(i), Dul));
570                   CIyz = FuncAdd(CIyz, FuncMul(IyzU(i), Dul));
571                 }
572               }//for: iL 
573             }//for: iGL
574
575             DimL(iLS) = FuncMul(CDim[0], lr);  
576             if(iGLEnd == 2)
577               ErrL(iLS) = FuncAdd(FuncMul(Abs(CDim[1]-CDim[0]),lr), ErrUL(iLS));
578
579             IxL(iLS) = FuncMul(CIx, lr);
580             IyL(iLS) = FuncMul(CIy, lr);
581             IzL(iLS) = FuncMul(CIz, lr); 
582             IxxL(iLS) = FuncMul(CIxx, lr);
583             IyyL(iLS) = FuncMul(CIyy, lr);
584             IzzL(iLS) = FuncMul(CIzz, lr); 
585             IxyL(iLS) = FuncMul(CIxy, lr);
586             IxzL(iLS) = FuncMul(CIxz, lr);
587             IyzL(iLS) = FuncMul(CIyz, lr);
588           }//for: (kL)iLS
589           //  Calculate/correct epsilon of computation by current value of Dim
590           //That is need for not spend time for 
591           if(JL == iLSubEnd)
592           {  
593             kLEnd = 2; 
594             Standard_Real DDim = 0.0;
595             for(i=1; i<=JL; i++)
596               DDim += DimL(i);
597
598             DDim = Abs(DDim*EpsDim);
599             if(DDim > Eps)
600             { 
601               Eps = DDim;
602               EpsL = 0.9*Eps;
603             }
604           }
605
606           if(kLEnd == 2)
607             ErrorL = ErrL(ErrL->Max());
608         }while((ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1);
609
610         for(i=1; i<=JL; i++)
611         {
612           Dim = FuncAdd(Dim, DimL(i)); 
613           Ix = FuncAdd(Ix, IxL(i));
614           Iy = FuncAdd(Iy, IyL(i));
615           Iz = FuncAdd(Iz, IzL(i));
616           Ixx = FuncAdd(Ixx, IxxL(i));
617           Iyy = FuncAdd(Iyy, IyyL(i));
618           Izz = FuncAdd(Izz, IzzL(i));
619           Ixy = FuncAdd(Ixy, IxyL(i));
620           Ixz = FuncAdd(Ixz, IxzL(i));
621           Iyz = FuncAdd(Iyz, IyzL(i));
622         }
623
624         ErrorLMax = Max(ErrorLMax, ErrorL);
625       }
626
627       if(isNaturalRestriction)
628         break;
629
630       D.Next();
631     }
632
633     if(Abs(Dim) >= EPS_DIM)
634     {
635       Ix /= Dim;
636       Iy /= Dim;
637       Iz /= Dim;
638       g.SetCoord (Ix, Iy, Iz);
639     }
640     else
641     {
642       Dim =0.0;
643       g.SetCoord (0., 0.,0.);
644     }
645
646     inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
647                       gp_XYZ (-Ixy, Iyy, -Iyz),
648                       gp_XYZ (-Ixz, -Iyz, Izz));
649
650     if(iGLEnd == 2)
651       Eps = Dim != 0.0? ErrorLMax/Abs(Dim): 0.0;
652     else
653       Eps = EpsDim;
654
655     return Eps;
656 }
657
658 static Standard_Real Compute(Face& S, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, 
659                              Standard_Real EpsDim) 
660 {
661   Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
662   Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
663   EpsDim = Abs(EpsDim);
664   Domain D;
665   return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
666 }
667
668 static Standard_Real Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, 
669                              Standard_Real EpsDim) 
670 {
671   Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
672   Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
673   EpsDim = Abs(EpsDim);
674   return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
675 }
676
677 static void Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& dim, gp_Pnt& g, gp_Mat& inertia)
678 {
679   Standard_Real  (*FuncAdd)(Standard_Real, Standard_Real);
680   Standard_Real  (*FuncMul)(Standard_Real, Standard_Real);
681
682   FuncAdd = Addition;
683   FuncMul = Multiplication;
684
685   Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
686   dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
687
688   Standard_Real x, y, z;
689   Standard_Integer NbCGaussgp_Pnts = 0;
690
691   Standard_Real l1, l2, lm, lr, l;   //boundary curve parametrization
692   Standard_Real v1, v2, vm, vr, v;   //Face parametrization in v direction
693   Standard_Real u1, u2, um, ur, u;
694   Standard_Real ds;                  //Jacobien (x, y, z) -> (u, v) = ||n||
695
696   gp_Pnt P;                    //On the Face
697   gp_Vec VNor;
698
699   gp_Pnt2d Puv;                //On the boundary curve u-v
700   gp_Vec2d Vuv;
701   Standard_Real Dul;                 // Dul = Du / Dl
702   Standard_Real CArea, CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz;
703   Standard_Real LocArea, LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy,
704     LocIxz, LocIyz;
705
706
707   S.Bounds (u1, u2, v1, v2);
708
709   if(Precision::IsInfinite(u1) || Precision::IsInfinite(u2) ||
710      Precision::IsInfinite(v1) || Precision::IsInfinite(v2))
711   {
712     FuncAdd = AdditionInf;
713     FuncMul = MultiplicationInf;
714   }
715
716
717   Standard_Integer NbUGaussgp_Pnts = Min(S.UIntegrationOrder (),
718     math::GaussPointsMax());
719   Standard_Integer NbVGaussgp_Pnts = Min(S.VIntegrationOrder (),
720     math::GaussPointsMax());
721
722   Standard_Integer NbGaussgp_Pnts = Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts);
723
724   //Number of Gauss points for the integration
725   //on the Face
726   math_Vector GaussSPV (1, NbGaussgp_Pnts);
727   math_Vector GaussSWV (1, NbGaussgp_Pnts);
728   math::GaussPoints  (NbGaussgp_Pnts,GaussSPV);
729   math::GaussWeights (NbGaussgp_Pnts,GaussSWV);
730
731
732   //location point used to compute the inertia
733   Standard_Real xloc, yloc, zloc;
734   loc.Coord (xloc, yloc, zloc);
735
736   while (D.More()) {
737
738     S.Load(D.Value());
739     NbCGaussgp_Pnts =  Min(S.IntegrationOrder (), math::GaussPointsMax());        
740
741     math_Vector GaussCP (1, NbCGaussgp_Pnts);
742     math_Vector GaussCW (1, NbCGaussgp_Pnts);
743     math::GaussPoints  (NbCGaussgp_Pnts,GaussCP);
744     math::GaussWeights (NbCGaussgp_Pnts,GaussCW);
745
746     CArea = 0.0;
747     CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
748     l1 = S.FirstParameter ();
749     l2 = S.LastParameter  ();
750     lm = 0.5 * (l2 + l1);         
751     lr = 0.5 * (l2 - l1);
752
753     Puv = S.Value2d (lm);
754     vm  = Puv.Y();
755     Puv = S.Value2d (lr);
756     vr  = Puv.Y();
757
758     for (Standard_Integer i = 1; i <= NbCGaussgp_Pnts; i++) {
759       l = lm + lr * GaussCP (i);
760       S.D12d(l, Puv, Vuv);
761       v   = Puv.Y();
762       u2  = Puv.X();
763       Dul = Vuv.Y();
764       Dul *= GaussCW (i);
765       um  = 0.5 * (u2 + u1);
766       ur  = 0.5 * (u2 - u1);
767       LocArea = LocIx  = LocIy  = LocIz = LocIxx = LocIyy = LocIzz = 
768         LocIxy  = LocIxz = LocIyz = 0.0;
769       for (Standard_Integer j = 1; j <= NbGaussgp_Pnts; j++) {
770         u = FuncAdd(um, FuncMul(ur, GaussSPV (j)));
771         S.Normal (u, v, P, VNor);
772         ds = VNor.Magnitude();    //normal.Magnitude
773         ds = FuncMul(ds, Dul) * GaussSWV (j); 
774         LocArea = FuncAdd(LocArea, ds); 
775         P.Coord (x, y, z);
776
777         x = FuncAdd(x, -xloc);
778         y = FuncAdd(y, -yloc);
779         z = FuncAdd(z, -zloc);
780
781         const Standard_Real XdS = FuncMul(x, ds);
782         const Standard_Real YdS = FuncMul(y, ds);
783         const Standard_Real ZdS = FuncMul(z, ds);
784         
785         LocIx = FuncAdd(LocIx, XdS);
786         LocIy = FuncAdd(LocIy, YdS);
787         LocIz = FuncAdd(LocIz, ZdS);
788         LocIxy = FuncAdd(LocIxy, FuncMul(x, YdS));
789         LocIyz = FuncAdd(LocIyz, FuncMul(y, ZdS));
790         LocIxz = FuncAdd(LocIxz, FuncMul(x, ZdS));
791         x = Precision::IsInfinite(x) ? Precision::Infinite() : x*x;
792         y = Precision::IsInfinite(y) ? Precision::Infinite() : y*y;
793         z = Precision::IsInfinite(z) ? Precision::Infinite() : z*z;
794         LocIxx = FuncAdd(LocIxx, FuncAdd(YdS, ZdS));
795         LocIyy = FuncAdd(LocIyy, FuncAdd(XdS, ZdS));
796         LocIzz = FuncAdd(LocIzz, FuncAdd(XdS, YdS));
797       }
798
799       CArea = FuncAdd(CArea, FuncMul(LocArea, ur));
800       CIx   = FuncAdd(CIx, FuncMul(LocIx, ur));
801       CIy   = FuncAdd(CIy, FuncMul(LocIy, ur));
802       CIz   = FuncAdd(CIz, FuncMul(LocIz, ur));
803       CIxx  = FuncAdd(CIxx, FuncMul(LocIxx, ur));
804       CIyy  = FuncAdd(CIyy, FuncMul(LocIyy, ur));
805       CIzz  = FuncAdd(CIzz, FuncMul(LocIzz, ur));
806       CIxy  = FuncAdd(CIxy, FuncMul(LocIxy, ur));
807       CIxz  = FuncAdd(CIxz, FuncMul(LocIxz, ur));
808       CIyz  = FuncAdd(CIyz, FuncMul(LocIyz, ur));
809     }
810
811     dim = FuncAdd(dim, FuncMul(CArea, lr));
812     Ix  = FuncAdd(Ix,  FuncMul(CIx, lr));
813     Iy  = FuncAdd(Iy,  FuncMul(CIy, lr));
814     Iz  = FuncAdd(Iz,  FuncMul(CIz, lr));
815     Ixx = FuncAdd(Ixx, FuncMul(CIxx, lr));
816     Iyy = FuncAdd(Iyy, FuncMul(CIyy, lr));
817     Izz = FuncAdd(Izz, FuncMul(CIzz, lr));
818     Ixy = FuncAdd(Ixy, FuncMul(CIxy, lr));
819     Ixz = FuncAdd(Iyz, FuncMul(CIxz, lr));
820     Iyz = FuncAdd(Ixz, FuncMul(CIyz, lr));
821     D.Next();
822   }
823
824   if (Abs(dim) >= EPS_DIM) {
825     Ix /= dim;
826     Iy /= dim;
827     Iz /= dim;
828     g.SetCoord (Ix, Iy, Iz);
829   }
830   else {
831     dim =0.;
832     g.SetCoord (0., 0.,0.);
833   }
834
835   inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
836                     gp_XYZ (-Ixy, Iyy, -Iyz),
837                     gp_XYZ (-Ixz, -Iyz, Izz));
838 }
839
840 static void Compute(const Face& S,
841                     const gp_Pnt& loc,
842                     Standard_Real& dim,
843                     gp_Pnt& g,
844                     gp_Mat& inertia)
845 {
846   Standard_Real  (*FuncAdd)(Standard_Real, Standard_Real);
847   Standard_Real  (*FuncMul)(Standard_Real, Standard_Real);
848
849   FuncAdd = Addition;
850   FuncMul = Multiplication;
851
852   Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
853   dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
854
855   Standard_Real LowerU, UpperU, LowerV, UpperV;
856   S.Bounds (LowerU, UpperU, LowerV, UpperV);
857
858   if(Precision::IsInfinite(LowerU) || Precision::IsInfinite(UpperU) ||
859      Precision::IsInfinite(LowerV) || Precision::IsInfinite(UpperV))
860   {
861     FuncAdd = AdditionInf;
862     FuncMul = MultiplicationInf;
863   }
864
865   Standard_Integer UOrder = Min(S.UIntegrationOrder (),
866     math::GaussPointsMax());
867   Standard_Integer VOrder = Min(S.VIntegrationOrder (),
868     math::GaussPointsMax());   
869   gp_Pnt P;          
870   gp_Vec VNor;   
871   Standard_Real dsi, ds;        
872   Standard_Real ur, um, u, vr, vm, v;
873   Standard_Real x, y, z; 
874   Standard_Real Ixi, Iyi, Izi, Ixxi, Iyyi, Izzi, Ixyi, Ixzi, Iyzi;
875   Standard_Real xloc, yloc, zloc;
876   loc.Coord (xloc, yloc, zloc);
877
878   Standard_Integer i, j;
879   math_Vector GaussPU (1, UOrder);     //gauss points and weights
880   math_Vector GaussWU (1, UOrder);
881   math_Vector GaussPV (1, VOrder);
882   math_Vector GaussWV (1, VOrder);
883
884   //Recuperation des points de Gauss dans le fichier GaussPoints.
885   math::GaussPoints  (UOrder,GaussPU);
886   math::GaussWeights (UOrder,GaussWU);
887   math::GaussPoints  (VOrder,GaussPV);
888   math::GaussWeights (VOrder,GaussWV);
889
890   // Calcul des integrales aux points de gauss :
891   um = 0.5 * FuncAdd(UpperU,  LowerU);
892   vm = 0.5 * FuncAdd(UpperV,  LowerV);
893   ur = 0.5 * FuncAdd(UpperU, -LowerU);
894   vr = 0.5 * FuncAdd(UpperV, -LowerV);
895
896   for (j = 1; j <= VOrder; j++) {
897     v = FuncAdd(vm, FuncMul(vr, GaussPV(j)));
898     dsi = Ixi = Iyi = Izi = Ixxi = Iyyi = Izzi = Ixyi = Ixzi = Iyzi = 0.0;
899
900     for (i = 1; i <= UOrder; i++) {
901       u = FuncAdd(um, FuncMul(ur, GaussPU (i)));
902       S.Normal (u, v, P, VNor); 
903       ds = FuncMul(VNor.Magnitude(), GaussWU (i));
904       P.Coord (x, y, z);
905
906       x = FuncAdd(x, -xloc);
907       y = FuncAdd(y, -yloc);
908       z = FuncAdd(z, -zloc);
909
910       dsi = FuncAdd(dsi, ds);
911
912       const Standard_Real XdS = FuncMul(x, ds);
913       const Standard_Real YdS = FuncMul(y, ds);
914       const Standard_Real ZdS = FuncMul(z, ds);
915
916       Ixi = FuncAdd(Ixi, XdS);
917       Iyi = FuncAdd(Iyi, YdS);
918       Izi = FuncAdd(Izi, ZdS);
919       Ixyi = FuncAdd(Ixyi, FuncMul(x, YdS));
920       Iyzi = FuncAdd(Iyzi, FuncMul(y, ZdS));
921       Ixzi = FuncAdd(Ixzi, FuncMul(x, ZdS));
922       x = Precision::IsInfinite(x) ? Precision::Infinite() : x*x;
923       y = Precision::IsInfinite(y) ? Precision::Infinite() : y*y;
924       z = Precision::IsInfinite(z) ? Precision::Infinite() : z*z;
925       Ixxi = FuncAdd(Ixxi, FuncAdd(YdS, ZdS));
926       Iyyi = FuncAdd(Iyyi, FuncAdd(XdS, ZdS));
927       Izzi = FuncAdd(Izzi, FuncAdd(XdS, YdS));
928     }
929
930     dim   = FuncAdd(dim, FuncMul(dsi, GaussWV (j)));
931     Ix    = FuncAdd(Ix,  FuncMul(Ixi, GaussWV (j)));
932     Iy    = FuncAdd(Iy,  FuncMul(Iyi, GaussWV (j)));
933     Iz    = FuncAdd(Iz,  FuncMul(Izi, GaussWV (j)));
934     Ixx   = FuncAdd(Ixx, FuncMul(Ixxi, GaussWV (j)));
935     Iyy   = FuncAdd(Iyy, FuncMul(Iyyi, GaussWV (j)));
936     Izz   = FuncAdd(Izz, FuncMul(Izzi, GaussWV (j)));
937     Ixy   = FuncAdd(Ixy, FuncMul(Ixyi, GaussWV (j)));
938     Iyz   = FuncAdd(Iyz, FuncMul(Iyzi, GaussWV (j)));
939     Ixz   = FuncAdd(Ixz, FuncMul(Ixzi, GaussWV (j)));
940   }
941
942   vr  = FuncMul(vr, ur);
943   Ixx = FuncMul(vr, Ixx);
944   Iyy = FuncMul(vr, Iyy);
945   Izz = FuncMul(vr, Izz);
946   Ixy = FuncMul(vr, Ixy);
947   Ixz = FuncMul(vr, Ixz);
948   Iyz = FuncMul(vr, Iyz);
949
950   if (Abs(dim) >= EPS_DIM)
951   {
952     Ix    /= dim;
953     Iy    /= dim;
954     Iz    /= dim;
955     dim   *= vr;
956     g.SetCoord (Ix, Iy, Iz);
957   }
958   else
959   {
960     dim =0.;
961     g.SetCoord (0.,0.,0.);
962   }
963
964   inertia = gp_Mat (gp_XYZ ( Ixx, -Ixy, -Ixz),
965                     gp_XYZ (-Ixy,  Iyy, -Iyz),
966                     gp_XYZ (-Ixz, -Iyz,  Izz));
967 }
968
969 GProp_SGProps::GProp_SGProps(){}
970
971 GProp_SGProps::GProp_SGProps (const Face&   S,
972                               const gp_Pnt& SLocation
973                               ) 
974 {
975   SetLocation(SLocation);
976   Perform(S);
977 }
978
979 GProp_SGProps::GProp_SGProps (Face&   S,
980                               Domain& D,
981                               const gp_Pnt& SLocation
982                               ) 
983 {
984   SetLocation(SLocation);
985   Perform(S,D);
986 }
987
988 GProp_SGProps::GProp_SGProps(Face& S, const gp_Pnt& SLocation, const Standard_Real Eps){
989   SetLocation(SLocation);
990   Perform(S, Eps);
991 }
992
993 GProp_SGProps::GProp_SGProps(Face& S, Domain& D, const gp_Pnt& SLocation, const Standard_Real Eps){
994   SetLocation(SLocation);
995   Perform(S, D, Eps);
996 }
997
998 void GProp_SGProps::SetLocation(const gp_Pnt& SLocation){
999   loc = SLocation;
1000 }
1001
1002 void GProp_SGProps::Perform(const Face& S){
1003   Compute(S,loc,dim,g,inertia);
1004   myEpsilon = 1.0;
1005   return;
1006 }
1007
1008 void GProp_SGProps::Perform(Face& S, Domain& D){
1009   Compute(S,D,loc,dim,g,inertia);
1010   myEpsilon = 1.0;
1011   return;
1012 }
1013
1014 Standard_Real GProp_SGProps::Perform(Face& S, const Standard_Real Eps){
1015   return myEpsilon = Compute(S,loc,dim,g,inertia,Eps);
1016 }
1017
1018 Standard_Real GProp_SGProps::Perform(Face& S, Domain& D, const Standard_Real Eps){
1019   return myEpsilon = Compute(S,D,loc,dim,g,inertia,Eps);
1020 }
1021
1022
1023 Standard_Real GProp_SGProps::GetEpsilon(){
1024   return myEpsilon;
1025 }