1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
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.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
16 #define No_Standard_RangeError
17 #define No_Standard_OutOfRange
18 #define No_Standard_DimensionError
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>
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()) {
35 Standard_Integer NbEqua = F.NbEquations() , NbVar = F.NbVariables();
37 Standard_Boolean IsOk;
38 math_Vector FVal1(1, NbEqua), FVal2(1, NbEqua), Tval(1, NbVar);
42 Standard_NotImplemented_Raise_if(
43 NbVar != 1 || Order.Value(Order.Lower()) > math::GaussPointsMax(),
44 "GaussSetIntegration ");
47 Done = Standard_False;
49 Standard_Real Xdeb = Lower.Value( Lower.Lower() );
50 Standard_Real Xfin = Upper.Value( Upper.Lower() );
51 Standard_Integer Ordre = Order.Value(Order.Lower());
53 math_Vector GaussP(1, Ordre), GaussW(1, Ordre);
55 // Recuperation des points de Gauss dans le fichier GaussPoints.
56 math::GaussPoints (Ordre, GaussP);
57 math::GaussWeights (Ordre, GaussW);
60 // Changement de variable pour la mise a l'echelle [Lower, Upper] :
61 Xm = 0.5 * (Xdeb + Xfin);
62 Xr = 0.5 * (Xfin - Xdeb);
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);
75 for (i=1; i<= ind; i++) {
76 Tval(1) = Xm + Xr * GaussP(i);
77 IsOk = F.Value(Tval, FVal1);
79 Tval(1) = Xm - Xr * GaussP(i);
80 IsOk = F.Value(Tval, FVal2);
91 void math_GaussSetIntegration::Dump(Standard_OStream& o) const
93 o <<"math_GaussSetIntegration ";
95 o << " Status = Done \n";
96 o << "Integration Value = " << Val<<"\n";
99 o << "Status = not Done \n";