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