1 // Created on: 2005-12-19
2 // Created by: Julia GERASIMOVA
3 // Copyright (c) 2005-2014 OPEN CASCADE SAS
5 // This file is part of Open CASCADE Technology software library.
7 // This library is free software; you can redistribute it and/or modify it under
8 // the terms of the GNU Lesser General Public License version 2.1 as published
9 // by the Free Software Foundation, with special exception defined in the file
10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11 // distribution for complete text of the license and disclaimer of any warranty.
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
17 #include <math_Array1OfValueAndWeight.hxx>
18 #include <math_ComputeGaussPointsAndWeights.hxx>
19 #include <math_EigenValuesSearcher.hxx>
20 #include <Standard_ErrorHandler.hxx>
23 math_ComputeGaussPointsAndWeights::math_ComputeGaussPointsAndWeights(const Standard_Integer Number)
25 myIsDone = Standard_False;
28 myPoints = new TColStd_HArray1OfReal(1, Number);
29 myWeights = new TColStd_HArray1OfReal(1, Number);
33 TColStd_Array1OfReal aDiag(1, Number);
34 TColStd_Array1OfReal aSubDiag(1, Number);
36 //Initialization of a real symmetric tridiagonal matrix for
37 //computation of Gauss quadrature.
39 for (i = 1; i <= Number; i++) {
45 Standard_Integer sqrIm1 = (i-1)*(i-1);
46 aSubDiag(i) = sqrIm1/(4.*sqrIm1 - 1);
47 aSubDiag(i) = Sqrt(aSubDiag(i));
51 // Compute eigen values.
52 math_EigenValuesSearcher EVsearch(aDiag, aSubDiag);
54 if (EVsearch.IsDone()) {
55 math_Array1OfValueAndWeight VWarray(1, Number);
56 for (i = 1; i <= Number; i++) {
57 math_Vector anEigenVector = EVsearch.EigenVector(i);
58 Standard_Real aWeight = anEigenVector(1);
59 aWeight = 2. * aWeight * aWeight;
60 math_ValueAndWeight EVW( EVsearch.EigenValue(i), aWeight );
64 std::sort (VWarray.begin(), VWarray.end());
66 for (i = 1; i <= Number; i++) {
67 myPoints->ChangeValue(i) = VWarray(i).Value();
68 myWeights->ChangeValue(i) = VWarray(i).Weight();
70 myIsDone = Standard_True;
72 } catch (Standard_Failure const&) {
76 Standard_Boolean math_ComputeGaussPointsAndWeights::IsDone() const
81 math_Vector math_ComputeGaussPointsAndWeights::Points() const
83 Standard_Integer Number = myPoints->Length();
84 math_Vector thePoints(1, Number);
85 for (Standard_Integer i = 1; i <= Number; i++)
86 thePoints(i) = myPoints->Value(i);
91 math_Vector math_ComputeGaussPointsAndWeights::Weights() const
93 Standard_Integer Number = myWeights->Length();
94 math_Vector theWeights(1, Number);
95 for (Standard_Integer i = 1; i <= Number; i++)
96 theWeights(i) = myWeights->Value(i);