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