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. |
b311480e |
14 | |
7fd59977 |
15 | // lpa le 8/08/91 |
16 | |
17 | // Ce programme utilise l algorithme d Uzawa pour resoudre un systeme |
18 | // dans le cas de contraintes. Si ce sont des contraintes d egalite, la |
19 | // resolution est directe. |
20 | // Le programme ci-dessous utilise la methode de Crout pour trouver |
21 | // l inverse d une matrice symetrique (Le gain est d environ 30% par |
22 | // rapport a Gauss.). Les calculs sur les matrices sont faits avec chaque |
23 | // coordonnee car il est plus long d utiliser les methodes deja ecrites |
24 | // de la classe Matrix avec un passage par valeur. |
25 | |
0797d9d3 |
26 | //#ifndef OCCT_DEBUG |
7fd59977 |
27 | #define No_Standard_RangeError |
28 | #define No_Standard_OutOfRange |
29 | #define No_Standard_DimensionError |
30 | //#endif |
31 | |
32 | #include <math_Uzawa.ixx> |
33 | #include <math_Crout.hxx> |
34 | #include <Standard_DimensionError.hxx> |
35 | |
36 | |
37 | math_Uzawa::math_Uzawa(const math_Matrix& Cont, const math_Vector& Secont, |
38 | const math_Vector& StartingPoint, |
39 | const Standard_Real EpsLix, const Standard_Real EpsLic, |
40 | const Standard_Integer NbIterations) : |
41 | Resul(1, Cont.ColNumber()), |
42 | Erruza(1, Cont.ColNumber()), |
43 | Errinit(1, Cont.ColNumber()), |
44 | Vardua(1, Cont.RowNumber()), |
45 | CTCinv(1, Cont.RowNumber(), |
46 | 1, Cont.RowNumber()){ |
47 | |
48 | Perform(Cont, Secont, StartingPoint, Cont.RowNumber(), 0, EpsLix, |
49 | EpsLic, NbIterations); |
50 | } |
51 | |
52 | math_Uzawa::math_Uzawa(const math_Matrix& Cont, const math_Vector& Secont, |
53 | const math_Vector& StartingPoint, |
54 | const Standard_Integer Nce, const Standard_Integer Nci, |
55 | const Standard_Real EpsLix, const Standard_Real EpsLic, |
56 | const Standard_Integer NbIterations) : |
57 | Resul(1, Cont.ColNumber()), |
58 | Erruza(1, Cont.ColNumber()), |
59 | Errinit(1, Cont.ColNumber()), |
60 | Vardua(1, Cont.RowNumber()), |
61 | CTCinv(1, Cont.RowNumber(), |
62 | 1, Cont.RowNumber()) { |
63 | |
64 | Perform(Cont, Secont, StartingPoint, Nce, Nci, EpsLix, EpsLic, NbIterations); |
65 | |
66 | } |
67 | |
68 | |
69 | void math_Uzawa::Perform(const math_Matrix& Cont, const math_Vector& Secont, |
70 | const math_Vector& StartingPoint, |
71 | const Standard_Integer Nce, const Standard_Integer Nci, |
72 | const Standard_Real EpsLix, const Standard_Real EpsLic, |
73 | const Standard_Integer NbIterations) { |
74 | |
75 | Standard_Integer i,j, k; |
76 | Standard_Real scale; |
77 | Standard_Real Normat, Normli, Xian, Xmax=0, Xmuian; |
78 | Standard_Real Rho, Err, Err1, ErrMax=0, Coef = 1./Sqrt(2.); |
79 | Standard_Integer Nlig = Cont.RowNumber(); |
80 | Standard_Integer Ncol = Cont.ColNumber(); |
81 | |
82 | Standard_DimensionError_Raise_if((Secont.Length() != Nlig) || |
83 | ((Nce+Nci) != Nlig), " "); |
84 | |
85 | // Calcul du vecteur Cont*X0 - D: (erreur initiale) |
86 | //================================================== |
87 | |
88 | for (i = 1; i<= Nlig; i++) { |
89 | Errinit(i) = Cont(i,1)*StartingPoint(1)-Secont(i); |
90 | for (j = 2; j <= Ncol; j++) { |
91 | Errinit(i) += Cont(i,j)*StartingPoint(j); |
92 | } |
93 | } |
94 | |
95 | if (Nci == 0) { // cas de resolution directe |
96 | NbIter = 1; //========================== |
97 | // Calcul de Cont*T(Cont) |
98 | for (i = 1; i <= Nlig; i++) { |
99 | for (j =1; j <= i; j++) { // a utiliser avec Crout. |
100 | // for (j = 1; j <= Nlig; j++) { // a utiliser pour Gauss. |
101 | CTCinv(i,j)= Cont(i,1)*Cont(j,1); |
102 | for (k = 2; k <= Ncol; k++) { |
103 | CTCinv(i,j) += Cont(i,k)*Cont(j,k); |
104 | } |
105 | } |
106 | } |
107 | // Calcul de l inverse de CTCinv : |
108 | //================================ |
109 | // CTCinv = CTCinv.Inverse(); // utilisation de Gauss. |
110 | math_Crout inv(CTCinv); // utilisation de Crout. |
111 | CTCinv = inv.Inverse(); |
112 | for (i =1; i <= Nlig; i++) { |
113 | scale = CTCinv(i,1)*Errinit(1); |
114 | for (j = 2; j <= i; j++) { |
115 | scale += CTCinv(i,j)*Errinit(j); |
116 | } |
117 | for (j =i+1; j <= Nlig; j++) { |
118 | scale += CTCinv(j,i)*Errinit(j); |
119 | } |
120 | Vardua(i) = scale; |
121 | } |
122 | for (i = 1; i <= Ncol; i++) { |
123 | Erruza(i) = -Cont(1,i)*Vardua(1); |
124 | for (j = 2; j <= Nlig; j++) { |
125 | Erruza(i) -= Cont(j,i)*Vardua(j); |
126 | } |
127 | } |
128 | // restitution des valeurs calculees: |
129 | //=================================== |
130 | Resul = StartingPoint + Erruza; |
131 | Done = Standard_True; |
132 | return; |
133 | } // Fin de la resolution directe. |
134 | //============================== |
135 | |
136 | else { |
137 | // Initialisation des variables duales. |
138 | //===================================== |
139 | for (i = 1; i <= Nlig; i++) { |
140 | if (i <= Nce) { |
141 | Vardua(i) = 0.0; |
142 | } |
143 | else { |
144 | Vardua(i) = 1.; |
145 | } |
146 | } |
147 | |
148 | // Calcul du coefficient Rho: |
149 | //=========================== |
150 | Normat = 0.0; |
151 | for (i = 1; i <= Nlig; i++) { |
152 | Normli = Cont(i,1)*Cont(i,1); |
153 | for (j = 2; j <= Ncol; j++) { |
154 | Normli += Cont(i,j)*Cont(i,j); |
155 | } |
156 | Normat += Normli; |
157 | } |
158 | Rho = Coef/Normat; |
159 | |
160 | // Boucle des iterations de la methode d Uzawa. |
161 | //============================================= |
162 | for (NbIter = 1; NbIter <= NbIterations; NbIter++) { |
163 | for (i = 1; i <= Ncol; i++) { |
164 | Xian = Erruza(i); |
165 | Erruza(i) = -Cont(1,i)*Vardua(1); |
166 | for(j =2; j <= Nlig; j++) { |
167 | Erruza(i) -= Cont(j,i)*Vardua(j); |
168 | } |
169 | if (NbIter > 1) { |
170 | if (i == 1) { |
171 | Xmax = Abs(Erruza(i) - Xian); |
172 | } |
173 | Xmax = Max(Xmax, Abs(Erruza(i)-Xian)); |
174 | } |
175 | } |
176 | |
177 | // Calcul de Xmu a l iteration NbIter et evaluation de l erreur sur |
178 | // la verification des contraintes. |
179 | //================================================================= |
180 | for (i = 1; i <= Nlig; i++) { |
181 | Err = Cont(i,1)*Erruza(1) + Errinit(i); |
182 | for (j = 2; j <= Ncol; j++) { |
183 | Err += Cont(i,j)*Erruza(j); |
184 | } |
185 | if (i <= Nce) { |
186 | Vardua(i) += Rho * Err; |
187 | Err1 = Abs(Rho*Err); |
188 | } |
189 | else { |
190 | Xmuian = Vardua(i); |
191 | Vardua(i) = Max(0.0, Vardua(i)+ Rho*Err); |
192 | Err1 = Abs(Vardua(i) - Xmuian); |
193 | } |
194 | if (i == 1) { |
195 | ErrMax = Err1; |
196 | } |
197 | ErrMax = Max(ErrMax, Err1); |
198 | } |
199 | |
200 | if (NbIter > 1) { |
201 | if (Xmax <= EpsLix) { |
202 | if (ErrMax <= EpsLic) { |
203 | // cout <<"Convergence atteinte dans Uzawa"<<endl; |
204 | Done = Standard_True; |
205 | } |
206 | else { |
207 | // cout <<"convergence non atteinte pour le probleme dual"<<endl; |
208 | Done = Standard_False; |
209 | return; |
210 | } |
211 | // Restitution des valeurs calculees |
212 | //================================== |
213 | Resul = StartingPoint + Erruza; |
214 | Done = Standard_True; |
215 | return; |
216 | } |
217 | } |
218 | |
219 | } // fin de la boucle d iterations. |
220 | Done = Standard_False; |
221 | } |
222 | } |
223 | |
224 | |
225 | void math_Uzawa::Duale(math_Vector& V) const |
226 | { |
227 | V = Vardua; |
228 | } |
229 | |
230 | void math_Uzawa::Dump(Standard_OStream& o) const { |
231 | |
232 | o << "math_Uzawa"; |
233 | if(Done) { |
234 | o << " Status = Done \n"; |
235 | o << " Number of iterations = " << NbIter << endl; |
236 | o << " The solution vector is: " << Resul << endl; |
237 | } |
238 | else { |
239 | o << " Status = not Done \n"; |
240 | } |
241 | } |
242 | |
243 | |
244 | |
245 | |
246 | |