src/math/math_ComputeGaussPointsAndWeights.cxx

   1 // Created on: 2005-12-19
   2 // Created by: Julia GERASIMOVA
   3 // Copyright (c) 2005-2014 OPEN CASCADE SAS
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   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.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement.
  15
  16
  17 #include <math_Array1OfValueAndWeight.hxx>
  18 #include <math_ComputeGaussPointsAndWeights.hxx>
  19 #include <math_EigenValuesSearcher.hxx>
  20 #include <Standard_ErrorHandler.hxx>
  21
  22 #include <algorithm>
  23 math_ComputeGaussPointsAndWeights::math_ComputeGaussPointsAndWeights(const Standard_Integer Number)
  24 {
  25   myIsDone = Standard_False;
  26
  27   try {
  28     myPoints  = new TColStd_HArray1OfReal(1, Number);
  29     myWeights = new TColStd_HArray1OfReal(1, Number);
  30
  31     Standard_Integer i;
  32
  33     TColStd_Array1OfReal aDiag(1, Number);
  34     TColStd_Array1OfReal aSubDiag(1, Number);
  35
  36     //Initialization of a real symmetric tridiagonal matrix for
  37     //computation of Gauss quadrature.
  38
  39     for (i = 1; i <= Number; i++) {
  40       aDiag(i) = 0.;
  41
  42       if (i == 1)
  43         aSubDiag(i) = 0.;
  44       else {
  45         Standard_Integer sqrIm1 = (i-1)*(i-1);
  46         aSubDiag(i) = sqrIm1/(4.*sqrIm1 - 1);
  47         aSubDiag(i) = Sqrt(aSubDiag(i));
  48       }
  49     }
  50
  51     // Compute eigen values.
  52     math_EigenValuesSearcher EVsearch(aDiag, aSubDiag);
  53
  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 );
  61         VWarray(i) = EVW;
  62       }
  63
  64       std::sort (VWarray.begin(), VWarray.end());
  65
  66       for (i = 1; i <= Number; i++) {
  67         myPoints->ChangeValue(i)  = VWarray(i).Value();
  68         myWeights->ChangeValue(i) = VWarray(i).Weight();
  69       }
  70       myIsDone = Standard_True;
  71     }
  72   } catch (Standard_Failure) {
  73   }
  74 }
  75
  76 Standard_Boolean math_ComputeGaussPointsAndWeights::IsDone() const
  77 {
  78   return myIsDone;
  79 }
  80
  81 math_Vector math_ComputeGaussPointsAndWeights::Points() const
  82 {
  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);
  87
  88   return thePoints;
  89 }
  90
  91 math_Vector math_ComputeGaussPointsAndWeights::Weights() const
  92 {
  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);
  97
  98   return theWeights;
  99 }