0022904: Clean up sccsid variables
[occt.git] / src / math / math_GaussMultipleIntegration.cxx
1 // File math_GaussMultipleIntegration.cxx
2
3
4 /*
5 Ce programme permet le calcul de l'integrale nieme d'une fonction. La 
6 resolution peut etre vue comme l'integrale de l'integrale de ....(n fois), ce
7 qui necessite l'emploi de fonction recursive.
8
9 Une integrale simple de Gauss est la somme de la fonction donnee sur les points
10 de Gauss affectee du poids de Gauss correspondant.
11
12 I = somme(F(GaussPt) * GaussW)   sur i = 1 jusqu'a l'ordre d'integration.
13
14 Une integrale multiple est une somme recursive sur toutes les variables.
15
16 Tout le calcul est effectue dans la methode recursive_iteration appelee n fois.
17  (n = nombre de variables)
18
19 Cette methode fait une sommation sur toutes les variables, sur tous les ordres
20 d'integration correspondants de la valeur de la fonction aux points et poids
21 de 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
38 class 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
77 Standard_Integer i, k;
78 math_IntegerVector inc(1, NVar);
79 inc.Init(0);
80 Fsav = &F;
81 NVarsav = NVar;
82 Ordsav = Ord;
83 Done = 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
109 Standard_Real IntegrationFunction::Value() {
110   return Val;
111 }
112
113 Standard_Boolean IntegrationFunction::IsDone() const{
114   return Done;
115 }
116
117 Standard_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
151 math_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
188 void 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