b311480e |
1 | // Copyright (c) 1997-1999 Matra Datavision |
973c2be1 |
2 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
b311480e |
3 | // |
973c2be1 |
4 | // This file is part of Open CASCADE Technology software library. |
b311480e |
5 | // |
d5f74e42 |
6 | // This library is free software; you can redistribute it and/or modify it under |
7 | // the terms of the GNU Lesser General Public License version 2.1 as published |
973c2be1 |
8 | // by the Free Software Foundation, with special exception defined in the file |
9 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
10 | // distribution for complete text of the license and disclaimer of any warranty. |
b311480e |
11 | // |
973c2be1 |
12 | // Alternatively, this file may be used under the terms of Open CASCADE |
13 | // commercial license or contractual agreement. |
7fd59977 |
14 | |
15 | /* |
16 | Ce programme permet le calcul de l'integrale nieme d'une fonction. La |
17 | resolution peut etre vue comme l'integrale de l'integrale de ....(n fois), ce |
18 | qui necessite l'emploi de fonction recursive. |
19 | |
20 | Une integrale simple de Gauss est la somme de la fonction donnee sur les points |
21 | de Gauss affectee du poids de Gauss correspondant. |
22 | |
23 | I = somme(F(GaussPt) * GaussW) sur i = 1 jusqu'a l'ordre d'integration. |
24 | |
25 | Une integrale multiple est une somme recursive sur toutes les variables. |
26 | |
27 | Tout le calcul est effectue dans la methode recursive_iteration appelee n fois. |
28 | (n = nombre de variables) |
29 | |
30 | Cette methode fait une sommation sur toutes les variables, sur tous les ordres |
31 | d'integration correspondants de la valeur de la fonction aux points et poids |
32 | de Gauss. |
33 | |
34 | */ |
35 | |
0797d9d3 |
36 | //#ifndef OCCT_DEBUG |
7fd59977 |
37 | #define No_Standard_RangeError |
38 | #define No_Standard_OutOfRange |
39 | #define No_Standard_DimensionError |
42cf5bc1 |
40 | |
7fd59977 |
41 | //#endif |
42 | |
7fd59977 |
43 | #include <math.hxx> |
42cf5bc1 |
44 | #include <math_GaussMultipleIntegration.hxx> |
7fd59977 |
45 | #include <math_IntegerVector.hxx> |
42cf5bc1 |
46 | #include <math_Matrix.hxx> |
7fd59977 |
47 | #include <math_MultipleVarFunction.hxx> |
42cf5bc1 |
48 | #include <math_Vector.hxx> |
49 | #include <StdFail_NotDone.hxx> |
7fd59977 |
50 | |
51 | class IntegrationFunction { |
52 | |
53 | // Cette classe sert a conserver dans les champs les valeurs de tests de |
54 | // la fonction recursive.(En particulier le pointeur sur la fonction F, |
55 | // heritant d'une classe abstraite, non traite en cdl.) |
56 | |
57 | math_MultipleVarFunction *Fsav; |
58 | math_IntegerVector Ordsav; |
59 | Standard_Integer NVarsav; |
60 | math_Vector xr; |
61 | math_Vector xm; |
62 | math_Matrix GaussPoint; |
63 | math_Matrix GaussWeight; |
64 | Standard_Real Val; |
65 | Standard_Boolean Done; |
66 | |
67 | public: |
68 | IntegrationFunction(math_MultipleVarFunction& F, |
69 | const Standard_Integer maxsav, const Standard_Integer NVar, |
70 | const math_IntegerVector& Ord, |
71 | const math_Vector& Lowsav, |
72 | const math_Vector& Uppsav); |
73 | |
74 | Standard_Real Value() ; |
75 | Standard_Boolean IsDone() const; |
76 | Standard_Boolean recursive_iteration(Standard_Integer& n, math_IntegerVector& inc); |
77 | }; |
78 | |
79 | IntegrationFunction::IntegrationFunction(math_MultipleVarFunction& F, |
80 | const Standard_Integer maxsav, const Standard_Integer NVar, |
81 | const math_IntegerVector& Ord, |
82 | const math_Vector& Lowsav,const math_Vector& Uppsav): |
83 | Ordsav(1, NVar), |
84 | xr(1, NVar), |
85 | xm(1, NVar), |
86 | GaussPoint(1, NVar, 1, maxsav), |
87 | GaussWeight(1, NVar, 1, maxsav) |
88 | { |
89 | |
90 | Standard_Integer i, k; |
91 | math_IntegerVector inc(1, NVar); |
92 | inc.Init(0); |
93 | Fsav = &F; |
94 | NVarsav = NVar; |
95 | Ordsav = Ord; |
96 | Done = Standard_False; |
97 | |
98 | //Recuperation des points et poids de Gauss dans le fichier GaussPoints |
99 | for (i =1; i<= NVarsav; i++) { |
100 | xm(i) = 0.5*(Lowsav(i) + Uppsav(i)); |
101 | xr(i) = 0.5*(Uppsav(i) - Lowsav(i)); |
102 | math_Vector GP(1,Ordsav(i)), GW(1,Ordsav(i)); |
103 | math::GaussPoints(Ordsav(i),GP); |
104 | math::GaussWeights(Ordsav(i),GW); |
105 | for (k = 1; k <= Ordsav(i); k++) { |
106 | GaussPoint(i,k) = GP(k); //kieme point et poids de |
107 | GaussWeight(i,k) = GW(k);//Gauss de la variable i. |
108 | } |
109 | } |
110 | Val = 0.0; |
111 | Standard_Integer Iterdeb = 1; |
112 | Standard_Boolean recur = recursive_iteration(Iterdeb, inc); |
113 | if (recur) { |
114 | //On ramene l'integration a la bonne echelle. |
115 | for ( i = 1; i <= NVarsav; i++) { |
116 | Val *= xr(i); |
117 | } |
118 | Done = Standard_True; |
119 | } |
120 | } |
121 | |
122 | Standard_Real IntegrationFunction::Value() { |
123 | return Val; |
124 | } |
125 | |
126 | Standard_Boolean IntegrationFunction::IsDone() const{ |
127 | return Done; |
128 | } |
129 | |
130 | Standard_Boolean IntegrationFunction::recursive_iteration(Standard_Integer& n, |
131 | math_IntegerVector& inc) { |
132 | |
133 | // Termination criterium : |
134 | // Calcul de la valeur de la fonction aux points de Gauss fixes: |
135 | int local; |
136 | if (n == (NVarsav+1)) { |
137 | math_Vector dx(1, NVarsav); |
138 | Standard_Integer j ; |
139 | for ( j = 1; j <= NVarsav; j++) { |
140 | dx(j) = xr(j)* GaussPoint(j, inc(j)); |
141 | } |
142 | Standard_Real F1; |
143 | Standard_Boolean Ok = Fsav->Value(xm + dx, F1); |
144 | if (!Ok) {return Standard_False;}; |
145 | Standard_Real Interm = 1; |
146 | for ( j = 1; j <= NVarsav; j++) { |
147 | Interm *= GaussWeight(j, inc(j)); |
148 | } |
149 | Val += Interm*F1; |
150 | return Standard_True; |
151 | } |
152 | |
153 | // Calcul recursif, pour chaque variable n de la valeur de la fonction aux |
154 | // Ordsav(n) points de Gauss. |
155 | |
156 | Standard_Boolean OK=Standard_False; |
157 | for (inc(n) = 1; inc(n) <= Ordsav(n); inc(n)++) { |
158 | local = n + 1; |
159 | OK = recursive_iteration(local, inc); |
160 | } |
161 | return OK; |
162 | } |
163 | |
164 | math_GaussMultipleIntegration:: |
165 | math_GaussMultipleIntegration(math_MultipleVarFunction& F, |
166 | const math_Vector& Lower, |
167 | const math_Vector& Upper, |
168 | const math_IntegerVector& Order) |
169 | { |
170 | Standard_Integer MaxOrder = math::GaussPointsMax(); |
171 | |
172 | Standard_Integer i, max =0; |
173 | Standard_Integer NVar = F.NbVariables(); |
174 | math_IntegerVector Ord(1, NVar); |
175 | math_Vector Lowsav(1, NVar); |
176 | math_Vector Uppsav(1, NVar); |
177 | Lowsav = Lower; |
178 | Uppsav = Upper; |
179 | // Ord = Order; |
180 | Done = Standard_False; |
181 | for (i = 1; i <= NVar; i++) { |
182 | if (Order(i) > MaxOrder) { |
183 | Ord(i) = MaxOrder; |
184 | } |
185 | else { |
186 | Ord(i) = Order(i); |
187 | } |
188 | if(Ord(i) >= max) {max = Ord(i);} |
189 | } |
190 | |
191 | //Calcul de l'integrale aux points de Gauss. |
192 | // Il s agit d une somme multiple sur le domaine [Lower, Upper]. |
193 | |
194 | IntegrationFunction Func(F, max, NVar, Ord, Lowsav, Uppsav); |
195 | if (Func.IsDone()) { |
196 | Val = Func.Value(); |
197 | Done = Standard_True; |
198 | } |
199 | } |
200 | |
201 | void math_GaussMultipleIntegration::Dump(Standard_OStream& o) const { |
202 | |
203 | o <<"math_GaussMultipleIntegration "; |
204 | if (Done) { |
205 | o << " Status = Done \n"; |
206 | o << " Integration value = " << Val <<"\n"; |
207 | } |
208 | else { |
209 | o << "Status = not Done \n"; |
210 | } |
211 | } |
212 | |
213 | |
214 | |