src/math/math_GaussMultipleIntegration.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 /*
  16 Ce programme permet le calcul de l'integrale nieme d'une fonction. La
  17 resolution peut etre vue comme l'integrale de l'integrale de ....(n fois), ce
  18 qui necessite l'emploi de fonction recursive.
  19
  20 Une integrale simple de Gauss est la somme de la fonction donnee sur les points
  21 de Gauss affectee du poids de Gauss correspondant.
  22
  23 I = somme(F(GaussPt) * GaussW)   sur i = 1 jusqu'a l'ordre d'integration.
  24
  25 Une integrale multiple est une somme recursive sur toutes les variables.
  26
  27 Tout le calcul est effectue dans la methode recursive_iteration appelee n fois.
  28  (n = nombre de variables)
  29
  30 Cette methode fait une sommation sur toutes les variables, sur tous les ordres
  31 d'integration correspondants de la valeur de la fonction aux points et poids
  32 de Gauss.
  33
  34 */
  35
  36 //#ifndef OCCT_DEBUG
  37 #define No_Standard_RangeError
  38 #define No_Standard_OutOfRange
  39 #define No_Standard_DimensionError
  40
  41 //#endif
  42
  43 #include <math.hxx>
  44 #include <math_GaussMultipleIntegration.hxx>
  45 #include <math_IntegerVector.hxx>
  46 #include <math_Matrix.hxx>
  47 #include <math_MultipleVarFunction.hxx>
  48 #include <math_Vector.hxx>
  49 #include <StdFail_NotDone.hxx>
  50
  51 class IntegrationFunction {
  52
  53 // Cette classe sert a conserver dans les champs les valeurs de tests de
  54 // la fonction recursive.(En particulier le pointeur sur la fonction F,
  55 // heritant d'une classe abstraite, non traite en cdl.)
  56
  57   math_MultipleVarFunction *Fsav;
  58   math_IntegerVector Ordsav;
  59   Standard_Integer NVarsav;
  60   math_Vector xr;
  61   math_Vector xm;
  62   math_Matrix GaussPoint;
  63   math_Matrix GaussWeight;
  64   Standard_Real Val;
  65   Standard_Boolean Done;
  66
  67  public:
  68     IntegrationFunction(math_MultipleVarFunction& F,
  69                         const Standard_Integer maxsav, const Standard_Integer NVar,
  70                         const math_IntegerVector& Ord,
  71                         const math_Vector& Lowsav,
  72                         const math_Vector& Uppsav);
  73
  74     Standard_Real Value() ;
  75     Standard_Boolean IsDone() const;
  76     Standard_Boolean recursive_iteration(Standard_Integer& n, math_IntegerVector& inc);
  77 };
  78
  79  IntegrationFunction::IntegrationFunction(math_MultipleVarFunction& F,
  80                         const Standard_Integer maxsav, const Standard_Integer NVar,
  81                         const math_IntegerVector& Ord,
  82                         const math_Vector& Lowsav,const  math_Vector& Uppsav):
  83                         Ordsav(1, NVar),
  84                         xr(1, NVar),
  85                         xm(1, NVar),
  86                         GaussPoint(1, NVar, 1, maxsav),
  87                         GaussWeight(1, NVar, 1, maxsav)
  88 {
  89
  90 Standard_Integer i, k;
  91 math_IntegerVector inc(1, NVar);
  92 inc.Init(0);
  93 Fsav = &F;
  94 NVarsav = NVar;
  95 Ordsav = Ord;
  96 Done = Standard_False;
  97
  98 //Recuperation des points et poids de Gauss dans le fichier GaussPoints
  99   for (i =1; i<= NVarsav; i++) {
 100     xm(i) = 0.5*(Lowsav(i) + Uppsav(i));
 101     xr(i) = 0.5*(Uppsav(i) - Lowsav(i));
 102     math_Vector GP(1,Ordsav(i)), GW(1,Ordsav(i));
 103     math::GaussPoints(Ordsav(i),GP);
 104     math::GaussWeights(Ordsav(i),GW);
 105     for (k = 1; k <= Ordsav(i); k++) {
 106       GaussPoint(i,k)  = GP(k);  //kieme point et poids de
 107       GaussWeight(i,k) = GW(k);//Gauss de la variable i.
 108     }
 109   }
 110   Val = 0.0;
 111   Standard_Integer Iterdeb = 1;
 112   Standard_Boolean recur = recursive_iteration(Iterdeb, inc);
 113   if (recur) {
 114 //On ramene l'integration a la bonne echelle.
 115      for ( i = 1; i <= NVarsav; i++) {
 116        Val *= xr(i);
 117      }
 118      Done = Standard_True;
 119   }
 120 }
 121
 122 Standard_Real IntegrationFunction::Value() {
 123   return Val;
 124 }
 125
 126 Standard_Boolean IntegrationFunction::IsDone() const{
 127   return Done;
 128 }
 129
 130 Standard_Boolean IntegrationFunction::recursive_iteration(Standard_Integer& n,
 131                                    math_IntegerVector& inc) {
 132
 133 // Termination criterium :
 134 // Calcul de la valeur de la fonction aux points de Gauss fixes:
 135         int local;
 136    if (n == (NVarsav+1)) {
 137      math_Vector dx(1, NVarsav);
 138      Standard_Integer j ;
 139      for ( j = 1; j <= NVarsav; j++) {
 140        dx(j) = xr(j)* GaussPoint(j, inc(j));
 141      }
 142      Standard_Real F1;
 143      Standard_Boolean Ok = Fsav->Value(xm + dx, F1);
 144      if (!Ok) {return Standard_False;};
 145      Standard_Real Interm = 1;
 146      for ( j = 1; j <= NVarsav; j++) {
 147        Interm *= GaussWeight(j, inc(j));
 148      }
 149      Val += Interm*F1;
 150      return Standard_True;
 151    }
 152
 153    // Calcul recursif, pour chaque variable n de la valeur de la fonction aux
 154    // Ordsav(n) points de Gauss.
 155
 156    Standard_Boolean OK=Standard_False;
 157    for (inc(n) = 1; inc(n) <= Ordsav(n); inc(n)++) {
 158         local = n + 1;
 159      OK = recursive_iteration(local, inc);
 160    }
 161    return OK;
 162 }
 163
 164 math_GaussMultipleIntegration::
 165                 math_GaussMultipleIntegration(math_MultipleVarFunction& F,
 166                                               const math_Vector& Lower,
 167                                               const math_Vector& Upper,
 168                                               const math_IntegerVector& Order)
 169 {
 170   Standard_Integer MaxOrder = math::GaussPointsMax();
 171
 172   Standard_Integer i,  max =0;
 173   Standard_Integer NVar = F.NbVariables();
 174   math_IntegerVector Ord(1, NVar);
 175   math_Vector Lowsav(1, NVar);
 176   math_Vector Uppsav(1, NVar);
 177   Lowsav = Lower;
 178   Uppsav = Upper;
 179 //  Ord = Order;
 180   Done = Standard_False;
 181   for (i = 1; i <= NVar; i++) {
 182     if (Order(i) > MaxOrder) {
 183       Ord(i) = MaxOrder;
 184     }
 185     else {
 186       Ord(i) = Order(i);
 187     }
 188     if(Ord(i) >= max) {max = Ord(i);}
 189   }
 190
 191 //Calcul de l'integrale aux points de Gauss.
 192 // Il s agit d une somme multiple sur le domaine [Lower, Upper].
 193
 194   IntegrationFunction Func(F, max, NVar, Ord, Lowsav, Uppsav);
 195   if (Func.IsDone()) {
 196     Val = Func.Value();
 197     Done = Standard_True;
 198   }
 199 }
 200
 201 void math_GaussMultipleIntegration::Dump(Standard_OStream& o) const {
 202
 203   o <<"math_GaussMultipleIntegration ";
 204    if (Done) {
 205      o << " Status = Done \n";
 206      o << " Integration value = " << Val <<"\n";
 207    }
 208    else {
 209      o << "Status = not Done \n";
 210    }
 211 }
 212
 213
 214