src/BRepGProp/BRepGProp_VinertGK.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_Domain.hxx>
  18 #include <BRepGProp_Face.hxx>
  19 #include <BRepGProp_TFunction.hxx>
  20 #include <BRepGProp_VinertGK.hxx>
  21 #include <gp_Pln.hxx>
  22 #include <gp_Pnt.hxx>
  23 #include <math_KronrodSingleIntegration.hxx>
  24 #include <TColStd_Array1OfBoolean.hxx>
  25 #include <TColStd_HArray1OfReal.hxx>
  26
  27 //==========================================================================
  28 //function : Constructor
  29 //==========================================================================
  30 BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face        &theSurface,
  31                                        const gp_Pnt          &theLocation,
  32                                        const Standard_Real    theTolerance,
  33                                        const Standard_Boolean theCGFlag,
  34                                        const Standard_Boolean theIFlag):
  35   myErrorReached(0.)
  36 {
  37   SetLocation(theLocation);
  38   Perform(theSurface, theTolerance, theCGFlag, theIFlag);
  39 }
  40
  41 //==========================================================================
  42 //function : Constructor
  43 //
  44 //==========================================================================
  45
  46 BRepGProp_VinertGK::BRepGProp_VinertGK(      BRepGProp_Face          &theSurface,
  47                                        const gp_Pnt        &thePoint,
  48                                        const gp_Pnt        &theLocation,
  49                                        const Standard_Real  theTolerance,
  50                                        const Standard_Boolean theCGFlag,
  51                                        const Standard_Boolean theIFlag):
  52   myErrorReached(0.)
  53 {
  54   SetLocation(theLocation);
  55   Perform(theSurface, thePoint, theTolerance, theCGFlag, theIFlag);
  56 }
  57
  58 //==========================================================================
  59 //function : Constructor
  60 //
  61 //==========================================================================
  62
  63 BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face        &theSurface,
  64                                        BRepGProp_Domain      &theDomain,
  65                                        const gp_Pnt          &theLocation,
  66                                        const Standard_Real    theTolerance,
  67                                        const Standard_Boolean theCGFlag,
  68                                        const Standard_Boolean theIFlag):
  69   myErrorReached(0.)
  70 {
  71   SetLocation(theLocation);
  72   Perform(theSurface, theDomain, theTolerance, theCGFlag, theIFlag);
  73 }
  74
  75 //==========================================================================
  76 //function : Constructor
  77 //
  78 //==========================================================================
  79
  80 BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face        &theSurface,
  81                                        BRepGProp_Domain      &theDomain,
  82                                        const gp_Pnt          &thePoint,
  83                                        const gp_Pnt          &theLocation,
  84                                        const Standard_Real    theTolerance,
  85                                        const Standard_Boolean theCGFlag,
  86                                        const Standard_Boolean theIFlag):
  87   myErrorReached(0.)
  88 {
  89   SetLocation(theLocation);
  90   Perform(theSurface, theDomain, thePoint, theTolerance, theCGFlag, theIFlag);
  91 }
  92
  93 //==========================================================================
  94 //function : Constructor
  95 //
  96 //==========================================================================
  97
  98 BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face        &theSurface,
  99                                        const gp_Pln          &thePlane,
 100                                        const gp_Pnt          &theLocation,
 101                                        const Standard_Real    theTolerance,
 102                                        const Standard_Boolean theCGFlag,
 103                                        const Standard_Boolean theIFlag):
 104   myErrorReached(0.)
 105 {
 106   SetLocation(theLocation);
 107   Perform(theSurface, thePlane, theTolerance, theCGFlag, theIFlag);
 108 }
 109
 110 //==========================================================================
 111 //function : Constructor
 112 //
 113 //==========================================================================
 114
 115 BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face        &theSurface,
 116                                        BRepGProp_Domain      &theDomain,
 117                                        const gp_Pln          &thePlane,
 118                                        const gp_Pnt          &theLocation,
 119                                        const Standard_Real    theTolerance,
 120                                        const Standard_Boolean theCGFlag,
 121                                        const Standard_Boolean theIFlag):
 122   myErrorReached(0.)
 123 {
 124   SetLocation(theLocation);
 125   Perform(theSurface, theDomain, thePlane, theTolerance, theCGFlag, theIFlag);
 126 }
 127
 128 //==========================================================================
 129 //function : Perform
 130 //           Compute the properties.
 131 //==========================================================================
 132
 133 Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face        &theSurface,
 134                                           const Standard_Real    theTolerance,
 135                                           const Standard_Boolean theCGFlag,
 136                                           const Standard_Boolean theIFlag)
 137
 138 {
 139   Standard_Real aShift[] = { 0., 0., 0. };
 140
 141   return PrivatePerform(theSurface, NULL, Standard_True, &aShift, theTolerance,
 142     theCGFlag, theIFlag);
 143 }
 144
 145 //==========================================================================
 146 //function : Perform
 147 //           Compute the properties.
 148 //==========================================================================
 149
 150 Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face        &theSurface,
 151                                           const gp_Pnt          &thePoint,
 152                                           const Standard_Real    theTolerance,
 153                                           const Standard_Boolean theCGFlag,
 154                                           const Standard_Boolean theIFlag)
 155
 156 {
 157   gp_XYZ        aXYZ(thePoint.XYZ().Subtracted(loc.XYZ()));
 158   Standard_Real aShift[3];
 159
 160   aXYZ.Coord(aShift[0], aShift[1], aShift[2]);
 161
 162   return PrivatePerform(theSurface, NULL, Standard_True, &aShift, theTolerance,
 163     theCGFlag, theIFlag);
 164 }
 165
 166 //==========================================================================
 167 //function : Perform
 168 //           Compute the properties.
 169 //==========================================================================
 170
 171 Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face        &theSurface,
 172                                           BRepGProp_Domain      &theDomain,
 173                                           const Standard_Real    theTolerance,
 174                                           const Standard_Boolean theCGFlag,
 175                                           const Standard_Boolean theIFlag)
 176
 177 {
 178   Standard_Real aShift[] = { 0., 0., 0. };
 179
 180   return PrivatePerform(theSurface, &theDomain,
 181     Standard_True, &aShift, theTolerance,
 182     theCGFlag, theIFlag);
 183 }
 184
 185 //==========================================================================
 186 //function : Perform
 187 //           Compute the properties.
 188 //==========================================================================
 189
 190 Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face        &theSurface,
 191                                           BRepGProp_Domain      &theDomain,
 192                                           const gp_Pnt          &thePoint,
 193                                           const Standard_Real    theTolerance,
 194                                           const Standard_Boolean theCGFlag,
 195                                           const Standard_Boolean theIFlag)
 196
 197 {
 198   gp_XYZ        aXYZ(thePoint.XYZ().Subtracted(loc.XYZ()));
 199   Standard_Real aShift[3];
 200
 201   aXYZ.Coord(aShift[0], aShift[1], aShift[2]);
 202
 203   return PrivatePerform(theSurface, &theDomain,
 204     Standard_True, &aShift, theTolerance,
 205     theCGFlag, theIFlag);
 206 }
 207
 208 //==========================================================================
 209 //function : Perform
 210 //           Compute the properties.
 211 //==========================================================================
 212
 213 Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face        &theSurface,
 214                                           const gp_Pln          &thePlane,
 215                                           const Standard_Real    theTolerance,
 216                                           const Standard_Boolean theCGFlag,
 217                                           const Standard_Boolean theIFlag)
 218
 219 {
 220   Standard_Real aCoeff[4];
 221   Standard_Real aXLoc;
 222   Standard_Real aYLoc;
 223   Standard_Real aZLoc;
 224
 225   loc.Coord(aXLoc, aYLoc, aZLoc);
 226   thePlane.Coefficients (aCoeff[0], aCoeff[1], aCoeff[2], aCoeff[3]);
 227   aCoeff[3] = aCoeff[3] - aCoeff[0]*aXLoc - aCoeff[1]*aYLoc - aCoeff[2]*aZLoc;
 228
 229   return PrivatePerform(theSurface, NULL,
 230     Standard_False, &aCoeff, theTolerance,
 231     theCGFlag, theIFlag);
 232 }
 233
 234 //==========================================================================
 235 //function : Perform
 236 //           Compute the properties.
 237 //==========================================================================
 238
 239 Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face        &theSurface,
 240                                           BRepGProp_Domain      &theDomain,
 241                                           const gp_Pln          &thePlane,
 242                                           const Standard_Real    theTolerance,
 243                                           const Standard_Boolean theCGFlag,
 244                                           const Standard_Boolean theIFlag)
 245
 246 {
 247   Standard_Real aCoeff[4];
 248   Standard_Real aXLoc;
 249   Standard_Real aYLoc;
 250   Standard_Real aZLoc;
 251
 252   loc.Coord(aXLoc, aYLoc, aZLoc);
 253   thePlane.Coefficients (aCoeff[0], aCoeff[1], aCoeff[2], aCoeff[3]);
 254   aCoeff[3] = aCoeff[3] - aCoeff[0]*aXLoc - aCoeff[1]*aYLoc - aCoeff[2]*aZLoc;
 255
 256   return PrivatePerform(theSurface, &theDomain,
 257     Standard_False, &aCoeff, theTolerance,
 258     theCGFlag, theIFlag);
 259 }
 260
 261 //==========================================================================
 262 //function : PrivatePerform
 263 //           Compute the properties.
 264 //==========================================================================
 265
 266 Standard_Real BRepGProp_VinertGK::PrivatePerform
 267 (BRepGProp_Face        &theSurface,
 268  const Standard_Address thePtrDomain,
 269  const Standard_Boolean IsByPoint,
 270  const Standard_Address theCoeffs,
 271  const Standard_Real    theTolerance,
 272  const Standard_Boolean theCGFlag,
 273  const Standard_Boolean theIFlag)
 274
 275 {
 276
 277   const Standard_Real aTTol = 1.e-9;
 278   Standard_Real *aCoeffs = (Standard_Real *)theCoeffs;
 279
 280   // Compute the number of 2d bounding curves of the face.
 281   BRepGProp_Domain           *aPDomain = NULL;
 282   Standard_Integer  aNbCurves = 0;
 283
 284   // If the pointer to the domain is NULL, there is only one curve to treat:
 285   // U isoline with the UMax parameter.
 286   if (thePtrDomain == NULL)
 287     aNbCurves = 1;
 288   else {
 289     aPDomain = (BRepGProp_Domain *)thePtrDomain;
 290
 291     for (aPDomain->Init(); aPDomain->More(); aPDomain->Next())
 292       aNbCurves++;
 293   }
 294
 295   if (aNbCurves == 0) {
 296     myErrorReached = -1.;
 297
 298     return myErrorReached;
 299   }
 300
 301   //Standard_Real    aCrvTol = 0.5*theTolerance/aNbCurves;
 302   Standard_Real    aCrvTol = 0.1*theTolerance;
 303   Standard_Real    aUMin;
 304   Standard_Real    aUMax;
 305   Standard_Real    aTMin;
 306   Standard_Real    aTMax;
 307   Standard_Integer aNbPnts;
 308   Standard_Integer aNbMaxIter = 1000;
 309   Standard_Integer aNbVal = 10;
 310   Standard_Integer k;
 311   math_Vector      aLocalValue(1, aNbVal);
 312   math_Vector      aLocalTolReached(1, aNbVal);
 313   math_Vector      aValue(1, aNbVal);
 314   math_Vector      aTolReached(1, aNbVal);
 315   TColStd_Array1OfBoolean CFlags(1, aNbVal);
 316   CFlags.Init(Standard_False);
 317   Standard_Boolean isMore;
 318
 319   //aNbVal = 1;
 320   aValue.Init(0.);
 321   aTolReached.Init(0.);
 322
 323   CFlags.Init(Standard_False);
 324   CFlags(1) = Standard_True;
 325
 326   if(theCGFlag || theIFlag) {
 327     Standard_Integer i;
 328     for(i = 2; i <= 4; ++i) {CFlags(i) = Standard_True;}
 329   }
 330
 331   if(theIFlag) {
 332     Standard_Integer i;
 333     for(i = 5; i <= 10; ++i) {CFlags(i) = Standard_True;}
 334   }
 335
 336   theSurface.Bounds(aUMin, aUMax, aTMin, aTMax);
 337
 338   if (thePtrDomain == NULL)
 339     isMore = Standard_True;
 340   else {
 341     aPDomain->Init();
 342     isMore = aPDomain->More();
 343   }
 344
 345   while(isMore) {
 346     // If the pointer to the domain is NULL, there is only one curve to treat:
 347     // U isoline with the UMax parameter.
 348
 349     if (thePtrDomain == NULL)
 350       theSurface.Load(Standard_False, GeomAbs_IsoU);
 351     else
 352       theSurface.Load(aPDomain->Value());
 353
 354     aTMin = theSurface.FirstParameter();
 355     aTMax = theSurface.LastParameter();
 356
 357
 358     // Get the spans on the curve.
 359     Handle(TColStd_HArray1OfReal) aTKnots;
 360     BRepGProp_TFunction               aTFunc(theSurface, loc, IsByPoint, theCoeffs,
 361       aUMin, aCrvTol);
 362
 363     theSurface.GetTKnots(aTMin, aTMax, aTKnots);
 364
 365     Standard_Integer              iU             = aTKnots->Upper();
 366     Standard_Integer              aNbTIntervals  = aTKnots->Length() - 1;
 367     //Standard_Real                 aTolSpan       = aCrvTol/aNbTIntervals;
 368     Standard_Real                 aTolSpan       = 0.9*theTolerance; //Relative error
 369     math_KronrodSingleIntegration anIntegral;
 370     GProp_ValueType               aValueType;
 371
 372
 373     // Empirical criterion.
 374     aNbPnts = Min(15, theSurface.IntegrationOrder()/aNbTIntervals + 1);
 375     aNbPnts = Max(5, aNbPnts);
 376     //     aNbPnts = theSurface.IntegrationOrder();
 377
 378     aLocalValue.Init(0.);
 379     aLocalTolReached.Init(0.);
 380
 381     for (k = 1; k <= aNbVal; k++) {
 382
 383       if(!CFlags(k)) continue;
 384
 385       Standard_Integer i = aTKnots->Lower();
 386
 387       switch (k) {
 388       case 1:   aValueType = GProp_Mass;        break;
 389       case 2:   aValueType = GProp_CenterMassX; break;
 390       case 3:   aValueType = GProp_CenterMassY; break;
 391       case 4:   aValueType = GProp_CenterMassZ; break;
 392       case 5:   aValueType = GProp_InertiaXX;   break;
 393       case 6:   aValueType = GProp_InertiaYY;   break;
 394       case 7:   aValueType = GProp_InertiaZZ;   break;
 395       case 8:   aValueType = GProp_InertiaXY;   break;
 396       case 9:   aValueType = GProp_InertiaXZ;   break;
 397       case 10:  aValueType = GProp_InertiaYZ;   break;
 398
 399       default: myErrorReached = -1.; return myErrorReached;
 400       }
 401       aTFunc.SetValueType(aValueType);
 402
 403       Standard_Real err1 = 0.;
 404       while (i < iU) {
 405
 406         //cout << "-------------- Span " << i << " nbp: " << aNbPnts << endl;
 407         Standard_Real aT1 = aTKnots->Value(i++);
 408         Standard_Real aT2 = aTKnots->Value(i);
 409
 410         if(aT2 - aT1 < aTTol) continue;
 411
 412         aTFunc.SetNbKronrodPoints(aNbPnts);
 413         aTFunc.Init();
 414         aTFunc.SetTolerance(aCrvTol/(aT2-aT1));
 415         anIntegral.Perform(aTFunc, aT1, aT2, aNbPnts, aTolSpan, aNbMaxIter);
 416
 417         if (!anIntegral.IsDone()) {
 418           myErrorReached = -1.;
 419
 420           return myErrorReached;
 421         }
 422
 423         aLocalValue(k)      += anIntegral.Value();
 424         err1 = aTFunc.AbsolutError()*(aT2 - aT1);
 425         //cout << "Errors: " << anIntegral.NbIterReached() << " " << anIntegral.AbsolutError() << " " << err1 << endl;
 426         aLocalTolReached(k) += anIntegral.AbsolutError() + err1;
 427         //cout << "--- Errors: " << anIntegral.NbIterReached() << " " << anIntegral.AbsolutError() << " " << err1 << endl;
 428       }
 429
 430       aValue(k)      += aLocalValue(k);
 431       aTolReached(k) += aLocalTolReached(k);
 432     }
 433
 434     // If the pointer to the domain is NULL, there is only one curve to treat:
 435     // U isoline with the UMax parameter.
 436     if (thePtrDomain == NULL)
 437       isMore = Standard_False;
 438     else {
 439       aPDomain->Next();
 440       isMore = aPDomain->More();
 441     }
 442   }
 443
 444   // Get volume value.
 445   dim            = aValue(1);
 446   myErrorReached = aTolReached(1);
 447   myAbsolutError = myErrorReached;
 448   Standard_Real anAbsDim = Abs(dim);
 449   Standard_Real aVolTol  = Epsilon(myAbsolutError);
 450   if(anAbsDim >= aVolTol) myErrorReached /= anAbsDim;
 451
 452   if(theCGFlag || theIFlag) {
 453     // Compute values of center of mass.
 454     if(anAbsDim >= aVolTol) {
 455       if (IsByPoint) {
 456         aValue(2) = aCoeffs[0] + aValue(2)/dim;
 457         aValue(3) = aCoeffs[1] + aValue(3)/dim;
 458         aValue(4) = aCoeffs[2] + aValue(4)/dim;
 459       } else {
 460         aValue(2) /= dim;
 461         aValue(3) /= dim;
 462         aValue(4) /= dim;
 463       }
 464     } else {
 465       aValue(2) = 0.;
 466       aValue(3) = 0.;
 467       aValue(4) = 0.;
 468       dim = 0.;
 469     }
 470     g.SetCoord(aValue(2), aValue(3), aValue(4));
 471   }
 472
 473   if(theIFlag) {
 474     // Fill the matrix of inertia.
 475     inertia.SetCols (gp_XYZ (aValue(5), aValue(8),  aValue(9)),
 476       gp_XYZ (aValue(8), aValue(6),  aValue(10)),
 477       gp_XYZ (aValue(9), aValue(10), aValue(7)));
 478   }
 479   //return myErrorReached;
 480   return myAbsolutError;
 481 }
 482