1 // Created on: 2005-12-09
2 // Created by: Sergey KHROMOV
3 // Copyright (c) 2005-2014 OPEN CASCADE SAS
5 // This file is part of Open CASCADE Technology software library.
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
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.
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
17 #include <BRepGProp_Domain.hxx>
18 #include <BRepGProp_Face.hxx>
19 #include <BRepGProp_TFunction.hxx>
20 #include <BRepGProp_VinertGK.hxx>
23 #include <math_KronrodSingleIntegration.hxx>
24 #include <TColStd_Array1OfBoolean.hxx>
25 #include <TColStd_HArray1OfReal.hxx>
27 //==========================================================================
28 //function : Constructor
29 //==========================================================================
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):
37 SetLocation(theLocation);
38 Perform(theSurface, theTolerance, theCGFlag, theIFlag);
41 //==========================================================================
42 //function : Constructor
44 //==========================================================================
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):
54 SetLocation(theLocation);
55 Perform(theSurface, thePoint, theTolerance, theCGFlag, theIFlag);
58 //==========================================================================
59 //function : Constructor
61 //==========================================================================
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):
71 SetLocation(theLocation);
72 Perform(theSurface, theDomain, theTolerance, theCGFlag, theIFlag);
75 //==========================================================================
76 //function : Constructor
78 //==========================================================================
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):
89 SetLocation(theLocation);
90 Perform(theSurface, theDomain, thePoint, theTolerance, theCGFlag, theIFlag);
93 //==========================================================================
94 //function : Constructor
96 //==========================================================================
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):
106 SetLocation(theLocation);
107 Perform(theSurface, thePlane, theTolerance, theCGFlag, theIFlag);
110 //==========================================================================
111 //function : Constructor
113 //==========================================================================
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):
124 SetLocation(theLocation);
125 Perform(theSurface, theDomain, thePlane, theTolerance, theCGFlag, theIFlag);
128 //==========================================================================
130 // Compute the properties.
131 //==========================================================================
133 Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face &theSurface,
134 const Standard_Real theTolerance,
135 const Standard_Boolean theCGFlag,
136 const Standard_Boolean theIFlag)
139 Standard_Real aShift[] = { 0., 0., 0. };
141 return PrivatePerform(theSurface, NULL, Standard_True, &aShift, theTolerance,
142 theCGFlag, theIFlag);
145 //==========================================================================
147 // Compute the properties.
148 //==========================================================================
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)
157 gp_XYZ aXYZ(thePoint.XYZ().Subtracted(loc.XYZ()));
158 Standard_Real aShift[3];
160 aXYZ.Coord(aShift[0], aShift[1], aShift[2]);
162 return PrivatePerform(theSurface, NULL, Standard_True, &aShift, theTolerance,
163 theCGFlag, theIFlag);
166 //==========================================================================
168 // Compute the properties.
169 //==========================================================================
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)
178 Standard_Real aShift[] = { 0., 0., 0. };
180 return PrivatePerform(theSurface, &theDomain,
181 Standard_True, &aShift, theTolerance,
182 theCGFlag, theIFlag);
185 //==========================================================================
187 // Compute the properties.
188 //==========================================================================
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)
198 gp_XYZ aXYZ(thePoint.XYZ().Subtracted(loc.XYZ()));
199 Standard_Real aShift[3];
201 aXYZ.Coord(aShift[0], aShift[1], aShift[2]);
203 return PrivatePerform(theSurface, &theDomain,
204 Standard_True, &aShift, theTolerance,
205 theCGFlag, theIFlag);
208 //==========================================================================
210 // Compute the properties.
211 //==========================================================================
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)
220 Standard_Real aCoeff[4];
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;
229 return PrivatePerform(theSurface, NULL,
230 Standard_False, &aCoeff, theTolerance,
231 theCGFlag, theIFlag);
234 //==========================================================================
236 // Compute the properties.
237 //==========================================================================
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)
247 Standard_Real aCoeff[4];
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;
256 return PrivatePerform(theSurface, &theDomain,
257 Standard_False, &aCoeff, theTolerance,
258 theCGFlag, theIFlag);
261 //==========================================================================
262 //function : PrivatePerform
263 // Compute the properties.
264 //==========================================================================
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)
277 const Standard_Real aTTol = 1.e-9;
278 Standard_Real *aCoeffs = (Standard_Real *)theCoeffs;
280 // Compute the number of 2d bounding curves of the face.
281 BRepGProp_Domain *aPDomain = NULL;
282 Standard_Integer aNbCurves = 0;
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)
289 aPDomain = (BRepGProp_Domain *)thePtrDomain;
291 for (aPDomain->Init(); aPDomain->More(); aPDomain->Next())
295 if (aNbCurves == 0) {
296 myErrorReached = -1.;
298 return myErrorReached;
301 //Standard_Real aCrvTol = 0.5*theTolerance/aNbCurves;
302 Standard_Real aCrvTol = 0.1*theTolerance;
307 Standard_Integer aNbPnts;
308 Standard_Integer aNbMaxIter = 1000;
309 Standard_Integer aNbVal = 10;
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;
321 aTolReached.Init(0.);
323 CFlags.Init(Standard_False);
324 CFlags(1) = Standard_True;
326 if(theCGFlag || theIFlag) {
328 for(i = 2; i <= 4; ++i) {CFlags(i) = Standard_True;}
333 for(i = 5; i <= 10; ++i) {CFlags(i) = Standard_True;}
336 theSurface.Bounds(aUMin, aUMax, aTMin, aTMax);
338 if (thePtrDomain == NULL)
339 isMore = Standard_True;
342 isMore = aPDomain->More();
346 // If the pointer to the domain is NULL, there is only one curve to treat:
347 // U isoline with the UMax parameter.
349 if (thePtrDomain == NULL)
350 theSurface.Load(Standard_False, GeomAbs_IsoU);
352 theSurface.Load(aPDomain->Value());
354 aTMin = theSurface.FirstParameter();
355 aTMax = theSurface.LastParameter();
358 // Get the spans on the curve.
359 Handle(TColStd_HArray1OfReal) aTKnots;
360 BRepGProp_TFunction aTFunc(theSurface, loc, IsByPoint, theCoeffs,
363 theSurface.GetTKnots(aTMin, aTMax, aTKnots);
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;
373 // Empirical criterion.
374 aNbPnts = Min(15, theSurface.IntegrationOrder()/aNbTIntervals + 1);
375 aNbPnts = Max(5, aNbPnts);
376 // aNbPnts = theSurface.IntegrationOrder();
378 aLocalValue.Init(0.);
379 aLocalTolReached.Init(0.);
381 for (k = 1; k <= aNbVal; k++) {
383 if(!CFlags(k)) continue;
385 Standard_Integer i = aTKnots->Lower();
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;
399 default: myErrorReached = -1.; return myErrorReached;
401 aTFunc.SetValueType(aValueType);
403 Standard_Real err1 = 0.;
406 //cout << "-------------- Span " << i << " nbp: " << aNbPnts << endl;
407 Standard_Real aT1 = aTKnots->Value(i++);
408 Standard_Real aT2 = aTKnots->Value(i);
410 if(aT2 - aT1 < aTTol) continue;
412 aTFunc.SetNbKronrodPoints(aNbPnts);
414 aTFunc.SetTolerance(aCrvTol/(aT2-aT1));
415 anIntegral.Perform(aTFunc, aT1, aT2, aNbPnts, aTolSpan, aNbMaxIter);
417 if (!anIntegral.IsDone()) {
418 myErrorReached = -1.;
420 return myErrorReached;
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;
430 aValue(k) += aLocalValue(k);
431 aTolReached(k) += aLocalTolReached(k);
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;
440 isMore = aPDomain->More();
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;
452 if(theCGFlag || theIFlag) {
453 // Compute values of center of mass.
454 if(anAbsDim >= aVolTol) {
456 aValue(2) = aCoeffs[0] + aValue(2)/dim;
457 aValue(3) = aCoeffs[1] + aValue(3)/dim;
458 aValue(4) = aCoeffs[2] + aValue(4)/dim;
470 g.SetCoord(aValue(2), aValue(3), aValue(4));
474 // Fill the matrix of inertia.
475 inertia.SetCols (gp_XYZ (aValue(5), aValue(8), aValue(9)),
476 gp_XYZ (aValue(8), aValue(6), aValue(10)),
477 gp_XYZ (aValue(9), aValue(10), aValue(7)));
479 //return myErrorReached;
480 return myAbsolutError;