0022904: Clean up sccsid variables
[occt.git] / src / math / math_GaussMultipleIntegration.cxx
CommitLineData
7fd59977 1// File math_GaussMultipleIntegration.cxx
2
3
4/*
5Ce programme permet le calcul de l'integrale nieme d'une fonction. La
6resolution peut etre vue comme l'integrale de l'integrale de ....(n fois), ce
7qui necessite l'emploi de fonction recursive.
8
9Une integrale simple de Gauss est la somme de la fonction donnee sur les points
10de Gauss affectee du poids de Gauss correspondant.
11
12I = somme(F(GaussPt) * GaussW) sur i = 1 jusqu'a l'ordre d'integration.
13
14Une integrale multiple est une somme recursive sur toutes les variables.
15
16Tout le calcul est effectue dans la methode recursive_iteration appelee n fois.
17 (n = nombre de variables)
18
19Cette methode fait une sommation sur toutes les variables, sur tous les ordres
20d'integration correspondants de la valeur de la fonction aux points et poids
21de Gauss.
22
23*/
24
25//#ifndef DEB
26#define No_Standard_RangeError
27#define No_Standard_OutOfRange
28#define No_Standard_DimensionError
29//#endif
30
31#include <math_GaussMultipleIntegration.ixx>
32#include <math.hxx>
33#include <math_Matrix.hxx>
34#include <math_Vector.hxx>
35#include <math_IntegerVector.hxx>
36#include <math_MultipleVarFunction.hxx>
37
38class IntegrationFunction {
39
40// Cette classe sert a conserver dans les champs les valeurs de tests de
41// la fonction recursive.(En particulier le pointeur sur la fonction F,
42// heritant d'une classe abstraite, non traite en cdl.)
43
44 math_MultipleVarFunction *Fsav;
45 math_IntegerVector Ordsav;
46 Standard_Integer NVarsav;
47 math_Vector xr;
48 math_Vector xm;
49 math_Matrix GaussPoint;
50 math_Matrix GaussWeight;
51 Standard_Real Val;
52 Standard_Boolean Done;
53
54 public:
55 IntegrationFunction(math_MultipleVarFunction& F,
56 const Standard_Integer maxsav, const Standard_Integer NVar,
57 const math_IntegerVector& Ord,
58 const math_Vector& Lowsav,
59 const math_Vector& Uppsav);
60
61 Standard_Real Value() ;
62 Standard_Boolean IsDone() const;
63 Standard_Boolean recursive_iteration(Standard_Integer& n, math_IntegerVector& inc);
64};
65
66 IntegrationFunction::IntegrationFunction(math_MultipleVarFunction& F,
67 const Standard_Integer maxsav, const Standard_Integer NVar,
68 const math_IntegerVector& Ord,
69 const math_Vector& Lowsav,const math_Vector& Uppsav):
70 Ordsav(1, NVar),
71 xr(1, NVar),
72 xm(1, NVar),
73 GaussPoint(1, NVar, 1, maxsav),
74 GaussWeight(1, NVar, 1, maxsav)
75{
76
77Standard_Integer i, k;
78math_IntegerVector inc(1, NVar);
79inc.Init(0);
80Fsav = &F;
81NVarsav = NVar;
82Ordsav = Ord;
83Done = Standard_False;
84
85//Recuperation des points et poids de Gauss dans le fichier GaussPoints
86 for (i =1; i<= NVarsav; i++) {
87 xm(i) = 0.5*(Lowsav(i) + Uppsav(i));
88 xr(i) = 0.5*(Uppsav(i) - Lowsav(i));
89 math_Vector GP(1,Ordsav(i)), GW(1,Ordsav(i));
90 math::GaussPoints(Ordsav(i),GP);
91 math::GaussWeights(Ordsav(i),GW);
92 for (k = 1; k <= Ordsav(i); k++) {
93 GaussPoint(i,k) = GP(k); //kieme point et poids de
94 GaussWeight(i,k) = GW(k);//Gauss de la variable i.
95 }
96 }
97 Val = 0.0;
98 Standard_Integer Iterdeb = 1;
99 Standard_Boolean recur = recursive_iteration(Iterdeb, inc);
100 if (recur) {
101//On ramene l'integration a la bonne echelle.
102 for ( i = 1; i <= NVarsav; i++) {
103 Val *= xr(i);
104 }
105 Done = Standard_True;
106 }
107}
108
109Standard_Real IntegrationFunction::Value() {
110 return Val;
111}
112
113Standard_Boolean IntegrationFunction::IsDone() const{
114 return Done;
115}
116
117Standard_Boolean IntegrationFunction::recursive_iteration(Standard_Integer& n,
118 math_IntegerVector& inc) {
119
120// Termination criterium :
121// Calcul de la valeur de la fonction aux points de Gauss fixes:
122 int local;
123 if (n == (NVarsav+1)) {
124 math_Vector dx(1, NVarsav);
125 Standard_Integer j ;
126 for ( j = 1; j <= NVarsav; j++) {
127 dx(j) = xr(j)* GaussPoint(j, inc(j));
128 }
129 Standard_Real F1;
130 Standard_Boolean Ok = Fsav->Value(xm + dx, F1);
131 if (!Ok) {return Standard_False;};
132 Standard_Real Interm = 1;
133 for ( j = 1; j <= NVarsav; j++) {
134 Interm *= GaussWeight(j, inc(j));
135 }
136 Val += Interm*F1;
137 return Standard_True;
138 }
139
140 // Calcul recursif, pour chaque variable n de la valeur de la fonction aux
141 // Ordsav(n) points de Gauss.
142
143 Standard_Boolean OK=Standard_False;
144 for (inc(n) = 1; inc(n) <= Ordsav(n); inc(n)++) {
145 local = n + 1;
146 OK = recursive_iteration(local, inc);
147 }
148 return OK;
149}
150
151math_GaussMultipleIntegration::
152 math_GaussMultipleIntegration(math_MultipleVarFunction& F,
153 const math_Vector& Lower,
154 const math_Vector& Upper,
155 const math_IntegerVector& Order)
156{
157 Standard_Integer MaxOrder = math::GaussPointsMax();
158
159 Standard_Integer i, max =0;
160 Standard_Integer NVar = F.NbVariables();
161 math_IntegerVector Ord(1, NVar);
162 math_Vector Lowsav(1, NVar);
163 math_Vector Uppsav(1, NVar);
164 Lowsav = Lower;
165 Uppsav = Upper;
166// Ord = Order;
167 Done = Standard_False;
168 for (i = 1; i <= NVar; i++) {
169 if (Order(i) > MaxOrder) {
170 Ord(i) = MaxOrder;
171 }
172 else {
173 Ord(i) = Order(i);
174 }
175 if(Ord(i) >= max) {max = Ord(i);}
176 }
177
178//Calcul de l'integrale aux points de Gauss.
179// Il s agit d une somme multiple sur le domaine [Lower, Upper].
180
181 IntegrationFunction Func(F, max, NVar, Ord, Lowsav, Uppsav);
182 if (Func.IsDone()) {
183 Val = Func.Value();
184 Done = Standard_True;
185 }
186}
187
188void math_GaussMultipleIntegration::Dump(Standard_OStream& o) const {
189
190 o <<"math_GaussMultipleIntegration ";
191 if (Done) {
192 o << " Status = Done \n";
193 o << " Integration value = " << Val <<"\n";
194 }
195 else {
196 o << "Status = not Done \n";
197 }
198}
199
200
201