diff --git a/src/ModelingData/TKGeomBase/GProp/GProp_SelGProps.cxx b/src/ModelingData/TKGeomBase/GProp/GProp_SelGProps.cxx index 56ba2164152..4d20585582d 100644 --- a/src/ModelingData/TKGeomBase/GProp/GProp_SelGProps.cxx +++ b/src/ModelingData/TKGeomBase/GProp/GProp_SelGProps.cxx @@ -122,9 +122,10 @@ void GProp_SelGProps::Perform(const gp_Cone& S, double Auxi1 = R + (Z2 + Z1) * Snt / 2.; double Auxi2 = (Z2 * Z2 + Z1 * Z2 + Z1 * Z1) / 3.; - dim = (Alpha2 - Alpha1) * Cnt * (Z2 - Z1) * Auxi1; + // The area element of gp_Cone is R + v sin(a) dv du. + dim = (Alpha2 - Alpha1) * (Z2 - Z1) * Auxi1; - double Ix = (R * R + R * (Z2 + Z1) * Snt + Snt * Auxi2) / Auxi1; + double Ix = (R * R + R * (Z2 + Z1) * Snt + Snt * Snt * Auxi2) / Auxi1; double Iy = Ix * (Cn1 - Cn2) / (Alpha2 - Alpha1); Ix = Ix * (Sn2 - Sn1) / (Alpha2 - Alpha1); double Iz = Cnt * (R * (Z2 + Z1) / 2. + Snt * Auxi2) / Auxi1; @@ -133,18 +134,21 @@ void GProp_SelGProps::Perform(const gp_Cone& S, Y0 + Ya1 * Ix + Ya2 * Iy + Ya3 * Iz, Z0 + Za1 * Ix + Za2 * Iy + Za3 * Iz); - double R1 = R + Z1 * Snt; - double R2 = R + Z2 * Snt; - double ZZ = (Z2 - Z1) * Cnt; - double IR2 = ZZ * Snt * (R1 * R1 * R1 + R1 * R1 * R2 + R1 * R2 * R2 + R2 * R2 * R2) / 4.; - double ICn2 = IR2 * (Alpha2 - Alpha1 + Cn2 * Sn2 - Cn1 * Sn1) / 2.; - double ISn2 = IR2 * (Alpha2 - Alpha1 + Cn2 * Sn2 - Cn1 * Sn1) / 2.; - double IZ2 = ZZ * Cnt * Cnt * (Z2 - Z1) * (Alpha2 - Alpha1) - * (R * Auxi2 + Snt * (Z2 * Z2 * Z2 + Z2 * Z2 * Z1 + Z2 * Z1 * Z1 + Z1 * Z1 * Z1)) - / 4.; - double ICnSn = IR2 * (Cn2 * Cn2 - Cn1 * Cn1); - double ICnz = Cnt * Snt * ZZ * (R * (Z1 + Z2) / 2. + Auxi2) * (Sn2 - Sn1); - double ISnz = Cnt * Snt * ZZ * (R * (Z1 + Z2) / 2. + Auxi2) * (Cn1 - Cn2); + double R1 = R + Z1 * Snt; + double R2 = R + Z2 * Snt; + double Z3 = Z2 * Z2 * Z2 + Z2 * Z2 * Z1 + Z2 * Z1 * Z1 + Z1 * Z1 * Z1; + double Z4 = (Z1 + Z2) * (Z1 * Z1 + Z2 * Z2); + // Per unit angle, r = R + v sin(a), z = v cos(a), dA = r dv: IR2 = int r^3, IZ2 = int z^2 r, + // IRZ = int r^2 z. + double IR2 = (Z2 - Z1) * (R1 * R1 * R1 + R1 * R1 * R2 + R1 * R2 * R2 + R2 * R2 * R2) / 4.; + double IRZ = + Cnt * (Z2 - Z1) * (R * R * (Z1 + Z2) / 2. + 2. * R * Snt * Auxi2 + Snt * Snt * Z4 / 4.); + double ICn2 = IR2 * (Alpha2 - Alpha1 + Cn2 * Sn2 - Cn1 * Sn1) / 2.; + double ISn2 = IR2 * (Alpha2 - Alpha1 - Cn2 * Sn2 + Cn1 * Sn1) / 2.; + double IZ2 = Cnt * Cnt * (Z2 - Z1) * (Alpha2 - Alpha1) * (R * Auxi2 + Snt * Z3 / 4.); + double ICnSn = IR2 * (Sn2 * Sn2 - Sn1 * Sn1) / 2.; + double ICnz = IRZ * (Sn2 - Sn1); + double ISnz = IRZ * (Cn1 - Cn2); math_Matrix Dm(1, 3, 1, 3); Dm(1, 1) = ISn2 + IZ2; @@ -154,32 +158,13 @@ void GProp_SelGProps::Perform(const gp_Cone& S, Dm(1, 3) = Dm(3, 1) = -ICnz; Dm(3, 2) = Dm(2, 3) = -ISnz; - math_Matrix Passage(1, 3, 1, 3); - Passage(1, 1) = Xa1; - Passage(1, 2) = Xa2; - Passage(1, 3) = Xa3; - Passage(2, 1) = Ya1; - Passage(2, 2) = Ya2; - Passage(2, 3) = Ya3; - Passage(3, 1) = Za1; - Passage(3, 2) = Za2; - Passage(3, 3) = Za3; - - math_Jacobi J(Dm); - math_Vector V1(1, 3), V2(1, 3), V3(1, 3); - J.Vector(1, V1); - V1.Multiply(Passage, V1); - V1.Multiply(J.Value(1)); - J.Vector(2, V2); - V2.Multiply(Passage, V2); - V2.Multiply(J.Value(2)); - J.Vector(3, V3); - V3.Multiply(Passage, V3); - V3.Multiply(J.Value(3)); - - inertia = - gp_Mat(gp_XYZ(V1(1), V2(1), V3(1)), gp_XYZ(V1(2), V2(2), V3(2)), gp_XYZ(V1(3), V2(3), V3(3))); + // Dm is about the cone location in the cone axes: take it to global axes, to g, then to loc. + gp_Mat Pm(Xa1, Xa2, Xa3, Ya1, Ya2, Ya3, Za1, Za2, Za3); + gp_Mat + Dg(Dm(1, 1), Dm(1, 2), Dm(1, 3), Dm(2, 1), Dm(2, 2), Dm(2, 3), Dm(3, 1), Dm(3, 2), Dm(3, 3)); gp_Mat Hop; + GProp::HOperator(g, S.Location(), dim, Hop); + inertia = Pm * Dg * Pm.Transposed() - Hop; GProp::HOperator(g, loc, dim, Hop); inertia = inertia + Hop; } diff --git a/src/ModelingData/TKGeomBase/GProp/GProp_VelGProps.cxx b/src/ModelingData/TKGeomBase/GProp/GProp_VelGProps.cxx index 67b3a5eab0a..755adec52eb 100644 --- a/src/ModelingData/TKGeomBase/GProp/GProp_VelGProps.cxx +++ b/src/ModelingData/TKGeomBase/GProp/GProp_VelGProps.cxx @@ -165,31 +165,46 @@ void GProp_VelGProps::Perform(const gp_Cone& S, double Sn1 = std::sin(Alpha1); double Cn2 = std::cos(Alpha2); double Cn1 = std::cos(Alpha1); - double ZZ = (Z2 - Z1) * (Z2 - Z1) * Cnt * Snt; - double Auxi1 = 2 * R + (Z2 + Z1) * Snt; + double Auxi2 = (Z2 * Z2 + Z1 * Z2 + Z1 * Z1) / 3.; - dim = ZZ * (Alpha2 - Alpha1) * Auxi1 / 2.; + // Volume swept between the axis and the cone, the frustum for a full turn. + dim = (Alpha2 - Alpha1) * Cnt * (Z2 - Z1) + * (3 * R * R + 3 * R * (Z2 + Z1) * Snt + (Z2 * Z2 + Z2 * Z1 + Z1 * Z1) * Snt * Snt) / 6.; - double R1 = R + Z1 * Snt; - double R2 = R + Z2 * Snt; + double R1 = R + Z1 * Snt; + double R2 = R + Z2 * Snt; + double Z4 = (Z1 + Z2) * (Z1 * Z1 + Z2 * Z2); + double Z5 = Z2 * Z2 * Z2 * Z2 + Z2 * Z2 * Z2 * Z1 + Z2 * Z2 * Z1 * Z1 + Z2 * Z1 * Z1 * Z1 + + Z1 * Z1 * Z1 * Z1; double Coef0 = (R1 * R1 + R1 * R2 + R2 * R2); - double Iz = Cnt * (R * (Z2 + Z1) + 2 * Snt * (Z1 * Z1 + Z1 * Z2 + Z2 * Z2) / 3.) / Auxi1; - double Ix = Coef0 * (Sn2 - Sn1) / (Alpha2 - Alpha1) / Auxi1; - double Iy = Coef0 * (Cn1 - Cn2) / (Alpha2 - Alpha1) / Auxi1; + double Coef1 = (R1 * R1 * R1 + R1 * R1 * R2 + R1 * R2 * R2 + R2 * R2 * R2); + // int v r^2 dv with r = R + v sin(a). + double Coef2 = R * R * (Z1 + Z2) / 2. + 2. * R * Snt * Auxi2 + Snt * Snt * Z4 / 4.; + double Iz = 3. * Cnt * Coef2 / Coef0; + double Ix = Coef1 * (Sn2 - Sn1) / (Alpha2 - Alpha1) / (2. * Coef0); + double Iy = Coef1 * (Cn1 - Cn2) / (Alpha2 - Alpha1) / (2. * Coef0); g.SetCoord(X0 + Xa1 * Ix + Xa2 * Iy + Xa3 * Iz, Y0 + Ya1 * Ix + Ya2 * Iy + Ya3 * Iz, Z0 + Za1 * Ix + Za2 * Iy + Za3 * Iz); - double IR2 = ZZ * (R2 * R2 * R2 + R2 * R2 * R1 + R1 * R1 * R2 + R1 * R1 * R1) / 4.; + // Per unit angle, r = R + v sin(a), z = v cos(a), dV = rho drho dz: IR2 = int r^4 / 4 dz, + // IZ2 = int z^2 r^2 / 2 dz, IRZ = int z r^3 / 3 dz. + double IR2 = Cnt * (Z2 - Z1) + * (R2 * R2 * R2 * R2 + R2 * R2 * R2 * R1 + R2 * R2 * R1 * R1 + R2 * R1 * R1 * R1 + + R1 * R1 * R1 * R1) + / 20.; + double IRZ = Cnt * Cnt * (Z2 - Z1) + * (R * R * R * (Z1 + Z2) / 2. + 3. * R * R * Snt * Auxi2 + + 3. * R * Snt * Snt * Z4 / 4. + Snt * Snt * Snt * Z5 / 5.) + / 3.; double ICn2 = IR2 * (Alpha2 - Alpha1 + Cn2 * Sn2 - Cn1 * Sn1) / 2.; - double ISn2 = IR2 * (Alpha2 - Alpha1 + Cn2 * Sn2 - Cn1 * Sn1) / 2.; - double IZ2 = ZZ * Cnt * Cnt * (Alpha2 - Alpha1) - * (Z1 * Z1 * (R / 3 + Z1 * Snt / 4) + Z2 * Z2 * (R / 3 + Z2 * Snt / 4) - + Z1 * Z2 * (R / 3 + Z1 * Snt / 4 + Z2 * Snt / 4)); - double ICnSn = IR2 * (Cn2 * Cn2 - Cn1 * Cn1); - double ICnz = (Z1 + Z2) * ZZ * Coef0 * (Sn2 - Sn1) / 3; - double ISnz = (Z1 + Z2) * ZZ * Coef0 * (Cn1 - Cn2) / 3; + double ISn2 = IR2 * (Alpha2 - Alpha1 - Cn2 * Sn2 + Cn1 * Sn1) / 2.; + double IZ2 = Cnt * Cnt * Cnt * (Z2 - Z1) * (Alpha2 - Alpha1) + * (R * R * Auxi2 + R * Snt * Z4 / 2. + Snt * Snt * Z5 / 5.) / 2.; + double ICnSn = IR2 * (Sn2 * Sn2 - Sn1 * Sn1) / 2.; + double ICnz = IRZ * (Sn2 - Sn1); + double ISnz = IRZ * (Cn1 - Cn2); math_Matrix Dm(1, 3, 1, 3); Dm(1, 1) = ISn2 + IZ2; @@ -199,32 +214,13 @@ void GProp_VelGProps::Perform(const gp_Cone& S, Dm(1, 3) = Dm(3, 1) = -ICnz; Dm(3, 2) = Dm(2, 3) = -ISnz; - math_Matrix Passage(1, 3, 1, 3); - Passage(1, 1) = Xa1; - Passage(1, 2) = Xa2; - Passage(1, 3) = Xa3; - Passage(2, 1) = Ya1; - Passage(2, 2) = Ya2; - Passage(2, 3) = Ya3; - Passage(3, 1) = Za1; - Passage(3, 2) = Za2; - Passage(3, 3) = Za3; - - math_Jacobi J(Dm); - math_Vector V1(1, 3), V2(1, 3), V3(1, 3); - J.Vector(1, V1); - V1.Multiply(Passage, V1); - V1.Multiply(J.Value(1)); - J.Vector(2, V2); - V2.Multiply(Passage, V2); - V2.Multiply(J.Value(2)); - J.Vector(3, V3); - V3.Multiply(Passage, V3); - V3.Multiply(J.Value(3)); - - inertia = - gp_Mat(gp_XYZ(V1(1), V2(1), V3(1)), gp_XYZ(V1(2), V2(2), V3(2)), gp_XYZ(V1(3), V2(3), V3(3))); + // Dm is about the cone location in the cone axes: take it to global axes, to g, then to loc. + gp_Mat Pm(Xa1, Xa2, Xa3, Ya1, Ya2, Ya3, Za1, Za2, Za3); + gp_Mat + Dg(Dm(1, 1), Dm(1, 2), Dm(1, 3), Dm(2, 1), Dm(2, 2), Dm(2, 3), Dm(3, 1), Dm(3, 2), Dm(3, 3)); gp_Mat Hop; + GProp::HOperator(g, S.Location(), dim, Hop); + inertia = Pm * Dg * Pm.Transposed() - Hop; GProp::HOperator(g, loc, dim, Hop); inertia = inertia + Hop; } diff --git a/src/ModelingData/TKGeomBase/GTests/FILES.cmake b/src/ModelingData/TKGeomBase/GTests/FILES.cmake index 01fe6fd2891..70f8efe023c 100644 --- a/src/ModelingData/TKGeomBase/GTests/FILES.cmake +++ b/src/ModelingData/TKGeomBase/GTests/FILES.cmake @@ -61,6 +61,8 @@ set(OCCT_TKGeomBase_GTests_FILES GeomLib_CheckCurveOnSurface_Test.cxx GProp_PEquation_Test.cxx GProp_PGProps_Test.cxx + GProp_SelGProps_Test.cxx + GProp_VelGProps_Test.cxx Hermit_Test.cxx IntAna_IntQuadQuad_Test.cxx ProjLib_Cone_Test.cxx diff --git a/src/ModelingData/TKGeomBase/GTests/GProp_SelGProps_Test.cxx b/src/ModelingData/TKGeomBase/GTests/GProp_SelGProps_Test.cxx new file mode 100644 index 00000000000..455e7234213 --- /dev/null +++ b/src/ModelingData/TKGeomBase/GTests/GProp_SelGProps_Test.cxx @@ -0,0 +1,219 @@ +// Copyright (c) 2026 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +namespace +{ +//! Mass, first moment and second moment (integral of p p^T) of a patch about the origin. +struct PatchMoments +{ + double Mass = 0.0; + gp_XYZ First; + double Second[3][3] = {}; +}; + +//! Gauss-Legendre integral over [theU1, theU2] x [theV1, theV2]. +//! theSample returns the point at (u, v) and the area element there. +PatchMoments integratePatch(const std::function& theSample, + const double theU1, + const double theU2, + const double theV1, + const double theV2) +{ + constexpr int aNbNodes = 16; + math_Vector aNodes(1, aNbNodes), aWeights(1, aNbNodes); + math::OrderedGaussPointsAndWeights(aNbNodes, aNodes, aWeights); + + PatchMoments aRes; + for (int i = 1; i <= aNbNodes; ++i) + { + for (int j = 1; j <= aNbNodes; ++j) + { + const double aU = 0.5 * (theU2 - theU1) * aNodes(i) + 0.5 * (theU2 + theU1); + const double aV = 0.5 * (theV2 - theV1) * aNodes(j) + 0.5 * (theV2 + theV1); + double aDm = 0.0; + const gp_XYZ aP = theSample(aU, aV, aDm); + const double aW = 0.25 * (theU2 - theU1) * (theV2 - theV1) * aWeights(i) * aWeights(j) * aDm; + aRes.Mass += aW; + aRes.First += aW * aP; + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + aRes.Second[a][b] += aW * aP.Coord(a + 1) * aP.Coord(b + 1); + } + } + } + } + return aRes; +} + +//! Matrix of inertia about the centre of mass from the moments about the origin. +gp_Mat inertiaAboutCentre(const PatchMoments& theMom, gp_XYZ& theCentre) +{ + theCentre = theMom.First / theMom.Mass; + double aS[3][3]; + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + aS[a][b] = + theMom.Second[a][b] - theMom.Mass * theCentre.Coord(a + 1) * theCentre.Coord(b + 1); + } + } + const double aTrace = aS[0][0] + aS[1][1] + aS[2][2]; + gp_Mat aRes; + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + aRes.SetValue(a + 1, b + 1, (a == b ? aTrace : 0.0) - aS[a][b]); + } + } + return aRes; +} + +//! Moment of inertia about the axis through theQ along the unit vector theD. +double momentAboutAxis(const PatchMoments& theMom, const gp_XYZ& theQ, const gp_XYZ& theD) +{ + double aTrace = 0.0, aDSD = 0.0; + for (int a = 0; a < 3; ++a) + { + aTrace += theMom.Second[a][a]; + for (int b = 0; b < 3; ++b) + { + aDSD += theD.Coord(a + 1) * theMom.Second[a][b] * theD.Coord(b + 1); + } + } + const double aQd = theQ.Dot(theD); + return (aTrace - 2.0 * theQ.Dot(theMom.First) + theMom.Mass * theQ.SquareModulus()) + - (aDSD - 2.0 * aQd * theD.Dot(theMom.First) + theMom.Mass * aQd * aQd); +} + +//! Global point of the local coordinates (theX, theY, theZ) in the frame theAx. +gp_XYZ toGlobal(const gp_Ax3& theAx, const double theX, const double theY, const double theZ) +{ + return theAx.Location().XYZ() + theX * theAx.XDirection().XYZ() + theY * theAx.YDirection().XYZ() + + theZ * theAx.Direction().XYZ(); +} + +struct ConeCase +{ + double SemiAngle; + double RefRadius; + double V1; + double V2; + double U1; + double U2; +}; + +const double kTol = 1.0e-9; +} // namespace + +TEST(GProp_SelGPropsTest, Cone_LateralArea_ClosedForm) +{ + // Area element of gp_Cone is R + v sin(a), no cos(a) factor. + const double aSemiAngle = M_PI / 6.0, aRadius = 5.0, aV1 = 0.0, aV2 = 10.0; + GProp_SelGProps aProps; + aProps.Perform(gp_Cone(gp_Ax3(gp::XOY()), aSemiAngle, aRadius), 0.0, 2.0 * M_PI, aV1, aV2); + + const double aExpected = + 2.0 * M_PI * (aV2 - aV1) * (aRadius + 0.5 * (aV1 + aV2) * std::sin(aSemiAngle)); + EXPECT_NEAR(aProps.Mass(), aExpected, kTol * aExpected); +} + +TEST(GProp_SelGPropsTest, Cone_MatchesQuadrature) +{ + const double aPi = M_PI; + const ConeCase aCases[] = {{aPi / 6.0, 5.0, 0.0, 10.0, 0.0, 2.0 * aPi}, + {aPi / 12.0, 8.0, 0.5, 6.0, 0.0, 2.0 * aPi}, + {aPi / 3.0, 8.0, 0.5, 6.0, 0.0, 2.0 * aPi}, + {-aPi / 6.0, 8.0, 0.5, 6.0, 0.0, 2.0 * aPi}, + {aPi / 6.0, 5.0, 1.0, 9.0, 0.3, 2.2}, + {-aPi / 4.0, 9.0, 0.5, 7.0, 0.3, 4.0}}; + const gp_Ax3 aFrames[] = { + gp_Ax3(gp::XOY()), + gp_Ax3(gp_Pnt(1.5, -2.0, 3.0), gp_Dir(1.0, 2.0, 3.0), gp_Dir(2.0, -1.0, 0.0))}; + const gp_Pnt aLoc(2.0, -1.0, 3.0); + const gp_XYZ aAxisDir = gp_XYZ(1.0, 2.0, 2.0) / 3.0; + + for (const gp_Ax3& aFrame : aFrames) + { + for (const ConeCase& aCase : aCases) + { + const double aSin = std::sin(aCase.SemiAngle), aCos = std::cos(aCase.SemiAngle); + const gp_Cone aCone(aFrame, aCase.SemiAngle, aCase.RefRadius); + const PatchMoments aMom = integratePatch( + [&](double theU, double theV, double& theDm) { + const double aR = aCase.RefRadius + theV * aSin; + theDm = aR; + return toGlobal(aFrame, aR * std::cos(theU), aR * std::sin(theU), theV * aCos); + }, + aCase.U1, + aCase.U2, + aCase.V1, + aCase.V2); + + GProp_SelGProps aProps; + aProps.Perform(aCone, aCase.U1, aCase.U2, aCase.V1, aCase.V2); + + gp_XYZ aCentre; + const gp_Mat aInertia = inertiaAboutCentre(aMom, aCentre); + double aScale = 0.0; + for (int a = 1; a <= 3; ++a) + { + for (int b = 1; b <= 3; ++b) + { + aScale = std::max(aScale, std::abs(aInertia.Value(a, b))); + } + } + + EXPECT_NEAR(aProps.Mass(), aMom.Mass, kTol * aMom.Mass); + EXPECT_NEAR(aProps.CentreOfMass().Distance(gp_Pnt(aCentre)), + 0.0, + kTol * (1.0 + aCentre.Modulus())); + const gp_Mat aMat = aProps.MatrixOfInertia(); + for (int a = 1; a <= 3; ++a) + { + for (int b = 1; b <= 3; ++b) + { + EXPECT_NEAR(aMat.Value(a, b), aInertia.Value(a, b), kTol * aScale); + } + } + + // The inertia about an axis through the location tests the Huyghens shift to it. + GProp_SelGProps aLocProps(aCone, aCase.U1, aCase.U2, aCase.V1, aCase.V2, aLoc); + const double aExpectedAxis = momentAboutAxis(aMom, aLoc.XYZ(), aAxisDir); + EXPECT_NEAR(aLocProps.MomentOfInertia(gp_Ax1(aLoc, gp_Dir(aAxisDir))), + aExpectedAxis, + kTol * std::abs(aExpectedAxis)); + } + } +} diff --git a/src/ModelingData/TKGeomBase/GTests/GProp_VelGProps_Test.cxx b/src/ModelingData/TKGeomBase/GTests/GProp_VelGProps_Test.cxx new file mode 100644 index 00000000000..1dfc2bde10d --- /dev/null +++ b/src/ModelingData/TKGeomBase/GTests/GProp_VelGProps_Test.cxx @@ -0,0 +1,241 @@ +// Copyright (c) 2026 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +namespace +{ +//! Mass, first moment and second moment (integral of p p^T) of a solid about the origin. +struct SolidMoments +{ + double Mass = 0.0; + gp_XYZ First; + double Second[3][3] = {}; +}; + +//! Gauss-Legendre integral over [theU1, theU2] x [theV1, theV2] x [0, 1]. +//! theSample returns the point at (u, v, t) and the volume element there. +SolidMoments integrateSolid(const std::function& theSample, + const double theU1, + const double theU2, + const double theV1, + const double theV2) +{ + constexpr int aNbNodes = 16; + math_Vector aNodes(1, aNbNodes), aWeights(1, aNbNodes); + math::OrderedGaussPointsAndWeights(aNbNodes, aNodes, aWeights); + + SolidMoments aRes; + for (int i = 1; i <= aNbNodes; ++i) + { + for (int j = 1; j <= aNbNodes; ++j) + { + for (int k = 1; k <= aNbNodes; ++k) + { + const double aU = 0.5 * (theU2 - theU1) * aNodes(i) + 0.5 * (theU2 + theU1); + const double aV = 0.5 * (theV2 - theV1) * aNodes(j) + 0.5 * (theV2 + theV1); + const double aT = 0.5 * aNodes(k) + 0.5; + double aDm = 0.0; + const gp_XYZ aP = theSample(aU, aV, aT, aDm); + const double aW = + 0.125 * (theU2 - theU1) * (theV2 - theV1) * aWeights(i) * aWeights(j) * aWeights(k) * aDm; + aRes.Mass += aW; + aRes.First += aW * aP; + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + aRes.Second[a][b] += aW * aP.Coord(a + 1) * aP.Coord(b + 1); + } + } + } + } + } + return aRes; +} + +//! Matrix of inertia about the centre of mass from the moments about the origin. +gp_Mat inertiaAboutCentre(const SolidMoments& theMom, gp_XYZ& theCentre) +{ + theCentre = theMom.First / theMom.Mass; + double aS[3][3]; + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + aS[a][b] = + theMom.Second[a][b] - theMom.Mass * theCentre.Coord(a + 1) * theCentre.Coord(b + 1); + } + } + const double aTrace = aS[0][0] + aS[1][1] + aS[2][2]; + gp_Mat aRes; + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + aRes.SetValue(a + 1, b + 1, (a == b ? aTrace : 0.0) - aS[a][b]); + } + } + return aRes; +} + +//! Moment of inertia about the axis through theQ along the unit vector theD. +double momentAboutAxis(const SolidMoments& theMom, const gp_XYZ& theQ, const gp_XYZ& theD) +{ + double aTrace = 0.0, aDSD = 0.0; + for (int a = 0; a < 3; ++a) + { + aTrace += theMom.Second[a][a]; + for (int b = 0; b < 3; ++b) + { + aDSD += theD.Coord(a + 1) * theMom.Second[a][b] * theD.Coord(b + 1); + } + } + const double aQd = theQ.Dot(theD); + return (aTrace - 2.0 * theQ.Dot(theMom.First) + theMom.Mass * theQ.SquareModulus()) + - (aDSD - 2.0 * aQd * theD.Dot(theMom.First) + theMom.Mass * aQd * aQd); +} + +//! Global point of the local coordinates (theX, theY, theZ) in the frame theAx. +gp_XYZ toGlobal(const gp_Ax3& theAx, const double theX, const double theY, const double theZ) +{ + return theAx.Location().XYZ() + theX * theAx.XDirection().XYZ() + theY * theAx.YDirection().XYZ() + + theZ * theAx.Direction().XYZ(); +} + +struct ConeCase +{ + double SemiAngle; + double RefRadius; + double V1; + double V2; + double U1; + double U2; +}; + +const double kTol = 1.0e-9; +} // namespace + +TEST(GProp_VelGPropsTest, Cone_Volume_IsTheFrustum) +{ + // Full turn: pi H (R1^2 + R1 R2 + R2^2) / 3 with H = (v2 - v1) cos(a). + const double aSemiAngle = M_PI / 6.0, aRadius = 5.0, aV1 = 0.0, aV2 = 10.0; + GProp_VelGProps aProps; + aProps.Perform(gp_Cone(gp_Ax3(gp::XOY()), aSemiAngle, aRadius), 0.0, 2.0 * M_PI, aV1, aV2); + + const double aR1 = aRadius + aV1 * std::sin(aSemiAngle); + const double aR2 = aRadius + aV2 * std::sin(aSemiAngle); + const double aH = (aV2 - aV1) * std::cos(aSemiAngle); + const double aExpected = M_PI * aH * (aR1 * aR1 + aR1 * aR2 + aR2 * aR2) / 3.0; + EXPECT_NEAR(aProps.Mass(), aExpected, kTol * aExpected); +} + +TEST(GProp_VelGPropsTest, Cone_NearCylinder_VolumeMatchesCylinder) +{ + GProp_VelGProps aCylProps; + aCylProps.Perform(gp_Cylinder(gp_Ax3(gp::XOY()), 5.0), 0.0, 2.0 * M_PI, 0.0, 10.0); + GProp_VelGProps aConeProps; + aConeProps.Perform(gp_Cone(gp_Ax3(gp::XOY()), 1.0e-6, 5.0), 0.0, 2.0 * M_PI, 0.0, 10.0); + + EXPECT_NEAR(aConeProps.Mass(), aCylProps.Mass(), 1.0e-4 * aCylProps.Mass()); +} + +TEST(GProp_VelGPropsTest, Cone_MatchesQuadrature) +{ + const double aPi = M_PI; + const ConeCase aCases[] = {{aPi / 6.0, 5.0, 0.0, 10.0, 0.0, 2.0 * aPi}, + {aPi / 12.0, 8.0, 0.5, 6.0, 0.0, 2.0 * aPi}, + {aPi / 3.0, 8.0, 0.5, 6.0, 0.0, 2.0 * aPi}, + {-aPi / 6.0, 8.0, 0.5, 6.0, 0.0, 2.0 * aPi}, + {aPi / 6.0, 5.0, 1.0, 9.0, 0.3, 2.2}, + {-aPi / 4.0, 9.0, 0.5, 7.0, 0.3, 4.0}}; + const gp_Ax3 aFrames[] = { + gp_Ax3(gp::XOY()), + gp_Ax3(gp_Pnt(1.5, -2.0, 3.0), gp_Dir(1.0, 2.0, 3.0), gp_Dir(2.0, -1.0, 0.0))}; + const gp_Pnt aLoc(2.0, -1.0, 3.0); + const gp_XYZ aAxisDir = gp_XYZ(1.0, 2.0, 2.0) / 3.0; + + for (const gp_Ax3& aFrame : aFrames) + { + for (const ConeCase& aCase : aCases) + { + const double aSin = std::sin(aCase.SemiAngle), aCos = std::cos(aCase.SemiAngle); + const gp_Cone aCone(aFrame, aCase.SemiAngle, aCase.RefRadius); + const SolidMoments aMom = integrateSolid( + [&](double theU, double theV, double theT, double& theDm) { + // Solid swept between the axis and the patch: rho = t r, dV = rho drho du dz. + const double aR = aCase.RefRadius + theV * aSin; + theDm = theT * aR * aR * aCos; + return toGlobal(aFrame, + theT * aR * std::cos(theU), + theT * aR * std::sin(theU), + theV * aCos); + }, + aCase.U1, + aCase.U2, + aCase.V1, + aCase.V2); + + GProp_VelGProps aProps; + aProps.Perform(aCone, aCase.U1, aCase.U2, aCase.V1, aCase.V2); + + gp_XYZ aCentre; + const gp_Mat aInertia = inertiaAboutCentre(aMom, aCentre); + double aScale = 0.0; + for (int a = 1; a <= 3; ++a) + { + for (int b = 1; b <= 3; ++b) + { + aScale = std::max(aScale, std::abs(aInertia.Value(a, b))); + } + } + + EXPECT_NEAR(aProps.Mass(), aMom.Mass, kTol * aMom.Mass); + EXPECT_NEAR(aProps.CentreOfMass().Distance(gp_Pnt(aCentre)), + 0.0, + kTol * (1.0 + aCentre.Modulus())); + const gp_Mat aMat = aProps.MatrixOfInertia(); + for (int a = 1; a <= 3; ++a) + { + for (int b = 1; b <= 3; ++b) + { + EXPECT_NEAR(aMat.Value(a, b), aInertia.Value(a, b), kTol * aScale); + } + } + + // The inertia about an axis through the location tests the Huyghens shift to it. + GProp_VelGProps aLocProps(aCone, aCase.U1, aCase.U2, aCase.V1, aCase.V2, aLoc); + const double aExpectedAxis = momentAboutAxis(aMom, aLoc.XYZ(), aAxisDir); + EXPECT_NEAR(aLocProps.MomentOfInertia(gp_Ax1(aLoc, gp_Dir(aAxisDir))), + aExpectedAxis, + kTol * std::abs(aExpectedAxis)); + } + } +}