src/math/math_GaussSetIntegration.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 //#ifndef OCCT_DEBUG
  16 #define No_Standard_RangeError
  17 #define No_Standard_OutOfRange
  18 #define No_Standard_DimensionError
  19
  20 //#endif
  21
  22 #include <math.hxx>
  23 #include <math_FunctionSet.hxx>
  24 #include <math_GaussSetIntegration.hxx>
  25 #include <math_Vector.hxx>
  26 #include <Standard_NotImplemented.hxx>
  27 #include <StdFail_NotDone.hxx>
  28
  29 math_GaussSetIntegration::math_GaussSetIntegration(math_FunctionSet& F,
  30                                                     const math_Vector& Lower,
  31                                                     const math_Vector& Upper,
  32                                                     const math_IntegerVector& Order)
  33                           : Val(1, F.NbEquations()) {
  34
  35    Standard_Integer NbEqua = F.NbEquations() , NbVar = F.NbVariables();
  36    Standard_Integer i;
  37    Standard_Boolean IsOk;
  38    math_Vector FVal1(1, NbEqua), FVal2(1, NbEqua), Tval(1, NbVar);
  39
  40
  41 // Verification
  42    Standard_NotImplemented_Raise_if(
  43             NbVar != 1 || Order.Value(Order.Lower()) > math::GaussPointsMax(),
  44             "GaussSetIntegration ");
  45
  46 // Initialisations
  47    Done = Standard_False;
  48
  49    Standard_Real Xdeb = Lower.Value( Lower.Lower() );
  50    Standard_Real Xfin = Upper.Value( Upper.Lower() );
  51    Standard_Integer Ordre = Order.Value(Order.Lower());
  52    Standard_Real Xm, Xr;
  53    math_Vector GaussP(1, Ordre), GaussW(1, Ordre);
  54
  55 // Recuperation des points de Gauss dans le fichier GaussPoints.
  56    math::GaussPoints  (Ordre, GaussP);
  57    math::GaussWeights (Ordre, GaussW);
  58
  59
  60 // Changement de variable pour la mise a l'echelle [Lower, Upper] :
  61    Xm = 0.5 * (Xdeb + Xfin);
  62    Xr = 0.5 * (Xfin - Xdeb);
  63
  64    Standard_Integer ind = Ordre/2, ind1 = (Ordre+1)/2;
  65    if(ind1 > ind) { // odder case
  66        Tval(1) =  Xm; // +  Xr * GaussP(ind1);
  67        IsOk = F.Value(Tval, Val);
  68        if (!IsOk) return;
  69        Val *= GaussW(ind1);
  70      }
  71    else {
  72      Val.Init(0);
  73    }
  74
  75    for (i=1; i<= ind; i++) {
  76        Tval(1) =  Xm +  Xr * GaussP(i);
  77        IsOk = F.Value(Tval, FVal1);
  78        if (!IsOk) return;
  79        Tval(1) =  Xm -  Xr * GaussP(i);
  80        IsOk = F.Value(Tval, FVal2);
  81        if (!IsOk) return;
  82        FVal1 += FVal2;
  83        FVal1 *=  GaussW(i);
  84        Val += FVal1;
  85      }
  86    Val *= Xr;
  87
  88    Done = Standard_True;
  89  }
  90
  91 void math_GaussSetIntegration::Dump(Standard_OStream& o) const
  92 {
  93   o <<"math_GaussSetIntegration ";
  94    if (Done) {
  95      o << " Status = Done \n";
  96      o << "Integration Value = " << Val<<"\n";
  97    }
  98    else {
  99      o << "Status = not Done \n";
 100    }
 101 }