src/math/math_Uzawa.cxx

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