0024002: Overall code and build procedure refactoring -- automatic
[occt.git] / src / BRepGProp / BRepGProp_VinertGK.cxx
1 // Created on: 2005-12-09
2 // Created by: Sergey KHROMOV
3 // Copyright (c) 2005-2014 OPEN CASCADE SAS
4 //
5 // This file is part of Open CASCADE Technology software library.
6 //
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.
12 //
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
15
16
17 #include <BRepGProp_Domain.hxx>
18 #include <BRepGProp_Face.hxx>
19 #include <BRepGProp_TFunction.hxx>
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>
26
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):
35   myErrorReached(0.)
36 {
37   SetLocation(theLocation);
38   Perform(theSurface, theTolerance, theCGFlag, theIFlag);
39 }
40
41 //==========================================================================
42 //function : Constructor
43 //           
44 //==========================================================================
45
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.)
53 {
54   SetLocation(theLocation);
55   Perform(theSurface, thePoint, theTolerance, theCGFlag, theIFlag);
56 }
57
58 //==========================================================================
59 //function : Constructor
60 //           
61 //==========================================================================
62
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.)
70 {
71   SetLocation(theLocation);
72   Perform(theSurface, theDomain, theTolerance, theCGFlag, theIFlag);
73 }
74
75 //==========================================================================
76 //function : Constructor
77 //           
78 //==========================================================================
79
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.)
88 {
89   SetLocation(theLocation);
90   Perform(theSurface, theDomain, thePoint, theTolerance, theCGFlag, theIFlag);
91 }
92
93 //==========================================================================
94 //function : Constructor
95 //           
96 //==========================================================================
97
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.)
105 {
106   SetLocation(theLocation);
107   Perform(theSurface, thePlane, theTolerance, theCGFlag, theIFlag);
108 }
109
110 //==========================================================================
111 //function : Constructor
112 //           
113 //==========================================================================
114
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.)
123 {
124   SetLocation(theLocation);
125   Perform(theSurface, theDomain, thePlane, theTolerance, theCGFlag, theIFlag);
126 }
127
128 //==========================================================================
129 //function : Perform
130 //           Compute the properties.
131 //==========================================================================
132
133 Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face        &theSurface,
134                                           const Standard_Real    theTolerance,
135                                           const Standard_Boolean theCGFlag,
136                                           const Standard_Boolean theIFlag)
137
138 {
139   Standard_Real aShift[] = { 0., 0., 0. };
140
141   return PrivatePerform(theSurface, NULL, Standard_True, &aShift, theTolerance, 
142     theCGFlag, theIFlag);
143 }
144
145 //==========================================================================
146 //function : Perform
147 //           Compute the properties.
148 //==========================================================================
149
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)
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, 
163     theCGFlag, theIFlag);
164 }
165
166 //==========================================================================
167 //function : Perform
168 //           Compute the properties.
169 //==========================================================================
170
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)
176
177 {
178   Standard_Real aShift[] = { 0., 0., 0. };
179
180   return PrivatePerform(theSurface, &theDomain,
181     Standard_True, &aShift, theTolerance, 
182     theCGFlag, theIFlag);
183 }
184
185 //==========================================================================
186 //function : Perform
187 //           Compute the properties.
188 //==========================================================================
189
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)
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,
204     Standard_True, &aShift, theTolerance, 
205     theCGFlag, theIFlag);
206 }
207
208 //==========================================================================
209 //function : Perform
210 //           Compute the properties.
211 //==========================================================================
212
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)
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,
230     Standard_False, &aCoeff, theTolerance, 
231     theCGFlag, theIFlag);
232 }
233
234 //==========================================================================
235 //function : Perform
236 //           Compute the properties.
237 //==========================================================================
238
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)
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,
257     Standard_False, &aCoeff, theTolerance, 
258     theCGFlag, theIFlag);
259 }
260
261 //==========================================================================
262 //function : PrivatePerform
263 //           Compute the properties.
264 //==========================================================================
265
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)
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.
281   BRepGProp_Domain           *aPDomain = NULL;
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 {
289     aPDomain = (BRepGProp_Domain *)thePtrDomain;
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;
360     BRepGProp_TFunction               aTFunc(theSurface, loc, IsByPoint, theCoeffs,
361       aUMin, aCrvTol);
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);
376     //     aNbPnts = theSurface.IntegrationOrder();
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
406         //cout << "-------------- Span " << i << " nbp: " << aNbPnts << endl;
407         Standard_Real aT1 = aTKnots->Value(i++);
408         Standard_Real aT2 = aTKnots->Value(i);
409
410         if(aT2 - aT1 < aTTol) continue;
411
412         aTFunc.SetNbKronrodPoints(aNbPnts);
413         aTFunc.Init();
414         aTFunc.SetTolerance(aCrvTol/(aT2-aT1));
415         anIntegral.Perform(aTFunc, aT1, aT2, aNbPnts, aTolSpan, aNbMaxIter);
416
417         if (!anIntegral.IsDone()) {
418           myErrorReached = -1.;
419
420           return myErrorReached;
421         }
422
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;
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) {
456         aValue(2) = aCoeffs[0] + aValue(2)/dim;
457         aValue(3) = aCoeffs[1] + aValue(3)/dim;
458         aValue(4) = aCoeffs[2] + aValue(4)/dim;
459       } else {
460         aValue(2) /= dim;
461         aValue(3) /= dim;
462         aValue(4) /= dim;
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)),
476       gp_XYZ (aValue(8), aValue(6),  aValue(10)),
477       gp_XYZ (aValue(9), aValue(10), aValue(7)));
478   }
479   //return myErrorReached;
480   return myAbsolutError;
481 }
482