src/BRepGProp/BRepGProp_Gauss.cxx

   1 // Copyright (c) 2008-2015 OPEN CASCADE SAS
   2 //
   3 // This file is part of Open CASCADE Technology software library.
   4 //
   5 // This library is free software; you can redistribute it and/or modify it under
   6 // the terms of the GNU Lesser General Public License version 2.1 as published
   7 // by the Free Software Foundation, with special exception defined in the file
   8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
   9 // distribution for complete text of the license and disclaimer of any warranty.
  10 //
  11 // Alternatively, this file may be used under the terms of Open CASCADE
  12 // commercial license or contractual agreement.
  13
  14 #include <math.hxx>
  15 #include <Precision.hxx>
  16 #include <TColStd_Array1OfReal.hxx>
  17 #include <Standard_Assert.hxx>
  18 #include <BRepGProp_Face.hxx>
  19 #include <BRepGProp_Domain.hxx>
  20 #include <BRepGProp_Gauss.hxx>
  21
  22 // If the following is defined the error of algorithm is calculated by static moments
  23 #define IS_MIN_DIM
  24
  25 namespace
  26 {
  27   // Minimal value of interval's range for computation | minimal value of "dim" | ...
  28   static const Standard_Real EPS_PARAM          = 1.e-12;
  29   static const Standard_Real EPS_DIM            = 1.e-30;
  30   static const Standard_Real ERROR_ALGEBR_RATIO = 2.0 / 3.0;
  31
  32   // Maximum of GaussPoints on a subinterval and maximum of subintervals
  33   static const Standard_Integer GPM        = math::GaussPointsMax();
  34   static const Standard_Integer SUBS_POWER = 32;
  35   static const Standard_Integer SM         = SUBS_POWER * GPM + 1;
  36
  37   // Auxiliary inner functions to perform arithmetic operations.
  38   static Standard_Real Add(const Standard_Real theA, const Standard_Real theB)
  39   {
  40     return theA + theB;
  41   }
  42
  43   static Standard_Real AddInf(const Standard_Real theA, const Standard_Real theB)
  44   {
  45     if (Precision::IsPositiveInfinite(theA))
  46     {
  47       if (Precision::IsNegativeInfinite(theB))
  48         return 0.0;
  49       else
  50         return Precision::Infinite();
  51     }
  52
  53     if (Precision::IsPositiveInfinite(theB))
  54     {
  55       if (Precision::IsNegativeInfinite(theA))
  56         return 0.0;
  57       else
  58         return Precision::Infinite();
  59     }
  60
  61     if (Precision::IsNegativeInfinite(theA))
  62     {
  63       if (Precision::IsPositiveInfinite(theB))
  64         return 0.0;
  65       else
  66         return -Precision::Infinite();
  67     }
  68
  69     if (Precision::IsNegativeInfinite(theB))
  70     {
  71       if (Precision::IsPositiveInfinite(theA))
  72         return 0.0;
  73       else
  74         return -Precision::Infinite();
  75     }
  76
  77     return theA + theB;
  78   }
  79
  80   static Standard_Real Mult(const Standard_Real theA, const Standard_Real theB)
  81   {
  82     return theA * theB;
  83   }
  84
  85   static Standard_Real MultInf(const Standard_Real theA, const Standard_Real theB)
  86   {
  87     if ((theA == 0.0) || (theB == 0.0))  //strictly zerro (without any tolerances)
  88       return 0.0;
  89
  90     if (Precision::IsPositiveInfinite(theA))
  91     {
  92       if (theB < 0.0)
  93         return -Precision::Infinite();
  94       else
  95         return Precision::Infinite();
  96     }
  97
  98     if (Precision::IsPositiveInfinite(theB))
  99     {
 100       if (theA < 0.0)
 101         return -Precision::Infinite();
 102       else
 103         return Precision::Infinite();
 104     }
 105
 106     if (Precision::IsNegativeInfinite(theA))
 107     {
 108       if (theB < 0.0)
 109         return +Precision::Infinite();
 110       else
 111         return -Precision::Infinite();
 112     }
 113
 114     if (Precision::IsNegativeInfinite(theB))
 115     {
 116       if (theA < 0.0)
 117         return +Precision::Infinite();
 118       else
 119         return -Precision::Infinite();
 120     }
 121
 122     return theA * theB;
 123   }
 124 }
 125
 126 //=======================================================================
 127 //function : BRepGProp_Gauss::Inert::Inert
 128 //purpose  : Constructor
 129 //=======================================================================
 130 BRepGProp_Gauss::Inertia::Inertia()
 131 : Mass(0.0),
 132   Ix  (0.0),
 133   Iy  (0.0),
 134   Iz  (0.0),
 135   Ixx (0.0),
 136   Iyy (0.0),
 137   Izz (0.0),
 138   Ixy (0.0),
 139   Ixz (0.0),
 140   Iyz (0.0)
 141 {
 142 }
 143
 144 //=======================================================================
 145 //function : Inertia::Reset
 146 //purpose  : Zeroes all values.
 147 //=======================================================================
 148 void BRepGProp_Gauss::Inertia::Reset()
 149 {
 150   memset(reinterpret_cast<void*>(this), 0, sizeof(BRepGProp_Gauss::Inertia));
 151 }
 152
 153 //=======================================================================
 154 //function : BRepGProp_Gauss
 155 //purpose  : Constructor
 156 //=======================================================================
 157 BRepGProp_Gauss::BRepGProp_Gauss(const BRepGProp_GaussType theType)
 158 : myType(theType)
 159 {
 160   add  = (::Add );
 161   mult = (::Mult);
 162 }
 163
 164 //=======================================================================
 165 //function : MaxSubs
 166 //purpose  :
 167 //=======================================================================
 168 Standard_Integer BRepGProp_Gauss::MaxSubs(const Standard_Integer theN,
 169                                           const Standard_Integer theCoeff)
 170 {
 171   return IntegerLast() / theCoeff < theN ?
 172          IntegerLast() : theN * theCoeff + 1;
 173 }
 174
 175 //=======================================================================
 176 //function : Init
 177 //purpose  :
 178 //=======================================================================
 179 void BRepGProp_Gauss::Init(NCollection_Handle<math_Vector>&         theOutVec,
 180                            const Standard_Real    theValue,
 181                            const Standard_Integer theFirst,
 182                            const Standard_Integer theLast)
 183 {
 184   if(theLast - theFirst == 0)
 185   {
 186     theOutVec->Init(theValue);
 187   }
 188   else
 189   {
 190     for (Standard_Integer i = theFirst; i <= theLast; ++i)
 191       theOutVec->Value(i) = theValue;
 192   }
 193 }
 194
 195 //=======================================================================
 196 //function : InitMass
 197 //purpose  :
 198 //=======================================================================
 199 void BRepGProp_Gauss::InitMass(const Standard_Real    theValue,
 200                                const Standard_Integer theFirst,
 201                                const Standard_Integer theLast,
 202                                InertiaArray&          theArray)
 203 {
 204   if (theArray.IsNull())
 205     return;
 206
 207   Standard_Integer aFirst = theFirst;
 208   Standard_Integer aLast  = theLast;
 209
 210   if (theLast - theFirst == 0)
 211   {
 212     aFirst = theArray->Lower();
 213     aLast  = theArray->Upper();
 214   }
 215
 216   for (Standard_Integer i = aFirst; i <= aLast; ++i)
 217       theArray->ChangeValue(i).Mass = theValue;
 218 }
 219
 220 //=======================================================================
 221 //function : FillIntervalBounds
 222 //purpose  :
 223 //=======================================================================
 224 Standard_Integer BRepGProp_Gauss::FillIntervalBounds(
 225   const Standard_Real         theA,
 226   const Standard_Real         theB,
 227   const TColStd_Array1OfReal& theKnots,
 228   const Standard_Integer      theNumSubs,
 229   InertiaArray&               theInerts,
 230   NCollection_Handle<math_Vector>&              theParam1,
 231   NCollection_Handle<math_Vector>&              theParam2,
 232   NCollection_Handle<math_Vector>&              theError,
 233   NCollection_Handle<math_Vector>&              theCommonError)
 234 {
 235   const Standard_Integer aSize =
 236     Max(theKnots.Upper(), MaxSubs(theKnots.Upper() - 1, theNumSubs));
 237
 238   if (aSize - 1 > theParam1->Upper())
 239   {
 240     theInerts = new NCollection_Array1<Inertia>(1, aSize);
 241     theParam1 = new math_Vector(1, aSize);
 242     theParam2 = new math_Vector(1, aSize);
 243     theError  = new math_Vector(1, aSize, 0.0);
 244
 245     if (theCommonError.IsNull() == Standard_False)
 246       theCommonError = new math_Vector(1, aSize, 0.0);
 247   }
 248
 249   Standard_Integer j = 1, k = 1;
 250   theParam1->Value(j++) = theA;
 251
 252   const Standard_Integer aLength = theKnots.Upper();
 253   for (Standard_Integer i = 1; i <= aLength; ++i)
 254   {
 255     const Standard_Real kn = theKnots(i);
 256     if (theA < kn)
 257     {
 258       if (kn < theB)
 259       {
 260         theParam1->Value(j++) = kn;
 261         theParam2->Value(k++) = kn;
 262       }
 263       else
 264         break;
 265     }
 266   }
 267
 268   theParam2->Value(k) = theB;
 269   return k;
 270 }
 271
 272
 273
 274 //=======================================================================
 275 //function : computeVInertiaOfElementaryPart
 276 //purpose  :
 277 //=======================================================================
 278 void BRepGProp_Gauss::computeVInertiaOfElementaryPart(
 279   const gp_Pnt&             thePoint,
 280   const gp_Vec&             theNormal,
 281   const gp_Pnt&             theLocation,
 282   const Standard_Real       theWeight,
 283   const Standard_Real       theCoeff[],
 284   const Standard_Boolean    theIsByPoint,
 285   BRepGProp_Gauss::Inertia& theOutInertia)
 286 {
 287   Standard_Real x = thePoint.X() - theLocation.X();
 288   Standard_Real y = thePoint.Y() - theLocation.Y();
 289   Standard_Real z = thePoint.Z() - theLocation.Z();
 290
 291   const Standard_Real xn = theNormal.X() * theWeight;
 292   const Standard_Real yn = theNormal.Y() * theWeight;
 293   const Standard_Real zn = theNormal.Z() * theWeight;
 294
 295   if (theIsByPoint)
 296   {
 297     /////////////////////                        ///////////////////////
 298     //    OFV code     //                        //    Initial code   //
 299     /////////////////////                        ///////////////////////
 300     // modified by APO
 301
 302     Standard_Real dv = x * xn + y * yn + z * zn; //xyz  = x * y * z;
 303     theOutInertia.Mass += dv / 3.0;              //Ixyi += zn * xyz;
 304     theOutInertia.Ix   += 0.25 * x * dv;         //Iyzi += xn * xyz;
 305     theOutInertia.Iy   += 0.25 * y * dv;         //Ixzi += yn * xyz;
 306     theOutInertia.Iz   += 0.25 * z * dv;         //xi = x * x * x * xn / 3.0;
 307     x -= theCoeff[0];                            //yi = y * y * y * yn / 3.0;
 308     y -= theCoeff[1];                            //zi = z * z * z * zn / 3.0;
 309     z -= theCoeff[2];                            //Ixxi += (yi + zi);
 310     dv *= 0.2;                                   //Iyyi += (xi + zi);
 311     theOutInertia.Ixy -= x * y * dv;             //Izzi += (xi + yi);
 312     theOutInertia.Iyz -= y * z * dv;             //x -= Coeff[0];
 313     theOutInertia.Ixz -= x * z * dv;             //y -= Coeff[1];
 314     x *= x;                                      //z -= Coeff[2];
 315     y *= y;                                      //dv = x * xn + y * yn + z * zn;
 316     z *= z;                                      //dvi +=  dv;
 317     theOutInertia.Ixx += (y + z) * dv;           //Ixi += x * dv;
 318     theOutInertia.Iyy += (x + z) * dv;           //Iyi += y * dv;
 319     theOutInertia.Izz += (x + y) * dv;           //Izi += z * dv;
 320   }
 321   else
 322   { // By plane
 323     const Standard_Real s = xn * theCoeff[0] + yn * theCoeff[1] + zn * theCoeff[2];
 324
 325     Standard_Real d1  = theCoeff[0] * x + theCoeff[1] * y + theCoeff[2] * z - theCoeff[3];
 326     Standard_Real d2  = d1 * d1;
 327     Standard_Real d3  = d1 * d2 / 3.0;
 328     Standard_Real dv  = s  * d1;
 329
 330     theOutInertia.Mass += dv;
 331     theOutInertia.Ix   += (x - (theCoeff[0] * d1 * 0.5)) * dv;
 332     theOutInertia.Iy   += (y - (theCoeff[1] * d1 * 0.5)) * dv;
 333     theOutInertia.Iz   += (z - (theCoeff[2] * d1 * 0.5)) * dv;
 334
 335     const Standard_Real px = x - theCoeff[0] * d1;
 336     const Standard_Real py = y - theCoeff[1] * d1;
 337     const Standard_Real pz = z - theCoeff[2] * d1;
 338
 339     x = px * px * d1 + px * theCoeff[0] * d2 + theCoeff[0] * theCoeff[0] * d3;
 340     y = py * py * d1 + py * theCoeff[1] * d2 + theCoeff[1] * theCoeff[1] * d3;
 341     z = pz * pz * d1 + pz * theCoeff[2] * d2 + theCoeff[2] * theCoeff[2] * d3;
 342
 343     theOutInertia.Ixx += (y + z) * s;
 344     theOutInertia.Iyy += (x + z) * s;
 345     theOutInertia.Izz += (x + y) * s;
 346
 347     d2 *= 0.5;
 348     x = (py * pz * d1) + (py * theCoeff[2] * d2) + (pz * theCoeff[1] * d2) + (theCoeff[1] * theCoeff[2] * d3);
 349     y = (px * pz * d1) + (pz * theCoeff[0] * d2) + (px * theCoeff[2] * d2) + (theCoeff[0] * theCoeff[2] * d3);
 350     z = (px * py * d1) + (px * theCoeff[1] * d2) + (py * theCoeff[0] * d2) + (theCoeff[0] * theCoeff[1] * d3);
 351
 352     theOutInertia.Ixy -= z * s;
 353     theOutInertia.Iyz -= x * s;
 354     theOutInertia.Ixz -= y * s;
 355   }
 356 }
 357
 358 //=======================================================================
 359 //function : computeSInertiaOfElementaryPart
 360 //purpose  :
 361 //=======================================================================
 362 void BRepGProp_Gauss::computeSInertiaOfElementaryPart(
 363   const gp_Pnt&             thePoint,
 364   const gp_Vec&             theNormal,
 365   const gp_Pnt&             theLocation,
 366   const Standard_Real       theWeight,
 367   BRepGProp_Gauss::Inertia& theOutInertia)
 368 {
 369   // ds - Jacobien (x, y, z) -> (u, v) = ||n||
 370   const Standard_Real ds = mult(theNormal.Magnitude(), theWeight);
 371   const Standard_Real x = add(thePoint.X(), -theLocation.X());
 372   const Standard_Real y = add(thePoint.Y(), -theLocation.Y());
 373   const Standard_Real z = add(thePoint.Z(), -theLocation.Z());
 374
 375   theOutInertia.Mass = add(theOutInertia.Mass, ds);
 376
 377   const Standard_Real XdS = mult(x, ds);
 378   const Standard_Real YdS = mult(y, ds);
 379   const Standard_Real ZdS = mult(z, ds);
 380
 381   theOutInertia.Ix  = add(theOutInertia.Ix,  XdS);
 382   theOutInertia.Iy  = add(theOutInertia.Iy,  YdS);
 383   theOutInertia.Iz  = add(theOutInertia.Iz,  ZdS);
 384   theOutInertia.Ixy = add(theOutInertia.Ixy, mult(x, YdS));
 385   theOutInertia.Iyz = add(theOutInertia.Iyz, mult(y, ZdS));
 386   theOutInertia.Ixz = add(theOutInertia.Ixz, mult(x, ZdS));
 387
 388   const Standard_Real XXdS = mult(x, XdS);
 389   const Standard_Real YYdS = mult(y, YdS);
 390   const Standard_Real ZZdS = mult(z, ZdS);
 391
 392   theOutInertia.Ixx = add(theOutInertia.Ixx, add(YYdS, ZZdS));
 393   theOutInertia.Iyy = add(theOutInertia.Iyy, add(XXdS, ZZdS));
 394   theOutInertia.Izz = add(theOutInertia.Izz, add(XXdS, YYdS));
 395 }
 396
 397 //=======================================================================
 398 //function : checkBounds
 399 //purpose  :
 400 //=======================================================================
 401 void BRepGProp_Gauss::checkBounds(const Standard_Real theU1,
 402                                   const Standard_Real theU2,
 403                                   const Standard_Real theV1,
 404                                   const Standard_Real theV2)
 405 {
 406   if (Precision::IsInfinite(theU1) || Precision::IsInfinite(theU2) ||
 407       Precision::IsInfinite(theV1) || Precision::IsInfinite(theV2))
 408   {
 409     add  = (::AddInf);
 410     mult = (::MultInf);
 411   }
 412 }
 413
 414 //=======================================================================
 415 //function : addAndRestoreInertia
 416 //purpose  :
 417 //=======================================================================
 418 void BRepGProp_Gauss::addAndRestoreInertia(
 419   const BRepGProp_Gauss::Inertia& theInInertia,
 420   BRepGProp_Gauss::Inertia&       theOutInertia)
 421 {
 422   theOutInertia.Mass = add(theOutInertia.Mass, theInInertia.Mass);
 423   theOutInertia.Ix   = add(theOutInertia.Ix,   theInInertia.Ix);
 424   theOutInertia.Iy   = add(theOutInertia.Iy,   theInInertia.Iy);
 425   theOutInertia.Iz   = add(theOutInertia.Iz,   theInInertia.Iz);
 426   theOutInertia.Ixx  = add(theOutInertia.Ixx,  theInInertia.Ixx);
 427   theOutInertia.Iyy  = add(theOutInertia.Iyy,  theInInertia.Iyy);
 428   theOutInertia.Izz  = add(theOutInertia.Izz,  theInInertia.Izz);
 429   theOutInertia.Ixy  = add(theOutInertia.Ixy,  theInInertia.Ixy);
 430   theOutInertia.Ixz  = add(theOutInertia.Ixz,  theInInertia.Ixz);
 431   theOutInertia.Iyz  = add(theOutInertia.Iyz,  theInInertia.Iyz);
 432 }
 433
 434 //=======================================================================
 435 //function : multAndRestoreInertia
 436 //purpose  :
 437 //=======================================================================
 438 void BRepGProp_Gauss::multAndRestoreInertia(
 439   const Standard_Real       theValue,
 440   BRepGProp_Gauss::Inertia& theInOutInertia)
 441 {
 442   theInOutInertia.Mass = mult(theInOutInertia.Mass, theValue);
 443   theInOutInertia.Ix   = mult(theInOutInertia.Ix,   theValue);
 444   theInOutInertia.Iy   = mult(theInOutInertia.Iy,   theValue);
 445   theInOutInertia.Iz   = mult(theInOutInertia.Iz,   theValue);
 446   theInOutInertia.Ixx  = mult(theInOutInertia.Ixx,  theValue);
 447   theInOutInertia.Iyy  = mult(theInOutInertia.Iyy,  theValue);
 448   theInOutInertia.Izz  = mult(theInOutInertia.Izz,  theValue);
 449   theInOutInertia.Ixy  = mult(theInOutInertia.Ixy,  theValue);
 450   theInOutInertia.Ixz  = mult(theInOutInertia.Ixz,  theValue);
 451   theInOutInertia.Iyz  = mult(theInOutInertia.Iyz,  theValue);
 452 }
 453
 454 //=======================================================================
 455 //function : convert
 456 //purpose  :
 457 //=======================================================================
 458 void BRepGProp_Gauss::convert(const BRepGProp_Gauss::Inertia& theInertia,
 459                               gp_Pnt&                         theOutGravityCenter,
 460                               gp_Mat&                         theOutMatrixOfInertia,
 461                               Standard_Real&                  theOutMass)
 462 {
 463   if (Abs(theInertia.Mass) >= EPS_DIM)
 464   {
 465     const Standard_Real anInvMass = 1.0 / theInertia.Mass;
 466     theOutGravityCenter.SetX(theInertia.Ix * anInvMass);
 467     theOutGravityCenter.SetY(theInertia.Iy * anInvMass);
 468     theOutGravityCenter.SetZ(theInertia.Iz * anInvMass);
 469
 470     theOutMass = theInertia.Mass;
 471   }
 472   else
 473   {
 474     theOutMass = 0.0;
 475     theOutGravityCenter.SetCoord(0.0, 0.0, 0.0);
 476   }
 477
 478   theOutMatrixOfInertia = gp_Mat(
 479     gp_XYZ ( theInertia.Ixx, -theInertia.Ixy, -theInertia.Ixz),
 480     gp_XYZ (-theInertia.Ixy,  theInertia.Iyy, -theInertia.Iyz),
 481     gp_XYZ (-theInertia.Ixz, -theInertia.Iyz,  theInertia.Izz));
 482 }
 483
 484 //=======================================================================
 485 //function : convert
 486 //purpose  :
 487 //=======================================================================
 488 void BRepGProp_Gauss::convert(const BRepGProp_Gauss::Inertia& theInertia,
 489                               const Standard_Real             theCoeff[],
 490                               const Standard_Boolean          theIsByPoint,
 491                               gp_Pnt&                         theOutGravityCenter,
 492                               gp_Mat&                         theOutMatrixOfInertia,
 493                               Standard_Real&                  theOutMass)
 494 {
 495   convert(theInertia, theOutGravityCenter, theOutMatrixOfInertia, theOutMass);
 496   if (Abs(theInertia.Mass) >= EPS_DIM && theIsByPoint)
 497   {
 498     const Standard_Real anInvMass = 1.0 / theInertia.Mass;
 499     if (theIsByPoint == Standard_True)
 500     {
 501       theOutGravityCenter.SetX(theCoeff[0] + theInertia.Ix * anInvMass);
 502       theOutGravityCenter.SetY(theCoeff[1] + theInertia.Iy * anInvMass);
 503       theOutGravityCenter.SetZ(theCoeff[2] + theInertia.Iz * anInvMass);
 504     }
 505     else
 506     {
 507       theOutGravityCenter.SetX(theInertia.Ix * anInvMass);
 508       theOutGravityCenter.SetY(theInertia.Iy * anInvMass);
 509       theOutGravityCenter.SetZ(theInertia.Iz * anInvMass);
 510     }
 511
 512     theOutMass = theInertia.Mass;
 513   }
 514   else
 515   {
 516     theOutMass = 0.0;
 517     theOutGravityCenter.SetCoord(0.0, 0.0, 0.0);
 518   }
 519
 520   theOutMatrixOfInertia = gp_Mat(
 521     gp_XYZ (theInertia.Ixx, theInertia.Ixy, theInertia.Ixz),
 522     gp_XYZ (theInertia.Ixy, theInertia.Iyy, theInertia.Iyz),
 523     gp_XYZ (theInertia.Ixz, theInertia.Iyz, theInertia.Izz));
 524 }
 525
 526 //=======================================================================
 527 //function : Compute
 528 //purpose  :
 529 //=======================================================================
 530 Standard_Real BRepGProp_Gauss::Compute(
 531   BRepGProp_Face&        theSurface,
 532   BRepGProp_Domain&      theDomain,
 533   const gp_Pnt&          theLocation,
 534   const Standard_Real    theEps,
 535   const Standard_Real    theCoeff[],
 536   const Standard_Boolean theIsByPoint,
 537   Standard_Real&         theOutMass,
 538   gp_Pnt&                theOutGravityCenter,
 539   gp_Mat&                theOutInertia)
 540 {
 541   const Standard_Boolean isErrorCalculation  =
 542     ( 0.0 > theEps || theEps < 0.001 ) ? Standard_True : Standard_False;
 543   const Standard_Boolean isVerifyComputation =
 544     ( 0.0 < theEps && theEps < 0.001 ) ? Standard_True : Standard_False;
 545
 546   Standard_Real anEpsilon= Abs(theEps);
 547
 548   BRepGProp_Gauss::Inertia anInertia;
 549   InertiaArray anInertiaL = new NCollection_Array1<Inertia>(1, SM);
 550   InertiaArray anInertiaU = new NCollection_Array1<Inertia>(1, SM);
 551
 552   // Prepare Gauss points and weights
 553   NCollection_Handle<math_Vector> LGaussP[2];
 554   NCollection_Handle<math_Vector> LGaussW[2];
 555   NCollection_Handle<math_Vector> UGaussP[2];
 556   NCollection_Handle<math_Vector> UGaussW[2];
 557
 558   const Standard_Integer aNbGaussPoint =
 559     RealToInt(Ceiling(ERROR_ALGEBR_RATIO * GPM));
 560
 561   LGaussP[0] = new math_Vector(1, GPM);
 562   LGaussP[1] = new math_Vector(1, aNbGaussPoint);
 563   LGaussW[0] = new math_Vector(1, GPM);
 564   LGaussW[1] = new math_Vector(1, aNbGaussPoint);
 565
 566   UGaussP[0] = new math_Vector(1, GPM);
 567   UGaussP[1] = new math_Vector(1, aNbGaussPoint);
 568   UGaussW[0] = new math_Vector(1, GPM);
 569   UGaussW[1] = new math_Vector(1, aNbGaussPoint);
 570
 571   NCollection_Handle<math_Vector> L1 = new math_Vector(1, SM);
 572   NCollection_Handle<math_Vector> L2 = new math_Vector(1, SM);
 573   NCollection_Handle<math_Vector> U1 = new math_Vector(1, SM);
 574   NCollection_Handle<math_Vector> U2 = new math_Vector(1, SM);
 575
 576   NCollection_Handle<math_Vector> ErrL  = new math_Vector(1, SM, 0.0);
 577   NCollection_Handle<math_Vector> ErrU  = new math_Vector(1, SM, 0.0);
 578   NCollection_Handle<math_Vector> ErrUL = new math_Vector(1, SM, 0.0);
 579
 580   // Face parametrization in U and V direction
 581   Standard_Real BV1, BV2, BU1, BU2;
 582   theSurface.Bounds(BU1, BU2, BV1, BV2);
 583   checkBounds(BU1, BU2, BV1, BV2);
 584
 585   //
 586   const Standard_Integer NumSubs = SUBS_POWER;
 587   const TopoDS_Face& aF = theSurface.GetFace();
 588   TopoDS_Iterator aWIter(aF);
 589   const Standard_Boolean isNaturalRestriction = !aWIter.More(); //theSurface.NaturalRestriction();
 590
 591   Standard_Real CIx, CIy, CIz, CIxy, CIxz, CIyz;
 592   Standard_Real CDim[2], CIxx[2], CIyy[2], CIzz[2];
 593
 594   // Boundary curve parametrization
 595   Standard_Real u1 = BU1, u2, l1, l2, lm, lr, l, v;
 596
 597   // On the boundary curve u-v
 598   gp_Pnt2d Puv;
 599   gp_Vec2d Vuv;
 600   Standard_Real Dul;  // Dul = Du / Dl
 601
 602   Standard_Integer iLS, iLSubEnd, iGL, iGLEnd, NbLGaussP[2], LRange[2], iL, kL, kLEnd, IL, JL;
 603   Standard_Integer i, iUSubEnd, NbUGaussP[2], URange[2], kU, kUEnd, IU, JU;
 604   Standard_Integer UMaxSubs, LMaxSubs;
 605
 606   Standard_Real ErrorU, ErrorL, ErrorLMax = 0.0, Eps = 0.0, EpsL = 0.0, EpsU = 0.0;
 607   iGLEnd = isErrorCalculation ? 2 : 1;
 608
 609   NbUGaussP[0] = theSurface.SIntOrder(anEpsilon);
 610   NbUGaussP[1] = RealToInt( Ceiling(ERROR_ALGEBR_RATIO * NbUGaussP[0]) );
 611
 612   math::GaussPoints (NbUGaussP[0], *UGaussP[0]);
 613   math::GaussWeights(NbUGaussP[0], *UGaussW[0]);
 614   math::GaussPoints (NbUGaussP[1], *UGaussP[1]);
 615   math::GaussWeights(NbUGaussP[1], *UGaussW[1]);
 616
 617   const Standard_Integer aNbUSubs = theSurface.SUIntSubs();
 618   TColStd_Array1OfReal UKnots(1, aNbUSubs + 1);
 619   theSurface.UKnots(UKnots);
 620
 621   while (isNaturalRestriction || theDomain.More())
 622   {
 623     if (isNaturalRestriction)
 624     {
 625       NbLGaussP[0] = Min(2 * NbUGaussP[0], math::GaussPointsMax());
 626     }
 627     else
 628     {
 629       theSurface.Load(theDomain.Value());
 630       NbLGaussP[0] = theSurface.LIntOrder(anEpsilon);
 631     }
 632
 633     NbLGaussP[1] = RealToInt( Ceiling(ERROR_ALGEBR_RATIO * NbLGaussP[0]) );
 634
 635     math::GaussPoints (NbLGaussP[0], *LGaussP[0]);
 636     math::GaussWeights(NbLGaussP[0], *LGaussW[0]);
 637     math::GaussPoints (NbLGaussP[1], *LGaussP[1]);
 638     math::GaussWeights(NbLGaussP[1], *LGaussW[1]);
 639
 640     const Standard_Integer aNbLSubs =
 641       isNaturalRestriction ? theSurface.SVIntSubs(): theSurface.LIntSubs();
 642     TColStd_Array1OfReal LKnots(1, aNbLSubs + 1);
 643
 644     if (isNaturalRestriction)
 645     {
 646       theSurface.VKnots(LKnots);
 647       l1 = BV1;
 648       l2 = BV2;
 649     }
 650     else
 651     {
 652       theSurface.LKnots(LKnots);
 653       l1 = theSurface.FirstParameter();
 654       l2 = theSurface.LastParameter();
 655     }
 656     ErrorL = 0.0;
 657     kLEnd = 1; JL = 0;
 658
 659     if (Abs(l2 - l1) > EPS_PARAM)
 660     {
 661       iLSubEnd = FillIntervalBounds(l1, l2, LKnots, NumSubs, anInertiaL, L1, L2, ErrL, ErrUL);
 662       LMaxSubs = BRepGProp_Gauss::MaxSubs(iLSubEnd);
 663
 664       if (LMaxSubs > SM)
 665         LMaxSubs = SM;
 666
 667       BRepGProp_Gauss::InitMass(0.0, 1, LMaxSubs, anInertiaL);
 668       BRepGProp_Gauss::Init(ErrL,  0.0, 1, LMaxSubs);
 669       BRepGProp_Gauss::Init(ErrUL, 0.0, 1, LMaxSubs);
 670
 671       do // while: L
 672       {
 673         if (++JL > iLSubEnd)
 674         {
 675           LRange[0] = IL = ErrL->Max();
 676           LRange[1] = JL;
 677           L1->Value(JL) = (L1->Value(IL) + L2->Value(IL)) * 0.5;
 678           L2->Value(JL) = L2->Value(IL);
 679           L2->Value(IL) = L1->Value(JL);
 680         }
 681         else
 682           LRange[0] = IL = JL;
 683
 684         if (JL == LMaxSubs || Abs(L2->Value(JL) - L1->Value(JL)) < EPS_PARAM)
 685         {
 686           if (kLEnd == 1)
 687           {
 688             anInertiaL->ChangeValue(JL).Reset();
 689             ErrL->Value(JL) = 0.0;
 690           }
 691           else
 692           {
 693             --JL;
 694             EpsL = ErrorL;
 695             Eps = EpsL / 0.9;
 696             break;
 697           }
 698         }
 699         else
 700           for (kL = 0; kL < kLEnd; kL++)
 701           {
 702             iLS = LRange[kL];
 703             lm = 0.5 * (L2->Value(iLS) + L1->Value(iLS));
 704             lr = 0.5 * (L2->Value(iLS) - L1->Value(iLS));
 705
 706             CIx = CIy = CIz = CIxy = CIxz = CIyz = 0.0;
 707
 708             for (iGL = 0; iGL < iGLEnd; ++iGL)
 709             {
 710               CDim[iGL] = CIxx[iGL] = CIyy[iGL] = CIzz[iGL] = 0.0;
 711
 712               for (iL = 1; iL <= NbLGaussP[iGL]; iL++)
 713               {
 714                 l = lm + lr * LGaussP[iGL]->Value(iL);
 715                 if (isNaturalRestriction)
 716                 {
 717                   v = l;
 718                   u2 = BU2;
 719                   Dul = LGaussW[iGL]->Value(iL);
 720                 }
 721                 else
 722                 {
 723                   theSurface.D12d (l, Puv, Vuv);
 724                   Dul = Vuv.Y() * LGaussW[iGL]->Value(iL);  // Dul = Du / Dl
 725
 726                   if (Abs(Dul) < EPS_PARAM)
 727                     continue;
 728
 729                   v  = Puv.Y();
 730                   u2 = Puv.X();
 731
 732                   // Check on cause out off bounds of value current parameter
 733                   if (v < BV1)
 734                     v = BV1;
 735                   else if (v > BV2)
 736                     v = BV2;
 737
 738                   if (u2 < BU1)
 739                     u2 = BU1;
 740                   else if (u2 > BU2)
 741                     u2 = BU2;
 742                 }
 743
 744                 ErrUL->Value(iLS) = 0.0;
 745                 kUEnd = 1;
 746                 JU = 0;
 747
 748                 if (Abs(u2 - u1) < EPS_PARAM)
 749                   continue;
 750
 751                 NCollection_Handle<math_Vector> aDummy;
 752                 iUSubEnd = FillIntervalBounds(u1, u2, UKnots, NumSubs, anInertiaU, U1, U2, ErrU, aDummy);
 753                 UMaxSubs = BRepGProp_Gauss::MaxSubs(iUSubEnd);
 754
 755                 if (UMaxSubs > SM)
 756                   UMaxSubs = SM;
 757
 758                 BRepGProp_Gauss::InitMass(0.0, 1, UMaxSubs, anInertiaU);
 759                 BRepGProp_Gauss::Init(ErrU, 0.0, 1, UMaxSubs);
 760                 ErrorU = 0.0;
 761
 762                 do
 763                 {//while: U
 764                   if (++JU > iUSubEnd)
 765                   {
 766                     URange[0] = IU = ErrU->Max();
 767                     URange[1] = JU;
 768
 769                     U1->Value(JU) = (U1->Value(IU) + U2->Value(IU)) * 0.5;
 770                     U2->Value(JU) = U2->Value(IU);
 771                     U2->Value(IU) = U1->Value(JU);
 772                   }
 773                   else
 774                     URange[0] = IU = JU;
 775
 776                   if (JU == UMaxSubs || Abs(U2->Value(JU) - U1->Value(JU)) < EPS_PARAM)
 777                     if (kUEnd == 1)
 778                     {
 779                       ErrU->Value(JU) = 0.0;
 780                       anInertiaU->ChangeValue(JU).Reset();
 781                     }
 782                     else
 783                     {
 784                       --JU;
 785                       EpsU = ErrorU;
 786                       Eps  = 10. * EpsU * Abs((u2 - u1) * Dul);
 787                       EpsL = 0.9 * Eps;
 788                       break;
 789                     }
 790                   else
 791                   {
 792                     gp_Pnt aPoint;
 793                     gp_Vec aNormal;
 794
 795                     for (kU = 0; kU < kUEnd; ++kU)
 796                     {
 797                       BRepGProp_Gauss::Inertia aLocal[2];
 798
 799                       Standard_Integer iUS = URange[kU];
 800                       const Standard_Integer aLength = iGLEnd - iGL;
 801
 802                       const Standard_Real um = 0.5 * (U2->Value(iUS) + U1->Value(iUS));
 803                       const Standard_Real ur = 0.5 * (U2->Value(iUS) - U1->Value(iUS));
 804
 805                       for (Standard_Integer iGU = 0; iGU < aLength; ++iGU)
 806                       {
 807                         for (Standard_Integer iU = 1; iU <= NbUGaussP[iGU]; ++iU)
 808                         {
 809                           Standard_Real w = UGaussW[iGU]->Value(iU);
 810                           const Standard_Real u = um + ur * UGaussP[iGU]->Value(iU);
 811
 812                           theSurface.Normal(u, v, aPoint, aNormal);
 813
 814                           if (myType == Vinert)
 815                           {
 816                             computeVInertiaOfElementaryPart(
 817                               aPoint, aNormal, theLocation, w, theCoeff, theIsByPoint, aLocal[iGU]);
 818                           }
 819                           else
 820                           {
 821                             if (iGU > 0)
 822                               aLocal[iGU].Mass += (w * aNormal.Magnitude());
 823                             else
 824                             {
 825                               computeSInertiaOfElementaryPart(
 826                                 aPoint, aNormal, theLocation, w, aLocal[iGU]);
 827                             }
 828                           }
 829                         }
 830                       }
 831
 832                       BRepGProp_Gauss::Inertia& anUI =
 833                         anInertiaU->ChangeValue(iUS);
 834
 835                       anUI.Mass = mult(aLocal[0].Mass, ur);
 836
 837                       if (myType == Vinert)
 838                       {
 839                         anUI.Ixx = mult(aLocal[0].Ixx, ur);
 840                         anUI.Iyy = mult(aLocal[0].Iyy, ur);
 841                         anUI.Izz = mult(aLocal[0].Izz, ur);
 842                       }
 843
 844                       if (iGL > 0)
 845                         continue;
 846
 847                       Standard_Real aDMass = Abs(aLocal[1].Mass - aLocal[0].Mass);
 848
 849                       if (myType == Vinert)
 850                       {
 851                         aLocal[1].Ixx = Abs(aLocal[1].Ixx - aLocal[0].Ixx);
 852                         aLocal[1].Iyy = Abs(aLocal[1].Iyy - aLocal[0].Iyy);
 853                         aLocal[1].Izz = Abs(aLocal[1].Izz - aLocal[0].Izz);
 854
 855                         anUI.Ix = mult(aLocal[0].Ix, ur);
 856                         anUI.Iy = mult(aLocal[0].Iy, ur);
 857                         anUI.Iz = mult(aLocal[0].Iz, ur);
 858
 859                         anUI.Ixy = mult(aLocal[0].Ixy, ur);
 860                         anUI.Ixz = mult(aLocal[0].Ixz, ur);
 861                         anUI.Iyz = mult(aLocal[0].Iyz, ur);
 862
 863                         #ifndef IS_MIN_DIM
 864                         aDMass = aLocal[1].Ixx + aLocal[1].Iyy + aLocal[1].Izz;
 865                         #endif
 866
 867                         ErrU->Value(iUS) = mult(aDMass, ur);
 868                       }
 869                       else
 870                       {
 871                         anUI.Ix  = mult(aLocal[0].Ix,  ur);
 872                         anUI.Iy  = mult(aLocal[0].Iy,  ur);
 873                         anUI.Iz  = mult(aLocal[0].Iz,  ur);
 874                         anUI.Ixx = mult(aLocal[0].Ixx, ur);
 875                         anUI.Iyy = mult(aLocal[0].Iyy, ur);
 876                         anUI.Izz = mult(aLocal[0].Izz, ur);
 877                         anUI.Ixy = mult(aLocal[0].Ixy, ur);
 878                         anUI.Ixz = mult(aLocal[0].Ixz, ur);
 879                         anUI.Iyz = mult(aLocal[0].Iyz, ur);
 880
 881                         ErrU->Value(iUS) = mult(aDMass, ur);
 882                       }
 883                     }
 884                   }
 885
 886                   if (JU == iUSubEnd)
 887                   {
 888                     kUEnd = 2;
 889                     ErrorU = ErrU->Value(ErrU->Max());
 890                   }
 891                 } while ( (ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1 );
 892
 893                 for (i = 1; i <= JU; ++i)
 894                 {
 895                   const BRepGProp_Gauss::Inertia& anIU =
 896                     anInertiaU->Value(i);
 897
 898                   CDim[iGL] = add(CDim[iGL], mult(anIU.Mass, Dul));
 899                   CIxx[iGL] = add(CIxx[iGL], mult(anIU.Ixx,  Dul));
 900                   CIyy[iGL] = add(CIyy[iGL], mult(anIU.Iyy,  Dul));
 901                   CIzz[iGL] = add(CIzz[iGL], mult(anIU.Izz,  Dul));
 902                 }
 903
 904                 if (iGL > 0)
 905                   continue;
 906
 907                 ErrUL->Value(iLS) = ErrorU * Abs((u2 - u1) * Dul);
 908
 909                 for (i = 1; i <= JU; ++i)
 910                 {
 911                   const BRepGProp_Gauss::Inertia& anIU =
 912                     anInertiaU->Value(i);
 913
 914                   CIx = add(CIx, mult(anIU.Ix, Dul));
 915                   CIy = add(CIy, mult(anIU.Iy, Dul));
 916                   CIz = add(CIz, mult(anIU.Iz, Dul));
 917
 918                   CIxy = add(CIxy, mult(anIU.Ixy, Dul));
 919                   CIxz = add(CIxz, mult(anIU.Ixz, Dul));
 920                   CIyz = add(CIyz, mult(anIU.Iyz, Dul));
 921                 }
 922               }//for: iL
 923             }//for: iGL
 924
 925             BRepGProp_Gauss::Inertia& aLI = anInertiaL->ChangeValue(iLS);
 926
 927             aLI.Mass = mult(CDim[0], lr);
 928             aLI.Ixx  = mult(CIxx[0], lr);
 929             aLI.Iyy  = mult(CIyy[0], lr);
 930             aLI.Izz  = mult(CIzz[0], lr);
 931
 932             if (iGLEnd == 2)
 933             {
 934               Standard_Real aSubDim = Abs(CDim[1] - CDim[0]);
 935
 936               if (myType == Vinert)
 937               {
 938                 ErrorU = ErrUL->Value(iLS);
 939
 940                 CIxx[1] = Abs(CIxx[1] - CIxx[0]);
 941                 CIyy[1] = Abs(CIyy[1] - CIyy[0]);
 942                 CIzz[1] = Abs(CIzz[1] - CIzz[0]);
 943
 944                 #ifndef IS_MIN_DIM
 945                 aSubDim = CIxx[1] + CIyy[1] + CIzz[1];
 946                 #endif
 947
 948                 ErrL->Value(iLS) = add(mult(aSubDim, lr), ErrorU);
 949               }
 950               else
 951               {
 952                 ErrL->Value(iLS) = add(mult(aSubDim, lr), ErrUL->Value(iLS));
 953               }
 954             }
 955
 956             aLI.Ix = mult(CIx, lr);
 957             aLI.Iy = mult(CIy, lr);
 958             aLI.Iz = mult(CIz, lr);
 959
 960             aLI.Ixy = mult(CIxy, lr);
 961             aLI.Ixz = mult(CIxz, lr);
 962             aLI.Iyz = mult(CIyz, lr);
 963           }//for: (kL)iLS
 964
 965           // Calculate/correct epsilon of computation by current value of dim
 966           // That is need for not spend time for
 967           if (JL == iLSubEnd)
 968           {
 969             kLEnd = 2;
 970
 971             Standard_Real DDim = 0.0;
 972             for (i = 1; i <= JL; ++i)
 973               DDim += anInertiaL->Value(i).Mass;
 974
 975             #ifndef IS_MIN_DIM
 976             {
 977               if (myType == Vinert)
 978               {
 979                 Standard_Real DIxx = 0.0, DIyy = 0.0, DIzz = 0.0;
 980                 for (i = 1; i <= JL; ++i)
 981                 {
 982                   const BRepGProp_Gauss::Inertia& aLocalL =
 983                     anInertiaL->Value(i);
 984
 985                   DIxx += aLocalL.Ixx;
 986                   DIyy += aLocalL.Iyy;
 987                   DIzz += aLocalL.Izz;
 988                 }
 989
 990                 DDim = Abs(DIxx) + Abs(DIyy) + Abs(DIzz);
 991               }
 992             }
 993             #endif
 994
 995             DDim = Abs(DDim * anEpsilon);
 996
 997             if (DDim > Eps)
 998             {
 999               Eps  = DDim;
1000               EpsL = 0.9 * Eps;
1001             }
1002           }
1003           if (kLEnd == 2)
1004           {
1005             ErrorL = ErrL->Value(ErrL->Max());
1006           }
1007       } while ( (ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1 );
1008
1009       for ( i = 1; i <= JL; i++ )
1010         addAndRestoreInertia(anInertiaL->Value(i), anInertia);
1011
1012       ErrorLMax = Max(ErrorLMax, ErrorL);
1013     }
1014
1015     if (isNaturalRestriction)
1016       break;
1017
1018     theDomain.Next();
1019   }
1020
1021   if (myType == Vinert)
1022     convert(anInertia, theCoeff, theIsByPoint, theOutGravityCenter, theOutInertia, theOutMass);
1023   else
1024     convert(anInertia, theOutGravityCenter, theOutInertia, theOutMass);
1025
1026   if (iGLEnd == 2)
1027   {
1028     if (theOutMass != 0.0)
1029     {
1030       Eps = ErrorLMax / Abs(theOutMass);
1031
1032       #ifndef IS_MIN_DIM
1033       {
1034         if (myType == Vinert)
1035           Eps = ErrorLMax / (Abs(anInertia.Ixx) +
1036                              Abs(anInertia.Iyy) +
1037                              Abs(anInertia.Izz));
1038       }
1039       #endif
1040     }
1041     else
1042     {
1043       Eps = 0.0;
1044     }
1045   }
1046   else
1047   {
1048     Eps = anEpsilon;
1049   }
1050
1051   return Eps;
1052 }
1053
1054 //=======================================================================
1055 //function : Compute
1056 //purpose  :
1057 //=======================================================================
1058 Standard_Real BRepGProp_Gauss::Compute(BRepGProp_Face&     theSurface,
1059                                        BRepGProp_Domain&   theDomain,
1060                                        const gp_Pnt&       theLocation,
1061                                        const Standard_Real theEps,
1062                                        Standard_Real&      theOutMass,
1063                                        gp_Pnt&             theOutGravityCenter,
1064                                        gp_Mat&             theOutInertia)
1065 {
1066   Standard_ASSERT_RAISE(myType == Sinert, "BRepGProp_Gauss: Incorrect type");
1067
1068   return Compute(theSurface,
1069                  theDomain,
1070                  theLocation,
1071                  theEps,
1072                  NULL,
1073                  Standard_True,
1074                  theOutMass,
1075                  theOutGravityCenter,
1076                  theOutInertia);
1077 }
1078
1079 //=======================================================================
1080 //function : Compute
1081 //purpose  :
1082 //=======================================================================
1083 void BRepGProp_Gauss::Compute(BRepGProp_Face&   theSurface,
1084                               BRepGProp_Domain& theDomain,
1085                               const gp_Pnt&     theLocation,
1086                               Standard_Real&    theOutMass,
1087                               gp_Pnt&           theOutGravityCenter,
1088                               gp_Mat&           theOutInertia)
1089 {
1090   Standard_ASSERT_RAISE(myType == Sinert, "BRepGProp_Gauss: Incorrect type");
1091
1092   Standard_Real u1, u2, v1, v2;
1093   theSurface.Bounds (u1, u2, v1, v2);
1094   checkBounds(u1, u2, v1, v2);
1095
1096   const Standard_Integer NbUGaussgp_Pnts =
1097     Min(theSurface.UIntegrationOrder(), math::GaussPointsMax());
1098
1099   const Standard_Integer NbVGaussgp_Pnts =
1100     Min(theSurface.VIntegrationOrder(), math::GaussPointsMax());
1101
1102   const Standard_Integer NbGaussgp_Pnts =
1103     Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts);
1104
1105   // Number of Gauss points for the integration on the face
1106   math_Vector GaussSPV (1, NbGaussgp_Pnts);
1107   math_Vector GaussSWV (1, NbGaussgp_Pnts);
1108   math::GaussPoints (NbGaussgp_Pnts, GaussSPV);
1109   math::GaussWeights(NbGaussgp_Pnts, GaussSWV);
1110
1111   BRepGProp_Gauss::Inertia anInertia;
1112   while (theDomain.More())
1113   {
1114     theSurface.Load(theDomain.Value());
1115
1116     Standard_Integer NbCGaussgp_Pnts =
1117       Min(theSurface.IntegrationOrder(), math::GaussPointsMax());
1118
1119     NbCGaussgp_Pnts = Max(NbCGaussgp_Pnts, NbGaussgp_Pnts);
1120
1121     math_Vector GaussCP(1, NbCGaussgp_Pnts);
1122     math_Vector GaussCW(1, NbCGaussgp_Pnts);
1123     math::GaussPoints (NbCGaussgp_Pnts, GaussCP);
1124     math::GaussWeights(NbCGaussgp_Pnts, GaussCW);
1125
1126
1127     const Standard_Real l1 = theSurface.FirstParameter();
1128     const Standard_Real l2 = theSurface.LastParameter ();
1129     const Standard_Real lm = 0.5 * (l2 + l1);
1130     const Standard_Real lr = 0.5 * (l2 - l1);
1131
1132     BRepGProp_Gauss::Inertia aCInertia;
1133     for (Standard_Integer i = 1; i <= NbCGaussgp_Pnts; ++i)
1134     {
1135       const Standard_Real l = lm + lr * GaussCP(i);
1136
1137       gp_Pnt2d Puv;
1138       gp_Vec2d Vuv;
1139       theSurface.D12d(l, Puv, Vuv);
1140
1141       const Standard_Real v = Puv.Y();
1142       u2  = Puv.X();
1143
1144       const Standard_Real Dul = Vuv.Y() * GaussCW(i);
1145       const Standard_Real um  = 0.5 * (u2 + u1);
1146       const Standard_Real ur  = 0.5 * (u2 - u1);
1147
1148       BRepGProp_Gauss::Inertia aLocalInertia;
1149       for (Standard_Integer j = 1; j <= NbGaussgp_Pnts; ++j)
1150       {
1151         const Standard_Real u = add(um, mult(ur, GaussSPV(j)));
1152         const Standard_Real aWeight = Dul * GaussSWV(j);
1153
1154         gp_Pnt aPoint;
1155         gp_Vec aNormal;
1156         theSurface.Normal (u, v, aPoint, aNormal);
1157
1158         computeSInertiaOfElementaryPart(aPoint, aNormal, theLocation, aWeight, aLocalInertia);
1159       }
1160
1161       multAndRestoreInertia(ur, aLocalInertia);
1162       addAndRestoreInertia (aLocalInertia, aCInertia);
1163     }
1164
1165     multAndRestoreInertia(lr, aCInertia);
1166     addAndRestoreInertia (aCInertia, anInertia);
1167
1168     theDomain.Next();
1169   }
1170
1171   convert(anInertia, theOutGravityCenter, theOutInertia, theOutMass);
1172 }
1173
1174 //=======================================================================
1175 //function : Compute
1176 //purpose  :
1177 //=======================================================================
1178 void BRepGProp_Gauss::Compute(BRepGProp_Face&         theSurface,
1179                               BRepGProp_Domain&       theDomain,
1180                               const gp_Pnt&           theLocation,
1181                               const Standard_Real     theCoeff[],
1182                               const Standard_Boolean  theIsByPoint,
1183                               Standard_Real&          theOutMass,
1184                               gp_Pnt&                 theOutGravityCenter,
1185                               gp_Mat&                 theOutInertia)
1186 {
1187   Standard_ASSERT_RAISE(myType == Vinert, "BRepGProp_Gauss: Incorrect type");
1188
1189   Standard_Real u1, v1, u2, v2;
1190   theSurface.Bounds (u1, u2, v1, v2);
1191   checkBounds(u1, u2, v1, v2);
1192
1193   Standard_Real _u2 = u2;  //OCC104
1194
1195   BRepGProp_Gauss::Inertia anInertia;
1196   while (theDomain.More())
1197   {
1198     theSurface.Load(theDomain.Value());
1199
1200     const Standard_Integer aVNbCGaussgp_Pnts =
1201       theSurface.VIntegrationOrder();
1202
1203     const Standard_Integer aNbGaussgp_Pnts =
1204       Min( Max(theSurface.IntegrationOrder(), aVNbCGaussgp_Pnts), math::GaussPointsMax() );
1205
1206     math_Vector GaussP(1, aNbGaussgp_Pnts);
1207     math_Vector GaussW(1, aNbGaussgp_Pnts);
1208     math::GaussPoints (aNbGaussgp_Pnts, GaussP);
1209     math::GaussWeights(aNbGaussgp_Pnts, GaussW);
1210
1211     const Standard_Real l1 = theSurface.FirstParameter();
1212     const Standard_Real l2 = theSurface.LastParameter();
1213     const Standard_Real lm = 0.5 * (l2 + l1);
1214     const Standard_Real lr = 0.5 * (l2 - l1);
1215
1216     BRepGProp_Gauss::Inertia aCInertia;
1217     for (Standard_Integer i = 1; i <= aNbGaussgp_Pnts; ++i)
1218     {
1219       const Standard_Real l = lm + lr * GaussP(i);
1220
1221       gp_Pnt2d Puv;
1222       gp_Vec2d Vuv;
1223
1224       theSurface.D12d(l, Puv, Vuv);
1225
1226       u2  = Puv.X();
1227       u2 = Min( Max(u1, u2), _u2 ); // OCC104
1228       const Standard_Real v = Min(Max(Puv.Y(), v1), v2);
1229
1230       const Standard_Real Dul = Vuv.Y() * GaussW(i);
1231       const Standard_Real um  = 0.5 * (u2 + u1);
1232       const Standard_Real ur  = 0.5 * (u2 - u1);
1233
1234       BRepGProp_Gauss::Inertia aLocalInertia;
1235       for (Standard_Integer j = 1; j <= aNbGaussgp_Pnts; ++j)
1236       {
1237         const Standard_Real u = um + ur * GaussP(j);
1238         const Standard_Real aWeight = Dul * GaussW(j);
1239
1240         gp_Pnt aPoint;
1241         gp_Vec aNormal;
1242
1243         theSurface.Normal(u, v, aPoint, aNormal);
1244
1245         computeVInertiaOfElementaryPart(
1246           aPoint,
1247           aNormal,
1248           theLocation,
1249           aWeight,
1250           theCoeff,
1251           theIsByPoint,
1252           aLocalInertia);
1253       }
1254
1255       multAndRestoreInertia(ur, aLocalInertia);
1256       addAndRestoreInertia (aLocalInertia, aCInertia);
1257     }
1258
1259     multAndRestoreInertia(lr,        aCInertia);
1260     addAndRestoreInertia (aCInertia, anInertia);
1261
1262     theDomain.Next();
1263   }
1264
1265   convert(anInertia, theCoeff, theIsByPoint, theOutGravityCenter, theOutInertia, theOutMass);
1266 }
1267
1268 //=======================================================================
1269 //function : Compute
1270 //purpose  :
1271 //=======================================================================
1272 void BRepGProp_Gauss::Compute(const BRepGProp_Face&  theSurface,
1273                               const gp_Pnt&          theLocation,
1274                               const Standard_Real    theCoeff[],
1275                               const Standard_Boolean theIsByPoint,
1276                               Standard_Real&         theOutMass,
1277                               gp_Pnt&                theOutGravityCenter,
1278                               gp_Mat&                theOutInertia)
1279 {
1280   Standard_Real LowerU, UpperU, LowerV, UpperV;
1281   theSurface.Bounds(LowerU, UpperU, LowerV, UpperV);
1282   checkBounds(LowerU, UpperU, LowerV, UpperV);
1283
1284   const Standard_Integer UOrder =
1285     Min(theSurface.UIntegrationOrder(), math::GaussPointsMax());
1286   const Standard_Integer VOrder =
1287     Min(theSurface.VIntegrationOrder(), math::GaussPointsMax());
1288
1289   // Gauss points and weights
1290   math_Vector GaussPU(1, UOrder);
1291   math_Vector GaussWU(1, UOrder);
1292   math_Vector GaussPV(1, VOrder);
1293   math_Vector GaussWV(1, VOrder);
1294
1295   math::GaussPoints (UOrder, GaussPU);
1296   math::GaussWeights(UOrder, GaussWU);
1297   math::GaussPoints (VOrder, GaussPV);
1298   math::GaussWeights(VOrder, GaussWV);
1299
1300   const Standard_Real um = 0.5 * add(UpperU,  LowerU);
1301   const Standard_Real vm = 0.5 * add(UpperV,  LowerV);
1302   Standard_Real ur = 0.5 * add(UpperU, -LowerU);
1303   Standard_Real vr = 0.5 * add(UpperV, -LowerV);
1304
1305   gp_Pnt aPoint;
1306   gp_Vec aNormal;
1307
1308   BRepGProp_Gauss::Inertia anInertia;
1309   for (Standard_Integer j = 1; j <= VOrder; ++j)
1310   {
1311     BRepGProp_Gauss::Inertia anInertiaOfElementaryPart;
1312     const Standard_Real v = add(vm, mult(vr, GaussPV(j)));
1313
1314     for (Standard_Integer i = 1; i <= UOrder; ++i)
1315     {
1316       const Standard_Real aWeight = GaussWU(i);
1317       const Standard_Real u = add(um, mult(ur, GaussPU (i)));
1318       theSurface.Normal(u, v, aPoint, aNormal);
1319
1320       if (myType == Vinert)
1321       {
1322         computeVInertiaOfElementaryPart(
1323           aPoint,
1324           aNormal,
1325           theLocation,
1326           aWeight,
1327           theCoeff,
1328           theIsByPoint,
1329           anInertiaOfElementaryPart);
1330       }
1331       else // Sinert
1332       {
1333         computeSInertiaOfElementaryPart(
1334           aPoint,
1335           aNormal,
1336           theLocation,
1337           aWeight,
1338           anInertiaOfElementaryPart);
1339       }
1340     }
1341
1342     multAndRestoreInertia(GaussWV(j), anInertiaOfElementaryPart);
1343     addAndRestoreInertia (anInertiaOfElementaryPart,  anInertia);
1344   }
1345   vr            = mult(vr, ur);
1346   anInertia.Ixx = mult(vr, anInertia.Ixx);
1347   anInertia.Iyy = mult(vr, anInertia.Iyy);
1348   anInertia.Izz = mult(vr, anInertia.Izz);
1349   anInertia.Ixy = mult(vr, anInertia.Ixy);
1350   anInertia.Ixz = mult(vr, anInertia.Ixz);
1351   anInertia.Iyz = mult(vr, anInertia.Iyz);
1352
1353   if (myType == Vinert)
1354   {
1355     convert(anInertia, theCoeff, theIsByPoint, theOutGravityCenter, theOutInertia, theOutMass);
1356   }
1357   else // Sinert
1358   {
1359     convert(anInertia, theOutGravityCenter, theOutInertia, theOutMass);
1360   }
1361
1362   theOutMass *= vr;
1363 }
1364
1365 //=======================================================================
1366 //function : Compute
1367 //purpose  :
1368 //=======================================================================
1369 void BRepGProp_Gauss::Compute(const BRepGProp_Face& theSurface,
1370                               const gp_Pnt&         theLocation,
1371                               Standard_Real&        theOutMass,
1372                               gp_Pnt&               theOutGravityCenter,
1373                               gp_Mat&               theOutInertia)
1374 {
1375   Standard_ASSERT_RAISE(myType == Sinert, "BRepGProp_Gauss: Incorrect type");
1376
1377   Compute(theSurface,
1378           theLocation,
1379           NULL,
1380           Standard_True,
1381           theOutMass,
1382           theOutGravityCenter,
1383           theOutInertia);
1384 }