projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0024428: Implementation of LGPL license
[occt.git] / src / GProp / GProp_SGProps.gxx
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//
973c2be1 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.
b311480e 11//
973c2be1 12// Alternatively, this file may be used under the terms of Open CASCADE
13// commercial license or contractual agreement.
b311480e 14
7fd59977 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
26class 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
48static Standard_Real EPS_PARAM = 1.e-12;
49static Standard_Real EPS_DIM = 1.e-20;
50static Standard_Real ERROR_ALGEBR_RATIO = 2.0/3.0;
51
52static Standard_Integer GPM = 61;
53static Standard_Integer SUBS_POWER = 32;
54static Standard_Integer SM = 1953;
55
56static math_Vector LGaussP0(1,GPM);
57static math_Vector LGaussW0(1,GPM);
58static math_Vector LGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
59static math_Vector LGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
60
61static math_Vector* LGaussP[] = {&LGaussP0,&LGaussP1};
62static math_Vector* LGaussW[] = {&LGaussW0,&LGaussW1};
63
64static HMath_Vector L1 = new math_Vector(1,SM,0.0);
65static HMath_Vector L2 = new math_Vector(1,SM,0.0);
66static HMath_Vector DimL = new math_Vector(1,SM,0.0);
67static HMath_Vector ErrL = new math_Vector(1,SM,0.0);
68static HMath_Vector ErrUL = new math_Vector(1,SM,0.0);
69static HMath_Vector IxL = new math_Vector(1,SM,0.0);
70static HMath_Vector IyL = new math_Vector(1,SM,0.0);
71static HMath_Vector IzL = new math_Vector(1,SM,0.0);
72static HMath_Vector IxxL = new math_Vector(1,SM,0.0);
73static HMath_Vector IyyL = new math_Vector(1,SM,0.0);
74static HMath_Vector IzzL = new math_Vector(1,SM,0.0);
75static HMath_Vector IxyL = new math_Vector(1,SM,0.0);
76static HMath_Vector IxzL = new math_Vector(1,SM,0.0);
77static HMath_Vector IyzL = new math_Vector(1,SM,0.0);
78
79static math_Vector UGaussP0(1,GPM);
80static math_Vector UGaussW0(1,GPM);
81static math_Vector UGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
82static math_Vector UGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
83
84static math_Vector* UGaussP[] = {&UGaussP0,&UGaussP1};
85static math_Vector* UGaussW[] = {&UGaussW0,&UGaussW1};
86
87static HMath_Vector U1 = new math_Vector(1,SM,0.0);
88static HMath_Vector U2 = new math_Vector(1,SM,0.0);
89static HMath_Vector DimU = new math_Vector(1,SM,0.0);
90static HMath_Vector ErrU = new math_Vector(1,SM,0.0);
91static HMath_Vector IxU = new math_Vector(1,SM,0.0);
92static HMath_Vector IyU = new math_Vector(1,SM,0.0);
93static HMath_Vector IzU = new math_Vector(1,SM,0.0);
94static HMath_Vector IxxU = new math_Vector(1,SM,0.0);
95static HMath_Vector IyyU = new math_Vector(1,SM,0.0);
96static HMath_Vector IzzU = new math_Vector(1,SM,0.0);
97static HMath_Vector IxyU = new math_Vector(1,SM,0.0);
98static HMath_Vector IxzU = new math_Vector(1,SM,0.0);
99static HMath_Vector IyzU = new math_Vector(1,SM,0.0);
100
c63628e8 101static 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
141static 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
178static inline Standard_Real Multiplication(Standard_Real theMA, Standard_Real theMB)
179{
180 return (theMA * theMB);
181}
182
183static inline Standard_Real Addition(Standard_Real theMA, Standard_Real theMB)
184{
185 return (theMA + theMB);
186}
187
7fd59977 188static Standard_Integer FillIntervalBounds(Standard_Real A,
189 Standard_Real B,
190 const TColStd_Array1OfReal& Knots,
c63628e8 191 HMath_Vector& VA,
7fd59977 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)
258ff83b 199 {
200 if(kn < B)
201 {
202 VA(j++) = VB(k++) = kn;
203 }
204 else
205 {
206 break;
207 }
208 }
7fd59977 209 }
210 VB(k) = B;
211 return k;
212}
213
214static inline Standard_Integer MaxSubs(Standard_Integer n, Standard_Integer coeff = SUBS_POWER){
c63628e8 215 // return n = IntegerLast()/coeff < n? IntegerLast(): n*coeff + 1;
7fd59977 216 return Min((n * coeff + 1),SM);
217}
218
219static Standard_Integer LFillIntervalBounds(Standard_Real A,
220 Standard_Real B,
221 const TColStd_Array1OfReal& Knots,
c63628e8 222 const Standard_Integer NumSubs)
7fd59977 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
245static Standard_Integer UFillIntervalBounds(Standard_Real A,
246 Standard_Real B,
247 const TColStd_Array1OfReal& Knots,
c63628e8 248 const Standard_Integer NumSubs)
7fd59977 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
270static 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,
c63628e8 276 const Standard_Real EpsDim,
277 const Standard_Boolean isErrorCalculation,
7fd59977 278 const Standard_Boolean isVerifyComputation)
279{
c63628e8 280 Standard_Real (*FuncAdd)(Standard_Real, Standard_Real);
281 Standard_Real (*FuncMul)(Standard_Real, Standard_Real);
282
283 FuncAdd = Addition;
284 FuncMul = Multiplication;
285
7fd59977 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;
c63628e8 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
7fd59977 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;
c63628e8 323
7fd59977 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);
c63628e8 339
7fd59977 340 NbUSubs = S.SUIntSubs();
341 TColStd_Array1OfReal UKnots(1,NbUSubs+1);
342 S.UKnots(UKnots);
7fd59977 343
c63628e8 344 while (isNaturalRestriction || D.More())
345 {
346 if(isNaturalRestriction)
347 {
7fd59977 348 NbLGaussP[0] = Min(2*NbUGaussP[0],math::GaussPointsMax());
c63628e8 349 }
350 else
351 {
7fd59977 352 S.Load(D.Value()); ++iD;
353 NbLGaussP[0] = S.LIntOrder(EpsDim);
354 }
355
7fd59977 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);
c63628e8 359
7fd59977 360 NbLSubs = isNaturalRestriction? S.SVIntSubs(): S.LIntSubs();
361
362 TColStd_Array1OfReal LKnots(1,NbLSubs+1);
c63628e8 363 if(isNaturalRestriction)
364 {
7fd59977 365 S.VKnots(LKnots);
366 l1 = BV1; l2 = BV2;
c63628e8 367 }
368 else
369 {
7fd59977 370 S.LKnots(LKnots);
371 l1 = S.FirstParameter(); l2 = S.LastParameter();
372 }
c63628e8 373
7fd59977 374 ErrorL = 0.0;
375 kLEnd = 1; JL = 0;
376 //OCC503(apo): if(Abs(l2-l1) < EPS_PARAM) continue;
c63628e8 377 if(Abs(l2-l1) > EPS_PARAM)
378 {
7fd59977 379 iLSubEnd = LFillIntervalBounds(l1, l2, LKnots, NumSubs);
380 LMaxSubs = MaxSubs(iLSubEnd);
c63628e8 381 if(LMaxSubs > DimL.Vector()->Upper())
382 LMaxSubs = DimL.Vector()->Upper();
383
7fd59977 384 DimL.Init(0.0,1,LMaxSubs); ErrL.Init(0.0,1,LMaxSubs); ErrUL.Init(0.0,1,LMaxSubs);
c63628e8 385
7fd59977 386 do{// while: L
c63628e8 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);
7fd59977 625 }
c63628e8 626
627 if(isNaturalRestriction)
628 break;
629
630 D.Next();
7fd59977 631 }
7fd59977 632
c63628e8 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;
7fd59977 656}
657
658static Standard_Real Compute(Face& S, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia,
c63628e8 659 Standard_Real EpsDim)
7fd59977 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
668static Standard_Real Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia,
c63628e8 669 Standard_Real EpsDim)
7fd59977 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
c63628e8 677static 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);
7fd59977 681
c63628e8 682 FuncAdd = Addition;
683 FuncMul = Multiplication;
7fd59977 684
c63628e8 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;
7fd59977 687
c63628e8 688 Standard_Real x, y, z;
689 Standard_Integer NbCGaussgp_Pnts = 0;
7fd59977 690
c63628e8 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||
7fd59977 695
c63628e8 696 gp_Pnt P; //On the Face
697 gp_Vec VNor;
7fd59977 698
c63628e8 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;
7fd59977 705
7fd59977 706
c63628e8 707 S.Bounds (u1, u2, v1, v2);
7fd59977 708
c63628e8 709 if(Precision::IsInfinite(u1) || Precision::IsInfinite(u2) ||
710 Precision::IsInfinite(v1) || Precision::IsInfinite(v2))
711 {
712 FuncAdd = AdditionInf;
713 FuncMul = MultiplicationInf;
714 }
7fd59977 715
716
c63628e8 717 Standard_Integer NbUGaussgp_Pnts = Min(S.UIntegrationOrder (),
718 math::GaussPointsMax());
719 Standard_Integer NbVGaussgp_Pnts = Min(S.VIntegrationOrder (),
720 math::GaussPointsMax());
7fd59977 721
c63628e8 722 Standard_Integer NbGaussgp_Pnts = Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts);
7fd59977 723
c63628e8 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 =
7fd59977 768 LocIxy = LocIxz = LocIyz = 0.0;
c63628e8 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));
7fd59977 797 }
c63628e8 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));
7fd59977 838}
839
c63628e8 840static 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 }
7fd59977 864
c63628e8 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));
7fd59977 967}
968
969GProp_SGProps::GProp_SGProps(){}
970
971GProp_SGProps::GProp_SGProps (const Face& S,
c63628e8 972 const gp_Pnt& SLocation
973 )
7fd59977 974{
c63628e8 975 SetLocation(SLocation);
976 Perform(S);
7fd59977 977}
978
979GProp_SGProps::GProp_SGProps (Face& S,
980 Domain& D,
c63628e8 981 const gp_Pnt& SLocation
982 )
7fd59977 983{
c63628e8 984 SetLocation(SLocation);
985 Perform(S,D);
7fd59977 986}
987
988GProp_SGProps::GProp_SGProps(Face& S, const gp_Pnt& SLocation, const Standard_Real Eps){
989 SetLocation(SLocation);
990 Perform(S, Eps);
991}
992
993GProp_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
998void GProp_SGProps::SetLocation(const gp_Pnt& SLocation){
999 loc = SLocation;
1000}
1001
1002void GProp_SGProps::Perform(const Face& S){
1003 Compute(S,loc,dim,g,inertia);
1004 myEpsilon = 1.0;
1005 return;
1006}
1007
1008void GProp_SGProps::Perform(Face& S, Domain& D){
1009 Compute(S,D,loc,dim,g,inertia);
1010 myEpsilon = 1.0;
1011 return;
1012}
1013
1014Standard_Real GProp_SGProps::Perform(Face& S, const Standard_Real Eps){
1015 return myEpsilon = Compute(S,loc,dim,g,inertia,Eps);
1016}
1017
1018Standard_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
1023Standard_Real GProp_SGProps::GetEpsilon(){
1024 return myEpsilon;
1025}
Open CASCADE Technology repository