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