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