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 | Par Gauss le calcul d'une integrale simple se transforme en sommation des |
17 | valeurs de la fonction donnee aux <Order> points de Gauss affectee des poids |
18 | de Gauss. |
19 | |
20 | Les points et poids de Gauss sont stockes dans GaussPoints.cxx. |
21 | Les points sont compris entre les valeurs -1 et +1, ce qui necessite un |
22 | changement de variable pour les faire varier dans l'intervalle [Lower, Upper]. |
23 | |
24 | |
25 | On veut calculer Integrale( f(u)* du) entre a et b. |
26 | |
27 | |
28 | Etapes du calcul: |
29 | |
30 | 1- calcul de la fonction au ieme point de Gauss (apres changement de variable). |
31 | |
32 | 2- multiplication de cette valeur par le ieme poids de Gauss. |
33 | |
34 | 3- sommation de toutes ces valeurs. |
35 | |
36 | 4- retour a l'intervalle [Lower, Upper] de notre integrale. |
37 | |
38 | */ |
39 | |
0797d9d3 |
40 | //#ifndef OCCT_DEBUG |
7fd59977 |
41 | #define No_Standard_RangeError |
42 | #define No_Standard_OutOfRange |
43 | #define No_Standard_DimensionError |
7fd59977 |
44 | |
42cf5bc1 |
45 | //#endif |
7fd59977 |
46 | |
47 | #include <math.hxx> |
7fd59977 |
48 | #include <math_Function.hxx> |
42cf5bc1 |
49 | #include <math_GaussSingleIntegration.hxx> |
50 | #include <math_Vector.hxx> |
51 | #include <StdFail_NotDone.hxx> |
7fd59977 |
52 | |
53 | math_GaussSingleIntegration::math_GaussSingleIntegration() : Done(Standard_False) |
54 | { |
55 | } |
56 | |
57 | math_GaussSingleIntegration:: |
58 | math_GaussSingleIntegration(math_Function& F, |
59 | const Standard_Real Lower, |
60 | const Standard_Real Upper, |
61 | const Standard_Integer Order) |
62 | { |
63 | Standard_Integer theOrder = Min(math::GaussPointsMax(), Order); |
64 | Perform(F, Lower, Upper, theOrder); |
65 | } |
66 | |
67 | math_GaussSingleIntegration:: |
68 | math_GaussSingleIntegration(math_Function& F, |
69 | const Standard_Real Lower, |
70 | const Standard_Real Upper, |
71 | const Standard_Integer Order, |
72 | const Standard_Real Tol) |
73 | { |
74 | Standard_Integer theOrder = Min(math::GaussPointsMax(), Order); |
75 | |
76 | const Standard_Integer IterMax = 13; // Max number of iteration |
77 | Standard_Integer NIter = 1; // current number of iteration |
78 | Standard_Integer NbInterval = 1; // current number of subintervals |
79 | Standard_Real dU,OldLen,Len; |
80 | |
81 | Perform(F, Lower, Upper, theOrder); |
82 | Len = Val; |
83 | do { |
84 | OldLen = Len; |
85 | Len = 0.; |
86 | NbInterval *= 2; |
87 | dU = (Upper-Lower)/NbInterval; |
88 | for (Standard_Integer i=1; i<=NbInterval; i++) { |
89 | Perform(F, Lower+(i-1)*dU, Lower+i*dU, theOrder); |
90 | if (!Done) return; |
91 | Len += Val; |
92 | } |
93 | NIter++; |
94 | } |
95 | while (fabs(OldLen-Len) > Tol && NIter <= IterMax); |
96 | |
97 | Val = Len; |
98 | } |
99 | |
100 | void math_GaussSingleIntegration::Perform(math_Function& F, |
101 | const Standard_Real Lower, |
102 | const Standard_Real Upper, |
103 | const Standard_Integer Order) |
104 | { |
105 | Standard_Real xr, xm, dx; |
106 | Standard_Integer j; |
107 | Standard_Real F1, F2; |
108 | Standard_Boolean Ok1; |
109 | math_Vector GaussP(1, Order); |
110 | math_Vector GaussW(1, Order); |
111 | Done = Standard_False; |
112 | |
113 | //Recuperation des points de Gauss dans le fichier GaussPoints. |
114 | math::GaussPoints(Order,GaussP); |
115 | math::GaussWeights(Order,GaussW); |
116 | |
117 | // Calcul de l'integrale aux points de Gauss : |
118 | |
119 | // Changement de variable pour la mise a l'echelle [Lower, Upper] : |
120 | xm = 0.5*(Upper + Lower); |
121 | xr = 0.5*(Upper - Lower); |
122 | Val = 0.; |
123 | |
124 | Standard_Integer ind = Order/2, ind1 = (Order+1)/2; |
125 | if(ind1 > ind) { // odder case |
126 | Ok1 = F.Value(xm, Val); |
127 | if (!Ok1) return; |
128 | Val *= GaussW(ind1); |
129 | } |
130 | // Sommation sur tous les points de Gauss: avec utilisation de la symetrie. |
131 | for (j = 1; j <= ind; j++) { |
132 | dx = xr*GaussP(j); |
133 | Ok1 = F.Value(xm-dx, F1); |
134 | if(!Ok1) return; |
135 | Ok1 = F.Value(xm+dx, F2); |
136 | if(!Ok1) return; |
137 | // Multiplication par les poids de Gauss. |
138 | Standard_Real FT = F1+F2; |
139 | Val += GaussW(j)*FT; |
140 | } |
141 | // Mise a l'echelle de l'intervalle [Lower, Upper] |
142 | Val *= xr; |
143 | Done = Standard_True; |
144 | } |
145 | |
146 | void math_GaussSingleIntegration::Dump(Standard_OStream& o) const { |
147 | |
148 | o <<"math_GaussSingleIntegration "; |
149 | if (Done) { |
150 | o << " Status = Done \n"; |
151 | o << "Integration Value = " << Val<<"\n"; |
152 | } |
153 | else { |
154 | o << "Status = not Done \n"; |
155 | } |
156 | } |