0024774: Convertation of the generic classes to the non-generic. Part 8
[occt.git] / src / BRepGProp / BRepGProp_TFunction.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 #include <BRepGProp_TFunction.ixx>
17
18 #include <TColStd_HArray1OfReal.hxx>
19 #include <math_KronrodSingleIntegration.hxx>
20 #include <Precision.hxx>
21
22 //=======================================================================
23 //function : Constructor.
24 //purpose  : 
25 //=======================================================================
26
27 BRepGProp_TFunction::BRepGProp_TFunction(const BRepGProp_Face   &theSurface,
28                                          const gp_Pnt           &theVertex,
29                                          const Standard_Boolean  IsByPoint,
30                                          const Standard_Address  theCoeffs,
31                                          const Standard_Real     theUMin,
32                                          const Standard_Real     theTolerance):
33   mySurface(theSurface),
34   myUFunction(mySurface, theVertex, IsByPoint, theCoeffs),
35   myUMin(theUMin),
36   myTolerance(theTolerance),
37   myTolReached(0.),
38   myErrReached(0.),
39   myAbsError(0.),
40   myValueType(GProp_Unknown),
41   myIsByPoint(IsByPoint),
42   myNbPntOuter(3)
43 {
44 }
45
46 //=======================================================================
47 //function : Init
48 //purpose  : 
49 //=======================================================================
50
51 void BRepGProp_TFunction::Init() 
52 {
53   myTolReached = 0.;
54   myErrReached = 0.;
55   myAbsError = 0.;
56 }
57
58 //=======================================================================
59 //function : Value
60 //purpose  : 
61 //=======================================================================
62
63 Standard_Boolean BRepGProp_TFunction::Value(const Standard_Real X,
64                                             Standard_Real &F)
65 {
66   const Standard_Real tolU = 1.e-9;
67
68   gp_Pnt2d                      aP2d;
69   gp_Vec2d                      aV2d;
70   Standard_Real                 aUMax;
71   Handle(TColStd_HArray1OfReal) anUKnots;
72
73   mySurface.D12d(X, aP2d, aV2d);
74   aUMax = aP2d.X();
75
76   if(aUMax - myUMin < tolU)
77   {
78     F = 0.;
79     return Standard_True;
80   }
81
82   mySurface.GetUKnots(myUMin, aUMax, anUKnots);
83   myUFunction.SetVParam(aP2d.Y());
84
85   // Compute the integral from myUMin to aUMax of myUFunction.
86   Standard_Integer i;
87   Standard_Real    aCoeff        = aV2d.Y();
88   //Standard_Integer aNbUIntervals = anUKnots->Length() - 1;
89   //Standard_Real    aTol          = myTolerance/aNbUIntervals;
90   Standard_Real    aTol          = myTolerance;
91
92   //aTol /= myNbPntOuter;
93
94   if (myValueType == GProp_Mass) {
95     if (myIsByPoint)
96       aCoeff /= 3.;
97   } else if (myValueType == GProp_CenterMassX ||
98     myValueType == GProp_CenterMassY ||
99     myValueType == GProp_CenterMassZ) {
100       if (myIsByPoint)
101         aCoeff *= 0.25;
102   } else if (myValueType == GProp_InertiaXX ||
103     myValueType == GProp_InertiaYY ||
104     myValueType == GProp_InertiaZZ ||
105     myValueType == GProp_InertiaXY ||
106     myValueType == GProp_InertiaXZ ||
107     myValueType == GProp_InertiaYZ) {
108       if (myIsByPoint)
109         aCoeff *= 0.2;
110   } else
111     return Standard_False;
112
113   Standard_Real aAbsCoeff = Abs(aCoeff);
114
115   if (aAbsCoeff <= Precision::Angular()) {
116     // No need to compute the integral. The value will be equal to 0.
117     F = 0.;
118     return Standard_True;
119   } 
120   //else if (aAbsCoeff > 10.*aTol)
121   //  aTol /= aAbsCoeff;
122   //else
123   //  aTol = 0.1;
124
125   Standard_Integer              iU           = anUKnots->Upper();
126   Standard_Integer              aNbPntsStart;
127   Standard_Integer              aNbMaxIter   = 1000;
128   math_KronrodSingleIntegration anIntegral;
129   Standard_Real                 aLocalErr    = 0.;
130
131   i            = anUKnots->Lower();
132   F            = 0.;
133
134   // Epmirical criterion
135   aNbPntsStart = Min(15, mySurface.UIntegrationOrder()/(anUKnots->Length() - 1)+1);
136   aNbPntsStart = Max(5, aNbPntsStart);
137
138   while (i < iU) {
139     Standard_Real aU1 = anUKnots->Value(i++);
140     Standard_Real aU2 = anUKnots->Value(i);
141
142     if(aU2 - aU1 < tolU) continue;
143
144     anIntegral.Perform(myUFunction, aU1, aU2, aNbPntsStart, aTol, aNbMaxIter);
145
146     if (!anIntegral.IsDone())
147       return Standard_False;
148
149     F         += anIntegral.Value();
150     aLocalErr += anIntegral.AbsolutError();
151     //cout << " TFunction : " << anIntegral.NbIterReached() << " " << anIntegral.AbsolutError() << endl;
152   }
153
154   F            *= aCoeff;
155   aLocalErr *= aAbsCoeff;
156   myAbsError = Max(myAbsError, aLocalErr);
157
158   myTolReached += aLocalErr;
159
160   if(Abs(F) > Epsilon(1.)) aLocalErr /= Abs(F);
161
162   myErrReached = Max(myErrReached, aLocalErr);
163
164
165   return Standard_True;
166 }
167
168 //=======================================================================
169 //function : GetStateNumber
170 //purpose  : 
171 //=======================================================================
172
173 Standard_Integer BRepGProp_TFunction::GetStateNumber()
174 {
175   //myErrReached  = myTolReached;
176   //myTolReached  = 0.;
177   //myNbPntOuter += RealToInt(0.5*myNbPntOuter);
178
179   //if (myNbPntOuter%2 == 0)
180   //myNbPntOuter++;
181
182   return 0;
183 }