1 // Copyright (c) 2008-2015 OPEN CASCADE SAS
2 //
3 // This file is part of Open CASCADE Technology software library.
4 //
5 // This library is free software; you can redistribute it and/or modify it under
6 // the terms of the GNU Lesser General Public License version 2.1 as published
7 // by the Free Software Foundation, with special exception defined in the file
8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
9 // distribution for complete text of the license and disclaimer of any warranty.
10 //
11 // Alternatively, this file may be used under the terms of Open CASCADE
12 // commercial license or contractual agreement.
13
14 #include <math.hxx>
15 #include <Precision.hxx>
16 #include <TColStd_Array1OfReal.hxx>
17 #include <Standard_Assert.hxx>
18 #include <BRepGProp_Face.hxx>
19 #include <BRepGProp_Domain.hxx>
20 #include <BRepGProp_Gauss.hxx>
21
22 // If the following is defined the error of algorithm is calculated by static moments
23 #define IS_MIN_DIM
24
25 namespace
26 {
27 // Minimal value of interval's range for computation | minimal value of "dim" | ...
28 static const Standard_Real EPS_PARAM = 1.e-12;
29 static const Standard_Real EPS_DIM = 1.e-30;
30 static const Standard_Real ERROR_ALGEBR_RATIO = 2.0 / 3.0;
31
32 // Maximum of GaussPoints on a subinterval and maximum of subintervals
33 static const Standard_Integer GPM = math::GaussPointsMax();
34 static const Standard_Integer SUBS_POWER = 32;
35 static const Standard_Integer SM = SUBS_POWER * GPM + 1;
36
37 // Auxiliary inner functions to perform arithmetic operations.
Add(const Standard_Real theA,const Standard_Real theB)38 static Standard_Real Add(const Standard_Real theA, const Standard_Real theB)
39 {
40 return theA + theB;
41 }
42
AddInf(const Standard_Real theA,const Standard_Real theB)43 static Standard_Real AddInf(const Standard_Real theA, const Standard_Real theB)
44 {
45 if (Precision::IsPositiveInfinite(theA))
46 {
47 if (Precision::IsNegativeInfinite(theB))
48 return 0.0;
49 else
50 return Precision::Infinite();
51 }
52
53 if (Precision::IsPositiveInfinite(theB))
54 {
55 if (Precision::IsNegativeInfinite(theA))
56 return 0.0;
57 else
58 return Precision::Infinite();
59 }
60
61 if (Precision::IsNegativeInfinite(theA))
62 {
63 if (Precision::IsPositiveInfinite(theB))
64 return 0.0;
65 else
66 return -Precision::Infinite();
67 }
68
69 if (Precision::IsNegativeInfinite(theB))
70 {
71 if (Precision::IsPositiveInfinite(theA))
72 return 0.0;
73 else
74 return -Precision::Infinite();
75 }
76
77 return theA + theB;
78 }
79
Mult(const Standard_Real theA,const Standard_Real theB)80 static Standard_Real Mult(const Standard_Real theA, const Standard_Real theB)
81 {
82 return theA * theB;
83 }
84
MultInf(const Standard_Real theA,const Standard_Real theB)85 static Standard_Real MultInf(const Standard_Real theA, const Standard_Real theB)
86 {
87 if ((theA == 0.0) || (theB == 0.0)) //strictly zerro (without any tolerances)
88 return 0.0;
89
90 if (Precision::IsPositiveInfinite(theA))
91 {
92 if (theB < 0.0)
93 return -Precision::Infinite();
94 else
95 return Precision::Infinite();
96 }
97
98 if (Precision::IsPositiveInfinite(theB))
99 {
100 if (theA < 0.0)
101 return -Precision::Infinite();
102 else
103 return Precision::Infinite();
104 }
105
106 if (Precision::IsNegativeInfinite(theA))
107 {
108 if (theB < 0.0)
109 return +Precision::Infinite();
110 else
111 return -Precision::Infinite();
112 }
113
114 if (Precision::IsNegativeInfinite(theB))
115 {
116 if (theA < 0.0)
117 return +Precision::Infinite();
118 else
119 return -Precision::Infinite();
120 }
121
122 return theA * theB;
123 }
124 }
125
126 //=======================================================================
127 //function : BRepGProp_Gauss::Inert::Inert
128 //purpose : Constructor
129 //=======================================================================
Inertia()130 BRepGProp_Gauss::Inertia::Inertia()
131 : Mass(0.0),
132 Ix (0.0),
133 Iy (0.0),
134 Iz (0.0),
135 Ixx (0.0),
136 Iyy (0.0),
137 Izz (0.0),
138 Ixy (0.0),
139 Ixz (0.0),
140 Iyz (0.0)
141 {
142 }
143
144 //=======================================================================
145 //function : Inertia::Reset
146 //purpose : Zeroes all values.
147 //=======================================================================
Reset()148 void BRepGProp_Gauss::Inertia::Reset()
149 {
150 memset(reinterpret_cast<void*>(this), 0, sizeof(BRepGProp_Gauss::Inertia));
151 }
152
153 //=======================================================================
154 //function : BRepGProp_Gauss
155 //purpose : Constructor
156 //=======================================================================
BRepGProp_Gauss(const BRepGProp_GaussType theType)157 BRepGProp_Gauss::BRepGProp_Gauss(const BRepGProp_GaussType theType)
158 : myType(theType)
159 {
160 add = (::Add );
161 mult = (::Mult);
162 }
163
164 //=======================================================================
165 //function : MaxSubs
166 //purpose :
167 //=======================================================================
MaxSubs(const Standard_Integer theN,const Standard_Integer theCoeff)168 Standard_Integer BRepGProp_Gauss::MaxSubs(const Standard_Integer theN,
169 const Standard_Integer theCoeff)
170 {
171 return IntegerLast() / theCoeff < theN ?
172 IntegerLast() : theN * theCoeff + 1;
173 }
174
175 //=======================================================================
176 //function : Init
177 //purpose :
178 //=======================================================================
Init(NCollection_Handle<math_Vector> & theOutVec,const Standard_Real theValue,const Standard_Integer theFirst,const Standard_Integer theLast)179 void BRepGProp_Gauss::Init(NCollection_Handle<math_Vector>& theOutVec,
180 const Standard_Real theValue,
181 const Standard_Integer theFirst,
182 const Standard_Integer theLast)
183 {
184 if(theLast - theFirst == 0)
185 {
186 theOutVec->Init(theValue);
187 }
188 else
189 {
190 for (Standard_Integer i = theFirst; i <= theLast; ++i)
191 theOutVec->Value(i) = theValue;
192 }
193 }
194
195 //=======================================================================
196 //function : InitMass
197 //purpose :
198 //=======================================================================
InitMass(const Standard_Real theValue,const Standard_Integer theFirst,const Standard_Integer theLast,InertiaArray & theArray)199 void BRepGProp_Gauss::InitMass(const Standard_Real theValue,
200 const Standard_Integer theFirst,
201 const Standard_Integer theLast,
202 InertiaArray& theArray)
203 {
204 if (theArray.IsNull())
205 return;
206
207 Standard_Integer aFirst = theFirst;
208 Standard_Integer aLast = theLast;
209
210 if (theLast - theFirst == 0)
211 {
212 aFirst = theArray->Lower();
213 aLast = theArray->Upper();
214 }
215
216 for (Standard_Integer i = aFirst; i <= aLast; ++i)
217 theArray->ChangeValue(i).Mass = theValue;
218 }
219
220 //=======================================================================
221 //function : FillIntervalBounds
222 //purpose :
223 //=======================================================================
FillIntervalBounds(const Standard_Real theA,const Standard_Real theB,const TColStd_Array1OfReal & theKnots,const Standard_Integer theNumSubs,InertiaArray & theInerts,NCollection_Handle<math_Vector> & theParam1,NCollection_Handle<math_Vector> & theParam2,NCollection_Handle<math_Vector> & theError,NCollection_Handle<math_Vector> & theCommonError)224 Standard_Integer BRepGProp_Gauss::FillIntervalBounds(
225 const Standard_Real theA,
226 const Standard_Real theB,
227 const TColStd_Array1OfReal& theKnots,
228 const Standard_Integer theNumSubs,
229 InertiaArray& theInerts,
230 NCollection_Handle<math_Vector>& theParam1,
231 NCollection_Handle<math_Vector>& theParam2,
232 NCollection_Handle<math_Vector>& theError,
233 NCollection_Handle<math_Vector>& theCommonError)
234 {
235 const Standard_Integer aSize =
236 Max(theKnots.Upper(), MaxSubs(theKnots.Upper() - 1, theNumSubs));
237
238 if (aSize - 1 > theParam1->Upper())
239 {
240 theInerts = new NCollection_Array1<Inertia>(1, aSize);
241 theParam1 = new math_Vector(1, aSize);
242 theParam2 = new math_Vector(1, aSize);
243 theError = new math_Vector(1, aSize, 0.0);
244
245 if (theCommonError.IsNull() == Standard_False)
246 theCommonError = new math_Vector(1, aSize, 0.0);
247 }
248
249 Standard_Integer j = 1, k = 1;
250 theParam1->Value(j++) = theA;
251
252 const Standard_Integer aLength = theKnots.Upper();
253 for (Standard_Integer i = 1; i <= aLength; ++i)
254 {
255 const Standard_Real kn = theKnots(i);
256 if (theA < kn)
257 {
258 if (kn < theB)
259 {
260 theParam1->Value(j++) = kn;
261 theParam2->Value(k++) = kn;
262 }
263 else
264 break;
265 }
266 }
267
268 theParam2->Value(k) = theB;
269 return k;
270 }
271
272
273
274 //=======================================================================
275 //function : computeVInertiaOfElementaryPart
276 //purpose :
277 //=======================================================================
computeVInertiaOfElementaryPart(const gp_Pnt & thePoint,const gp_Vec & theNormal,const gp_Pnt & theLocation,const Standard_Real theWeight,const Standard_Real theCoeff[],const Standard_Boolean theIsByPoint,BRepGProp_Gauss::Inertia & theOutInertia)278 void BRepGProp_Gauss::computeVInertiaOfElementaryPart(
279 const gp_Pnt& thePoint,
280 const gp_Vec& theNormal,
281 const gp_Pnt& theLocation,
282 const Standard_Real theWeight,
283 const Standard_Real theCoeff[],
284 const Standard_Boolean theIsByPoint,
285 BRepGProp_Gauss::Inertia& theOutInertia)
286 {
287 Standard_Real x = thePoint.X() - theLocation.X();
288 Standard_Real y = thePoint.Y() - theLocation.Y();
289 Standard_Real z = thePoint.Z() - theLocation.Z();
290
291 const Standard_Real xn = theNormal.X() * theWeight;
292 const Standard_Real yn = theNormal.Y() * theWeight;
293 const Standard_Real zn = theNormal.Z() * theWeight;
294
295 if (theIsByPoint)
296 {
297 ///////////////////// ///////////////////////
298 // OFV code // // Initial code //
299 ///////////////////// ///////////////////////
300 // modified by APO
301
302 Standard_Real dv = x * xn + y * yn + z * zn; //xyz = x * y * z;
303 theOutInertia.Mass += dv / 3.0; //Ixyi += zn * xyz;
304 theOutInertia.Ix += 0.25 * x * dv; //Iyzi += xn * xyz;
305 theOutInertia.Iy += 0.25 * y * dv; //Ixzi += yn * xyz;
306 theOutInertia.Iz += 0.25 * z * dv; //xi = x * x * x * xn / 3.0;
307 x -= theCoeff[0]; //yi = y * y * y * yn / 3.0;
308 y -= theCoeff[1]; //zi = z * z * z * zn / 3.0;
309 z -= theCoeff[2]; //Ixxi += (yi + zi);
310 dv *= 0.2; //Iyyi += (xi + zi);
311 theOutInertia.Ixy -= x * y * dv; //Izzi += (xi + yi);
312 theOutInertia.Iyz -= y * z * dv; //x -= Coeff[0];
313 theOutInertia.Ixz -= x * z * dv; //y -= Coeff[1];
314 x *= x; //z -= Coeff[2];
315 y *= y; //dv = x * xn + y * yn + z * zn;
316 z *= z; //dvi += dv;
317 theOutInertia.Ixx += (y + z) * dv; //Ixi += x * dv;
318 theOutInertia.Iyy += (x + z) * dv; //Iyi += y * dv;
319 theOutInertia.Izz += (x + y) * dv; //Izi += z * dv;
320 }
321 else
322 { // By plane
323 const Standard_Real s = xn * theCoeff[0] + yn * theCoeff[1] + zn * theCoeff[2];
324
325 Standard_Real d1 = theCoeff[0] * x + theCoeff[1] * y + theCoeff[2] * z - theCoeff[3];
326 Standard_Real d2 = d1 * d1;
327 Standard_Real d3 = d1 * d2 / 3.0;
328 Standard_Real dv = s * d1;
329
330 theOutInertia.Mass += dv;
331 theOutInertia.Ix += (x - (theCoeff[0] * d1 * 0.5)) * dv;
332 theOutInertia.Iy += (y - (theCoeff[1] * d1 * 0.5)) * dv;
333 theOutInertia.Iz += (z - (theCoeff[2] * d1 * 0.5)) * dv;
334
335 const Standard_Real px = x - theCoeff[0] * d1;
336 const Standard_Real py = y - theCoeff[1] * d1;
337 const Standard_Real pz = z - theCoeff[2] * d1;
338
339 x = px * px * d1 + px * theCoeff[0] * d2 + theCoeff[0] * theCoeff[0] * d3;
340 y = py * py * d1 + py * theCoeff[1] * d2 + theCoeff[1] * theCoeff[1] * d3;
341 z = pz * pz * d1 + pz * theCoeff[2] * d2 + theCoeff[2] * theCoeff[2] * d3;
342
343 theOutInertia.Ixx += (y + z) * s;
344 theOutInertia.Iyy += (x + z) * s;
345 theOutInertia.Izz += (x + y) * s;
346
347 d2 *= 0.5;
348 x = (py * pz * d1) + (py * theCoeff[2] * d2) + (pz * theCoeff[1] * d2) + (theCoeff[1] * theCoeff[2] * d3);
349 y = (px * pz * d1) + (pz * theCoeff[0] * d2) + (px * theCoeff[2] * d2) + (theCoeff[0] * theCoeff[2] * d3);
350 z = (px * py * d1) + (px * theCoeff[1] * d2) + (py * theCoeff[0] * d2) + (theCoeff[0] * theCoeff[1] * d3);
351
352 theOutInertia.Ixy -= z * s;
353 theOutInertia.Iyz -= x * s;
354 theOutInertia.Ixz -= y * s;
355 }
356 }
357
358 //=======================================================================
359 //function : computeSInertiaOfElementaryPart
360 //purpose :
361 //=======================================================================
computeSInertiaOfElementaryPart(const gp_Pnt & thePoint,const gp_Vec & theNormal,const gp_Pnt & theLocation,const Standard_Real theWeight,BRepGProp_Gauss::Inertia & theOutInertia)362 void BRepGProp_Gauss::computeSInertiaOfElementaryPart(
363 const gp_Pnt& thePoint,
364 const gp_Vec& theNormal,
365 const gp_Pnt& theLocation,
366 const Standard_Real theWeight,
367 BRepGProp_Gauss::Inertia& theOutInertia)
368 {
369 // ds - Jacobien (x, y, z) -> (u, v) = ||n||
370 const Standard_Real ds = mult(theNormal.Magnitude(), theWeight);
371 const Standard_Real x = add(thePoint.X(), -theLocation.X());
372 const Standard_Real y = add(thePoint.Y(), -theLocation.Y());
373 const Standard_Real z = add(thePoint.Z(), -theLocation.Z());
374
375 theOutInertia.Mass = add(theOutInertia.Mass, ds);
376
377 const Standard_Real XdS = mult(x, ds);
378 const Standard_Real YdS = mult(y, ds);
379 const Standard_Real ZdS = mult(z, ds);
380
381 theOutInertia.Ix = add(theOutInertia.Ix, XdS);
382 theOutInertia.Iy = add(theOutInertia.Iy, YdS);
383 theOutInertia.Iz = add(theOutInertia.Iz, ZdS);
384 theOutInertia.Ixy = add(theOutInertia.Ixy, mult(x, YdS));
385 theOutInertia.Iyz = add(theOutInertia.Iyz, mult(y, ZdS));
386 theOutInertia.Ixz = add(theOutInertia.Ixz, mult(x, ZdS));
387
388 const Standard_Real XXdS = mult(x, XdS);
389 const Standard_Real YYdS = mult(y, YdS);
390 const Standard_Real ZZdS = mult(z, ZdS);
391
392 theOutInertia.Ixx = add(theOutInertia.Ixx, add(YYdS, ZZdS));
393 theOutInertia.Iyy = add(theOutInertia.Iyy, add(XXdS, ZZdS));
394 theOutInertia.Izz = add(theOutInertia.Izz, add(XXdS, YYdS));
395 }
396
397 //=======================================================================
398 //function : checkBounds
399 //purpose :
400 //=======================================================================
checkBounds(const Standard_Real theU1,const Standard_Real theU2,const Standard_Real theV1,const Standard_Real theV2)401 void BRepGProp_Gauss::checkBounds(const Standard_Real theU1,
402 const Standard_Real theU2,
403 const Standard_Real theV1,
404 const Standard_Real theV2)
405 {
406 if (Precision::IsInfinite(theU1) || Precision::IsInfinite(theU2) ||
407 Precision::IsInfinite(theV1) || Precision::IsInfinite(theV2))
408 {
409 add = (::AddInf);
410 mult = (::MultInf);
411 }
412 }
413
414 //=======================================================================
415 //function : addAndRestoreInertia
416 //purpose :
417 //=======================================================================
addAndRestoreInertia(const BRepGProp_Gauss::Inertia & theInInertia,BRepGProp_Gauss::Inertia & theOutInertia)418 void BRepGProp_Gauss::addAndRestoreInertia(
419 const BRepGProp_Gauss::Inertia& theInInertia,
420 BRepGProp_Gauss::Inertia& theOutInertia)
421 {
422 theOutInertia.Mass = add(theOutInertia.Mass, theInInertia.Mass);
423 theOutInertia.Ix = add(theOutInertia.Ix, theInInertia.Ix);
424 theOutInertia.Iy = add(theOutInertia.Iy, theInInertia.Iy);
425 theOutInertia.Iz = add(theOutInertia.Iz, theInInertia.Iz);
426 theOutInertia.Ixx = add(theOutInertia.Ixx, theInInertia.Ixx);
427 theOutInertia.Iyy = add(theOutInertia.Iyy, theInInertia.Iyy);
428 theOutInertia.Izz = add(theOutInertia.Izz, theInInertia.Izz);
429 theOutInertia.Ixy = add(theOutInertia.Ixy, theInInertia.Ixy);
430 theOutInertia.Ixz = add(theOutInertia.Ixz, theInInertia.Ixz);
431 theOutInertia.Iyz = add(theOutInertia.Iyz, theInInertia.Iyz);
432 }
433
434 //=======================================================================
435 //function : multAndRestoreInertia
436 //purpose :
437 //=======================================================================
multAndRestoreInertia(const Standard_Real theValue,BRepGProp_Gauss::Inertia & theInOutInertia)438 void BRepGProp_Gauss::multAndRestoreInertia(
439 const Standard_Real theValue,
440 BRepGProp_Gauss::Inertia& theInOutInertia)
441 {
442 theInOutInertia.Mass = mult(theInOutInertia.Mass, theValue);
443 theInOutInertia.Ix = mult(theInOutInertia.Ix, theValue);
444 theInOutInertia.Iy = mult(theInOutInertia.Iy, theValue);
445 theInOutInertia.Iz = mult(theInOutInertia.Iz, theValue);
446 theInOutInertia.Ixx = mult(theInOutInertia.Ixx, theValue);
447 theInOutInertia.Iyy = mult(theInOutInertia.Iyy, theValue);
448 theInOutInertia.Izz = mult(theInOutInertia.Izz, theValue);
449 theInOutInertia.Ixy = mult(theInOutInertia.Ixy, theValue);
450 theInOutInertia.Ixz = mult(theInOutInertia.Ixz, theValue);
451 theInOutInertia.Iyz = mult(theInOutInertia.Iyz, theValue);
452 }
453
454 //=======================================================================
455 //function : convert
456 //purpose :
457 //=======================================================================
convert(const BRepGProp_Gauss::Inertia & theInertia,gp_Pnt & theOutGravityCenter,gp_Mat & theOutMatrixOfInertia,Standard_Real & theOutMass)458 void BRepGProp_Gauss::convert(const BRepGProp_Gauss::Inertia& theInertia,
459 gp_Pnt& theOutGravityCenter,
460 gp_Mat& theOutMatrixOfInertia,
461 Standard_Real& theOutMass)
462 {
463 if (Abs(theInertia.Mass) >= EPS_DIM)
464 {
465 const Standard_Real anInvMass = 1.0 / theInertia.Mass;
466 theOutGravityCenter.SetX(theInertia.Ix * anInvMass);
467 theOutGravityCenter.SetY(theInertia.Iy * anInvMass);
468 theOutGravityCenter.SetZ(theInertia.Iz * anInvMass);
469
470 theOutMass = theInertia.Mass;
471 }
472 else
473 {
474 theOutMass = 0.0;
475 theOutGravityCenter.SetCoord(0.0, 0.0, 0.0);
476 }
477
478 theOutMatrixOfInertia = gp_Mat(
479 gp_XYZ ( theInertia.Ixx, -theInertia.Ixy, -theInertia.Ixz),
480 gp_XYZ (-theInertia.Ixy, theInertia.Iyy, -theInertia.Iyz),
481 gp_XYZ (-theInertia.Ixz, -theInertia.Iyz, theInertia.Izz));
482 }
483
484 //=======================================================================
485 //function : convert
486 //purpose :
487 //=======================================================================
convert(const BRepGProp_Gauss::Inertia & theInertia,const Standard_Real theCoeff[],const Standard_Boolean theIsByPoint,gp_Pnt & theOutGravityCenter,gp_Mat & theOutMatrixOfInertia,Standard_Real & theOutMass)488 void BRepGProp_Gauss::convert(const BRepGProp_Gauss::Inertia& theInertia,
489 const Standard_Real theCoeff[],
490 const Standard_Boolean theIsByPoint,
491 gp_Pnt& theOutGravityCenter,
492 gp_Mat& theOutMatrixOfInertia,
493 Standard_Real& theOutMass)
494 {
495 convert(theInertia, theOutGravityCenter, theOutMatrixOfInertia, theOutMass);
496 if (Abs(theInertia.Mass) >= EPS_DIM && theIsByPoint)
497 {
498 const Standard_Real anInvMass = 1.0 / theInertia.Mass;
499 if (theIsByPoint == Standard_True)
500 {
501 theOutGravityCenter.SetX(theCoeff[0] + theInertia.Ix * anInvMass);
502 theOutGravityCenter.SetY(theCoeff[1] + theInertia.Iy * anInvMass);
503 theOutGravityCenter.SetZ(theCoeff[2] + theInertia.Iz * anInvMass);
504 }
505 else
506 {
507 theOutGravityCenter.SetX(theInertia.Ix * anInvMass);
508 theOutGravityCenter.SetY(theInertia.Iy * anInvMass);
509 theOutGravityCenter.SetZ(theInertia.Iz * anInvMass);
510 }
511
512 theOutMass = theInertia.Mass;
513 }
514 else
515 {
516 theOutMass = 0.0;
517 theOutGravityCenter.SetCoord(0.0, 0.0, 0.0);
518 }
519
520 theOutMatrixOfInertia = gp_Mat(
521 gp_XYZ (theInertia.Ixx, theInertia.Ixy, theInertia.Ixz),
522 gp_XYZ (theInertia.Ixy, theInertia.Iyy, theInertia.Iyz),
523 gp_XYZ (theInertia.Ixz, theInertia.Iyz, theInertia.Izz));
524 }
525
526 //=======================================================================
527 //function : Compute
528 //purpose :
529 //=======================================================================
Compute(BRepGProp_Face & theSurface,BRepGProp_Domain & theDomain,const gp_Pnt & theLocation,const Standard_Real theEps,const Standard_Real theCoeff[],const Standard_Boolean theIsByPoint,Standard_Real & theOutMass,gp_Pnt & theOutGravityCenter,gp_Mat & theOutInertia)530 Standard_Real BRepGProp_Gauss::Compute(
531 BRepGProp_Face& theSurface,
532 BRepGProp_Domain& theDomain,
533 const gp_Pnt& theLocation,
534 const Standard_Real theEps,
535 const Standard_Real theCoeff[],
536 const Standard_Boolean theIsByPoint,
537 Standard_Real& theOutMass,
538 gp_Pnt& theOutGravityCenter,
539 gp_Mat& theOutInertia)
540 {
541 const Standard_Boolean isErrorCalculation =
542 ( 0.0 > theEps || theEps < 0.001 ) ? Standard_True : Standard_False;
543 const Standard_Boolean isVerifyComputation =
544 ( 0.0 < theEps && theEps < 0.001 ) ? Standard_True : Standard_False;
545
546 Standard_Real anEpsilon= Abs(theEps);
547
548 BRepGProp_Gauss::Inertia anInertia;
549 InertiaArray anInertiaL = new NCollection_Array1<Inertia>(1, SM);
550 InertiaArray anInertiaU = new NCollection_Array1<Inertia>(1, SM);
551
552 // Prepare Gauss points and weights
553 NCollection_Handle<math_Vector> LGaussP[2];
554 NCollection_Handle<math_Vector> LGaussW[2];
555 NCollection_Handle<math_Vector> UGaussP[2];
556 NCollection_Handle<math_Vector> UGaussW[2];
557
558 const Standard_Integer aNbGaussPoint =
559 RealToInt(Ceiling(ERROR_ALGEBR_RATIO * GPM));
560
561 LGaussP[0] = new math_Vector(1, GPM);
562 LGaussP[1] = new math_Vector(1, aNbGaussPoint);
563 LGaussW[0] = new math_Vector(1, GPM);
564 LGaussW[1] = new math_Vector(1, aNbGaussPoint);
565
566 UGaussP[0] = new math_Vector(1, GPM);
567 UGaussP[1] = new math_Vector(1, aNbGaussPoint);
568 UGaussW[0] = new math_Vector(1, GPM);
569 UGaussW[1] = new math_Vector(1, aNbGaussPoint);
570
571 NCollection_Handle<math_Vector> L1 = new math_Vector(1, SM);
572 NCollection_Handle<math_Vector> L2 = new math_Vector(1, SM);
573 NCollection_Handle<math_Vector> U1 = new math_Vector(1, SM);
574 NCollection_Handle<math_Vector> U2 = new math_Vector(1, SM);
575
576 NCollection_Handle<math_Vector> ErrL = new math_Vector(1, SM, 0.0);
577 NCollection_Handle<math_Vector> ErrU = new math_Vector(1, SM, 0.0);
578 NCollection_Handle<math_Vector> ErrUL = new math_Vector(1, SM, 0.0);
579
580 // Face parametrization in U and V direction
581 Standard_Real BV1, BV2, BU1, BU2;
582 theSurface.Bounds(BU1, BU2, BV1, BV2);
583 checkBounds(BU1, BU2, BV1, BV2);
584
585 //
586 const Standard_Integer NumSubs = SUBS_POWER;
587 const TopoDS_Face& aF = theSurface.GetFace();
588 const Standard_Boolean isNaturalRestriction = (aF.NbChildren () == 0); //theSurface.NaturalRestriction();
589
590 Standard_Real CIx, CIy, CIz, CIxy, CIxz, CIyz;
591 Standard_Real CDim[2], CIxx[2], CIyy[2], CIzz[2];
592
593 // Boundary curve parametrization
594 Standard_Real u1 = BU1, u2, l1, l2, lm, lr, l, v;
595
596 // On the boundary curve u-v
597 gp_Pnt2d Puv;
598 gp_Vec2d Vuv;
599 Standard_Real Dul; // Dul = Du / Dl
600
601 Standard_Integer iLS, iLSubEnd, iGL, iGLEnd, NbLGaussP[2], LRange[2], iL, kL, kLEnd, IL, JL;
602 Standard_Integer i, iUSubEnd, NbUGaussP[2], URange[2], kU, kUEnd, IU, JU;
603 Standard_Integer UMaxSubs, LMaxSubs;
604
605 Standard_Real ErrorU, ErrorL, ErrorLMax = 0.0, Eps = 0.0, EpsL = 0.0, EpsU = 0.0;
606 iGLEnd = isErrorCalculation ? 2 : 1;
607
608 NbUGaussP[0] = theSurface.SIntOrder(anEpsilon);
609 NbUGaussP[1] = RealToInt( Ceiling(ERROR_ALGEBR_RATIO * NbUGaussP[0]) );
610
611 math::GaussPoints (NbUGaussP[0], *UGaussP[0]);
612 math::GaussWeights(NbUGaussP[0], *UGaussW[0]);
613 math::GaussPoints (NbUGaussP[1], *UGaussP[1]);
614 math::GaussWeights(NbUGaussP[1], *UGaussW[1]);
615
616 const Standard_Integer aNbUSubs = theSurface.SUIntSubs();
617 TColStd_Array1OfReal UKnots(1, aNbUSubs + 1);
618 theSurface.UKnots(UKnots);
619
620 while (isNaturalRestriction || theDomain.More())
621 {
622 if (isNaturalRestriction)
623 {
624 NbLGaussP[0] = Min(2 * NbUGaussP[0], math::GaussPointsMax());
625 }
626 else
627 {
628 if (!theSurface.Load(theDomain.Value()))
629 {
630 return Precision::Infinite();
631 }
632 NbLGaussP[0] = theSurface.LIntOrder(anEpsilon);
633 }
634
635 NbLGaussP[1] = RealToInt( Ceiling(ERROR_ALGEBR_RATIO * NbLGaussP[0]) );
636
637 math::GaussPoints (NbLGaussP[0], *LGaussP[0]);
638 math::GaussWeights(NbLGaussP[0], *LGaussW[0]);
639 math::GaussPoints (NbLGaussP[1], *LGaussP[1]);
640 math::GaussWeights(NbLGaussP[1], *LGaussW[1]);
641
642 const Standard_Integer aNbLSubs =
643 isNaturalRestriction ? theSurface.SVIntSubs(): theSurface.LIntSubs();
644 TColStd_Array1OfReal LKnots(1, aNbLSubs + 1);
645
646 if (isNaturalRestriction)
647 {
648 theSurface.VKnots(LKnots);
649 l1 = BV1;
650 l2 = BV2;
651 }
652 else
653 {
654 theSurface.LKnots(LKnots);
655 l1 = theSurface.FirstParameter();
656 l2 = theSurface.LastParameter();
657 }
658 ErrorL = 0.0;
659 kLEnd = 1; JL = 0;
660
661 if (Abs(l2 - l1) > EPS_PARAM)
662 {
663 iLSubEnd = FillIntervalBounds(l1, l2, LKnots, NumSubs, anInertiaL, L1, L2, ErrL, ErrUL);
664 LMaxSubs = BRepGProp_Gauss::MaxSubs(iLSubEnd);
665
666 if (LMaxSubs > SM)
667 {
668 LMaxSubs = SM;
669 }
670
671 BRepGProp_Gauss::InitMass(0.0, 1, LMaxSubs, anInertiaL);
672 BRepGProp_Gauss::Init(ErrL, 0.0, 1, LMaxSubs);
673 BRepGProp_Gauss::Init(ErrUL, 0.0, 1, LMaxSubs);
674
675 do // while: L
676 {
677 if (++JL > iLSubEnd)
678 {
679 LRange[0] = IL = ErrL->Max();
680 LRange[1] = JL;
681 L1->Value(JL) = (L1->Value(IL) + L2->Value(IL)) * 0.5;
682 L2->Value(JL) = L2->Value(IL);
683 L2->Value(IL) = L1->Value(JL);
684 }
685 else
686 {
687 LRange[0] = IL = JL;
688 }
689
690 if (JL == LMaxSubs || Abs(L2->Value(JL) - L1->Value(JL)) < EPS_PARAM)
691 {
692 if (kLEnd == 1)
693 {
694 anInertiaL->ChangeValue(JL).Reset();
695 ErrL->Value(JL) = 0.0;
696 }
697 else
698 {
699 --JL;
700 EpsL = ErrorL;
701 Eps = EpsL / 0.9;
702 break;
703 }
704 }
705 else
706 {
707 for (kL = 0; kL < kLEnd; kL++)
708 {
709 iLS = LRange[kL];
710 lm = 0.5 * (L2->Value(iLS) + L1->Value(iLS));
711 lr = 0.5 * (L2->Value(iLS) - L1->Value(iLS));
712
713 CIx = CIy = CIz = CIxy = CIxz = CIyz = 0.0;
714
715 for (iGL = 0; iGL < iGLEnd; ++iGL)
716 {
717 CDim[iGL] = CIxx[iGL] = CIyy[iGL] = CIzz[iGL] = 0.0;
718
719 for (iL = 1; iL <= NbLGaussP[iGL]; iL++)
720 {
721 l = lm + lr * LGaussP[iGL]->Value(iL);
722 if (isNaturalRestriction)
723 {
724 v = l;
725 u2 = BU2;
726 Dul = LGaussW[iGL]->Value(iL);
727 }
728 else
729 {
730 theSurface.D12d (l, Puv, Vuv);
731 Dul = Vuv.Y() * LGaussW[iGL]->Value(iL); // Dul = Du / Dl
732
733 if (Abs(Dul) < EPS_PARAM)
734 continue;
735
736 v = Puv.Y();
737 u2 = Puv.X();
738
739 // Check on cause out off bounds of value current parameter
740 if (v < BV1)
741 v = BV1;
742 else if (v > BV2)
743 v = BV2;
744
745 if (u2 < BU1)
746 u2 = BU1;
747 else if (u2 > BU2)
748 u2 = BU2;
749 }
750
751 ErrUL->Value(iLS) = 0.0;
752 kUEnd = 1;
753 JU = 0;
754
755 if (Abs(u2 - u1) < EPS_PARAM)
756 continue;
757
758 NCollection_Handle<math_Vector> aDummy;
759 iUSubEnd = FillIntervalBounds(u1, u2, UKnots, NumSubs, anInertiaU, U1, U2, ErrU, aDummy);
760 UMaxSubs = BRepGProp_Gauss::MaxSubs(iUSubEnd);
761
762 if (UMaxSubs > SM)
763 UMaxSubs = SM;
764
765 BRepGProp_Gauss::InitMass(0.0, 1, UMaxSubs, anInertiaU);
766 BRepGProp_Gauss::Init(ErrU, 0.0, 1, UMaxSubs);
767 ErrorU = 0.0;
768
769 do
770 {//while: U
771 if (++JU > iUSubEnd)
772 {
773 URange[0] = IU = ErrU->Max();
774 URange[1] = JU;
775
776 U1->Value(JU) = (U1->Value(IU) + U2->Value(IU)) * 0.5;
777 U2->Value(JU) = U2->Value(IU);
778 U2->Value(IU) = U1->Value(JU);
779 }
780 else
781 URange[0] = IU = JU;
782
783 if (JU == UMaxSubs || Abs(U2->Value(JU) - U1->Value(JU)) < EPS_PARAM)
784 if (kUEnd == 1)
785 {
786 ErrU->Value(JU) = 0.0;
787 anInertiaU->ChangeValue(JU).Reset();
788 }
789 else
790 {
791 --JU;
792 EpsU = ErrorU;
793 Eps = 10. * EpsU * Abs((u2 - u1) * Dul);
794 EpsL = 0.9 * Eps;
795 break;
796 }
797 else
798 {
799 gp_Pnt aPoint;
800 gp_Vec aNormal;
801
802 for (kU = 0; kU < kUEnd; ++kU)
803 {
804 BRepGProp_Gauss::Inertia aLocal[2];
805
806 Standard_Integer iUS = URange[kU];
807 const Standard_Integer aLength = iGLEnd - iGL;
808
809 const Standard_Real um = 0.5 * (U2->Value(iUS) + U1->Value(iUS));
810 const Standard_Real ur = 0.5 * (U2->Value(iUS) - U1->Value(iUS));
811
812 for (Standard_Integer iGU = 0; iGU < aLength; ++iGU)
813 {
814 for (Standard_Integer iU = 1; iU <= NbUGaussP[iGU]; ++iU)
815 {
816 Standard_Real w = UGaussW[iGU]->Value(iU);
817 const Standard_Real u = um + ur * UGaussP[iGU]->Value(iU);
818
819 theSurface.Normal(u, v, aPoint, aNormal);
820
821 if (myType == Vinert)
822 {
823 computeVInertiaOfElementaryPart(
824 aPoint, aNormal, theLocation, w, theCoeff, theIsByPoint, aLocal[iGU]);
825 }
826 else
827 {
828 if (iGU > 0)
829 aLocal[iGU].Mass += (w * aNormal.Magnitude());
830 else
831 {
832 computeSInertiaOfElementaryPart(
833 aPoint, aNormal, theLocation, w, aLocal[iGU]);
834 }
835 }
836 }
837 }
838
839 BRepGProp_Gauss::Inertia& anUI =
840 anInertiaU->ChangeValue(iUS);
841
842 anUI.Mass = mult(aLocal[0].Mass, ur);
843
844 if (myType == Vinert)
845 {
846 anUI.Ixx = mult(aLocal[0].Ixx, ur);
847 anUI.Iyy = mult(aLocal[0].Iyy, ur);
848 anUI.Izz = mult(aLocal[0].Izz, ur);
849 }
850
851 if (iGL > 0)
852 continue;
853
854 Standard_Real aDMass = Abs(aLocal[1].Mass - aLocal[0].Mass);
855
856 if (myType == Vinert)
857 {
858 aLocal[1].Ixx = Abs(aLocal[1].Ixx - aLocal[0].Ixx);
859 aLocal[1].Iyy = Abs(aLocal[1].Iyy - aLocal[0].Iyy);
860 aLocal[1].Izz = Abs(aLocal[1].Izz - aLocal[0].Izz);
861
862 anUI.Ix = mult(aLocal[0].Ix, ur);
863 anUI.Iy = mult(aLocal[0].Iy, ur);
864 anUI.Iz = mult(aLocal[0].Iz, ur);
865
866 anUI.Ixy = mult(aLocal[0].Ixy, ur);
867 anUI.Ixz = mult(aLocal[0].Ixz, ur);
868 anUI.Iyz = mult(aLocal[0].Iyz, ur);
869
870 #ifndef IS_MIN_DIM
871 aDMass = aLocal[1].Ixx + aLocal[1].Iyy + aLocal[1].Izz;
872 #endif
873
874 ErrU->Value(iUS) = mult(aDMass, ur);
875 }
876 else
877 {
878 anUI.Ix = mult(aLocal[0].Ix, ur);
879 anUI.Iy = mult(aLocal[0].Iy, ur);
880 anUI.Iz = mult(aLocal[0].Iz, ur);
881 anUI.Ixx = mult(aLocal[0].Ixx, ur);
882 anUI.Iyy = mult(aLocal[0].Iyy, ur);
883 anUI.Izz = mult(aLocal[0].Izz, ur);
884 anUI.Ixy = mult(aLocal[0].Ixy, ur);
885 anUI.Ixz = mult(aLocal[0].Ixz, ur);
886 anUI.Iyz = mult(aLocal[0].Iyz, ur);
887
888 ErrU->Value(iUS) = mult(aDMass, ur);
889 }
890 }
891 }
892
893 if (JU == iUSubEnd)
894 {
895 kUEnd = 2;
896 ErrorU = ErrU->Value(ErrU->Max());
897 }
898 } while ( (ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1 );
899
900 for (i = 1; i <= JU; ++i)
901 {
902 const BRepGProp_Gauss::Inertia& anIU =
903 anInertiaU->Value(i);
904
905 CDim[iGL] = add(CDim[iGL], mult(anIU.Mass, Dul));
906 CIxx[iGL] = add(CIxx[iGL], mult(anIU.Ixx, Dul));
907 CIyy[iGL] = add(CIyy[iGL], mult(anIU.Iyy, Dul));
908 CIzz[iGL] = add(CIzz[iGL], mult(anIU.Izz, Dul));
909 }
910
911 if (iGL > 0)
912 continue;
913
914 ErrUL->Value(iLS) = ErrorU * Abs((u2 - u1) * Dul);
915
916 for (i = 1; i <= JU; ++i)
917 {
918 const BRepGProp_Gauss::Inertia& anIU =
919 anInertiaU->Value(i);
920
921 CIx = add(CIx, mult(anIU.Ix, Dul));
922 CIy = add(CIy, mult(anIU.Iy, Dul));
923 CIz = add(CIz, mult(anIU.Iz, Dul));
924
925 CIxy = add(CIxy, mult(anIU.Ixy, Dul));
926 CIxz = add(CIxz, mult(anIU.Ixz, Dul));
927 CIyz = add(CIyz, mult(anIU.Iyz, Dul));
928 }
929 }//for: iL
930 }//for: iGL
931
932 BRepGProp_Gauss::Inertia& aLI = anInertiaL->ChangeValue(iLS);
933
934 aLI.Mass = mult(CDim[0], lr);
935 aLI.Ixx = mult(CIxx[0], lr);
936 aLI.Iyy = mult(CIyy[0], lr);
937 aLI.Izz = mult(CIzz[0], lr);
938
939 if (iGLEnd == 2)
940 {
941 Standard_Real aSubDim = Abs(CDim[1] - CDim[0]);
942
943 if (myType == Vinert)
944 {
945 ErrorU = ErrUL->Value(iLS);
946
947 CIxx[1] = Abs(CIxx[1] - CIxx[0]);
948 CIyy[1] = Abs(CIyy[1] - CIyy[0]);
949 CIzz[1] = Abs(CIzz[1] - CIzz[0]);
950
951 #ifndef IS_MIN_DIM
952 aSubDim = CIxx[1] + CIyy[1] + CIzz[1];
953 #endif
954
955 ErrL->Value(iLS) = add(mult(aSubDim, lr), ErrorU);
956 }
957 else
958 {
959 ErrL->Value(iLS) = add(mult(aSubDim, lr), ErrUL->Value(iLS));
960 }
961 }
962
963 aLI.Ix = mult(CIx, lr);
964 aLI.Iy = mult(CIy, lr);
965 aLI.Iz = mult(CIz, lr);
966
967 aLI.Ixy = mult(CIxy, lr);
968 aLI.Ixz = mult(CIxz, lr);
969 aLI.Iyz = mult(CIyz, lr);
970 }//for: (kL)iLS
971 }
972
973 // Calculate/correct epsilon of computation by current value of dim
974 // That is need for not spend time for
975 if (JL == iLSubEnd)
976 {
977 kLEnd = 2;
978
979 Standard_Real DDim = 0.0;
980 for (i = 1; i <= JL; ++i)
981 DDim += anInertiaL->Value(i).Mass;
982
983 #ifndef IS_MIN_DIM
984 {
985 if (myType == Vinert)
986 {
987 Standard_Real DIxx = 0.0, DIyy = 0.0, DIzz = 0.0;
988 for (i = 1; i <= JL; ++i)
989 {
990 const BRepGProp_Gauss::Inertia& aLocalL =
991 anInertiaL->Value(i);
992
993 DIxx += aLocalL.Ixx;
994 DIyy += aLocalL.Iyy;
995 DIzz += aLocalL.Izz;
996 }
997
998 DDim = Abs(DIxx) + Abs(DIyy) + Abs(DIzz);
999 }
1000 }
1001 #endif
1002
1003 DDim = Abs(DDim * anEpsilon);
1004
1005 if (DDim > Eps)
1006 {
1007 Eps = DDim;
1008 EpsL = 0.9 * Eps;
1009 }
1010 }
1011 if (kLEnd == 2)
1012 {
1013 ErrorL = ErrL->Value(ErrL->Max());
1014 }
1015 } while ( (ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1 );
1016
1017 for ( i = 1; i <= JL; i++ )
1018 {
1019 addAndRestoreInertia(anInertiaL->Value(i), anInertia);
1020 }
1021
1022 ErrorLMax = Max(ErrorLMax, ErrorL);
1023 }
1024
1025 if (isNaturalRestriction)
1026 break;
1027
1028 theDomain.Next();
1029 }
1030
1031 if (myType == Vinert)
1032 convert(anInertia, theCoeff, theIsByPoint, theOutGravityCenter, theOutInertia, theOutMass);
1033 else
1034 convert(anInertia, theOutGravityCenter, theOutInertia, theOutMass);
1035
1036 if (iGLEnd == 2)
1037 {
1038 if (theOutMass != 0.0)
1039 {
1040 Eps = ErrorLMax / Abs(theOutMass);
1041
1042 #ifndef IS_MIN_DIM
1043 {
1044 if (myType == Vinert)
1045 Eps = ErrorLMax / (Abs(anInertia.Ixx) +
1046 Abs(anInertia.Iyy) +
1047 Abs(anInertia.Izz));
1048 }
1049 #endif
1050 }
1051 else
1052 {
1053 Eps = 0.0;
1054 }
1055 }
1056 else
1057 {
1058 Eps = anEpsilon;
1059 }
1060
1061 return Eps;
1062 }
1063
1064 //=======================================================================
1065 //function : Compute
1066 //purpose :
1067 //=======================================================================
Compute(BRepGProp_Face & theSurface,BRepGProp_Domain & theDomain,const gp_Pnt & theLocation,const Standard_Real theEps,Standard_Real & theOutMass,gp_Pnt & theOutGravityCenter,gp_Mat & theOutInertia)1068 Standard_Real BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface,
1069 BRepGProp_Domain& theDomain,
1070 const gp_Pnt& theLocation,
1071 const Standard_Real theEps,
1072 Standard_Real& theOutMass,
1073 gp_Pnt& theOutGravityCenter,
1074 gp_Mat& theOutInertia)
1075 {
1076 Standard_ASSERT_RAISE(myType == Sinert, "BRepGProp_Gauss: Incorrect type");
1077
1078 return Compute(theSurface,
1079 theDomain,
1080 theLocation,
1081 theEps,
1082 NULL,
1083 Standard_True,
1084 theOutMass,
1085 theOutGravityCenter,
1086 theOutInertia);
1087 }
1088
1089 //=======================================================================
1090 //function : Compute
1091 //purpose :
1092 //=======================================================================
Compute(BRepGProp_Face & theSurface,BRepGProp_Domain & theDomain,const gp_Pnt & theLocation,Standard_Real & theOutMass,gp_Pnt & theOutGravityCenter,gp_Mat & theOutInertia)1093 void BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface,
1094 BRepGProp_Domain& theDomain,
1095 const gp_Pnt& theLocation,
1096 Standard_Real& theOutMass,
1097 gp_Pnt& theOutGravityCenter,
1098 gp_Mat& theOutInertia)
1099 {
1100 Standard_ASSERT_RAISE(myType == Sinert, "BRepGProp_Gauss: Incorrect type");
1101
1102 Standard_Real u1, u2, v1, v2;
1103 theSurface.Bounds (u1, u2, v1, v2);
1104 checkBounds(u1, u2, v1, v2);
1105
1106 const Standard_Integer NbUGaussgp_Pnts =
1107 Min(theSurface.UIntegrationOrder(), math::GaussPointsMax());
1108
1109 const Standard_Integer NbVGaussgp_Pnts =
1110 Min(theSurface.VIntegrationOrder(), math::GaussPointsMax());
1111
1112 const Standard_Integer NbGaussgp_Pnts =
1113 Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts);
1114
1115 // Number of Gauss points for the integration on the face
1116 math_Vector GaussSPV (1, NbGaussgp_Pnts);
1117 math_Vector GaussSWV (1, NbGaussgp_Pnts);
1118 math::GaussPoints (NbGaussgp_Pnts, GaussSPV);
1119 math::GaussWeights(NbGaussgp_Pnts, GaussSWV);
1120
1121 BRepGProp_Gauss::Inertia anInertia;
1122 for (; theDomain.More(); theDomain.Next())
1123 {
1124 if (!theSurface.Load(theDomain.Value()))
1125 {
1126 return;
1127 }
1128
1129 Standard_Integer NbCGaussgp_Pnts =
1130 Min(theSurface.IntegrationOrder(), math::GaussPointsMax());
1131
1132 NbCGaussgp_Pnts = Max(NbCGaussgp_Pnts, NbGaussgp_Pnts);
1133
1134 math_Vector GaussCP(1, NbCGaussgp_Pnts);
1135 math_Vector GaussCW(1, NbCGaussgp_Pnts);
1136 math::GaussPoints (NbCGaussgp_Pnts, GaussCP);
1137 math::GaussWeights(NbCGaussgp_Pnts, GaussCW);
1138
1139
1140 const Standard_Real l1 = theSurface.FirstParameter();
1141 const Standard_Real l2 = theSurface.LastParameter ();
1142 const Standard_Real lm = 0.5 * (l2 + l1);
1143 const Standard_Real lr = 0.5 * (l2 - l1);
1144
1145 BRepGProp_Gauss::Inertia aCInertia;
1146 for (Standard_Integer i = 1; i <= NbCGaussgp_Pnts; ++i)
1147 {
1148 const Standard_Real l = lm + lr * GaussCP(i);
1149
1150 gp_Pnt2d Puv;
1151 gp_Vec2d Vuv;
1152 theSurface.D12d(l, Puv, Vuv);
1153
1154 const Standard_Real v = Puv.Y();
1155 u2 = Puv.X();
1156
1157 const Standard_Real Dul = Vuv.Y() * GaussCW(i);
1158 const Standard_Real um = 0.5 * (u2 + u1);
1159 const Standard_Real ur = 0.5 * (u2 - u1);
1160
1161 BRepGProp_Gauss::Inertia aLocalInertia;
1162 for (Standard_Integer j = 1; j <= NbGaussgp_Pnts; ++j)
1163 {
1164 const Standard_Real u = add(um, mult(ur, GaussSPV(j)));
1165 const Standard_Real aWeight = Dul * GaussSWV(j);
1166
1167 gp_Pnt aPoint;
1168 gp_Vec aNormal;
1169 theSurface.Normal (u, v, aPoint, aNormal);
1170
1171 computeSInertiaOfElementaryPart(aPoint, aNormal, theLocation, aWeight, aLocalInertia);
1172 }
1173
1174 multAndRestoreInertia(ur, aLocalInertia);
1175 addAndRestoreInertia (aLocalInertia, aCInertia);
1176 }
1177
1178 multAndRestoreInertia(lr, aCInertia);
1179 addAndRestoreInertia (aCInertia, anInertia);
1180 }
1181
1182 convert(anInertia, theOutGravityCenter, theOutInertia, theOutMass);
1183 }
1184
1185 //=======================================================================
1186 //function : Compute
1187 //purpose :
1188 //=======================================================================
Compute(BRepGProp_Face & theSurface,BRepGProp_Domain & theDomain,const gp_Pnt & theLocation,const Standard_Real theCoeff[],const Standard_Boolean theIsByPoint,Standard_Real & theOutMass,gp_Pnt & theOutGravityCenter,gp_Mat & theOutInertia)1189 void BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface,
1190 BRepGProp_Domain& theDomain,
1191 const gp_Pnt& theLocation,
1192 const Standard_Real theCoeff[],
1193 const Standard_Boolean theIsByPoint,
1194 Standard_Real& theOutMass,
1195 gp_Pnt& theOutGravityCenter,
1196 gp_Mat& theOutInertia)
1197 {
1198 Standard_ASSERT_RAISE(myType == Vinert, "BRepGProp_Gauss: Incorrect type");
1199
1200 Standard_Real u1, v1, u2, v2;
1201 theSurface.Bounds (u1, u2, v1, v2);
1202 checkBounds(u1, u2, v1, v2);
1203
1204 Standard_Real _u2 = u2; //OCC104
1205
1206 BRepGProp_Gauss::Inertia anInertia;
1207 for (; theDomain.More(); theDomain.Next())
1208 {
1209 if (!theSurface.Load(theDomain.Value()))
1210 {
1211 return;
1212 }
1213
1214 const Standard_Integer aVNbCGaussgp_Pnts =
1215 theSurface.VIntegrationOrder();
1216
1217 const Standard_Integer aNbGaussgp_Pnts =
1218 Min( Max(theSurface.IntegrationOrder(), aVNbCGaussgp_Pnts), math::GaussPointsMax() );
1219
1220 math_Vector GaussP(1, aNbGaussgp_Pnts);
1221 math_Vector GaussW(1, aNbGaussgp_Pnts);
1222 math::GaussPoints (aNbGaussgp_Pnts, GaussP);
1223 math::GaussWeights(aNbGaussgp_Pnts, GaussW);
1224
1225 const Standard_Real l1 = theSurface.FirstParameter();
1226 const Standard_Real l2 = theSurface.LastParameter();
1227 const Standard_Real lm = 0.5 * (l2 + l1);
1228 const Standard_Real lr = 0.5 * (l2 - l1);
1229
1230 BRepGProp_Gauss::Inertia aCInertia;
1231 for (Standard_Integer i = 1; i <= aNbGaussgp_Pnts; ++i)
1232 {
1233 const Standard_Real l = lm + lr * GaussP(i);
1234
1235 gp_Pnt2d Puv;
1236 gp_Vec2d Vuv;
1237
1238 theSurface.D12d(l, Puv, Vuv);
1239
1240 u2 = Puv.X();
1241 u2 = Min( Max(u1, u2), _u2 ); // OCC104
1242 const Standard_Real v = Min(Max(Puv.Y(), v1), v2);
1243
1244 const Standard_Real Dul = Vuv.Y() * GaussW(i);
1245 const Standard_Real um = 0.5 * (u2 + u1);
1246 const Standard_Real ur = 0.5 * (u2 - u1);
1247
1248 BRepGProp_Gauss::Inertia aLocalInertia;
1249 for (Standard_Integer j = 1; j <= aNbGaussgp_Pnts; ++j)
1250 {
1251 const Standard_Real u = um + ur * GaussP(j);
1252 const Standard_Real aWeight = Dul * GaussW(j);
1253
1254 gp_Pnt aPoint;
1255 gp_Vec aNormal;
1256
1257 theSurface.Normal(u, v, aPoint, aNormal);
1258
1259 computeVInertiaOfElementaryPart(
1260 aPoint,
1261 aNormal,
1262 theLocation,
1263 aWeight,
1264 theCoeff,
1265 theIsByPoint,
1266 aLocalInertia);
1267 }
1268
1269 multAndRestoreInertia(ur, aLocalInertia);
1270 addAndRestoreInertia (aLocalInertia, aCInertia);
1271 }
1272
1273 multAndRestoreInertia(lr, aCInertia);
1274 addAndRestoreInertia (aCInertia, anInertia);
1275 }
1276
1277 convert(anInertia, theCoeff, theIsByPoint, theOutGravityCenter, theOutInertia, theOutMass);
1278 }
1279
1280 //=======================================================================
1281 //function : Compute
1282 //purpose :
1283 //=======================================================================
Compute(const BRepGProp_Face & theSurface,const gp_Pnt & theLocation,const Standard_Real theCoeff[],const Standard_Boolean theIsByPoint,Standard_Real & theOutMass,gp_Pnt & theOutGravityCenter,gp_Mat & theOutInertia)1284 void BRepGProp_Gauss::Compute(const BRepGProp_Face& theSurface,
1285 const gp_Pnt& theLocation,
1286 const Standard_Real theCoeff[],
1287 const Standard_Boolean theIsByPoint,
1288 Standard_Real& theOutMass,
1289 gp_Pnt& theOutGravityCenter,
1290 gp_Mat& theOutInertia)
1291 {
1292 Standard_Real LowerU, UpperU, LowerV, UpperV;
1293 theSurface.Bounds(LowerU, UpperU, LowerV, UpperV);
1294 checkBounds(LowerU, UpperU, LowerV, UpperV);
1295
1296 const Standard_Integer UOrder =
1297 Min(theSurface.UIntegrationOrder(), math::GaussPointsMax());
1298 const Standard_Integer VOrder =
1299 Min(theSurface.VIntegrationOrder(), math::GaussPointsMax());
1300
1301 // Gauss points and weights
1302 math_Vector GaussPU(1, UOrder);
1303 math_Vector GaussWU(1, UOrder);
1304 math_Vector GaussPV(1, VOrder);
1305 math_Vector GaussWV(1, VOrder);
1306
1307 math::GaussPoints (UOrder, GaussPU);
1308 math::GaussWeights(UOrder, GaussWU);
1309 math::GaussPoints (VOrder, GaussPV);
1310 math::GaussWeights(VOrder, GaussWV);
1311
1312 const Standard_Real um = 0.5 * add(UpperU, LowerU);
1313 const Standard_Real vm = 0.5 * add(UpperV, LowerV);
1314 Standard_Real ur = 0.5 * add(UpperU, -LowerU);
1315 Standard_Real vr = 0.5 * add(UpperV, -LowerV);
1316
1317 gp_Pnt aPoint;
1318 gp_Vec aNormal;
1319
1320 BRepGProp_Gauss::Inertia anInertia;
1321 for (Standard_Integer j = 1; j <= VOrder; ++j)
1322 {
1323 BRepGProp_Gauss::Inertia anInertiaOfElementaryPart;
1324 const Standard_Real v = add(vm, mult(vr, GaussPV(j)));
1325
1326 for (Standard_Integer i = 1; i <= UOrder; ++i)
1327 {
1328 const Standard_Real aWeight = GaussWU(i);
1329 const Standard_Real u = add(um, mult(ur, GaussPU (i)));
1330 theSurface.Normal(u, v, aPoint, aNormal);
1331
1332 if (myType == Vinert)
1333 {
1334 computeVInertiaOfElementaryPart(
1335 aPoint,
1336 aNormal,
1337 theLocation,
1338 aWeight,
1339 theCoeff,
1340 theIsByPoint,
1341 anInertiaOfElementaryPart);
1342 }
1343 else // Sinert
1344 {
1345 computeSInertiaOfElementaryPart(
1346 aPoint,
1347 aNormal,
1348 theLocation,
1349 aWeight,
1350 anInertiaOfElementaryPart);
1351 }
1352 }
1353
1354 multAndRestoreInertia(GaussWV(j), anInertiaOfElementaryPart);
1355 addAndRestoreInertia (anInertiaOfElementaryPart, anInertia);
1356 }
1357 vr = mult(vr, ur);
1358 anInertia.Ixx = mult(vr, anInertia.Ixx);
1359 anInertia.Iyy = mult(vr, anInertia.Iyy);
1360 anInertia.Izz = mult(vr, anInertia.Izz);
1361 anInertia.Ixy = mult(vr, anInertia.Ixy);
1362 anInertia.Ixz = mult(vr, anInertia.Ixz);
1363 anInertia.Iyz = mult(vr, anInertia.Iyz);
1364
1365 if (myType == Vinert)
1366 {
1367 convert(anInertia, theCoeff, theIsByPoint, theOutGravityCenter, theOutInertia, theOutMass);
1368 }
1369 else // Sinert
1370 {
1371 convert(anInertia, theOutGravityCenter, theOutInertia, theOutMass);
1372 }
1373
1374 theOutMass *= vr;
1375 }
1376
1377 //=======================================================================
1378 //function : Compute
1379 //purpose :
1380 //=======================================================================
Compute(const BRepGProp_Face & theSurface,const gp_Pnt & theLocation,Standard_Real & theOutMass,gp_Pnt & theOutGravityCenter,gp_Mat & theOutInertia)1381 void BRepGProp_Gauss::Compute(const BRepGProp_Face& theSurface,
1382 const gp_Pnt& theLocation,
1383 Standard_Real& theOutMass,
1384 gp_Pnt& theOutGravityCenter,
1385 gp_Mat& theOutInertia)
1386 {
1387 Standard_ASSERT_RAISE(myType == Sinert, "BRepGProp_Gauss: Incorrect type");
1388
1389 Compute(theSurface,
1390 theLocation,
1391 NULL,
1392 Standard_True,
1393 theOutMass,
1394 theOutGravityCenter,
1395 theOutInertia);
1396 }
1397