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