mirror of
https://git.dev.opencascade.org/repos/occt.git
synced 2025-04-10 18:51:21 +03:00
9bd59d1c5b
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.
1381 lines
46 KiB
C++
1381 lines
46 KiB
C++
// 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);
|
|
}
|