b311480e |
1 | // Created on: 2005-12-09 |
2 | // Created by: Sergey KHROMOV |
973c2be1 |
3 | // Copyright (c) 2005-2014 OPEN CASCADE SAS |
b311480e |
4 | // |
973c2be1 |
5 | // This file is part of Open CASCADE Technology software library. |
b311480e |
6 | // |
d5f74e42 |
7 | // This library is free software; you can redistribute it and/or modify it under |
8 | // the terms of the GNU Lesser General Public License version 2.1 as published |
973c2be1 |
9 | // by the Free Software Foundation, with special exception defined in the file |
10 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
11 | // distribution for complete text of the license and disclaimer of any warranty. |
b311480e |
12 | // |
973c2be1 |
13 | // Alternatively, this file may be used under the terms of Open CASCADE |
14 | // commercial license or contractual agreement. |
7fd59977 |
15 | |
424cd6bb |
16 | |
42cf5bc1 |
17 | #include <BRepGProp_Domain.hxx> |
18 | #include <BRepGProp_Face.hxx> |
424cd6bb |
19 | #include <BRepGProp_TFunction.hxx> |
42cf5bc1 |
20 | #include <BRepGProp_VinertGK.hxx> |
21 | #include <gp_Pln.hxx> |
22 | #include <gp_Pnt.hxx> |
23 | #include <math_KronrodSingleIntegration.hxx> |
24 | #include <TColStd_Array1OfBoolean.hxx> |
25 | #include <TColStd_HArray1OfReal.hxx> |
7fd59977 |
26 | |
27 | //========================================================================== |
28 | //function : Constructor |
7fd59977 |
29 | //========================================================================== |
424cd6bb |
30 | BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face &theSurface, |
31 | const gp_Pnt &theLocation, |
32 | const Standard_Real theTolerance, |
33 | const Standard_Boolean theCGFlag, |
34 | const Standard_Boolean theIFlag): |
35 | myErrorReached(0.) |
7fd59977 |
36 | { |
37 | SetLocation(theLocation); |
38 | Perform(theSurface, theTolerance, theCGFlag, theIFlag); |
39 | } |
40 | |
41 | //========================================================================== |
42 | //function : Constructor |
43 | // |
44 | //========================================================================== |
45 | |
424cd6bb |
46 | BRepGProp_VinertGK::BRepGProp_VinertGK( BRepGProp_Face &theSurface, |
47 | const gp_Pnt &thePoint, |
48 | const gp_Pnt &theLocation, |
49 | const Standard_Real theTolerance, |
50 | const Standard_Boolean theCGFlag, |
51 | const Standard_Boolean theIFlag): |
52 | myErrorReached(0.) |
7fd59977 |
53 | { |
54 | SetLocation(theLocation); |
55 | Perform(theSurface, thePoint, theTolerance, theCGFlag, theIFlag); |
56 | } |
57 | |
58 | //========================================================================== |
59 | //function : Constructor |
60 | // |
61 | //========================================================================== |
62 | |
424cd6bb |
63 | BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face &theSurface, |
64 | BRepGProp_Domain &theDomain, |
65 | const gp_Pnt &theLocation, |
66 | const Standard_Real theTolerance, |
67 | const Standard_Boolean theCGFlag, |
68 | const Standard_Boolean theIFlag): |
69 | myErrorReached(0.) |
7fd59977 |
70 | { |
71 | SetLocation(theLocation); |
72 | Perform(theSurface, theDomain, theTolerance, theCGFlag, theIFlag); |
73 | } |
74 | |
75 | //========================================================================== |
76 | //function : Constructor |
77 | // |
78 | //========================================================================== |
79 | |
424cd6bb |
80 | BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face &theSurface, |
81 | BRepGProp_Domain &theDomain, |
82 | const gp_Pnt &thePoint, |
83 | const gp_Pnt &theLocation, |
84 | const Standard_Real theTolerance, |
85 | const Standard_Boolean theCGFlag, |
86 | const Standard_Boolean theIFlag): |
87 | myErrorReached(0.) |
7fd59977 |
88 | { |
89 | SetLocation(theLocation); |
90 | Perform(theSurface, theDomain, thePoint, theTolerance, theCGFlag, theIFlag); |
91 | } |
92 | |
93 | //========================================================================== |
94 | //function : Constructor |
95 | // |
96 | //========================================================================== |
97 | |
424cd6bb |
98 | BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face &theSurface, |
99 | const gp_Pln &thePlane, |
100 | const gp_Pnt &theLocation, |
101 | const Standard_Real theTolerance, |
102 | const Standard_Boolean theCGFlag, |
103 | const Standard_Boolean theIFlag): |
104 | myErrorReached(0.) |
7fd59977 |
105 | { |
106 | SetLocation(theLocation); |
107 | Perform(theSurface, thePlane, theTolerance, theCGFlag, theIFlag); |
108 | } |
109 | |
110 | //========================================================================== |
111 | //function : Constructor |
112 | // |
113 | //========================================================================== |
114 | |
424cd6bb |
115 | BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face &theSurface, |
116 | BRepGProp_Domain &theDomain, |
117 | const gp_Pln &thePlane, |
118 | const gp_Pnt &theLocation, |
119 | const Standard_Real theTolerance, |
120 | const Standard_Boolean theCGFlag, |
121 | const Standard_Boolean theIFlag): |
122 | myErrorReached(0.) |
7fd59977 |
123 | { |
124 | SetLocation(theLocation); |
125 | Perform(theSurface, theDomain, thePlane, theTolerance, theCGFlag, theIFlag); |
126 | } |
127 | |
128 | //========================================================================== |
129 | //function : Perform |
130 | // Compute the properties. |
131 | //========================================================================== |
132 | |
424cd6bb |
133 | Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face &theSurface, |
134 | const Standard_Real theTolerance, |
135 | const Standard_Boolean theCGFlag, |
136 | const Standard_Boolean theIFlag) |
7fd59977 |
137 | |
138 | { |
139 | Standard_Real aShift[] = { 0., 0., 0. }; |
140 | |
141 | return PrivatePerform(theSurface, NULL, Standard_True, &aShift, theTolerance, |
424cd6bb |
142 | theCGFlag, theIFlag); |
7fd59977 |
143 | } |
144 | |
145 | //========================================================================== |
146 | //function : Perform |
147 | // Compute the properties. |
148 | //========================================================================== |
149 | |
424cd6bb |
150 | Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face &theSurface, |
151 | const gp_Pnt &thePoint, |
152 | const Standard_Real theTolerance, |
153 | const Standard_Boolean theCGFlag, |
154 | const Standard_Boolean theIFlag) |
7fd59977 |
155 | |
156 | { |
157 | gp_XYZ aXYZ(thePoint.XYZ().Subtracted(loc.XYZ())); |
158 | Standard_Real aShift[3]; |
159 | |
160 | aXYZ.Coord(aShift[0], aShift[1], aShift[2]); |
161 | |
162 | return PrivatePerform(theSurface, NULL, Standard_True, &aShift, theTolerance, |
424cd6bb |
163 | theCGFlag, theIFlag); |
7fd59977 |
164 | } |
165 | |
166 | //========================================================================== |
167 | //function : Perform |
168 | // Compute the properties. |
169 | //========================================================================== |
170 | |
424cd6bb |
171 | Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face &theSurface, |
172 | BRepGProp_Domain &theDomain, |
173 | const Standard_Real theTolerance, |
174 | const Standard_Boolean theCGFlag, |
175 | const Standard_Boolean theIFlag) |
7fd59977 |
176 | |
177 | { |
178 | Standard_Real aShift[] = { 0., 0., 0. }; |
179 | |
180 | return PrivatePerform(theSurface, &theDomain, |
424cd6bb |
181 | Standard_True, &aShift, theTolerance, |
182 | theCGFlag, theIFlag); |
7fd59977 |
183 | } |
184 | |
185 | //========================================================================== |
186 | //function : Perform |
187 | // Compute the properties. |
188 | //========================================================================== |
189 | |
424cd6bb |
190 | Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face &theSurface, |
191 | BRepGProp_Domain &theDomain, |
192 | const gp_Pnt &thePoint, |
193 | const Standard_Real theTolerance, |
194 | const Standard_Boolean theCGFlag, |
195 | const Standard_Boolean theIFlag) |
7fd59977 |
196 | |
197 | { |
198 | gp_XYZ aXYZ(thePoint.XYZ().Subtracted(loc.XYZ())); |
199 | Standard_Real aShift[3]; |
200 | |
201 | aXYZ.Coord(aShift[0], aShift[1], aShift[2]); |
202 | |
203 | return PrivatePerform(theSurface, &theDomain, |
424cd6bb |
204 | Standard_True, &aShift, theTolerance, |
205 | theCGFlag, theIFlag); |
7fd59977 |
206 | } |
207 | |
208 | //========================================================================== |
209 | //function : Perform |
210 | // Compute the properties. |
211 | //========================================================================== |
212 | |
424cd6bb |
213 | Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face &theSurface, |
214 | const gp_Pln &thePlane, |
215 | const Standard_Real theTolerance, |
216 | const Standard_Boolean theCGFlag, |
217 | const Standard_Boolean theIFlag) |
7fd59977 |
218 | |
219 | { |
220 | Standard_Real aCoeff[4]; |
221 | Standard_Real aXLoc; |
222 | Standard_Real aYLoc; |
223 | Standard_Real aZLoc; |
224 | |
225 | loc.Coord(aXLoc, aYLoc, aZLoc); |
226 | thePlane.Coefficients (aCoeff[0], aCoeff[1], aCoeff[2], aCoeff[3]); |
227 | aCoeff[3] = aCoeff[3] - aCoeff[0]*aXLoc - aCoeff[1]*aYLoc - aCoeff[2]*aZLoc; |
228 | |
229 | return PrivatePerform(theSurface, NULL, |
424cd6bb |
230 | Standard_False, &aCoeff, theTolerance, |
231 | theCGFlag, theIFlag); |
7fd59977 |
232 | } |
233 | |
234 | //========================================================================== |
235 | //function : Perform |
236 | // Compute the properties. |
237 | //========================================================================== |
238 | |
424cd6bb |
239 | Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face &theSurface, |
240 | BRepGProp_Domain &theDomain, |
241 | const gp_Pln &thePlane, |
242 | const Standard_Real theTolerance, |
243 | const Standard_Boolean theCGFlag, |
244 | const Standard_Boolean theIFlag) |
7fd59977 |
245 | |
246 | { |
247 | Standard_Real aCoeff[4]; |
248 | Standard_Real aXLoc; |
249 | Standard_Real aYLoc; |
250 | Standard_Real aZLoc; |
251 | |
252 | loc.Coord(aXLoc, aYLoc, aZLoc); |
253 | thePlane.Coefficients (aCoeff[0], aCoeff[1], aCoeff[2], aCoeff[3]); |
254 | aCoeff[3] = aCoeff[3] - aCoeff[0]*aXLoc - aCoeff[1]*aYLoc - aCoeff[2]*aZLoc; |
255 | |
256 | return PrivatePerform(theSurface, &theDomain, |
424cd6bb |
257 | Standard_False, &aCoeff, theTolerance, |
258 | theCGFlag, theIFlag); |
7fd59977 |
259 | } |
260 | |
261 | //========================================================================== |
262 | //function : PrivatePerform |
263 | // Compute the properties. |
264 | //========================================================================== |
265 | |
424cd6bb |
266 | Standard_Real BRepGProp_VinertGK::PrivatePerform |
267 | (BRepGProp_Face &theSurface, |
268 | const Standard_Address thePtrDomain, |
269 | const Standard_Boolean IsByPoint, |
270 | const Standard_Address theCoeffs, |
271 | const Standard_Real theTolerance, |
272 | const Standard_Boolean theCGFlag, |
273 | const Standard_Boolean theIFlag) |
7fd59977 |
274 | |
275 | { |
276 | |
277 | const Standard_Real aTTol = 1.e-9; |
278 | Standard_Real *aCoeffs = (Standard_Real *)theCoeffs; |
279 | |
280 | // Compute the number of 2d bounding curves of the face. |
424cd6bb |
281 | BRepGProp_Domain *aPDomain = NULL; |
7fd59977 |
282 | Standard_Integer aNbCurves = 0; |
283 | |
284 | // If the pointer to the domain is NULL, there is only one curve to treat: |
285 | // U isoline with the UMax parameter. |
286 | if (thePtrDomain == NULL) |
287 | aNbCurves = 1; |
288 | else { |
424cd6bb |
289 | aPDomain = (BRepGProp_Domain *)thePtrDomain; |
7fd59977 |
290 | |
291 | for (aPDomain->Init(); aPDomain->More(); aPDomain->Next()) |
292 | aNbCurves++; |
293 | } |
294 | |
295 | if (aNbCurves == 0) { |
296 | myErrorReached = -1.; |
297 | |
298 | return myErrorReached; |
299 | } |
300 | |
301 | //Standard_Real aCrvTol = 0.5*theTolerance/aNbCurves; |
302 | Standard_Real aCrvTol = 0.1*theTolerance; |
303 | Standard_Real aUMin; |
304 | Standard_Real aUMax; |
305 | Standard_Real aTMin; |
306 | Standard_Real aTMax; |
307 | Standard_Integer aNbPnts; |
308 | Standard_Integer aNbMaxIter = 1000; |
309 | Standard_Integer aNbVal = 10; |
310 | Standard_Integer k; |
311 | math_Vector aLocalValue(1, aNbVal); |
312 | math_Vector aLocalTolReached(1, aNbVal); |
313 | math_Vector aValue(1, aNbVal); |
314 | math_Vector aTolReached(1, aNbVal); |
315 | TColStd_Array1OfBoolean CFlags(1, aNbVal); |
316 | CFlags.Init(Standard_False); |
317 | Standard_Boolean isMore; |
318 | |
319 | //aNbVal = 1; |
320 | aValue.Init(0.); |
321 | aTolReached.Init(0.); |
322 | |
323 | CFlags.Init(Standard_False); |
324 | CFlags(1) = Standard_True; |
325 | |
326 | if(theCGFlag || theIFlag) { |
327 | Standard_Integer i; |
328 | for(i = 2; i <= 4; ++i) {CFlags(i) = Standard_True;} |
329 | } |
330 | |
331 | if(theIFlag) { |
332 | Standard_Integer i; |
333 | for(i = 5; i <= 10; ++i) {CFlags(i) = Standard_True;} |
334 | } |
335 | |
336 | theSurface.Bounds(aUMin, aUMax, aTMin, aTMax); |
337 | |
338 | if (thePtrDomain == NULL) |
339 | isMore = Standard_True; |
340 | else { |
341 | aPDomain->Init(); |
342 | isMore = aPDomain->More(); |
343 | } |
344 | |
345 | while(isMore) { |
346 | // If the pointer to the domain is NULL, there is only one curve to treat: |
347 | // U isoline with the UMax parameter. |
348 | |
349 | if (thePtrDomain == NULL) |
350 | theSurface.Load(Standard_False, GeomAbs_IsoU); |
351 | else |
352 | theSurface.Load(aPDomain->Value()); |
353 | |
354 | aTMin = theSurface.FirstParameter(); |
355 | aTMax = theSurface.LastParameter(); |
356 | |
357 | |
358 | // Get the spans on the curve. |
359 | Handle(TColStd_HArray1OfReal) aTKnots; |
424cd6bb |
360 | BRepGProp_TFunction aTFunc(theSurface, loc, IsByPoint, theCoeffs, |
361 | aUMin, aCrvTol); |
7fd59977 |
362 | |
363 | theSurface.GetTKnots(aTMin, aTMax, aTKnots); |
364 | |
365 | Standard_Integer iU = aTKnots->Upper(); |
366 | Standard_Integer aNbTIntervals = aTKnots->Length() - 1; |
367 | //Standard_Real aTolSpan = aCrvTol/aNbTIntervals; |
368 | Standard_Real aTolSpan = 0.9*theTolerance; //Relative error |
369 | math_KronrodSingleIntegration anIntegral; |
370 | GProp_ValueType aValueType; |
371 | |
372 | |
373 | // Empirical criterion. |
374 | aNbPnts = Min(15, theSurface.IntegrationOrder()/aNbTIntervals + 1); |
375 | aNbPnts = Max(5, aNbPnts); |
424cd6bb |
376 | // aNbPnts = theSurface.IntegrationOrder(); |
7fd59977 |
377 | |
378 | aLocalValue.Init(0.); |
379 | aLocalTolReached.Init(0.); |
380 | |
381 | for (k = 1; k <= aNbVal; k++) { |
382 | |
383 | if(!CFlags(k)) continue; |
384 | |
385 | Standard_Integer i = aTKnots->Lower(); |
386 | |
387 | switch (k) { |
388 | case 1: aValueType = GProp_Mass; break; |
389 | case 2: aValueType = GProp_CenterMassX; break; |
390 | case 3: aValueType = GProp_CenterMassY; break; |
391 | case 4: aValueType = GProp_CenterMassZ; break; |
392 | case 5: aValueType = GProp_InertiaXX; break; |
393 | case 6: aValueType = GProp_InertiaYY; break; |
394 | case 7: aValueType = GProp_InertiaZZ; break; |
395 | case 8: aValueType = GProp_InertiaXY; break; |
396 | case 9: aValueType = GProp_InertiaXZ; break; |
397 | case 10: aValueType = GProp_InertiaYZ; break; |
398 | |
399 | default: myErrorReached = -1.; return myErrorReached; |
400 | } |
401 | aTFunc.SetValueType(aValueType); |
402 | |
403 | Standard_Real err1 = 0.; |
404 | while (i < iU) { |
405 | |
424cd6bb |
406 | //cout << "-------------- Span " << i << " nbp: " << aNbPnts << endl; |
407 | Standard_Real aT1 = aTKnots->Value(i++); |
408 | Standard_Real aT2 = aTKnots->Value(i); |
7fd59977 |
409 | |
424cd6bb |
410 | if(aT2 - aT1 < aTTol) continue; |
7fd59977 |
411 | |
424cd6bb |
412 | aTFunc.SetNbKronrodPoints(aNbPnts); |
413 | aTFunc.Init(); |
414 | aTFunc.SetTolerance(aCrvTol/(aT2-aT1)); |
415 | anIntegral.Perform(aTFunc, aT1, aT2, aNbPnts, aTolSpan, aNbMaxIter); |
7fd59977 |
416 | |
424cd6bb |
417 | if (!anIntegral.IsDone()) { |
418 | myErrorReached = -1.; |
7fd59977 |
419 | |
424cd6bb |
420 | return myErrorReached; |
421 | } |
7fd59977 |
422 | |
424cd6bb |
423 | aLocalValue(k) += anIntegral.Value(); |
424 | err1 = aTFunc.AbsolutError()*(aT2 - aT1); |
425 | //cout << "Errors: " << anIntegral.NbIterReached() << " " << anIntegral.AbsolutError() << " " << err1 << endl; |
426 | aLocalTolReached(k) += anIntegral.AbsolutError() + err1; |
427 | //cout << "--- Errors: " << anIntegral.NbIterReached() << " " << anIntegral.AbsolutError() << " " << err1 << endl; |
7fd59977 |
428 | } |
429 | |
430 | aValue(k) += aLocalValue(k); |
431 | aTolReached(k) += aLocalTolReached(k); |
432 | } |
433 | |
434 | // If the pointer to the domain is NULL, there is only one curve to treat: |
435 | // U isoline with the UMax parameter. |
436 | if (thePtrDomain == NULL) |
437 | isMore = Standard_False; |
438 | else { |
439 | aPDomain->Next(); |
440 | isMore = aPDomain->More(); |
441 | } |
442 | } |
443 | |
444 | // Get volume value. |
445 | dim = aValue(1); |
446 | myErrorReached = aTolReached(1); |
447 | myAbsolutError = myErrorReached; |
448 | Standard_Real anAbsDim = Abs(dim); |
449 | Standard_Real aVolTol = Epsilon(myAbsolutError); |
450 | if(anAbsDim >= aVolTol) myErrorReached /= anAbsDim; |
451 | |
452 | if(theCGFlag || theIFlag) { |
453 | // Compute values of center of mass. |
454 | if(anAbsDim >= aVolTol) { |
455 | if (IsByPoint) { |
424cd6bb |
456 | aValue(2) = aCoeffs[0] + aValue(2)/dim; |
457 | aValue(3) = aCoeffs[1] + aValue(3)/dim; |
458 | aValue(4) = aCoeffs[2] + aValue(4)/dim; |
7fd59977 |
459 | } else { |
424cd6bb |
460 | aValue(2) /= dim; |
461 | aValue(3) /= dim; |
462 | aValue(4) /= dim; |
7fd59977 |
463 | } |
464 | } else { |
465 | aValue(2) = 0.; |
466 | aValue(3) = 0.; |
467 | aValue(4) = 0.; |
468 | dim = 0.; |
469 | } |
470 | g.SetCoord(aValue(2), aValue(3), aValue(4)); |
471 | } |
472 | |
473 | if(theIFlag) { |
474 | // Fill the matrix of inertia. |
475 | inertia.SetCols (gp_XYZ (aValue(5), aValue(8), aValue(9)), |
424cd6bb |
476 | gp_XYZ (aValue(8), aValue(6), aValue(10)), |
477 | gp_XYZ (aValue(9), aValue(10), aValue(7))); |
7fd59977 |
478 | } |
479 | //return myErrorReached; |
480 | return myAbsolutError; |
481 | } |
424cd6bb |
482 | |