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.
--- /dev/null
+// 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);
+}
--- /dev/null
+// 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
// 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;
}
// 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;
}
--- /dev/null
+BRepGProp_Gauss.hxx
+BRepGProp_Gauss.cxx
\ No newline at end of file