0025630: Possible memory leaks in BRepGProp_Vinert and BRepGProp_Sinert
authorazn <azn@opencascade.com>
Thu, 19 Mar 2015 12:50:09 +0000 (15:50 +0300)
committerbugmaster <bugmaster@opencascade.com>
Thu, 19 Mar 2015 14:08:13 +0000 (17:08 +0300)
Code refactoring of BRepGProp_Sinert and BRepGProp_Vinert classes.
- All static variables have been removed.
- Common functionality connected with Gauss integration has beem moved from BRepGProp_Sinert and BRepGProp_Vinert classes to the new BRepGProp_Gauss class.

Slight changes in the comments.

Fix compilation error.

Fix Sinert errors. Rebased on new master.

Elimination of constant conditional expression warnings.

Small fix in comment.

src/BRepGProp/BRepGProp_Gauss.cxx [new file with mode: 0644]
src/BRepGProp/BRepGProp_Gauss.hxx [new file with mode: 0644]
src/BRepGProp/BRepGProp_Sinert.cxx
src/BRepGProp/BRepGProp_Vinert.cxx
src/BRepGProp/FILES [new file with mode: 0644]

diff --git a/src/BRepGProp/BRepGProp_Gauss.cxx b/src/BRepGProp/BRepGProp_Gauss.cxx
new file mode 100644 (file)
index 0000000..d905f75
--- /dev/null
@@ -0,0 +1,1380 @@
+// Copyright (c) 2008-2015 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#include <math.hxx>
+#include <Precision.hxx>
+#include <TColStd_Array1OfReal.hxx>
+#include <Standard_Assert.hxx>
+#include <BRepGProp_Face.hxx>
+#include <BRepGProp_Domain.hxx>
+#include <BRepGProp_Gauss.hxx>
+
+// If the following is defined the error of algorithm is calculated by static moments
+#define IS_MIN_DIM
+
+namespace
+{
+  // Minimal value of interval's range for computation | minimal value of "dim" | ... 
+  static const Standard_Real EPS_PARAM          = 1.e-12;
+  static const Standard_Real EPS_DIM            = 1.e-30;
+  static const Standard_Real ERROR_ALGEBR_RATIO = 2.0 / 3.0;
+
+  // Maximum of GaussPoints on a subinterval and maximum of subintervals
+  static const Standard_Integer GPM        = math::GaussPointsMax();
+  static const Standard_Integer SUBS_POWER = 32;
+  static const Standard_Integer SM         = SUBS_POWER * GPM + 1;
+
+  // Auxiliary inner functions to perform arithmetic operations.
+  static Standard_Real Add(const Standard_Real theA, const Standard_Real theB)
+  {
+    return theA + theB;
+  }
+
+  static Standard_Real AddInf(const Standard_Real theA, const Standard_Real theB)
+  {
+    if (Precision::IsPositiveInfinite(theA))
+    {
+      if (Precision::IsNegativeInfinite(theB))
+        return 0.0;
+      else
+        return Precision::Infinite();
+    }
+
+    if (Precision::IsPositiveInfinite(theB))
+    {
+      if (Precision::IsNegativeInfinite(theA))
+        return 0.0;
+      else
+        return Precision::Infinite();
+    }
+
+    if (Precision::IsNegativeInfinite(theA))
+    {
+      if (Precision::IsPositiveInfinite(theB))
+        return 0.0;
+      else
+        return -Precision::Infinite();
+    }
+
+    if (Precision::IsNegativeInfinite(theB))
+    {
+      if (Precision::IsPositiveInfinite(theA))
+        return 0.0;
+      else
+        return -Precision::Infinite();
+    }
+
+    return theA + theB;
+  }
+
+  static Standard_Real Mult(const Standard_Real theA, const Standard_Real theB)
+  {
+    return theA * theB;
+  }
+
+  static Standard_Real MultInf(const Standard_Real theA, const Standard_Real theB)
+  {
+    if ((theA == 0.0) || (theB == 0.0))  //strictly zerro (without any tolerances)
+      return 0.0;
+
+    if (Precision::IsPositiveInfinite(theA))
+    {
+      if (theB < 0.0)
+        return -Precision::Infinite();
+      else
+        return Precision::Infinite();
+    }
+
+    if (Precision::IsPositiveInfinite(theB))
+    {
+      if (theA < 0.0)
+        return -Precision::Infinite();
+      else
+        return Precision::Infinite();
+    }
+
+    if (Precision::IsNegativeInfinite(theA))
+    {
+      if (theB < 0.0)
+        return +Precision::Infinite();
+      else
+        return -Precision::Infinite();
+    }
+
+    if (Precision::IsNegativeInfinite(theB))
+    {
+      if (theA < 0.0)
+        return +Precision::Infinite();
+      else
+        return -Precision::Infinite();
+    }
+
+    return theA * theB;
+  }
+}
+
+//=======================================================================
+//function : BRepGProp_Gauss::Inert::Inert
+//purpose  : Constructor
+//=======================================================================
+BRepGProp_Gauss::Inertia::Inertia()
+: Mass(0.0),
+  Ix  (0.0),
+  Iy  (0.0),
+  Iz  (0.0),
+  Ixx (0.0),
+  Iyy (0.0),
+  Izz (0.0),
+  Ixy (0.0),
+  Ixz (0.0),
+  Iyz (0.0)
+{
+}
+
+//=======================================================================
+//function : Inertia::Reset
+//purpose  : Zeroes all values.
+//=======================================================================
+void BRepGProp_Gauss::Inertia::Reset()
+{
+  memset(reinterpret_cast<void*>(this), 0, sizeof(BRepGProp_Gauss::Inertia));
+}
+
+//=======================================================================
+//function : BRepGProp_Gauss
+//purpose  : Constructor
+//=======================================================================
+BRepGProp_Gauss::BRepGProp_Gauss(const BRepGProp_GaussType theType)
+: myType(theType)
+{
+  add  = (::Add );
+  mult = (::Mult);
+}
+
+//=======================================================================
+//function : MaxSubs
+//purpose  : 
+//=======================================================================
+Standard_Integer BRepGProp_Gauss::MaxSubs(const Standard_Integer theN,
+                                          const Standard_Integer theCoeff)
+{
+  return IntegerLast() / theCoeff < theN ?
+         IntegerLast() : theN * theCoeff + 1;
+}
+
+//=======================================================================
+//function : Init
+//purpose  : 
+//=======================================================================
+void BRepGProp_Gauss::Init(Handle_Vector&         theOutVec,
+                           const Standard_Real    theValue,
+                           const Standard_Integer theFirst,
+                           const Standard_Integer theLast)
+{
+  if(theLast - theFirst == 0)
+  {
+    theOutVec->Init(theValue);
+  }
+  else
+  {
+    for (Standard_Integer i = theFirst; i <= theLast; ++i)
+      theOutVec->Value(i) = theValue;
+  }
+}
+
+//=======================================================================
+//function : InitMass
+//purpose  : 
+//=======================================================================
+void BRepGProp_Gauss::InitMass(const Standard_Real    theValue,
+                               const Standard_Integer theFirst,
+                               const Standard_Integer theLast,
+                               InertiaArray&          theArray)
+{
+  if (theArray.IsNull())
+    return;
+
+  Standard_Integer aFirst = theFirst;
+  Standard_Integer aLast  = theLast;
+
+  if (theLast - theFirst == 0)
+  {
+    aFirst = theArray->Lower();
+    aLast  = theArray->Upper();
+  }
+
+  for (Standard_Integer i = aFirst; i <= aLast; ++i)
+      theArray->ChangeValue(i).Mass = theValue;
+}
+
+//=======================================================================
+//function : FillIntervalBounds
+//purpose  : 
+//=======================================================================
+Standard_Integer BRepGProp_Gauss::FillIntervalBounds(
+  const Standard_Real         theA,
+  const Standard_Real         theB,
+  const TColStd_Array1OfReal& theKnots,
+  const Standard_Integer      theNumSubs,
+  InertiaArray&               theInerts,
+  Handle_Vector&              theParam1,
+  Handle_Vector&              theParam2,
+  Handle_Vector&              theError,
+  Handle_Vector&              theCommonError)
+{
+  const Standard_Integer aSize =
+    Max(theKnots.Upper(), MaxSubs(theKnots.Upper() - 1, theNumSubs));
+
+  if (aSize - 1 > theParam1->Upper())
+  {
+    theInerts = new NCollection_Array1<Inertia>(1, aSize);
+    theParam1 = new math_Vector(1, aSize);
+    theParam2 = new math_Vector(1, aSize);
+    theError  = new math_Vector(1, aSize, 0.0);
+
+    if (theCommonError.IsNull() == Standard_False)
+      theCommonError = new math_Vector(1, aSize, 0.0);
+  }
+
+  Standard_Integer j = 1, k = 1;
+  theParam1->Value(j++) = theA;
+
+  const Standard_Integer aLength = theKnots.Upper();
+  for (Standard_Integer i = 1; i <= aLength; ++i)
+  {
+    const Standard_Real kn = theKnots(i);
+    if (theA < kn)
+    {
+      if (kn < theB) 
+      {
+        theParam1->Value(j++) = kn;
+        theParam2->Value(k++) = kn; 
+      }
+      else
+        break;
+    }
+  }
+
+  theParam2->Value(k) = theB;
+  return k;
+}
+
+
+
+//=======================================================================
+//function : computeVInertiaOfElementaryPart
+//purpose  : 
+//=======================================================================
+void BRepGProp_Gauss::computeVInertiaOfElementaryPart(
+  const gp_Pnt&             thePoint,
+  const gp_Vec&             theNormal,
+  const gp_Pnt&             theLocation,
+  const Standard_Real       theWeight,
+  const Standard_Real       theCoeff[],
+  const Standard_Boolean    theIsByPoint,
+  BRepGProp_Gauss::Inertia& theOutInertia)
+{
+  Standard_Real x = thePoint.X() - theLocation.X();
+  Standard_Real y = thePoint.Y() - theLocation.Y();
+  Standard_Real z = thePoint.Z() - theLocation.Z();
+
+  const Standard_Real xn = theNormal.X() * theWeight;
+  const Standard_Real yn = theNormal.Y() * theWeight;
+  const Standard_Real zn = theNormal.Z() * theWeight;
+
+  if (theIsByPoint)
+  {
+    /////////////////////                        ///////////////////////
+    //    OFV code     //                        //    Initial code   //
+    /////////////////////                        ///////////////////////
+    // modified by APO
+
+    Standard_Real dv = x * xn + y * yn + z * zn; //xyz  = x * y * z;
+    theOutInertia.Mass += dv / 3.0;              //Ixyi += zn * xyz; 
+    theOutInertia.Ix   += 0.25 * x * dv;         //Iyzi += xn * xyz;
+    theOutInertia.Iy   += 0.25 * y * dv;         //Ixzi += yn * xyz; 
+    theOutInertia.Iz   += 0.25 * z * dv;         //xi = x * x * x * xn / 3.0;
+    x -= theCoeff[0];                            //yi = y * y * y * yn / 3.0;
+    y -= theCoeff[1];                            //zi = z * z * z * zn / 3.0;
+    z -= theCoeff[2];                            //Ixxi += (yi + zi);
+    dv *= 0.2;                                   //Iyyi += (xi + zi);
+    theOutInertia.Ixy -= x * y * dv;             //Izzi += (xi + yi);
+    theOutInertia.Iyz -= y * z * dv;             //x -= Coeff[0];
+    theOutInertia.Ixz -= x * z * dv;             //y -= Coeff[1];
+    x *= x;                                      //z -= Coeff[2];
+    y *= y;                                      //dv = x * xn + y * yn + z * zn;
+    z *= z;                                      //dvi +=  dv; 
+    theOutInertia.Ixx += (y + z) * dv;           //Ixi += x * dv;
+    theOutInertia.Iyy += (x + z) * dv;           //Iyi += y * dv;
+    theOutInertia.Izz += (x + y) * dv;           //Izi += z * dv;
+  }
+  else
+  { // By plane
+    const Standard_Real s = xn * theCoeff[0] + yn * theCoeff[1] + zn * theCoeff[2];
+
+    Standard_Real d1  = theCoeff[0] * x + theCoeff[1] * y + theCoeff[2] * z - theCoeff[3];
+    Standard_Real d2  = d1 * d1;
+    Standard_Real d3  = d1 * d2 / 3.0;
+    Standard_Real dv  = s  * d1;
+
+    theOutInertia.Mass += dv;
+    theOutInertia.Ix   += (x - (theCoeff[0] * d1 * 0.5)) * dv;
+    theOutInertia.Iy   += (y - (theCoeff[1] * d1 * 0.5)) * dv;
+    theOutInertia.Iz   += (z - (theCoeff[2] * d1 * 0.5)) * dv;
+
+    const Standard_Real px = x - theCoeff[0] * d1;
+    const Standard_Real py = y - theCoeff[1] * d1;
+    const Standard_Real pz = z - theCoeff[2] * d1;
+
+    x = px * px * d1 + px * theCoeff[0] * d2 + theCoeff[0] * theCoeff[0] * d3;
+    y = py * py * d1 + py * theCoeff[1] * d2 + theCoeff[1] * theCoeff[1] * d3;
+    z = pz * pz * d1 + pz * theCoeff[2] * d2 + theCoeff[2] * theCoeff[2] * d3;
+
+    theOutInertia.Ixx += (y + z) * s;
+    theOutInertia.Iyy += (x + z) * s;
+    theOutInertia.Izz += (x + y) * s;
+
+    d2 *= 0.5;
+    x = (py * pz * d1) + (py * theCoeff[2] * d2) + (pz * theCoeff[1] * d2) + (theCoeff[1] * theCoeff[2] * d3);
+    y = (px * pz * d1) + (pz * theCoeff[0] * d2) + (px * theCoeff[2] * d2) + (theCoeff[0] * theCoeff[2] * d3);
+    z = (px * py * d1) + (px * theCoeff[1] * d2) + (py * theCoeff[0] * d2) + (theCoeff[0] * theCoeff[1] * d3);
+
+    theOutInertia.Ixy -= z * s;
+    theOutInertia.Iyz -= x * s;
+    theOutInertia.Ixz -= y * s;
+  }
+}
+
+//=======================================================================
+//function : computeSInertiaOfElementaryPart
+//purpose  : 
+//=======================================================================
+void BRepGProp_Gauss::computeSInertiaOfElementaryPart(
+  const gp_Pnt&             thePoint,
+  const gp_Vec&             theNormal,
+  const gp_Pnt&             theLocation,
+  const Standard_Real       theWeight,
+  BRepGProp_Gauss::Inertia& theOutInertia)
+{
+  // ds - Jacobien (x, y, z) -> (u, v) = ||n||
+  const Standard_Real ds = mult(theNormal.Magnitude(), theWeight);
+  const Standard_Real x = add(thePoint.X(), -theLocation.X());
+  const Standard_Real y = add(thePoint.Y(), -theLocation.Y());
+  const Standard_Real z = add(thePoint.Z(), -theLocation.Z());
+
+  theOutInertia.Mass = add(theOutInertia.Mass, ds);
+
+  const Standard_Real XdS = mult(x, ds);
+  const Standard_Real YdS = mult(y, ds);
+  const Standard_Real ZdS = mult(z, ds);
+
+  theOutInertia.Ix  = add(theOutInertia.Ix,  XdS);
+  theOutInertia.Iy  = add(theOutInertia.Iy,  YdS);
+  theOutInertia.Iz  = add(theOutInertia.Iz,  ZdS);
+  theOutInertia.Ixy = add(theOutInertia.Ixy, mult(x, YdS));
+  theOutInertia.Iyz = add(theOutInertia.Iyz, mult(y, ZdS));
+  theOutInertia.Ixz = add(theOutInertia.Ixz, mult(x, ZdS));
+
+  const Standard_Real XXdS = mult(x, XdS);
+  const Standard_Real YYdS = mult(y, YdS);
+  const Standard_Real ZZdS = mult(z, ZdS);
+
+  theOutInertia.Ixx = add(theOutInertia.Ixx, add(YYdS, ZZdS));
+  theOutInertia.Iyy = add(theOutInertia.Iyy, add(XXdS, ZZdS));
+  theOutInertia.Izz = add(theOutInertia.Izz, add(XXdS, YYdS));
+}
+
+//=======================================================================
+//function : checkBounds
+//purpose  : 
+//=======================================================================
+void BRepGProp_Gauss::checkBounds(const Standard_Real theU1,
+                                  const Standard_Real theU2,
+                                  const Standard_Real theV1,
+                                  const Standard_Real theV2)
+{
+  if (Precision::IsInfinite(theU1) || Precision::IsInfinite(theU2) ||
+      Precision::IsInfinite(theV1) || Precision::IsInfinite(theV2))
+  {
+    add  = (::AddInf);
+    mult = (::MultInf);
+  }
+}
+
+//=======================================================================
+//function : addAndRestoreInertia
+//purpose  : 
+//=======================================================================
+void BRepGProp_Gauss::addAndRestoreInertia(
+  const BRepGProp_Gauss::Inertia& theInInertia,
+  BRepGProp_Gauss::Inertia&       theOutInertia)
+{
+  theOutInertia.Mass = add(theOutInertia.Mass, theInInertia.Mass);
+  theOutInertia.Ix   = add(theOutInertia.Ix,   theInInertia.Ix);
+  theOutInertia.Iy   = add(theOutInertia.Iy,   theInInertia.Iy);
+  theOutInertia.Iz   = add(theOutInertia.Iz,   theInInertia.Iz);
+  theOutInertia.Ixx  = add(theOutInertia.Ixx,  theInInertia.Ixx);
+  theOutInertia.Iyy  = add(theOutInertia.Iyy,  theInInertia.Iyy);
+  theOutInertia.Izz  = add(theOutInertia.Izz,  theInInertia.Izz);
+  theOutInertia.Ixy  = add(theOutInertia.Ixy,  theInInertia.Ixy);
+  theOutInertia.Ixz  = add(theOutInertia.Ixz,  theInInertia.Ixz);
+  theOutInertia.Iyz  = add(theOutInertia.Iyz,  theInInertia.Iyz);
+}
+
+//=======================================================================
+//function : multAndRestoreInertia
+//purpose  : 
+//=======================================================================
+void BRepGProp_Gauss::multAndRestoreInertia(
+  const Standard_Real       theValue,
+  BRepGProp_Gauss::Inertia& theInOutInertia)
+{
+  theInOutInertia.Mass = mult(theInOutInertia.Mass, theValue);
+  theInOutInertia.Ix   = mult(theInOutInertia.Ix,   theValue);
+  theInOutInertia.Iy   = mult(theInOutInertia.Iy,   theValue);
+  theInOutInertia.Iz   = mult(theInOutInertia.Iz,   theValue);
+  theInOutInertia.Ixx  = mult(theInOutInertia.Ixx,  theValue);
+  theInOutInertia.Iyy  = mult(theInOutInertia.Iyy,  theValue);
+  theInOutInertia.Izz  = mult(theInOutInertia.Izz,  theValue);
+  theInOutInertia.Ixy  = mult(theInOutInertia.Ixy,  theValue);
+  theInOutInertia.Ixz  = mult(theInOutInertia.Ixz,  theValue);
+  theInOutInertia.Iyz  = mult(theInOutInertia.Iyz,  theValue);
+}
+
+//=======================================================================
+//function : convert
+//purpose  : 
+//=======================================================================
+void BRepGProp_Gauss::convert(const BRepGProp_Gauss::Inertia& theInertia,
+                              gp_Pnt&                         theOutGravityCenter,
+                              gp_Mat&                         theOutMatrixOfInertia,
+                              Standard_Real&                  theOutMass)
+{
+  if (Abs(theInertia.Mass) >= EPS_DIM)
+  {
+    const Standard_Real anInvMass = 1.0 / theInertia.Mass;
+    theOutGravityCenter.SetX(theInertia.Ix * anInvMass);
+    theOutGravityCenter.SetY(theInertia.Iy * anInvMass);
+    theOutGravityCenter.SetZ(theInertia.Iz * anInvMass);
+
+    theOutMass = theInertia.Mass;
+  }
+  else
+  {
+    theOutMass = 0.0;
+    theOutGravityCenter.SetCoord(0.0, 0.0, 0.0);
+  }
+
+  theOutMatrixOfInertia = gp_Mat(
+    gp_XYZ ( theInertia.Ixx, -theInertia.Ixy, -theInertia.Ixz),
+    gp_XYZ (-theInertia.Ixy,  theInertia.Iyy, -theInertia.Iyz),
+    gp_XYZ (-theInertia.Ixz, -theInertia.Iyz,  theInertia.Izz));
+}
+
+//=======================================================================
+//function : convert
+//purpose  : 
+//=======================================================================
+void BRepGProp_Gauss::convert(const BRepGProp_Gauss::Inertia& theInertia,
+                              const Standard_Real             theCoeff[],
+                              const Standard_Boolean          theIsByPoint,
+                              gp_Pnt&                         theOutGravityCenter,
+                              gp_Mat&                         theOutMatrixOfInertia,
+                              Standard_Real&                  theOutMass)
+{
+  convert(theInertia, theOutGravityCenter, theOutMatrixOfInertia, theOutMass);
+  if (Abs(theInertia.Mass) >= EPS_DIM && theIsByPoint)
+  {
+    const Standard_Real anInvMass = 1.0 / theInertia.Mass;
+    if (theIsByPoint == Standard_True)
+    {
+      theOutGravityCenter.SetX(theCoeff[0] + theInertia.Ix * anInvMass);
+      theOutGravityCenter.SetY(theCoeff[1] + theInertia.Iy * anInvMass);
+      theOutGravityCenter.SetZ(theCoeff[2] + theInertia.Iz * anInvMass);
+    }
+    else
+    {
+      theOutGravityCenter.SetX(theInertia.Ix * anInvMass);
+      theOutGravityCenter.SetY(theInertia.Iy * anInvMass);
+      theOutGravityCenter.SetZ(theInertia.Iz * anInvMass);
+    }
+
+    theOutMass = theInertia.Mass;
+  }
+  else
+  {
+    theOutMass = 0.0;
+    theOutGravityCenter.SetCoord(0.0, 0.0, 0.0);
+  }
+
+  theOutMatrixOfInertia = gp_Mat(
+    gp_XYZ (theInertia.Ixx, theInertia.Ixy, theInertia.Ixz),
+    gp_XYZ (theInertia.Ixy, theInertia.Iyy, theInertia.Iyz),
+    gp_XYZ (theInertia.Ixz, theInertia.Iyz, theInertia.Izz));
+}
+
+//=======================================================================
+//function : Compute
+//purpose  : 
+//=======================================================================
+Standard_Real BRepGProp_Gauss::Compute(
+  BRepGProp_Face&        theSurface,
+  BRepGProp_Domain&      theDomain,
+  const gp_Pnt&          theLocation,
+  const Standard_Real    theEps,
+  const Standard_Real    theCoeff[],
+  const Standard_Boolean theIsByPoint,
+  Standard_Real&         theOutMass,
+  gp_Pnt&                theOutGravityCenter,
+  gp_Mat&                theOutInertia)
+{
+  const Standard_Boolean isErrorCalculation  =
+    ( 0.0 > theEps || theEps < 0.001 ) ? Standard_True : Standard_False;
+  const Standard_Boolean isVerifyComputation =
+    ( 0.0 < theEps && theEps < 0.001 ) ? Standard_True : Standard_False;
+
+  Standard_Real anEpsilon= Abs(theEps);
+
+  BRepGProp_Gauss::Inertia anInertia;
+  InertiaArray anInertiaL = new NCollection_Array1<Inertia>(1, SM);
+  InertiaArray anInertiaU = new NCollection_Array1<Inertia>(1, SM);
+
+  // Prepare Gauss points and weights
+  Handle_Vector LGaussP[2];
+  Handle_Vector LGaussW[2];
+  Handle_Vector UGaussP[2];
+  Handle_Vector UGaussW[2];
+
+  const Standard_Integer aNbGaussPoint =
+    RealToInt(Ceiling(ERROR_ALGEBR_RATIO * GPM));
+
+  LGaussP[0] = new math_Vector(1, GPM);
+  LGaussP[1] = new math_Vector(1, aNbGaussPoint);
+  LGaussW[0] = new math_Vector(1, GPM);
+  LGaussW[1] = new math_Vector(1, aNbGaussPoint);
+
+  UGaussP[0] = new math_Vector(1, GPM);
+  UGaussP[1] = new math_Vector(1, aNbGaussPoint);
+  UGaussW[0] = new math_Vector(1, GPM);
+  UGaussW[1] = new math_Vector(1, aNbGaussPoint);
+
+  Handle_Vector L1 = new math_Vector(1, SM);
+  Handle_Vector L2 = new math_Vector(1, SM);
+  Handle_Vector U1 = new math_Vector(1, SM);
+  Handle_Vector U2 = new math_Vector(1, SM);
+
+  Handle_Vector ErrL  = new math_Vector(1, SM, 0.0);
+  Handle_Vector ErrU  = new math_Vector(1, SM, 0.0);
+  Handle_Vector ErrUL = new math_Vector(1, SM, 0.0);
+
+  // Face parametrization in U and V direction
+  Standard_Real BV1, BV2, BU1, BU2;
+  theSurface.Bounds(BU1, BU2, BV1, BV2);
+  checkBounds(BU1, BU2, BV1, BV2);
+
+  //
+  const Standard_Integer NumSubs = SUBS_POWER;
+  const Standard_Boolean isNaturalRestriction = theSurface.NaturalRestriction();
+
+  Standard_Real CIx, CIy, CIz, CIxy, CIxz, CIyz;
+  Standard_Real CDim[2], CIxx[2], CIyy[2], CIzz[2];
+
+  // Boundary curve parametrization
+  Standard_Real u1 = BU1, u2, l1, l2, lm, lr, l, v;
+
+  // On the boundary curve u-v
+  gp_Pnt2d Puv;
+  gp_Vec2d Vuv;
+  Standard_Real Dul;  // Dul = Du / Dl
+
+  Standard_Integer iLS, iLSubEnd, iGL, iGLEnd, NbLGaussP[2], LRange[2], iL, kL, kLEnd, IL, JL;
+  Standard_Integer i, iUSubEnd, NbUGaussP[2], URange[2], kU, kUEnd, IU, JU;
+  Standard_Integer UMaxSubs, LMaxSubs;
+
+  Standard_Real ErrorU, ErrorL, ErrorLMax = 0.0, Eps = 0.0, EpsL = 0.0, EpsU = 0.0;
+  iGLEnd = isErrorCalculation ? 2 : 1;
+
+  NbUGaussP[0] = theSurface.SIntOrder(anEpsilon);
+  NbUGaussP[1] = RealToInt( Ceiling(ERROR_ALGEBR_RATIO * NbUGaussP[0]) );
+
+  math::GaussPoints (NbUGaussP[0], *UGaussP[0]);
+  math::GaussWeights(NbUGaussP[0], *UGaussW[0]);
+  math::GaussPoints (NbUGaussP[1], *UGaussP[1]);
+  math::GaussWeights(NbUGaussP[1], *UGaussW[1]);
+
+  const Standard_Integer aNbUSubs = theSurface.SUIntSubs();
+  TColStd_Array1OfReal UKnots(1, aNbUSubs + 1);
+  theSurface.UKnots(UKnots);
+
+  while (isNaturalRestriction || theDomain.More())
+  {
+    if (isNaturalRestriction)
+    { 
+      NbLGaussP[0] = Min(2 * NbUGaussP[0], math::GaussPointsMax());
+    }
+    else
+    {
+      theSurface.Load(theDomain.Value());
+      NbLGaussP[0] = theSurface.LIntOrder(anEpsilon);
+    }
+
+    NbLGaussP[1] = RealToInt( Ceiling(ERROR_ALGEBR_RATIO * NbLGaussP[0]) );
+
+    math::GaussPoints (NbLGaussP[0], *LGaussP[0]);
+    math::GaussWeights(NbLGaussP[0], *LGaussW[0]);
+    math::GaussPoints (NbLGaussP[1], *LGaussP[1]);
+    math::GaussWeights(NbLGaussP[1], *LGaussW[1]);
+
+    const Standard_Integer aNbLSubs =
+      isNaturalRestriction ? theSurface.SVIntSubs(): theSurface.LIntSubs();
+    TColStd_Array1OfReal LKnots(1, aNbLSubs + 1);
+
+    if (isNaturalRestriction)
+    {
+      theSurface.VKnots(LKnots);
+      l1 = BV1;
+      l2 = BV2;
+    }
+    else
+    {
+      theSurface.LKnots(LKnots);
+      l1 = theSurface.FirstParameter();
+      l2 = theSurface.LastParameter();
+    }
+    ErrorL = 0.0;
+    kLEnd = 1; JL = 0;
+
+    if (Abs(l2 - l1) > EPS_PARAM)
+    {
+      iLSubEnd = FillIntervalBounds(l1, l2, LKnots, NumSubs, anInertiaL, L1, L2, ErrL, ErrUL);
+      LMaxSubs = BRepGProp_Gauss::MaxSubs(iLSubEnd);
+
+      if (LMaxSubs > SM)
+        LMaxSubs = SM;
+
+      BRepGProp_Gauss::InitMass(0.0, 1, LMaxSubs, anInertiaL);
+      BRepGProp_Gauss::Init(ErrL,  0.0, 1, LMaxSubs);
+      BRepGProp_Gauss::Init(ErrUL, 0.0, 1, LMaxSubs);
+
+      do // while: L
+      {
+        if (++JL > iLSubEnd)
+        {
+          LRange[0] = IL = ErrL->Max();
+          LRange[1] = JL;
+          L1->Value(JL) = (L1->Value(IL) + L2->Value(IL)) * 0.5;
+          L2->Value(JL) = L2->Value(IL);
+          L2->Value(IL) = L1->Value(JL);
+        }
+        else
+          LRange[0] = IL = JL;
+
+        if (JL == LMaxSubs || Abs(L2->Value(JL) - L1->Value(JL)) < EPS_PARAM)
+        {
+          if (kLEnd == 1)
+          {
+            anInertiaL->ChangeValue(JL).Reset();
+            ErrL->Value(JL) = 0.0;
+          }
+          else
+          {
+            --JL;
+            EpsL = ErrorL;
+            Eps = EpsL / 0.9;
+            break;
+          }
+        }
+        else
+          for (kL = 0; kL < kLEnd; kL++)
+          {
+            iLS = LRange[kL];
+            lm = 0.5 * (L2->Value(iLS) + L1->Value(iLS));
+            lr = 0.5 * (L2->Value(iLS) - L1->Value(iLS));
+
+            CIx = CIy = CIz = CIxy = CIxz = CIyz = 0.0;
+
+            for (iGL = 0; iGL < iGLEnd; ++iGL)
+            {
+              CDim[iGL] = CIxx[iGL] = CIyy[iGL] = CIzz[iGL] = 0.0;
+
+              for (iL = 1; iL <= NbLGaussP[iGL]; iL++)
+              {
+                l = lm + lr * LGaussP[iGL]->Value(iL);
+                if (isNaturalRestriction)
+                { 
+                  v = l;
+                  u2 = BU2;
+                  Dul = LGaussW[iGL]->Value(iL);
+                }
+                else
+                {
+                  theSurface.D12d (l, Puv, Vuv);
+                  Dul = Vuv.Y() * LGaussW[iGL]->Value(iL);  // Dul = Du / Dl
+
+                  if (Abs(Dul) < EPS_PARAM)
+                    continue;
+
+                  v  = Puv.Y();
+                  u2 = Puv.X();
+
+                  // Check on cause out off bounds of value current parameter
+                  if (v < BV1)
+                    v = BV1;
+                  else if (v > BV2)
+                    v = BV2;
+
+                  if (u2 < BU1)
+                    u2 = BU1;
+                  else if (u2 > BU2)
+                    u2 = BU2; 
+                }
+
+                ErrUL->Value(iLS) = 0.0;
+                kUEnd = 1;
+                JU = 0;
+
+                if (Abs(u2 - u1) < EPS_PARAM)
+                  continue;
+
+                Handle_Vector aDummy;
+                iUSubEnd = FillIntervalBounds(u1, u2, UKnots, NumSubs, anInertiaU, U1, U2, ErrU, aDummy);
+                UMaxSubs = BRepGProp_Gauss::MaxSubs(iUSubEnd);
+
+                if (UMaxSubs > SM)
+                  UMaxSubs = SM;
+
+                BRepGProp_Gauss::InitMass(0.0, 1, UMaxSubs, anInertiaU);
+                BRepGProp_Gauss::Init(ErrU, 0.0, 1, UMaxSubs);
+                ErrorU = 0.0;
+
+                do
+                {//while: U
+                  if (++JU > iUSubEnd)
+                  {
+                    URange[0] = IU = ErrU->Max();
+                    URange[1] = JU;
+
+                    U1->Value(JU) = (U1->Value(IU) + U2->Value(IU)) * 0.5;
+                    U2->Value(JU) = U2->Value(IU);
+                    U2->Value(IU) = U1->Value(JU);
+                  }
+                  else
+                    URange[0] = IU = JU;
+
+                  if (JU == UMaxSubs || Abs(U2->Value(JU) - U1->Value(JU)) < EPS_PARAM)
+                    if (kUEnd == 1)
+                    {
+                      ErrU->Value(JU) = 0.0;
+                      anInertiaU->ChangeValue(JU).Reset();
+                    }
+                    else
+                    {
+                      --JU;
+                      EpsU = ErrorU;
+                      Eps  = 10. * EpsU * Abs((u2 - u1) * Dul);
+                      EpsL = 0.9 * Eps;
+                      break;
+                    }
+                  else
+                  {
+                    gp_Pnt aPoint;
+                    gp_Vec aNormal;
+
+                    for (kU = 0; kU < kUEnd; ++kU)
+                    {
+                      BRepGProp_Gauss::Inertia aLocal[2];
+
+                      Standard_Integer iUS = URange[kU];
+                      const Standard_Integer aLength = iGLEnd - iGL;
+
+                      const Standard_Real um = 0.5 * (U2->Value(iUS) + U1->Value(iUS));
+                      const Standard_Real ur = 0.5 * (U2->Value(iUS) - U1->Value(iUS));
+
+                      for (Standard_Integer iGU = 0; iGU < aLength; ++iGU)
+                      {
+                        for (Standard_Integer iU = 1; iU <= NbUGaussP[iGU]; ++iU)
+                        {
+                          Standard_Real w = UGaussW[iGU]->Value(iU);
+                          const Standard_Real u = um + ur * UGaussP[iGU]->Value(iU);
+
+                          theSurface.Normal(u, v, aPoint, aNormal);
+
+                          if (myType == Vinert)
+                          {
+                            computeVInertiaOfElementaryPart(
+                              aPoint, aNormal, theLocation, w, theCoeff, theIsByPoint, aLocal[iGU]);
+                          }
+                          else
+                          {
+                            if (iGU > 0)
+                              aLocal[iGU].Mass += (w * aNormal.Magnitude());
+                            else
+                            {
+                              computeSInertiaOfElementaryPart(
+                                aPoint, aNormal, theLocation, w, aLocal[iGU]);
+                            }
+                          }
+                        }
+                      }
+
+                      BRepGProp_Gauss::Inertia& anUI =
+                        anInertiaU->ChangeValue(iUS);
+
+                      anUI.Mass = mult(aLocal[0].Mass, ur);
+
+                      if (myType == Vinert)
+                      {
+                        anUI.Ixx = mult(aLocal[0].Ixx, ur);
+                        anUI.Iyy = mult(aLocal[0].Iyy, ur);
+                        anUI.Izz = mult(aLocal[0].Izz, ur);
+                      }
+
+                      if (iGL > 0)
+                        continue;
+
+                      Standard_Real aDMass = Abs(aLocal[1].Mass - aLocal[0].Mass);
+
+                      if (myType == Vinert)
+                      {
+                        aLocal[1].Ixx = Abs(aLocal[1].Ixx - aLocal[0].Ixx);
+                        aLocal[1].Iyy = Abs(aLocal[1].Iyy - aLocal[0].Iyy);
+                        aLocal[1].Izz = Abs(aLocal[1].Izz - aLocal[0].Izz);
+
+                        anUI.Ix = mult(aLocal[0].Ix, ur);
+                        anUI.Iy = mult(aLocal[0].Iy, ur);
+                        anUI.Iz = mult(aLocal[0].Iz, ur);
+
+                        anUI.Ixy = mult(aLocal[0].Ixy, ur);
+                        anUI.Ixz = mult(aLocal[0].Ixz, ur);
+                        anUI.Iyz = mult(aLocal[0].Iyz, ur);
+
+                        #ifndef IS_MIN_DIM
+                        aDMass = aLocal[1].Ixx + aLocal[1].Iyy + aLocal[1].Izz;
+                        #endif
+
+                        ErrU->Value(iUS) = mult(aDMass, ur);
+                      }
+                      else
+                      {
+                        anUI.Ix  = mult(aLocal[0].Ix,  ur);
+                        anUI.Iy  = mult(aLocal[0].Iy,  ur);
+                        anUI.Iz  = mult(aLocal[0].Iz,  ur);
+                        anUI.Ixx = mult(aLocal[0].Ixx, ur);
+                        anUI.Iyy = mult(aLocal[0].Iyy, ur);
+                        anUI.Izz = mult(aLocal[0].Izz, ur);
+                        anUI.Ixy = mult(aLocal[0].Ixy, ur);
+                        anUI.Ixz = mult(aLocal[0].Ixz, ur);
+                        anUI.Iyz = mult(aLocal[0].Iyz, ur);
+
+                        ErrU->Value(iUS) = mult(aDMass, ur);
+                      }
+                    }
+                  }
+
+                  if (JU == iUSubEnd)
+                  {
+                    kUEnd = 2;
+                    ErrorU = ErrU->Value(ErrU->Max());
+                  }
+                } while ( (ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1 );
+
+                for (i = 1; i <= JU; ++i)
+                {
+                  const BRepGProp_Gauss::Inertia& anIU =
+                    anInertiaU->Value(i);
+
+                  CDim[iGL] = add(CDim[iGL], mult(anIU.Mass, Dul));
+                  CIxx[iGL] = add(CIxx[iGL], mult(anIU.Ixx,  Dul));
+                  CIyy[iGL] = add(CIyy[iGL], mult(anIU.Iyy,  Dul));
+                  CIzz[iGL] = add(CIzz[iGL], mult(anIU.Izz,  Dul));
+                }
+
+                if (iGL > 0)
+                  continue;
+
+                ErrUL->Value(iLS) = ErrorU * Abs((u2 - u1) * Dul);
+
+                for (i = 1; i <= JU; ++i)
+                {
+                  const BRepGProp_Gauss::Inertia& anIU =
+                    anInertiaU->Value(i);
+
+                  CIx = add(CIx, mult(anIU.Ix, Dul));
+                  CIy = add(CIy, mult(anIU.Iy, Dul));
+                  CIz = add(CIz, mult(anIU.Iz, Dul));
+
+                  CIxy = add(CIxy, mult(anIU.Ixy, Dul));
+                  CIxz = add(CIxz, mult(anIU.Ixz, Dul));
+                  CIyz = add(CIyz, mult(anIU.Iyz, Dul));
+                }
+              }//for: iL 
+            }//for: iGL
+
+            BRepGProp_Gauss::Inertia& aLI = anInertiaL->ChangeValue(iLS);
+
+            aLI.Mass = mult(CDim[0], lr);
+            aLI.Ixx  = mult(CIxx[0], lr);
+            aLI.Iyy  = mult(CIyy[0], lr);
+            aLI.Izz  = mult(CIzz[0], lr);
+
+            if (iGLEnd == 2)
+            {
+              Standard_Real aSubDim = Abs(CDim[1] - CDim[0]);
+
+              if (myType == Vinert)
+              {
+                ErrorU = ErrUL->Value(iLS);
+
+                CIxx[1] = Abs(CIxx[1] - CIxx[0]);
+                CIyy[1] = Abs(CIyy[1] - CIyy[0]);
+                CIzz[1] = Abs(CIzz[1] - CIzz[0]);
+
+                #ifndef IS_MIN_DIM
+                aSubDim = CIxx[1] + CIyy[1] + CIzz[1];
+                #endif
+
+                ErrL->Value(iLS) = add(mult(aSubDim, lr), ErrorU);
+              }
+              else
+              {
+                ErrL->Value(iLS) = add(mult(aSubDim, lr), ErrUL->Value(iLS));
+              }
+            }
+
+            aLI.Ix = mult(CIx, lr);
+            aLI.Iy = mult(CIy, lr);
+            aLI.Iz = mult(CIz, lr);
+
+            aLI.Ixy = mult(CIxy, lr);
+            aLI.Ixz = mult(CIxz, lr);
+            aLI.Iyz = mult(CIyz, lr);
+          }//for: (kL)iLS
+
+          // Calculate/correct epsilon of computation by current value of dim
+          // That is need for not spend time for 
+          if (JL == iLSubEnd)
+          {
+            kLEnd = 2;
+
+            Standard_Real DDim = 0.0;
+            for (i = 1; i <= JL; ++i)
+              DDim += anInertiaL->Value(i).Mass;
+
+            #ifndef IS_MIN_DIM
+            {
+              if (myType == Vinert)
+              {
+                Standard_Real DIxx = 0.0, DIyy = 0.0, DIzz = 0.0;
+                for (i = 1; i <= JL; ++i)
+                {
+                  const BRepGProp_Gauss::Inertia& aLocalL =
+                    anInertiaL->Value(i);
+
+                  DIxx += aLocalL.Ixx;
+                  DIyy += aLocalL.Iyy;
+                  DIzz += aLocalL.Izz;
+                }
+
+                DDim = Abs(DIxx) + Abs(DIyy) + Abs(DIzz);
+              }
+            }
+            #endif
+
+            DDim = Abs(DDim * anEpsilon);
+
+            if (DDim > Eps)
+            {
+              Eps  = DDim;
+              EpsL = 0.9 * Eps;
+            }
+          }
+          if (kLEnd == 2)
+          {
+            ErrorL = ErrL->Value(ErrL->Max());
+          }
+      } while ( (ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1 );
+
+      for ( i = 1; i <= JL; i++ )
+        addAndRestoreInertia(anInertiaL->Value(i), anInertia);
+
+      ErrorLMax = Max(ErrorLMax, ErrorL);
+    }
+
+    if (isNaturalRestriction)
+      break;
+
+    theDomain.Next();
+  }
+
+  if (myType == Vinert)
+    convert(anInertia, theCoeff, theIsByPoint, theOutGravityCenter, theOutInertia, theOutMass);
+  else
+    convert(anInertia, theOutGravityCenter, theOutInertia, theOutMass);
+
+  if (iGLEnd == 2)
+  {
+    if (theOutMass != 0.0)
+    {
+      Eps = ErrorLMax / Abs(theOutMass);
+
+      #ifndef IS_MIN_DIM
+      {
+        if (myType == Vinert)
+          Eps = ErrorLMax / (Abs(anInertia.Ixx) +
+                             Abs(anInertia.Iyy) +
+                             Abs(anInertia.Izz));
+      }
+      #endif
+    }
+    else
+    {
+      Eps = 0.0;
+    }
+  }
+  else
+  {
+    Eps = anEpsilon;
+  }
+
+  return Eps;
+}
+
+//=======================================================================
+//function : Compute
+//purpose  : 
+//=======================================================================
+Standard_Real BRepGProp_Gauss::Compute(BRepGProp_Face&     theSurface,
+                                       BRepGProp_Domain&   theDomain,
+                                       const gp_Pnt&       theLocation,
+                                       const Standard_Real theEps,
+                                       Standard_Real&      theOutMass,
+                                       gp_Pnt&             theOutGravityCenter,
+                                       gp_Mat&             theOutInertia)
+{
+  Standard_ASSERT_RAISE(myType == Sinert, "BRepGProp_Gauss: Incorrect type");
+
+  return Compute(theSurface,
+                 theDomain,
+                 theLocation,
+                 theEps,
+                 NULL,
+                 Standard_True,
+                 theOutMass,
+                 theOutGravityCenter,
+                 theOutInertia);
+}
+
+//=======================================================================
+//function : Compute
+//purpose  : 
+//=======================================================================
+void BRepGProp_Gauss::Compute(BRepGProp_Face&   theSurface,
+                              BRepGProp_Domain& theDomain,
+                              const gp_Pnt&     theLocation,
+                              Standard_Real&    theOutMass,
+                              gp_Pnt&           theOutGravityCenter,
+                              gp_Mat&           theOutInertia)
+{
+  Standard_ASSERT_RAISE(myType == Sinert, "BRepGProp_Gauss: Incorrect type");
+
+  Standard_Real u1, u2, v1, v2;
+  theSurface.Bounds (u1, u2, v1, v2);
+  checkBounds(u1, u2, v1, v2);
+
+  const Standard_Integer NbUGaussgp_Pnts =
+    Min(theSurface.UIntegrationOrder(), math::GaussPointsMax());
+
+  const Standard_Integer NbVGaussgp_Pnts =
+    Min(theSurface.VIntegrationOrder(), math::GaussPointsMax());
+
+  const Standard_Integer NbGaussgp_Pnts =
+    Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts);
+
+  // Number of Gauss points for the integration on the face
+  math_Vector GaussSPV (1, NbGaussgp_Pnts);
+  math_Vector GaussSWV (1, NbGaussgp_Pnts);
+  math::GaussPoints (NbGaussgp_Pnts, GaussSPV);
+  math::GaussWeights(NbGaussgp_Pnts, GaussSWV);
+
+  BRepGProp_Gauss::Inertia anInertia;
+  while (theDomain.More())
+  {
+    theSurface.Load(theDomain.Value());
+
+    const Standard_Integer NbCGaussgp_Pnts =
+      Min(theSurface.IntegrationOrder(), math::GaussPointsMax());
+
+    math_Vector GaussCP(1, NbCGaussgp_Pnts);
+    math_Vector GaussCW(1, NbCGaussgp_Pnts);
+    math::GaussPoints (NbCGaussgp_Pnts, GaussCP);
+    math::GaussWeights(NbCGaussgp_Pnts, GaussCW);
+
+    
+    const Standard_Real l1 = theSurface.FirstParameter();
+    const Standard_Real l2 = theSurface.LastParameter ();
+    const Standard_Real lm = 0.5 * (l2 + l1);
+    const Standard_Real lr = 0.5 * (l2 - l1);
+
+    BRepGProp_Gauss::Inertia aCInertia;
+    for (Standard_Integer i = 1; i <= NbCGaussgp_Pnts; ++i)
+    {
+      const Standard_Real l = lm + lr * GaussCP(i);
+
+      gp_Pnt2d Puv;
+      gp_Vec2d Vuv;
+      theSurface.D12d(l, Puv, Vuv);
+
+      const Standard_Real v = Puv.Y();
+      u2  = Puv.X();
+
+      const Standard_Real Dul = Vuv.Y() * GaussCW(i);
+      const Standard_Real um  = 0.5 * (u2 + u1);
+      const Standard_Real ur  = 0.5 * (u2 - u1);
+
+      BRepGProp_Gauss::Inertia aLocalInertia;
+      for (Standard_Integer j = 1; j <= NbGaussgp_Pnts; ++j)
+      {
+        const Standard_Real u = add(um, mult(ur, GaussSPV(j)));
+        const Standard_Real aWeight = Dul * GaussSWV(j);
+
+        gp_Pnt aPoint;
+        gp_Vec aNormal;
+        theSurface.Normal (u, v, aPoint, aNormal);
+
+        computeSInertiaOfElementaryPart(aPoint, aNormal, theLocation, aWeight, aLocalInertia);
+      }
+
+      multAndRestoreInertia(ur, aLocalInertia);
+      addAndRestoreInertia (aLocalInertia, aCInertia);
+    }
+
+    multAndRestoreInertia(lr, aCInertia);
+    addAndRestoreInertia (aCInertia, anInertia);
+
+    theDomain.Next();
+  }
+
+  convert(anInertia, theOutGravityCenter, theOutInertia, theOutMass);
+}
+
+//=======================================================================
+//function : Compute
+//purpose  : 
+//=======================================================================
+void BRepGProp_Gauss::Compute(BRepGProp_Face&         theSurface,
+                              BRepGProp_Domain&       theDomain,
+                              const gp_Pnt&           theLocation,
+                              const Standard_Real     theCoeff[],
+                              const Standard_Boolean  theIsByPoint,
+                              Standard_Real&          theOutMass,
+                              gp_Pnt&                 theOutGravityCenter,
+                              gp_Mat&                 theOutInertia)
+{
+  Standard_ASSERT_RAISE(myType == Vinert, "BRepGProp_Gauss: Incorrect type");
+
+  Standard_Real u1, v1, u2, v2;
+  theSurface.Bounds (u1, u2, v1, v2);
+  checkBounds(u1, u2, v1, v2);
+
+  Standard_Real _u2 = u2;  //OCC104
+
+  BRepGProp_Gauss::Inertia anInertia;
+  while (theDomain.More())
+  {
+    theSurface.Load(theDomain.Value());
+
+    const Standard_Integer aVNbCGaussgp_Pnts =
+      theSurface.VIntegrationOrder();
+
+    const Standard_Integer aNbGaussgp_Pnts =
+      Min( Max(theSurface.IntegrationOrder(), aVNbCGaussgp_Pnts), math::GaussPointsMax() );
+
+    math_Vector GaussP(1, aNbGaussgp_Pnts);
+    math_Vector GaussW(1, aNbGaussgp_Pnts);
+    math::GaussPoints (aNbGaussgp_Pnts, GaussP);
+    math::GaussWeights(aNbGaussgp_Pnts, GaussW);
+
+    const Standard_Real l1 = theSurface.FirstParameter();
+    const Standard_Real l2 = theSurface.LastParameter();
+    const Standard_Real lm = 0.5 * (l2 + l1);
+    const Standard_Real lr = 0.5 * (l2 - l1);
+
+    BRepGProp_Gauss::Inertia aCInertia;
+    for (Standard_Integer i = 1; i <= aNbGaussgp_Pnts; ++i)
+    {
+      const Standard_Real l = lm + lr * GaussP(i);
+
+      gp_Pnt2d Puv;
+      gp_Vec2d Vuv;
+
+      theSurface.D12d(l, Puv, Vuv);
+
+      u2  = Puv.X();
+      u2 = Min( Max(u1, u2), _u2 ); // OCC104
+      const Standard_Real v = Min(Max(Puv.Y(), v1), v2);
+
+      const Standard_Real Dul = Vuv.Y() * GaussW(i);
+      const Standard_Real um  = 0.5 * (u2 + u1);
+      const Standard_Real ur  = 0.5 * (u2 - u1);
+
+      BRepGProp_Gauss::Inertia aLocalInertia;
+      for (Standard_Integer j = 1; j <= aNbGaussgp_Pnts; ++j)
+      {
+        const Standard_Real u = um + ur * GaussP(j);
+        const Standard_Real aWeight = Dul * GaussW(j);
+
+        gp_Pnt aPoint;
+        gp_Vec aNormal;
+
+        theSurface.Normal(u, v, aPoint, aNormal);
+
+        computeVInertiaOfElementaryPart(
+          aPoint,
+          aNormal,
+          theLocation,
+          aWeight,
+          theCoeff,
+          theIsByPoint,
+          aLocalInertia);
+      }
+
+      multAndRestoreInertia(ur, aLocalInertia);
+      addAndRestoreInertia (aLocalInertia, aCInertia);
+    }
+
+    multAndRestoreInertia(lr,        aCInertia);
+    addAndRestoreInertia (aCInertia, anInertia);
+
+    theDomain.Next();
+  }
+
+  convert(anInertia, theCoeff, theIsByPoint, theOutGravityCenter, theOutInertia, theOutMass);
+}
+
+//=======================================================================
+//function : Compute
+//purpose  : 
+//=======================================================================
+void BRepGProp_Gauss::Compute(const BRepGProp_Face&  theSurface,
+                              const gp_Pnt&          theLocation,
+                              const Standard_Real    theCoeff[],
+                              const Standard_Boolean theIsByPoint,
+                              Standard_Real&         theOutMass,
+                              gp_Pnt&                theOutGravityCenter,
+                              gp_Mat&                theOutInertia)
+{
+  Standard_Real LowerU, UpperU, LowerV, UpperV;
+  theSurface.Bounds(LowerU, UpperU, LowerV, UpperV);
+  checkBounds(LowerU, UpperU, LowerV, UpperV);
+
+  const Standard_Integer UOrder =
+    Min(theSurface.UIntegrationOrder(), math::GaussPointsMax());
+  const Standard_Integer VOrder =
+    Min(theSurface.VIntegrationOrder(), math::GaussPointsMax());
+
+  // Gauss points and weights
+  math_Vector GaussPU(1, UOrder);
+  math_Vector GaussWU(1, UOrder);
+  math_Vector GaussPV(1, VOrder);
+  math_Vector GaussWV(1, VOrder);
+
+  math::GaussPoints (UOrder, GaussPU);
+  math::GaussWeights(UOrder, GaussWU);
+  math::GaussPoints (VOrder, GaussPV);
+  math::GaussWeights(VOrder, GaussWV);
+
+  const Standard_Real um = 0.5 * add(UpperU,  LowerU);
+  const Standard_Real vm = 0.5 * add(UpperV,  LowerV);
+  Standard_Real ur = 0.5 * add(UpperU, -LowerU);
+  Standard_Real vr = 0.5 * add(UpperV, -LowerV);
+
+  gp_Pnt aPoint;
+  gp_Vec aNormal;
+
+  BRepGProp_Gauss::Inertia anInertia;
+  for (Standard_Integer j = 1; j <= VOrder; ++j)
+  {
+    BRepGProp_Gauss::Inertia anInertiaOfElementaryPart;
+    const Standard_Real v = add(vm, mult(vr, GaussPV(j)));
+
+    for (Standard_Integer i = 1; i <= UOrder; ++i)
+    {
+      const Standard_Real aWeight = GaussWU(i);
+      const Standard_Real u = add(um, mult(ur, GaussPU (i)));
+      theSurface.Normal(u, v, aPoint, aNormal);
+
+      if (myType == Vinert)
+      {
+        computeVInertiaOfElementaryPart(
+          aPoint,
+          aNormal,
+          theLocation,
+          aWeight,
+          theCoeff,
+          theIsByPoint,
+          anInertiaOfElementaryPart);
+      }
+      else // Sinert
+      {
+        computeSInertiaOfElementaryPart(
+          aPoint,
+          aNormal,
+          theLocation,
+          aWeight,
+          anInertiaOfElementaryPart);
+      }
+    }
+
+    multAndRestoreInertia(GaussWV(j), anInertiaOfElementaryPart);
+    addAndRestoreInertia (anInertiaOfElementaryPart,  anInertia);
+  }
+  vr            = mult(vr, ur);
+  anInertia.Ixx = mult(vr, anInertia.Ixx);
+  anInertia.Iyy = mult(vr, anInertia.Iyy);
+  anInertia.Izz = mult(vr, anInertia.Izz);
+  anInertia.Ixy = mult(vr, anInertia.Ixy);
+  anInertia.Ixz = mult(vr, anInertia.Ixz);
+  anInertia.Iyz = mult(vr, anInertia.Iyz);
+
+  if (myType == Vinert)
+  {
+    convert(anInertia, theCoeff, theIsByPoint, theOutGravityCenter, theOutInertia, theOutMass);
+  }
+  else // Sinert
+  {
+    convert(anInertia, theOutGravityCenter, theOutInertia, theOutMass);
+  }
+
+  theOutMass *= vr;
+}
+
+//=======================================================================
+//function : Compute
+//purpose  : 
+//=======================================================================
+void BRepGProp_Gauss::Compute(const BRepGProp_Face& theSurface,
+                              const gp_Pnt&         theLocation,
+                              Standard_Real&        theOutMass,
+                              gp_Pnt&               theOutGravityCenter,
+                              gp_Mat&               theOutInertia)
+{
+  Standard_ASSERT_RAISE(myType == Sinert, "BRepGProp_Gauss: Incorrect type");
+
+  Compute(theSurface,
+          theLocation,
+          NULL,
+          Standard_True,
+          theOutMass,
+          theOutGravityCenter,
+          theOutInertia);
+}
diff --git a/src/BRepGProp/BRepGProp_Gauss.hxx b/src/BRepGProp/BRepGProp_Gauss.hxx
new file mode 100644 (file)
index 0000000..6ded4b8
--- /dev/null
@@ -0,0 +1,284 @@
+// Copyright (c) 2015 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _BRepGProp_Gauss_HeaderFile
+#define _BRepGProp_Gauss_HeaderFile
+
+#include <NCollection_Handle.hxx>
+#include <NCollection_Array1.hxx>
+
+class math_Vector;
+
+//! Class performs computing of the global inertia properties
+//! of geometric object in 3D space by adaptive and non-adaptive
+//! 2D Gauss integration algorithms.
+class BRepGProp_Gauss
+{
+  //! Auxiliary structure for storing of inertial moments.
+  struct Inertia
+  {
+    //! Mass of the current system (without density).
+    //! May correspond to: length, area, volume.
+    Standard_Real Mass;
+
+    //! Static moments of inertia.
+    Standard_Real Ix;
+    Standard_Real Iy;
+    Standard_Real Iz;
+
+    //! Quadratic moments of inertia.
+    Standard_Real Ixx;
+    Standard_Real Iyy;
+    Standard_Real Izz;
+    Standard_Real Ixy;
+    Standard_Real Ixz;
+    Standard_Real Iyz;
+
+    //! Default constructor.
+    Inertia();
+
+    //! Zeroes all values.
+    void Reset();
+  };
+
+  typedef NCollection_Handle< NCollection_Array1<Inertia> > InertiaArray;
+  typedef NCollection_Handle<math_Vector>                   Handle_Vector;
+  typedef Standard_Real(*BRepGProp_GaussFunc)(const Standard_Real, const Standard_Real);
+
+public: //! @name public API
+
+  //! Describes types of geometric objects.
+  //! - Vinert is 3D closed region of space delimited with:
+  //! -- Surface;
+  //! -- Point and Surface;
+  //! -- Plane and Surface.
+  //! - Sinert is face in 3D space.
+  typedef enum { Vinert = 0, Sinert } BRepGProp_GaussType;
+
+  //! Constructor
+  Standard_EXPORT explicit BRepGProp_Gauss(const BRepGProp_GaussType theType);
+
+  //! Computes the global properties of a solid region of 3D space which can be
+  //! delimited by the surface and point or surface and plane. Surface can be closed.
+  //! The method is quick and its precision is enough for many cases of analytical surfaces.
+  //! Non-adaptive 2D Gauss integration with predefined numbers of Gauss points
+  //! is used. Numbers of points depend on types of surfaces and curves.
+  //! Error of the computation is not calculated.
+  //! @param theSurface - bounding surface of the region;
+  //! @param theLocation - location of the point or the plane;
+  //! @param theCoeff - plane coefficients;
+  //! @param theIsByPoint - flag of restricition (point/plane);
+  //! @param theOutMass[out] - mass (volume) of region;
+  //! @param theOutGravityCenter[out] - garvity center of region;
+  //! @param theOutInertia[out] - matrix of inertia;
+  Standard_EXPORT void Compute(
+    const BRepGProp_Face&  theSurface,
+    const gp_Pnt&          theLocation,
+    const Standard_Real    theCoeff[],
+    const Standard_Boolean theIsByPoint,
+    Standard_Real&         theOutMass,
+    gp_Pnt&                theOutGravityCenter,
+    gp_Mat&                theOutInertia);
+
+  //! Computes the global properties of a surface. Surface can be closed.
+  //! The method is quick and its precision is enough for many cases of analytical surfaces.
+  //! Non-adaptive 2D Gauss integration with predefined numbers of Gauss points
+  //! is used. Numbers of points depend on types of surfaces and curves.
+  //! Error of the computation is not calculated.
+  //! @param theSurface - bounding surface of the region;
+  //! @param theLocation - surface location;
+  //! @param theOutMass[out] - mass (volume) of region;
+  //! @param theOutGravityCenter[out] - garvity center of region;
+  //! @param theOutInertia[out] - matrix of inertia;
+  Standard_EXPORT void Compute(
+    const BRepGProp_Face&  theSurface,
+    const gp_Pnt&          theLocation,
+    Standard_Real&         theOutMass,
+    gp_Pnt&                theOutGravityCenter,
+    gp_Mat&                theOutInertia);
+
+  //! Computes the global properties of a region of 3D space which can be
+  //! delimited by the surface and point or surface and plane. Surface can be closed.
+  //! The method is quick and its precision is enough for many cases of analytical surfaces.
+  //! Non-adaptive 2D Gauss integration with predefined numbers of Gauss points is used.
+  //! Numbers of points depend on types of surfaces and curves.
+  //! Error of the computation is not calculated.
+  //! @param theSurface - bounding surface of the region;
+  //! @param theDomain - surface boundings;
+  //! @param theLocation - location of the point or the plane;
+  //! @param theCoeff - plane coefficients;
+  //! @param theIsByPoint - flag of restricition (point/plane);
+  //! @param theOutMass[out] - mass (volume) of region;
+  //! @param theOutGravityCenter[out] - garvity center of region;
+  //! @param theOutInertia[out] - matrix of inertia;
+  Standard_EXPORT void Compute(
+    BRepGProp_Face&        theSurface,
+    BRepGProp_Domain&      theDomain,
+    const gp_Pnt&          theLocation,
+    const Standard_Real    theCoeff[],
+    const Standard_Boolean theIsByPoint,
+    Standard_Real&         theOutMass,
+    gp_Pnt&                theOutGravityCenter,
+    gp_Mat&                theOutInertia);
+
+  //! Computes the global properties of a surface. Surface can be closed.
+  //! The method is quick and its precision is enough for many cases of analytical surfaces.
+  //! Non-adaptive 2D Gauss integration with predefined numbers of Gauss points
+  //! is used. Numbers of points depend on types of surfaces and curves.
+  //! Error of the computation is not calculated.
+  //! @param theSurface - bounding surface of the region;
+  //! @param theDomain - surface boundings;
+  //! @param theLocation - surface location;
+  //! @param theOutMass[out] - mass (volume) of region;
+  //! @param theOutGravityCenter[out] - garvity center of region;
+  //! @param theOutInertia[out] - matrix of inertia;
+  Standard_EXPORT void Compute(
+    BRepGProp_Face&        theSurface,
+    BRepGProp_Domain&      theDomain,
+    const gp_Pnt&          theLocation,
+    Standard_Real&         theOutMass,
+    gp_Pnt&                theOutGravityCenter,
+    gp_Mat&                theOutInertia);
+
+  //! Computes the global properties of the region of 3D space which can be
+  //! delimited by the surface and point or surface and plane.
+  //! Adaptive 2D Gauss integration is used.
+  //! If Epsilon more than 0.001 then algorithm performs non-adaptive integration.
+  //! @param theSurface - bounding surface of the region;
+  //! @param theDomain - surface boundings;
+  //! @param theLocation - location of the point or the plane;
+  //! @param theEps - maximal relative error of computed mass (volume) for face;
+  //! @param theCoeff - plane coefficients;
+  //! @param theIsByPoint - flag of restricition (point/plane);
+  //! @param theOutMass[out] - mass (volume) of region;
+  //! @param theOutGravityCenter[out] - garvity center of region;
+  //! @param theOutInertia[out] - matrix of inertia;
+  //! @return value of error which is calculated as
+  //! Abs((M(i+1)-M(i))/M(i+1)), M(i+1) and M(i) are values
+  //! for two successive steps of adaptive integration.
+  Standard_EXPORT Standard_Real Compute(
+    BRepGProp_Face&        theSurface,
+    BRepGProp_Domain&      theDomain,
+    const gp_Pnt&          theLocation,
+    const Standard_Real    theEps,
+    const Standard_Real    theCoeff[],
+    const Standard_Boolean theByPoint,
+    Standard_Real&         theOutMass,
+    gp_Pnt&                theOutGravityCenter,
+    gp_Mat&                theOutInertia);
+
+  //! Computes the global properties of the face. Adaptive 2D Gauss integration is used.
+  //! If Epsilon more than 0.001 then algorithm performs non-adaptive integration.
+  //! @param theSurface - bounding surface of the region;
+  //! @param theDomain - surface boundings;
+  //! @param theLocation - surface location;
+  //! @param theEps - maximal relative error of computed mass (square) for face;
+  //! @param theOutMass[out] - mass (volume) of region;
+  //! @param theOutGravityCenter[out] - garvity center of region;
+  //! @param theOutInertia[out] - matrix of inertia;
+  //! @return value of error which is calculated as
+  //! Abs((M(i+1)-M(i))/M(i+1)), M(i+1) and M(i) are values
+  //! for two successive steps of adaptive integration.
+  Standard_EXPORT Standard_Real Compute(
+    BRepGProp_Face&        theSurface,
+    BRepGProp_Domain&      theDomain,
+    const gp_Pnt&          theLocation,
+    const Standard_Real    theEps,
+    Standard_Real&         theOutMass,
+    gp_Pnt&                theOutGravityCenter,
+    gp_Mat&                theOutInertia);
+
+private: //! @name private methods
+
+  BRepGProp_Gauss(BRepGProp_Gauss const&);
+  BRepGProp_Gauss& operator= (BRepGProp_Gauss const&);
+
+  void computeVInertiaOfElementaryPart(
+    const gp_Pnt&             thePoint,
+    const gp_Vec&             theNormal,
+    const gp_Pnt&             theLocation,
+    const Standard_Real       theWeight,
+    const Standard_Real       theCoeff[],
+    const Standard_Boolean    theIsByPoint,
+    BRepGProp_Gauss::Inertia& theOutInertia);
+
+  void computeSInertiaOfElementaryPart(
+    const gp_Pnt&             thePoint,
+    const gp_Vec&             theNormal,
+    const gp_Pnt&             theLocation,
+    const Standard_Real       theWeight,
+    BRepGProp_Gauss::Inertia& theOutInertia);
+
+  void checkBounds(
+    const Standard_Real theU1,
+    const Standard_Real theU2,
+    const Standard_Real theV1,
+    const Standard_Real theV2);
+
+  void addAndRestoreInertia(
+    const BRepGProp_Gauss::Inertia& theInInertia,
+    BRepGProp_Gauss::Inertia&       theOutInertia);
+
+  void multAndRestoreInertia(
+    const Standard_Real       theValue,
+    BRepGProp_Gauss::Inertia& theInertia);
+
+  void convert(
+    const BRepGProp_Gauss::Inertia& theInertia,
+    gp_Pnt&                         theOutGravityCenter,
+    gp_Mat&                         theOutMatrixOfInertia,
+    Standard_Real&                  theOutMass);
+
+  void convert(
+    const BRepGProp_Gauss::Inertia& theInertia,
+    const Standard_Real             theCoeff[],
+    const Standard_Boolean          theIsByPoint,
+    gp_Pnt&                         theOutGravityCenter,
+    gp_Mat&                         theOutMatrixOfInertia,
+    Standard_Real&                  theOutMass);
+
+  static Standard_Integer MaxSubs(
+    const Standard_Integer theN,
+    const Standard_Integer theCoeff = 32);
+
+  static void Init(
+    Handle_Vector&         theOutVec,
+    const Standard_Real    theValue,
+    const Standard_Integer theFirst = 0,
+    const Standard_Integer theLast  = 0);
+
+  static void InitMass(
+    const Standard_Real    theValue,
+    const Standard_Integer theFirst,
+    const Standard_Integer theLast,
+    InertiaArray&          theArray);
+
+  static Standard_Integer FillIntervalBounds(
+    const Standard_Real         theA,
+    const Standard_Real         theB,
+    const TColStd_Array1OfReal& theKnots,
+    const Standard_Integer      theNumSubs,
+    InertiaArray&               theInerts,
+    Handle_Vector&              theParam1,
+    Handle_Vector&              theParam2,
+    Handle_Vector&              theError,
+    Handle_Vector&              theCommonError);
+
+private: //! @name private fields
+
+  BRepGProp_GaussType myType; //!< Type of geometric object
+  BRepGProp_GaussFunc add;    //!< Pointer on the add function
+  BRepGProp_GaussFunc mult;   //!< Pointer on the mult function
+};
+
+#endif
\ No newline at end of file
index 24b8877..b1928c1 100644 (file)
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
-#include <BRepGProp_Sinert.ixx>
-
-#include <math.hxx>
-#include <TColStd_Array1OfReal.hxx>
-#include <Precision.hxx>
-
-class HMath_Vector{
-  math_Vector *pvec;
-  void operator=(const math_Vector&){}
- public:
-  HMath_Vector(){ pvec = 0;}
-  HMath_Vector(math_Vector* pv){ pvec = pv;}
-  ~HMath_Vector(){ if(pvec != 0) delete pvec;}
-  void operator=(math_Vector* pv){ if(pvec != pv && pvec != 0) delete pvec;  pvec = pv;}
-  Standard_Real& operator()(Standard_Integer i){ return (*pvec).operator()(i);}
-  const Standard_Real& operator()(Standard_Integer i) const{ return (*pvec).operator()(i);}
-  const math_Vector* operator->() const{ return pvec;}
-  math_Vector* operator->(){ return pvec;}
-  math_Vector* Vector(){ return pvec;}
-  math_Vector* Init(Standard_Real v, Standard_Integer i = 0, Standard_Integer iEnd = 0){ 
-    if(pvec == 0) return pvec;
-    if(iEnd - i == 0) pvec->Init(v);
-    else { Standard_Integer End = (iEnd <= pvec->Upper()) ? iEnd : pvec->Upper();
-           for(; i <= End; i++) pvec->operator()(i) = v; }
-    return pvec;
-  }
-};
-
-static Standard_Real EPS_PARAM          = 1.e-12;
-static Standard_Real EPS_DIM            = 1.e-20;
-static Standard_Real ERROR_ALGEBR_RATIO = 2.0/3.0;
-
-static Standard_Integer  GPM        = 61;
-static Standard_Integer  SUBS_POWER = 32;
-static Standard_Integer  SM         = 1953;
-
-static math_Vector LGaussP0(1,GPM);
-static math_Vector LGaussW0(1,GPM);
-static math_Vector LGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
-static math_Vector LGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
-
-static math_Vector* LGaussP[] = {&LGaussP0,&LGaussP1};
-static math_Vector* LGaussW[] = {&LGaussW0,&LGaussW1};
-
-static HMath_Vector L1    = new math_Vector(1,SM,0.0);
-static HMath_Vector L2    = new math_Vector(1,SM,0.0);
-static HMath_Vector DimL  = new math_Vector(1,SM,0.0);
-static HMath_Vector ErrL  = new math_Vector(1,SM,0.0);
-static HMath_Vector ErrUL = new math_Vector(1,SM,0.0);
-static HMath_Vector IxL   = new math_Vector(1,SM,0.0);
-static HMath_Vector IyL   = new math_Vector(1,SM,0.0);
-static HMath_Vector IzL   = new math_Vector(1,SM,0.0);
-static HMath_Vector IxxL  = new math_Vector(1,SM,0.0);
-static HMath_Vector IyyL  = new math_Vector(1,SM,0.0);
-static HMath_Vector IzzL  = new math_Vector(1,SM,0.0);
-static HMath_Vector IxyL  = new math_Vector(1,SM,0.0);
-static HMath_Vector IxzL  = new math_Vector(1,SM,0.0);
-static HMath_Vector IyzL  = new math_Vector(1,SM,0.0);
-
-static math_Vector UGaussP0(1,GPM);
-static math_Vector UGaussW0(1,GPM);
-static math_Vector UGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
-static math_Vector UGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
-
-static math_Vector* UGaussP[] = {&UGaussP0,&UGaussP1};
-static math_Vector* UGaussW[] = {&UGaussW0,&UGaussW1};
-
-static HMath_Vector U1   = new math_Vector(1,SM,0.0);
-static HMath_Vector U2   = new math_Vector(1,SM,0.0);
-static HMath_Vector DimU = new math_Vector(1,SM,0.0);
-static HMath_Vector ErrU = new math_Vector(1,SM,0.0);
-static HMath_Vector IxU  = new math_Vector(1,SM,0.0);
-static HMath_Vector IyU  = new math_Vector(1,SM,0.0);
-static HMath_Vector IzU  = new math_Vector(1,SM,0.0);
-static HMath_Vector IxxU = new math_Vector(1,SM,0.0);
-static HMath_Vector IyyU = new math_Vector(1,SM,0.0);
-static HMath_Vector IzzU = new math_Vector(1,SM,0.0);
-static HMath_Vector IxyU = new math_Vector(1,SM,0.0);
-static HMath_Vector IxzU = new math_Vector(1,SM,0.0);
-static HMath_Vector IyzU = new math_Vector(1,SM,0.0);
-
-static inline Standard_Real MultiplicationInf(Standard_Real theMA, Standard_Real theMB)
-{
-  if((theMA == 0.0) || (theMB == 0.0))  //strictly zerro (without any tolerances)
-    return 0.0;
-
-  if(Precision::IsPositiveInfinite(theMA))
-  {
-    if(theMB < 0.0)
-      return -Precision::Infinite();
-    else
-      return Precision::Infinite();
-  }
-
-  if(Precision::IsPositiveInfinite(theMB))
-  {
-    if(theMA < 0.0)
-      return -Precision::Infinite();
-    else
-      return Precision::Infinite();
-  }
-
-  if(Precision::IsNegativeInfinite(theMA))
-  {
-    if(theMB < 0.0)
-      return +Precision::Infinite();
-    else
-      return -Precision::Infinite();
-  }
-
-  if(Precision::IsNegativeInfinite(theMB))
-  {
-    if(theMA < 0.0)
-      return +Precision::Infinite();
-    else
-      return -Precision::Infinite();
-  }
-
-  return (theMA * theMB);
-}
-
-static inline Standard_Real AdditionInf(Standard_Real theMA, Standard_Real theMB)
+#include<BRepGProp_Sinert.ixx>
+#include<BRepGProp_Gauss.hxx>
+
+//=======================================================================
+//function : BRepGProp_Sinert
+//purpose  : Constructor
+//=======================================================================
+BRepGProp_Sinert::BRepGProp_Sinert()
 {
-  if(Precision::IsPositiveInfinite(theMA))
-  {
-    if(Precision::IsNegativeInfinite(theMB))
-      return 0.0;
-    else
-      return Precision::Infinite();
-  }
-
-  if(Precision::IsPositiveInfinite(theMB))
-  {
-    if(Precision::IsNegativeInfinite(theMA))
-      return 0.0;
-    else
-      return Precision::Infinite();
-  }
-
-  if(Precision::IsNegativeInfinite(theMA))
-  {
-    if(Precision::IsPositiveInfinite(theMB))
-      return 0.0;
-    else
-      return -Precision::Infinite();
-  }
-
-  if(Precision::IsNegativeInfinite(theMB))
-  {
-    if(Precision::IsPositiveInfinite(theMA))
-      return 0.0;
-    else
-      return -Precision::Infinite();
-  }
-
-  return (theMA + theMB);
 }
 
-static inline Standard_Real Multiplication(Standard_Real theMA, Standard_Real theMB)
+//=======================================================================
+//function : BRepGProp_Sinert
+//purpose  : Constructor
+//=======================================================================
+BRepGProp_Sinert::BRepGProp_Sinert (const BRepGProp_Face& theSurface,
+                                    const gp_Pnt&         theLocation)
 {
-  return (theMA * theMB);
+  SetLocation(theLocation);
+  Perform(theSurface);
 }
 
-static inline Standard_Real Addition(Standard_Real theMA, Standard_Real theMB)
+//=======================================================================
+//function : BRepGProp_Sinert
+//purpose  : Constructor
+//=======================================================================
+BRepGProp_Sinert::BRepGProp_Sinert(BRepGProp_Face&   theSurface,
+                                   BRepGProp_Domain& theDomain,
+                                   const gp_Pnt&     theLocation)
 {
-  return (theMA + theMB);
+  SetLocation(theLocation);
+  Perform(theSurface, theDomain);
 }
 
-static Standard_Integer FillIntervalBounds(Standard_Real               A,
-                                           Standard_Real               B,
-                                           const TColStd_Array1OfReal& Knots, 
-                                           HMath_Vector&               VA,
-                                           HMath_Vector&               VB)
+//=======================================================================
+//function : BRepGProp_Sinert
+//purpose  : Constructor
+//=======================================================================
+BRepGProp_Sinert::BRepGProp_Sinert(BRepGProp_Face&     theSurface,
+                                   const gp_Pnt&       theLocation,
+                                   const Standard_Real theEps)
 {
-  Standard_Integer i = 1, iEnd = Knots.Upper(), j = 1, k = 1;
-  VA(j++) = A;
-  for(; i <= iEnd; i++){
-    Standard_Real kn = Knots(i);
-    if(A < kn)
-    {
-      if(kn < B) 
-      {
-        VA(j++) = VB(k++) = kn;
-      }
-      else 
-      {
-        break;
-      }
-    }
-  }
-  VB(k) = B;
-  return k;
-}
-
-static inline Standard_Integer MaxSubs(Standard_Integer n, Standard_Integer coeff = SUBS_POWER){
-  //  return n = IntegerLast()/coeff < n? IntegerLast(): n*coeff + 1;
-  return Min((n * coeff + 1),SM);
-}
-
-static Standard_Integer LFillIntervalBounds(Standard_Real               A,
-                                            Standard_Real               B,
-                                            const TColStd_Array1OfReal& Knots, 
-                                            const Standard_Integer      NumSubs)
+  SetLocation(theLocation);
+  Perform(theSurface, theEps);
+}
+
+//=======================================================================
+//function : BRepGProp_Sinert
+//purpose  : Constructor
+//=======================================================================
+BRepGProp_Sinert::BRepGProp_Sinert(BRepGProp_Face&     theSurface,
+                                   BRepGProp_Domain&   theDomain,
+                                   const gp_Pnt&       theLocation,
+                                   const Standard_Real theEps)
 {
-  Standard_Integer iEnd = MaxSubs(Knots.Upper()-1, NumSubs), jEnd = L1->Upper();
-  iEnd = Max(iEnd, Knots.Upper());
-  if(iEnd - 1 > jEnd){
-    L1    = new math_Vector(1,iEnd);
-    L2    = new math_Vector(1,iEnd);
-    DimL  = new math_Vector(1,iEnd);
-    ErrL  = new math_Vector(1,iEnd,0.0);
-    ErrUL = new math_Vector(1,iEnd,0.0);
-    IxL   = new math_Vector(1,iEnd);
-    IyL   = new math_Vector(1,iEnd);
-    IzL   = new math_Vector(1,iEnd);
-    IxxL  = new math_Vector(1,iEnd);
-    IyyL  = new math_Vector(1,iEnd);
-    IzzL  = new math_Vector(1,iEnd);
-    IxyL  = new math_Vector(1,iEnd);
-    IxzL  = new math_Vector(1,iEnd);
-    IyzL  = new math_Vector(1,iEnd);
-  }
-  return FillIntervalBounds(A, B, Knots, L1, L2);
+  SetLocation(theLocation);
+  Perform(theSurface, theDomain, theEps);
 }
 
-static Standard_Integer UFillIntervalBounds(Standard_Real               A,
-                                            Standard_Real               B,
-                                            const TColStd_Array1OfReal& Knots, 
-                                            const Standard_Integer      NumSubs)
+//=======================================================================
+//function : SetLocation
+//purpose  : 
+//=======================================================================
+void BRepGProp_Sinert::SetLocation(const gp_Pnt& theLocation)
 {
-  Standard_Integer iEnd = MaxSubs(Knots.Upper()-1, NumSubs), jEnd = U1->Upper();
-  iEnd = Max(iEnd, Knots.Upper());
-  if(iEnd - 1 > jEnd){
-    U1   = new math_Vector(1,iEnd);
-    U2   = new math_Vector(1,iEnd);
-    DimU = new math_Vector(1,iEnd);
-    ErrU = new math_Vector(1,iEnd,0.0);
-    IxU  = new math_Vector(1,iEnd);
-    IyU  = new math_Vector(1,iEnd);
-    IzU  = new math_Vector(1,iEnd);
-    IxxU = new math_Vector(1,iEnd);
-    IyyU = new math_Vector(1,iEnd);
-    IzzU = new math_Vector(1,iEnd);
-    IxyU = new math_Vector(1,iEnd);
-    IxzU = new math_Vector(1,iEnd);
-    IyzU = new math_Vector(1,iEnd);
-  }
-  return FillIntervalBounds(A, B, Knots, U1, U2);
+  loc = theLocation;
 }
 
-static Standard_Real CCompute(BRepGProp_Face&                  S,
-                              BRepGProp_Domain&                D,
-                              const gp_Pnt&          loc,
-                              Standard_Real&         Dim,
-                              gp_Pnt&                g,
-                              gp_Mat&                inertia,
-                              const Standard_Real    EpsDim,
-                              const Standard_Boolean isErrorCalculation,
-                              const Standard_Boolean isVerifyComputation) 
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void BRepGProp_Sinert::Perform(const BRepGProp_Face& theSurface)
 {
-  Standard_Real  (*FuncAdd)(Standard_Real, Standard_Real);
-  Standard_Real  (*FuncMul)(Standard_Real, Standard_Real);
-
-  FuncAdd = Addition;
-  FuncMul = Multiplication;
-
-  Standard_Boolean isNaturalRestriction = S.NaturalRestriction();
-
-  Standard_Integer NumSubs = SUBS_POWER;
-
-  Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
-  Dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
-  Standard_Real x, y, z;
-  //boundary curve parametrization
-  Standard_Real l1, l2, lm, lr, l;   
-  //BRepGProp_Face parametrization in U and V direction
-  Standard_Real BV1, BV2, v;         
-  Standard_Real BU1, BU2, u1, u2, um, ur, u;
-  S.Bounds (BU1, BU2, BV1, BV2);
-  u1 = BU1;
-
-  if(Precision::IsInfinite(BU1) || Precision::IsInfinite(BU2) ||
-     Precision::IsInfinite(BV1) || Precision::IsInfinite(BV2))
-  {
-    FuncAdd = AdditionInf;
-    FuncMul = MultiplicationInf;
-  }
-
-
-  //location point used to compute the inertia
-  Standard_Real xloc, yloc, zloc;
-  loc.Coord (xloc, yloc, zloc); // use member of parent class
-  //Jacobien (x, y, z) -> (u, v) = ||n||
-  Standard_Real ds;                  
-  //On the BRepGProp_Face
-  gp_Pnt Ps;                    
-  gp_Vec VNor;
-  //On the boundary curve u-v
-  gp_Pnt2d Puv;                
-  gp_Vec2d Vuv;
-  Standard_Real Dul;  // Dul = Du / Dl
-  Standard_Real CDim[2], CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz;
-  Standard_Real LocDim[2], LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy, LocIxz, LocIyz;
-
-  Standard_Real ErrorU, ErrorL, ErrorLMax = 0.0, Eps=0.0, EpsL=0.0, EpsU=0.0;
-
-  Standard_Integer iD = 0, NbLSubs, iLS, iLSubEnd, iGL, iGLEnd, NbLGaussP[2], LRange[2], iL, kL, kLEnd, IL, JL;
-  Standard_Integer i, NbUSubs, iUS, iUSubEnd, iGU, iGUEnd, NbUGaussP[2], URange[2], iU, kU, kUEnd, IU, JU;
-  Standard_Integer UMaxSubs, LMaxSubs;
-  iGLEnd = isErrorCalculation? 2: 1; 
-  for(i = 0; i < 2; i++) {
-    LocDim[i] = 0.0;
-    CDim[i] = 0.0;
-  }
-
-  NbUGaussP[0] = S.SIntOrder(EpsDim);  
-  NbUGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbUGaussP[0]));
-  math::GaussPoints(NbUGaussP[0],UGaussP0);  math::GaussWeights(NbUGaussP[0],UGaussW0);
-  math::GaussPoints(NbUGaussP[1],UGaussP1);  math::GaussWeights(NbUGaussP[1],UGaussW1);
-
-  NbUSubs = S.SUIntSubs();
-  TColStd_Array1OfReal UKnots(1,NbUSubs+1);
-  S.UKnots(UKnots);
-
-  while (isNaturalRestriction || D.More())
-  {
-    if(isNaturalRestriction)
-    { 
-      NbLGaussP[0] = Min(2*NbUGaussP[0],math::GaussPointsMax());
-    }
-    else
-    {
-      S.Load(D.Value());  ++iD;
-      NbLGaussP[0] = S.LIntOrder(EpsDim);  
-    }
-
-    NbLGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbLGaussP[0]));
-    math::GaussPoints(NbLGaussP[0],LGaussP0);  math::GaussWeights(NbLGaussP[0],LGaussW0);
-    math::GaussPoints(NbLGaussP[1],LGaussP1);  math::GaussWeights(NbLGaussP[1],LGaussW1);
-
-    NbLSubs = isNaturalRestriction? S.SVIntSubs(): S.LIntSubs();
-
-    TColStd_Array1OfReal LKnots(1,NbLSubs+1);
-    if(isNaturalRestriction)
-    {
-      S.VKnots(LKnots); 
-      l1 = BV1; l2 = BV2;
-    }
-    else
-    {
-      S.LKnots(LKnots);
-      l1 = S.FirstParameter();  l2 = S.LastParameter();
-    }
-
-    ErrorL = 0.0;
-    kLEnd = 1; JL = 0;
-    //OCC503(apo): if(Abs(l2-l1) < EPS_PARAM) continue;
-    if(Abs(l2-l1) > EPS_PARAM)
-    {
-      iLSubEnd = LFillIntervalBounds(l1, l2, LKnots, NumSubs);
-      LMaxSubs = MaxSubs(iLSubEnd);
-      if(LMaxSubs > DimL.Vector()->Upper())
-        LMaxSubs = DimL.Vector()->Upper();
-
-      DimL.Init(0.0,1,LMaxSubs);  ErrL.Init(0.0,1,LMaxSubs);  ErrUL.Init(0.0,1,LMaxSubs);
-      
-      do{// while: L
-        if(++JL > iLSubEnd)
-        {
-          LRange[0] = IL = ErrL->Max();
-          LRange[1] = JL;
-          L1(JL) = (L1(IL) + L2(IL))/2.0;
-          L2(JL) = L2(IL);
-          L2(IL) = L1(JL);
-        }
-        else
-          LRange[0] = IL = JL;
-        
-        if(JL == LMaxSubs || Abs(L2(JL) - L1(JL)) < EPS_PARAM)
-          if(kLEnd == 1)
-          {
-            DimL(JL) = ErrL(JL) = IxL(JL) = IyL(JL) = IzL(JL) = 
-              IxxL(JL) = IyyL(JL) = IzzL(JL) = IxyL(JL) = IxzL(JL) = IyzL(JL) = 0.0;
-          }else{
-            JL--;
-            EpsL = ErrorL;  Eps = EpsL/0.9;
-            break;
-          }
-        else
-          for(kL=0; kL < kLEnd; kL++)
-          {
-            iLS = LRange[kL];
-            lm = 0.5*(L2(iLS) + L1(iLS));         
-            lr = 0.5*(L2(iLS) - L1(iLS));
-            CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
-            
-            for(iGL=0; iGL < iGLEnd; iGL++)
-            {
-              CDim[iGL] = 0.0;
-              for(iL=1; iL<=NbLGaussP[iGL]; iL++)
-              {
-                l = lm + lr*(*LGaussP[iGL])(iL);
-                if(isNaturalRestriction)
-                {
-                  v = l; u2 = BU2; Dul = (*LGaussW[iGL])(iL);
-                }
-                else
-                {
-                  S.D12d (l, Puv, Vuv);
-                  Dul = Vuv.Y()*(*LGaussW[iGL])(iL);  // Dul = Du / Dl
-                  if(Abs(Dul) < EPS_PARAM)
-                    continue;
-
-                  v  = Puv.Y();
-                  u2 = Puv.X();
-                  //Check on cause out off bounds of value current parameter
-                  if(v < BV1)
-                    v = BV1;
-                  else if(v > BV2)
-                    v = BV2;
-
-                  if(u2 < BU1)
-                    u2 = BU1;
-                  else if(u2 > BU2)
-                    u2 = BU2;
-                }
-
-                ErrUL(iLS) = 0.0;
-                kUEnd = 1;
-                JU = 0;
-
-                if(Abs(u2-u1) < EPS_PARAM)
-                  continue;
-
-                iUSubEnd = UFillIntervalBounds(u1, u2, UKnots, NumSubs);
-                UMaxSubs = MaxSubs(iUSubEnd);
-                if(UMaxSubs > DimU.Vector()->Upper())
-                  UMaxSubs = DimU.Vector()->Upper();
-
-                DimU.Init(0.0,1,UMaxSubs);  ErrU.Init(0.0,1,UMaxSubs);  ErrorU = 0.0;
-                
-                do{//while: U
-                  if(++JU > iUSubEnd)
-                  {
-                    URange[0] = IU = ErrU->Max();
-                    URange[1] = JU;
-                    U1(JU) = (U1(IU)+U2(IU))/2.0;
-                    U2(JU) = U2(IU);
-                    U2(IU) = U1(JU);
-                  }
-                  else
-                    URange[0] = IU = JU;
-
-                  if(JU == UMaxSubs || Abs(U2(JU) - U1(JU)) < EPS_PARAM)
-                    if(kUEnd == 1)
-                    {
-                      DimU(JU) = ErrU(JU) = IxU(JU) = IyU(JU) = IzU(JU) = 
-                        IxxU(JU) = IyyU(JU) = IzzU(JU) = IxyU(JU) = IxzU(JU) = IyzU(JU) = 0.0;
-                    }else
-                    {
-                      JU--;  
-                      EpsU = ErrorU;  Eps = EpsU*Abs((u2-u1)*Dul)/0.1;  EpsL = 0.9*Eps;
-                      break;
-                    }
-                  else
-                    for(kU=0; kU < kUEnd; kU++)
-                    {
-                      iUS = URange[kU];
-                      um = 0.5*(U2(iUS) + U1(iUS));
-                      ur = 0.5*(U2(iUS) - U1(iUS));
-                      LocIx = LocIy = LocIz = LocIxx = LocIyy = LocIzz = LocIxy = LocIxz = LocIyz = 0.0;
-                      iGUEnd = iGLEnd - iGL;
-                      for(iGU=0; iGU < iGUEnd; iGU++)
-                      {
-                        LocDim[iGU] = 0.0;
-                        for(iU=1; iU <= NbUGaussP[iGU]; iU++)
-                        {
-                          u = um + ur*(*UGaussP[iGU])(iU);
-                          S.Normal(u, v, Ps, VNor);
-                          ds = VNor.Magnitude();    //Jacobien(x,y,z) -> (u,v)=||n||
-                          ds *= (*UGaussW[iGU])(iU); 
-                          LocDim[iGU] += ds; 
-
-                          if(iGU > 0)
-                            continue;
-
-                          Ps.Coord(x, y, z);  
-                          x = FuncAdd(x, -xloc);
-                          y = FuncAdd(y, -yloc);
-                          z = FuncAdd(z, -zloc);
-
-                          const Standard_Real XdS = FuncMul(x, ds);
-                          const Standard_Real YdS = FuncMul(y, ds);
-                          const Standard_Real ZdS = FuncMul(z, ds);
-
-                          LocIx = FuncAdd(LocIx, XdS);
-                          LocIy = FuncAdd(LocIy, YdS);
-                          LocIz = FuncAdd(LocIz, ZdS);
-                          LocIxy = FuncAdd(LocIxy, FuncMul(x, YdS));
-                          LocIyz = FuncAdd(LocIyz, FuncMul(y, ZdS));
-                          LocIxz = FuncAdd(LocIxz, FuncMul(x, ZdS));
-
-                          const Standard_Real XXdS = FuncMul(x, XdS);
-                          const Standard_Real YYdS = FuncMul(y, YdS);
-                          const Standard_Real ZZdS = FuncMul(z, ZdS);
-
-                          LocIxx = FuncAdd(LocIxx, FuncAdd(YYdS, ZZdS));
-                          LocIyy = FuncAdd(LocIyy, FuncAdd(XXdS, ZZdS));
-                          LocIzz = FuncAdd(LocIzz, FuncAdd(XXdS, YYdS));
-                        }//for: iU
-                      }//for: iGU
-
-                      DimU(iUS) = FuncMul(LocDim[0],ur);
-                      if(iGL > 0)
-                        continue;
-
-                      ErrU(iUS) = FuncMul(Abs(LocDim[1]-LocDim[0]), ur);
-                      IxU(iUS) = FuncMul(LocIx, ur);
-                      IyU(iUS) = FuncMul(LocIy, ur);
-                      IzU(iUS) = FuncMul(LocIz, ur);
-                      IxxU(iUS) = FuncMul(LocIxx, ur);
-                      IyyU(iUS) = FuncMul(LocIyy, ur);
-                      IzzU(iUS) = FuncMul(LocIzz, ur);
-                      IxyU(iUS) = FuncMul(LocIxy, ur);
-                      IxzU(iUS) = FuncMul(LocIxz, ur);
-                      IyzU(iUS) = FuncMul(LocIyz, ur);
-                    }//for: kU (iUS)
-
-                    if(JU == iUSubEnd)
-                      kUEnd = 2;
-
-                    if(kUEnd == 2)
-                      ErrorU = ErrU(ErrU->Max());
-                }while((ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1);
-
-                for(i=1; i<=JU; i++)
-                  CDim[iGL] = FuncAdd(CDim[iGL], FuncMul(DimU(i), Dul));
-
-                if(iGL > 0)
-                  continue;
-
-                ErrUL(iLS) = ErrorU*Abs((u2-u1)*Dul);
-                for(i=1; i<=JU; i++)
-                {
-                  CIx = FuncAdd(CIx, FuncMul(IxU(i), Dul));
-                  CIy = FuncAdd(CIy, FuncMul(IyU(i), Dul));
-                  CIz = FuncAdd(CIz, FuncMul(IzU(i), Dul));
-                  CIxx = FuncAdd(CIxx, FuncMul(IxxU(i), Dul));
-                  CIyy = FuncAdd(CIyy, FuncMul(IyyU(i), Dul));
-                  CIzz = FuncAdd(CIzz, FuncMul(IzzU(i), Dul));
-                  CIxy = FuncAdd(CIxy, FuncMul(IxyU(i), Dul));
-                  CIxz = FuncAdd(CIxz, FuncMul(IxzU(i), Dul));
-                  CIyz = FuncAdd(CIyz, FuncMul(IyzU(i), Dul));
-                }
-              }//for: iL 
-            }//for: iGL
-
-            DimL(iLS) = FuncMul(CDim[0], lr);  
-            if(iGLEnd == 2)
-              ErrL(iLS) = FuncAdd(FuncMul(Abs(CDim[1]-CDim[0]),lr), ErrUL(iLS));
-
-            IxL(iLS) = FuncMul(CIx, lr);
-            IyL(iLS) = FuncMul(CIy, lr);
-            IzL(iLS) = FuncMul(CIz, lr); 
-            IxxL(iLS) = FuncMul(CIxx, lr);
-            IyyL(iLS) = FuncMul(CIyy, lr);
-            IzzL(iLS) = FuncMul(CIzz, lr); 
-            IxyL(iLS) = FuncMul(CIxy, lr);
-            IxzL(iLS) = FuncMul(CIxz, lr);
-            IyzL(iLS) = FuncMul(CIyz, lr);
-          }//for: (kL)iLS
-          //  Calculate/correct epsilon of computation by current value of Dim
-          //That is need for not spend time for 
-          if(JL == iLSubEnd)
-          {  
-            kLEnd = 2; 
-            Standard_Real DDim = 0.0;
-            for(i=1; i<=JL; i++)
-              DDim += DimL(i);
-
-            DDim = Abs(DDim*EpsDim);
-            if(DDim > Eps)
-            { 
-              Eps = DDim;
-              EpsL = 0.9*Eps;
-            }
-          }
-
-          if(kLEnd == 2)
-            ErrorL = ErrL(ErrL->Max());
-        }while((ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1);
-
-        for(i=1; i<=JL; i++)
-        {
-          Dim = FuncAdd(Dim, DimL(i)); 
-          Ix = FuncAdd(Ix, IxL(i));
-          Iy = FuncAdd(Iy, IyL(i));
-          Iz = FuncAdd(Iz, IzL(i));
-          Ixx = FuncAdd(Ixx, IxxL(i));
-          Iyy = FuncAdd(Iyy, IyyL(i));
-          Izz = FuncAdd(Izz, IzzL(i));
-          Ixy = FuncAdd(Ixy, IxyL(i));
-          Ixz = FuncAdd(Ixz, IxzL(i));
-          Iyz = FuncAdd(Iyz, IyzL(i));
-        }
-
-        ErrorLMax = Max(ErrorLMax, ErrorL);
-      }
-
-      if(isNaturalRestriction)
-        break;
-
-      D.Next();
-    }
-
-    if(Abs(Dim) >= EPS_DIM)
-    {
-      Ix /= Dim;
-      Iy /= Dim;
-      Iz /= Dim;
-      g.SetCoord (Ix, Iy, Iz);
-    }
-    else
-    {
-      Dim =0.0;
-      g.SetCoord (0., 0.,0.);
-    }
-
-    inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
-                      gp_XYZ (-Ixy, Iyy, -Iyz),
-                      gp_XYZ (-Ixz, -Iyz, Izz));
-
-    if(iGLEnd == 2)
-      Eps = Dim != 0.0? ErrorLMax/Abs(Dim): 0.0;
-    else
-      Eps = EpsDim;
-
-    return Eps;
-}
-
-static Standard_Real Compute(BRepGProp_Face& S, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, 
-                             Standard_Real EpsDim) 
-{
-  Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
-  Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
-  EpsDim = Abs(EpsDim);
-  BRepGProp_Domain D;
-  return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
-}
+  myEpsilon = 1.0;
 
-static Standard_Real Compute(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, 
-                             Standard_Real EpsDim) 
-{
-  Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
-  Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
-  EpsDim = Abs(EpsDim);
-  return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
+  BRepGProp_Gauss aGauss(BRepGProp_Gauss::Sinert);
+  aGauss.Compute(theSurface, loc, dim, g, inertia);
 }
 
-static void Compute(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pnt& loc, Standard_Real& dim, gp_Pnt& g, gp_Mat& inertia)
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void BRepGProp_Sinert::Perform(BRepGProp_Face&   theSurface,
+                               BRepGProp_Domain& theDomain)
 {
-  Standard_Real  (*FuncAdd)(Standard_Real, Standard_Real);
-  Standard_Real  (*FuncMul)(Standard_Real, Standard_Real);
-
-  FuncAdd = Addition;
-  FuncMul = Multiplication;
-
-  Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
-  dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
-
-  Standard_Real x, y, z;
-  Standard_Integer NbCGaussgp_Pnts = 0;
-
-  Standard_Real l1, l2, lm, lr, l;   //boundary curve parametrization
-  Standard_Real v1, v2,         v;   //BRepGProp_Face parametrization in v direction
-  Standard_Real u1, u2, um, ur, u;
-  Standard_Real ds;                  //Jacobien (x, y, z) -> (u, v) = ||n||
-
-  gp_Pnt P;                    //On the BRepGProp_Face
-  gp_Vec VNor;
-
-  gp_Pnt2d Puv;                //On the boundary curve u-v
-  gp_Vec2d Vuv;
-  Standard_Real Dul;                 // Dul = Du / Dl
-  Standard_Real CArea, CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz;
-  Standard_Real LocArea, LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy,
-    LocIxz, LocIyz;
-
-
-  S.Bounds (u1, u2, v1, v2);
-
-  if(Precision::IsInfinite(u1) || Precision::IsInfinite(u2) ||
-     Precision::IsInfinite(v1) || Precision::IsInfinite(v2))
-  {
-    FuncAdd = AdditionInf;
-    FuncMul = MultiplicationInf;
-  }
-
-
-  Standard_Integer NbUGaussgp_Pnts = Min(S.UIntegrationOrder (),
-    math::GaussPointsMax());
-  Standard_Integer NbVGaussgp_Pnts = Min(S.VIntegrationOrder (),
-    math::GaussPointsMax());
-
-  Standard_Integer NbGaussgp_Pnts = Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts);
-
-  //Number of Gauss points for the integration
-  //on the BRepGProp_Face
-  math_Vector GaussSPV (1, NbGaussgp_Pnts);
-  math_Vector GaussSWV (1, NbGaussgp_Pnts);
-  math::GaussPoints  (NbGaussgp_Pnts,GaussSPV);
-  math::GaussWeights (NbGaussgp_Pnts,GaussSWV);
-
-
-  //location point used to compute the inertia
-  Standard_Real xloc, yloc, zloc;
-  loc.Coord (xloc, yloc, zloc);
-
-  while (D.More()) {
-
-    S.Load(D.Value());
-    NbCGaussgp_Pnts =  Min(S.IntegrationOrder (), math::GaussPointsMax());        
-
-    math_Vector GaussCP (1, NbCGaussgp_Pnts);
-    math_Vector GaussCW (1, NbCGaussgp_Pnts);
-    math::GaussPoints  (NbCGaussgp_Pnts,GaussCP);
-    math::GaussWeights (NbCGaussgp_Pnts,GaussCW);
-
-    CArea = 0.0;
-    CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
-    l1 = S.FirstParameter ();
-    l2 = S.LastParameter  ();
-    lm = 0.5 * (l2 + l1);         
-    lr = 0.5 * (l2 - l1);
-
-    for (Standard_Integer i = 1; i <= NbCGaussgp_Pnts; i++) {
-      l = lm + lr * GaussCP (i);
-      S.D12d(l, Puv, Vuv);
-      v   = Puv.Y();
-      u2  = Puv.X();
-      Dul = Vuv.Y();
-      Dul *= GaussCW (i);
-      um  = 0.5 * (u2 + u1);
-      ur  = 0.5 * (u2 - u1);
-      LocArea = LocIx  = LocIy  = LocIz = LocIxx = LocIyy = LocIzz = 
-        LocIxy  = LocIxz = LocIyz = 0.0;
-      for (Standard_Integer j = 1; j <= NbGaussgp_Pnts; j++) {
-        u = FuncAdd(um, FuncMul(ur, GaussSPV (j)));
-        S.Normal (u, v, P, VNor);
-        ds = VNor.Magnitude();    //normal.Magnitude
-        ds = FuncMul(ds, Dul) * GaussSWV (j); 
-        LocArea = FuncAdd(LocArea, ds); 
-        P.Coord (x, y, z);
-
-        x = FuncAdd(x, -xloc);
-        y = FuncAdd(y, -yloc);
-        z = FuncAdd(z, -zloc);
-
-        const Standard_Real XdS = FuncMul(x, ds);
-        const Standard_Real YdS = FuncMul(y, ds);
-        const Standard_Real ZdS = FuncMul(z, ds);
-        
-        LocIx = FuncAdd(LocIx, XdS);
-        LocIy = FuncAdd(LocIy, YdS);
-        LocIz = FuncAdd(LocIz, ZdS);
-        LocIxy = FuncAdd(LocIxy, FuncMul(x, YdS));
-        LocIyz = FuncAdd(LocIyz, FuncMul(y, ZdS));
-        LocIxz = FuncAdd(LocIxz, FuncMul(x, ZdS));
-
-        const Standard_Real XXdS = FuncMul(x, XdS);
-        const Standard_Real YYdS = FuncMul(y, YdS);
-        const Standard_Real ZZdS = FuncMul(z, ZdS);
-        
-        LocIxx = FuncAdd(LocIxx, FuncAdd(YYdS, ZZdS));
-        LocIyy = FuncAdd(LocIyy, FuncAdd(XXdS, ZZdS));
-        LocIzz = FuncAdd(LocIzz, FuncAdd(XXdS, YYdS));
-      }
-
-      CArea = FuncAdd(CArea, FuncMul(LocArea, ur));
-      CIx   = FuncAdd(CIx, FuncMul(LocIx, ur));
-      CIy   = FuncAdd(CIy, FuncMul(LocIy, ur));
-      CIz   = FuncAdd(CIz, FuncMul(LocIz, ur));
-      CIxx  = FuncAdd(CIxx, FuncMul(LocIxx, ur));
-      CIyy  = FuncAdd(CIyy, FuncMul(LocIyy, ur));
-      CIzz  = FuncAdd(CIzz, FuncMul(LocIzz, ur));
-      CIxy  = FuncAdd(CIxy, FuncMul(LocIxy, ur));
-      CIxz  = FuncAdd(CIxz, FuncMul(LocIxz, ur));
-      CIyz  = FuncAdd(CIyz, FuncMul(LocIyz, ur));
-    }
-
-    dim = FuncAdd(dim, FuncMul(CArea, lr));
-    Ix  = FuncAdd(Ix,  FuncMul(CIx, lr));
-    Iy  = FuncAdd(Iy,  FuncMul(CIy, lr));
-    Iz  = FuncAdd(Iz,  FuncMul(CIz, lr));
-    Ixx = FuncAdd(Ixx, FuncMul(CIxx, lr));
-    Iyy = FuncAdd(Iyy, FuncMul(CIyy, lr));
-    Izz = FuncAdd(Izz, FuncMul(CIzz, lr));
-    Ixy = FuncAdd(Ixy, FuncMul(CIxy, lr));
-    Ixz = FuncAdd(Iyz, FuncMul(CIxz, lr));
-    Iyz = FuncAdd(Ixz, FuncMul(CIyz, lr));
-    D.Next();
-  }
-
-  if (Abs(dim) >= EPS_DIM) {
-    Ix /= dim;
-    Iy /= dim;
-    Iz /= dim;
-    g.SetCoord (Ix, Iy, Iz);
-  }
-  else {
-    dim =0.;
-    g.SetCoord (0., 0.,0.);
-  }
+  myEpsilon = 1.0;
 
-  inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
-                    gp_XYZ (-Ixy, Iyy, -Iyz),
-                    gp_XYZ (-Ixz, -Iyz, Izz));
+  BRepGProp_Gauss aGauss(BRepGProp_Gauss::Sinert);
+  aGauss.Compute(theSurface, theDomain, loc, dim, g, inertia);
 }
 
-static void Compute(const BRepGProp_Face& S,
-                    const gp_Pnt& loc,
-                    Standard_Real& dim,
-                    gp_Pnt& g,
-                    gp_Mat& inertia)
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+Standard_Real BRepGProp_Sinert::Perform(BRepGProp_Face&     theSurface,
+                                        const Standard_Real theEps)
 {
-  Standard_Real  (*FuncAdd)(Standard_Real, Standard_Real);
-  Standard_Real  (*FuncMul)(Standard_Real, Standard_Real);
-
-  FuncAdd = Addition;
-  FuncMul = Multiplication;
-
-  Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
-  dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
-
-  Standard_Real LowerU, UpperU, LowerV, UpperV;
-  S.Bounds (LowerU, UpperU, LowerV, UpperV);
-
-  if(Precision::IsInfinite(LowerU) || Precision::IsInfinite(UpperU) ||
-     Precision::IsInfinite(LowerV) || Precision::IsInfinite(UpperV))
-  {
-    FuncAdd = AdditionInf;
-    FuncMul = MultiplicationInf;
-  }
-
-  Standard_Integer UOrder = Min(S.UIntegrationOrder (),
-    math::GaussPointsMax());
-  Standard_Integer VOrder = Min(S.VIntegrationOrder (),
-    math::GaussPointsMax());   
-  gp_Pnt P;          
-  gp_Vec VNor;   
-  Standard_Real dsi, ds;        
-  Standard_Real ur, um, u, vr, vm, v;
-  Standard_Real x, y, z; 
-  Standard_Real Ixi, Iyi, Izi, Ixxi, Iyyi, Izzi, Ixyi, Ixzi, Iyzi;
-  Standard_Real xloc, yloc, zloc;
-  loc.Coord (xloc, yloc, zloc);
-
-  Standard_Integer i, j;
-  math_Vector GaussPU (1, UOrder);     //gauss points and weights
-  math_Vector GaussWU (1, UOrder);
-  math_Vector GaussPV (1, VOrder);
-  math_Vector GaussWV (1, VOrder);
-
-  //Recuperation des points de Gauss dans le fichier GaussPoints.
-  math::GaussPoints  (UOrder,GaussPU);
-  math::GaussWeights (UOrder,GaussWU);
-  math::GaussPoints  (VOrder,GaussPV);
-  math::GaussWeights (VOrder,GaussWV);
-
-  // Calcul des integrales aux points de gauss :
-  um = 0.5 * FuncAdd(UpperU,  LowerU);
-  vm = 0.5 * FuncAdd(UpperV,  LowerV);
-  ur = 0.5 * FuncAdd(UpperU, -LowerU);
-  vr = 0.5 * FuncAdd(UpperV, -LowerV);
-
-  for (j = 1; j <= VOrder; j++) {
-    v = FuncAdd(vm, FuncMul(vr, GaussPV(j)));
-    dsi = Ixi = Iyi = Izi = Ixxi = Iyyi = Izzi = Ixyi = Ixzi = Iyzi = 0.0;
-
-    for (i = 1; i <= UOrder; i++) {
-      u = FuncAdd(um, FuncMul(ur, GaussPU (i)));
-      S.Normal (u, v, P, VNor); 
-      ds = FuncMul(VNor.Magnitude(), GaussWU (i));
-      P.Coord (x, y, z);
-
-      x = FuncAdd(x, -xloc);
-      y = FuncAdd(y, -yloc);
-      z = FuncAdd(z, -zloc);
-
-      dsi = FuncAdd(dsi, ds);
-
-      const Standard_Real XdS = FuncMul(x, ds);
-      const Standard_Real YdS = FuncMul(y, ds);
-      const Standard_Real ZdS = FuncMul(z, ds);
-
-      Ixi = FuncAdd(Ixi, XdS);
-      Iyi = FuncAdd(Iyi, YdS);
-      Izi = FuncAdd(Izi, ZdS);
-      Ixyi = FuncAdd(Ixyi, FuncMul(x, YdS));
-      Iyzi = FuncAdd(Iyzi, FuncMul(y, ZdS));
-      Ixzi = FuncAdd(Ixzi, FuncMul(x, ZdS));
-
-      const Standard_Real XXdS = FuncMul(x, XdS);
-      const Standard_Real YYdS = FuncMul(y, YdS);
-      const Standard_Real ZZdS = FuncMul(z, ZdS);
-
-      Ixxi = FuncAdd(Ixxi, FuncAdd(YYdS, ZZdS));
-      Iyyi = FuncAdd(Iyyi, FuncAdd(XXdS, ZZdS));
-      Izzi = FuncAdd(Izzi, FuncAdd(XXdS, YYdS));
-    }
-
-    dim   = FuncAdd(dim, FuncMul(dsi, GaussWV (j)));
-    Ix    = FuncAdd(Ix,  FuncMul(Ixi, GaussWV (j)));
-    Iy    = FuncAdd(Iy,  FuncMul(Iyi, GaussWV (j)));
-    Iz    = FuncAdd(Iz,  FuncMul(Izi, GaussWV (j)));
-    Ixx   = FuncAdd(Ixx, FuncMul(Ixxi, GaussWV (j)));
-    Iyy   = FuncAdd(Iyy, FuncMul(Iyyi, GaussWV (j)));
-    Izz   = FuncAdd(Izz, FuncMul(Izzi, GaussWV (j)));
-    Ixy   = FuncAdd(Ixy, FuncMul(Ixyi, GaussWV (j)));
-    Iyz   = FuncAdd(Iyz, FuncMul(Iyzi, GaussWV (j)));
-    Ixz   = FuncAdd(Ixz, FuncMul(Ixzi, GaussWV (j)));
-  }
-
-  vr  = FuncMul(vr, ur);
-  Ixx = FuncMul(vr, Ixx);
-  Iyy = FuncMul(vr, Iyy);
-  Izz = FuncMul(vr, Izz);
-  Ixy = FuncMul(vr, Ixy);
-  Ixz = FuncMul(vr, Ixz);
-  Iyz = FuncMul(vr, Iyz);
-
-  if (Abs(dim) >= EPS_DIM)
-  {
-    Ix    /= dim;
-    Iy    /= dim;
-    Iz    /= dim;
-    dim   *= vr;
-    g.SetCoord (Ix, Iy, Iz);
-  }
-  else
-  {
-    dim =0.;
-    g.SetCoord (0.,0.,0.);
-  }
-
-  inertia = gp_Mat (gp_XYZ ( Ixx, -Ixy, -Ixz),
-                    gp_XYZ (-Ixy,  Iyy, -Iyz),
-                    gp_XYZ (-Ixz, -Iyz,  Izz));
+  BRepGProp_Domain anEmptyDomian;
+  return Perform(theSurface, anEmptyDomian, theEps);
 }
 
-BRepGProp_Sinert::BRepGProp_Sinert(){}
-
-BRepGProp_Sinert::BRepGProp_Sinert (const BRepGProp_Face&   S,
-                              const gp_Pnt& SLocation
-                              ) 
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+Standard_Real BRepGProp_Sinert::Perform(BRepGProp_Face&     theSurface,
+                                        BRepGProp_Domain&   theDomain,
+                                        const Standard_Real theEps)
 {
-  SetLocation(SLocation);
-  Perform(S);
+  BRepGProp_Gauss  aGauss(BRepGProp_Gauss::Sinert);
+  return myEpsilon = aGauss.Compute(theSurface, theDomain, loc, theEps, dim, g, inertia);
 }
 
-BRepGProp_Sinert::BRepGProp_Sinert (BRepGProp_Face&   S,
-                              BRepGProp_Domain& D,
-                              const gp_Pnt& SLocation
-                              ) 
+//=======================================================================
+//function : GetEpsilon
+//purpose  : 
+//=======================================================================
+Standard_Real BRepGProp_Sinert::GetEpsilon()
 {
-  SetLocation(SLocation);
-  Perform(S,D);
-}
-
-BRepGProp_Sinert::BRepGProp_Sinert(BRepGProp_Face& S, const gp_Pnt& SLocation, const Standard_Real Eps){
-  SetLocation(SLocation);
-  Perform(S, Eps);
-}
-
-BRepGProp_Sinert::BRepGProp_Sinert(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pnt& SLocation, const Standard_Real Eps){
-  SetLocation(SLocation);
-  Perform(S, D, Eps);
-}
-
-void BRepGProp_Sinert::SetLocation(const gp_Pnt& SLocation){
-  loc = SLocation;
-}
-
-void BRepGProp_Sinert::Perform(const BRepGProp_Face& S){
-  Compute(S,loc,dim,g,inertia);
-  myEpsilon = 1.0;
-  return;
-}
-
-void BRepGProp_Sinert::Perform(BRepGProp_Face& S, BRepGProp_Domain& D){
-  Compute(S,D,loc,dim,g,inertia);
-  myEpsilon = 1.0;
-  return;
-}
-
-Standard_Real BRepGProp_Sinert::Perform(BRepGProp_Face& S, const Standard_Real Eps){
-  return myEpsilon = Compute(S,loc,dim,g,inertia,Eps);
-}
-
-Standard_Real BRepGProp_Sinert::Perform(BRepGProp_Face& S, BRepGProp_Domain& D, const Standard_Real Eps){
-  return myEpsilon = Compute(S,D,loc,dim,g,inertia,Eps);
-}
-
-
-Standard_Real BRepGProp_Sinert::GetEpsilon(){
   return myEpsilon;
 }
index 843f83f..524f1f2 100644 (file)
 // commercial license or contractual agreement.
 
 #include <BRepGProp_Vinert.ixx>
+#include <BRepGProp_Gauss.hxx>
 
-#include <math.hxx>
-#include <TColStd_Array1OfReal.hxx>
-
-class HMath_Vector{
-  math_Vector *pvec;
-  void operator=(const math_Vector&){}
-public:
-  HMath_Vector(){ pvec = 0;}
-  HMath_Vector(math_Vector* pv){ pvec = pv;}
-  ~HMath_Vector(){ if(pvec != 0) delete pvec;}
-  void operator=(math_Vector* pv){ if(pvec != pv && pvec != 0) delete pvec;  pvec = pv;}
-  Standard_Real& operator()(Standard_Integer i){ return (*pvec).operator()(i);}
-  const Standard_Real& operator()(Standard_Integer i) const{ return (*pvec).operator()(i);}
-  const math_Vector* operator->() const{ return pvec;}
-  math_Vector* operator->(){ return pvec;}
-  math_Vector* Init(Standard_Real v, Standard_Integer i = 0, Standard_Integer iEnd = 0){ 
-    if(pvec == 0) return pvec;
-    if(iEnd - i == 0) pvec->Init(v);
-    else for(; i <= iEnd; i++) pvec->operator()(i) = v;
-    return pvec;
-  }
-};
-
-//Minimal value of interval's range for computation | minimal value of "dim" | ... 
-static Standard_Real     EPS_PARAM = Precision::Angular(), EPS_DIM = 1.E-30, ERROR_ALGEBR_RATIO = 2.0/3.0;
-//Maximum of GaussPoints on a subinterval and maximum of subintervals
-static Standard_Integer  GPM = math::GaussPointsMax(), SUBS_POWER = 32, SM = SUBS_POWER*GPM + 1; 
-static Standard_Boolean  IS_MIN_DIM = 1; // if the value equal 0 error of algorithm calculted by static moments
-
-static math_Vector  LGaussP0(1,GPM), LGaussW0(1,GPM),
-LGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM))), LGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
-static HMath_Vector  L1 = new math_Vector(1,SM), L2 = new math_Vector(1,SM),
-DimL = new math_Vector(1,SM), ErrL = new math_Vector(1,SM), ErrUL = new math_Vector(1,SM,0.0),
-IxL = new math_Vector(1,SM), IyL = new math_Vector(1,SM), IzL = new math_Vector(1,SM),
-IxxL = new math_Vector(1,SM), IyyL = new math_Vector(1,SM), IzzL = new math_Vector(1,SM),
-IxyL = new math_Vector(1,SM), IxzL = new math_Vector(1,SM), IyzL = new math_Vector(1,SM);
-
-static math_Vector* LGaussP[] = {&LGaussP0,&LGaussP1};
-static math_Vector* LGaussW[] = {&LGaussW0,&LGaussW1};
-
-static math_Vector  UGaussP0(1,GPM), UGaussW0(1,GPM),
-UGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM))), UGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
-static HMath_Vector  U1 = new math_Vector(1,SM), U2 = new math_Vector(1,SM), 
-DimU = new math_Vector(1,SM), ErrU = new math_Vector(1,SM,0.0),
-IxU = new math_Vector(1,SM), IyU = new math_Vector(1,SM), IzU = new math_Vector(1,SM),
-IxxU = new math_Vector(1,SM), IyyU = new math_Vector(1,SM), IzzU = new math_Vector(1,SM),
-IxyU = new math_Vector(1,SM), IxzU = new math_Vector(1,SM), IyzU = new math_Vector(1,SM);
-
-static math_Vector* UGaussP[] = {&UGaussP0,&UGaussP1};
-static math_Vector* UGaussW[] = {&UGaussW0,&UGaussW1};
-
-static Standard_Integer FillIntervalBounds(Standard_Real A, Standard_Real B, const TColStd_Array1OfReal& Knots, 
-                                           HMath_Vector& VA, HMath_Vector& VB)
+//=======================================================================
+//function : BRepGProp_Vinert
+//purpose  : Constructor
+//=======================================================================
+BRepGProp_Vinert::BRepGProp_Vinert()
 {
-  Standard_Integer i = 1, iEnd = Knots.Upper(), j = 1, k = 1;
-  VA(j++) = A;
-  for(; i <= iEnd; i++){
-    Standard_Real kn = Knots(i);
-    if(A < kn)
-    {
-      if(kn < B) 
-      {
-        VA(j++) = VB(k++) = kn; 
-      }
-      else
-      {
-        break;
-      }
-    }
-  }
-  VB(k) = B;
-  return k;
 }
 
-static inline Standard_Integer MaxSubs(Standard_Integer n, Standard_Integer coeff = SUBS_POWER){
-  return n = IntegerLast()/coeff < n? IntegerLast(): n*coeff + 1; 
-}
-
-static Standard_Integer LFillIntervalBounds(Standard_Real A, Standard_Real B, const TColStd_Array1OfReal& Knots, 
-                                            const Standard_Integer NumSubs)
+//=======================================================================
+//function : BRepGProp_Vinert
+//purpose  : 
+//=======================================================================
+BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face&     theSurface,
+                                   const gp_Pnt&       theLocation,
+                                   const Standard_Real theEps)
 {
-  Standard_Integer iEnd = Knots.Upper(), jEnd = L1->Upper();
-
-  //  Modified by Sergey KHROMOV - Wed Mar 26 11:22:50 2003
-  iEnd = Max(iEnd, MaxSubs(iEnd-1,NumSubs));
-  if(iEnd - 1 > jEnd){ 
-    //     iEnd = MaxSubs(iEnd-1,NumSubs); 
-    //  Modified by Sergey KHROMOV - Wed Mar 26 11:22:51 2003
-    L1 = new math_Vector(1,iEnd);  L2 = new math_Vector(1,iEnd);
-    DimL = new math_Vector(1,iEnd); ErrL = new math_Vector(1,iEnd,0.0); ErrUL = new math_Vector(1,iEnd,0.0);
-    IxL = new math_Vector(1,iEnd); IyL = new math_Vector(1,iEnd); IzL = new math_Vector(1,iEnd);
-    IxxL = new math_Vector(1,iEnd); IyyL = new math_Vector(1,iEnd); IzzL = new math_Vector(1,iEnd);
-    IxyL = new math_Vector(1,iEnd); IxzL = new math_Vector(1,iEnd); IyzL = new math_Vector(1,iEnd);
-  }
-  return FillIntervalBounds(A, B, Knots, L1, L2);
-}
-
-static Standard_Integer UFillIntervalBounds(Standard_Real A, Standard_Real B, const TColStd_Array1OfReal& Knots,
-                                            const Standard_Integer NumSubs)
+  SetLocation(theLocation);
+  Perform(theSurface, theEps);
+}
+
+//=======================================================================
+//function : BRepGProp_Vinert
+//purpose  : 
+//=======================================================================
+BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face&     theSurface,
+                                   BRepGProp_Domain&   theDomain,
+                                   const gp_Pnt&       theLocation,
+                                   const Standard_Real theEps)
 {
-  Standard_Integer iEnd = Knots.Upper(), jEnd = U1->Upper();
-
-  //  Modified by Sergey KHROMOV - Wed Mar 26 11:22:50 2003
-  iEnd = Max(iEnd, MaxSubs(iEnd-1,NumSubs));
-  if(iEnd - 1 > jEnd){
-    //     iEnd = MaxSubs(iEnd-1,NumSubs); 
-    //  Modified by Sergey KHROMOV - Wed Mar 26 11:22:51 2003
-    U1 = new math_Vector(1,iEnd);  U2 = new math_Vector(1,iEnd);
-    DimU = new math_Vector(1,iEnd); ErrU = new math_Vector(1,iEnd,0.0);
-    IxU = new math_Vector(1,iEnd); IyU = new math_Vector(1,iEnd); IzU = new math_Vector(1,iEnd);
-    IxxU = new math_Vector(1,iEnd); IyyU = new math_Vector(1,iEnd); IzzU = new math_Vector(1,iEnd);
-    IxyU = new math_Vector(1,iEnd); IxzU = new math_Vector(1,iEnd); IyzU = new math_Vector(1,iEnd);
-  }
-  return FillIntervalBounds(A, B, Knots, U1, U2);
+  SetLocation(theLocation);
+  Perform(theSurface, theDomain, theEps);
 }
 
-static Standard_Real CCompute(BRepGProp_Face& S, BRepGProp_Domain& D, const Standard_Boolean ByPoint, const Standard_Real Coeff[],
-                              const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, 
-                              const Standard_Real EpsDim, 
-                              const Standard_Boolean isErrorCalculation, const Standard_Boolean isVerifyComputation) 
+//=======================================================================
+//function : BRepGProp_Vinert
+//purpose  : 
+//=======================================================================
+BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face&   theSurface,
+                                   BRepGProp_Domain& theDomain,
+                                   const gp_Pnt&     theLocation)
 {
-  Standard_Boolean isNaturalRestriction = S.NaturalRestriction();
-
-  Standard_Integer NumSubs = SUBS_POWER;
-  Standard_Boolean isMinDim = IS_MIN_DIM;
-
-  Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
-  Dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
-  //boundary curve parametrization
-  Standard_Real l1, l2, lm, lr, l;   
-  //BRepGProp_Face parametrization in U and V direction
-  Standard_Real BV1, BV2, v;         
-  Standard_Real BU1, BU2, u1, u2, um, ur, u;
-  S.Bounds (BU1, BU2, BV1, BV2);  u1 = BU1;
-  //location point used to compute the inertia
-  Standard_Real xloc, yloc, zloc;
-  loc.Coord (xloc, yloc, zloc);
-  //location point used to compute the inertiard (xloc, yloc, zloc);
-  //Jacobien (x, y, z) -> (u, v) = ||n||
-  Standard_Real xn, yn, zn, s, ds, dDim;
-  Standard_Real x, y, z, xi, px, py, pz, yi, zi, d1, d2, d3;
-  //On the BRepGProp_Face
-  gp_Pnt Ps;                    
-  gp_Vec VNor;
-  //On the boundary curve u-v
-  gp_Pnt2d Puv;                
-  gp_Vec2d Vuv;
-  Standard_Real Dul;  // Dul = Du / Dl
-  Standard_Real CDim[2], CIx, CIy, CIz, CIxx[2], CIyy[2], CIzz[2], CIxy, CIxz, CIyz;
-  Standard_Real LocDim[2], LocIx[2], LocIy[2], LocIz[2], LocIxx[2], LocIyy[2], LocIzz[2], LocIxy[2], LocIxz[2], LocIyz[2];
-
-  Standard_Integer iD = 0, NbLSubs, iLS, iLSubEnd, iGL, iGLEnd, NbLGaussP[2], LRange[2], iL, kL, kLEnd, IL, JL;
-  Standard_Integer i, NbUSubs, iUS, iUSubEnd, iGU, iGUEnd, NbUGaussP[2], URange[2], iU, kU, kUEnd, IU, JU;
-  Standard_Integer UMaxSubs, LMaxSubs;
-
-  Standard_Real ErrorU, ErrorL, ErrorLMax = 0.0, Eps=0.0, EpsL=0.0, EpsU=0.0;
-  iGLEnd = isErrorCalculation? 2: 1; 
-
-  for(i = 0; i < 2; i++) {
-    LocDim[i] = 0.0;
-    LocIx[i] = 0.0;
-    LocIy[i] = 0.0;
-    LocIz[i] = 0.0;
-    LocIxx[i] = 0.0;
-    LocIyy[i] = 0.0;
-    LocIzz[i] = 0.0;
-    LocIxy[i] = 0.0;
-    LocIyz[i] = 0.0;
-    LocIxz[i] = 0.0;
-  }
-
-  NbUGaussP[0] = S.SIntOrder(EpsDim);  
-  NbUGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbUGaussP[0]));
-  math::GaussPoints(NbUGaussP[0],UGaussP0);  math::GaussWeights(NbUGaussP[0],UGaussW0);
-  math::GaussPoints(NbUGaussP[1],UGaussP1);  math::GaussWeights(NbUGaussP[1],UGaussW1);
-
-  NbUSubs = S.SUIntSubs();
-  TColStd_Array1OfReal UKnots(1,NbUSubs+1);
-  S.UKnots(UKnots);
-
-  while (isNaturalRestriction || D.More()) {
-    if(isNaturalRestriction){ 
-      NbLGaussP[0] = Min(2*NbUGaussP[0],math::GaussPointsMax());
-    }else{
-      S.Load(D.Value());  ++iD;
-      NbLGaussP[0] = S.LIntOrder(EpsDim);  
-    }
-    NbLGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbLGaussP[0]));
-    math::GaussPoints(NbLGaussP[0],LGaussP0);  math::GaussWeights(NbLGaussP[0],LGaussW0);
-    math::GaussPoints(NbLGaussP[1],LGaussP1);  math::GaussWeights(NbLGaussP[1],LGaussW1);
-
-    NbLSubs = isNaturalRestriction? S.SVIntSubs(): S.LIntSubs();
-    TColStd_Array1OfReal LKnots(1,NbLSubs+1);
-    if(isNaturalRestriction){
-      S.VKnots(LKnots); 
-      l1 = BV1; l2 = BV2;
-    }else{
-      S.LKnots(LKnots);
-      l1 = S.FirstParameter();  l2 = S.LastParameter();
-    }
-    ErrorL = 0.0;
-    kLEnd = 1; JL = 0;
-    //OCC503(apo): if(Abs(l2-l1) < EPS_PARAM) continue;
-    if(Abs(l2-l1) > EPS_PARAM) {
-      iLSubEnd = LFillIntervalBounds(l1, l2, LKnots, NumSubs);
-      LMaxSubs = MaxSubs(iLSubEnd);
-      //-- exception avoiding
-      if(LMaxSubs > SM) LMaxSubs = SM;
-      DimL.Init(0.0,1,LMaxSubs);  ErrL.Init(0.0,1,LMaxSubs);  ErrUL.Init(0.0,1,LMaxSubs);
-      do{// while: L
-        if(++JL > iLSubEnd){
-          LRange[0] = IL = ErrL->Max();  LRange[1] = JL;
-          L1(JL) = (L1(IL) + L2(IL))/2.0;  L2(JL) = L2(IL);  L2(IL) = L1(JL);
-        }else  LRange[0] = IL = JL;
-        if(JL == LMaxSubs || Abs(L2(JL) - L1(JL)) < EPS_PARAM)
-          if(kLEnd == 1){
-            DimL(JL) = ErrL(JL) = IxL(JL) = IyL(JL) = IzL(JL) = 
-              IxxL(JL) = IyyL(JL) = IzzL(JL) = IxyL(JL) = IxzL(JL) = IyzL(JL) = 0.0;
-          }else{
-            JL--;
-            EpsL = ErrorL;  Eps = EpsL/0.9;
-            break;
-          }
-        else
-          for(kL=0; kL < kLEnd; kL++){
-            iLS = LRange[kL];
-            lm = 0.5*(L2(iLS) + L1(iLS));         
-            lr = 0.5*(L2(iLS) - L1(iLS));
-            CIx = CIy = CIz = CIxy = CIxz = CIyz = 0.0;
-            for(iGL=0; iGL < iGLEnd; iGL++){//
-              CDim[iGL] = CIxx[iGL] = CIyy[iGL] = CIzz[iGL] = 0.0;
-              for(iL=1; iL<=NbLGaussP[iGL]; iL++){
-                l = lm + lr*(*LGaussP[iGL])(iL);
-                if(isNaturalRestriction){ 
-                  v = l; u2 = BU2; Dul = (*LGaussW[iGL])(iL);
-                }else{
-                  S.D12d (l, Puv, Vuv);
-                  Dul = Vuv.Y()*(*LGaussW[iGL])(iL);  // Dul = Du / Dl
-                  if(Abs(Dul) < EPS_PARAM) continue;
-                  v  = Puv.Y();  u2 = Puv.X();
-                  //Check on cause out off bounds of value current parameter
-                  if(v < BV1) v = BV1; else if(v > BV2) v = BV2;
-                  if(u2 < BU1) u2 = BU1; else if(u2 > BU2) u2 = BU2; 
-                }
-                ErrUL(iLS) = 0.0;
-                kUEnd = 1; JU = 0;
-                if(Abs(u2-u1) < EPS_PARAM) continue;
-                iUSubEnd = UFillIntervalBounds(u1, u2, UKnots, NumSubs);
-                UMaxSubs = MaxSubs(iUSubEnd);
-                //-- exception avoiding
-                if(UMaxSubs > SM) UMaxSubs = SM;
-                DimU.Init(0.0,1,UMaxSubs);  ErrU.Init(0.0,1,UMaxSubs);  ErrorU = 0.0;
-                do{//while: U
-                  if(++JU > iUSubEnd){
-                    URange[0] = IU = ErrU->Max();  URange[1] = JU;  
-                    U1(JU) = (U1(IU)+U2(IU))/2.0;  U2(JU) = U2(IU);  U2(IU) = U1(JU);
-                  }else  URange[0] = IU = JU;
-                  if(JU == UMaxSubs || Abs(U2(JU) - U1(JU)) < EPS_PARAM)
-                    if(kUEnd == 1){
-                      DimU(JU) = ErrU(JU) = IxU(JU) = IyU(JU) = IzU(JU) = 
-                        IxxU(JU) = IyyU(JU) = IzzU(JU) = IxyU(JU) = IxzU(JU) = IyzU(JU) = 0.0;
-                    }else{
-                      JU--;  
-                      EpsU = ErrorU;  Eps = EpsU*Abs((u2-u1)*Dul)/0.1;  EpsL = 0.9*Eps;
-                      break;
-                    }
-                  else
-                    for(kU=0; kU < kUEnd; kU++){
-                      iUS = URange[kU];
-                      um = 0.5*(U2(iUS) + U1(iUS));
-                      ur = 0.5*(U2(iUS) - U1(iUS));
-                      iGUEnd = iGLEnd - iGL;
-                      for(iGU=0; iGU < iGUEnd; iGU++){//
-                        LocDim[iGU] = 
-                          LocIxx[iGU] = LocIyy[iGU] = LocIzz[iGU] = 
-                          LocIx[iGU] = LocIy[iGU] = LocIz[iGU] = 
-                          LocIxy[iGU] = LocIxz[iGU] = LocIyz[iGU] = 0.0;
-                        for(iU=1; iU<=NbUGaussP[iGU]; iU++){
-                          u = um + ur*(*UGaussP[iGU])(iU);
-                          S.Normal(u, v, Ps, VNor);
-                          VNor.Coord(xn, yn, zn);
-                          Ps.Coord(x, y, z);  
-                          x -= xloc;  y -= yloc;  z -= zloc;
-                          xn *= (*UGaussW[iGU])(iU);
-                          yn *= (*UGaussW[iGU])(iU);
-                          zn *= (*UGaussW[iGU])(iU);
-                          if(ByPoint){
-                            //volume of elementary cone
-                            dDim = (x*xn+y*yn+z*zn)/3.0;
-                            //coordinates of cone's center mass 
-                            px = 0.75*x;  py = 0.75*y;  pz = 0.75*z;  
-                            LocDim[iGU] +=  dDim;
-                            //if(iGU > 0) continue;
-                            LocIx[iGU] += px*dDim;  
-                            LocIy[iGU] += py*dDim;  
-                            LocIz[iGU] += pz*dDim;
-                            x -= Coeff[0];  y -= Coeff[1];  z -= Coeff[2];
-                            dDim *= 3.0/5.0;
-                            LocIxy[iGU] -= x*y*dDim;  
-                            LocIyz[iGU] -= y*z*dDim;  
-                            LocIxz[iGU] -= x*z*dDim;
-                            xi = x*x;  yi = y*y;  zi = z*z;
-                            LocIxx[iGU] += (yi+zi)*dDim;  
-                            LocIyy[iGU] += (xi+zi)*dDim;  
-                            LocIzz[iGU] += (xi+yi)*dDim;
-                          }else{ // by plane
-                            s = xn*Coeff[0] + yn*Coeff[1] + zn*Coeff[2];
-                            d1 = Coeff[0]*x + Coeff[1]*y + Coeff[2]*z - Coeff[3];
-                            d2 = d1*d1;
-                            d3 = d1*d2/3.0;
-                            ds = s*d1;
-                            LocDim[iGU] += ds;
-                            //if(iGU > 0) continue;
-                            LocIx[iGU] += (x - Coeff[0]*d1/2.0) * ds;  
-                            LocIy[iGU] += (y - Coeff[1]*d1/2.0) * ds;
-                            LocIz[iGU] += (z - Coeff[2]*d1/2.0) * ds;
-                            px = x-Coeff[0]*d1;  py = y-Coeff[1]*d1;  pz = z-Coeff[2]*d1;
-                            xi = px*px*d1 + px*Coeff[0]*d2 + Coeff[0]*Coeff[0]*d3;
-                            yi = py*py*d1 + py*Coeff[1]*d2 + Coeff[1]*Coeff[1]*d3;
-                            zi = pz*pz*d1 + pz*Coeff[2]*d2 + Coeff[2]*Coeff[2]*d3;
-                            LocIxx[iGU] += (yi+zi)*s;  
-                            LocIyy[iGU] += (xi+zi)*s;  
-                            LocIzz[iGU] += (xi+yi)*s;
-                            d2 /= 2.0;
-                            xi = py*pz*d1 + py*Coeff[2]*d2 + pz*Coeff[1]*d2 + Coeff[1]*Coeff[2]*d3;
-                            yi = px*pz*d1 + pz*Coeff[0]*d2 + px*Coeff[2]*d2 + Coeff[0]*Coeff[2]*d3;
-                            zi = px*py*d1 + px*Coeff[1]*d2 + py*Coeff[0]*d2 + Coeff[0]*Coeff[1]*d3;
-                            LocIxy[iGU] -= zi*s;  LocIyz[iGU] -= xi*s;  LocIxz[iGU] -= yi*s;
-                          }
-                        }//for: iU
-                      }//for: iGU
-                      DimU(iUS) = LocDim[0]*ur;
-                      IxxU(iUS) = LocIxx[0]*ur; IyyU(iUS) = LocIyy[0]*ur; IzzU(iUS) = LocIzz[0]*ur;
-                      if(iGL > 0) continue;
-                      LocDim[1] = Abs(LocDim[1]-LocDim[0]); 
-                      LocIxx[1] = Abs(LocIxx[1]-LocIxx[0]); 
-                      LocIyy[1] = Abs(LocIyy[1]-LocIyy[0]); 
-                      LocIzz[1] = Abs(LocIzz[1]-LocIzz[0]);
-                      ErrU(iUS) = isMinDim? LocDim[1]*ur: (LocIxx[1] + LocIyy[1] + LocIzz[1])*ur;
-                      IxU(iUS) = LocIx[0]*ur; IyU(iUS) = LocIy[0]*ur; IzU(iUS) = LocIz[0]*ur;
-                      IxyU(iUS) = LocIxy[0]*ur; IxzU(iUS) = LocIxz[0]*ur; IyzU(iUS) = LocIyz[0]*ur;
-                    }//for: kU (iUS)
-                    if(JU == iUSubEnd)  kUEnd = 2; 
-                    if(kUEnd == 2)  {
-                      Standard_Integer imax = ErrU->Max();
-                      if(imax > 0) ErrorU = ErrU(imax);
-                      else ErrorU = 0.0;
-                    }
-                }while((ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1);
-                for(i=1; i<=JU; i++) {
-                  CDim[iGL] += DimU(i)*Dul;
-                  CIxx[iGL] += IxxU(i)*Dul; CIyy[iGL] += IyyU(i)*Dul; CIzz[iGL] += IzzU(i)*Dul;
-                }
-                if(iGL > 0) continue;
-                ErrUL(iLS) = ErrorU*Abs((u2-u1)*Dul);
-                for(i=1; i<=JU; i++){
-                  CIx += IxU(i)*Dul; CIy += IyU(i)*Dul; CIz += IzU(i)*Dul;
-                  //CIxx += IxxU(i)*Dul; CIyy += IyyU(i)*Dul; CIzz += IzzU(i)*Dul;
-                  CIxy += IxyU(i)*Dul; CIxz += IxzU(i)*Dul; CIyz += IyzU(i)*Dul;
-                }
-              }//for: iL 
-            }//for: iGL
-            DimL(iLS) = CDim[0]*lr;  
-            IxxL(iLS) = CIxx[0]*lr; IyyL(iLS) = CIyy[0]*lr; IzzL(iLS) = CIzz[0]*lr; 
-            if(iGLEnd == 2) {
-              //ErrL(iLS) = Abs(CDim[1]-CDim[0])*lr + ErrUL(iLS);
-              CDim[1] = Abs(CDim[1]-CDim[0]);
-              CIxx[1] = Abs(CIxx[1]-CIxx[0]); CIyy[1] = Abs(CIyy[1]-CIyy[0]); CIzz[1] = Abs(CIzz[1]-CIzz[0]);
-              ErrorU = ErrUL(iLS);
-              ErrL(iLS) = (isMinDim? CDim[1]: (CIxx[1] + CIyy[1] + CIzz[1]))*lr + ErrorU;
-            }
-            IxL(iLS) = CIx*lr; IyL(iLS) = CIy*lr; IzL(iLS) = CIz*lr; 
-            //IxxL(iLS) = CIxx*lr; IyyL(iLS) = CIyy*lr; IzzL(iLS) = CIzz*lr; 
-            IxyL(iLS) = CIxy*lr; IxzL(iLS) = CIxz*lr; IyzL(iLS) = CIyz*lr; 
-          }//for: (kL)iLS
-          //  Calculate/correct epsilon of computation by current value of Dim
-          //That is need for not spend time for 
-          if(JL == iLSubEnd){  
-            kLEnd = 2; 
-            Standard_Real DDim = 0.0, DIxx = 0.0, DIyy = 0.0, DIzz = 0.0;
-            for(i=1; i<=JL; i++) {
-              DDim += DimL(i);
-              DIxx += IxxL(i);  DIyy += IyyL(i);  DIzz += IzzL(i);
-            }
-            DDim = isMinDim? Abs(DDim): Abs(DIxx) + Abs(DIyy) + Abs(DIzz);
-            DDim = Abs(DDim*EpsDim);
-            if(DDim > Eps) { 
-              Eps = DDim;  EpsL = 0.9*Eps;
-            }
-          }
-          if(kLEnd == 2) {
-            Standard_Integer imax = ErrL->Max();
-            if(imax > 0) ErrorL = ErrL(imax);
-            else ErrorL = 0.0;
-          }
-      }while((ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1);
-      for(i=1; i<=JL; i++){
-        Dim += DimL(i); 
-        Ix += IxL(i); Iy += IyL(i); Iz += IzL(i);
-        Ixx += IxxL(i); Iyy += IyyL(i); Izz += IzzL(i);
-        Ixy += IxyL(i); Ixz += IxzL(i); Iyz += IyzL(i);
-      }
-      ErrorLMax = Max(ErrorLMax, ErrorL);
-    }
-    if(isNaturalRestriction) break;
-    D.Next();
-  }
-  if(Abs(Dim) >= EPS_DIM){
-    if(ByPoint){
-      Ix = Coeff[0] + Ix/Dim;
-      Iy = Coeff[1] + Iy/Dim;
-      Iz = Coeff[2] + Iz/Dim;
-    }else{
-      Ix /= Dim;
-      Iy /= Dim;
-      Iz /= Dim;
-    }
-    g.SetCoord (Ix, Iy, Iz);
-  }else{
-    Dim =0.;
-    g.SetCoord(0.,0.,0.);
-  }
-  inertia.SetCols (gp_XYZ (Ixx, Ixy, Ixz),
-    gp_XYZ (Ixy, Iyy, Iyz),
-    gp_XYZ (Ixz, Iyz, Izz));
-  if(iGLEnd == 2) 
-    Eps = Dim != 0.0? ErrorLMax/(isMinDim? Abs(Dim): (Abs(Ixx) + Abs(Iyy) + Abs(Izz))): 0.0;
-  else Eps = EpsDim;
-  return Eps;
+  SetLocation(theLocation);
+  Perform(theSurface, theDomain);
 }
 
-static Standard_Real Compute(BRepGProp_Face& S, const Standard_Boolean ByPoint, const Standard_Real Coeff[],
-                             const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, Standard_Real EpsDim) 
+//=======================================================================
+//function : BRepGProp_Vinert
+//purpose  : 
+//=======================================================================
+BRepGProp_Vinert::BRepGProp_Vinert(const BRepGProp_Face& theSurface,
+                                   const gp_Pnt&         theLocation)
 {
-  Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
-  Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
-  EpsDim = Abs(EpsDim);
-  BRepGProp_Domain D;
-  return CCompute(S,D,ByPoint,Coeff,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
+  SetLocation(theLocation);
+  Perform(theSurface);
+}
+
+//=======================================================================
+//function : BRepGProp_Vinert
+//purpose  : 
+//=======================================================================
+BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face&     theSurface,
+                                   const gp_Pnt&       theOrigin,
+                                   const gp_Pnt&       theLocation,
+                                   const Standard_Real theEps)
+{
+  SetLocation(theLocation);
+  Perform(theSurface, theOrigin, theEps);
+}
+
+//=======================================================================
+//function : BRepGProp_Vinert
+//purpose  : 
+//=======================================================================
+BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face&     theSurface,
+                                   BRepGProp_Domain&   theDomain,
+                                   const gp_Pnt&       theOrigin,
+                                   const gp_Pnt&       theLocation,
+                                   const Standard_Real theEps)
+{
+  SetLocation(theLocation);
+  Perform(theSurface, theDomain, theOrigin, theEps);
 }
 
-static Standard_Real Compute(BRepGProp_Face& S, BRepGProp_Domain& D, const Standard_Boolean ByPoint, const Standard_Real Coeff[],
-                             const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, Standard_Real EpsDim) 
+//=======================================================================
+//function : BRepGProp_Vinert
+//purpose  : 
+//=======================================================================
+BRepGProp_Vinert::BRepGProp_Vinert(const BRepGProp_Face& theSurface,
+                                   const gp_Pnt&         theOrigin,
+                                   const gp_Pnt&         theLocation)
+{
+  SetLocation(theLocation);
+  Perform(theSurface, theOrigin);
+}
+
+//=======================================================================
+//function : BRepGProp_Vinert
+//purpose  : 
+//=======================================================================
+BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face&   theSurface,
+                                   BRepGProp_Domain& theDomain,
+                                   const gp_Pnt&     theOrigin,
+                                   const gp_Pnt&     theLocation)
+{
+  SetLocation(theLocation);
+  Perform(theSurface, theDomain, theOrigin);
+}
+
+//=======================================================================
+//function : BRepGProp_Vinert
+//purpose  : 
+//=======================================================================
+BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face&     theSurface,
+                                   const gp_Pln&       thePlane,
+                                   const gp_Pnt&       theLocation,
+                                   const Standard_Real theEps)
 {
-  Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
-  Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
-  EpsDim = Abs(EpsDim);
-  return CCompute(S,D,ByPoint,Coeff,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
+  SetLocation(theLocation);
+  Perform(theSurface, thePlane, theEps);
+}
+
+//=======================================================================
+//function : BRepGProp_Vinert
+//purpose  : 
+//=======================================================================
+BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face&     theSurface,
+                                   BRepGProp_Domain&   theDomain,
+                                   const gp_Pln&       thePlane,
+                                   const gp_Pnt&       theLocation,
+                                   const Standard_Real theEps)
+{
+  SetLocation(theLocation);
+  Perform(theSurface, theDomain, thePlane, theEps);
 }
 
-static void Compute(const BRepGProp_Face& S,
-                    const Standard_Boolean ByPoint,
-                    const Standard_Real  Coeff[],
-                    const gp_Pnt& Loc,
-                    Standard_Real& Volu,
-                    gp_Pnt& G,
-                    gp_Mat& Inertia) 
+//=======================================================================
+//function : BRepGProp_Vinert
+//purpose  : 
+//=======================================================================
+BRepGProp_Vinert::BRepGProp_Vinert(const BRepGProp_Face& theSurface,
+                                   const gp_Pln&         thePlane,
+                                   const gp_Pnt&         theLocation)
 {
-
-  gp_Pnt P;            
-  gp_Vec VNor;     
-  Standard_Real dvi, dv;
-  Standard_Real ur, um, u, vr, vm, v;
-  Standard_Real x, y, z, xn, yn, zn, xi, yi, zi;
-  //  Standard_Real x, y, z, xn, yn, zn, xi, yi, zi, xyz;
-  Standard_Real px,py,pz,s,d1,d2,d3;
-  Standard_Real Ixi, Iyi, Izi, Ixxi, Iyyi, Izzi, Ixyi, Ixzi, Iyzi;
-  Standard_Real xloc, yloc, zloc;
-  Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
-
-  Volu = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
-  Loc.Coord (xloc, yloc, zloc);
-
-  Standard_Real LowerU, UpperU, LowerV, UpperV;
-  S.Bounds ( LowerU, UpperU, LowerV, UpperV);
-  Standard_Integer UOrder = Min(S.UIntegrationOrder (),
-    math::GaussPointsMax());
-  Standard_Integer VOrder = Min(S.VIntegrationOrder (),
-    math::GaussPointsMax());
-
-  Standard_Integer i, j;
-  math_Vector GaussPU (1, UOrder);     //gauss points and weights
-  math_Vector GaussWU (1, UOrder);
-  math_Vector GaussPV (1, VOrder);
-  math_Vector GaussWV (1, VOrder);
-
-  math::GaussPoints  (UOrder,GaussPU);
-  math::GaussWeights (UOrder,GaussWU);
-  math::GaussPoints  (VOrder,GaussPV);
-  math::GaussWeights (VOrder,GaussWV);
-
-  um = 0.5 * (UpperU + LowerU);
-  vm = 0.5 * (UpperV + LowerV);
-  ur = 0.5 * (UpperU - LowerU);
-  vr = 0.5 * (UpperV - LowerV);
-
-  for (j = 1; j <= VOrder; j++) {
-    v = vm + vr * GaussPV (j);
-    dvi = Ixi = Iyi = Izi = Ixxi = Iyyi = Izzi = Ixyi = Ixzi = Iyzi = 0.0;
-
-    for (i = 1; i <= UOrder; i++) {
-      u = um + ur * GaussPU (i);
-      S.Normal (u, v, P, VNor); 
-      VNor.Coord (xn, yn, zn);
-      P.Coord (x, y, z);
-      x -= xloc;  y -= yloc;  z -= zloc;
-      xn *= GaussWU (i);  yn *= GaussWU (i);  zn *= GaussWU (i);
-      if (ByPoint) {   
-        /////////////////////           ///////////////////////
-        //    OFV code     //           //    Initial code   //
-        /////////////////////           ///////////////////////
-        // modified by APO 
-        dv = (x*xn+y*yn+z*zn)/3.0;     //xyz  = x * y * z;
-        dvi += dv;                     //Ixyi += zn * xyz; 
-        Ixi += 0.75*x*dv;              //Iyzi += xn * xyz;
-        Iyi += 0.75*y*dv;              //Ixzi += yn * xyz; 
-        Izi += 0.75*z*dv;              //xi = x * x * x * xn / 3.0;
-        x -= Coeff[0];                 //yi = y * y * y * yn / 3.0;
-        y -= Coeff[1];                 //zi = z * z * z * zn / 3.0;
-        z -= Coeff[2];                 //Ixxi += (yi + zi);
-        dv *= 3.0/5.0;                 //Iyyi += (xi + zi);
-        Ixyi -= x*y*dv;                //Izzi += (xi + yi);
-        Iyzi -= y*z*dv;                //x -= Coeff[0];
-        Ixzi -= x*z*dv;                //y -= Coeff[1];
-        xi = x*x;                      //z -= Coeff[2];
-        yi = y*y;                      //dv = x * xn + y * yn + z * zn;
-        zi = z*z;                      //dvi +=  dv; 
-        Ixxi += (yi + zi)*dv;          //Ixi += x * dv;
-        Iyyi += (xi + zi)*dv;          //Iyi += y * dv;
-        Izzi += (xi + yi)*dv;          //Izi += z * dv;
-      }
-      else { // by plane
-        s   = xn * Coeff[0] + yn * Coeff[1] + zn * Coeff[2];
-        d1  = Coeff[0] * x + Coeff[1] * y + Coeff[2] * z - Coeff[3];
-        d2  = d1 * d1;
-        d3  = d1 * d2 / 3.0;
-        dv  = s * d1;
-        dvi += dv;
-        Ixi += (x - (Coeff[0] * d1 / 2.0)) * dv;
-        Iyi += (y - (Coeff[1] * d1 / 2.0)) * dv;
-        Izi += (z - (Coeff[2] * d1 / 2.0)) * dv;
-        px  = x - Coeff[0] * d1;
-        py  = y - Coeff[1] * d1;
-        pz  = z - Coeff[2] * d1;
-        xi  = px * px * d1 + px * Coeff[0]* d2 + Coeff[0] * Coeff[0] * d3;
-        yi  = py * py * d1 + py * Coeff[1] * d2 + Coeff[1] * Coeff[1] * d3;
-        zi  = pz * pz * d1 + pz * Coeff[2] * d2 + Coeff[2] * Coeff[2] * d3;
-        Ixxi += (yi + zi) * s;
-        Iyyi += (xi + zi) * s;
-        Izzi += (xi + yi) * s;
-        d2  /= 2.0;
-        xi  = (py * pz * d1) + (py * Coeff[2] * d2) + (pz * Coeff[1] * d2) + (Coeff[1] * Coeff[2] * d3);
-        yi  = (px * pz * d1) + (pz * Coeff[0] * d2) + (px * Coeff[2] * d2) + (Coeff[0] * Coeff[2] * d3);
-        zi  = (px * py * d1) + (px * Coeff[1] * d2) + (py * Coeff[0] * d2) + (Coeff[0] * Coeff[1] * d3);
-        Ixyi -=  zi * s;
-        Iyzi -=  xi * s;
-        Ixzi -=  yi * s;
-      }
-    }
-    Volu += dvi  * GaussWV (j);
-    Ix   += Ixi  * GaussWV (j);
-    Iy   += Iyi  * GaussWV (j);
-    Iz   += Izi  * GaussWV (j);
-    Ixx  += Ixxi * GaussWV (j);
-    Iyy  += Iyyi * GaussWV (j);
-    Izz  += Izzi * GaussWV (j);
-    Ixy  += Ixyi * GaussWV (j);
-    Ixz  += Ixzi * GaussWV (j);
-    Iyz  += Iyzi * GaussWV (j);
-  }
-  vr  *= ur;
-  Ixx *= vr;
-  Iyy *= vr;
-  Izz *= vr;
-  Ixy *= vr;
-  Ixz *= vr;
-  Iyz *= vr;
-  if (Abs(Volu) >= EPS_DIM ) {
-    if (ByPoint) {
-      Ix  = Coeff[0] + Ix/Volu;
-      Iy  = Coeff[1] + Iy/Volu;
-      Iz  = Coeff[2] + Iz/Volu;
-      Volu *= vr;
-    }
-    else { //by plane
-      Ix /= Volu;
-      Iy /= Volu;
-      Iz /= Volu;
-      Volu *= vr;
-    }
-    G.SetCoord (Ix, Iy, Iz);
-  } 
-  else {
-    G.SetCoord(0.,0.,0.);
-    Volu =0.;
-  }
-  Inertia.SetCols (gp_XYZ (Ixx, Ixy, Ixz),
-    gp_XYZ (Ixy, Iyy, Iyz),
-    gp_XYZ (Ixz, Iyz, Izz));
-
+  SetLocation(theLocation);
+  Perform(theSurface, thePlane);
+}
+
+//=======================================================================
+//function : BRepGProp_Vinert
+//purpose  : 
+//=======================================================================
+BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face&   theSurface,
+                                   BRepGProp_Domain& theDomain,
+                                   const gp_Pln&     thePlane,
+                                   const gp_Pnt&     theLocation)
+{
+  SetLocation(theLocation);
+  Perform(theSurface, theDomain, thePlane);
 }
 
-// Last modified by OFV 5.2001:
-// 1). surface and edge integration order is equal now
-// 2). "by point" works now rathre correctly (it looks so...)
-static void Compute(BRepGProp_Face& S, BRepGProp_Domain& D, const Standard_Boolean ByPoint, const Standard_Real  Coeff[],
-                    const gp_Pnt& Loc, Standard_Real& Volu, gp_Pnt& G, gp_Mat& Inertia) 
-
+//=======================================================================
+//function : SetLocation
+//purpose  : 
+//=======================================================================
+void BRepGProp_Vinert::SetLocation(const gp_Pnt& theLocation)
 {
-  Standard_Real x, y, z, xi, yi, zi, l1, l2, lm, lr, l, v1, v2, v, u1, u2, um, ur, u, ds, Dul, xloc, yloc, zloc;
-  Standard_Real LocVolu, LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy, LocIxz, LocIyz;
-  Standard_Real CVolu, CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz, Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
-  Standard_Real xn, yn, zn, px, py, pz, s, d1, d2, d3, dSigma;
-  Standard_Integer i, j, vio, sio, max, NbGaussgp_Pnts;
-
-  gp_Pnt Ps;
-  gp_Vec VNor; 
-  gp_Pnt2d Puv;
-  gp_Vec2d Vuv;
-
-  Loc.Coord (xloc, yloc, zloc);
-  Volu = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
-  S.Bounds (u1, u2, v1, v2);
-  Standard_Real _u2 = u2;  //OCC104
-  vio = S.VIntegrationOrder ();
-
-  while (D.More())
-  {
-    S.Load(D.Value());
-    sio = S.IntegrationOrder ();
-    max = Max(vio,sio);
-    NbGaussgp_Pnts = Min(max,math::GaussPointsMax());
-
-    math_Vector GaussP (1, NbGaussgp_Pnts);
-    math_Vector GaussW (1, NbGaussgp_Pnts);
-    math::GaussPoints  (NbGaussgp_Pnts,GaussP);
-    math::GaussWeights (NbGaussgp_Pnts,GaussW);
-
-    CVolu = CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
-    l1  = S.FirstParameter();
-    l2  = S.LastParameter();
-    lm  = 0.5 * (l2 + l1);         
-    lr  = 0.5 * (l2 - l1);
-
-    for (i=1; i<=NbGaussgp_Pnts; i++)
-    {
-      l = lm + lr * GaussP(i);
-      S.D12d (l, Puv, Vuv);
-      v   = Puv.Y();
-      u2  = Puv.X();
-
-      //OCC104
-      v = v < v1? v1: v; 
-      v = v > v2? v2: v; 
-      u2 = u2 < u1? u1: u2; 
-      u2 = u2 > _u2? _u2: u2; 
-
-      Dul = Vuv.Y() * GaussW(i);
-      um  = 0.5 * (u2 + u1);
-      ur  = 0.5 * (u2 - u1);
-      LocVolu = LocIx = LocIy = LocIz = LocIxx = LocIyy = LocIzz = LocIxy = LocIxz = LocIyz = 0.0;
-
-      for (j=1; j<=NbGaussgp_Pnts; j++)
-      {
-        u = um + ur * GaussP(j);
-        S.Normal (u, v, Ps, VNor);
-        VNor.Coord (xn, yn, zn);
-        Ps.Coord (x, y, z);
-        x      -= xloc; 
-        y      -= yloc;
-        z      -= zloc;
-        xn = xn * Dul * GaussW(j);
-        yn = yn * Dul * GaussW(j);
-        zn = zn * Dul * GaussW(j);
-        if(ByPoint)
-        {
-          dSigma = (x*xn+y*yn+z*zn)/3.0;
-          LocVolu +=  dSigma; 
-          LocIx   += 0.75*x*dSigma;  
-          LocIy   += 0.75*y*dSigma;
-          LocIz   += 0.75*z*dSigma;
-          x      -= Coeff[0];
-          y      -= Coeff[1];
-          z      -= Coeff[2];
-          dSigma *= 3.0/5.0; 
-          LocIxy -= x*y*dSigma;
-          LocIyz -= y*z*dSigma;
-          LocIxz -= x*z*dSigma;
-          xi      = x*x;
-          yi      = y*y;
-          zi      = z*z;
-          LocIxx += (yi + zi)*dSigma;
-          LocIyy += (xi + zi)*dSigma;
-          LocIzz += (xi + yi)*dSigma;
-        }
-        else
-        {
-          s   = xn * Coeff[0] + yn * Coeff[1] + zn * Coeff[2];
-          d1  = Coeff[0] * x + Coeff[1] * y + Coeff[2] * z;
-          d2  = d1 * d1;
-          d3  = d1 * d2 / 3.0;
-          ds  = s * d1;
-          LocVolu += ds;
-          LocIx += (x - Coeff[0] * d1 / 2.0) * ds;
-          LocIy += (y - Coeff[1] * d1 / 2.0) * ds;
-          LocIz += (z - Coeff[2] * d1 / 2.0) * ds;
-          px  = x - Coeff[0] * d1;
-          py  = y - Coeff[1] * d1;
-          pz  = z - Coeff[2] * d1;
-          xi  = (px * px * d1) + (px * Coeff[0]* d2) + (Coeff[0] * Coeff[0] * d3);
-          yi  = (py * py * d1) + (py * Coeff[1] * d2) + (Coeff[1] * Coeff[1] * d3);
-          zi  = pz * pz * d1 + pz * Coeff[2] * d2 + (Coeff[2] * Coeff[2] * d3);
-          LocIxx += (yi + zi) * s;
-          LocIyy += (xi + zi) * s;
-          LocIzz += (xi + yi) * s;
-          d2  /= 2.0;
-          xi  = (py * pz * d1) + (py * Coeff[2] * d2) + (pz * Coeff[1] * d2) + (Coeff[1] * Coeff[2] * d3);
-          yi  = (px * pz * d1) + (pz * Coeff[0] * d2) + (px * Coeff[2] * d2) + (Coeff[0] * Coeff[2] * d3);
-          zi  = (px * py * d1) + (px * Coeff[1] * d2) + (py * Coeff[0] * d2) + (Coeff[0] * Coeff[1] * d3);
-          LocIxy -=  zi * s;
-          LocIyz -=  xi * s;
-          LocIxz -=  yi * s;
-        }
-      }
-      CVolu += LocVolu * ur;
-      CIx   += LocIx   * ur;
-      CIy   += LocIy   * ur;
-      CIz   += LocIz   * ur;
-      CIxx  += LocIxx  * ur;
-      CIyy  += LocIyy  * ur;
-      CIzz  += LocIzz  * ur;
-      CIxy  += LocIxy  * ur;
-      CIxz  += LocIxz  * ur;
-      CIyz  += LocIyz  * ur;
-    }
-    Volu   += CVolu * lr;
-    Ix     += CIx   * lr;
-    Iy     += CIy   * lr;
-    Iz     += CIz   * lr;
-    Ixx    += CIxx  * lr;
-    Iyy    += CIyy  * lr;
-    Izz    += CIzz  * lr;
-    Ixy    += CIxy  * lr;
-    Ixz    += CIxz  * lr;
-    Iyz    += CIyz  * lr;
-    D.Next();
-  }
-
-  if(Abs(Volu) >= EPS_DIM)
-  {
-    if(ByPoint)
-    {
-      Ix = Coeff[0] + Ix/Volu;
-      Iy = Coeff[1] + Iy/Volu;
-      Iz = Coeff[2] + Iz/Volu;
-    }
-    else
-    {
-      Ix /= Volu;
-      Iy /= Volu;
-      Iz /= Volu;
-    }
-    G.SetCoord (Ix, Iy, Iz);
-  }
-  else
-  {
-    Volu =0.;
-    G.SetCoord(0.,0.,0.);
-  }
-
-  Inertia.SetCols (gp_XYZ (Ixx, Ixy, Ixz),
-    gp_XYZ (Ixy, Iyy, Iyz),
-    gp_XYZ (Ixz, Iyz, Izz));
-
+  loc = theLocation;
 }
 
-BRepGProp_Vinert::BRepGProp_Vinert(){}
-
-BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face& S, const gp_Pnt& VLocation, const Standard_Real Eps){
-  SetLocation(VLocation);
-  Perform(S,Eps);
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+Standard_Real BRepGProp_Vinert::Perform(BRepGProp_Face&     theSurface,
+                                        const Standard_Real theEps)
+{
+  BRepGProp_Domain anEmptyDomain;
+  return Perform(theSurface, anEmptyDomain, theEps);
 }
 
-BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pnt& VLocation, const Standard_Real Eps){
-  SetLocation(VLocation);
-  Perform(S,D,Eps);
-}
+//=======================================================================
+//function : 
+//purpose  : 
+//=======================================================================
+Standard_Real BRepGProp_Vinert::Perform(BRepGProp_Face&     theSurface,
+                                        BRepGProp_Domain&   theDomain,
+                                        const Standard_Real theEps)
+{
+  const Standard_Real aCoeff[] = {0.0, 0.0, 0.0};
+  BRepGProp_Gauss     aGauss(BRepGProp_Gauss::Vinert);
 
-BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pnt& VLocation){
-  SetLocation(VLocation);
-  Perform(S,D);
+  return myEpsilon =
+    aGauss.Compute(theSurface, theDomain, loc, theEps,
+                   aCoeff, Standard_True, dim, g, inertia);
 }
 
-BRepGProp_Vinert::BRepGProp_Vinert(const BRepGProp_Face& S, const gp_Pnt& VLocation){
-  SetLocation(VLocation);
-  Perform(S);
-}
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void BRepGProp_Vinert::Perform(const BRepGProp_Face& theSurface)
+{
+  const Standard_Real aCoeff[] = {0.0, 0.0, 0.0};
+  BRepGProp_Gauss     aGauss(BRepGProp_Gauss::Vinert);
 
-BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face& S, const gp_Pnt& O, const gp_Pnt& VLocation, const Standard_Real Eps){
-  SetLocation(VLocation);
-  Perform(S,O,Eps);
+  myEpsilon = 1.0;
+  aGauss.Compute(theSurface, loc, aCoeff, Standard_True, dim, g, inertia);
 }
 
-BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pnt& O, const gp_Pnt& VLocation, const Standard_Real Eps){
-  SetLocation(VLocation);
-  Perform(S,D,O,Eps);
-}
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void BRepGProp_Vinert::Perform(BRepGProp_Face&   theSurface,
+                               BRepGProp_Domain& theDomain)
+{
+  const Standard_Real aCoeff[] = {0.0, 0.0, 0.0};
+  BRepGProp_Gauss     aGauss(BRepGProp_Gauss::Vinert);
 
-BRepGProp_Vinert::BRepGProp_Vinert(const BRepGProp_Face& S, const gp_Pnt& O, const gp_Pnt& VLocation){
-  SetLocation(VLocation);
-  Perform(S,O);
+  myEpsilon = 1.0;
+  aGauss.Compute(theSurface, theDomain, loc, aCoeff, Standard_True, dim, g, inertia);
 }
 
-BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pnt& O, const gp_Pnt& VLocation){
-  SetLocation(VLocation);
-  Perform(S,D,O);
-}
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+Standard_Real BRepGProp_Vinert::Perform(BRepGProp_Face&     theSurface,
+                                        const gp_Pnt&       theOrigin,
+                                        const Standard_Real theEps)
+{
+  BRepGProp_Domain anEmptyDomain;
+  return Perform(theSurface, anEmptyDomain, theOrigin, theEps);
+}
+
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+Standard_Real BRepGProp_Vinert::Perform(BRepGProp_Face&     theSurface,
+                                        BRepGProp_Domain&   theDomain,
+                                        const gp_Pnt&       theOrigin,
+                                        const Standard_Real theEps)
+{
+  const Standard_Real aCoeff[] =
+  {
+    theOrigin.X() - loc.X(),
+    theOrigin.Y() - loc.Y(),
+    theOrigin.Z() - loc.Z()
+  };
 
-BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face& S, const gp_Pln& Pl, const gp_Pnt& VLocation, const Standard_Real Eps){
-  SetLocation(VLocation);
-  Perform(S,Pl,Eps);
-}
+  BRepGProp_Gauss  aGauss(BRepGProp_Gauss::Vinert);
 
-BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pln& Pl, const gp_Pnt& VLocation, const Standard_Real Eps){
-  SetLocation(VLocation);
-  Perform(S,D,Pl,Eps);
+  return myEpsilon =
+    aGauss.Compute(theSurface, theDomain, loc, theEps,
+                   aCoeff, Standard_True, dim, g, inertia);
 }
 
-BRepGProp_Vinert::BRepGProp_Vinert(const BRepGProp_Face& S, const gp_Pln& Pl, const gp_Pnt& VLocation){
-  SetLocation(VLocation);
-  Perform(S,Pl);
-}
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void BRepGProp_Vinert::Perform(const BRepGProp_Face& theSurface,
+                               const gp_Pnt&         theOrigin)
+{
+  BRepGProp_Gauss     aGauss(BRepGProp_Gauss::Vinert);
+  const Standard_Real aCoeff[] =
+  {
+    theOrigin.X() - loc.X(),
+    theOrigin.Y() - loc.Y(),
+    theOrigin.Z() - loc.Z()
+  };
 
-BRepGProp_Vinert::BRepGProp_Vinert(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pln& Pl, const gp_Pnt& VLocation){
-  SetLocation(VLocation);
-  Perform(S,D,Pl);
+  myEpsilon = 1.0;
+  aGauss.Compute(theSurface, loc, aCoeff, Standard_True, dim, g, inertia);
 }
 
-void BRepGProp_Vinert::SetLocation(const gp_Pnt& VLocation){
-  loc = VLocation;
-}
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void BRepGProp_Vinert::Perform(BRepGProp_Face&   theSurface,
+                               BRepGProp_Domain& theDomain,
+                               const gp_Pnt&     theOrigin)
+{
+  BRepGProp_Gauss     aGauss(BRepGProp_Gauss::Vinert);
+  const Standard_Real aCoeff[] =
+  {
+    theOrigin.X() - loc.X(),
+    theOrigin.Y() - loc.Y(),
+    theOrigin.Z() - loc.Z()
+  };
 
-Standard_Real BRepGProp_Vinert::Perform(BRepGProp_Face& S, const Standard_Real Eps){
-  Standard_Real Coeff[] = {0., 0., 0.};
-  return myEpsilon = Compute(S,Standard_True,Coeff,loc,dim,g,inertia,Eps);
+  myEpsilon = 1.0;
+  aGauss.Compute(theSurface, theDomain, loc, aCoeff, Standard_True, dim, g, inertia);
 }
 
-Standard_Real BRepGProp_Vinert::Perform(BRepGProp_Face& S, BRepGProp_Domain& D, const Standard_Real Eps){
-  Standard_Real Coeff[] = {0., 0., 0.};
-  return myEpsilon = Compute(S,D,Standard_True,Coeff,loc,dim,g,inertia,Eps);
-}
+//=======================================================================
+//function : 
+//purpose  : 
+//=======================================================================
+Standard_Real BRepGProp_Vinert::Perform(BRepGProp_Face&     theSurface,
+                                        const gp_Pln&       thePlane,
+                                        const Standard_Real theEps)
+{
+  BRepGProp_Domain anEmptyDomain;
+  return Perform(theSurface, anEmptyDomain, thePlane, theEps);
+}
+
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+Standard_Real BRepGProp_Vinert::Perform(BRepGProp_Face&     theSurface,
+                                        BRepGProp_Domain&   theDomain,
+                                        const gp_Pln&       thePlane,
+                                        const Standard_Real theEps)
+{
+  Standard_Real aCoeff[4];
+  thePlane.Coefficients(aCoeff[0], aCoeff[1], aCoeff[2], aCoeff[3]);
+  aCoeff[3] = aCoeff[3] - aCoeff[0] * loc.X()
+                        - aCoeff[1] * loc.Y()
+                        - aCoeff[2] * loc.Z();
 
-void BRepGProp_Vinert::Perform(const BRepGProp_Face& S){
-  Standard_Real Coeff[] = {0., 0., 0.};
-  Compute(S,Standard_True,Coeff,loc,dim,g,inertia);
-  myEpsilon = 1.0;
-  return;
-}
+  BRepGProp_Gauss  aGauss(BRepGProp_Gauss::Vinert);
 
-void BRepGProp_Vinert::Perform(BRepGProp_Face& S, BRepGProp_Domain& D){
-  Standard_Real Coeff[] = {0., 0., 0.};
-  Compute(S,D,Standard_True,Coeff,loc,dim,g,inertia);
-  myEpsilon = 1.0;
-  return;
+  return myEpsilon =
+    aGauss.Compute(theSurface, theDomain, loc, theEps,
+                   aCoeff, Standard_False, dim, g, inertia);
 }
 
-Standard_Real BRepGProp_Vinert::Perform(BRepGProp_Face& S, const gp_Pnt&  O, const Standard_Real Eps){
-  Standard_Real xloc, yloc, zloc;
-  loc.Coord(xloc, yloc, zloc);
-  Standard_Real Coeff[3];
-  O.Coord (Coeff[0], Coeff[1], Coeff[2]);
-  Coeff[0] -= xloc;  Coeff[1] -= yloc;   Coeff[2] -= zloc;
-  return myEpsilon = Compute(S,Standard_True,Coeff,loc,dim,g,inertia,Eps);
-}
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void BRepGProp_Vinert::Perform(const BRepGProp_Face& theSurface,
+                               const gp_Pln&         thePlane)
+{
+  BRepGProp_Gauss aGauss(BRepGProp_Gauss::Vinert);
+  Standard_Real   aCoeff[4];
 
-Standard_Real BRepGProp_Vinert::Perform(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pnt&  O, const Standard_Real Eps){
-  Standard_Real xloc, yloc, zloc;
-  loc.Coord(xloc, yloc, zloc);
-  Standard_Real Coeff[3];
-  O.Coord (Coeff[0], Coeff[1], Coeff[2]);
-  Coeff[0] -= xloc;  Coeff[1] -= yloc;   Coeff[2] -= zloc;
-  return myEpsilon = Compute(S,D,Standard_True,Coeff,loc,dim,g,inertia,Eps);
-}
+  thePlane.Coefficients  (aCoeff[0],
+                          aCoeff[1],
+                          aCoeff[2],
+                          aCoeff[3]);
 
-void BRepGProp_Vinert::Perform(const BRepGProp_Face& S, const gp_Pnt&  O){
-  Standard_Real xloc, yloc, zloc;
-  loc.Coord(xloc, yloc, zloc);
-  Standard_Real Coeff[3];
-  O.Coord (Coeff[0], Coeff[1], Coeff[2]);
-  Coeff[0] -= xloc;  Coeff[1] -= yloc;   Coeff[2] -= zloc;
-  Compute(S,Standard_True,Coeff,loc,dim,g,inertia);
-  myEpsilon = 1.0;
-  return;
-}
+  aCoeff[3] = aCoeff[3] - aCoeff[0] * loc.X()
+                        - aCoeff[1] * loc.Y()
+                        - aCoeff[2] * loc.Z();
 
-void BRepGProp_Vinert::Perform(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pnt&  O){
-  Standard_Real xloc, yloc, zloc;
-  loc.Coord(xloc, yloc, zloc);
-  Standard_Real Coeff[3];
-  O.Coord (Coeff[0], Coeff[1], Coeff[2]);
-  Coeff[0] -= xloc;  Coeff[1] -= yloc;   Coeff[2] -= zloc;
-  Compute(S,D,Standard_True,Coeff,loc,dim,g,inertia);
   myEpsilon = 1.0;
-  return;
+  aGauss.Compute(theSurface, loc, aCoeff, Standard_False, dim, g, inertia);
 }
 
-Standard_Real BRepGProp_Vinert::Perform(BRepGProp_Face& S, const gp_Pln& Pl, const Standard_Real Eps){
-  Standard_Real xloc, yloc, zloc;
-  loc.Coord (xloc, yloc, zloc);
-  Standard_Real Coeff[4];
-  Pl.Coefficients (Coeff[0], Coeff[1],Coeff[2],Coeff[3]); 
-  Coeff[3] = Coeff[3] - Coeff[0]*xloc - Coeff[1]*yloc - Coeff[2]*zloc;
-  return myEpsilon = Compute(S,Standard_False,Coeff,loc,dim,g,inertia,Eps);
-}
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void BRepGProp_Vinert::Perform(BRepGProp_Face&   theSurface,
+                               BRepGProp_Domain& theDomain,
+                               const gp_Pln&     thePlane)
+{
+  BRepGProp_Gauss aGauss(BRepGProp_Gauss::Vinert);
+  Standard_Real   aCoeff[4];
 
-Standard_Real BRepGProp_Vinert::Perform(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pln& Pl, const Standard_Real Eps){
-  Standard_Real xloc, yloc, zloc;
-  loc.Coord (xloc, yloc, zloc);
-  Standard_Real Coeff[4];
-  Pl.Coefficients (Coeff[0], Coeff[1],Coeff[2],Coeff[3]); 
-  Coeff[3] = Coeff[3] - Coeff[0]*xloc - Coeff[1]*yloc - Coeff[2]*zloc;
-  return myEpsilon = Compute(S,D,Standard_False,Coeff,loc,dim,g,inertia,Eps);
-}
+  thePlane.Coefficients  (aCoeff[0],
+                          aCoeff[1],
+                          aCoeff[2],
+                          aCoeff[3]);
 
-void BRepGProp_Vinert::Perform(const BRepGProp_Face& S, const gp_Pln& Pl){
-  Standard_Real xloc, yloc, zloc;
-  loc.Coord (xloc, yloc, zloc);
-  Standard_Real Coeff[4];
-  Pl.Coefficients (Coeff[0], Coeff[1],Coeff[2],Coeff[3]); 
-  Coeff[3] = Coeff[3] - Coeff[0]*xloc - Coeff[1]*yloc - Coeff[2]*zloc;
-  Compute(S,Standard_False,Coeff,loc,dim,g,inertia);
-  myEpsilon = 1.0;
-  return;
-}
+  aCoeff[3] = aCoeff[3] - aCoeff[0] * loc.X()
+                        - aCoeff[1] * loc.Y()
+                        - aCoeff[2] * loc.Z();
 
-void BRepGProp_Vinert::Perform(BRepGProp_Face& S, BRepGProp_Domain& D, const gp_Pln& Pl){
-  Standard_Real xloc, yloc, zloc;
-  loc.Coord (xloc, yloc, zloc);
-  Standard_Real Coeff[4];
-  Pl.Coefficients (Coeff[0], Coeff[1],Coeff[2],Coeff[3]); 
-  Coeff[3] = Coeff[3] - Coeff[0]*xloc - Coeff[1]*yloc - Coeff[2]*zloc;
-  Compute(S,D,Standard_False,Coeff,loc,dim,g,inertia);
   myEpsilon = 1.0;
-  return;
+  aGauss.Compute(theSurface, theDomain, loc, aCoeff, Standard_False, dim, g, inertia);
 }
 
-Standard_Real BRepGProp_Vinert::GetEpsilon(){
+//=======================================================================
+//function : GetEpsilon
+//purpose  : 
+//=======================================================================
+Standard_Real BRepGProp_Vinert::GetEpsilon()
+{
   return myEpsilon;
 }
diff --git a/src/BRepGProp/FILES b/src/BRepGProp/FILES
new file mode 100644 (file)
index 0000000..b47e4d6
--- /dev/null
@@ -0,0 +1,2 @@
+BRepGProp_Gauss.hxx
+BRepGProp_Gauss.cxx
\ No newline at end of file