0025934: Review PERF test group
[occt.git] / src / BRepGProp / BRepGProp_Sinert.cxx
CommitLineData
b311480e 1// Copyright (c) 1995-1999 Matra Datavision
973c2be1 2// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 3//
973c2be1 4// This file is part of Open CASCADE Technology software library.
b311480e 5//
d5f74e42 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
973c2be1 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.
b311480e 11//
973c2be1 12// Alternatively, this file may be used under the terms of Open CASCADE
13// commercial license or contractual agreement.
b311480e 14
424cd6bb 15#include <BRepGProp_Sinert.ixx>
7fd59977 16
424cd6bb 17#include <math.hxx>
7fd59977 18#include <TColStd_Array1OfReal.hxx>
19#include <Precision.hxx>
20
21class 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
43static Standard_Real EPS_PARAM = 1.e-12;
44static Standard_Real EPS_DIM = 1.e-20;
45static Standard_Real ERROR_ALGEBR_RATIO = 2.0/3.0;
46
47static Standard_Integer GPM = 61;
48static Standard_Integer SUBS_POWER = 32;
49static Standard_Integer SM = 1953;
50
51static math_Vector LGaussP0(1,GPM);
52static math_Vector LGaussW0(1,GPM);
53static math_Vector LGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
54static math_Vector LGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
55
56static math_Vector* LGaussP[] = {&LGaussP0,&LGaussP1};
57static math_Vector* LGaussW[] = {&LGaussW0,&LGaussW1};
58
59static HMath_Vector L1 = new math_Vector(1,SM,0.0);
60static HMath_Vector L2 = new math_Vector(1,SM,0.0);
61static HMath_Vector DimL = new math_Vector(1,SM,0.0);
62static HMath_Vector ErrL = new math_Vector(1,SM,0.0);
63static HMath_Vector ErrUL = new math_Vector(1,SM,0.0);
64static HMath_Vector IxL = new math_Vector(1,SM,0.0);
65static HMath_Vector IyL = new math_Vector(1,SM,0.0);
66static HMath_Vector IzL = new math_Vector(1,SM,0.0);
67static HMath_Vector IxxL = new math_Vector(1,SM,0.0);
68static HMath_Vector IyyL = new math_Vector(1,SM,0.0);
69static HMath_Vector IzzL = new math_Vector(1,SM,0.0);
70static HMath_Vector IxyL = new math_Vector(1,SM,0.0);
71static HMath_Vector IxzL = new math_Vector(1,SM,0.0);
72static HMath_Vector IyzL = new math_Vector(1,SM,0.0);
73
74static math_Vector UGaussP0(1,GPM);
75static math_Vector UGaussW0(1,GPM);
76static math_Vector UGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
77static math_Vector UGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
78
79static math_Vector* UGaussP[] = {&UGaussP0,&UGaussP1};
80static math_Vector* UGaussW[] = {&UGaussW0,&UGaussW1};
81
82static HMath_Vector U1 = new math_Vector(1,SM,0.0);
83static HMath_Vector U2 = new math_Vector(1,SM,0.0);
84static HMath_Vector DimU = new math_Vector(1,SM,0.0);
85static HMath_Vector ErrU = new math_Vector(1,SM,0.0);
86static HMath_Vector IxU = new math_Vector(1,SM,0.0);
87static HMath_Vector IyU = new math_Vector(1,SM,0.0);
88static HMath_Vector IzU = new math_Vector(1,SM,0.0);
89static HMath_Vector IxxU = new math_Vector(1,SM,0.0);
90static HMath_Vector IyyU = new math_Vector(1,SM,0.0);
91static HMath_Vector IzzU = new math_Vector(1,SM,0.0);
92static HMath_Vector IxyU = new math_Vector(1,SM,0.0);
93static HMath_Vector IxzU = new math_Vector(1,SM,0.0);
94static HMath_Vector IyzU = new math_Vector(1,SM,0.0);
95
c63628e8 96static 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
136static 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
173static inline Standard_Real Multiplication(Standard_Real theMA, Standard_Real theMB)
174{
175 return (theMA * theMB);
176}
177
178static inline Standard_Real Addition(Standard_Real theMA, Standard_Real theMB)
179{
180 return (theMA + theMB);
181}
182
7fd59977 183static Standard_Integer FillIntervalBounds(Standard_Real A,
184 Standard_Real B,
185 const TColStd_Array1OfReal& Knots,
c63628e8 186 HMath_Vector& VA,
7fd59977 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)
258ff83b 194 {
195 if(kn < B)
196 {
197 VA(j++) = VB(k++) = kn;
198 }
199 else
200 {
201 break;
202 }
203 }
7fd59977 204 }
205 VB(k) = B;
206 return k;
207}
208
209static inline Standard_Integer MaxSubs(Standard_Integer n, Standard_Integer coeff = SUBS_POWER){
c63628e8 210 // return n = IntegerLast()/coeff < n? IntegerLast(): n*coeff + 1;
7fd59977 211 return Min((n * coeff + 1),SM);
212}
213
214static Standard_Integer LFillIntervalBounds(Standard_Real A,
215 Standard_Real B,
216 const TColStd_Array1OfReal& Knots,
c63628e8 217 const Standard_Integer NumSubs)
7fd59977 218{
e97a6598 219 Standard_Integer iEnd = MaxSubs(Knots.Upper()-1, NumSubs), jEnd = L1->Upper();
220 iEnd = Max(iEnd, Knots.Upper());
7fd59977 221 if(iEnd - 1 > jEnd){
7fd59977 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
240static Standard_Integer UFillIntervalBounds(Standard_Real A,
241 Standard_Real B,
242 const TColStd_Array1OfReal& Knots,
c63628e8 243 const Standard_Integer NumSubs)
7fd59977 244{
e97a6598 245 Standard_Integer iEnd = MaxSubs(Knots.Upper()-1, NumSubs), jEnd = U1->Upper();
246 iEnd = Max(iEnd, Knots.Upper());
7fd59977 247 if(iEnd - 1 > jEnd){
7fd59977 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
424cd6bb 265static Standard_Real CCompute(BRepGProp_Face& S,
266 BRepGProp_Domain& D,
7fd59977 267 const gp_Pnt& loc,
268 Standard_Real& Dim,
269 gp_Pnt& g,
270 gp_Mat& inertia,
c63628e8 271 const Standard_Real EpsDim,
272 const Standard_Boolean isErrorCalculation,
7fd59977 273 const Standard_Boolean isVerifyComputation)
274{
c63628e8 275 Standard_Real (*FuncAdd)(Standard_Real, Standard_Real);
276 Standard_Real (*FuncMul)(Standard_Real, Standard_Real);
277
278 FuncAdd = Addition;
279 FuncMul = Multiplication;
280
7fd59977 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;
424cd6bb 290 //BRepGProp_Face parametrization in U and V direction
7fd59977 291 Standard_Real BV1, BV2, v;
292 Standard_Real BU1, BU2, u1, u2, um, ur, u;
c63628e8 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
7fd59977 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;
424cd6bb 309 //On the BRepGProp_Face
7fd59977 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;
c63628e8 318
7fd59977 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);
c63628e8 334
7fd59977 335 NbUSubs = S.SUIntSubs();
336 TColStd_Array1OfReal UKnots(1,NbUSubs+1);
337 S.UKnots(UKnots);
7fd59977 338
c63628e8 339 while (isNaturalRestriction || D.More())
340 {
341 if(isNaturalRestriction)
342 {
7fd59977 343 NbLGaussP[0] = Min(2*NbUGaussP[0],math::GaussPointsMax());
c63628e8 344 }
345 else
346 {
7fd59977 347 S.Load(D.Value()); ++iD;
348 NbLGaussP[0] = S.LIntOrder(EpsDim);
349 }
350
7fd59977 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);
c63628e8 354
7fd59977 355 NbLSubs = isNaturalRestriction? S.SVIntSubs(): S.LIntSubs();
356
357 TColStd_Array1OfReal LKnots(1,NbLSubs+1);
c63628e8 358 if(isNaturalRestriction)
359 {
7fd59977 360 S.VKnots(LKnots);
361 l1 = BV1; l2 = BV2;
c63628e8 362 }
363 else
364 {
7fd59977 365 S.LKnots(LKnots);
366 l1 = S.FirstParameter(); l2 = S.LastParameter();
367 }
c63628e8 368
7fd59977 369 ErrorL = 0.0;
370 kLEnd = 1; JL = 0;
371 //OCC503(apo): if(Abs(l2-l1) < EPS_PARAM) continue;
c63628e8 372 if(Abs(l2-l1) > EPS_PARAM)
373 {
7fd59977 374 iLSubEnd = LFillIntervalBounds(l1, l2, LKnots, NumSubs);
375 LMaxSubs = MaxSubs(iLSubEnd);
c63628e8 376 if(LMaxSubs > DimL.Vector()->Upper())
377 LMaxSubs = DimL.Vector()->Upper();
378
7fd59977 379 DimL.Init(0.0,1,LMaxSubs); ErrL.Init(0.0,1,LMaxSubs); ErrUL.Init(0.0,1,LMaxSubs);
c63628e8 380
7fd59977 381 do{// while: L
c63628e8 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));
a24c75d9 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));
c63628e8 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);
7fd59977 622 }
c63628e8 623
624 if(isNaturalRestriction)
625 break;
626
627 D.Next();
7fd59977 628 }
7fd59977 629
c63628e8 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;
7fd59977 653}
654
424cd6bb 655static Standard_Real Compute(BRepGProp_Face& S, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia,
c63628e8 656 Standard_Real EpsDim)
7fd59977 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);
424cd6bb 661 BRepGProp_Domain D;
7fd59977 662 return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
663}
664
424cd6bb 665static Standard_Real Compute(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia,
c63628e8 666 Standard_Real EpsDim)
7fd59977 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
424cd6bb 674static void Compute(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pnt& loc, Standard_Real& dim, gp_Pnt& g, gp_Mat& inertia)
c63628e8 675{
676 Standard_Real (*FuncAdd)(Standard_Real, Standard_Real);
677 Standard_Real (*FuncMul)(Standard_Real, Standard_Real);
7fd59977 678
c63628e8 679 FuncAdd = Addition;
680 FuncMul = Multiplication;
7fd59977 681
c63628e8 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;
7fd59977 684
c63628e8 685 Standard_Real x, y, z;
686 Standard_Integer NbCGaussgp_Pnts = 0;
7fd59977 687
c63628e8 688 Standard_Real l1, l2, lm, lr, l; //boundary curve parametrization
424cd6bb 689 Standard_Real v1, v2, v; //BRepGProp_Face parametrization in v direction
c63628e8 690 Standard_Real u1, u2, um, ur, u;
691 Standard_Real ds; //Jacobien (x, y, z) -> (u, v) = ||n||
7fd59977 692
424cd6bb 693 gp_Pnt P; //On the BRepGProp_Face
c63628e8 694 gp_Vec VNor;
7fd59977 695
c63628e8 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;
7fd59977 702
7fd59977 703
c63628e8 704 S.Bounds (u1, u2, v1, v2);
7fd59977 705
c63628e8 706 if(Precision::IsInfinite(u1) || Precision::IsInfinite(u2) ||
707 Precision::IsInfinite(v1) || Precision::IsInfinite(v2))
708 {
709 FuncAdd = AdditionInf;
710 FuncMul = MultiplicationInf;
711 }
7fd59977 712
713
c63628e8 714 Standard_Integer NbUGaussgp_Pnts = Min(S.UIntegrationOrder (),
715 math::GaussPointsMax());
716 Standard_Integer NbVGaussgp_Pnts = Min(S.VIntegrationOrder (),
717 math::GaussPointsMax());
7fd59977 718
c63628e8 719 Standard_Integer NbGaussgp_Pnts = Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts);
7fd59977 720
c63628e8 721 //Number of Gauss points for the integration
424cd6bb 722 //on the BRepGProp_Face
c63628e8 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
c63628e8 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 =
7fd59977 760 LocIxy = LocIxz = LocIyz = 0.0;
c63628e8 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));
a24c75d9 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));
7fd59977 791 }
c63628e8 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));
7fd59977 832}
833
424cd6bb 834static void Compute(const BRepGProp_Face& S,
c63628e8 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 }
7fd59977 858
c63628e8 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));
a24c75d9 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));
c63628e8 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));
7fd59977 963}
964
424cd6bb 965BRepGProp_Sinert::BRepGProp_Sinert(){}
7fd59977 966
424cd6bb 967BRepGProp_Sinert::BRepGProp_Sinert (const BRepGProp_Face& S,
c63628e8 968 const gp_Pnt& SLocation
969 )
7fd59977 970{
c63628e8 971 SetLocation(SLocation);
972 Perform(S);
7fd59977 973}
974
424cd6bb 975BRepGProp_Sinert::BRepGProp_Sinert (BRepGProp_Face& S,
976 BRepGProp_Domain& D,
c63628e8 977 const gp_Pnt& SLocation
978 )
7fd59977 979{
c63628e8 980 SetLocation(SLocation);
981 Perform(S,D);
7fd59977 982}
983
424cd6bb 984BRepGProp_Sinert::BRepGProp_Sinert(BRepGProp_Face& S, const gp_Pnt& SLocation, const Standard_Real Eps){
7fd59977 985 SetLocation(SLocation);
986 Perform(S, Eps);
987}
988
424cd6bb 989BRepGProp_Sinert::BRepGProp_Sinert(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pnt& SLocation, const Standard_Real Eps){
7fd59977 990 SetLocation(SLocation);
991 Perform(S, D, Eps);
992}
993
424cd6bb 994void BRepGProp_Sinert::SetLocation(const gp_Pnt& SLocation){
7fd59977 995 loc = SLocation;
996}
997
424cd6bb 998void BRepGProp_Sinert::Perform(const BRepGProp_Face& S){
7fd59977 999 Compute(S,loc,dim,g,inertia);
1000 myEpsilon = 1.0;
1001 return;
1002}
1003
424cd6bb 1004void BRepGProp_Sinert::Perform(BRepGProp_Face& S, BRepGProp_Domain& D){
7fd59977 1005 Compute(S,D,loc,dim,g,inertia);
1006 myEpsilon = 1.0;
1007 return;
1008}
1009
424cd6bb 1010Standard_Real BRepGProp_Sinert::Perform(BRepGProp_Face& S, const Standard_Real Eps){
7fd59977 1011 return myEpsilon = Compute(S,loc,dim,g,inertia,Eps);
1012}
1013
424cd6bb 1014Standard_Real BRepGProp_Sinert::Perform(BRepGProp_Face& S, BRepGProp_Domain& D, const Standard_Real Eps){
7fd59977 1015 return myEpsilon = Compute(S,D,loc,dim,g,inertia,Eps);
1016}
1017
1018
424cd6bb 1019Standard_Real BRepGProp_Sinert::GetEpsilon(){
7fd59977 1020 return myEpsilon;
1021}