src/BRepGProp/BRepGProp_TFunction.cxx

   1 // Created on: 2005-12-09
   2 // Created by: Sergey KHROMOV
   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 <BRepGProp_Face.hxx>
  18 #include <BRepGProp_TFunction.hxx>
  19 #include <gp_Pnt.hxx>
  20 #include <math_KronrodSingleIntegration.hxx>
  21 #include <Precision.hxx>
  22 #include <TColStd_HArray1OfReal.hxx>
  23
  24 //=======================================================================
  25 //function : Constructor.
  26 //purpose  :
  27 //=======================================================================
  28 BRepGProp_TFunction::BRepGProp_TFunction(const BRepGProp_Face   &theSurface,
  29                                          const gp_Pnt           &theVertex,
  30                                          const Standard_Boolean  IsByPoint,
  31                                          const Standard_Address  theCoeffs,
  32                                          const Standard_Real     theUMin,
  33                                          const Standard_Real     theTolerance):
  34   mySurface(theSurface),
  35   myUFunction(mySurface, theVertex, IsByPoint, theCoeffs),
  36   myUMin(theUMin),
  37   myTolerance(theTolerance),
  38   myTolReached(0.),
  39   myErrReached(0.),
  40   myAbsError(0.),
  41   myValueType(GProp_Unknown),
  42   myIsByPoint(IsByPoint),
  43   myNbPntOuter(3)
  44 {
  45 }
  46
  47 //=======================================================================
  48 //function : Init
  49 //purpose  :
  50 //=======================================================================
  51
  52 void BRepGProp_TFunction::Init()
  53 {
  54   myTolReached = 0.;
  55   myErrReached = 0.;
  56   myAbsError = 0.;
  57 }
  58
  59 //=======================================================================
  60 //function : Value
  61 //purpose  :
  62 //=======================================================================
  63
  64 Standard_Boolean BRepGProp_TFunction::Value(const Standard_Real X,
  65                                             Standard_Real &F)
  66 {
  67   const Standard_Real tolU = 1.e-9;
  68
  69   gp_Pnt2d                      aP2d;
  70   gp_Vec2d                      aV2d;
  71   Standard_Real                 aUMax;
  72   Handle(TColStd_HArray1OfReal) anUKnots;
  73
  74   mySurface.D12d(X, aP2d, aV2d);
  75   aUMax = aP2d.X();
  76
  77   if(aUMax - myUMin < tolU)
  78   {
  79     F = 0.;
  80     return Standard_True;
  81   }
  82
  83   mySurface.GetUKnots(myUMin, aUMax, anUKnots);
  84   myUFunction.SetVParam(aP2d.Y());
  85
  86   // Compute the integral from myUMin to aUMax of myUFunction.
  87   Standard_Integer i;
  88   Standard_Real    aCoeff        = aV2d.Y();
  89   //Standard_Integer aNbUIntervals = anUKnots->Length() - 1;
  90   //Standard_Real    aTol          = myTolerance/aNbUIntervals;
  91   Standard_Real    aTol          = myTolerance;
  92
  93   //aTol /= myNbPntOuter;
  94
  95   if (myValueType == GProp_Mass) {
  96     if (myIsByPoint)
  97       aCoeff /= 3.;
  98   } else if (myValueType == GProp_CenterMassX ||
  99     myValueType == GProp_CenterMassY ||
 100     myValueType == GProp_CenterMassZ) {
 101       if (myIsByPoint)
 102         aCoeff *= 0.25;
 103   } else if (myValueType == GProp_InertiaXX ||
 104     myValueType == GProp_InertiaYY ||
 105     myValueType == GProp_InertiaZZ ||
 106     myValueType == GProp_InertiaXY ||
 107     myValueType == GProp_InertiaXZ ||
 108     myValueType == GProp_InertiaYZ) {
 109       if (myIsByPoint)
 110         aCoeff *= 0.2;
 111   } else
 112     return Standard_False;
 113
 114   Standard_Real aAbsCoeff = Abs(aCoeff);
 115
 116   if (aAbsCoeff <= Precision::Angular()) {
 117     // No need to compute the integral. The value will be equal to 0.
 118     F = 0.;
 119     return Standard_True;
 120   }
 121   //else if (aAbsCoeff > 10.*aTol)
 122   //  aTol /= aAbsCoeff;
 123   //else
 124   //  aTol = 0.1;
 125
 126   Standard_Integer              iU           = anUKnots->Upper();
 127   Standard_Integer              aNbPntsStart;
 128   Standard_Integer              aNbMaxIter   = 1000;
 129   math_KronrodSingleIntegration anIntegral;
 130   Standard_Real                 aLocalErr    = 0.;
 131
 132   i            = anUKnots->Lower();
 133   F            = 0.;
 134
 135   // Epmirical criterion
 136   aNbPntsStart = Min(15, mySurface.UIntegrationOrder()/(anUKnots->Length() - 1)+1);
 137   aNbPntsStart = Max(5, aNbPntsStart);
 138
 139   while (i < iU) {
 140     Standard_Real aU1 = anUKnots->Value(i++);
 141     Standard_Real aU2 = anUKnots->Value(i);
 142
 143     if(aU2 - aU1 < tolU) continue;
 144
 145     anIntegral.Perform(myUFunction, aU1, aU2, aNbPntsStart, aTol, aNbMaxIter);
 146
 147     if (!anIntegral.IsDone())
 148       return Standard_False;
 149
 150     F         += anIntegral.Value();
 151     aLocalErr += anIntegral.AbsolutError();
 152     //std::cout << " TFunction : " << anIntegral.NbIterReached() << " " << anIntegral.AbsolutError() << std::endl;
 153   }
 154
 155   F            *= aCoeff;
 156   aLocalErr *= aAbsCoeff;
 157   myAbsError = Max(myAbsError, aLocalErr);
 158
 159   myTolReached += aLocalErr;
 160
 161   if(Abs(F) > Epsilon(1.)) aLocalErr /= Abs(F);
 162
 163   myErrReached = Max(myErrReached, aLocalErr);
 164
 165
 166   return Standard_True;
 167 }
 168
 169 //=======================================================================
 170 //function : GetStateNumber
 171 //purpose  :
 172 //=======================================================================
 173
 174 Standard_Integer BRepGProp_TFunction::GetStateNumber()
 175 {
 176   //myErrReached  = myTolReached;
 177   //myTolReached  = 0.;
 178   //myNbPntOuter += RealToInt(0.5*myNbPntOuter);
 179
 180   //if (myNbPntOuter%2 == 0)
 181   //myNbPntOuter++;
 182
 183   return 0;
 184 }