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