b311480e |
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 | |
7fd59977 |
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 | |
30 | class 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 | |
52 | static Standard_Real EPS_PARAM = 1.e-12; |
53 | static Standard_Real EPS_DIM = 1.e-20; |
54 | static Standard_Real ERROR_ALGEBR_RATIO = 2.0/3.0; |
55 | |
56 | static Standard_Integer GPM = 61; |
57 | static Standard_Integer SUBS_POWER = 32; |
58 | static Standard_Integer SM = 1953; |
59 | |
60 | static math_Vector LGaussP0(1,GPM); |
61 | static math_Vector LGaussW0(1,GPM); |
62 | static math_Vector LGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM))); |
63 | static math_Vector LGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM))); |
64 | |
65 | static math_Vector* LGaussP[] = {&LGaussP0,&LGaussP1}; |
66 | static math_Vector* LGaussW[] = {&LGaussW0,&LGaussW1}; |
67 | |
68 | static HMath_Vector L1 = new math_Vector(1,SM,0.0); |
69 | static HMath_Vector L2 = new math_Vector(1,SM,0.0); |
70 | static HMath_Vector DimL = new math_Vector(1,SM,0.0); |
71 | static HMath_Vector ErrL = new math_Vector(1,SM,0.0); |
72 | static HMath_Vector ErrUL = new math_Vector(1,SM,0.0); |
73 | static HMath_Vector IxL = new math_Vector(1,SM,0.0); |
74 | static HMath_Vector IyL = new math_Vector(1,SM,0.0); |
75 | static HMath_Vector IzL = new math_Vector(1,SM,0.0); |
76 | static HMath_Vector IxxL = new math_Vector(1,SM,0.0); |
77 | static HMath_Vector IyyL = new math_Vector(1,SM,0.0); |
78 | static HMath_Vector IzzL = new math_Vector(1,SM,0.0); |
79 | static HMath_Vector IxyL = new math_Vector(1,SM,0.0); |
80 | static HMath_Vector IxzL = new math_Vector(1,SM,0.0); |
81 | static HMath_Vector IyzL = new math_Vector(1,SM,0.0); |
82 | |
83 | static math_Vector UGaussP0(1,GPM); |
84 | static math_Vector UGaussW0(1,GPM); |
85 | static math_Vector UGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM))); |
86 | static math_Vector UGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM))); |
87 | |
88 | static math_Vector* UGaussP[] = {&UGaussP0,&UGaussP1}; |
89 | static math_Vector* UGaussW[] = {&UGaussW0,&UGaussW1}; |
90 | |
91 | static HMath_Vector U1 = new math_Vector(1,SM,0.0); |
92 | static HMath_Vector U2 = new math_Vector(1,SM,0.0); |
93 | static HMath_Vector DimU = new math_Vector(1,SM,0.0); |
94 | static HMath_Vector ErrU = new math_Vector(1,SM,0.0); |
95 | static HMath_Vector IxU = new math_Vector(1,SM,0.0); |
96 | static HMath_Vector IyU = new math_Vector(1,SM,0.0); |
97 | static HMath_Vector IzU = new math_Vector(1,SM,0.0); |
98 | static HMath_Vector IxxU = new math_Vector(1,SM,0.0); |
99 | static HMath_Vector IyyU = new math_Vector(1,SM,0.0); |
100 | static HMath_Vector IzzU = new math_Vector(1,SM,0.0); |
101 | static HMath_Vector IxyU = new math_Vector(1,SM,0.0); |
102 | static HMath_Vector IxzU = new math_Vector(1,SM,0.0); |
103 | static HMath_Vector IyzU = new math_Vector(1,SM,0.0); |
104 | |
c63628e8 |
105 | static 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 | |
145 | static 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 | |
182 | static inline Standard_Real Multiplication(Standard_Real theMA, Standard_Real theMB) |
183 | { |
184 | return (theMA * theMB); |
185 | } |
186 | |
187 | static inline Standard_Real Addition(Standard_Real theMA, Standard_Real theMB) |
188 | { |
189 | return (theMA + theMB); |
190 | } |
191 | |
7fd59977 |
192 | static Standard_Integer FillIntervalBounds(Standard_Real A, |
193 | Standard_Real B, |
194 | const TColStd_Array1OfReal& Knots, |
c63628e8 |
195 | HMath_Vector& VA, |
7fd59977 |
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) |
258ff83b |
203 | { |
204 | if(kn < B) |
205 | { |
206 | VA(j++) = VB(k++) = kn; |
207 | } |
208 | else |
209 | { |
210 | break; |
211 | } |
212 | } |
7fd59977 |
213 | } |
214 | VB(k) = B; |
215 | return k; |
216 | } |
217 | |
218 | static inline Standard_Integer MaxSubs(Standard_Integer n, Standard_Integer coeff = SUBS_POWER){ |
c63628e8 |
219 | // return n = IntegerLast()/coeff < n? IntegerLast(): n*coeff + 1; |
7fd59977 |
220 | return Min((n * coeff + 1),SM); |
221 | } |
222 | |
223 | static Standard_Integer LFillIntervalBounds(Standard_Real A, |
224 | Standard_Real B, |
225 | const TColStd_Array1OfReal& Knots, |
c63628e8 |
226 | const Standard_Integer NumSubs) |
7fd59977 |
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 | |
249 | static Standard_Integer UFillIntervalBounds(Standard_Real A, |
250 | Standard_Real B, |
251 | const TColStd_Array1OfReal& Knots, |
c63628e8 |
252 | const Standard_Integer NumSubs) |
7fd59977 |
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 | |
274 | static 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, |
c63628e8 |
280 | const Standard_Real EpsDim, |
281 | const Standard_Boolean isErrorCalculation, |
7fd59977 |
282 | const Standard_Boolean isVerifyComputation) |
283 | { |
c63628e8 |
284 | Standard_Real (*FuncAdd)(Standard_Real, Standard_Real); |
285 | Standard_Real (*FuncMul)(Standard_Real, Standard_Real); |
286 | |
287 | FuncAdd = Addition; |
288 | FuncMul = Multiplication; |
289 | |
7fd59977 |
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; |
c63628e8 |
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 | |
7fd59977 |
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; |
c63628e8 |
327 | |
7fd59977 |
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); |
c63628e8 |
343 | |
7fd59977 |
344 | NbUSubs = S.SUIntSubs(); |
345 | TColStd_Array1OfReal UKnots(1,NbUSubs+1); |
346 | S.UKnots(UKnots); |
7fd59977 |
347 | |
c63628e8 |
348 | while (isNaturalRestriction || D.More()) |
349 | { |
350 | if(isNaturalRestriction) |
351 | { |
7fd59977 |
352 | NbLGaussP[0] = Min(2*NbUGaussP[0],math::GaussPointsMax()); |
c63628e8 |
353 | } |
354 | else |
355 | { |
7fd59977 |
356 | S.Load(D.Value()); ++iD; |
357 | NbLGaussP[0] = S.LIntOrder(EpsDim); |
358 | } |
359 | |
7fd59977 |
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); |
c63628e8 |
363 | |
7fd59977 |
364 | NbLSubs = isNaturalRestriction? S.SVIntSubs(): S.LIntSubs(); |
365 | |
366 | TColStd_Array1OfReal LKnots(1,NbLSubs+1); |
c63628e8 |
367 | if(isNaturalRestriction) |
368 | { |
7fd59977 |
369 | S.VKnots(LKnots); |
370 | l1 = BV1; l2 = BV2; |
c63628e8 |
371 | } |
372 | else |
373 | { |
7fd59977 |
374 | S.LKnots(LKnots); |
375 | l1 = S.FirstParameter(); l2 = S.LastParameter(); |
376 | } |
c63628e8 |
377 | |
7fd59977 |
378 | ErrorL = 0.0; |
379 | kLEnd = 1; JL = 0; |
380 | //OCC503(apo): if(Abs(l2-l1) < EPS_PARAM) continue; |
c63628e8 |
381 | if(Abs(l2-l1) > EPS_PARAM) |
382 | { |
7fd59977 |
383 | iLSubEnd = LFillIntervalBounds(l1, l2, LKnots, NumSubs); |
384 | LMaxSubs = MaxSubs(iLSubEnd); |
c63628e8 |
385 | if(LMaxSubs > DimL.Vector()->Upper()) |
386 | LMaxSubs = DimL.Vector()->Upper(); |
387 | |
7fd59977 |
388 | DimL.Init(0.0,1,LMaxSubs); ErrL.Init(0.0,1,LMaxSubs); ErrUL.Init(0.0,1,LMaxSubs); |
c63628e8 |
389 | |
7fd59977 |
390 | do{// while: L |
c63628e8 |
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); |
7fd59977 |
629 | } |
c63628e8 |
630 | |
631 | if(isNaturalRestriction) |
632 | break; |
633 | |
634 | D.Next(); |
7fd59977 |
635 | } |
7fd59977 |
636 | |
c63628e8 |
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; |
7fd59977 |
660 | } |
661 | |
662 | static Standard_Real Compute(Face& S, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, |
c63628e8 |
663 | Standard_Real EpsDim) |
7fd59977 |
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 | |
672 | static Standard_Real Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, |
c63628e8 |
673 | Standard_Real EpsDim) |
7fd59977 |
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 | |
c63628e8 |
681 | static 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); |
7fd59977 |
685 | |
c63628e8 |
686 | FuncAdd = Addition; |
687 | FuncMul = Multiplication; |
7fd59977 |
688 | |
c63628e8 |
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; |
7fd59977 |
691 | |
c63628e8 |
692 | Standard_Real x, y, z; |
693 | Standard_Integer NbCGaussgp_Pnts = 0; |
7fd59977 |
694 | |
c63628e8 |
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|| |
7fd59977 |
699 | |
c63628e8 |
700 | gp_Pnt P; //On the Face |
701 | gp_Vec VNor; |
7fd59977 |
702 | |
c63628e8 |
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; |
7fd59977 |
709 | |
7fd59977 |
710 | |
c63628e8 |
711 | S.Bounds (u1, u2, v1, v2); |
7fd59977 |
712 | |
c63628e8 |
713 | if(Precision::IsInfinite(u1) || Precision::IsInfinite(u2) || |
714 | Precision::IsInfinite(v1) || Precision::IsInfinite(v2)) |
715 | { |
716 | FuncAdd = AdditionInf; |
717 | FuncMul = MultiplicationInf; |
718 | } |
7fd59977 |
719 | |
720 | |
c63628e8 |
721 | Standard_Integer NbUGaussgp_Pnts = Min(S.UIntegrationOrder (), |
722 | math::GaussPointsMax()); |
723 | Standard_Integer NbVGaussgp_Pnts = Min(S.VIntegrationOrder (), |
724 | math::GaussPointsMax()); |
7fd59977 |
725 | |
c63628e8 |
726 | Standard_Integer NbGaussgp_Pnts = Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts); |
7fd59977 |
727 | |
c63628e8 |
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 = |
7fd59977 |
772 | LocIxy = LocIxz = LocIyz = 0.0; |
c63628e8 |
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)); |
7fd59977 |
801 | } |
c63628e8 |
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)); |
7fd59977 |
842 | } |
843 | |
c63628e8 |
844 | static 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 | } |
7fd59977 |
868 | |
c63628e8 |
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)); |
7fd59977 |
971 | } |
972 | |
973 | GProp_SGProps::GProp_SGProps(){} |
974 | |
975 | GProp_SGProps::GProp_SGProps (const Face& S, |
c63628e8 |
976 | const gp_Pnt& SLocation |
977 | ) |
7fd59977 |
978 | { |
c63628e8 |
979 | SetLocation(SLocation); |
980 | Perform(S); |
7fd59977 |
981 | } |
982 | |
983 | GProp_SGProps::GProp_SGProps (Face& S, |
984 | Domain& D, |
c63628e8 |
985 | const gp_Pnt& SLocation |
986 | ) |
7fd59977 |
987 | { |
c63628e8 |
988 | SetLocation(SLocation); |
989 | Perform(S,D); |
7fd59977 |
990 | } |
991 | |
992 | GProp_SGProps::GProp_SGProps(Face& S, const gp_Pnt& SLocation, const Standard_Real Eps){ |
993 | SetLocation(SLocation); |
994 | Perform(S, Eps); |
995 | } |
996 | |
997 | GProp_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 | |
1002 | void GProp_SGProps::SetLocation(const gp_Pnt& SLocation){ |
1003 | loc = SLocation; |
1004 | } |
1005 | |
1006 | void GProp_SGProps::Perform(const Face& S){ |
1007 | Compute(S,loc,dim,g,inertia); |
1008 | myEpsilon = 1.0; |
1009 | return; |
1010 | } |
1011 | |
1012 | void GProp_SGProps::Perform(Face& S, Domain& D){ |
1013 | Compute(S,D,loc,dim,g,inertia); |
1014 | myEpsilon = 1.0; |
1015 | return; |
1016 | } |
1017 | |
1018 | Standard_Real GProp_SGProps::Perform(Face& S, const Standard_Real Eps){ |
1019 | return myEpsilon = Compute(S,loc,dim,g,inertia,Eps); |
1020 | } |
1021 | |
1022 | Standard_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 | |
1027 | Standard_Real GProp_SGProps::GetEpsilon(){ |
1028 | return myEpsilon; |
1029 | } |